Blowup when trying to reproduce MITgcm Barotropic Ocean Gyre

[hdrake]

I am trying to reproduce the barotropic gyre verification experiment in the MITgcm: https://mitgcm.readthedocs.io/en/latest/examples/barotropic_gyre/barotropic_gyre.html

using Oceananigans
using Printf

Nx = 60
Ny = 60

Lx = 1200kilometers # east-west extent [m]
Ly = 1200kilometers # north-south extent [m]
Lz = 5kilometers    # depth [m]

grid = RectilinearGrid(size = (Nx, Ny, 1),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = (-Lz, 0),
                       topology = (Bounded, Bounded, Bounded))

coriolis = BetaPlane(f₀=1e-4, β=1e-11)

# Forcing by zonal wind stress
@show surface_wind_stress_parameters = (τ₀ = 0.1, Lφ = grid.Ly)
@inline surface_wind_stress(x, y, t, p) = -p.τ₀ * cos(π * y / p.Lφ)
surface_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress,
                                               parameters = surface_wind_stress_parameters)
u_bcs = FieldBoundaryConditions(top = surface_wind_stress_bc)

# Friction due to constant horizontal viscosity
@show νₕ = 400
constant_horizontal_diffusivity = HorizontalScalarDiffusivity(ν = νₕ)

@show model = HydrostaticFreeSurfaceModel(grid = grid,
                                    coriolis = coriolis,
                                    boundary_conditions = (u=u_bcs,),
                                    closure = constant_horizontal_diffusivity)

