Adapt cell_advection_timescale for ShallowWaterModel

[francispoulin]

I will look into this soon as I really want to avoid NaN's, as I'm getting today.

[glwagner]

What is the error you are getting when you try to use it?

The only function you'd need to extend is update_Δt! I think, which is probably a small piece of work if we know what error we're getting.

https://github.com/CliMA/Oceananigans.jl/blob/35f749ef09fe021d13e8bad6433d889e8c39d6ca/src/Simulations/time_step_wizard.jl#L32

[francispoulin]

I had not tried it before but I just tried it now.

The error is that ShallowWaterModel has no field velocities.

cell_advection_timescale.jl defined here meeds to know the velocities, u,v,w. In ShallowWaterModel we know uh, vh, and h, from which we can compute the velocities.

Is this a time for dispatching?

[glwagner]

Hmm, I think then cell_advection_timescale and cell_diffusion_timescale need to be extended for ShallowWaterModel, after which the TimeStepWizard will probably work.

Perhaps change the name of this issue to "Adapt cell_advection_timescale for ShallowWaterModel"?

[glwagner]

And just to clarify, when I say "extend" I'm referring to "method extension" --- in other words, defining new methods for a function that dispatch on the number or arguments and argument types.

[francispoulin]

And just to clarify, when I say "extend" I'm referring to "method extension" --- in other words, defining new methods for a function that dispatch on the number or arguments and argument types.

That makes sense and it should be easy enough to define ShallowWaterModel versions of these. I will work on the advection case and let you know what I come up with.

[francispoulin]

This is my first attempt at making a version of cell_advection_timescale for ShallowWaterModel. Note that we are dividing elements of mass flux and height, which are not defined at the same cell points, but if we are looking for estimates I wonder if this will be sufficient.

I have not tested this as I'm not sure how to integrate this into the rest of the code, but I thought we could try that after we agree on what the script should look like.

"Returns the time-scale for advection on a regular grid across a single grid cell 
 for ShallowWaterModel."
function cell_advection_timescale(uh, vh, h, grid)
    umax = maximum(abs, uh / h)
    vmax = maximum(abs, vh / h)

    Δx = grid.Δx
    Δy = grid.Δy

    return min(Δx/umax, Δy/vmax)
end

cell_advection_timescale(model) =
    cell_advection_timescale(model.solution.uh.data.parent,
                             model.solution.vh.data.parent,
                             model.solution.h.data.parent,
                             model.grid)

[glwagner]

You could maybe form a conservative estimate using

uhmax = maximum(abs, uh)
vhmax = maximum(abs, vh)
hmax = maximum(abs, h)

umax = max(uhmax / hmax, vhmax / hmax)

but this would be overly conservative for very nonlinear problems, I suppose.

One issue is that writing uh ./ h will allocate memory. Note that you need ./ (with the .), not just /.

You'll also need to dispatch on model::ShallowWaterModel.

[francispoulin]

Good point. I would rather be conservative then get Nan's so I will make that fix, using hmin. It occurs to me that if/when we have vanishing layer depths, this could get more complicated, but will save that concern for later, when it's relevant.

I belive these lines of code will fix it.

    uhmax = maximum(abs, uh)
    vhmax = maximum(abs, vh)
     hmin = minimum(abs,  h)

     umax = uhmax/hmin
     vmax = vhmax/hmin

[glwagner]

using hmin

ah yes good catch.

I belive these lines of code will fix it.

You also need to dispatch on ShallowWaterModel.

[francispoulin]

Agreed. Unfortunately, I'm not entirely clear how to dispatch. It will presumably require I also change something in IncompressibleModel,

If can give me high level descriptions of what needs to be done where, I can give it a try later today.

[glwagner]

Unfortunately, I'm not entirely clear how to dispatch

Rather than

cell_advection_timescale(model) =
    cell_advection_timescale(model.solution.uh.data.parent,
                             model.solution.vh.data.parent,
                             model.solution.h.data.parent,
                             model.grid)

as posted abvove, you'll need to write

cell_advection_timescale(model::ShallowWaterModel) =
    cell_advection_timescale(model.solution.uh.data.parent,
                             model.solution.vh.data.parent,
                             model.solution.h.data.parent,
                             model.grid)

[francispoulin]

thanks. And presumably a similar change needs to be done in the original script to include IncompressibleModel. That I think I can find and do. The last thing I can think of is where to define the function before there is a dispatch that occurs depending on the model.

[glwagner]

And presumably a similar change needs to be done in the original script to include IncompressibleModel.

No, you don't need to do that. But you can if you want.

[francispoulin]

The fewer changes the better so happy to keep this as is and will give it a try.

[francispoulin]

I figure you need to include this file somewhere. I moved the file into ShallowWaterModels folder and then modified this file here

