Scaling ShallowWaterModel causes static issue

[rbenamran]

Hi,
I have been working with a shallow water model simulation that consists simply of a stochastic height field as the initial condition. It's been working well on a 2pi by 2pi domain without using the Units module.

However, once I tried to scale this same simulation to a 2000km by 2000km domain, I encountered some issues. Specifically, it seems the simulation becomes static. The output I'm viewing is just the height field animation. I also apologize if the initial condition code is confusing and added an image of the initial height function.

I scaled the initial height function to match the domain so that should not be a source of bugs, theoretically. I've been testing and trying to fix this issue but am not sure what to do.

Thanks for any help or advice!

Code:

using Oceananigans
using Oceananigans.Models: ShallowWaterModel
using Oceananigans.Units

#Parameters
Nx, Ny = 128, 128
Lx = 1000kilometer
Ly = 1000kilometer

grid = RectilinearGrid(size = (Nx, Ny),
                            x=(-Lx, Lx), y=(-Lx, Ly),
                            topology = (Periodic, Periodic, Flat))

model = ShallowWaterModel(
    grid=grid,
    gravitational_acceleration=9.81,
    coriolis=FPlane(f=1),
    )

#SMOOTH RANDOM FUNCTION
using Random
Random.seed!(17)

n = 16
m = 16
len = n * m * 4

global alea
alea = rand(-1:1, (n, m, 4))

function rand_smth(x::Float64, y::Float64)
    X = x / 1000000
    Y = y / 1000000
    val = 0
    for i in 1:n
        for j in 1:m
            ar = alea[i,j,1]
            br = alea[i,j,2]
            cr = alea[i,j,3]
            dr = alea[i,j,4]
            val += ar * cos(i * x) * cos(j * y)
            val += br * cos(i * x) * sin(j * y)
            val += cr * sin(i * x) * cos(j * y)
            val += dr * sin(i * x) * sin(j * y)
        end
    end
    return val
end
#End of smooth random function

H = 1kilometer # unperturbed layer depth

h₀(x, y, t) = H + rand_smth(x,y)

set!(model, h=h₀)

uh, vh, h = model.solution

simulation = Simulation(model, Δt=0.002, stop_iteration=10000)


using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:fields] =
    JLD2OutputWriter(
        model,
        model.solution,
        filename = "simkm1.jld2",
        schedule=IterationInterval(4),
        overwrite_existing = 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 = 1.5meter
    h_levels = range(H - h_lim, stop = H + h_lim, length=40)

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

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

    plot(h_plot, title="Height", layout=(1), size=(600, 500))
end

mp4(anim, "simkm1.mp4", fps = 64) # hide

Initial random height function:

Output (small domain):

sim_reg.mp4

Output (large domain):

simkm1.mp4

[glwagner]

@rbenamran do you mind if I convert this to a discussion? This would fit in the "computational science" category: https://github.com/CliMA/Oceananigans.jl/discussions/categories/computational-science

The issue here is the function random_smth which uses wavenumbers i, j that range from 1 to 16:

function rand_smth(x::Float64, y::Float64)
    X = x / 1000000
    Y = y / 1000000
    val = 0
    for i in 1:n
        for j in 1:m
            ar = alea[i,j,1]
            br = alea[i,j,2]
            cr = alea[i,j,3]
            dr = alea[i,j,4]
            val += ar * cos(i * x) * cos(j * y)
            val += br * cos(i * x) * sin(j * y)
            val += cr * sin(i * x) * cos(j * y)
            val += dr * sin(i * x) * sin(j * y)
        end
    end
    return val
end

These represent long wavelengths on domains of size 2pi by 2pi, but very short wavelengths on a domain that's 1000 kilometers square. A quick change is to define random_smth as

function rand_smth(x::Float64, y::Float64)
    X = x / 1000000
    Y = y / 1000000
    val = 0
    for i in 1:n
        for j in 1:m
            
            # Dimensional wavenumbers
            k = 2π / Lx
            l = 2π / Lx

            ar = alea[i,j,1]
            br = alea[i,j,2]
            cr = alea[i,j,3]
            dr = alea[i,j,4]
            val += ar * cos(k * x) * cos(l * y)
            val += br * cos(k * x) * sin(l * y)
            val += cr * sin(k * x) * cos(l * y)
            val += dr * sin(k * x) * sin(l * y)
        end
    end
    return val
end

[simone-silvestri]

Hmmm, the velocities in your system are in the order of 0.005 m/s (maximum(uh) ~ 0.5).

simkm_9.mp4

This means that in 10000 s you move by roughly 50 meters. Since the cell is 7.8 km large, in that time span it is impossible to see a difference from the initial field. Are all the parameters correct?

If you look at geostrophic balance it kind of makes sense since u ~ ∂y(h)/f where h is relatively small compared to the size of your domain (which means very small pressure force) and f is quite large

[rbenamran]

Interesting, thanks for the analysis! That does make sense, I hadn't thought it would be so low. Changing the parameters to lower Coriolis and higher h definitely increases the velocity, though still a little slow. I also wanted to avoid too many nonlinear effects by increasing h, and need at least a decent Coriolis effect. Beyond these parameters, is there anything else I can change to increase the velocities within this large domain with realistic gravity without doing something like forcing?

Thanks again for the help in figuring out what the problem was!

[glwagner]

If you fix the horizontal and vertical size of the domain, Coriolis parameter, gravity, and the nonlinearity / amplitude of the initial condition, I believe the problem is completely defined, correct?

[rbenamran]

Yup, just checking to see if I missed anything!

[glwagner]

Also just to clarify --- the time-scale in this problem is set by the wave propagation time-scale, not the amplitude of the velocities, because the dynamics are nearly linear (or only weakly nonlinear). The phase and group speed is c = sqrt(g * h); so the time-scale for changes to appear on a particular length scale L is T = L / c. With L as the grid spacing, this also scales the time-step you have to use for numerical stability.

[francispoulin]

Also, I noticed that you have f=1, when you define the model. I wonder if you might want to use something like f=1e-4, to get closer to geophysical values? If yes, then you probably want to increase the time scale appropriately.

[glwagner]

Good point! The appropriate time step / model evolution time should also scale.