Fixing tendencies in shallow water model and bringing back regression tests

[francispoulin]

Following up on #3394, this PR tries to do the first two tasks, fix the tendencies in the shallow water model and bring back the regression tests.

[francispoulin]

I'm trying to run the near global shallow water model and having difficulties with the latitude longitude grid. The line that calls it is here

https://github.com/CliMA/Oceananigans.jl/blob/5ba82e5786d1b7725d6bfa3e45072d2f6fad487a/validation/shallow_water_model/near_global_shallow_water_quarter_degree.jl#L95

And the error that I get is the following,

julia> include("near_global_shallow_water_quarter_degree.jl")
┌ Warning: Over-writing registration of the datadep
│   name = "quarter_degree_near_global_lat_lon"
└ @ DataDeps ~/.julia/packages/DataDeps/Y2lje/src/registration.jl:15
ERROR: LoadError: UndefVarError: `Nλ` not defined
Stacktrace:
 [1] validate_lat_lon_grid_args(topology::Tuple{…}, size::Tuple{…}, halo::Tuple{…}, FT::Type, latitude::Tuple{…}, longitude::Tuple{…}, z::Nothing, precompute_metrics::Bool)
   @ Oceananigans.Grids ~/Software/Oceananigans.jl/src/Grids/latitude_longitude_grid.jl:257
 [2] LatitudeLongitudeGrid(architecture::GPU, FT::DataType; size::Tuple{…}, longitude::Tuple{…}, latitude::Tuple{…}, z::Nothing, radius::Float64, topology::Tuple{…}, precompute_metrics::Bool, halo::Tuple{…})
   @ Oceananigans.Grids ~/Software/Oceananigans.jl/src/Grids/latitude_longitude_grid.jl:189
 [3] macro expansion
   @ show.jl:1181 [inlined]
 [4] top-level scope
   @ ~/Software/Oceananigans.jl/validation/shallow_water_model/near_global_shallow_water_quarter_degree.jl:95
 [5] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [6] top-level scope
   @ REPL[70]:1
 [7] top-level scope
   @ ~/.julia/packages/CUDA/nbRJk/src/initialization.jl:205
in expression starting at /home/fpoulin/Software/Oceananigans.jl/validation/shallow_water_model/near_global_shallow_water_quarter_degree.jl:95
Some type information was truncated. Use `show(err)` to see complete types.

Should the lat-lon grid be able to handle flat in the vertical? @simone-silvestri, maybe this is something you tried to fix the last time we chatted?

[simone-silvestri]

It definitely should. Indeed the previous near global simulation with the shallow water model was done on a lat lon grid flat in the vertical. Something must have changed. I will take a look

[navidcy]

on #main:

julia> grid = LatitudeLongitudeGrid(size = (10, 10), latitude = (-60, 60), longitude = (-180, 180), topology = (Periodic, Bounded, Flat))
10×10×1 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Flat} on CPU with 3×3×0 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [-180.0, 180.0)  regularly spaced with Δλ=36.0
├── latitude:  Bounded  φ ∈ [-60.0, 60.0]    regularly spaced with Δφ=12.0
└── z:         Flat z

but on this branch:

julia> grid = LatitudeLongitudeGrid(size = (10, 10), latitude = (-60, 60), longitude = (-180, 180), topology = (Periodic, Bounded, Flat))
ERROR: UndefVarError: `Nλ` not defined
Stacktrace:
 [1] validate_lat_lon_grid_args(topology::Tuple{…}, size::Tuple{…}, halo::Nothing, FT::Type, latitude::Tuple{…}, longitude::Tuple{…}, z::Nothing, precompute_metrics::Bool)
   @ Oceananigans.Grids ~/Research/OC10.jl/src/Grids/latitude_longitude_grid.jl:257
 [2] LatitudeLongitudeGrid(architecture::CPU, FT::DataType; size::Tuple{…}, longitude::Tuple{…}, latitude::Tuple{…}, z::Nothing, radius::Float64, topology::Tuple{…}, precompute_metrics::Bool, halo::Nothing)
   @ Oceananigans.Grids ~/Research/OC10.jl/src/Grids/latitude_longitude_grid.jl:189
 [3] top-level scope
   @ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.

So it has something to do with the changes here.

[francispoulin]

grid = LatitudeLongitudeGrid(size = (10, 10), latitude = (-60, 60), longitude = (-180, 180), topology = (Periodic, Bounded, Flat))

Yes. If you change the line to add a 1 to the size as follows, it does run. However, ShallowWaterModel does no tallow for a 3-tuple for size so this is not an option here.

