======================== 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 :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. The two share only the words "slope limiter". .. 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 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 :math:`\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 :math:`\tau_0` (see :doc:`../theory/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 :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 elevation-gradient construction the GWCE right-hand-side assembler uses (:doc:`../theory/gwce`), including its 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. 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 :math:`A_e`, each corner contributes with its lumped area weight, three times :math:`\int \phi_i`, which is the weight ``active_element_area`` is built with (a triangle's three times :math:`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: 1. GWCE continuity solve produces :math:`\zeta^{n+1}`. 2. 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). 3. Elemental slope limiter pass. 4. A ghost halo exchange refreshes the overwritten :math:`\zeta^{n+1}`. 5. Momentum solve consumes the limited :math:`\zeta^{n+1}`. Configuration ------------- The limiter is controlled by a single key, ``physics.elemental_slope_limiter``, which takes one of three forms: .. 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` - 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. 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. 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.