rtx-cfd: mask hysteresis — sticky cell classification against a reference mask

The measured FSI3 killer is the bistable mask: one geometry (|d|
identical to 4 digits) samples two load branches (60 vs 120 kN), and
the traced s=1 death is the secant walking a 20x load cliff
(68,886 -> 1,307,938 N over a 1e-3 candidate change). A sticky band
makes the load map single-valued at the crossing: formerly-fluid
cells flip only at phi < -band, formerly-solid at phi > band,
classified against the mask held at rebuild time — in a coupling
loop, the restored committed step-start mask.

Band 0 is structurally bit-identical (phi > -0.0 <=> phi > 0.0) and
verified digit-for-digit on both committed defaults (FSI2 and FSI3,
every physics digit). Measured cost on the translating-circle MMS at
band 0.25h: +0.5% field error (u/p ratios 1.17/2.12 vs the
no-hysteresis moving levels 1.16/2.11); flip delay = band/(v dt),
deterministic. Exposed as RTX_FSI{2,3}_HYST in multiples of h_min.

Co-Authored-By: Claude Fable 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded.rs

@@ -119,6 +119,8 @@ pub struct EmbeddedPisoSolver {
     body: Option<EmbeddedBody>,
     mask: Option<EmbeddedMask>,
     moving: bool,
+    /// Mask hysteresis band in multiples of the min cell size (0 = off).
+    mask_hysteresis: f64,
     time: f64,
     initialized: bool,
 }
@@ -135,11 +137,26 @@ impl EmbeddedPisoSolver {
             body: None,
             mask: None,
             moving: false,
+            mask_hysteresis: 0.0,
             time: 0.0,
             initialized: false,
         })
     }
 
+    /// Mask hysteresis for the moving-body rebuild, as a fraction of the
+    /// min cell size (default 0, exactly the plain rebuild). With a band,
+    /// a cell within `band * h_min` of the surface keeps the
+    /// classification it has in the mask held at rebuild time — in a
+    /// coupling loop that restores a [`Self::snapshot`] before each
+    /// subiteration, that is the committed step-start mask, so every pass
+    /// of a step classifies against ONE reference and candidate geometries
+    /// within the band all see the SAME mask (the pass map stops flipping
+    /// cells on sub-band candidate differences). The cost is the effective
+    /// wall lagging the true surface by up to the band.
+    pub fn set_mask_hysteresis(&mut self, band_in_h: f64) {
+        self.mask_hysteresis = band_in_h;
+    }
+
     /// Volumetric momentum source `(x, y, t) -> (f_x, f_y)` per unit volume.
     pub fn set_momentum_source<F>(&mut self, f: F)
     where