grid = LatitudeLongitudeGrid(size = (10, 10, 1), latitude = (-60, 60), longitude = (-180, 180), topology = (Periodic, Bounded, Flat))

[francispoulin]

I am able to run the validatoin example. Unfortunately, it produces NaNs after just over an hour. Sigh.

Even if I reduce the time step from 20 seconds to 1 second, it still goes unstable at about the same time. I need to look more closely as to why this might be happening. I would share an animation but one hour doesn't seem to be enough time for much to happen so I don't see any changes.

[navidcy]

I remember there is a problematic issue with initialisation and bathymetry. You need to be extra careful to put initial condition that is consistent with the bathymetry. We should actually fix that because it was so counterintuitive. That’s what I was trying to understand back in the day and fix but I lost steam. I couldn’t set up very simple things, let alone global configuration.

[francispoulin]

Very good point @navidcy!

Instead of trying to do this first, it makes more sense to include simple topography.

I'll start with the example and include the affect of shelflike topography, as Glenn and I did in this paper. It should be easy to set up and will teach us how to include topography in a much simpler setting.

https://journals.ametsoc.org/view/journals/phoc/35/5/jpo2719.1.xml?tab_body=abstract-display

[simone-silvestri]

There must be a simple bug to fix. I remember the simulation was quite stable. I think it has to do with how bathymetry is defined in the code. A simple configuration might snoop out the issue

[simone-silvestri]

I think in this PR we changed the definition of bathymetry to how it was defined before. Also, the initial height cannot be zero (because of wind stress forcing begin $\tau / h$). On main, it looks like it runs stably (correcting the momentum advection scheme that changed in the meantime)

[francispoulin]

I have created a new validation script called shallow_jet_topography.jl that evolves the barotropically unstable jet (like in the shallow water example) but now over a shelflike topography. It can be found in valiation/shallow_water.

Some preliminary results are shown below.

First, this is the case with a flat bottom and it does produce the same results as what we have in the example.

Second, this is an example of prograde topography with the shallow water on the right of the propogating jet.

Third, this is an example of retrograde topopgraphy with the shallow water on the right of the propogating jet.

The growth rates differ but I want to update my linear stability code to compute the growth rates so that we can validate that the results are quantitatively close.

But for this case here the prograde case is slightly more unstable with the rettrograde case, which is qualitatively correct.

[francispoulin]

This is a sample plot of the growth rate in the case of a flat bottom, and it agrees with what I expecte to see. Also the growth rate is 0.139, which matches with the theory to at least three significant digits.

This is what I get with some topography, it doesn't matter which direction, they are both very similar to this. The large oscillations at the beginning are a cluse that something is wrong. There is a rapid initial growth that shouldn't be there either.

This is using the conservative form. I'll try the vector invariant form and see if this happens in both cases or not.

[navidcy]

I have created a new validation script called shallow_jet_topography.jl that evolves the barotropically unstable jet (like in the shallow water example) but now over a shelflike topography. It can be found in valiation/shallow_water.

I don't see it; did you push?

[francispoulin]

I have created a new validation script called shallow_jet_topography.jl that evolves the barotropically unstable jet (like in the shallow water example) but now over a shelflike topography. It can be found in valiation/shallow_water.

I don't see it; did you push?

Sorry! I have added it now.

https://github.com/CliMA/Oceananigans.jl/blob/fp-ss/shallow-water-version2/validation/shallow_water_model/shallow_water_jet_topography.jl

One of my mitakes is that I didn't define the bathymetry. I need to figure out how to do that. Thanks for the very helpful feedback!

[simone-silvestri]

you need to pass it to the model as a keyword argument. In this case with bathymetry = b

[francispoulin]

you need to pass it to the model as a keyword argument. In this case with bathymetry = b

Nice! I will give it a try today.

[francispoulin]

This is the case of prograde topography, which has a growth rate that looks reasoanble and is 0.147 > 0.139 as we have in the flat bottom case. This is consistent with theory.

This is the case of retrograde topography, which has a growth rate that looks reasoanble and is 0.115 < 0.139 as we have in the flat bottom case. This is consistent with theory.

This is all very promising and I will confirm we get the same values with the vector invariant form.

[francispoulin]

I am happy to say that the VectorInvariant formulation yields similar growth rates 0.144 and 0.123. Given the numerical methods are different I'm not bothered by this difference. They should converge with increased resolution.

Any more questions/concerns on this simple validation problem before I return to the global one?

[navidcy]

