Implement communicationless water borrowing limiters

[Zhengyu-Huang]

Let $q$ denote the tracer we want to remain between $q_{\min}$ and $q_{\max}$. In our implementation, we have prognostic variables: density $\rho$, tracer mass $Q = \rho q$, and we have quadrature weight $\omega$ at each quadrature point.

The basic idea of the limiter_clip_and_sum subroutine is to update $q$ to $\tilde{q}$ in the constraint set $[q_{\min}, q_{\max}]$, such that $\Vert \tilde{q}\rho w - Q w\Vert_1$ is minimal.

There are two possible implementations, which are listed below. A. The first one applies limiter_clip_and_sum for all quadrature points, but requires communication; B. The second one loops each element and applies limiter_clip_and_sum for quadrature points in the element, it does not require communication, but might fail. (This is the implementation used in E3SM).

Implementation A

Click for details. Not implementing due to communication requirement

• Compute $q = Q/\rho$
• Loop each quadrature point, do clipping and compute added mass.
  • $addmass = 0$
  • Loop each quadrature point
    • if $q > q_{\max}$ set $q = q_{\max}$ and compute $addmass += (q - q_{\max}) \rho w$;
    • if $q < q_{\min}$ set $q = q_{\min}$ and compute $addmass += (q - q_{\min}) \rho w$
• Now all q are between $[q_{\min}, q_{\max}]$, and we need to handle the added mass $addmass$. If $addmass > 0$, we need to add mass, otherwise, we need to remove mass.
  • if $addmass > 0$, $dq = q_{\max} - q$ otherwise $dq = q - q_{\min}$ (We know here $dq \geq 0$)
  • $denominator = \sum dq \rho w$
  • if $denominator > 0$, $q = q + \frac{addmass}{denominator} dq$ (When q is clipped before $dq = 0$, and updated $q$ are between $[q_{\min}, q_{\max}]$.)
• Update $Q = \rho q$

Implementation B

It is basically applying algorithm A in each element, but $q_{\max}$ and $q_{min}$ need to be adjusted, and hence here the summation $\sum$ is for the summation in the element

• In each element
  • Compute $q = Q/\rho$
  • Adjust $q_{\max}$ and $q_{min}$ in this element (The enforcement of the constraint is soft)
    • if $\sum q\rho w > q_{\max} \sum q\rho w$, $q_{\max} = \frac{\sum q\rho w } {\sum \rho w}$
    • if $\sum q\rho w < q_{\min} \sum q\rho w$, $q_{\min} = \frac{\sum q\rho w } {\sum \rho w}$
  • Loop each quadrature point, do clipping and compute added mass.
    • $addmass = 0$
    • Loop each quadrature point
      • if $q > q_{\max}$ set $q = q_{\max}$ and compute $addmass += (q - q_{\max}) \rho w$;
      • if $q < q_{\min}$ set $q = q_{\min}$ and compute $addmass += (q - q_{\min}) \rho w$
  • Now all q are between $[q_{\min}, q_{\max}]$, and we need to handle the added mass $addmass$. If $addmass > 0$, we need to add mass, otherwise, we need to remove mass.
    • if $addmass > 0$, $dq = q_{\max} - q$ otherwise $dq = q - q_{\min}$ (We know here $dq \geq 0$)
    • $denominator = \sum dq \rho w$
    • if $denominator > 0$, $q = q + \frac{addmass}{denominator} dq$ (When q is clipped before $dq = 0$, and updated $q$ are between $[q_{\min}, q_{\max}]$.)
  • Update $Q = \rho q$

This function can be called after each RK stage (E3SM implementation), or after each time step

[Zhengyu-Huang]

@valeriabarra @akshaysridhar @simonbyrne @charleskawczynski When you have time could you review this issue? Let me know whether the clipping method makes sense for you.

[valeriabarra]

Thanks Daniel. I think in the very first paragraph, you meant $Q$, tracer mass. Correct? Not tracer only

[OsKnoth]

The quadrature weights should also contain the Jacobian |dXdx| including the vertical distance.

[simonbyrne]

@valeriabarra can you review this and see how it differs from what we currently have in the limiters?

[simonbyrne]

@Zhengyu-Huang how does this integrate with the timestepper? i.e. at what point do they apply these?

[charleskawczynski]

For implementation A, what is the sum over for dqρw?

[Zhengyu-Huang]

Thanks Daniel. I think in the very first paragraph, you meant Q, tracer mass. Correct? Not tracer only

Yes, thanks!

[Zhengyu-Huang]

The quadrature weights should also contain the Jacobian |dXdx| including the vertical distance.

Yes, I agree! I think the algorithm preserves the total mass, and tracer mass, so it should have vertical distance in it

[Zhengyu-Huang]

@Zhengyu-Huang how does this integrate with the timestepper? i.e. at what point do they apply these?

They apply it after Euler step, I think it corresponds to each RK stage

[Zhengyu-Huang]

