Help with shallow water initial height conditions and NaN error.

[rbenamran]

Hi,
I've been trying to learn to use the ShallowWaterModel and tried to set a basic example with a centered gaussian height field for the initial conditions. I'm having trouble understanding the results I get. I've changed and tested many different parameters in the grid and model, timestepping and decreasing the initial height condition. However, it always seems that the height field falls normally at first, then rises quickly until the simulation aborts due to a NaN error in for the fields (perhaps the numbers become too large?).

Here is the animation :

Simple.Gaussian.Height3.mp4

and the code:

using Oceananigans
using Oceananigans.Models: ShallowWaterModel

Nx, Ny = 16, 16
Lx, Ly = 16, 16

grid = RegularRectilinearGrid(size = (Nx, Ny),
                            x=(0, Lx), y=(0, Ly),
                            topology = (Periodic, Periodic, Flat))

model = ShallowWaterModel(
    architecture=CPU(),
    grid=grid,
    advection=nothing,
    gravitational_acceleration=10,
    coriolis=nothing,
    )

x₀, y₀ = Lx/2, Ly/2

h₀(x,y,z) = 0.2 * exp( - ((x - x₀)^2 + (y - y₀)^2)/2 )

set!(model, h=h₀)


uh, vh, h = model.solution
simulation = Simulation(model, Δt=0.1, stop_iteration=60)

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:fields] =
    JLD2OutputWriter(
        model,
        model.solution,
        prefix = "Simple Gaussian Height3",
        schedule=IterationInterval(1),
        force = true)

run!(simulation)

using JLD2

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

xh, yh, zh = nodes(h)

using Plots

@info "Making a neat movie of height..."

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    h_snapshot = file["timeseries/h/$iteration"][:, :, 1]

    h_lim = 2.0
    h_levels = range(-h_lim, stop=h_lim, length=20)

    kwargs = (xlabel="x", ylabel="y", aspectratio=1, linewidth=0, colorbar=true,
        xlims=(0, model.grid.Lx), ylims=(0, model.grid.Ly))

    h_plot = contourf(xh, yh, clamp.(h_snapshot', -h_lim, h_lim);
                             color = :balance,
                            levels = h_levels,
                             clims = (-h_lim, h_lim),
                            kwargs...)

    plot(h_plot, title="Height", layout=(1), size=(600, 500))
end
mp4(anim, "Simple Gaussian Height3.mp4", fps = 8) # hide

Thanks in advance for any help understanding this!

[francispoulin]

If you don't have a Coriolis force initially, then you will always get that an initial bump in height will fall. If you use a coriolis force and choose the velocity in geostrophic balance, then you can get a vortex, that will be a steady solution to the nonlinear steady equations.

As for the NaN, have you tried to reduce the time step to see if the problem is with the CFL condition, i.e. time step is too big?

[glwagner]

@rbenamran I think @francispoulin is suggesting to use something like

H = 10 # unperturbed layer depth
h₀(x, y, z) = H + 0.2 * exp( - ((x - x₀)^2 + (y - y₀)^2)/2 )

for the initial condition. In other words h is the total layer depth, not just the perturbation around some mean / constant value.

[rbenamran]

Understood, it indeed works now by both adding the unperturbed value as well as adapting the time step to the model conditions. I will keep in mind these factors in the future to avoid the numerical issues.

Thanks a lot to both of you for your help as well as for the development of Oceananigans, I'm really enjoying using it!

(Here is one of the results I got, for info, with Δx=0.25, Δt=0.01, H=10):

Simple.Gaussian.Height4.mp4

[francispoulin]

Very happy that it worked out for you and that you are enjoying it. Lots of fun can be had wtih a shallow water model!

[williamtiberio]

yeah

[williamtiberio]

this is fine