Climatemachine.jl: Filters need to be applied correctly to state variables

[bischtob]

Description

When running a solid-body rotation test case, the flow does not remain at rest and a hydrostatic state shows spurious warming and cooling at the top of the domain and at the element edges. The problem is more severe when diffusion or exponential filters are applied and the problem gets more severe as the filter exponent is increased. The problem occurs for explicit and IMEX-type solvers alike, but the 1DIMEX is what we need to use for the GCM for now. In addition, right now, on the sphere, the simulations will no stably integrate without exponential filters. The problem occurs in Float32 and Float64.

This most likely causes significant issues in the Held-Suarez experiment case.

Steps to reproduce

Steps to reproduce the behavior:

1. Run the solid-body test from here.
2. Use VizCLIMA to genereate plots using this script here.

Expected behavior

In a solid-body test case we should expect the state variables to remain unchanged to high precision. This does not seem to be the case.

Screenshots

Example 1: No sponges, No forcing, No diffusion, No hyperdiffusion, Isothermal reference and initial temperature profile, exponential filter exponent: 64

The figure shows vertical slices of the difference between the final state at 2hrs and the initial state. Top-left: zonal-mean zonal velocity difference, top-right: zonal-mean meridional velocity difference, bottom-left: zonal-mean vertical velocity difference, bottom-right: zonal-mean temperature difference. All fields are evaluated at the equator of the planet.

Example 2: No sponges, No forcing, No diffusion, No hyperdiffusion, Isothermal reference and initial temperature profile, exponential filter exponent: 16

The figure shows vertical slices of the difference between the final state at 2hrs and the initial state. Top-left: zonal-mean zonal velocity difference, top-right: zonal-mean meridional velocity difference, bottom-left: zonal-mean vertical velocity difference, bottom-right: zonal-mean temperature difference. All fields are evaluated at the equator of the planet.

Example 3: No sponges, No forcing, No diffusion, No hyperdiffusion, Stably stratified reference and initial temperature profile, exponential filter exponent: 64

The figure shows vertical slices of the final state at 1hr at various locations on the sphere. Top-left: zonal-mean zonal velocity, top-right: zonal-mean meridional velocity, bottom-left: zonal-mean vertical velocity, bottom-right: zonal-mean temperature. Note the impact on the temperature field at the top of the domain (the small kink).

Environment and error messages

No error messages, just physical inconsistencies.

Additional context

Add any other context about the problem here.

For CLIMA developers

• [x] The bug can be reproduced
• [x] Labels have been added and developers assigned
• [ ] The bug can be fixed (if not, please explain why not in a comment below)
• [ ] There is a timeline to fix the bug
• [ ] The bug has been fixed (Please link the PR)

[smarras79]

This is exactly what @slifer50 is experiencing with the tropical cyclone test in a LES domain.

I'll create a balanced flow test for the flow in a box with topography. You can then export it to the sphere.
This test verifies that the discrete solution is hydrostatically balanced to machine precision. I.e. given a still atmosphere over a mountain, this should stay still (w=1e-16 m/s) forever.
If we don't pass this test, it means that the background state it is either not correctly setup or the "numerical solution is not balanced". l'll try push it by today.

[jkozdon]

Is the situation any better with Float64?

[bischtob]

Is the situation any better with Float64?

The problem occurs in Float32 and Float64.

[jkozdon]

The problem occurs in Float32 and Float64.

Thanks for checking.

[thomasgibson]

@mwarusz aptly pointed out that the issue is perhaps related to the filter and I reproduced his experiment. Below is the same solid body case at 2hr with 10 vertical elements, Float 64, and no filtering applied.

So the deviations are much smaller. However, we need filters to keep the numerics stable. So a working hypothesis is that the filter should not be applied to the full state; instead, we think it should be applied to the perturbation: \Delta q := q^{n+1} - q^{n}.

We will keep this issue updated as we find out more.

[bischtob]

Great work @mwarusz and @thomasgibson, I can also reproduce those results as well with Float64 and no filters. Without filtering this setup goes unstable for me at about 0.15 days. See below:

[smarras79]

@bischtob what diffusion method are you using to overcome the lack of filter?

[bischtob]

@bischtob what diffusion method are you using to overcome the lack of filter?

None. We need to make sure everything is right before adding another term to the equation, I think.

[smarras79]

@bischtob what diffusion method are you using to overcome the lack of filter?

None. We need to make sure everything is right before adding another term to the equation, I think.

You have gravity waves and no filter. without a stabilization mechanism it will blow up. This is DG which is not dissipative enough to preserve numerical stability with important motion

[bischtob]

@bischtob what diffusion method are you using to overcome the lack of filter?

None. We need to make sure everything is right before adding another term to the equation, I think.

You have gravity waves and no filter. without a stabilization mechanism it will blow up. This is DG which is not dissipative enough to preserve numerical stability with important motion

Totally correct :)! The point of this issue is to understand if and when filtering stabilizes the model without adding other sources. As you can see above, the filtering does not stabilize the situation (figures at the top). Finding out what is going on, aka whether it is the filters that are wrong, or something else, is what I think we should do before I add another term to the equations. Let me know if that makes any sense.

[thomasgibson]

@mwarusz has a potential fix for how the filter is applied to the state. The commit can be found here: https://github.com/CliMA/ClimateMachine.jl/commit/54b3b31c3a954f091f931cc48080a5e49d8418ef

Basically what this fix does is it changes how the filters are applied to the state vector. Rather than applying it to q directly, the reference state q_0 is subtracted out, leaving just the perturbation: \Delta q, where q = q_0 + \Delta q. The filter is applied to \Delta q rather than q = q_0 + \Delta q. Afterward, the q_0 data is added back.

With this fix, the solid body test case looks much better. Here is the vertical slice deviation (using Float64 and running for 1000 days):

Here is the kinetic energy profile (the spike is there but note the magnitude is still 10^{-14}:

All of these results were generated using VizCLIMA.jl and diagnostics were generated every 100 steps.

[bischtob]

@thomasgibson, @mwarusz, this looks MUCH better and it makes a lot more sense theoretically, I think, to apply the filter to the stare-ref_state.

[bischtob]

@mwarusz, @thomasgibson, once the reference state differs from the initial condition (BOTH being in hydrostatic balance) the simulation does not remain balanced anymore. This may be unrelated to the filtering, but may be related to #1169. @jkozdon, any ideas?