A quasi-monotone flux limiter for isopycnal diffusion (Redi)

[mark-petersen]

These improvements were originally written by @dengwirda, as described in E3SM-Ocean-Discussion#53. Please refer to that PR for additional test results and discussion.

This PR addresses a range of issues concerning the isopycnal diffusion operator (Redi):

1. Updating the use of the triad slopes in one of the Redi fluxes.
2. Correcting the implementation of the Redi fluxes.
3. Adding a new new quasi-monotone flux limiting strategy that attempts to prevent the development of non-monotone tracer values.

[non-BFB]
[CC]

[mark-petersen]

@dengwirda thank you for your detailed investigation and testing that resulted in the PR. Feel free to make further comments here, as time permits with your new schedule.

[mark-petersen]

Copying over documentation and testing from E3SM-Ocean-Discussion#53 so that this PR is complete.

Formulation

Isopycnal diffusion of a tracer $\psi$ is achieved via a rotated diffusion operator

$\partial_{t} \psi = \dots + \nabla \cdot \big( \kappa \mathbf{K} \nabla \psi \big)$

where $\kappa$ is the diffusivity and $\mathbf{K}$ is a tensor built from the isopycnal slopes $S$. Integrating over layers via finite-volumes

$\partial_{t} h \psi = \dots + \nabla \cdot (h F_{x}) + \partial_{z} F_{z}$

$F_{x} = \kappa \nabla \psi + \kappa S \partial_{z} \psi, \quad (\text{horizontal flux})$

$F_{z} = \kappa S \cdot \nabla \psi + \kappa |S|^2 \partial_{z} \psi, \quad (\text{vertical flux})$

which leads to the 'four-term' expression that is implemented in MPAS-O

$\partial_{t} h \psi = \dots + \nabla \cdot (\kappa h \nabla \psi) + \nabla \cdot (\kappa h S \partial_{z} \psi) + \partial_{z} (\kappa S \cdot \nabla \psi) + \partial_{z} (\kappa |S|^2 \partial_{z} \psi).$

These four terms are denoted $\tau_{1,2,3,4}$, with the first three evaluated as part of the horizontal operator and the fourth as part of the implicit vertical mixing.

Evaluating $\tau_{3}$

Evaluating the cross-term fluxes $\tau_{2,3}$ is tricky, and a 'slope-triad' formulation is often used to deal with interactions between the isopycnal slopes and the (active) tracer gradients. MPAS-O uses the triad formulation currently.

This PR updates the way the triad slopes are used to evaluate the $\tau_{3}$ term, which needs to be computed at layer interfaces. The new code averages $S \cdot \nabla \psi$ (with up-down triads) using the two layers either side of an interface, compared to the one-sided approach used previously. Centering the evaluation at the layer interface appears to improve the accuracy of the scheme, helping to maintain the shape of tracer profiles being diffused along isopycnals steeply-sloped with respect to the model grid.

Remapping $\tau_{3}$

The $\tau_{3}$ term is evaluated first on edges, and then remapped onto cells. This PR updates the remapping weights, which previously included an extra factor of 2.0. This was leading to serious errors which would generate spurious 'stripe' instabilities, and contribute a kind of 'anti-diffusion' effect.

Flux limiter

It's well known that the cross-terms $\tau_{2,3}$ can generate non-monotone behaviour under certain conditions, so this PR also adds a cross-term flux limiter that attempts to prevent the development of non-monotone tracer values. Compared to the existing approach that enforces a global lower bound (typically zero) as well as allowing diffusion to be turned off in the upper layers of the ocean (config_Redi_min_layers_diag_terms), the new approach assembles a min/max bound for every cell in the ocean and deactivates the (potentially non-monotone) cross-term fluxes if a tracer value would violate these bounds otherwise.

The limiter bounds are taken as the min/max tracer values over the stencil of the operator --- all cells that are edge neighbours of a given cell over the upper, self, and lower layers.

The limiter may be configured using the new config_Redi_quasi_monotone_safety_factor option, which enables cross-term fluxes to be ramped toward zero proportional to the size of the monotonicity violation, with smaller safety_factors $\gamma$ ramping more aggressively.

$\psi^{*} = \psi^{n} + \frac{\Delta t}{h} (\tau_1 + \tau_2 + \tau_3) \quad \text{(predictor)}$

$\beta = \gamma (1 - \frac{|\psi^{*} - \psi_{bnds}|}{|\Delta \psi|})$

$\tau_{2}* = \beta \tau_{2}, \tau_{3}* = \beta \tau_{3} \quad \text{(cross-term scaling)}$

$\psi^{n+1} = \psi^{n} + \frac{\Delta t}{h} (\tau_1 + \tau_{2}* + \tau_{3}*) \quad \text{(corrector)}$

A hard on/off limiter may be obtained with config_Redi_quasi_monotone_safety_factor = 0.0.

Global spin-up

Using the new approach, a 1-yr MPAS-O QU30km standalone spin-up shows effectively no departure from initial temperature and salinity bounds, with e.g. the maximum SST remaining > -2 deg and <= 31 deg throughout. This run used the new config_Redi_use_quasi_monotone_limiter = .true. and did not disable any surface layers, with config_Redi_min_layers_diag_terms = 0. Previous spin-ups have shown a strong dependence on the number of surface layers disabled, with runs often developing SSTs in excess of 70 deg, exhibiting spurious temperature hot-spots in shelf regions, and ultimately crashing due to strong T/S gradients.

