Don't separate the pressure into hydrostatic and nonhydrostatic in NonhydrostaticModel

[tomchor]

This PR changes NonhydrostaticModel to use a single pressure variable. Prior to this PR, pressure is decomposed into a "hydrostatic anomaly" and a "non-hydrostatic" component. The vertical momentum equation in NonhydrostaticModel is

$$ \partial_t w = - \partial_z p + b + G_w' $$

where $G_w'$ collects the contributions to the vertical momentum tendency other than pressure gradient and buoyancy. Prior to this PR, the kinematic pressure $p$ is decomposed into (note that in this Boussinesq formulation, a hydrostatic component has already been removed from $p$)

$$ p = p_h' + p_n $$

with the hydrostatic anomaly $p_h'$ defined by

$$ p_h' = - \int^0_z b \mathrm{d} z $$

Using the hydrostatic pressure anomaly eliminates buoyancy from the vertical momentum equation, such that

$$ \partial_t w = - \partial_z p_n + G_w' $$

Buoyancy forces enter the dynamics via horizontal gradients of the hydrostatic pressure anomaly. For example, x-gradient of the kinematic pressure becomes

$$ \partial_x p = \partial_x p_h' + \partial_x p_n $$

This decomposition is advantageous for two reasons. First, in a hydrostatic model the vertical momentum equation reduces to the equation for $p_h'$. This means that switching from a hydrostatic to non-hydrostatic model is particularly simple given this decomposition.

Second --- and this must be evaluated --- it's possible to carefully control the evaluation of the hydrostatic integral so that a resting stratified fluid remains at rest, even in the presence of complex bathymetry. When we use the "MITgcm algorithm" we achieve this perfectly, even with partial cell bathymetry.

This PR proposes to eliminate the pressure decomposition so that there is only one pressure. In principle, this has a computational advantage because the hydrostatic pressure integral does not need to be evaluated (in practice, this computation has a negligible cost). It also reduces the number of memory loads that take place in the momentum advection kernels (though these are typically domained by advection scheme, so this may not matter except for centered advection schemes). Also in principle, it would allow 3D domain decompositions for distributed computations, in addition to 2D (but again, these are rarely used because typical ocean domains are shallow and wide, rather than deep and narrow). Having a single pressure also simplifies diagnostics. Finally, and perhaps most importantly, we can avoid allocating memory for an additional 3D field. In the absolute best case scenario of a model with no tracers and pure implicit dissipation, this means we go from 14 3D fields (9 for prognostic momentum + tendencies, 4 (?) for nonhydrostatic pressure including scratch variables for FFTs, and 1 for hydrostatic pressure) to 13 3D fields. So it saves about 7%. In more typical situations with LES closure and one active tracer, the savings is more marginal: we go from 19 3D fields to 18 3D fields, and thus have 5% more memory. Note also that more scratch variables are needed for FFTs in domains that are not horizontally-periodic. (This memory bookkeeping should be double checked.)

In summary, the main advantages and risks of this PR are:

Risk: loss of accuracy in scenarios that are dominated by hydrostatic balance. This issues may be exacerbated by experimental or not-yet-existing features, such as: cut cell and partial cell bathymetry, curvilinear grids, nonlinear free surfaces, reduced precision arithemetic, etc.

Advantage: memory savings of at most 7% but more typically 1-5%, and a cleaner code and user interface.

Note, there was another attempt to coalesce the pressures in #1910. However, buoyancy was not reconstructed properly in the momentum equations (buoyancy is at tracer points; thus the buoyancy force has to be interpolated to by included in the vertical momentum balance). and thus the discretization was incorrect and produced spurious dynamics. This bug was fixed by #3079.

In all of our tests so far, the "dynamics seem clean". However, it's not clear whether there are unforeseen issues in scenarios that we haven't tested, or rather are impossible to test because the feature does not exist yet (such as accurate reduced precision algorithm or nonlinear free surface). Thus we should consider this PR carefully.

[tomchor]

Docs are previewing here: https://clima.github.io/OceananigansDocumentation/previews/PR3080/

I checked all the examples with NonhdyrostaticModel and they all look the same as they do on the stable branch.

Furthermore, the few tests failures that we have are all something like

JLD2 output writer [CPU]: Test Failed at /var/lib/buildkite-agent/builds/tartarus-1/clima/oceananigans/test/test_jld2_output_writer.jl:131
--
  | Expression: wu == zero(FT)
  | Evaluated: -3.009265538105056e-35 == 0.0

