error from twoband.jl with high vertical resolution

[johnryantaylor]

I am getting an error from twoband.jl which seems to occur when I run with high (<1m) vertical resolution. When I run eady.jl and use 128 gridpoints in the vertical, I get the error message below. I am using Julia 1.9.2 on my macbook pro. I also get a very similar error message in a high resolution version of a column model (where I first encountered this error).

[ Info: Initializing simulation...
i: 0, sim time: 0 seconds, wall time: 0 seconds, Δt: 16.500 minutes, CFL: 1.84e-01
[ Info: ... simulation initialization complete (9.006 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (54.796 seconds).
ERROR: DomainError with -2.618155629978065e21:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(x::Float64)
@ Base.Math ./math.jl:37
[2] ^(x::Float64, y::Float64)
@ Base.Math ./math.jl:1123
[3] macro expansion
@ ~/.julia/packages/OceanBioME/mii75/src/Light/2band.jl:30 [inlined]
[4] macro expansion
@ ~/.julia/packages/KernelAbstractions/cWlFz/src/extras/loopinfo.jl:26 [inlined]
[5] macro expansion
@ ~/.julia/packages/OceanBioME/mii75/src/Light/2band.jl:28 [inlined]
[6] cpu_update_TwoBandPhotosyntheticallyActiveRadiation!
@ ~/.julia/packages/KernelAbstractions/cWlFz/src/macros.jl:276 [inlined]
[7] cpu_update_TwoBandPhotosyntheticallyActiveRadiation!(ctx::Ker

[johnryantaylor]

This appears to be due to P becoming negative since I can prevent it by putting P inside an abs() in twoband.jl. I'm not sure why P is negative when the resolution is high and not otherwise (maybe something to do with the halo?).

[johnryantaylor]

I would be interested to know if anyone else gets the same error by running the eady.jl example with 128 gridpoints in the vertical direction.

[jagoosw]

This error usually occurs because of negative P. Perhaps we should try and catch it somehow and make it produce a better error message.

For the Eady example, I also get this error with the default timestep, but if I reduce the initial timestep to 40s and make the CFL conditions both 0.5 it runs fine (and settles at about Δt = 2 minutes). I think the negative P is just from the CFL condition, but I also think from this and other things that we may have an issue with the negative tracer scaling.

I think that here P becomes too negative for it to scale the tracers, the callback throws an error but it does not get raised properly (like the issue with Ki's code the other day), and then the negative P does not get resolved before it tries to integrate the light. I've also had problems with TendencyCallsite callbacks not correctly modifying tendencies in Oceananigans but haven't worked out what's going on with that yet, so this might be it.

I have found that callbacks that do this can be fixed by moving them to update_tendencies!, so when I've finished #134 I think I might be able to tidy the negativity protection up a lot and move it to there so it works properly. I've also found

[johnryantaylor]

Yes, you're right that it seems to be down to the timestep, thanks! Actually, it seems to be fixed by simply reducing the size of the first timestep and leaving the CFL condition unchanged. I'll reduce the size of the initial timestep in the Eady example - since the CFL condition isn't checked for the first timestep, its better to start out more conservative and let the timestep increase from there.

[johnryantaylor]

I am going to close this issue, but I agree with @jagoosw that we should try to make the negative value prevention more robust (particular with respect to the light functions).