Implemented semi-lagrangian advection solver

[sdeastham]

The Eulerian advection scheme caused massive diffusion, slow simulations, and had multiple bugs (see #109 and #110). This PR replaces it with a semi-Lagrangian approach, shifting the local vector by floor(CFL) and then performing a single forward Euler step for the remaining fractional timestep. Credit to @coco-yeung for the implementation on which this is based!

This is NOT zero-diff. It dramatically reduces numerical diffusion, which has the side effect of quite severely changing the simulation; in a reference case (issl_rhi140) the lifetime is decreased (see attached). It also increases run speed by 1.2-1.5x in terms of simulation hours per real-time hour.

The reason for this decrease in lifetime seems to be that, with the high degree of numerical diffusion in the Eulerian code, the contrail core is artificially preserved:

A review of the code would be greatly appreciated (nominating @lrobion if possible, and input from @Calebsakhtar would also be greatly appreciated). I believe that this code change does resolve several outstanding issues, but it would be good to verify that the new behaviour is expected.

[speth]

A couple thoughts:

• If we adopt this, I think we need some clear documentation of what the scheme is. This becomes more important as the implementation evolves and neither the Fritz et al (2021) paper nor the Xu (2024) thesis provide a good description of how the code works any more. I'd say the minimal version would be an extended docstring for the AdvDiffSystem class, though I realize we don't even build Doxygen docs for this project.
• Could we run this before/after comparison with different grid resolution parameters as a rough convergence study? It would help to see that the two methods are approaching the same solution on a sufficiently fine grid (though perhaps the current method isn't even consistent in the numerical analysis sense, given the known issues).
• If the old scheme was to blame for "slow simulations" but this makes it slower per simulated hour, what does that mean? Of course, correctness beats speed.

[sdeastham]

@speth Thanks for this! On each:

1. A docstring makes total sense. I'll add something ASAP. Once the code is a little cleaner I would argue the entire project would benefit from improved documentation, perhaps rising to the level of a model development paper (or at least a public manual).
2. A resolution study is long, long overdue. I've done some thinking on convergence w.r.t. time step - see Process subcycling and time-averaged diffusion #126 - but I can extend that to include a grid resolution study. I was a bit reluctant to do so but I see the point that this is necessary.
3. To clarify - this scheme is faster per simulated hour than the old one, but I think the main benefit is really the suppression of numerical diffusion which was causing artificial sublimation of ice crystals.

[lrobion]

Hi Seb,
Thanks for this! As a I review this, can you point me to what paper this implementation is based on so that I have a reference. It looks a bit like a flux form semi-Lagrangian approach like Lin & Rood 1996. Is that the one or is there another paper you used?

I agree with Ray that we need new doc, and a before / after (maybe the before should have the fixes to #109 and #110) would be super helpful.

[sdeastham]

Thanks @lrobion ! The implementation is much, much simpler than that. Frankly, calling it "semi-Lagrangian" isn't strictly correct but I'm not aware of a better term. The basic idea is:

• Calculate the CFL, c.
• If c exceeds 1, then take the integer part N (e.g. if c = 11.7, take N = 11) and shift the entire vector by N.
• Perform a 1st-order Eulerian flux calculation for a timestep of dt * (c-N)/c.

This only works because the speed is always constant in a given row or column. This is more or less described by Ritchie (1986) (https://doi.org/10.1175/1520-0493(1986)114%3C0135:ETIAWT%3E2.0.CO;2). Hopefully that (somewhat) helps!

Will get started on that doc. I'm on leave next week so that might be a good time.. 😂

[lrobion]

I think there might be a small mistake in the forward Euler step because it does not always preserve monotonicity even though we use a flux limiter for that purpose. If you add this test to the test suite this will fail (the first failure being that we create an artificial negative value).

// In test_adv_diff_solver.cpp
    TEST_CASE("Semi-Lagrangian 1D Advection preserves monotonicity", "[advection]"){

        // Monotone non-decreasing profile
        const std::vector<double> initial = {0.0, 1.0, 3.0, 4.0, 5.0, 5.0};
        const double dt = 1.0;
        const double ds = 1.0;
        const double bc_left = 0.0;
        const double bc_right = 5.0;

        auto checkMonotonicity = [&](const std::vector<double>& slice){
            // 1. Values should stay in the range spanned by the initial data and the BCs.
            for (std::size_t m = 0; m < slice.size(); m++) {
                INFO("cell " << m << " = " << slice[m]);
                REQUIRE(slice[m] >= 0.0 - 1e-12);
                REQUIRE(slice[m] <= 5.0 + 1e-12);
            }

            // 2. A monotone profile must stay monotone
            for (std::size_t m = 1; m < slice.size(); m++) {
                INFO("cells " << m - 1 << ", " << m << " = " << slice[m-1] << ", " << slice[m]);
                REQUIRE(slice[m] >= slice[m-1] - 1e-12);
            }

            // 3. Total variation must not increase
            double tv_initial = 0.0;
            double tv_final = 0.0;
            for (std::size_t m = 1; m < slice.size(); m++) {
                tv_initial += std::abs(initial[m] - initial[m-1]);
                tv_final += std::abs(slice[m] - slice[m-1]);
            }
            INFO("TV before = " << tv_initial << ", TV after = " << tv_final);
            REQUIRE(tv_final <= tv_initial + 1e-12);
        };

        SECTION("Fractional CFL 0.80, above the 2/3 TVD limit"){
            std::vector<double> slice = initial;
            double velocity = 0.8; // dt = ds = 1 -> no integer shift, fractional CFL = 0.8
            semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right);
            checkMonotonicity(slice);
        }

        SECTION("Fractional CFL 0.60, below the 2/3 TVD limit"){
            std::vector<double> slice = initial;
            double velocity = 0.6; // dt = ds = 1 -> no integer shift, fractional CFL = 0.6
            semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right);
            checkMonotonicity(slice);
        }
    }

I think this is because the face flux for cell $m$ e.g. $F_{m + 1/2}$ is incorrectly reconstructed (see line 591 of AdvDiffSystem.cpp).

For a cell $m$ assuming positive wind, when we update its value $\phi_m$ with a forward Euler step we do

$\phi_m^{t+\Delta t} = \phi_m^{t} - c* (F_{m + 1/2} - F_{m - 1/2})$

which is correct ($c = u\Delta t / \Delta s = CFL)$. But $F_{m + 1/2}$ here represents the time average value of the flux between $t$ and $t + \Delta t$.

We can derive that value by considering the mass that crosses the cell face $m+1/2$ between $t$ and $t + \Delta t$. This mass is exactly $\phi(x)$ for $x \in [x_{m + 1/2 - u\Delta t}, x_{m + 1/2}]$ (mass that gets pushed through by the wind during one timestep).

$m = \int_{x_{m + 1/2 - u\Delta t}}^{x_{m + 1/2}} \phi(x)dx$

Because $c < 1$ by construction (we use forward Euler for the remaining substep) all of this mass lies within a single MUSCL reconstruction (piecewise linear reconstruction of the cell values from the cell average, by definition: $\phi(x) = \phi_m + \sigma_m \xi$ where , $\xi = (x − x_m)/ \Delta s$ and $\xi \in [−1/2, 1/2]$ and $\sigma_m$ is the slope we reconstruct (we use minmod in APCEMM).

Then to get $F_{m + 1/2} = m / \Delta t = m / (c * \Delta s)$, we use the previous equation:

$F_{m + 1/2} =1 / (c \Delta s) \int_{x_{m + 1/2 - u\Delta t}}^{x_{m + 1/2}} \phi(x)dx$

and substituting $\phi(x)$ with the cell reconstruction + changing variables to $\xi$:

$F_{m + 1/2} =1 / c \int_{1/2 - c}^{1/2} (\phi_m + \sigma_m \xi) d\xi$

$F_{m + 1/2} =1 / c * (c \phi_m + \frac{\sigma_m}{2} [\xi^2]_{1/2 -c}^{1/2}) $

$F_{m + 1/2} =\phi_m + \frac{\sigma_m}{2} * (1 - c)$

This is different from line 591 (and other places of the SL solver in AdvDiffSystem.cpp) which just do $F_{m + 1/2} =\phi_m + \frac{\sigma_m}{2}$ .

Rerunning the test with this fix makes the tests pass. I am not an expert in this so I maybe be wrong, this cropped up because I wanted to add tests for properties we know about the solver here.

[lrobion]

This looks good, I am happy to discuss any thing I've flagged. Some of the allocation stuff is just for performance and I made not sure how much of a difference it would make to hoist them out of the hot loops.

[sdeastham]

@lrobion - thanks so much for your review! I've tried to make limited-scope edits which I think address the issues you raised, but would value your thoughts.

@speth - there's now a docstring, but the goal is to generate a much more thorough set of documentation which will try to cover recent major changes (including this one). I'll get that resolution study together too, and will make that available for you to review before going any further with the PR.

[lrobion]

From my end I think we resolved everything for the PR.
Two notes:

1. I think it'd be worth adding the test I commented here as part of the normal test suite
2. optional: if we're in the business of speed there are some optimizations we can do because the hot loops do not vectorize right now (checked with -fopt-info) due to branching (ternaries and minmod). But I think transport duration is now dominated by the SOR solve so the gains are likely minor.

[sdeastham]

Fantastic - thank you so much @lrobion ! I've added your proposed test (which I'm happy to say the code passes). I'm now working on documentation. I am also performing a resolution sweep, although for that I'm using the code in #126 - it's a PR which will be on top of this one, but which has a couple of additional improvements. They're zero-diff for the standard case though. I'll put the results of the resolution testing there, and I'm happy to hold off on merging both PRs until review is complete for that.

[lrobion]

That works for me. #126 looks very promising, but I'll look into it in detail once you're done with the resolution sweep!