Gravity Wave Speed Expectation

[rbenamran]

Hi,

This question may be due to a misunderstanding of mine, but I have been using the same code from #2672, that is to say a shallow water model sim with a random height field initial condition, and have been seeing slower speeds than expected.

I turned off the Coriolis to test this, as I then expected to observe gravity wave speeds of sqrt(gH). But, for example with the precise parameters from #2672, the gravity waves should be around 31m/s while the max speeds I obtain throughout the iterations are only between 1 and 2 m/s.

Is there something I’m missing to justify the results? I've kept in kind geostropic balance when setting different parameters but I am still always about an order of magnitude less. I've also increased time steps so duration shouldn't be the cause. This is important as I’m hoping to do some spectral analysis on Poincare waves further on (with Coriolis turned back on of course). Thanks

[glwagner]

There's many possibilities: nonlinearity, dissipation, a bug in the script, or a bug in the source code. Can you provide code to reproduce your result?

I would start by looking at an initial value problem in a one-dimensional model (eg propagate a very weak and therefore linear Gaussian disturbance) to see if the wave speed is wrong there. If is, we have our work cut out for us. If it's not, we can go to 2D.

[rbenamran]

Hi, so sorry for the long absence! I've been coming back to this project and reran the code with the new versions of Julia and Oceananigans, and still found the wave speed discrepancy with my expectations. However, I think I may be analyzing the JDL2 data wrong somehow.

Here are some tests I ran, using the logistic function to try to make the video clearer, but this worked with the Gaussian as well. For these tests, in 1D or 2D, so far as I have seen, the gravity-wave speed is reproduced correctly (set up so that the waves reach the boundary at the final simulation time).

FermiDiracX.Wavefront.mp4

Gaussian.Wavefront.mp4

My trouble right now is that, when taking these sims for which I can see the wave speed, I still find wrong values when analyzing the JLD2 files. I tried both the Conservative and VectorInvariant formulation to check if division by h was wrong. The simple code for the Gaussian is below, could you let me know why I can't find a value around 10m/s directly from the model.solution arrays or the output writers?

Thanks for your input!

using Oceananigans
using Oceananigans.Models: ShallowWaterModel

Nx, Ny = 64, 64   #grid nodes
Lx, Ly, Lz = 60, 60, 10  #grid size

g = 9.81 #gravitational_acceleration       

c = sqrt(g * Lz) #surface gravity wave speed

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

model = ShallowWaterModel(
    grid = grid,
    gravitational_acceleration = g,
    coriolis = nothing,
    formulation = ConservativeFormulation(),
    # momentum_advection = VectorInvariant() #needed for other formulation
)

#Surface layer η
#wavelength to be long compared to lattice spacing and small compared to system.
A = 0.2   #amplitude
k = 1   #wavenumbers tested for cos condition

x₀ = Lx/2
y₀ = Ly/2
#Initial condition
h₀(x, y, z) = Lz + A * exp( - ((x - x₀)^2 + (y - y₀)^2)/4 )
# h₀(x, y, z) = Lz + A / (exp(5 * (sqrt((x-x₀)^2 + (y-y₀)^2) - 4)) + 1)

set!(model, h=h₀)

uh, vh, h = model.solution
u = uh/h

#CFL condition plus matching the coriolis parameter for frequency dispersion
simulation = Simulation(model, Δt=0.001, stop_iteration=3000)


#Plot.
using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:fields] =
    JLD2OutputWriter(
        model,
        model.solution,
        filename = "test.jld2",
        schedule=IterationInterval(1)
)

run!(simulation)

using JLD2, Printf

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_base = Lz
    h_lim = 0.4
    h_levels = range(h_base - h_lim, stop = h_base + h_lim, length=80)

    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_base - h_lim, h_base + h_lim);
                             color = :balance,
                            levels = h_levels,
                             clims = (h_base - h_lim, h_base + h_lim),
                            kwargs...)

    plot(h_plot, layout=(1), size=(600, 500))
    title!(@sprintf("Height at t = %.2f, H=%.0f, c=%.2f, maxu = %.2f", t, Lz, c, maximum(u)))
end
mp4(anim, "test.mp4", fps = 8) # hide