rtx-fsi: the noise floor interrogated — smoothing refuted, IQN-ILS lands, tight coupling reopened

The tenth session ended on "lower the interface noise floor". This
builds the levers and measures them, and the measurements overturn the
diagnosis:

- smooth_tractions: arclength moving average over the wetted surface,
  area-weighted, smooth normal-similarity factor so corners do not mix
  and the smoothed load stays continuous in the geometry. Nine unit
  tests. MEASURED NEGATIVE RESULT: the flip-scan floor is unchanged to
  0.2% at radii 1-3h — the flip's load jump is coherent through the
  fluid field (mask rebuild shifts the pressure around the flipped
  cell), and a surface average preserves coherent shifts. Default off;
  the probe pins the attribution so nobody re-reaches for this lever.

- IqnIls: interface quasi-Newton with inverse least squares (Degroote
  2009) — filtered MGS least squares over secant columns (filter
  RELATIVE to column norm), cross-step history reuse, per-step
  set_tolerance. Model-map tests: exact on anisotropic linear maps
  within dim+2 passes (scalar Aitken provably cannot be), scale
  invariant, history reuse shortens the next step, stalls at the noise
  scale instead of diverging (fixture lesson: per-pass noise, not
  state-dependent noise — the latter has a genuine fixed point).

- tests/fsi2_harness/: the FSI2 machinery extracted shared; verified
  pure code motion (committed release response reproduced to every
  printed digit). March gains RTX_FSI2_SMOOTH / RTX_FSI2_COUPLER=iqn /
  RTX_FSI2_REUSE knobs; pinned bands guard the default configuration.

- tests/fsi2_interface_noise.rs: the probe. Flip-scan floor at
  subcycle 8: 3.05e-5 (pinned); smoothing attribution pinned; the
  cross-subcycle scan recorded but unpinned (the fixed geometry
  increment's wall-velocity trend, increment/dt_c, swamps the flip
  signal at small dt_c — a dt_c^2 scaling hypothesis died in that
  operationalization). THE OPERATIONAL FLOOR — the real release step
  subiterated at tolerance 1e-9 with residuals traced — converges DEEP
  at both subcycles: s8 aitken 3.4e-9 / iqn 1.6e-9, s2 both ~6.4e-10
  in 5-6 passes. The flip jumps are events at specific geometries, not
  a floor under every step: the tenth session's subcycle-2 blowup was
  tolerance mis-budgeting (2e-4 held fixed while dt_c shrank), not an
  impassable floor. Probe bug found and fixed on the way: stale shared
  geometry leaked a 4.5e-5 phantom first residual into the first stall
  run; every measurement now resets the geometry on entry.

All 924+17 tests green: lib 44 (was 27), piston 2, curved edge 1,
FSI1, the committed FSI2 march (release response identical), the probe.

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

crates/specialized/rtx-fsi/tests/fsi2_harness/mod.rs

