Egregious errors in `NonhydrostaticModel` simulations on `ImmersedBoundaryGrid` when nonhydrostatic pressure is not separated about oceananigans.jl HOT 21 OPEN

@glwagner, am I correct that this warning is only needed when 1) hydrostatic_pressure_anomaly = nothing and 2) immersed boundaries are being used?

The pressure solver is incorrect even with hydrostatic_pressure_anomaly = CenterField(grid), just that apparently only the nonhydrostatic part is incorrect so the error is smaller than with hydrostatic_pressure_anomaly = nothing.
#3188 should solve this issue

from oceananigans.jl.

I've just been playing around with immersed boundaries and found that with hydrostatic_pressure_anomaly = nothing I wasn't able to maintain an ocean at rest. It's almost like buoyant fluid was erupting out of the immersed boundary. T and S being zero inside the immersed boundary are probably causing huge gradients to show up somewhere?

I haven't dug into why this issue is so bad in my case, but I definitely second @tomchor's suggestion to make hydrostatic_pressure_anomaly = CenterField(grid) the default!

hydrostatic_pressure_anomaly = nothing

hydrostatic_pressure_anomaly = CenterField(grid)

from oceananigans.jl.

@xkykai so you're saying the fix for the instability exists already and it's just a matter of implementing it in #3188? Is that right?

Yes, and it’s here at Yixiao-Zhang@c7983b8

This should work and fix the numerical issues, but has not been extensively validated yet.

from oceananigans.jl.

It's hard to imagine buoyancy is globally conserved in the first of Isba's movies above, although I guess a divergent velocity field could create a widespread tracer extrema like that while still technically conserving buoyancy...

We'll run some simulations with flux BCs and check conservation.

from oceananigans.jl.

It's almost like buoyant fluid was erupting out of the immersed boundary.

Wow, this looks like a much more egregious version of the spurious buoyancy sources that @ikeshwani was seeing! I think it probably has to specifically do with buoyancy gradients across the boundary. @ikeshwani's initial condition happened to be b=0, which is the same as the fill value inside the immersed boundary, so the gradients only show up over time as the system adjusts.

from oceananigans.jl.

Right and that's a good succinct description. We could even have two warnings: (i) it's wrong and (ii) it's catastrophically wrong.

#3188 doesn't solve the problem, it actually introduces a totally new solver. After #3188, someone doing simulations in complex domains will have the choice of using an iterative solver with adjustable tolerance / iterations, or the old "incorrect" solver (which will be cheaper, and it seems that some people have been able to get mileage out of it).

from oceananigans.jl.

@xkykai would you like help with #3188 ?

Yes I will work on it.

Would you like any help?

I will take a look at the state of it and see what is required to complete it, and will let you know when I need any help! Thanks a lot.

from oceananigans.jl.

Likely introduced by #3574.

@tomchor recommends changing the default setting of NonHydrostaticModel to hydrostatic_pressure_anomaly = CenterField(grid) (at least until issue has been addressed).

from oceananigans.jl.

Could the false tracer extrema be due to the advection scheme, which exhibits more errors in the case that the pressure is not separated? Note that the pressure solver is approximate in either case. But when the pressure is separated, I think the pressure solver may be more accurate.

I'm not sure there is any way that the pressure separation could directly impact the conservation of tracer, so I'm hypothesizing that the effect is indirect and occurs through the differences in the emergent dynamics.

from oceananigans.jl.

Sounds like a good idea to change the default for sure!

from oceananigans.jl.

I'm not sure there is any way that the pressure separation could directly impact the conservation of tracer, so I'm hypothesizing that the effect is indirect and occurs through the differences in the emergent dynamics.

Yes, I think the tracer errors must be arising due to differences in the emergent dynamics. Isn't there something about the errors in the pressure solver meaning that the velocity field is partially divergent and that the no-normal-flow boundary condition is not perfectly satisfied?

from oceananigans.jl.

@xkykai can answer more thoroughly I think but in general yes, because the pressure solver is only approximate, we are left with a choice between 1) a divergent velocity field or 2) leakage through the boundary.

I believe though that we have selected option (1), which means that tracer advection is also approximate but also that tracers are conserved (unlike option 2, where tracers would leak through the boundary but, because the velocity is non-divergent as our advection formulation assumes, tracer advection is "accurate").

Might be good to check whether tracer is globally conserved.

from oceananigans.jl.

Another test that reveals a problem with non-separated pressure is a simple horizontal wall in a 2D setup.

from oceananigans.jl.

I'm not sure this is a bug by the way. The code is correct -- it may simply be that the pressure solver numerical method itself is wrong (or "approximate" if you want to be generous).

A pressure solver that is at least theoretically correct for complex domains is/was being developed here: #3188

from oceananigans.jl.

Indeed the choice of divergent boundary flow has been made so that there is no tracer leakage through the boundary. One can try using #3188 which satisfies both no divergent flow and tracer leakage conditions, but is slower than using the FFT solver itself. Numerical instability has been observed in this solver occasionally due to numerical errors of the conjugate gradient solver occasionally giving zero eigenvalues which creates the instability. I was trying to find the fix for it but I think it's in @Yixiao-Zhang 's fork of Oceananigans.

from oceananigans.jl.

@xkykai so you're saying the fix for the instability exists already and it's just a matter of implementing it in #3188? Is that right?

from oceananigans.jl.

@xkykai would you like help with #3188 ?

It seems that we should also change the default treatment of hydrostatic pressure. I also think it would be appropriate to add a warning remarking about the approximate nature of the current pressure solver (at least until #3188 is merged). Would anyone like to take the lead on this?

from oceananigans.jl.

I also think it would be appropriate to add a warning remarking about the approximate nature of the current pressure solver (at least until #3188 is merged).

@glwagner, am I correct that this warning is only needed when 1) hydrostatic_pressure_anomaly = nothing and 2) immersed boundaries are being used?

from oceananigans.jl.

@xkykai would you like help with #3188 ?

Yes I will work on it.

from oceananigans.jl.

@xkykai would you like help with #3188 ?

Yes I will work on it.

Would you like any help?

from oceananigans.jl.

I made a PR that adds a warning and changes the default for immersed boundary grids: #3692

I don't think this issue can really be closed. There will continue to be errors on immersed boundary grid for non-separated pressure.

Perhaps, once we have a new solver that we are confident reduces the chance of egregious errors, we can convert this to a discussion.

from oceananigans.jl.