======================== 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 :math:`\alpha_L`) scales down elevation gradients in the GWCE stiffness matrix everywhere in shallow water (see :doc:`../theory/wetdry`). 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. .. contents:: On This Page :local: :depth: 2 Overview -------- After each GWCE continuity solve produces a new elevation field :math:`\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 :math:`\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 :math:`\tau_0` (see :doc:`../theory/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 :math:`e` with nodes :math:`i,j,k`, area :math:`A_e`, and the stored shape-function derivative coefficients :math:`b_i` (contributing to :math:`\partial/\partial x`) and :math:`a_i` (contributing to :math:`\partial/\partial y`), the two components of the elemental elevation gradient are .. math:: \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 :math:`S_{mx,e}` and :math:`S_{my,e}` are the per-element momentum scale factors that carry the CPP/spherical map correction (see :doc:`../theory/coordinate_systems`). The gradient magnitude tested at the node is .. math:: \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 (:doc:`../theory/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 :math:`T_n` be the threshold at node :math:`n`. During the element pass a node is flagged active when .. math:: \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 .. math:: \bar{\zeta}_e = \frac{\zeta_i + \zeta_j + \zeta_k}{3}, the limiter accumulates, per node, .. math:: \text{patchSum}(n) = \sum_{e \ni n,\; e\ \text{wet}} A_e\, \bar{\zeta}_e , and divides by the node's current wet area :math:`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 .. math:: \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: 1. GWCE continuity solve produces :math:`\zeta^{n+1}`. 2. Wet/dry classification updates element/node status. 3. Ghost elevations and wet/dry status are exchanged (so boundary-element gradients read consistent neighbor elevations). 4. **Elemental slope limiter pass** (this feature). 5. A ghost halo exchange refreshes the overwritten :math:`\zeta^{n+1}`. 6. Momentum solve consumes the limited :math:`\zeta^{n+1}`. Configuration ------------- The limiter is controlled by a single key, ``physics.elemental_slope_limiter``. It is **tri-state**: .. list-table:: :header-rows: 1 :widths: 25 75 * - 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 ^^^^^^^^^^^^^^^^^^ .. code-block:: yaml physics: elemental_slope_limiter: 0.001 # constant m/m threshold at every node Per-node thresholds from the mesh ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: yaml 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 :doc:`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: .. code-block:: yaml physics: manning_n: 0.025 tau0: 0.005 # elemental_slope_limiter not present -> disabled Sentinels --------- The activation test is exactly :math:`\lvert \nabla \zeta \rvert_e \ge T_n`, and the threshold must be non-negative: .. list-table:: :header-rows: 1 :widths: 18 82 * - Threshold - Behavior * - :math:`T_n > 0` - **Normal limiting.** The node is limited on any step where an incident wet element's gradient magnitude meets or exceeds :math:`T_n`. * - :math:`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. See also -------- - :doc:`../theory/gwce`: the GWCE continuity solve the limiter post-processes, and the elevation-gradient construction it mirrors. - :doc:`../theory/wetdry`: the *distinct* wet/dry :math:`\alpha_L` slope limiter (``wet_dry_slope_limiter``) and the wet-element area bookkeeping. - :doc:`mesh_preparation`: storing per-node values as mesh nodal attributes.