software/Oceananigans.jl/src/Models/ShallowWaterModels/ShallowWaterModels.jl

to have the following line

include("cell_advection_timescale_shallow_water.jl")

unfortunately, it still complains that ShallowWaterModel does not have a field velocities.

[francispoulin]

Also, do we prefer the name cell_advection_timescale_shallow_water.jl or shallow_water_cell_advection_timescale.jl?

[glwagner]

I prefer shallow_water_cell_advection_timescale.jl.

At the top of the file, write

import Oceananigans.Utils: cell_advection_timescale

[francispoulin]

Thanks. I agree that I need this line and it creates cell_advection_timescale but thenI need to include the file I created. Should I do that in this same file right after the import?

[glwagner]

You should write

import Oceananigans.Utils: cell_advection_timescale

at the top of the file shallow_water_cell_advection_timescale.jl.

This "imports" the name cell_advection_timescale into your present scope such that the function can be extended with new methods. Next, you want to define a new method for cell_advection_timescale with the signature

cell_advection_timescale(model::ShallowWaterModel) = # definition here

Note that the function signature

cell_advection_timescale(uh, vh, h, grid)

is already taken by one of the existing methods for cell_advection_timescale:

https://github.com/CliMA/Oceananigans.jl/blob/2016e730350e0b6b625a5cf85a68269a370686fd/src/Utils/cell_advection_timescale.jl#L2-L12

This means that you need to either use type annotations to distinguish the version that takes arguments uh, vh, h, grid from the method that takes u, v, w, grid, or (probably better) use a unique name. A good name for the method that takes the arguments uh, vh, h, grid could be

shallow_water_cell_advection_timescale(uh, vh, h, grid) = # definition here

Remember that with this new name, the function cell_advection_timescale(model::ShallowWaterModel) must call shallow_water_cell_advection_timescale; eg

cell_advection_timescale(model::ShallowWaterModel) =
    shallow_water_cell_advection_timescale(model.solution.uh.data.parent,
                             model.solution.vh.data.parent,
                             model.solution.h.data.parent,
                             model.grid)

when this file is put together, you should include it in ShallowWaterModels.jl.

This should allow you to run! a Simulation with TimeStepWizard. Perhaps it makes sense to add a test for this in this PR?

If this does not work, feel free to copy/paste the error here so that we can debug the issue.

Note: I was confused earlier --- I definitely prefer shallow_water_cell_advection_timescale.jl for the filename. But anything will do.

[francispoulin]

Thanks @glwagner I'm happy to say this worked! Sorry that I was slow in putting the pieces together.

I agree that a test would no be a bad idea. I glanced in test_time_stepping.jl and I did not see anything obvious that tested the time stepping. I modfiied test_shallow_water.jl to use the wizard but not it uses it for everything. Not necessarily bad, but not sure if this is what we want? Not sure if we want to have two of these functions.

If there are not any tests for the wizard in IncompresibleModel maybe something should be done there and ShallowWaterModel can parallel that?

function time_stepping_shallow_water_model_works(arch, topo, coriolis)
    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(2π, 2π, 2π), topology=topo)
    model = ShallowWaterModel(grid=grid, gravitational_acceleration=1, architecture=arch, coriolis=coriolis)
    set!(model, h=1)

    wizard = TimeStepWizard(cfl=1.0, Δt=1.0, max_change=1.1, max_Δt=10)

    simulation = Simulation(model, Δt=wizard, stop_iteration=1)
    run!(simulation)

    return model.clock.iteration == 1
end

[glwagner]

This looks like a great simple test. Perhaps time_step_wizard_shallow_water_model_works is a good name for the function? You may want stop_iteration=2 to ensure that the update_dt works as expected.

[francispoulin]

Done! There is a test for the time stepping wizard that loops over archs, but not topologies. I can certainly test those as well but didn't think that was necessary. Happy to fix it if there is any concern.

If there aren't then I'm happy to close this issue and then focus on the accuracy of the solution, and comparing the results with the linear stability results that I posed yesterday.

[glwagner]

Is there a PR that resolves this issue?

[francispoulin]

Not quite. It is part of the Bickey Jet example. I am sorry that I didn't make a separrate PR, that would have been better. If you want I can make two separate PR's, but I suspect that's going to be more trouble.

[glwagner]

It's ok --- it just means that we can't easily review your changes.

You want to link the PR to this issue rather than closing this issue manually. To do this write "Resolves #1307" in the comments on that PR.

[francispoulin]

Thanks @glwagner , and I just did that. If I need to do something else to close this, just let me know.

[glwagner]

The issue should be closed automatically when #1305 is merged!

[francispoulin]

Great, and glad that github is that smart!

[ali-ramadhan]

It wasn't lol so I'll close it manually now.