Oceananigans.jl
Oceananigans.jl copied to clipboard
CliMA
Egregious errors in `NonhydrostaticModel` simulations on `ImmersedBoundaryGrid` when nonhydrostatic pressure is not separated
I have noticed an issue with immersed boundaries in the latest version of Oceananigans. At the immersed boundaries we expect the default tracer flux boundary conditions to be zero. While doing energy analysis, I noticed that something seemed off with my simulation. After increasing the colorrange in my animation, I noticed that near my immersed boundaries, the minimum buoyancy was drifting to values much less than the minimum value prescribed by the surface value boundary conditions, suggesting that the tracer flux at the immersed boundary is non-zero.
For context, I'm running a 2D Horizontal Convection simulation, in which the buoyancy is initialized as zero everywhere to start. I apply a buoyancy gradient at the surface, cooling half of the surface and warming the other half. The boundaries are insulated and I have a pair of gaussian hills at the bottom defined using the immersed boundary function. Link to my simulation setup.
Here is an animation of topographically-constrained horizontal convection where the nonhydrostatic pressure is not separated:
https://github.com/user-attachments/assets/9fd054c9-f1da-4b8d-9d04-b99479dbe348
Notice there seems to be a source of dense fluid in the basin between the hills.
To address this issue @hdrake and I separated the hydrostatic and nonhydrostatic pressure components in the NonHydrostaticModel by changing https://github.com/CliMA/Oceananigans.jl/blob/406eb9c5c7a9fc86947747116128c8c1ba4c93d4/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl#L70.
To hydrostatic_pressure_anomaly = CenterField(grid)
Here is the animation of topographically-constrained horizontal convection when the pressure components are separated:
https://github.com/user-attachments/assets/249e9814-335f-49bb-999e-73a6f95fcf37
Likely introduced by https://github.com/CliMA/Oceananigans.jl/pull/3574.
@tomchor recommends changing the default setting of NonHydrostaticModel to hydrostatic_pressure_anomaly = CenterField(grid) (at least until issue has been addressed).
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.
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?
@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.
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.
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
https://github.com/user-attachments/assets/fef6e145-d7d9-4e28-bbb0-e22181344dd8
hydrostatic_pressure_anomaly = CenterField(grid)
https://github.com/user-attachments/assets/6bd55ad5-5974-488e-a06c-341849df2e43
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.
Another test that reveals a problem with non-separated pressure is a simple horizontal wall in a 2D setup.
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: https://github.com/CliMA/Oceananigans.jl/pull/3188
Indeed the choice of divergent boundary flow has been made so that there is no tracer leakage through the boundary. One can try using https://github.com/CliMA/Oceananigans.jl/pull/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.
@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?
@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 https://github.com/Yixiao-Zhang/Oceananigans.jl/commit/c7983b8002b91cd5939018a7c999eae77e2105ac
This should work and fix the numerical issues, but has not been extensively validated yet.
@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?
I also think it would be appropriate to add a warning remarking about the approximate nature of the current pressure solver (at least until https://github.com/CliMA/Oceananigans.jl/pull/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?
@glwagner, am I correct that this warning is only needed when 1)
hydrostatic_pressure_anomaly = nothingand 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
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).
@xkykai would you like help with #3188 ?
Yes I will work on it.
Would you like any help?
@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.
I made a PR that adds a warning and changes the default for immersed boundary grids: https://github.com/CliMA/Oceananigans.jl/pull/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.
Just wanted to comment here for (self)awareness that I also stumbled into this issue while using an immersed boundary for an ice shelf (which seems especially bad when your immersed boundary has a buoyancy forcing!) and the comments above were particularly useful. The hydrostatic_pressure_anomaly = CenterField(grid) fix did significantly reduce the boundary divergence (excluding melt) although I think the conjugate gradient pressure solve may be needed in this particular case for improving the accuracy of the solution. Thanks for the help!
Just wanted to comment here for (self)awareness that I also stumbled into this issue while using an immersed boundary for an ice shelf (which seems especially bad when your immersed boundary has a buoyancy forcing!) and the comments above were particularly useful. The hydrostatic_pressure_anomaly = CenterField(grid) fix did significantly reduce the boundary divergence (excluding melt) although I think the conjugate gradient pressure solve may be needed in this particular case for improving the accuracy of the solution. Thanks for the help!
of course! Thanks for jumping into this conversation!
The hydrostatic_pressure_anomaly = CenterField(grid) fix did significantly reduce the boundary divergence (excluding melt)
This is default now right? So a fix should not be needed (just checking...)
We might be able to close this issue
No, "hydrostatic_pressure_anomaly = nothing," seems to be the default in v0.94, so if I were to run the default, the hydrostatic_pressure_anomaly would not be computed (separate from the "nonhydrostatic_pressure = CenterField(grid)" and seems to lead to the spurious immersed boundary buoyancy source. Let me know if I misunderstood something.
Here's the logic:
https://github.com/CliMA/Oceananigans.jl/blob/2c955dac657ffa3760ffb46c5ff7a36c0215c59b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl#L143-L160
Note that we also noticed there is an issue even not on immersed boundary grid (which is quite subtle), so we decided to separate the pressure computation any time there is buoyancy. That's documented on #3795