Don't precompute hydrostatic pressure in NonhydrostaticModel

[glwagner]

This PR inserts the vertical buoyant acceleration directly into the vertical momentum equation in NonhydrostaticModel, rather than integrating to find the hydrostatic pressure and inserting its gradient into the horizontal momentum equations.

I just did the easiest thing right because I'm curious if regression tests pass. If they do, we can refactor NonhydrostaticModel to eliminate hydrostatic pressure, speed up the model, and reduce its memory footprint without too much pain (🎉)

If the regression tests don't pass, we will unfortunately have slightly more pain ahead of us in refactoring the regression tests. In addition to the above advantages, we also need to eliminate vertical integrals in NonhydrostaticModel to permit 2D distributed memory parallelization. This is because PencilFFTs only allows parallelization across dimensions higher than 1 (or in other words, we cannot parallelize in x). Thus with a vertical integral, we can only parallelize easily across y. If we eliminate the vertical integral, we'll be able to parallelize in y, z.

[glwagner]

Heh, well the stratified_fluid_remains_at_rest_with_tilted_gravity_temperature_tracer test fails.

The log suggests that the problem is confined to a few grid points, since the mean is correct:

[2021/07/29 16:10:56.038] INFO  Simulation is stopping. Model time 1 hour has hit or exceeded simulation stop time 1 hour.
--
  | [2021/07/29 16:10:57.288] INFO  ∂T∂z * g̃[2] = 0.0001065658790393914, mean(∂y_T) = 0.0001065658790360618, Δ = 3.32959842423014e-15 at t = 1 hour with θ=1°
  | [2021/07/29 16:10:57.288] INFO  ∂T∂z * g̃[3] = 0.006105155121314884, mean(∂z_T) = 0.006105155121314366, Δ = 5.178149575790769e-16 at t = 1 hour with θ=1°
  | Tilted gravity: Test Failed at /var/lib/buildkite-agent/builds/tartarus-5/clima/oceananigans/test/test_dynamics.jl:290
  | Expression: all(∂T∂z * g̃[2] .≈ interior(∂y_T))

Oceananigans.jl/test/test_dynamics.jl

Line 241 in 9ecddac

function stratified_fluid_remains_at_rest_with_tilted_gravity_temperature_tracer(arch, FT; N=32, L=2000, θ=60, N²=1e-5)

I am skeptical about the value of this test. In general, discretizations do not necessarily represent the balance states of the continuous systems they approximate. Thus in general this test should fail! It will succeed of course in the case that the discrete system shares a balanced state with the continuous state.

I think this is flawed as a test because, while certainly formulating discrete systems that share balanced states with their continuous systems is desirable, I'm not sure it's something we want to guarantee via CI. We should be allowed to violate this prescription if the discrete system has other more important desirable properties.

It could be a nice validation test though.

In this particular case, the nonhydrostatic useful is most often used for simulations of vigorous turbulence in which weak flows induced by a violation of hydrostatic balance are likely irrelevant.

[glwagner]

The change in this PR might conceivably produce changes in the dynamics because it changes the temporal discretization of the system. In particular, the pressure correction method corresponds to a first-order backward Euler discretziation of the pressure gradient term. On the other hand, since we calculate the hydrostatic pressure analytically coincident with other terms on the right hand side, it's treated with an explicit discretization (though the discretization is still first-order for the quasi-Adams-Bashforth scheme, it can have different error properties).

Which time discretization is better? I don't know...

[glwagner]

@tomchor since you're the main person using the tilted gravity feature, I'm wondering if you can help provide some insight into this "stratified fluid at rest" test. The main issue is that the dynamics can be "correct" but the test can fail. I feel its a bad test for this reason.

[tomchor]

@tomchor since you're the main person using the tilted gravity feature, I'm wondering if you can help provide some insight into this "stratified fluid at rest" test. The main issue is that the dynamics can be "correct" but the test can fail. I feel its a bad test for this reason.

I don't have much to add to the discussion. I agree with you that a balanced state in a continuous system doesn't necessarily translate exactly to a discrete one. When I (or Ali?) came up with this test I figured this translation error would be small enough to be acceptable, and when the test actually passed I was happy enough with that.

So if this translation error is indeed large enough with the new solver that the tests don't pass I'm very much okay with changing the test. I can't, for the moment, think of another simple test to replace it though. My best guess is to do something similar to what I did for the rotated Coriolis: solve a system with gravity pointing upwards and then the same system with gravity pointing to the x or y direction and see if they match after the proper rotation.

[glwagner]

@tomchor since you're the main person using the tilted gravity feature, I'm wondering if you can help provide some insight into this "stratified fluid at rest" test. The main issue is that the dynamics can be "correct" but the test can fail. I feel its a bad test for this reason.

I don't have much to add to the discussion. I agree with you that a balanced state in a continuous system doesn't necessarily translate exactly to a discrete one. When I (or Ali?) came up with this test I figured this translation error would be small enough to be acceptable, and when the test actually passed I was happy enough with that.

So if this translation error is indeed large enough with the new solver that the tests don't pass I'm very much okay with changing the test. I can't, for the moment, think of another simple test to replace it though. My best guess is to do something similar to what I did for the rotated Coriolis: solve a system with gravity pointing upwards and then the same system with gravity pointing to the x or y direction and see if they match after the proper rotation.

Can't we just test directly that the output of x_dot_b is as expected (as well as the others)? Along with a test that the constructor works without error this seems sufficient.

More complicated integration tests, like testing that the discrete system has a balanced state analogous to the continuous one, seem better suited for a validation test, I think.

[tomchor]

Can't we just test directly that the output of x_dot_b is as expected (as well as the others)? Along with a test that the constructor works without error this seems sufficient.

