The “NaN error" problem

[zwei961120]

Hi all,

Conducting LES by using Oceananigans is very enjoyable and efficient, but recently I usually meet the NaN error problem and I can not figure out what happens. The case I try to simulate is about the evolution of surface and bottom boundary layers under the conditions of surface heat loss, geostrophic balance, and tidal forcing. At first, I just consider the heat loss and geostrophic balance, and everything goes well, but when I try to add the tidal forcing, there is always a NaN error when the simulation time arrives at about two days. I have checked out the discussions and found that @glwagner put forward possible reasons for the NaN error problem in discussion #1436. I have tried to change the time step, but it seems doesn't work.

Could you please help me figure out the reasons? Thanks in advance! My script is attached here :-D

using Random
using Printf
using JLD2
using Oceananigans
using Oceananigans:fields
using Oceananigans.Grids: on_architecture
using Oceananigans.ImmersedBoundaries: GridFittedBottom, GridFittedBoundary
using Oceananigans.Architectures: device
using Oceananigans.Units: minute, minutes, hour, hours

const H = 45
grid = RectilinearGrid(GPU();
                       size = (128, 128, 72),
                       halo = (3, 3, 3),
                          x = (0, 2π*H),
                          y = (0, 2π*H),
                          z = (-H, 0))
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4))

u₁₀ = 10    # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 1.25e-3 # dimensionless drag coefficient
ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean
Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

Qʰ = 200.0  # W m⁻², surface _heat_ flux
cᴾ = 3991.0 # J K⁻¹ kg⁻¹, typical heat capacity for seawater

Qᵀ = Qʰ / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ))

z₀ = 0.01 # m (roughness length)
κ = 0.4 # von Karman constant
z₁ = grid.Δzᵃᵃᶠ/2 # half the grid cell thickness
cᴰ = (κ / log(z₁ / z₀))^2 # Drag coefficient

 @inline drag_u(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + v^2) * u
 @inline drag_v(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + v^2) * v

 drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰ))
 drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰ))

 u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom = drag_bc_u)
 v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

f₀ = 1e-4
coriolis = FPlane(f=f₀)

U∞ = 0.25 # m s⁻¹

using Oceananigans.BuoyancyModels: g_Earth

const amplitude = 0.8 # m
const wavelength = 60  # m
const wavenumber = 2π / wavelength # m⁻¹
const frequency = sqrt(g_Earth * wavenumber * tanh(wavenumber * H)) # s⁻¹

const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹

∂z_uˢ(z, t) = Uˢ * wavenumber * sinh(2 * wavenumber * (z + H)) / (sinh(wavenumber * H))^2


const pgrad = f₀ * U∞
@inline pressure_gradient(x, y, z, t) = pgrad
v_forcing = Forcing(pressure_gradient)

const U₀ = 0.015 # m s⁻¹
const ω₂ = 1.4e-4 # rad s⁻¹, corresponding to 12.4 hours M2 tidal period
@inline pressure_gradient(x, y, z, t) = - U₀ * ω₂ * cos(ω₂ * t)
u_forcing = Forcing(pressure_gradient)

model = NonhydrostaticModel(; grid,
                advection = WENOFifthOrder(),
              timestepper = :RungeKutta3,
                  tracers = (:T, :S),
                  coriolis = coriolis,
                 buoyancy = buoyancy,
                  closure = AnisotropicMinimumDissipation(),
             stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
              boundary_conditions = (;T=T_bcs, v=v_bcs, u=u_bcs),
                  forcing = (;u=u_forcing, v=v_forcing))

S₀ = 35.0
T₀ = 20.0
dTdz = 0.1 # K m⁻¹
T0(z) = T₀ + dTdz * z

Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

Tᵢ(x, y, z) = T0(z) + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

uᵢ(x, y, z) = U∞ + 1e-5 * Ξ(z)

wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)


set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=S₀)

simulation = Simulation(model, Δt=5.0, stop_time=240hours)

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s, advective CFL: %.2e, diffusive CFL: %.2e\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w),
                                prettytime(sim.run_wall_time),
                                AdvectiveCFL(sim.Δt)(sim.model),
                                DiffusiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ)

simulation.output_writers[:slices] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity),
                         filename = "sliceslrtideheatloss.jld2",
                          indices = (:, grid.Ny/2, :),
                         schedule = TimeInterval(10minutes),
               overwrite_existing = true)

fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivity_fields.νₑ,))
simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity),
                     filename = "fieldslrtideheatloss.jld2",
                     schedule = TimeInterval(10minutes),
           overwrite_existing = true)