i.e. very small approximation errors that aren't indicative of any significant errors in the model.

In other words, I think this is working well! I vote we simplify the model and get rid of the hydrostatic separation.

@glwagner as you mentioned, this isn't a trivial change. If you wanna move forward with it, feel free to push to this PR or close this one and open another. I can also help if you want, just lmk what I should focus on.

PS: Just like we did in #1910 we might need/want to replace the stratified_fluid_remains_at_rest_with_tilted_gravity_buoyancy_tracer() test for something simpler.

[glwagner]

annoying that the tests need a little refactoring too, but maybe its all for the best.

[tomchor]

With this last commit both the hydrostatic and nonhydrostatic models are working locally for me and they only have a model.pressure field; no more NamedTuples of pressure. But let's see if tests are passing.

The PressureField functions (formerly PressureFields) should probably be moved from field_tuples.jl since they don't return tuples anymore. Any suggestion as to where they might be moved?

[glwagner]

Do we need this at all, or is it simpler to put it in the model constructor?

[tomchor]

If that's the only place that method is called, then I think better to put it directly in the model constructor. But I'm not sure that's the case.

[tomchor]

What's called in the constructor is

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl

Line 176 in 6c38841

pressure = PressureField(pressure, grid, boundary_conditions)

with pressure==nothing by default. So I guess the reason we don't call PressureField(grid, bcs) directly is to make it easier for users to pass a pressure directly to the user-facing model constructor and override the default.

[tomchor]

The issue was that buoyancy was not reconstructed properly in the vertical. (There was also a problem with reconstructing buoyancy in the horizontal, but this only affects tilted domains.) So there was a bug and the discretization was incorrect.

