Fix Spherical Grid capability (#78)

[CyprienBosserelle]

Fix Spherical grid correction

Spherical grid capability has been shaky!. Added to a previous branch (#78) still had a serious bug. This may have been solved in the development branch.
This branch tris to bring it all together and resolve any remaining bugs.

Basework

• Modify Calcres for spherical grids
• Add cm and fmv correction to Kurganov scheme
• Add cm and fmv correction to Buttinger scheme
• Add spherical coordinate correction to the advance scheme
• Add fix for adaptive BUQ
• Add delta reference to new gradient
• Make spherical/geo keyword a bool input
• Double check validity of Gradients
• Check which engine/forcing can/cannot use spherical grid
• Add limiter for latitude to stay between -90/90
• Hotstarting Pressure and adding inv. baro pressure to boundary water level forcing

Test

• Lake-at-rest with spherical adaptive grid

Demo

• Tsunami Comparison of projected model and spherical model
• storm surgeComparison of projected model and spherical model
• rain Comparison of projected model and spherical model
• rivers Comparison of projected model and spherical model

[CyprienBosserelle]

Check forcing in spherical mode

• Rain is OK because it affects h and is area independent
• River Needs work to include cm
• Wind?
• Atm pressure

[CyprienBosserelle]

Does the Spherical code and the projected code produce similar results?

BG_Flood can run on a spherical grid or a crtesian (project) grid. One of the sanity check for validating the code is to run a similar simulation both in projected coordinate and in spherical coordinates

Comparing apples and oranges

because the calculations are somewhat different on a sphere then on a projected grid we can't really assume that both model produce exactly identical results. However we can check if thay are sensibly similar.

Samoa Tsunami test

The first test is using a constant resolution: 100 m in cartesian and 0.001 degree in spherical wich is 10ish % larger. As you can see theyr are comparable resolution but not exactly the same.
secondly we use an adaptive run with the same criteria but that will lead to inconsistent grid resolution.

Adaptation = Inrange,-10.0,-5.0,zb;

Cartesian domain

UTM model

Spherical model

Adaptive domain

Here you can see that the small differences in the grid layout lead to different adaptation even though the original data is different

UTM model

Spherical model

Also working for Buttinger reconstruction

Runtime consideration

Using the example above we can compare whether there is an advantage of using the spherical vs the projected coordinate.

All the test are done on my NVidia Quadro P620.

Simulation	Resolution	Number of Blocks / nodes	Runtime	rt/block [s/block]
Projected cartesian (UTM)	200m	1250 / 320,000	310	0.25
Spherical cartesian	0.002 deg	1058 / 270,848	286	0.27
UTM adaptive	100-800m	433 / 110,848	132	0.30
Spherical adaptive	0.001-0.008 deg	375 / 96,000	143	0.38

As you can see the computational cost of using spherical grid is about 8% for a cartesian grid, which isn't much compared with the adaptivity cost (~20%). But the combined cost (adaptivity+spherical) is much more with 52% slower per block. That is because there are additional steps/calculation for adaptive run in spherical mode.

Summary

Corrections for spherical grids does what is expected for tsunami and essenytially produce the same results.

However, when running small domain in spherical grid you might end up spending more time in computation (>25%) then if you project your DEM in projected coordinate.

In many case, the spherical solution is the only way to simulate the hazard (e.g. trans-pacific tsunami). in this case, make sure you think about your adaptation strategy. if you can reduce the number of nodes to 40% less than a cartesian grid your model will run faster!

[CyprienBosserelle]

Verification of Atm Press forcing

Running a UTM and spherical model with similar forcing (projected to UTM to be as similar as possible) show that both models are essentially similar. Minor difference are due to refinement difference and possibly how we treat the roughness in spherical mode (?)

Spherical:

UTM:

Verification of Wind forcing

As above they are both essentially similar:

Spherical:

UTM:

[CyprienBosserelle]

River check

Spherical implementation should also work . We can compare the volume discharge from a river in a UTM grid vs from a Spherical grid. The table show that the simulated injected water is very similar to the theoretical volume. The errors here are likely due to the steep slope the water is injected to which causes mass to be imbalanced (see pic below).

Theoretical volume injected	UTM simulated volume	Spherical simulated volume
1.2e6	1.1999e6 (-0.003%)	1.2002 (+0.02%)

[CyprienBosserelle]

Warm-starting with atmospheric pressure

During the investigation in veryfying atmospheric pressure we found some inconsistencies in the first few step in the model above. This is because the model is initialised with flat water and at the first step lifted/pushed by atmospheric pressure. In addition, boundaries are fighting against the atmospheric pressure forcing for a balance of boundary.

warm stating water level to match inverse Barometer forcing

I order to warm up the model one shouldn't use flat ocean as a warm up but instead use an inverse barometer calculation.

Testing on steady state convergence

This seem easy enough and should improve the model convergence. But it doesn't! This is because the inverse barometer may be suggesting that water level is significantly higher/lower than the boundary is trying to acheive. After 10 or 20 minute the model ends up in a worse shape than the cold start.
In the example below is modelled water level after 600s. The top panel is warm-started (see above image) but the boundaries left alone; in the lower panel the model is initiated with zero water level and boundaries are left alone.

Inverse Barometer forcing on the boundary

If we warm start the model and modify the boundary to match the inverse barometer behaviour we do get rid of boundary generated waves and the model warm start brings it to pretty-much steady state situation. Below picture is as above where the top panel is the model after 600s with the inverse barometer forcing the warm start and the boundary. lower panel is the model after 600s with not inverse barometer forcing the warm start or the boundary.

Automating warm start and boundary forcing makes sense

Letting BG_Flood deal with the inverse barometer warm start and boundary adjustment makes a lot of sense because in most cases the modeller would have to generate those. There are no prescribed switch to turn this feature off at this stage but a modeller could hack their pressure forcing to avoid the problem. Modifying the value in a 3d Netcdf file is far easier than making inverse barometer adjusted boundaries.

[AliceHarang]

This is a very good and useful development. It is well documented and tested.
As discussed, a continuous integration test might be added for the spherical grid and perhaps for the hotstart too (with pressure and without?)... but it can also go in the "completing the tests" part for a later PR!