Elemental Slope Limiter

The elemental slope limiter is an optional, per-node numerical-stabilization step applied to the water surface elevation immediately after the GWCE continuity solve. It watches the per-element elevation gradient and, at nodes where that gradient exceeds a user-supplied threshold, replaces the elevation with a local area-weighted patch average, clamping a runaway gradient at a handful of problem nodes without globally smoothing the solution. It is the Cocoa analogue of ADCIRC’s elemental_slope_limiter fort.13 nodal attribute.

The feature is disabled by default and costs nothing when it is not configured: nothing is allocated and the time-step pass is skipped entirely. It is available on the CG (continuous-Galerkin) solver only.

Note

This is a different feature from the wet/dry slope limiter. The wet_dry_slope_limiter configuration key (ADCIRC SLIM, internally \(\alpha_L\)) scales down elevation gradients in the GWCE stiffness matrix everywhere in shallow water (see Wetting and Drying). The elemental slope limiter documented here is a separate, selective post-solve correction keyed off a per-node threshold. The two share only the words “slope limiter”.

Overview 

After each GWCE continuity solve produces a new elevation field \(\zeta^{n+1}\), and after wet/dry status and ghost elevations have been brought up to date, the limiter runs a single element pass over the current wet-element list. For every wet element it computes the magnitude of the elemental elevation gradient. Each of the element’s nodes (three on a triangle, four on a quadrilateral) carries a threshold; a node is flagged active for this step when the gradient meets or exceeds its threshold. Active nodes then have their elevation overwritten by the wet-element-area-weighted average of the surrounding elements’ mean elevations. The pass runs before the momentum solve consumes \(\zeta\), and a ghost halo exchange refreshes the overwritten values across partition boundaries.

A small number of nodes can develop spuriously steep elevation gradients (often at complex bathymetry or wet/dry fronts) and destabilize the run. The limiter smooths only those nodes, and only on the steps where the gradient exceeds the threshold, instead of adding global damping.

When to use it 

Use the elemental slope limiter when a run is stable almost everywhere but develops a localized elevation blow-up at a few identifiable nodes. It does not make up for an inadequate time step, a poorly conditioned mesh, or an inappropriate \(\tau_0\) (see Generalized Wave Continuity Equation (GWCE)). Because it overwrites the solved elevation, a threshold set too low smooths real signal; tune it so that it engages only at the problem nodes.

Theory 

Elemental elevation gradient 

For a linear (P1) triangular element \(e\) with nodes \(i,j,k\), area \(A_e\), and the stored shape-function derivative coefficients \(b_i\) (contributing to \(\partial/\partial x\)) and \(a_i\) (contributing to \(\partial/\partial y\)), the two components of the elemental elevation gradient are

\[\frac{\partial \zeta}{\partial x}\bigg|_e = \frac{S_{mx,e}}{2 A_e} \sum_{p \in \{i,j,k\}} \zeta_p\, b_p, \qquad \frac{\partial \zeta}{\partial y}\bigg|_e = \frac{S_{my,e}}{2 A_e} \sum_{p \in \{i,j,k\}} \zeta_p\, a_p,\]

where \(S_{mx,e}\) and \(S_{my,e}\) are the per-element momentum scale factors that carry the CPP/spherical map correction (see Coordinate Systems). The gradient magnitude tested at the node is

\[\lvert \nabla \zeta \rvert_e = \sqrt{\left(\frac{\partial \zeta}{\partial x}\bigg|_e\right)^2 + \left(\frac{\partial \zeta}{\partial y}\bigg|_e\right)^2}.\]

This is the elevation-gradient construction the GWCE right-hand-side assembler uses (Generalized Wave Continuity Equation (GWCE)), including its scale-factor convention, so the two agree to floating-point tolerance.

Activation rule 

Let \(T_n\) be the threshold at node \(n\). During the element pass a node is flagged active when

\[\lvert \nabla \zeta \rvert_e \ge T_n .\]

Because a node touches several elements, it is flagged active if any incident wet element’s gradient meets the threshold.

Patch-average replacement 

Each wet element contributes its area-weighted mean elevation to each of its nodes. With the element-mean elevation

\[\bar{\zeta}_e = \frac{\zeta_i + \zeta_j + \zeta_k}{3},\]

the limiter accumulates, per node,

\[\text{patchSum}(n) = \sum_{e \ni n,\; e\ \text{wet}} A_e\, \bar{\zeta}_e ,\]

and divides by the node’s current wet area \(A^{\text{wet}}(n) = \sum_{e \ni n,\, e\ \text{wet}} A_e\) (the dynamic active_element_area). At an active node with positive wet area the elevation is replaced by

\[\zeta(n) \;\leftarrow\; \frac{\text{patchSum}(n)}{A^{\text{wet}}(n)} .\]

Nodes that are not active, or whose surrounding wet area is zero, are left unchanged. The patch sum and a per-node count of exceedances (a node is active when its count is positive) accumulate in one pass over the owned wet-element list; in a distributed run the contributions are summed into the owning ranks before the node update.

