Maintaining discrete hydrostatic balance

[mwarusz]

Description

Numerics should maintain hydrostatically balanced states to a very high precision. For that to happen the pressure gradient term must be in discrete balance with the gravity source term and the numerical flux shouldn't break this balance. Calculating rho -= 1 / g dp/dz using DG results in discrete balance but Rusanov flux breaks it due to the penalty term wavespeed * (rho+ - rho-) and discontinuous rho.

Ideas:

• calculate p from rho using discrete numerical integration that uses projection to a space that is one order lower to make the state continuous
• use a different numerical flux (current strategy is 1) to test central fluxes and split in IMEX between Rusanov & central and 2) implement Roe fluxes later)
• approach based on https://arxiv.org/pdf/1511.08739

Additional context

Add any other reasons why this feature will help with the goals of CLIMA.

For CLIMA Developers

• [ ] This feature can be added (if it cannot be, explain why in a comment below -- lack of technical expertise, not relevant to the scope of this project, too ambitious)
• [ ] There is a timeline for when this feature can be implemented
• [ ] The feature has been (or will be) implemented (Please link the PR)

[tapios]

A couple important things:

• It’s important for a geophysical model not to generate flow spontaneously from a hydrostatic rest state, even when density is discontinuous in the rest state. This probably means we need to rethink numerical fluxes.
• The continuity equation is a tracer equation for a constant tracer (e.g., ‘chi=1’). The equations (and numerical fluxes) for density and tracers should be treated consistently.

[simonbyrne]

• It’s important for a geophysical model not to generate flow spontaneously from a hydrostatic rest state, even when density is discontinuous in the rest state. This probably means we need to rethink numerical fluxes.

I would argue that if density is discontinuous, then it isn't hydrostatic.

[tapios]

I would argue that if density is discontinuous, then it isn't hydrostatic.

Why? We use models with discontinuous density in hydrostatic balance all the time. The classical models of atmospheric flow are hydrostatic with discontinuous densities (e.g., Phillips two-layer model, https://www.whoi.edu/cms/files/adoucette/2006/4/Pedlosky_12_824_Ch10_9446.pdf).

[christophernhill]

@mwarusz does the Rusanov term cause the solution to grow in an unbounded way? If so maybe some tests with a smaller Rusanov wave speed is worth trying - some ramblings below...

Numerically we are all used to having small truncations that need to be stabilized. Rusanov does that, but my impression is that we need to "lie" a little about the wave speed otherwise everything falls apart. In the ocean we have settled on a Rusanov term that has a speed slower than the external gravity wave speed. It still helps stabilize Rieman problem "jumps" but without crazy CFL restrictions. Its not entirely physical, but we aren't representing that physics because our tilmestep is too long?

Some of the initial thinking goes back to Godunov, but there my impression is that the notion was some realistic physics would compensate these jumps. In theory in the atmosphere that is acoustic waves, but we don't really have those. So maybe a test with a more modest wave speed is worth trying?

We are still working towards proving this fully works for ocean of course!!

[mwarusz]

@christophernhill My experience has been that the Rusanov term doesn't destabilize the solution in this case. It simply generates unphysical flow. But in general I agree that we need a more careful look at wavespeeds when we have unresolved acoustic waves that we don't really care about.

[simonbyrne]

Ah, I misunderstood: it seems there is an inherent tension between what is mathematically possible, and what is representible in a DG context.

[tapios]

@lcw Can we use Roe fluxes? From the discussion this morning, it wasn‘t clear to me if that is what you did with the tracers. If so, does it work with mass? It seems like a straightforward choice for a consistent flux that vanishes for vanishing velocities (which we need).

[lcw]

@lcw Can we use Roe fluxes? From the discussion this morning, it wasn‘t clear to me if that is what you did with the tracers. If so, does it work with mass? It seems like a straightforward choice for a consistent flux that vanishes for vanishing velocities (which we need).

I have never used a Roe flux in DG for Euler. The code definitely supports defining it. I guess one would have to try it and see how it works.

@fxgiraldo Do you have any experience with a Roe flux in NUMA?

[tapios]

Seems like a Roe flux would be worth a try, not only in terms of preserving hydrostatic balance, but also in terms of stability: https://www.sciencedirect.com/science/article/pii/S0021999116305642

[fxgiraldo]

@lcw Can we use Roe fluxes? From the discussion this morning, it wasn‘t clear to me if that is what you did with the tracers. If so, does it work with mass? It seems like a straightforward choice for a consistent flux that vanishes for vanishing velocities (which we need). I have never used a Roe flux in DG for Euler. The code definitely supports defining it. I guess one would have to try it and see how it works. @fxgiraldo Do you have any experience with a Roe flux in NUMA?

@lcw Not in NUMA since we have never seen a need for it. However, this is one of the fluxes I used in my PhD work. There are pros and cons to using this flux (expense being one disadvantage and the other being that it admits expansion shocks but there is an easy fix and since we are not supersonic then we may not need it. Will have to refresh my memory on this flux). Bottom line, there are many fluxes to choose from, like Lucas said at the monday meeting: Rusanov is just the simplest one and that's what we started with.

[simonbyrne]

To make this concrete (since it isn't stated in the issue), we are defining "discrete hydrostatic balance" to such that the DGModel operator returns zero tendencies (up to machine precision).

The aims in order of priority:

1. Our current reference state should satisfy this property
2. Be able to construct a reference state from an arbitrary continuous temperature profile
3. Be able to construct a reference state from an arbitrary discontinuous temperature profile

[bischtob]

@mwarusz, @thomasgibson and I put together a deep and shallow baroclinic-wave here. Please make sure we test against this one as well and compare to NUMA's test cases if possible.

[tapios]

The discussion in section 5 of https://arxiv.org/pdf/2005.02317.pdf may be relevant here (h/t @OsKnoth). It makes a case for a split form of the flux divergence that is both stable and free-stream preserving (i.e., discretely preserving a constant velocity field).

[mwarusz]

In context of IMEX splitting a simple way to have a numerical flux that deals properly with density discontinuities is to use central flux on the linear part and Rusanov flux on the nonlinear part (using only advective wavespeed). There is still a couple of issues

• stability, is this stable for long-term realistic simulations (to check we need a working test case for the GCM, idealized tests seem to indicate that it is stable but requires more dissipation)
• for 1D IMEX we would need PDE level splitting that includes directions which is addressed by #1037

However, this is promising since using central flux on the linear part is also required to construct the Schur complement and it would be an easy fix.