@@ -913,7 +930,16 @@ impl EmbeddedPisoSolver {
         if self.moving {
             if let Some(body) = &self.body {
                 let (nx, ny, dx, dy) = field.grid_info();
-                let new_mask = EmbeddedMask::build(body, nx, ny, dx, dy, t_new)?;
+                let new_mask = EmbeddedMask::build_with_reference(
+                    body,
+                    nx,
+                    ny,
+                    dx,
+                    dy,
+                    t_new,
+                    self.mask.as_ref(),
+                    self.mask_hysteresis * dx.min(dy),
+                )?;
                 if let Some(old_mask) = &self.mask {
                     for j in 0..ny {
                         for i in 0..nx {

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded_body.rs

@@ -442,6 +442,31 @@ impl EmbeddedMask {
         dy: f64,
         t: f64,
     ) -> CfdResult<Self> {
+        Self::build_with_reference(body, nx, ny, dx, dy, t, None, 0.0)
+    }
+
+    /// [`Self::build`] with mask hysteresis: a cell whose centre lies
+    /// within `band` (metres) of the surface keeps the classification it
+    /// has in `reference`, flipping only once `phi` crosses `band` on the
+    /// far side — a formerly-fluid cell goes solid only at `phi < -band`,
+    /// a formerly-solid cell goes fluid only at `phi > band`. This makes
+    /// the classification a single-valued function of geometry around the
+    /// reference: two candidate geometries within the band produce the
+    /// SAME mask, at the cost of the effective wall lagging the true
+    /// surface by up to `band`. With `band = 0` or no reference (or a
+    /// reference of different dimensions) this is exactly [`Self::build`].
+    #[allow(clippy::too_many_arguments)]
+    pub fn build_with_reference(
+        body: &EmbeddedBody,
+        nx: usize,
+        ny: usize,
+        dx: f64,
+        dy: f64,
+        t: f64,
+        reference: Option<&EmbeddedMask>,
+        band: f64,
+    ) -> CfdResult<Self> {
+        let sticky = reference.filter(|m| band > 0.0 && m.nx == nx && m.ny == ny);
         let xc = |i: usize| (i as f64 + 0.5) * dx;
         let yc = |j: usize| (j as f64 + 0.5) * dy;
         let mut cell_fluid = vec![true; nx * ny];
@@ -449,7 +474,12 @@ impl EmbeddedMask {
         let mut anchor = None;
         for j in 0..ny {
             for i in 0..nx {
-                let fluid = body.phi(xc(i), yc(j), t) > 0.0;
+                let phi = body.phi(xc(i), yc(j), t);
+                let fluid = match sticky.map(|m| m.cell_fluid[j * nx + i]) {
+                    Some(true) => phi > -band,
+                    Some(false) => phi > band,
+                    None => phi > 0.0,
+                };
                 cell_fluid[j * nx + i] = fluid;
                 if fluid {
                     fluid_cells += 1;
@@ -1170,6 +1200,63 @@ mod tests {
         }
     }
 
+    /// Mask hysteresis: within the band every cell keeps the reference
+    /// classification (two geometries within the band produce the SAME
+    /// mask); past the band cells flip; band 0 with a reference is exactly
+    /// the plain build.
+    #[test]
+    fn hysteresis_keeps_the_reference_classification_within_the_band() {
+        let n = 32;
+        let h = 1.0 / n as f64;
+        let band = 0.5 * h;
+        let circle_at = |cx: f64| EmbeddedBody::circle(cx, 0.5, 0.2);
+        let reference = EmbeddedMask::build(&circle_at(0.5), n, n, h, h, 0.0).unwrap();
+        let flips = |a: &EmbeddedMask, b: &EmbeddedMask| {
+            let mut count = 0;
+            for j in 0..n {
+                for i in 0..n {
+                    if a.is_fluid_cell(j, i) != b.is_fluid_cell(j, i) {
+                        count += 1;
+                    }
+                }
+            }
+            count
+        };
+
+        // A shift inside the band: the plain build flips cells, the sticky
+        // build must equal the reference cell-for-cell.
+        let shifted = circle_at(0.5 + 0.4 * h);
+        let plain = EmbeddedMask::build(&shifted, n, n, h, h, 0.0).unwrap();
+        let sticky =
+            EmbeddedMask::build_with_reference(&shifted, n, n, h, h, 0.0, Some(&reference), band)
+                .unwrap();
+        assert!(
+            flips(&plain, &reference) > 0,
+            "a 0.4h shift flips no cells — the test is vacuous"
+        );
+        assert_eq!(
+            flips(&sticky, &reference),
+            0,
+            "cells flipped inside the hysteresis band"
+        );
+
+        // A shift past the band flips cells even with hysteresis.
+        let far = circle_at(0.5 + 2.0 * h);
+        let sticky_far =
+            EmbeddedMask::build_with_reference(&far, n, n, h, h, 0.0, Some(&reference), band)
+                .unwrap();
+        assert!(
+            flips(&sticky_far, &reference) > 0,
+            "the band froze the mask against a 2h shift"
+        );
+
+        // Band 0 with a reference is exactly the plain build.
+        let zero =
+            EmbeddedMask::build_with_reference(&shifted, n, n, h, h, 0.0, Some(&reference), 0.0)
+                .unwrap();
+        assert_eq!(flips(&zero, &plain), 0);
+    }
+
     #[test]
     fn interface_velocity_interpolates_along_the_nearest_edge() {
         // Unit square, CCW; each vertex carries a distinct velocity.