Quadrilaterals 

On a hybrid mesh the same pass runs over the wet quadrilaterals, into the same two accumulators, so a node’s patch spans all of its incident elements of either type. A quadrilateral’s gradient is the area-averaged gradient from its own four corner coefficients, by the formula above, and its mean elevation is the four-corner average. In place of \(A_e\), each corner contributes with its lumped area weight, three times \(\int \phi_i\), which is the weight active_element_area is built with (a triangle’s three times \(A_e/3\) is its area), so the division still recovers an average.

When it runs in the time step 

The limiter is applied once per time step, in this order:

GWCE continuity solve produces \(\zeta^{n+1}\).
Wet/dry classification updates element/node status, combining every element-to-node result over all ranks, so every ghost copy leaves it with its owner’s elevation and status (boundary-element gradients read consistent neighbor elevations).
Elemental slope limiter pass.
A ghost halo exchange refreshes the overwritten \(\zeta^{n+1}\).
Momentum solve consumes the limited \(\zeta^{n+1}\).

Configuration 

The limiter is controlled by a single key, physics.elemental_slope_limiter, which takes one of three forms:

Value	Meaning
absent / `null`	Disabled (the default). The attribute is never registered, nothing is allocated, and the time-step pass is skipped.
a number (e.g. `0.001`)	Enabled with a single constant threshold applied at every node.
`mesh`	Enabled with per-node thresholds read from the mesh NetCDF nodal-attribute `elemental_slope_limiter`.

The threshold has units of m/m (a dimensionless elevation slope: meters of elevation change per meter of horizontal distance).

Constant threshold 

physics:
  elemental_slope_limiter: 0.001    # constant m/m threshold at every node

Per-node thresholds from the mesh 

physics:
  elemental_slope_limiter: mesh     # per-node values from the mesh file

With mesh, the mesh NetCDF must contain a nodal-attribute variable named elemental_slope_limiter (in the /nodal_attributes group, or the root group in the legacy flat layout); a missing variable fails loud at startup. See Mesh Preparation for how nodal attributes are stored. Per-node values must be non-negative (zero is the always-active sentinel below); a negative threshold, a NaN, or an inf in the mesh field fails loud, naming the variable and the offending node count.

Disabled (default)

Omitting the key, or setting it to null, disables the feature:

physics:
  manning_n: 0.025
  tau0: 0.005
  # elemental_slope_limiter not present -> disabled

Sentinels 

The activation test is exactly \(\lvert \nabla \zeta \rvert_e \ge T_n\), and the threshold must be non-negative:

Threshold	Behavior
\(T_n > 0\)	The node is limited on any step where an incident wet element’s gradient magnitude meets or exceeds \(T_n\).
\(T_n = 0\)	Always active: a zero gradient always meets a zero threshold, so the node is limited (patch-averaged) on every step it is wet.

ADCIRC suggests 0.001 (m/m) as a starting value; the right threshold is problem-dependent and must be tuned.

Deviations from ADCIRC 

Cocoa’s elemental slope limiter matches the ADCIRC feature’s intent but differs from it in several respects. Each difference is deliberate and pinned by a test.

Negative thresholds 

ADCIRC treats a negative threshold as a diagnostic-only sentinel: it logs the offending node and leaves the elevation untouched. Cocoa has no such mode, so a negative threshold from either the config or the mesh is rejected at startup instead of being reinterpreted.

Stateless per-step activation 

ADCIRC’s documentation states that once a node is activated it “remains active for the remainder of the run.” ADCIRC’s code does not do this: it recomputes activation from scratch every step. Cocoa follows the code and recomputes activation every step. A node that limits on one step is limited again on a later step only if its gradient still meets the threshold then; when the gradient subsides, the node stops being limited. Only ADCIRC’s logging flag is sticky (“ever activated”), which Cocoa does not reproduce as a persistent field.

Dynamic area normalizer 

The patch average divides patchSum (built from the currently wet elements) by the node’s wet area. Cocoa uses the dynamic active_element_area, recomputed every step from the current wet set, so the numerator and denominator always share the same set of elements. ADCIRC instead divides by a static TotalArea computed once at cold start from the initial wet/dry state. On a node whose surrounding elements have partially dried since cold start, ADCIRC’s numerator and denominator disagree and bias the average. Cocoa’s patch average can therefore differ from ADCIRC’s on drying patches.

Finite-amplitude assumption 

ADCIRC multiplies the elevation contribution by its finite-amplitude flag IFNLFA (0 or 1). Cocoa runs finite-amplitude nonlinearity as its standard path and treats this factor as 1; there is no linear-elevation mode.

No reporting or hotstart persistence 

ADCIRC writes an ESLNodes.63 file recording every ever-activated node and reloads it on hotstart to resume limited-node state. Cocoa reports nothing at run time: there is no activation count on the screen log, no per-node output field, and no activation state in a checkpoint. Because activation is stateless (above), there is no run-spanning state to persist.