Repository navigation
How to determine the maximum stable time step for biogeochemistry dynamics? #237
Description
I was just playing around a LOBSTER box model (basically the first example in the docs) and noticed that I was getting some super negative values. Turns out it's probably numerical instability due to the time step because reducing the time step eventually leads to convergence (although sometimes blowup at intermediate time step sizes).
I'm curious if there's a way to determine the maximum stable time step for biogeochemistry dynamics (e.g. LOBSTER). This way we could take it into account when computing time step sizes for adaptive time-stepping.
If not, maybe the solution is an equivalent of the NaNChecker that stops a simulation if a tracer becomes too negative?
cc @suki-atdepthmrv @almacarolina
MWE:
using CairoMakie
using Oceananigans
using Oceananigans.Units
using Oceananigans.Fields: FunctionField
using OceanBioME
const year = years = 365day
grid = BoxModelGrid()
clock = Clock(time = zero(grid))
PAR⁰(t) = 100 * (1 - cos((t + 150days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
const z = -10
PAR_func(t) = PAR⁰(t) * exp(0.2z)
PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)
biogeochemistry = LOBSTER(;
grid,
light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation(PAR),
carbonates = true,
oxygen = true,
variable_redfield = true,
scale_negatives = false
)
model = BoxModel(; biogeochemistry, clock)
set!(model,
NO₃ = 50,
DIC = 2300,
Alk = 2400,
O₂ = 250,
NH₄ = 0.1,
DOC = 100,
DON = 10,
sPOC = 10,
bPOC = 1,
sPON = 1,
bPON = 1,
P = 5,
Z = 2
)
simulation = Simulation(model; Δt = 10seconds, stop_time = 30days)
simulation.output_writers[:fields] =
JLD2OutputWriter(model, model.fields;
filename = "box.jld2",
schedule = TimeInterval(1hour),
overwrite_existing = true
)
prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(simulation.run_wall_time))"
simulation.callbacks[:progress] = Callback(prog, TimeInterval(1day))
@info "Running the model..."
run!(simulation)
times = FieldTimeSeries("box.jld2", "P").times
time_series = NamedTuple{keys(model.fields)}(FieldTimeSeries("box.jld2", "$field")[1, 1, 1, :] for field in keys(model.fields))
fig = Figure(size = (1200, 1200), fontsize = 24)
axs = []
for (name, tracer) in pairs(time_series)
idx = (length(axs))
push!(axs, Axis(fig[floor(Int, idx/2), Int(idx%2)], ylabel = "$name", xlabel = "Time (days)"))
lines!(axs[end], times / day, tracer, linewidth = 3)
end
save("box_model.png", fig)With a time step of 15 mins (super negative NO3, Alk, and DIC):
With a time step of 5 mins (crash and NaNs after ~0.5 days):
With a time step of 1 min (seems fine):
With a time step of 10 seconds (seems converged):
Environment:
OceanBioME.jl main branch with
Julia Version 1.10.7
Commit 4976d05258e (2024-11-26 15:57 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 48 × AMD Ryzen Threadripper 7960X 24-Cores
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 16 default, 0 interactive, 8 GC (on 48 virtual cores)
Environment:
LD_PRELOAD = /usr/NX/lib/libnxegl.so
Activity
jagoosw commentedon Jan 20, 2025
Hi Ali,
Sorry, I completely missed this before.
As far as I'm aware it's not generally possible to determine a stable timestep and it would be very involved to compute it.
This is why I came up with the
ScaleNegativeTracersmodifier, but I think in a box model it makes sense to error on negatives. Maybe it would be good to have a callback for that in the source code.