Elemental Slope Limiter

The elemental slope limiter is an optional, per-node numerical-stabilization device 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 imposes zero cost 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 not the same thing as 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. They share the word “slope limiter” and nothing else; the internal code names the two differently (ElevationSlopeLimiter versus slope_limiter) to avoid the collision.

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 three nodes 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.

The intent is targeted numerical stabilization: a small number of nodes develop spuriously steep elevation gradients (often at complex bathymetry or wet/dry fronts) and can destabilize the run. Rather than adding global damping, the limiter smooths only those nodes, and only on the steps where the gradient actually exceeds the threshold.

When to use it 

Reach for the elemental slope limiter when a run is stable almost everywhere but develops a localized elevation blow-up at a few identifiable nodes. It is a targeted tool, not a substitute for an adequate time step, a well-conditioned mesh, or an appropriate \(\tau_0\) (see Generalized Wave Continuity Equation (GWCE)). Because it overwrites the solved elevation, an over-eager threshold silently smooths real signal; tune the threshold so that it engages only at the genuinely pathological 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 same elevation-gradient construction the GWCE right-hand-side assembler already uses (Generalized Wave Continuity Equation (GWCE)); the limiter mirrors it exactly rather than re-deriving the 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. Both the accumulation and the count are element-parallel atomic scatters fused into one kernel pass over the owned wet-element list; in a distributed run the scattered contributions are summed into the owning ranks before the node update.

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.
Ghost elevations and wet/dry status are exchanged (so boundary-element gradients read consistent neighbor elevations).
Elemental slope limiter pass (this feature).
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. It is tri-state:

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\)	Normal limiting. 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.

A reasonable starting value is 0.001 (m/m), the value ADCIRC suggests. This is only a starting point: the appropriate threshold is problem-dependent and must be tuned.

Reporting 

The limiter reports nothing at run time: there is no activation count on the screen log, no per-node output field, and no persistence of activation state across a hotstart (see Deviations from ADCIRC).

Deviations from ADCIRC 

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

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 rather than silently 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, so activation is recomputed 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 are consistent by construction. Cocoa’s patch average therefore can differ from ADCIRC’s on drying patches, a deliberate, physically consistent choice.

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 per-node output 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 does not produce an equivalent output field and does not persist activation state across a hotstart. Instead it logs a per-interval count of limited nodes (see Reporting). Because activation is stateless (above), there is no run-spanning state to persist.

Condensed-node (VEW) skip omitted 

ADCIRC skips condensed (vertical element wall / VEW) nodes during the gradient check. That feature is not present on this branch, so the skip is omitted; every wet element participates.