DomainError when changing grid size in simple column model

[andrewwatford]

I am interested if anyone is able to point me in the right direction with an error I am encountering. I am trying to modify the simple column example here to have a different grid (instead of 3D, it is 1D, and the number of points is greater). The error is ERROR: DomainError with -0.000840929789711402: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. I think it has something to do with the two-band light attenuation, but not exactly sure what is going wrong. Any help would be appreciated, thanks! Here's a minimal example that reproduces the error:

using Oceananigans.Fields: FunctionField, ConstantField
using Oceananigans.Units

const year = years = 365days;

Nz = 100
Lz = 200meters
grid = RectilinearGrid(size = Nz, extent = Lz, topology = (Flat, Flat, Bounded));

@inline H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)

@inline fmld1(t) = H(t, 50days, year) * (1 / (1 + exp(-(t - 100days) / 5days))) * (1 / (1 + exp((t - 330days) / 25days)))

@inline MLD(t) = - (10 + 340 * (1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) - fmld1(mod(t, year))))

@inline κₜ(z, t) = 1e-2 * (1 + tanh((z - MLD(t)) / 10)) / 2 + 1e-4

biogeochemistry = LOBSTER(; grid)
initial_conditions = (; P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05);

clock = Clock(; time = 0.0) 

model = NonhydrostaticModel(; grid,
                            clock = clock,
                            closure = ScalarDiffusivity(ν = κₜ, κ = κₜ),
                            biogeochemistry = biogeochemistry)

set!(model; initial_conditions...)

simulation = Simulation(model, Δt = 3minutes, stop_time = 5days)
run!(simulation)

[francispoulin]

Thanks in advance to anyone who can help us understand this problem!

[jagoosw]

Hi @andrewwatford and @francispoulin,

Have you tried making the timestep smaller? Or maybe setting scale_negatives=true in LOBSTER. This error usually comes up when numerical error has made the phytoplankton/chlorophyll go negative. The diffusive CFL is only around 0.45 now the grid is finer which is probably too high for AB2.

If this doesn't work I can look into it more tomorrow!

[francispoulin]

Thanks! When I reduced the time step I did get a stable result. I am sorry that we didn't try this before posting the code. I had forgotten that the fields can become negative, and that would cause problems.

[jagoosw]

No problem!

[jagoosw]

We should really try to tell the user when this is happening but you can't throw warnings/errors in the normal way on the GPU so it might take some work

[ali-ramadhan]

We should really try to tell the user when this is happening but you can't throw warnings/errors in the normal way on the GPU so it might take some work

One easy solution would be to add a simulation callback to check for this and either print a warning or throw an error to stop the simulation. Similar to how Oceananigans adds a NaNChecker.

[jagoosw]

I guess this would work as an UpdateStateCallsite callback because that happens before the update_biogeochemical_state!

[jagoosw]

Also, we have the "negativity protection" things as modifiers in OceanBioME Biogeochemistry which would also stop this.

https://github.com/OceanBioME/OceanBioME.jl/blob/main/src/Utils/negative_tracers.jl

OceanBioME.jl/src/Models/AdvectedPopulations/NPZD.jl

Line 202 in 54cfd8a

scaler = ScaleNegativeTracers(underlying_biogeochemistry, grid; invalid_fill_value)