Skamarock-Klemp Hydrostatic is broken

[tommbendall]

Unfortunately the case studies haven't been tested for a few months, and it now appears that the Skamarock-Klemp hydrostatic gravity wave test is broken. It crashes on the first timestep, having produced a ridiculous explicit forcing.

I am making a gusto issue rather than gusto_case_studies issue because this seems to be a gusto issue.

This test is a "vertical slice" but on a 3D mesh which is 1 cell wide in the y-direction, and involves the Coriolis force which should be balanced by an "extra forcing term" that we add. However, the crash seems to be caused by a stupidly large forcing in the y-direction, generating unphysically fast winds.

By process of elimination I have worked out that this wrong forcing is surprisingly coming from the pressure gradient term, and specifically this part:

jump(theta_v*w, n)*avg(exner)*dS_v_qp

I can break this into horizontal and vertical components and keep only the vertical:

        wk = dot(w, domain.k)*domain.k
        wh = w - wk

        pressure_gradient_form = ...
                + jump(theta_v*wk, n)*avg(exner)*dS_v_qp

which seems to lead to sensible behaviour.

So why is jump(theta_v*wh, n)*avg(exner)*dS_v_qp not zero? I expect n to be vertical only but it seems to have some horizontal component... Has this changed for this kind of mesh without us noticing? I can't see a Gusto change that could explain this behaviour in the last few months.

[JHopeCollins]

I expect n to be vertical only but it seems to have some horizontal component.

Maybe I have misunderstood the form here, but isn't this the expected behaviour? dS_v is integrating over vertical facets, which will only have a normal component in the horizontal direction. Is w here the velocity?

But that would mean that your modified form would always evaluate to zero, so maybe there is something else going on here.

I remember @colinjcotter saying that jump(..., n) had some unintuitive behaviour for certain cases but I forget exactly what. Maybe there's something going on there.

[tommbendall]

Ah yes, I forgot which way around dS_v and dS_h are, sorry! So yes we should only expect this to be zero if theta_v is constant (which it is).

The full form is actually

        pressure_gradient_form = pressure_gradient(subject(prognostic(
            cp*(-div(theta_v*w)*exner*dx_qp
                + jump(theta_v*w, n)*avg(exner)*dS_v_qp), 'u'), self.X))

where w is the test function in the velocity space

[JHopeCollins]

Ah, so you mean if theta_v is constant then you expect jump(theta_v*w, n) to be zero always? Because the horizontal velocity is continuous in the horizontal direction, and the vertical velocity obviously has zero component against that facet?

[tommbendall]

Yes that's right! So this is suggesting that something about jump(theta_v*w, n) isn't what I expect on the horizontal faces

[JHopeCollins]

What's up with the horizontal faces? Aren't the dS_h integrals manually elided?

[JHopeCollins]

Personal bugbear, but I have always struggled with the Firedrake style of casting facet fluxes using jump instead of writing the flux expression explicitly and then writing the contribution to each element.
I always have to translate it to a finite volume type form before I can tell if it's doing what I think it is!

In FV speak, this is a central difference flux approximation to the pressure gradient at the face right? But having the jump operator in here always makes me initially think there's some unwind-ness to it.
EDIT: oh it isn't quite a standard central difference because the potential temperature is from the local side of the facet rather than being included in the avg.