@@ -0,0 +1,539 @@
+//! Shared harness for the Turek–Hron FSI2 tests: the benchmark geometry,
+//! the flag mesh, the wetted-interface bookkeeping, the embedded fluid
+//! configuration, and the load sampling (spike clamp + optional surface
+//! smoothing). `turek_hron_fsi2.rs` runs the coupled march on it;
+//! `fsi2_interface_noise.rs` measures the continuity of one coupling pass
+//! on the same machinery.
+//!
+//! Everything here is code motion from the tenth-session FSI2 test —
+//! the physics and defaults are unchanged unless a test says otherwise.
+
+#![allow(dead_code)] // two test crates share this; each uses a subset
+
+use std::cell::Cell;
+use std::sync::{Arc, RwLock};
+
+use nalgebra::Vector3;
+use rtx_cfd::CfdConfig;
+use rtx_cfd::solvers::incompressible::{
+    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
+    FlowField, PoissonSolverKind, SideBoundary, polygon_interface_velocity,
+    polygon_signed_distance,
+};
+use rtx_fea::assembly::dof_mapping::DofComponent;
+use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use rtx_fsi::{FluidFace, WettedSurface, smooth_tractions};
+
+pub const L: f64 = 2.5;
+pub const H: f64 = 0.41;
+pub const RHO_F: f64 = 1000.0;
+pub const NU_F: f64 = 1e-3;
+pub const U_MEAN: f64 = 1.0;
+pub const RHO_S: f64 = 10_000.0;
+pub const E_S: f64 = 1.4e6;
+pub const NU_S: f64 = 0.4;
+
+pub const FLAG_X0: f64 = 0.25;
+pub const FLAG_X1: f64 = 0.6;
+pub const FLAG_Y0: f64 = 0.19;
+pub const FLAG_Y1: f64 = 0.21;
+
+pub fn circle_sdf(x: f64, y: f64) -> f64 {
+    ((x - 0.2).powi(2) + (y - 0.2).powi(2)).sqrt() - 0.05
+}
+
+/// The ramped parabolic inflow of the benchmark definition.
+pub fn inflow(y: f64, t: f64) -> f64 {
+    let ramp = if t < 2.0 {
+        0.5 * (1.0 - (std::f64::consts::PI * t / 2.0).cos())
+    } else {
+        1.0
+    };
+    ramp * 1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
+}
+
+pub fn env_or(name: &str, default: f64) -> f64 {
+    std::env::var(name)
+        .map(|v| v.parse().expect(name))
+        .unwrap_or(default)
+}
+
+/// The flag's Quad8 mesh (as in FSI1 and the CSM tests).
+pub fn flag_mesh(nx: usize, ny: usize) -> Mesh {
+    let mut mesh = Mesh::new(2).unwrap();
+    let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
+    let mut grid = vec![vec![None; ly]; lx];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            if i % 2 == 1 && j % 2 == 1 {
+                continue;
+            }
+            let x = FLAG_X0 + (FLAG_X1 - FLAG_X0) * i as f64 / (2 * nx) as f64;
+            let y = FLAG_Y0 + (FLAG_Y1 - FLAG_Y0) * j as f64 / (2 * ny) as f64;
+            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
+        }
+    }
+    for i in 0..nx {
+        for j in 0..ny {
+            let (a, b) = (2 * i, 2 * j);
+            let nodes = vec![
+                grid[a][b].unwrap(),
+                grid[a + 2][b].unwrap(),
+                grid[a + 2][b + 2].unwrap(),
+                grid[a][b + 2].unwrap(),
+                grid[a + 1][b].unwrap(),
+                grid[a + 2][b + 1].unwrap(),
+                grid[a + 1][b + 2].unwrap(),
+                grid[a][b + 1].unwrap(),
+            ];
+            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    mesh
+}
+
+/// The wetted-interface bookkeeping (FSI1's, plus vertex velocities).
+pub struct Interface {
+    pub wetted: Vec<NodeId>,
+    pub reference: Vec<(f64, f64)>,
+    /// Ordered boundary walk: indices into `wetted` (`usize::MAX` marks
+    /// the fixed anchor vertices inside the cylinder / at the clamp).
+    pub walk: Vec<(usize, (f64, f64))>,
+}
+
+impl Interface {
+    pub fn build(mesh: &Mesh) -> Self {
+        let eps = 1e-9;
+        let on_bottom = |p: Vector3<f64>| (p.y - FLAG_Y0).abs() < eps;
+        let on_top = |p: Vector3<f64>| (p.y - FLAG_Y1).abs() < eps;
+        let on_tip = |p: Vector3<f64>| (p.x - FLAG_X1).abs() < eps;
+        let clamped = |p: Vector3<f64>| (p.x - FLAG_X0).abs() < eps;
+
+        let mut wetted: Vec<(NodeId, (f64, f64))> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, node)| {
+                let p = node.position();
+                (on_bottom(p) || on_top(p) || on_tip(p)) && !clamped(p)
+            })
+            .map(|(&id, node)| (id, (node.position().x, node.position().y)))
+            .collect();
+        wetted.sort_by_key(|(id, _)| *id);
+        let index_of = |id: NodeId| wetted.iter().position(|(w, _)| *w == id).unwrap();
+
+        let mut bottom: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| on_bottom(n.position()) && !clamped(n.position()))
+            .map(|(&id, n)| (id, n.position().x))
+            .collect();
+        bottom.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
+        let mut tip: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| {
+                let p = n.position();
+                on_tip(p) && !on_bottom(p) && !on_top(p)
+            })
+            .map(|(&id, n)| (id, n.position().y))
+            .collect();
+        tip.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
+        let mut top: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| on_top(n.position()) && !clamped(n.position()))
+            .map(|(&id, n)| (id, n.position().x))
+            .collect();
+        top.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap());
+
+        let mut walk: Vec<(usize, (f64, f64))> = Vec::new();
+        walk.push((usize::MAX, (0.22, FLAG_Y0)));
+        walk.push((usize::MAX, (FLAG_X0, FLAG_Y0)));
+        for (id, _) in &bottom {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        for (id, _) in &tip {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        for (id, _) in &top {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        walk.push((usize::MAX, (FLAG_X0, FLAG_Y1)));
+        walk.push((usize::MAX, (0.22, FLAG_Y1)));
+
+        let reference = wetted.iter().map(|(_, p)| *p).collect();
+        Self {
+            wetted: wetted.into_iter().map(|(id, _)| id).collect(),
+            reference,
+            walk,
+        }
+    }
+
+    /// Deformed polygon vertices for the interface vector `d`.
+    pub fn polygon(&self, d: &[f64]) -> Vec<(f64, f64)> {
+        self.walk
+            .iter()
+            .map(|&(k, anchor)| {
+                if k == usize::MAX {
+                    anchor
+                } else {
+                    let (x0, y0) = self.reference[k];
+                    (x0 + d[2 * k], y0 + d[2 * k + 1])
+                }
+            })
+            .collect()
+    }
+
+    /// Per-vertex velocities for the interface velocity vector `ddot`
+    /// (anchors do not move).
+    pub fn walk_velocities(&self, ddot: &[f64]) -> Vec<(f64, f64)> {
+        self.walk
+            .iter()
+            .map(|&(k, _)| {
+                if k == usize::MAX {
+                    (0.0, 0.0)
+                } else {
+                    (ddot[2 * k], ddot[2 * k + 1])
+                }
+            })
+            .collect()
+    }
+
+    /// Deformed wetted node positions for the transfer.
+    pub fn deformed_nodes(&self, d: &[f64]) -> Vec<Vector3<f64>> {
+        self.reference
+            .iter()
+            .enumerate()
+            .map(|(k, &(x0, y0))| Vector3::new(x0 + d[2 * k], y0 + d[2 * k + 1], 0.0))
+            .collect()
+    }
+}
+
+pub fn clamp_left(mesh: &Mesh) -> BoundaryConditionSet {
+    let clamped: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| (node.position().x - FLAG_X0).abs() < 1e-9)
+        .map(|(&id, _)| id)
+        .collect();
+    let mut set = BoundaryConditionSet::new();
+    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
+        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
+            nodes: clamped.clone(),
+            components: vec![component],
+            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
+            time_range: None,
+            ramping_factor: 1.0,
+            gradual_enforcement: false,
+        }));
+    }
+    set
+}
+
+pub fn mid_amp(series: &[f64]) -> (f64, f64) {
+    let max = series.iter().copied().fold(f64::MIN, f64::max);
+    let min = series.iter().copied().fold(f64::MAX, f64::min);
+    (0.5 * (max + min), 0.5 * (max - min))
+}
+
+/// Frequency from linearly interpolated upward crossings of the mean.
+pub fn crossing_frequency(times: &[f64], series: &[f64]) -> Option<f64> {
+    let (mean, _) = mid_amp(series);
+    let mut crossings: Vec<f64> = Vec::new();
+    for k in 1..series.len() {
+        let (a, b) = (series[k - 1] - mean, series[k] - mean);
+        if a < 0.0 && b >= 0.0 {
+            crossings.push(times[k - 1] + (a / (a - b)) * (times[k] - times[k - 1]));
+        }
+    }
+    if crossings.len() < 3 {
+        return None;
+    }
+    Some((crossings.len() - 1) as f64 / (crossings.last().unwrap() - crossings.first().unwrap()))
+}
+
+/// The fluid + interface machinery every FSI2 test shares: the embedded
+/// solver configured for the benchmark channel, the deformable-geometry
+/// lock, and the load sampling with its spike clamp and (optional)
+/// surface smoothing.
+pub struct Fsi2Harness {
+    pub mesh: Mesh,
+    pub interface: Interface,
+    pub a_node: NodeId,
+    pub ny: usize,
+    pub nx: usize,
+    pub h: f64,
+    pub mu: f64,
+    pub dt_fluid: f64,
+    /// Traction smoothing radius along the surface, in metres
+    /// (0 disables). Set from `RTX_FSI2_SMOOTH` (in multiples of `h`).
+    pub smooth_radius: f64,
+    /// The deformable geometry AND its velocity, behind one lock: the
+    /// fluid's per-step mask rebuild reads the polygon; the no-slip
+    /// closure reads both.
+    pub shared: Arc<RwLock<(Vec<(f64, f64)>, Vec<(f64, f64)>)>>,
+    pub spiked_total: Cell<usize>,
+}
+
+impl Fsi2Harness {
+    /// Build the harness plus the configured solver and an at-rest field.
+    pub fn build(
+        ny: usize,
+        flag_nx: usize,
+        smooth_in_h: f64,
+    ) -> (Self, EmbeddedPisoSolver, FlowField) {
+        let h = H / ny as f64;
+        let nx = (L / h).round() as usize;
+        let mu = RHO_F * NU_F;
+        let u_peak = 1.5 * 1.5 * U_MEAN;
+        // The fluid's explicit limit; the coupling (and the flag's
+        // Newmark) run at `subcycle` fluid steps per coupled step.
+        let dt_fluid = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
+
+        let mesh = flag_mesh(flag_nx, 2);
+        let interface = Interface::build(&mesh);
+        let a_node = mesh
+            .nodes
+            .iter()
+            .find(|(_, n)| {
+                (n.position().x - 0.6).abs() < 1e-9 && (n.position().y - 0.2).abs() < 1e-9
+            })
+            .map(|(&id, _)| id)
+            .expect("point A");
+
+        let zero_d = vec![0.0; 2 * interface.wetted.len()];
+        let shared = Arc::new(RwLock::new((
+            interface.polygon(&zero_d),
+            interface.walk_velocities(&zero_d),
+        )));
+        let sdf_shared = shared.clone();
+        let vel_shared = shared.clone();
+
+        let config = CfdConfig::new()
+            .with_density(RHO_F)
+            .with_viscosity(mu)
+            .with_reference_velocity(U_MEAN)
+            .with_reference_length(0.1);
+        let params = EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-7,
+            boundaries: AleBoundaries {
+                left: SideBoundary::Velocity,
+                right: SideBoundary::PressureOutlet,
+                bottom: SideBoundary::Velocity,
+                top: SideBoundary::Velocity,
+            },
+            poisson_solver: PoissonSolverKind::Multigrid,
+            // TVD, deliberately: FSI2 marches in time and needs the
+            // shedding physics upwind's numerical viscosity killed on
+            // these grids (CFD3's finding). The limiter chatter that
+            // defeats steady fixed points (FSI1's finding) is harmless
+            // here — each step's fixed point is the interface
+            // displacement of THAT step, not a steady load.
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+        };
+        let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
+        solver.set_boundary_velocity(|x, y, t| {
+            if x <= 0.0 {
+                (inflow(y, t), 0.0)
+            } else {
+                (0.0, 0.0)
+            }
+        });
+        solver.set_moving_body(
+            EmbeddedBody::from_sdf(move |x, y, _| {
+                let geometry = sdf_shared.read().unwrap();
+                circle_sdf(x, y).min(polygon_signed_distance(&geometry.0, x, y))
+            })
+            .with_surface_velocity(move |x, y, _| {
+                let geometry = vel_shared.read().unwrap();
+                if circle_sdf(x, y) <= polygon_signed_distance(&geometry.0, x, y) {
+                    (0.0, 0.0)
+                } else {
+                    polygon_interface_velocity(&geometry.0, &geometry.1, x, y)
+                }
+            }),
+        );
+
+        // Start at rest; the ramp brings the inflow up from zero.
+        let mut field = FlowField::new(nx, ny, h, h).unwrap();
+        solver.initialize(&mut field).unwrap();
+
+        let harness = Self {
+            mesh,
+            interface,
+            a_node,
+            ny,
+            nx,
+            h,
+            mu,
+            dt_fluid,
+            smooth_radius: smooth_in_h * h,
+            shared,
+            spiked_total: Cell::new(0),
+        };
+        (harness, solver, field)
+    }
+
+    /// Publish an interface geometry (+ velocity) to the fluid.
+    pub fn set_geometry(&self, d: &[f64], ddot: &[f64]) {
+        let mut geometry = self.shared.write().unwrap();
+        geometry.0 = self.interface.polygon(d);
+        geometry.1 = self.interface.walk_velocities(ddot);
+    }
+
+    /// Surface drag and lift on cylinder + flag at the current geometry.
+    pub fn measure_force(&self, solver: &EmbeddedPisoSolver, field: &FlowField) -> (f64, f64) {
+        let mask = solver.mask().unwrap();
+        let body = solver.body().unwrap();
+        let vertices = self.shared.read().unwrap().0.clone();
+        let mut drag = 0.0;
+        let mut lift = 0.0;
+        let poly_probe = EmbeddedBody::polygon(vertices.clone());
+        for s in poly_probe.surface_samples(0.5 * self.h) {
+            if circle_sdf(s.x, s.y) < 1e-9 {
+                continue;
+            }
+            if let Some((tx, ty)) = mask.traction_at(
+                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                drag += tx * s.ds;
+                lift += ty * s.ds;
+            }
+        }
+        let circle_probe = EmbeddedBody::circle(0.2, 0.2, 0.05);
+        for s in circle_probe.surface_samples(0.5 * self.h) {
+            if polygon_signed_distance(&vertices, s.x, s.y) < 1e-9 {
+                continue;
+            }
+            if let Some((tx, ty)) = mask.traction_at(
+                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                drag += tx * s.ds;
+                lift += ty * s.ds;
+            }
+        }
+        (drag, lift)
+    }
+
+    /// Tractions on the flag's wetted surface for a given geometry, from
+    /// the solver's current field/mask; returns the transferred nodal
+    /// forces, the conservation defect, and the samples dropped (probe
+    /// failures plus spike rejections).
+    pub fn sample_load(
+        &self,
+        solver: &EmbeddedPisoSolver,
+        field: &FlowField,
+        d: &[f64],
+    ) -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
+        let vertices = self.interface.polygon(d);
+        let poly_probe = EmbeddedBody::polygon(vertices);
+        let mask = solver.mask().unwrap();
+        let body = solver.body().unwrap();
+        let mut faces = Vec::new();
+        let mut tractions: Vec<Vector3<f64>> = Vec::new();
+        let mut skipped = 0usize;
+        for s in poly_probe.surface_samples(0.5 * self.h) {
+            if circle_sdf(s.x, s.y) < 1e-9 {
+                continue; // buried in the cylinder
+            }
+            match mask.traction_at(
+                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                Some((tx, ty)) => {
+                    faces.push(FluidFace {
+                        centroid: Vector3::new(s.x, s.y, 0.0),
+                        normal: Vector3::new(s.nx, s.ny, 0.0),
+                        area: s.ds,
+                    });
+                    tractions.push(Vector3::new(tx, ty, 0.0));
+                }
+                None => skipped += 1,
+            }
+        }
+        // Spike guard: a near-degenerate reconstruction can return a
+        // finite but wild traction (the linear-fit condition sits just
+        // above its truncation threshold at concave junctions). CLAMP
+        // samples to 20x the median magnitude, keeping their direction —
+        // the physical load varies smoothly along the surface — and COUNT
+        // them: a non-zero count is a measurement of the pathology, not a
+        // silent repair. Clamping, not dropping: a hard drop threshold
+        // makes the coupling pass discontinuous in the candidate geometry
+        // (a boundary sample flips in/out of the kept set between
+        // subiterations, and the load jumps by the spike magnitude —
+        // measured as a residual bouncing at the scale of the step
+        // increment); the clamp is continuous.
+        let mut magnitudes: Vec<f64> = tractions.iter().map(nalgebra::Vector3::norm).collect();
+        magnitudes.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        let median = magnitudes.get(magnitudes.len() / 2).copied().unwrap_or(0.0);
+        if median > 0.0 {
+            let cap = 20.0 * median;
+            for traction in &mut tractions {
+                let norm = traction.norm();
+                if norm > cap {
+                    *traction *= cap / norm;
+                    self.spiked_total.set(self.spiked_total.get() + 1);
+                }
+            }
+        }
+        // Surface smoothing (after the clamp: the clamp kills the wild
+        // outliers, the smoothing spreads what remains over the stencil
+        // the cell resolution can actually support — this is the
+        // interface-noise-floor lever, measured by
+        // `fsi2_interface_noise.rs`).
+        if self.smooth_radius > 0.0 {
+            tractions = smooth_tractions(&faces, &tractions, self.smooth_radius).unwrap();
+        }
+        let nodes_now = self.interface.deformed_nodes(d);
+        let surface = WettedSurface::build(&faces, &nodes_now).expect("transfer build");
+        let nodal = surface.transfer_load(&faces, &tractions).unwrap();
+        let total_sampled: Vector3<f64> =
+            faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
+        let total_nodal: Vector3<f64> = nodal.iter().sum();
+        let conservation = (total_nodal - total_sampled).norm() / total_sampled.norm().max(1e-30);
+        (
+            self.interface
+                .wetted
+                .iter()
+                .zip(nodal)
+                .map(|(&id, f)| (id, f))
+                .collect(),
+            conservation,
+            skipped,
+        )
+    }
+
+    /// Run `subcycle` fluid substeps from the current solver/field state,
+    /// interpolating the interface geometry from `d_n` to `d_candidate`
+    /// across the substeps with the candidate's constant interface
+    /// velocity `(d_candidate - d_n) / dt` — one coupling pass's fluid
+    /// half, exactly as the coupled march runs it.
+    pub fn advance_subcycled(
+        &self,
+        solver: &mut EmbeddedPisoSolver,
+        field: &mut FlowField,
+        d_n: &[f64],
+        d_candidate: &[f64],
+        subcycle: usize,
+    ) {
+        let dt = self.dt_fluid * subcycle as f64;
+        let ddot: Vec<f64> = d_candidate
+            .iter()
+            .zip(d_n)
+            .map(|(new, old)| (new - old) / dt)
+            .collect();
+        for m in 1..=subcycle {
+            let fraction = m as f64 / subcycle as f64;
+            let d_sub: Vec<f64> = d_n
+                .iter()
+                .zip(d_candidate)
+                .map(|(old, new)| old + fraction * (new - old))
+                .collect();
+            self.set_geometry(&d_sub, &ddot);
+            futures::executor::block_on(solver.advance(field, self.dt_fluid)).unwrap();
+        }
+    }
+}