simulation = Simulation(model, Δt = 20minutes, stop_time = 365days)
output_fields = merge(model.velocities, (η=model.free_surface.η,))
output_prefix = "barotropic_gyre"

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    @info @sprintf("Time: %s, iteration: %d, max(u): %.2e m s⁻¹, wall time: %s",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(sim.model.velocities.u),
                   prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(Progress(time_ns()), IterationInterval(20))


simulation.output_writers[:fields] = JLD2Writer(model, output_fields,
                                                schedule = TimeInterval(1hour),
                                                filename = output_prefix,
                                                overwrite_existing = true)

@show run!(simulation)

The output is:

surface_wind_stress_parameters = (τ₀ = 0.1, Lφ = grid.Ly) = (τ₀ = 0.1, Lφ = 1.2e6)
νₕ = 400 = 400
model = HydrostaticFreeSurfaceModel(grid = grid, coriolis = coriolis, boundary_conditions = (u = u_bcs,), closure = constant_horizontal_diffusivity) = HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 60×60×1 RectilinearGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×1 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=400.0)
├── buoyancy: Nothing
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme: 
│   └── momentum: Vector Invariant, Dimension-by-dimension reconstruction
└── coriolis: BetaPlane{Float64}
run!(simulation) = nothing
[ Info: Initializing simulation...
[ Info: Time: 0 seconds, iteration: 0, max(u): 0.00e+00 m s⁻¹, wall time: 73.100 ms
[ Info:     ... simulation initialization complete (54.822 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.005 ms).
[ Info: Time: 6.667 hours, iteration: 20, max(u): 2.88e-01 m s⁻¹, wall time: 40.664 ms
[ Info: Time: 13.333 hours, iteration: 40, max(u): 5.74e-01 m s⁻¹, wall time: 32.157 ms
[ Info: Time: 20 hours, iteration: 60, max(u): 8.62e-01 m s⁻¹, wall time: 29.913 ms
[ Info: Time: 1.111 days, iteration: 80, max(u): 1.15e+00 m s⁻¹, wall time: 31.254 ms
[ Info: Time: 1.389 days, iteration: 100, max(u): 1.44e+00 m s⁻¹, wall time: 30.519 ms
[ Info: Time: 1.667 days, iteration: 120, max(u): 1.74e+00 m s⁻¹, wall time: 93.118 ms
[ Info: Time: 1.944 days, iteration: 140, max(u): 2.06e+00 m s⁻¹, wall time: 29.967 ms
[ Info: Time: 2.222 days, iteration: 160, max(u): 2.39e+00 m s⁻¹, wall time: 30.089 ms
[ Info: Time: 2.500 days, iteration: 180, max(u): 2.75e+00 m s⁻¹, wall time: 28.808 ms
[ Info: Time: 2.778 days, iteration: 200, max(u): 3.16e+00 m s⁻¹, wall time: 30.130 ms
[ Info: Time: 3.056 days, iteration: 220, max(u): 3.61e+00 m s⁻¹, wall time: 32.050 ms
[ Info: Time: 3.333 days, iteration: 240, max(u): 4.32e+00 m s⁻¹, wall time: 30.466 ms
[ Info: Time: 3.611 days, iteration: 260, max(u): 5.14e+00 m s⁻¹, wall time: 32.652 ms
[ Info: Time: 3.889 days, iteration: 280, max(u): 5.82e+00 m s⁻¹, wall time: 31.829 ms
[ Info: Time: 4.167 days, iteration: 300, max(u): 6.40e+00 m s⁻¹, wall time: 31.153 ms
[ Info: Time: 4.444 days, iteration: 320, max(u): 6.90e+00 m s⁻¹, wall time: 35.594 ms
[ Info: Time: 4.722 days, iteration: 340, max(u): 7.35e+00 m s⁻¹, wall time: 38.299 ms
[ Info: Time: 5 days, iteration: 360, max(u): 7.86e+00 m s⁻¹, wall time: 42.543 ms
[ Info: Time: 5.278 days, iteration: 380, max(u): 8.32e+00 m s⁻¹, wall time: 33.223 ms
[ Info: Time: 5.556 days, iteration: 400, max(u): 8.72e+00 m s⁻¹, wall time: 32.344 ms
[ Info: Time: 5.833 days, iteration: 420, max(u): 9.07e+00 m s⁻¹, wall time: 31.446 ms
[ Info: Time: 6.111 days, iteration: 440, max(u): 9.39e+00 m s⁻¹, wall time: 34.026 ms
[ Info: Time: 6.389 days, iteration: 460, max(u): 9.68e+00 m s⁻¹, wall time: 38.273 ms
[ Info: Time: 6.667 days, iteration: 480, max(u): 9.95e+00 m s⁻¹, wall time: 34.761 ms
[ Info: Time: 6.944 days, iteration: 500, max(u): 1.02e+01 m s⁻¹, wall time: 37.534 ms
[ Info: Time: 7.222 days, iteration: 520, max(u): 1.05e+01 m s⁻¹, wall time: 34.562 ms
[ Info: Time: 7.500 days, iteration: 540, max(u): 1.08e+01 m s⁻¹, wall time: 30.470 ms
[ Info: Time: 7.778 days, iteration: 560, max(u): 1.11e+01 m s⁻¹, wall time: 33.772 ms
[ Info: Time: 8.056 days, iteration: 580, max(u): 1.14e+01 m s⁻¹, wall time: 34.482 ms
[ Info: Time: 8.333 days, iteration: 600, max(u): 1.20e+01 m s⁻¹, wall time: 30.379 ms
[ Info: Time: 8.611 days, iteration: 620, max(u): 1.70e+01 m s⁻¹, wall time: 33.998 ms
[ Info: Time: 8.889 days, iteration: 640, max(u): 3.50e+14 m s⁻¹, wall time: 33.991 ms
[ Info: Time: 9.167 days, iteration: 660, max(u): NaN m s⁻¹, wall time: 32.250 ms
[ Info: Time: 9.444 days, iteration: 680, max(u): NaN m s⁻¹, wall time: 109.758 ms
[ Info: time = 840000.0, iteration = 700: NaN found in field u. Stopping simulation.
[ Info: Time: 9.722 days, iteration: 700, max(u): NaN m s⁻¹, wall time: 34.036 ms

The simulations blows up after a few days due to a numerical instability related to grid-scale noise in the vertical velocity:

I've tried changing the advection scheme, timestepper, timestep, and viscosity, but the blow-up persisted in all of my experiments.

[glwagner]

Can you try

@inline surface_wind_stress(x, y, t, p) = p.τ₀ * cos(π * y / p.Lφ)

I believe our sign convention is opposite MITgcm. Not sure that would make the difference though.

[glwagner]

Oh!

You also need

@show surface_wind_stress_parameters = (τ₀ = 0.1 / 1000, Lφ = grid.Ly)

because the reference density in MITgcm is set to 1000. I think.

[hdrake]

My hero–thank you! It was just a difference in the convention for the units for $\tau$!

[glwagner]

Nice!

Do you mind if I convert this to a discussion? I want it to be more visible to folks looking for answers to simulation setups

[glwagner]

My hero–thank you! It was just a difference in the convention for the units for τ !

Also a comment, we don't even have a reference density at all for this simulation. We could introduce one, but as you can infer, it's only purpose in this kind of setup would be to change the units of user inputs. I personally think that is not a desirable design (the behavior of the simulation depends on hidden variables that may be irrelevant in many contexts, eg without wind stress), so we don't do that right now. It does have the downside of being different than what is typical in other modeling systems though.