For implementation A, what is the sum over for dqρw?

It it the sum of dq*ρ*w, so that the same of q*ρ*w will be addmass in the last step

[valeriabarra]

@valeriabarra can you review this and see how it differs from what we currently have in the limiters?

At a very first glance, a simple difference I can spot is that the existing limiter uses the least square norm to determine the constraints, whereas this one uses the L-1 norm, so it seems to be less accurate in that sense. I think this is what they mean in HOMME where in this comment they say "its solution quality is not established".

I'll read more.

[OsKnoth]

In homme the method is used in connection with a special SSP-RK method (the classical Shu-Osher method). For this method the SSP form has one Euler step in each RK-stage. In case of a general explicit SSP-RK method the limiter should be applied for each tendency computation with the maximal Euler time step for this tendency in the SSP form, so that it can be applied safely in further stages. The norm for the distance to the intermediate solution should not play a role for the positivity or monotonicity.

[Zhengyu-Huang]

Mark Taylor also mentioned, E3SM no longer use SE transportation, finite volumn + 3D Semi lagrangian. But they used this clipping method before

[Zhengyu-Huang]

@Zhengyu-Huang how does this integrate with the timestepper? i.e. at what point do they apply these?

If we use ARK scheme, I think we can apply it after each time step, since the clipping is applied on the state instead of tendency

[valeriabarra]

@valeriabarra can you review this and see how it differs from what we currently have in the limiters?

This was put on-hold at the Atmos meeting, so I did not get to this until now.

The structure of the limiter_clip_and_sum algorithm is very similar to the limiter_optim_iter_full one.

• The calculation of the variable addmass (which we translated in Δtracer_mass) is the same
• The main difference stems in what we do in the case that this addmass is > or < 0:
  • Instead of working on the average mass (weightssum), the limiter_clip_and_sum algorithm seems to be working on the tracer mass ρq at the quadrature points but with this difference according to the sign

I jotted down some code/prototype in this branch valeria/clip-and-sum-limiter (which hopefully is easier to follow than this short summary above, and can also be taken for reference for the implementation)

[simonbyrne]

Following discussion with @valeriabarra:

1. From what we can tell, E3SM only rearranges mass within an element.
2. The computation of $q_\min$ and $q_\max$ remains the same (i.e. can be calculated either from neighboring element limits, or simply use $q_\min = 0, q_\max = \infty$ for basic postivity enforcement.

Both limiters are based on a basic thresholding of the tracer concentration, and computing the mass to be redistributed:

$$ \tilde{q} = \text{clamp}(q, q_\min, q_\max) $$

$$ \Delta \rho q = \sum \rho (q - \tilde q) W J $$

For what follows below, we only consider the case $\Delta \rho q > 0$ (the same applies with $q_\min$ if $\Delta \rho q < 0$):

Our current limiter then uses an iterative scheme: compute the "available total mass":

$$ \Delta \rho^* = \sum_{q < q_\max} \rho W J $$

and redistributes the tracer according to that

$$ q^* = \begin{cases} \tilde{q} + \frac{\Delta \rho q}{\Delta \rho^*} \text{if } q < q_\max\ \tilde{q} \end{cases} $$

This may not obey the threshold, so it may require multiple iterations. This I think minimizes the L2 error.

The clip-and-limit approach instead simply computes the "available tracer mass"

$$ \Delta \rho q^* = \sum \rho (q_\max - \tilde{q}) WJ $$

and then distributes the required mass proportionally:

$$ q^* = \tilde{q} + \frac{\Delta \rho q}{\Delta \rho q^*} (q_\max - \tilde{q}) $$

[charleskawczynski]

Jotting some notes down here:

From the paper that @OsKnoth shared:

Appendix B: The MassBorrow fixer

The mass borrower borrows tracer mass from an adjacent layer. It conserves the mass and can avoid negative tracers. At level k, it will first borrow the mass from the layer k+1 (lower level). If the mass is not sufficient in layer k + 1, it will borrow mass from layer k + 2. The borrower will proceed with this process until the bottom layer is reached. If the tracer mass in the bottom layer goes negative, it will repeat the process from the bottom to the top. The mass borrower works for any shape of mass profiles, as long as the column-integrated mass is positive.

The code is adapted from the tracer mass borrower implemented in the global aerosol-climate model ECHAM–HAM2 (Stier et al., 2005; Zhang et al., 2012).

is looks like it shares between potentially multiple layers. Also, based on figure 1, it looks like it's called multiple time.

[charleskawczynski]

The aerosol-climate model ECHAM5-HAM: https://acp.copernicus.org/articles/5/1125/2005/ doc: https://acp.copernicus.org/articles/5/1125/2005/acp-5-1125-2005.pdf

The global aerosol-climate model ECHAM-HAM, version 2: sensitivity to improvements in process representations: https://acp.copernicus.org/articles/12/8911/2012/, doc: https://acp.copernicus.org/articles/12/8911/2012/acp-12-8911-2012.pdf