I'm aware of that. My point was more that I don't quite remember how the issue was manifested in the dynamics (the recognition of which was what prompted us to abandon #1910).

@tomchor I updated the PR description. Feel free to edit it further.

Thanks, that's a great detailed description. I only added one item to the advantages: simpler code and user interface.

[glwagner]

My point was more that I don't quite remember how the issue was manifested in the dynamics

That's ok. The dynamics associated with incorrect numerics are not important.

[glwagner]

What do you mean by cleaner code? You mean update_state!? We need all of these functions still for the hydrostatic model so I don't think on the whole there's much of a change to the source code.

[tomchor]

What do you mean by cleaner code? You mean update_state!? We need all of these functions still for the hydrostatic model so I don't think on the whole there's much of a change to the source code.

I'm referring to the fact that dealing with just one pressure makes for cleaner/simpler code in general. But feel free to remove that statement if you don't agree.

[glwagner]

Just trying to understand what was meant by that since it's being mentioned as one of the major advantages (as opposed to the many minor advantages like computational efficiency)

[glwagner]

I think it makes sense to merge #3188 before this PR --- if we can show that the new pressure solver in #3188 is unaffected by this PR for some examples of interest (maybe including cases dominated by hydrostatic pressure) then I think we can move forward.

[tomchor]

I think it makes sense to merge #3188 before this PR --- if we can show that the new pressure solver in #3188 is unaffected by this PR for some examples of interest (maybe including cases dominated by hydrostatic pressure) then I think we can move forward.

Sounds good. Do you have a sense of how close #3188 is to being ready for merging?

[tomchor]

I'll leave this here for the record. I'm currently experiencing the first significant dynamical difference I've seen so far between the model with and without the pressure separation.

In a simulation where I'm studying flow past an obstacle (therefore with immersed boundaries) the simulation runs fine on this branch, but (everything else being the same) doesn't finish running on main. On main the velocities keep increasing, leading to a progressive decrease in Δt to satisfy CFL condition, but the velocities keep increasing despite that, with smaller and smaller Δts. (I believe that's called a slow blow-up?)

So this is a case where the simulation fails on main, but is successful in this branch.

The simulation is far too complex to post here, but I'll try to come up with a MWE that reproduces it.

[glwagner]

I'll leave this here for the record. I'm currently experiencing the first significant dynamical difference I've seen so far between the model with and without the pressure separation.

In a simulation where I'm studying flow past an obstacle (therefore with immersed boundaries) the simulation runs fine on this branch, but (everything else being the same) doesn't finish running on main. On main the velocities keep increasing, leading to a progressive decrease in Δt to satisfy CFL condition, but the velocities keep increasing despite that, with smaller and smaller Δts. (I believe that's called a slow blow-up?)

So this is a case where the simulation fails on main, but is successful in this branch.

The simulation is far too complex to post here, but I'll try to come up with a MWE that reproduces it.

Interesting that is with immersed boundaries. The pressure solver isn't correct in that case so I'm not sure how to interpret that. The error should be smaller on main (because it's only the correction to the hydrostatic anomaly that is wrong on main).

The main uncertainty is how this PR will interact with #3188. We could explore using the new immersed pressure solver on this branch to test that out.

[tomchor]

The main uncertainty is how this PR will interact with #3188. We could explore using the new immersed pressure solver on this branch to test that out.

Is that branch working? i.e. do I just have to merge both branches and run? If so, I can do that and report back.

[glwagner]

I haven't run the code myself, but @xkykai has produced quite a few animations with that new pressure solver...

[glwagner]

Hmm, there are a lot of conflicts here, probably because of the MPI communication PR. Should we re-open this or resolve the conflicts? @tomchor @simone-silvestri

[tomchor]

Hmm, there are a lot of conflicts here, probably because of the MPI communication PR. Should we re-open this or resolve the conflicts? @tomchor @simone-silvestri

I think either option only makes sense if we're at the point of actually considering merging this. Are we? Otherwise it's probably best to just let this PR sit around then.

[glwagner]

@tomchor I think we should revisit this and get it merged. It doesn't seem like the immersed Poisson solver is progressing right now, so I don't think it makes sense to keep this on hold for it --- especially because features like triply periodic simulations require this change. Happy to help resolve merge conflicts or opening a new PR whichever makes more sense.

[tomchor]

@tomchor I think we should revisit this and get it merged. It doesn't seem like the immersed Poisson solver is progressing right now, so I don't think it makes sense to keep this on hold for it --- especially because features like triply periodic simulations require this change. Happy to help resolve merge conflicts or opening a new PR whichever makes more sense.

I'm leaning towards opening another PR and making the non-separation optional. The rationale is that even though the immersed Poisson solver isn't progressing right now (and I agree we shouldn't wait for it), when it does progress it may turn out that it works better with the original (i.e. hydrostatic pressure separation) algorithm. Plus if we keep the current Poisson solver as default, then we don't need to re-run the regression tests.

What do you think?

[glwagner]

I think that's a fine strategy. We can add a kwarg to NonhydrostaticModel called hydrostatic_pressure_anomaly. We can set it to CenterField(grid) to preserve existing behavior, or set it to nothing to avoid the separation. And we should probably make nothing default so that triply periodic problems can be done out of the box. Then we don't have to re-do the regression tests either because we preserve existing behavior...

I think that's also a less invasive change than this PR because we don't have to change pressures to pressure everywhere, hmm.

Since you've done most of the legwork I think you have prerogative to open a new PR if you like (and I can help once you do).

[tomchor]

I think that's a fine strategy. We can add a kwarg to NonhydrostaticModel called hydrostatic_pressure_anomaly. We can set it to CenterField(grid) to preserve existing behavior, or set it to nothing to avoid the separation. And we should probably make nothing default so that triply periodic problems can be done out of the box. Then we don't have to re-do the regression tests either because we preserve existing behavior...

I think that's also a less invasive change than this PR because we don't have to change pressures to pressure everywhere, hmm.

Agreed.

Since you've done most of the legwork I think you have prerogative to open a new PR if you like (and I can help once you do).

Thanks, but I unfortunately won't have time to dedicate to this for at least a few weeks. So please feel free to start the PR!

[glwagner]

I think that's a fine strategy. We can add a kwarg to NonhydrostaticModel called hydrostatic_pressure_anomaly. We can set it to CenterField(grid) to preserve existing behavior, or set it to nothing to avoid the separation. And we should probably make nothing default so that triply periodic problems can be done out of the box. Then we don't have to re-do the regression tests either because we preserve existing behavior...
I think that's also a less invasive change than this PR because we don't have to change pressures to pressure everywhere, hmm.

Agreed.

Since you've done most of the legwork I think you have prerogative to open a new PR if you like (and I can help once you do).

Thanks, but I unfortunately won't have time to dedicate to this for at least a few weeks. So please feel free to start the PR!

Great. I think the PR is nearly finished over at #3574. Give it a look over when you have time.

[tomchor]

I'm closing this since we merged #3574