More complicated integration tests, like testing that the discrete system has a balanced state analogous to the continuous one, seem better suited for a validation test, I think.

Yeah, if you think that's enough for CI then I'm definitely okay with that!

[tomchor]

Removed the "fluid at rest" test that was failing and included another (simpler) one for tilted buoyancy along with constructors. The new test checks that x_dot_g_b, etc., match for two different domains: one with z-aligned gravity and another with y-aligned gravity.

The only downside is that I only included the x_dot_g_b test for model=BuoyancyTracer(). Another separate test is needed for model=SeawaterBuoyancy() if we think that's necessary.

(Constructors are tested for both models, though.)

[glwagner]

Removed the "fluid at rest" test that was failing and included another (simpler) one for tilted buoyancy along with constructors. The new test checks that x_dot_g_b, etc., match for two different domains: one with z-aligned gravity and another with y-aligned gravity.

The only downside is that I only included the x_dot_g_b test for model=BuoyancyTracer(). Another separate test is needed for model=SeawaterBuoyancy() if we think that's necessary.

(Constructors are tested for both models, though.)

Thanks @tomchor !!

[tomchor]

@glwagner I'm trying to use this branch for my research so I tried resolving the conflicts as best as I could. I hope you don't mind. Otherwise I can always revert back!

[tomchor]

The errors of the cpu-simulation-tests are of this kind:

JLD2 output writer [CPU]: Test Failed at /var/lib/buildkite-agent/builds/tartarus-2/clima/oceananigans/test/test_jld2_output_writer.jl:135
--
  | Expression: wT == zero(FT)
  | Evaluated: -1.2037062f-34 == 0.0

which is something that should be solved by replacing == with \approx, but I wonder why this used to pass before but isn't passing now. In other words: I'm not sure if merely replacing equality with approximate equality is the right way to solve this test.

[tomchor]

Just a note: it seems like this PR is currently changing the dynamics. For example comparing the convecting plankton example in this PR has significantly more noise than the same example for the main branch. I haven't checked every example, but this also appears to happen in the Langmuir turbulence example, and especially in the internal wave one, which looks pretty bad...

It doesn't show up in all of them though. For example the wind mixing one for some reason.

I'm not sure if I did something wrong when merging the main branch, or if the effects were already there before (I don't remember checking before; @glwagner did you?)

[johncmarshall54]

Is this a GPU problem? Is it still there when running on CPUs? We had some very early issues with this. Ali may remember.

…

On Thu, Dec 16, 2021, 2:56 PM Tomas Chor ***@***.***> wrote: Just a note: it seems like this PR is currently changing the dynamics. For example comparing the convecting plankton example in this PR <https://clima.github.io/OceananigansDocumentation/previews/PR1910/generated/convecting_plankton/> has significantly more noise than the same example for the main branch <https://clima.github.io/OceananigansDocumentation/stable/generated/convecting_plankton/>. I haven't checked every example, but this also appears to happen in the Langmuir turbulence example <https://clima.github.io/OceananigansDocumentation/previews/PR1910/generated/langmuir_turbulence/>, and *especially* in the internal wave one <https://clima.github.io/OceananigansDocumentation/previews/PR1910/generated/internal_wave/>, which looks pretty bad... It doesn't show up in all of them though. For example the wind mixing one <https://clima.github.io/OceananigansDocumentation/previews/PR1910/generated/ocean_wind_mixing_and_convection/> for some reason. I'm not sure if I did something wrong when merging the main branch, or if the effects were already there before (I don't remember checking before; @glwagner <https://github.com/glwagner> did you?) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1910 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKXUEQU53Y5LV3FCQU5CWUTURHVZHANCNFSM5BHC6H6A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[tomchor]

Is this a GPU problem? Is it still there when running on CPUs? We had some very early issues with this. Ali may remember.
…

The equality failing in the tests happens for both CPUs and GPUs. The artifacts that appear in some simulations were only tested in CPUs I think. So it seems that GPUs aren't to blame here.Both issues (

[glwagner]

Is this a GPU problem? Is it still there when running on CPUs? We had some very early issues with this. Ali may remember.

I may not have been around for those tests! Here we're trying a different numerical algorithm that doesn't separate the hydrostatic and non-hydrostatic pressure. The numerics is different this way --- the non-hydrostatic pressure is in some sense treated "implicitly" (eg computed at t^{n+1} in order to advance the velocity field from n to n+1). The hydrostatic pressure, before this PR, was computed explicitly, at time-step t^n. The difference between the two numerical schemes could explain the issue @tomchor's observed (or there could be another bug).

[glwagner]

Unfortunately the PR previews are no longer showing so I'll try to update this PR...

[glwagner]

Maybe we can add a keyword argument to the NonhydrostaticModel constructor precompute_hydrostatic_pressure_anomaly = false or something like that. If it needs to be true by default (for numerical reasons if @tomchor's observations hold out) we can try to switch it to false if arch isa MultiArch and perhaps emit a warning.

[glwagner]

Looks like this is a bad idea, so I'm closing.

[tomchor]

Looks like this is a bad idea, so I'm closing.

I'm not sure it's a bad idea. I've been having some issues with spurious waves in rotated domains and I'm thinking this might be the cause. (Since we're modifying the direction gravity is acting on, but we're still integrating the pressure in the model's z direction.)

I still haven't been able to test my rotated domain on this branch, so I'm not sure. But if indeed this is the culprit, it might be a good idea to have a flag that turns the hydrostatic separation off for rotated domains.

[glwagner]

Hmm, ok I see that.

I think this PR is a little stale but when you find time to document the problem, maybe open an issue and we can proceed from there?