I have a bit of trouble and confusion regarding conventions for setting up the bathymetry. I'll put a few comments on the script and hopefully that will showcase where my confusion is.

[francispoulin]

Thanks everyone for the feedback. I did some debugging and I'm happy to say that bathymetry is now specified as we expected, no negative sign, and it gives results that are consistent with previous findings: prograde is more unstable and retrograde is more stable.

You can see the code here.

Any other thoughts on this before I move on to the global model?

[simone-silvestri]

So in the end bathymetry is defined as a negative distance from the surface?

[francispoulin]

So in the end bathymetry is defined as a negative distance from the surface?

Yes, exactly that.

I still think that we should change h_plus_hB in the tendencies to free_surface or something, as it is exactly that and more descriptive.

[simone-silvestri]

So in the end bathymetry is defined as a negative distance from the surface?

Yes, exactly that.

I still think that we should change h_plus_hB in the tendencies to free_surface or something, as it is exactly that and more descriptive.

Sure!

[francispoulin]

Lots of the checks have failed and I don't understand why. Anyone have an idea as to what's going wrong?

[francispoulin]

Thanks @glwagner for merging main into this branch.

I have returned to the global simulations and have learned something.

Using all the same parameters as before, I get NaNs after just over an hour.

Upon closer inspection, I see that the height becomes negative after 56 minutes or so. I have tried reducing the time step by a factor of 20 and that didn't help.

Do you think I should play around with other advection schemes, say some high order upwinding schemes?

[navidcy]

Thanks @glwagner for merging main into this branch.

I have returned to the global simulations and have learned something.

Using all the same parameters as before, I get NaNs after just over an hour.

Upon closer inspection, I see that the height becomes negative after 56 minutes or so. I have tried reducing the time step by a factor of 20 and that didn't help.

Do you think I should play around with other advection schemes, say some high order upwinding schemes?

Seems to me (but I don't know for sure) that the problem is more basic than the advection scheme... You sure that the initial condition is consistent with bathymetry?

[simone-silvestri]

Thanks @glwagner for merging main into this branch.

I have returned to the global simulations and have learned something.

Using all the same parameters as before, I get NaNs after just over an hour.

Upon closer inspection, I see that the height becomes negative after 56 minutes or so. I have tried reducing the time step by a factor of 20 and that didn't help.

Do you think I should play around with other advection schemes, say some high order upwinding schemes?

You need to make sure that h is never zero, that will give you a NaN because the forcing is divided by h. It is kind of annoying but you have to set a "minimum height" in the initial conditions (cannot start with h = 0 somewhere). I remember setting hmin = 10 was enough to get a stable simulation

[francispoulin]

Thanks @glwagner for merging main into this branch. I have returned to the global simulations and have learned something. Using all the same parameters as before, I get NaNs after just over an hour. Upon closer inspection, I see that the height becomes negative after 56 minutes or so. I have tried reducing the time step by a factor of 20 and that didn't help. Do you think I should play around with other advection schemes, say some high order upwinding schemes?

You need to make sure that h is never zero, that will give you a NaN because the forcing is divided by h. It is kind of annoying but you have to set a "minimum height" in the initial conditions (cannot start with h = 0 somewhere). I remember setting hmin = 10 was enough to get a stable simulation

Thanks @simone-silvestri.

The validation script was set up with hmin = 10 and that did go unstable. I tried 20 and 30 and they both go unstable. I'm tring 40 and it seems better. I thought if this at least runs, and shows us the ACC going on the right direction, thsi would be a big leap forward. Then we can worry abou thow to minimize the minimum h, if we like.

I will keep you posted.

Update: 50 is still unstable after just under 3 hours. I'm now going to try 100. It seems like a pretty big cushion but I am curious what it takes.

And 100 yields negative height after less than 4 hours. Not good.

[simone-silvestri]

Does the script crash on main?

[navidcy]

Smells like there might be something else happening also here @francispoulin...

[francispoulin]

Good question @simone-silvestri !

On main it doesn't run as is. There are two issues.

First, there is device!(2) in line 30 but I commented that out easily.

Second, on line 159 we load VorticityStencil, which is not defined. I deleted that and changed the momentum_advection to VectorInvariant(), but it occurs to me that this is not what was done before. Before it was using a vorticity stencil. How would I change to use the other formulation?

However, when I try running this on main it goes far beyond what I see in the PR. So I guess the version on main does run and is stable. Pasat a day easy. Hmm...

If we have a problem with the height going negative, I wonder if this is a sign that the conservation of mass equation is not quite right? I will look at that today.