distributed_hydrostatic_turbulence.jl yields NaNs

[francispoulin]

I have tried running all the scripts in distributed_simulations and they all failed for me. I thought I would point out the problems I have found here and we could clean them up. First, let's start with this one.

distributed_hydrostatic_turbulence.jl

It starts fine but then I get NaN, which suggests to me that the time step is too large. I reduced the cfl parameter from 0.2 to 0.1 and instead of dying at iteration 200 it died at 6100. Better but not great.

I am going to try 0.05 but is it a concern that it ran several months ago with these parameters and now it doesn't? Also, why does the cfl need to be so small? I would think that a cfl of 0.2 should be great for pretty much any simulation.

[glwagner]

I don't think you can actually adapt the time step right now. Try with a fixed time step.

Also the setup looks weird to me:

    model = HydrostaticFreeSurfaceModel(; grid,
                                        momentum_advection = VectorInvariant(vorticity_scheme=WENO(order=9)),
                                        free_surface = SplitExplicitFreeSurface(grid, substeps=10),
                                        tracer_advection = WENO(),
                                        buoyancy = nothing,
                                        coriolis = FPlane(f = 1),
                                        tracers = :c)

This is WENO for vorticity but second order for everything else? And no other closure. Can you modify the physics so that we have a hope the simulation will run? You want to use WENOVectorInvariant. Also 10 substeps for split explicit is too few probably.

[francispoulin]

Very good idea! I just noticed that it was proceeding with a time step of 1e-84 before it produced an NaN.

[glwagner]

Might make sense to try it in serial and make sure the setup runs before trying to distribute it.

[francispoulin]

This is my serial version of the hydrostatic script. It fails for me after 100 iterations with NaN in the u field.

Maybe the suggestions that @glwagner made will fix this up and then we can go to the distributed version?

using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
using Printf
using Statistics
using Oceananigans.BoundaryConditions
using Random
using JLD2
using Oceananigans.ImmersedBoundaries: ActiveCellsIBG, active_interior_map

function run_simulation(nx, ny; topology = (Periodic, Periodic, Bounded))
    grid = RectilinearGrid(CPU(); topology, size = (Nx, Ny, 10), extent=(4π, 4π, 0.5), halo=(8, 8, 8))
    
    bottom(x, y) = (x > π && x < 3π/2 && y > π/2 && y < 3π/2) ? 1.0 : - grid.Lz - 1.0
    grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom); active_cells_map = true)

    model = HydrostaticFreeSurfaceModel(; grid,
                                        momentum_advection = VectorInvariant(vorticity_scheme=WENO(order=9)),
                                        free_surface = SplitExplicitFreeSurface(grid, substeps=10),
                                        tracer_advection = WENO(),
                                        buoyancy = nothing,
                                        coriolis = FPlane(f = 1),
                                        tracers = :c)

    Random.seed!(1234)

    set!(model, u = (x, y, z) -> 1-2rand(), v = (x, y, z) -> 1-2rand())
    
    mask(x, y, z) = x > 3π/2 && x < 5π/2 && y > 3π/2 && y < 5π/2
    c = model.tracers.c

    set!(c, mask)

    u, v, _ = model.velocities
    η = model.free_surface.η
    outputs = merge(model.velocities, model.tracers)

    progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time: $(sim.model.clock.time), Δt: $(sim.Δt)"
    simulation = Simulation(model, Δt=0.02, stop_time=100.0)

    simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))
 
    filepath = "hydrostatic_turbulence"
    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, outputs, filename=filepath, schedule=TimeInterval(0.1),
                         overwrite_existing=true)

    run!(simulation)
end

Nx, Ny = 32, 32

run_simulation(Nx, Ny)

[simone-silvestri]

Oh, this validation is a little old and not up to date. I ll open a PR to correct it.

[francispoulin]

Thank you @simone-silvestri !

[glwagner]

For reference the script is here: https://github.com/CliMA/Oceananigans.jl/blob/main/validation/distributed_simulations/distributed_hydrostatic_turbulence.jl