Advection in `ShallowWaterModel`

[francispoulin]

Inspired by trying to resolve #1866, I have been looking over shallow_water_advection_operators.jl.

The first cocern is with line 22. . We are trying to discretize a term the term

$\partial_x (h u u)$

Using _advective_momentum_flux_Uu is a good idea, since then we can use a variety of advections schemes. However, we are advecting solution[1], $uh$, with solution[1], $uh$, and then dividing by $h$. I fear that the dividing by $h$ afterward is going to eliminate any higher order that we might achieve by say WENO5. Anyone know whether this is true?

A few ideas come to mind of how we could do things differently.

1. We could advect solution[1]/solution[3], $uh/h$, by solution[1], $uh$.
2. We could advect solution[1], $uh$, by solution[1]/solution[3], $uh/h$.
3. We could store $u$ and $v$ as solution[4], and solution[5] and then change either of the solution[1]/solution[3] above to `solution[4].

I don't know which is better but is this something we should consider?

Thoughts @glwagner @simone-silvestri ?

[glwagner]

I suggest creating a model property called "velocities" and storing objects that represent the velocity components there. For ConservativeFormulation these will be AbstractOperations.

[glwagner]

Should we discuss this on #1866?

[simone-silvestri]

I was just now trying to assess the impact of the different formulations and advection scheme orders in the shallow_water_bickley_jet example. I ll let you know the results, but it seems like the conservative formulation as it is implemented right now is quite nice. The problem with dividing the advecting velocities before it's that they are at a different location, by reconstructing the flux before and then diving you preserve location for at least two fluxes

[simone-silvestri]

By the way, WENO (and other advection schemes) are always going to be at most second order as they are implemented right now. So the order might not be the correct metric to look at the performance, the truncation error is probably what we want to look at

[francispoulin]

Thanks for the quick replies!

Sorry if I should have continued in #1866. The code has changed a lot since then so I was unsure as to how to start. I am happy to shift there if people like?

I like the idea of creating properties called "velocities", as they are used a lot for diagnostics and probably other things. I guess the drawback is more storage?

Probably a good idea to wait for @simone-silvestri 's assessement before we do anything. Good to know that the conservative scheme looks good so far. Also good to know that we can't do better than second order. I would be interested to know if there is a better approach in terms of stability and possibily accuracy.

[simone-silvestri]

I am actually looking into implementing a multidimensional sweep approach that will allow us to break the 2nd order limit

https://d-nb.info/1124132775/34 https://www.sciencedirect.com/science/article/pii/S0021999104002281

I ll see if it is too difficult or unfeasible in terms of performance

[simone-silvestri]

Maybe just as a preliminary example, the vortex merger experiment (at N =400) from https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2021MS002663

Conservative formulation WENO(order=5)

Vector Invariant Formulation WENO(order=5, vector_invariant=VorticityStencil())

Vector Invariant Formulation WENO(order=5, vector_invariant=VelocityStencil())

So the conservative formulation does not seem to be bad at all (I ll be more quantitative with the growth rate of the bickley jet). There is a problem near the boundaries but it is just because of the interpolation of h which will be solved with the new operators

[francispoulin]

Those look great! I agree that they are very similar, yet all a bit different. What resolution are you using? It would be interesting to compare with the results in the paper.

Also, just to point out that the height deviations in this example is 20% of the whole domain, which is pretty big. This is a good example that should illustrate problems with significant changes in depth.

[simone-silvestri]

For a single gaussian vortex in geostrophic balance integrated to t=10 the results are the following (N^2 = 64^2):

╭──────────────┬───────┬─────────────┬───────────────────┬─────────────┬───────────────────╮
│  Formulation │ Order │ L₂ error(h) │ L₂ error(u or uh) │ L∞ error(h) │ L∞ error(u or uh) │
├──────────────┼───────┼─────────────┼───────────────────┼─────────────┼───────────────────┤
│ Conservative │     3 │    1.66e-03 │          9.43e-03 │    5.05e-02 │          3.86e-01 │
│ VI_Vorticity │     3 │    6.80e-04 │          2.50e-03 │    1.63e-02 │          8.11e-02 │
│  VI_Velocity │     3 │    6.20e-04 │          2.30e-03 │    1.49e-02 │          7.33e-02 │
│ Conservative │     5 │    1.63e-03 │          9.03e-03 │    4.95e-02 │          3.73e-01 │
│ VI_Vorticity │     5 │    5.30e-04 │          1.95e-03 │    1.03e-02 │          5.63e-02 │
│  VI_Velocity │     5 │    4.93e-04 │          1.94e-03 │    9.35e-03 │          5.04e-02 │
│ Conservative │     7 │    1.63e-03 │          8.87e-03 │    4.93e-02 │          3.67e-01 │
│ VI_Vorticity │     7 │    4.99e-04 │          1.86e-03 │    9.17e-03 │          5.24e-02 │
│  VI_Velocity │     7 │    4.84e-04 │          1.86e-03 │    8.93e-03 │          4.93e-02 │
│ Conservative │     9 │    1.64e-03 │          8.63e-03 │    4.96e-02 │          3.62e-01 │
│ VI_Vorticity │     9 │    4.66e-04 │          1.73e-03 │    7.81e-03 │          4.52e-02 │
│  VI_Velocity │     9 │    4.55e-04 │          1.74e-03 │    7.57e-03 │          4.61e-02 │
╰──────────────┴───────┴─────────────┴───────────────────┴─────────────┴───────────────────╯

So, indeed, conservative formulation is not as good as the vector invariant formulation (at least in this case), especially when looking at momentum.

@francispoulin the resolution of the figures above is 400^2

[francispoulin]

Excellent to see the results.

I agree that the conservative scheme is not getting any better with increasing the order of accuracy.

The vector invariant scheme is slightly better, but there is not a bit difference when we increase the order of accuracy. I don't suppose it's the temporal error that is drowning this out? Maybe it's the pressure terms, which are of low order?

[glwagner]

Isn't that because the conservative implementation of advection is incorrect as pointed out on #1866 ? Or is there something else?

[francispoulin]

From looking at #1866 , there are a few things that were pointed out. I'll try and summarize them here.

1. Advection operators don't look right, in particular div_hUu because we don't multiply the advective fluxes by the cell area.
2. Tracer flux is not correct, because of missing area's, see div_Uc.
3. Functions that calculate velocities are not correct because the height is not averaged to the correct locations.
4. Momentum flux operators aren't right because they don't produce the correct high-order advection stencil, example momentum_flux_huu.

Current state fo affairs:

1 and 2. We now divide by area's insead fo volumes but that doesn't change anything sinze shallow water has \Delta z = 1. However, I don't believe that we need the lengths because they are already contained on _advective_momentum_flux_Uu and the like. For example the following line shows it for a centered scheme.

@inline advective_momentum_flux_Uu(i, j, k, grid, scheme::Centered, U, u) = Axᶜᶜᶜ(i, j, k, grid) * _symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, U) * _symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, u)

3. Inline defintions for u and v no longer appear. They were wrong before, but now they don't appear. so this problem has gone away.

4. I think this gets back to the discussion that I raised here, that we should be including the height in either the advecting or the advected quantity.

[glwagner]

Why don't we just fix the bugs?

For example the following line shows it for a centered scheme.

Can you link to the code rather than copy/pasting?

[francispoulin]

I am including the link from above here.

This shows that areas, or lenghts for shallow water, do appear but they are at the last step of the inbedded functions.

It's not clear to me that points 1 and 2 are a bug since the areas do appear in the calculations. The evidence to support that is the above plots where the conservative scheme is very similar to the vector invariant scheme.

Point 3 no long is relevant to the code.

Point 4 is not a bug but a matter of maybe reducing the accuracy to second order, which might be the highest we can do anyhow.

[simone-silvestri]

I would have rather left the volume in the flux divergence to avoid this type of confusion. In a finite volume scheme, the advection tendency is $\frac{\sum_k F_k}{V_c}$ , where there is no mention of areas because they should be embedded in the fluxes calculation.

Then the cell volume $V_c$ is a 1D metric in a 1D problem, a 2D metric in a 2D problem and a 3D metric in a 3D problem.

[francispoulin]

I see your point @simone-silvestri . I thought these changes would make things clearer but perhaps they have done the opposite. My apologies.

If you wanted to change them back to volumes, I'm perfectly happy with that.

Also, can you confirm that we don't need any extra metrics on the calculation of the momentum, as @glwagner originaly asked about it #1866 ?

[simone-silvestri]

yeah, the area is always added in the _advective_momentum_fluxs and the _advective_tracer_fluxs

[francispoulin]

Thanks for confirming!

So the two outstanding issues are:

• [ ] Do we want to divide by volumes instead of areas to be more general? (Cosmetic but maybe worth fixing)
• [ ] Figure out how the height field should be included to achieve higher accuracy, which neither formulation is able to do.
• [ ] Do we evolve the tracer or thickness weighted tracers?

[simone-silvestri]

Another outstanding issue is the formulation of tracer evolution. At the moment we are evolving tracers but maybe we want to evolve thickness weighted tracers

[francispoulin]

Thanks @simone-silvestri . I agree. I changed my list to a checked boxes, to help with monitoring the progress.

I can also change the main script to go back to volumes, if we decide this is what we perfer.

[glwagner]

Should we close this?