While the previous method + analysis has focused on surface effects / layers, non-monotone T/S behaviour should be prevented throughout the ocean in general. Disabling surface layers may also be undesirable in regions where isopycnals outcrop.

QU30km 1-yr standalone-spinup SST:

Idealised tests

Behaviour in the idealised 2d test case has been significantly improved, with arbitrary, sharp tracer distributions now being diffused smoothly along steeply tilted isopycnals for both linear and nonlinear equations-of-state (EoS = jm is used here, with Redi applied to T/S + passive tracers):

This compares to previous behaviour, in which issues with the $\tau_{3}$ flux led to instabilities:

Even with the fixes to $\tau_{3}$ it appears the flux limiter is still needed to prevent spurious min/max violations and instabilities. The updated code with config_Redi_use_quasi_monotone_limiter = .false. showing instabilities in T/S:

[dengwirda]

Thanks @mark-petersen, @vanroekel and all for your testing and work on this --- please let me know if any questions come up as this goes through.

@jonbob and all re: E3SM defaults:

• I'm hoping this PR should eliminate need for config_Redi_min_layers_diag_terms in any config., with config_Redi_use_quasi_monotone_limiter = .true. used in its place. This is expected to better control spurious tracer extrema while allowing Redi to be active throughout the full column.
• Additional tuning on config_Redi_quasi_monotone_safety_factor could be worthwhile though. config_Redi_quasi_monotone_safety_factor = 0.0 will do the maximal amount of limiting and config_Redi_quasi_monotone_safety_factor = 1.0 the minimum. I expect one may begin to flirt with non-monotonicity as config_Redi_quasi_monotone_safety_factor ==> 1.0 in some cases, so perhaps a safe value may be config_Redi_quasi_monotone_safety_factor = 0.8 or so? Initial standalone + G-case testing has not revealed strong sensitivity to this setting re: climate metrics.

[mark-petersen]

Testing three PRs together: #5946, #5947, #5945 (this one). I merged them as follows:

git reset --hard origin/master
git merge vanroekel/vanroekel/ocean/fix-richardson-number-calculation #5946 
git merge vanroekel/vanroekel/ocean/change-redi-slopes-and-limiter #5947
git merge -s recursive -X theirs mark-petersen/ocean/redi-flux-limiter #5945

I then had to remove variables kappaRediEdgeVal and RediKappaScaling, which were removed in #5947 but had a few strays come back in #5945 (this one).

Specifically, those changes are as follows:

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
index 9974a1386c..f58851168f 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
@@ -310,7 +310,7 @@ contains
                if (.not.use_Redi_diag_terms) cycle

                ! term 2: div( h S d/dz psi )
-               flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* &
+               flux_term2 = coef*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* &
                             0.25_RKIND* &
                             (slopeTriadDown(k, 1, iEdge)*(tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1)) &
                              /(zMid(k, cell1) - zMid(k + 1, cell1)) &
@@ -348,7 +348,7 @@ contains
                   if (.not.use_Redi_diag_terms) cycle

                   ! term 2: div( h S d/dz psi )
-                  flux_term2 = coef*RediKappa(iEdge)*kappaRediEdgeVal*layerThickEdgeMean(k, iEdge)* &
+                  flux_term2 = coef*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* &
                                0.25_RKIND* &
                                (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) &
                                 /(zMid(k - 1, cell1) - zMid(k, cell1)) &
@@ -428,9 +428,7 @@ contains
             do k = minLevelCell(iCell), maxLevelCell(iCell)
                do iTr = 1, nTracers
                   flux_term3 = rediKappaCell(iCell) * ( &
-                     RediKappaScaling(k, iCell)* &
                      redi_term3_topOfCell(iTr, k, iCell) - &
-                     RediKappaScaling(k + 1, iCell)* &
                      redi_term3_topOfCell(iTr, k + 1, iCell))

                   redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3
@@ -530,9 +528,7 @@ contains

                   tend(iTr, k, iCell) = tend(iTr, k, iCell) + &
                      rediKappaCell(iCell) * ( &
-                        RediKappaScaling(k, iCell)* &
                         scalUp*redi_term3_topOfCell(iTr, k, iCell) - &
-                        RediKappaScaling(k + 1, iCell)* &
                         scalDn*redi_term3_topOfCell(iTr, k + 1, iCell))
                end do
             end do

Alternatively, you could do this:

git reset --hard origin/master
git merge vanroekel/vanroekel/ocean/fix-richardson-number-calculation #5946 
git merge vanroekel/vanroekel/ocean/change-redi-slopes-and-limiter #5947
git merge mark-petersen/ocean/redi-flux-limiter #5945

and then choose the second block of the conflicted code for all cases, and remove kappaRediEdgeVal and RediKappaScaling lines as shown above. Then complete the merge with

git add components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
git merge --continue

With these changes, the three-PR merge passes the nighty suite with gnu and intel, both debug and optimized on chrysalis.

[dengwirda]

And config_Redi_quasi_monotone_safety_factor can be set as above @mark-petersen.

[mark-petersen]

Thanks. We are currently testing a G case with config_Redi_quasi_monotone_safety_factor=0.9, so I set that default here so we don't forget. This PR passes the nightly suite with gnu debug using config_Redi_use_quasi_monotone_limiter=true.

[vanroekel]

Approved by visual inspection and extensive testing in idealized, G-, and B-cases

[jonbob]

merged to next

[jonbob]

Re-merged to next after fix for threading issues

[jonbob]

merged to master