u, v, w = model.velocities
U = Average(u, dims=(1, 2))
V = Average(v, dims=(1, 2))
T = Average(model.tracers.T, dims=(1, 2))
wt = Average(w * model.tracers.T, dims=(1, 2))
wu = Average(w * u, dims=(1, 2))
wv = Average(w * v, dims=(1, 2))
ww = Average(w * w, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2OutputWriter(model, (u=U, v=V, T=T, wt=wt, wu=wu, wv=wv, ww=ww),
                     schedule = TimeInterval(10minutes),
                     filename = "averageslrtideheatloss.jld2",
           overwrite_existing = true)

run!(simulation)

The Hovmöller Diagram of XY-Plane-Averaged-U-velocity before the time simulation stopped is shown below:

And the error information is:

Iteration: 25900, time: 1.781 days, Δt: 6.972 seconds, max(|w|) = 6.1e-02 ms⁻¹, wall time: 1.153 hours, advective CFL: 1.00e+00, diffusive CFL: 1.37e-01
Iteration: 26000, time: 1.789 days, Δt: 7.140 seconds, max(|w|) = 6.2e-02 ms⁻¹, wall time: 1.158 hours, advective CFL: 1.00e+00, diffusive CFL: 1.52e-01
Iteration: 26100, time: 1.797 days, Δt: 7.150 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 1.162 hours, advective CFL: 1.00e+00, diffusive CFL: 1.93e-01
Iteration: 26200, time: 1.805 days, Δt: 7.123 seconds, max(|w|) = 6.7e-02 ms⁻¹, wall time: 1.166 hours, advective CFL: 1.00e+00, diffusive CFL: 1.63e-01
Iteration: 26300, time: 1.813 days, Δt: 6.937 seconds, max(|w|) = 8.0e-02 ms⁻¹, wall time: 1.170 hours, advective CFL: 1.00e+00, diffusive CFL: 1.78e-01
Iteration: 26400, time: 1.821 days, Δt: 6.752 seconds, max(|w|) = 8.8e-02 ms⁻¹, wall time: 1.175 hours, advective CFL: 1.00e+00, diffusive CFL: 1.66e-01
Iteration: 26500, time: 1.829 days, Δt: 6.825 seconds, max(|w|) = 8.8e-02 ms⁻¹, wall time: 1.179 hours, advective CFL: 1.00e+00, diffusive CFL: 1.78e-01
Iteration: 26600, time: 1.837 days, Δt: 6.878 seconds, max(|w|) = 9.0e-02 ms⁻¹, wall time: 1.183 hours, advective CFL: 1.00e+00, diffusive CFL: 2.17e-01
Iteration: 26700, time: 1.845 days, Δt: 6.762 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 1.187 hours, advective CFL: 1.00e+00, diffusive CFL: 2.11e-01
Iteration: 26800, time: 1.853 days, Δt: 6.340 seconds, max(|w|) = 9.9e-02 ms⁻¹, wall time: 1.192 hours, advective CFL: 1.00e+00, diffusive CFL: 1.91e-01
Iteration: 26900, time: 1.860 days, Δt: 6.971 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 1.196 hours, advective CFL: 1.00e+00, diffusive CFL: 2.30e-01
Iteration: 27000, time: 1.868 days, Δt: 7.114 seconds, max(|w|) = 7.6e-02 ms⁻¹, wall time: 1.200 hours, advective CFL: 1.00e+00, diffusive CFL: 2.09e-01
Iteration: 27100, time: 1.877 days, Δt: 7.216 seconds, max(|w|) = 8.7e-02 ms⁻¹, wall time: 1.205 hours, advective CFL: 1.00e+00, diffusive CFL: 2.69e-01
Iteration: 27200, time: 1.885 days, Δt: 6.969 seconds, max(|w|) = 9.0e-02 ms⁻¹, wall time: 1.209 hours, advective CFL: 1.00e+00, diffusive CFL: 1.86e-01
Iteration: 27300, time: 1.893 days, Δt: 7.518 seconds, max(|w|) = 8.3e-02 ms⁻¹, wall time: 1.213 hours, advective CFL: 1.00e+00, diffusive CFL: 2.21e-01
Iteration: 27400, time: 1.903 days, Δt: 8.321 seconds, max(|w|) = 6.5e-02 ms⁻¹, wall time: 1.218 hours, advective CFL: 1.00e+00, diffusive CFL: 2.36e-01
[ Info: time = 165029.77834255935, iteration = 27500: **NaN found in field u. Stopping simulation.**
Iteration: 27500, time: 1.910 days, Δt: 8.189 seconds, max(|w|) = 0.0e+00 ms⁻¹, wall time: 1.224 hours, advective CFL: 0.00e+00, diffusive CFL: 0.00e+00

[glwagner]

@zwei961120 I would try reducing the cfl in the TimeStepWizard to cfl=0.5. It looks like your Δt is increasing just prior to blow up:

Iteration: 27200, time: 1.885 days, Δt: 6.969 seconds
Iteration: 27300, time: 1.893 days, Δt: 7.518 seconds
Iteration: 27400, time: 1.903 days, Δt: 8.321 seconds

A cfl=1.0 is at the very edge of the "practical CFL" for RK3. I'd recommend going no higher than cfl=0.8. But cfl=0.5 (or even just cfl=0.1, for sanity's sake) should provide a conservative buffer.

Also, to diagnose the problem I think you should look at slices of w, rather than the horizontal average or u (which is rather smoothed).

Other than that, nothing else pops out when I read your script. The setup seems sound!

[zwei961120]

Hi @glwagner, thanks for your kind reply and suggestions, I tried to reduce the cfl number to 0.7 and everything goes well!