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/core/rtx-backend-metal/src/lib.rs

@@ -388,6 +388,23 @@ impl Backend for MetalBackend {
         ops::conv::avg_pool2d(&input, kernel_size, stride, padding, count_include_pad)
     }
 
+    // ==================== Index Operations ====================
+
+    fn index_select<const D: usize>(
+        tensor: Self::TensorPrimitive<D>,
+        indices: &[usize],
+    ) -> Self::TensorPrimitive<D> {
+        ops::index::index_select(&tensor, indices)
+    }
+
+    fn index_add<const D: usize>(
+        tensor: Self::TensorPrimitive<D>,
+        indices: &[usize],
+        num_rows: usize,
+    ) -> Self::TensorPrimitive<D> {
+        ops::index::index_add(&tensor, indices, num_rows)
+    }
+
     // ==================== Device Management ====================
 
     fn device<const D: usize>(tensor: &Self::TensorPrimitive<D>) -> Self::Device {

crates/core/rtx-backend-metal/src/ops/conv.rs

@@ -1,6 +1,10 @@
-//! Convolution and pooling operations using Metal GPU shaders.
+//! Convolution and pooling operations.
 //!
-//! These operations dispatch to Metal compute kernels for GPU-accelerated execution.
+//! `conv2d` dispatches to the Metal `conv2d_f32` kernel when the parameters
+//! fit its signature (symmetric stride/padding, dilation 1, groups 1); all
+//! other cases — and both pooling ops, which have no Metal kernel in
+//! rtx-metal — are computed on host (same fallback pattern as
+//! `reduction::sum_dim`).
 
 use crate::MetalTensorPrimitive;
 use rtx_metal::{MetalBuffer, MetalBufferUsage};
@@ -8,7 +12,7 @@ use rtx_metal::{MetalBuffer, MetalBufferUsage};
 #[cfg(target_os = "macos")]
 use rtx_metal::ops::nn;
 
-/// 2D Convolution using Metal GPU kernel.
+/// 2D Convolution.
 pub fn conv2d(
     input: &MetalTensorPrimitive<4>,
     weight: &MetalTensorPrimitive<4>,
@@ -18,8 +22,8 @@ pub fn conv2d(
     dilation: [usize; 2],
     groups: usize,
 ) -> MetalTensorPrimitive<4> {
-    let [batch, _in_channels, in_h, in_w] = input.shape;
+    let [batch, in_channels, in_h, in_w] = input.shape;
-    let [out_channels, _in_channels_per_group, kernel_h, kernel_w] = weight.shape;
+    let [out_channels, in_channels_per_group, kernel_h, kernel_w] = weight.shape;
 
     let out_h = (in_h + 2 * padding[0] - dilation[0] * (kernel_h - 1) - 1) / stride[0] + 1;
     let out_w = (in_w + 2 * padding[1] - dilation[1] * (kernel_w - 1) - 1) / stride[1] + 1;
@@ -28,40 +32,93 @@ pub fn conv2d(
     let out_numel: usize = out_shape.iter().product();
 
     let device = input.device.metal_device();
-    let mut output = MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
-        .expect("Failed to allocate output buffer for conv2d");
 
+    // The Metal kernel only supports symmetric stride/padding, dilation 1,
+    // groups 1. Dispatch to it when possible.
     #[cfg(target_os = "macos")]
     {
-        nn::conv2d(
+        if stride[0] == stride[1]
-            device,
+            && padding[0] == padding[1]
-            input.data(),
+            && dilation == [1, 1]
-            weight.data(),
+            && groups == 1
-            bias.map(|b| b.data()),
+        {
-            &mut output,
+            let mut output = MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
-            input.shape,
+                .expect("Failed to allocate output buffer for conv2d");
-            weight.shape,
+            nn::conv2d(
-            stride,
+                device,
-            padding,
+                input.data(),
-            dilation,
+                weight.data(),
-            groups,
+                bias.map(|b| b.data()),
-        )
+                &mut output,
-        .expect("Failed to execute Metal conv2d kernel");
+                batch,
+                in_channels,
+                out_channels,
+                in_h,
+                in_w,
+                kernel_h,
+                kernel_w,
+                stride[0],
+                padding[0],
+            )
+            .expect("Failed to execute Metal conv2d kernel");
+            return MetalTensorPrimitive::new(output, out_shape, input.device.clone());
+        }
     }
 
-    #[cfg(not(target_os = "macos"))]
+    // Host fallback: naive direct convolution with full stride / padding /
-    {
+    // dilation / groups support.
-        let _input_data = input.to_vec();
+    let input_data = input.to_vec();
-        let _weight_data = weight.to_vec();
+    let weight_data = weight.to_vec();
-        let _bias_data = bias.map(|b| b.to_vec());
+    let bias_data = bias.map(|b| b.to_vec());
-        let result = vec![0.0f32; out_numel];
+    let mut result = vec![0.0f32; out_numel];
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
+
+    let out_channels_per_group = out_channels / groups;
+
+    for b in 0..batch {
+        for oc in 0..out_channels {
+            let group = oc / out_channels_per_group;
+            for oh in 0..out_h {
+                for ow in 0..out_w {
+                    let mut acc = bias_data.as_ref().map_or(0.0, |bd| bd[oc]);
+                    for ic in 0..in_channels_per_group {
+                        let ic_global = group * in_channels_per_group + ic;
+                        for kh in 0..kernel_h {
+                            for kw in 0..kernel_w {
+                                let ih = (oh * stride[0] + kh * dilation[0]) as isize
+                                    - padding[0] as isize;
+                                let iw = (ow * stride[1] + kw * dilation[1]) as isize
+                                    - padding[1] as isize;
+                                if ih >= 0
+                                    && iw >= 0
+                                    && (ih as usize) < in_h
+                                    && (iw as usize) < in_w
+                                {
+                                    let in_idx = b * in_channels * in_h * in_w
+                                        + ic_global * in_h * in_w
+                                        + ih as usize * in_w
+                                        + iw as usize;
+                                    let w_idx = oc * in_channels_per_group * kernel_h * kernel_w
+                                        + ic * kernel_h * kernel_w
+                                        + kh * kernel_w
+                                        + kw;
+                                    acc += input_data[in_idx] * weight_data[w_idx];
+                                }
+                            }
+                        }
+                    }
+                    let out_idx =
+                        b * out_channels * out_h * out_w + oc * out_h * out_w + oh * out_w + ow;
+                    result[out_idx] = acc;
+                }
+            }
+        }
     }
 
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
     MetalTensorPrimitive::new(output, out_shape, input.device.clone())
 }
 
-/// 2D max pooling using Metal GPU kernel.
+/// 2D max pooling (host fallback — no Metal kernel in rtx-metal).
 pub fn max_pool2d(
     input: &MetalTensorPrimitive<4>,
     kernel_size: [usize; 2],
@@ -77,67 +134,40 @@ pub fn max_pool2d(
     let out_numel: usize = out_shape.iter().product();
 
     let device = input.device.metal_device();
-    let mut output = MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
+    let input_data = input.to_vec();
-        .expect("Failed to allocate output buffer for max_pool2d");
+    let mut result = vec![f32::NEG_INFINITY; out_numel];
 
-    #[cfg(target_os = "macos")]
+    for b in 0..batch {
-    {
+        for c in 0..channels {
-        nn::max_pool2d(
+            for oh in 0..out_h {
-            device,
+                for ow in 0..out_w {
-            input.data(),
+                    let mut max_val = f32::NEG_INFINITY;
-            &mut output,
+                    for kh in 0..kernel_size[0] {
-            input.shape,
+                        for kw in 0..kernel_size[1] {
-            kernel_size,
+                            let ih = (oh * stride[0] + kh) as isize - padding[0] as isize;
-            stride,
+                            let iw = (ow * stride[1] + kw) as isize - padding[1] as isize;
-            padding,
+                            if ih >= 0 && iw >= 0 && (ih as usize) < in_h && (iw as usize) < in_w {
-        )
+                                let in_idx = b * channels * in_h * in_w
-        .expect("Failed to execute Metal max_pool2d kernel");
+                                    + c * in_h * in_w
-    }
+                                    + ih as usize * in_w
-
+                                    + iw as usize;
-    #[cfg(not(target_os = "macos"))]
+                                max_val = max_val.max(input_data[in_idx]);
-    {
-        let input_data = input.to_vec();
-        let mut result = vec![f32::NEG_INFINITY; out_numel];
-
-        for b in 0..batch {
-            for c in 0..channels {
-                for oh in 0..out_h {
-                    for ow in 0..out_w {
-                        let mut max_val = f32::NEG_INFINITY;
-                        for kh in 0..kernel_size[0] {
-                            for kw in 0..kernel_size[1] {
-                                let ih = oh * stride[0] + kh;
-                                let iw = ow * stride[1] + kw;
-                                if ih >= padding[0]
-                                    && iw >= padding[1]
-                                    && ih < in_h + padding[0]
-                                    && iw < in_w + padding[1]
-                                {
-                                    let ih_actual = ih - padding[0];
-                                    let iw_actual = iw - padding[1];
-                                    let in_idx = b * channels * in_h * in_w
-                                        + c * in_h * in_w
-                                        + ih_actual * in_w
-                                        + iw_actual;
-                                    max_val = max_val.max(input_data[in_idx]);
-                                }
                             }
                         }
-                        let out_idx =
-                            b * channels * out_h * out_w + c * out_h * out_w + oh * out_w + ow;
-                        result[out_idx] = max_val;
                     }
+                    let out_idx =
+                        b * channels * out_h * out_w + c * out_h * out_w + oh * out_w + ow;
+                    result[out_idx] = max_val;
                 }
             }
         }
-
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
     }
 
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
     MetalTensorPrimitive::new(output, out_shape, input.device.clone())
 }
 
-/// 2D average pooling using Metal GPU kernel.
+/// 2D average pooling (host fallback — no Metal kernel in rtx-metal).
 pub fn avg_pool2d(
     input: &MetalTensorPrimitive<4>,
     kernel_size: [usize; 2],
@@ -154,67 +184,39 @@ pub fn avg_pool2d(
     let out_numel: usize = out_shape.iter().product();
 
     let device = input.device.metal_device();
-    let mut output = MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
+    let input_data = input.to_vec();
-        .expect("Failed to allocate output buffer for avg_pool2d");
+    let mut result = vec![0.0f32; out_numel];
 
-    #[cfg(target_os = "macos")]
+    for b in 0..batch {
-    {
+        for c in 0..channels {
-        nn::avg_pool2d(
+            for oh in 0..out_h {
-            device,
+                for ow in 0..out_w {
-            input.data(),
+                    let mut sum = 0.0f32;
-            &mut output,
+                    let mut count = 0usize;
-            input.shape,
+                    for kh in 0..kernel_size[0] {
-            kernel_size,
+                        for kw in 0..kernel_size[1] {
-            stride,
+                            let ih = (oh * stride[0] + kh) as isize - padding[0] as isize;
-            padding,
+                            let iw = (ow * stride[1] + kw) as isize - padding[1] as isize;
-            count_include_pad,
+                            if ih >= 0 && iw >= 0 && (ih as usize) < in_h && (iw as usize) < in_w {
-        )
+                                let in_idx = b * channels * in_h * in_w
-        .expect("Failed to execute Metal avg_pool2d kernel");
+                                    + c * in_h * in_w
-    }
+                                    + ih as usize * in_w
-
+                                    + iw as usize;
-    #[cfg(not(target_os = "macos"))]
+                                sum += input_data[in_idx];
-    {
+                                count += 1;
-        let input_data = input.to_vec();
+                            } else if count_include_pad {
-        let mut result = vec![0.0f32; out_numel];
+                                count += 1;
-
-        for b in 0..batch {
-            for c in 0..channels {
-                for oh in 0..out_h {
-                    for ow in 0..out_w {
-                        let mut sum = 0.0f32;
-                        let mut count = 0;
-                        for kh in 0..kernel_size[0] {
-                            for kw in 0..kernel_size[1] {
-                                let ih = oh * stride[0] + kh;
-                                let iw = ow * stride[1] + kw;
-                                if ih >= padding[0]
-                                    && iw >= padding[1]
-                                    && ih < in_h + padding[0]
-                                    && iw < in_w + padding[1]
-                                {
-                                    let ih_actual = ih - padding[0];
-                                    let iw_actual = iw - padding[1];
-                                    let in_idx = b * channels * in_h * in_w
-                                        + c * in_h * in_w
-                                        + ih_actual * in_w
-                                        + iw_actual;
-                                    sum += input_data[in_idx];
-                                    count += 1;
-                                } else if count_include_pad {
-                                    count += 1;
-                                }
                             }
                         }
-                        let out_idx =
-                            b * channels * out_h * out_w + c * out_h * out_w + oh * out_w + ow;
-                        result[out_idx] = if count > 0 { sum / count as f32 } else { 0.0 };
                     }
+                    let out_idx =
+                        b * channels * out_h * out_w + c * out_h * out_w + oh * out_w + ow;
+                    result[out_idx] = if count > 0 { sum / count as f32 } else { 0.0 };
                 }
             }
         }
-
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
     }
 
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
     MetalTensorPrimitive::new(output, out_shape, input.device.clone())
 }

crates/core/rtx-backend-metal/src/ops/index.rs

@@ -0,0 +1,282 @@
+//! Index operations (`index_select` / `index_add`) via sparse one-hot SpMM.
+//!
+//! Gathering rows `out[i] = x[indices[i]]` is exactly `S @ X` where `S` is the
+//! `[E x N]` one-hot selection matrix with `S[i, indices[i]] = 1`. Scatter-add
+//! (`index_add`) is the adjoint, `S^T @ X`, whose CSR form is built directly
+//! with a counting sort (so the GPU kernel — one thread per output element,
+//! looping a row's nonzeros — needs no atomics even with duplicate indices).
+//!
+//! CSR matrices are cached per thread, keyed by the exact index list and
+//! matrix dimensions, so repeated calls with a static graph topology (the GNN
+//! message-passing case) rebuild nothing.
+
+use crate::MetalTensorPrimitive;
+
+#[cfg(target_os = "macos")]
+use rtx_metal::{MetalBuffer, MetalBufferUsage};
+
+#[cfg(target_os = "macos")]
+use rtx_metal::sparse::{spmm_csr, CsrMatrix};
+
+#[cfg(target_os = "macos")]
+use std::cell::RefCell;
+#[cfg(target_os = "macos")]
+use std::collections::HashMap;
+#[cfg(target_os = "macos")]
+use std::rc::Rc;
+
+/// Maximum number of cached CSR matrices per thread before the cache is
+/// cleared. Each entry holds `O(nnz)` GPU memory, so keep this small.
+#[cfg(target_os = "macos")]
+const CSR_CACHE_CAP: usize = 32;
+
+#[cfg(target_os = "macos")]
+#[derive(PartialEq, Eq, Hash, Clone)]
+struct CsrCacheKey {
+    /// False: selection matrix `S` (E x N). True: adjoint `S^T` (num_rows x E).
+    transposed: bool,
+    /// The dense dimension (N for select, num_rows for add).
+    dim: usize,
+    /// The exact index list (collision-proof; hashing is cheap relative to
+    /// the SpMM itself and only rebuilt entries pay the CSR construction).
+    indices: Vec<usize>,
+}
+
+#[cfg(target_os = "macos")]
+thread_local! {
+    static CSR_CACHE: RefCell<HashMap<CsrCacheKey, Rc<CsrMatrix<f32>>>> =
+        RefCell::new(HashMap::new());
+}
+
+/// Get or build the `[E x N]` one-hot selection CSR for `indices`.
+#[cfg(target_os = "macos")]
+fn cached_select_csr(
+    device: &rtx_metal::MetalDevice,
+    indices: &[usize],
+    num_src_rows: usize,
+) -> Option<Rc<CsrMatrix<f32>>> {
+    let key = CsrCacheKey {
+        transposed: false,
+        dim: num_src_rows,
+        indices: indices.to_vec(),
+    };
+    CSR_CACHE.with(|cache| {
+        let mut cache = cache.borrow_mut();
+        if let Some(csr) = cache.get(&key) {
+            return Some(csr.clone());
+        }
+        let e = indices.len();
+        let row_ptr: Vec<i32> = (0..=e as i32).collect();
+        let col_indices: Vec<i32> = indices.iter().map(|&i| i as i32).collect();
+        let values = vec![1.0f32; e];
+        let csr =
+            CsrMatrix::new(device, e, num_src_rows, &row_ptr, &col_indices, &values).ok()?;
+        let csr = Rc::new(csr);
+        if cache.len() >= CSR_CACHE_CAP {
+            cache.clear();
+        }
+        cache.insert(key, csr.clone());
+        Some(csr)
+    })
+}
+
+/// Get or build the `[num_rows x E]` transposed one-hot CSR (counting sort).
+#[cfg(target_os = "macos")]
+fn cached_scatter_csr(
+    device: &rtx_metal::MetalDevice,
+    indices: &[usize],
+    num_rows: usize,
+) -> Option<Rc<CsrMatrix<f32>>> {
+    let key = CsrCacheKey {
+        transposed: true,
+        dim: num_rows,
+        indices: indices.to_vec(),
+    };
+    CSR_CACHE.with(|cache| {
+        let mut cache = cache.borrow_mut();
+        if let Some(csr) = cache.get(&key) {
+            return Some(csr.clone());
+        }
+        let e = indices.len();
+        // Counting sort: row r of S^T holds the input positions i with
+        // indices[i] == r.
+        let mut row_ptr = vec![0i32; num_rows + 1];
+        for &r in indices {
+            row_ptr[r + 1] += 1;
+        }
+        for r in 0..num_rows {
+            row_ptr[r + 1] += row_ptr[r];
+        }
+        let mut next: Vec<i32> = row_ptr[..num_rows].to_vec();
+        let mut col_indices = vec![0i32; e];
+        for (i, &r) in indices.iter().enumerate() {
+            col_indices[next[r] as usize] = i as i32;
+            next[r] += 1;
+        }
+        let values = vec![1.0f32; e];
+        let csr = CsrMatrix::new(device, num_rows, e, &row_ptr, &col_indices, &values).ok()?;
+        let csr = Rc::new(csr);
+        if cache.len() >= CSR_CACHE_CAP {
+            cache.clear();
+        }
+        cache.insert(key, csr.clone());
+        Some(csr)
+    })
+}
+
+/// Host fallback mirroring the `Backend::index_select` default body.
+fn host_index_select<const D: usize>(
+    tensor: &MetalTensorPrimitive<D>,
+    indices: &[usize],
+    row_len: usize,
+    out_shape: [usize; D],
+) -> MetalTensorPrimitive<D> {
+    let src = tensor.to_vec();
+    let mut out = Vec::with_capacity(indices.len() * row_len);
+    for &idx in indices {
+        out.extend_from_slice(&src[idx * row_len..(idx + 1) * row_len]);
+    }
+    super::creation::from_data(&out, out_shape, &tensor.device)
+}
+
+/// Host fallback mirroring the `Backend::index_add` default body.
+fn host_index_add<const D: usize>(
+    tensor: &MetalTensorPrimitive<D>,
+    indices: &[usize],
+    num_rows: usize,
+    row_len: usize,
+    out_shape: [usize; D],
+) -> MetalTensorPrimitive<D> {
+    let src = tensor.to_vec();
+    let mut out = vec![0.0f32; num_rows * row_len];
+    for (i, &idx) in indices.iter().enumerate() {
+        let dst = &mut out[idx * row_len..(idx + 1) * row_len];
+        let row = &src[i * row_len..(i + 1) * row_len];
+        for (d, &s) in dst.iter_mut().zip(row) {
+            *d += s;
+        }
+    }
+    super::creation::from_data(&out, out_shape, &tensor.device)
+}
+
+/// Gather rows along dim 0: `out[i, ..] = tensor[indices[i], ..]`.
+///
+/// GPU path: one-hot CSR `[E x N]` times the dense tensor via `spmm_csr`.
+/// Falls back to the host implementation for degenerate shapes or any
+/// sparse-pipeline failure.
+///
+/// # Panics
+/// Panics if any index is `>= shape[0]`, if `D == 0`, or if the tensor is
+/// not contiguous.
+pub fn index_select<const D: usize>(
+    tensor: &MetalTensorPrimitive<D>,
+    indices: &[usize],
+) -> MetalTensorPrimitive<D> {
+    assert!(D >= 1, "index_select requires at least one dimension");
+    assert!(
+        tensor.is_contiguous(),
+        "index_select: tensor must be contiguous"
+    );
+    let shape = tensor.shape;
+    let num_rows = shape[0];
+    let row_len: usize = shape[1..].iter().product();
+    for &idx in indices {
+        assert!(
+            idx < num_rows,
+            "index_select: index {idx} out of range for {num_rows} rows"
+        );
+    }
+    let mut out_shape = shape;
+    out_shape[0] = indices.len();
+
+    #[cfg(target_os = "macos")]
+    {
+        if row_len > 0
+            && !indices.is_empty()
+            && num_rows <= i32::MAX as usize
+            && indices.len() < i32::MAX as usize
+        {
+            let device = tensor.device.metal_device();
+            if let Some(csr) = cached_select_csr(device, indices, num_rows) {
+                let out_numel = indices.len() * row_len;
+                if let Ok(mut output) =
+                    MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
+                {
+                    if spmm_csr(device, &csr, tensor.data(), &mut output, row_len).is_ok() {
+                        return MetalTensorPrimitive::new(
+                            output,
+                            out_shape,
+                            tensor.device.clone(),
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    host_index_select(tensor, indices, row_len, out_shape)
+}
+
+/// Scatter-add rows along dim 0 into a zero tensor with `num_rows` rows:
+/// `out[indices[i], ..] += tensor[i, ..]`. Adjoint of [`index_select`].
+///
+/// GPU path: transposed one-hot CSR `[num_rows x E]` (built by counting
+/// sort — duplicate indices land in the same CSR row, so the kernel
+/// accumulates them without atomics) times the dense tensor via `spmm_csr`.
+///
+/// # Panics
+/// Panics if `indices.len() != shape[0]`, if any index is `>= num_rows`,
+/// if `D == 0`, or if the tensor is not contiguous.
+pub fn index_add<const D: usize>(
+    tensor: &MetalTensorPrimitive<D>,
+    indices: &[usize],
+    num_rows: usize,
+) -> MetalTensorPrimitive<D> {
+    assert!(D >= 1, "index_add requires at least one dimension");
+    assert!(
+        tensor.is_contiguous(),
+        "index_add: tensor must be contiguous"
+    );
+    let shape = tensor.shape;
+    assert_eq!(
+        indices.len(),
+        shape[0],
+        "index_add: indices.len() must equal the number of input rows"
+    );
+    let row_len: usize = shape[1..].iter().product();
+    for &idx in indices {
+        assert!(
+            idx < num_rows,
+            "index_add: index {idx} out of range for {num_rows} rows"
+        );
+    }
+    let mut out_shape = shape;
+    out_shape[0] = num_rows;
+
+    #[cfg(target_os = "macos")]
+    {
+        if row_len > 0
+            && num_rows > 0
+            && indices.len() < i32::MAX as usize
+            && num_rows <= i32::MAX as usize
+        {
+            let device = tensor.device.metal_device();
+            if let Some(csr) = cached_scatter_csr(device, indices, num_rows) {
+                let out_numel = num_rows * row_len;
+                if let Ok(mut output) =
+                    MetalBuffer::new(device, out_numel, MetalBufferUsage::Shared)
+                {
+                    if spmm_csr(device, &csr, tensor.data(), &mut output, row_len).is_ok() {
+                        return MetalTensorPrimitive::new(
+                            output,
+                            out_shape,
+                            tensor.device.clone(),
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    host_index_add(tensor, indices, num_rows, row_len, out_shape)
+}

crates/core/rtx-backend-metal/src/ops/mod.rs

@@ -10,6 +10,7 @@ pub mod conv;
 pub mod creation;
 pub mod device;
 pub mod gemm;
+pub mod index;
 pub mod normalization;
 pub mod reduction;
 pub mod shape;

crates/core/rtx-backend-metal/src/ops/reduction.rs

@@ -157,35 +157,24 @@ pub fn min<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimi
     MetalTensorPrimitive::new(output, [1], tensor.device.clone())
 }
 
-/// Variance of all elements using Metal GPU kernel.
+/// Variance of all elements (population variance, divisor `n`).
+///
+/// rtx-metal has no dedicated variance kernel; computed on host
+/// (same fallback pattern as `sum_dim`).
 pub fn var<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimitive<1> {
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, 1, MetalBufferUsage::Shared)
+    let n = data.len() as f32;
-        .expect("Failed to allocate output buffer for var");
+    let mean: f32 = data.iter().sum::<f32>() / n;
-
+    let variance: f32 = data
-    #[cfg(target_os = "macos")]
+        .iter()
-    {
+        .map(|&x| {
-        tensor_ops::var(device, tensor.data(), &mut output)
+            let diff = x - mean;
-            .expect("Failed to execute Metal var kernel");
+            diff * diff
-    }
+        })
-
+        .sum::<f32>()
-    #[cfg(not(target_os = "macos"))]
+        / n;
-    {
+    let output = MetalBuffer::from_slice(device, &[variance]).expect("Failed to create buffer");
-        let data = tensor.to_vec();
-        let n = data.len() as f32;
-        let mean: f32 = data.iter().sum::<f32>() / n;
-        let variance: f32 = data
-            .iter()
-            .map(|&x| {
-                let diff = x - mean;
-                diff * diff
-            })
-            .sum::<f32>()
-            / n;
-        output = MetalBuffer::from_slice(device, &[variance]).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, [1], tensor.device.clone())
 }
 

crates/core/rtx-backend-metal/src/ops/shape.rs

@@ -1,4 +1,9 @@
 //! Shape manipulation operations.
+//!
+//! All ops here return **contiguous** tensors. The rest of the backend
+//! (elementwise kernels, MPS matmul, `to_vec`) reads raw buffers and ignores
+//! strides, so a stride-swapped "view" would silently corrupt every
+//! downstream op — `swap_dims` therefore physically permutes the data.
 
 use crate::MetalTensorPrimitive;
 use rtx_metal::MetalBuffer;
@@ -12,6 +17,10 @@ pub fn reshape<const D1: usize, const D2: usize>(
     let old_numel: usize = tensor.shape.iter().product();
     let new_numel: usize = shape.iter().product();
     assert_eq!(old_numel, new_numel, "Total elements must remain the same");
+    assert!(
+        tensor.is_contiguous(),
+        "reshape requires a contiguous tensor (raw-buffer reinterpretation)"
+    );
 
     let strides = MetalTensorPrimitive::<D2>::compute_strides(&shape);
 
@@ -24,7 +33,7 @@ pub fn reshape<const D1: usize, const D2: usize>(
     }
 }
 
-/// Transpose the last two dimensions.
+/// Transpose the last two dimensions (physical, contiguous result).
 pub fn transpose<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimitive<D> {
     if D < 2 {
         return tensor.clone();
@@ -32,7 +41,8 @@ pub fn transpose<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTenso
     swap_dims(tensor, D - 2, D - 1)
 }
 
-/// Swap two dimensions.
+/// Swap two dimensions, physically permuting the data so the result is
+/// contiguous in the new shape.
 pub fn swap_dims<const D: usize>(
     tensor: &MetalTensorPrimitive<D>,
     dim1: usize,
@@ -45,30 +55,34 @@ pub fn swap_dims<const D: usize>(
     }
 
     let mut new_shape = tensor.shape;
-    let mut new_strides = tensor.strides;
     new_shape.swap(dim1, dim2);
-    new_strides.swap(dim1, dim2);
+    let out_strides = MetalTensorPrimitive::<D>::compute_strides(&new_shape);
 
-    let data = tensor.to_vec();
+    let src = tensor.to_vec();
     let numel = tensor.numel();
     let mut result = vec![0.0f32; numel];
 
-    for i in 0..numel {
+    for (out_pos, r) in result.iter_mut().enumerate() {
-        let mut indices = [0usize; D];
+        // Decompose the contiguous output position into a new_shape
-        let mut remaining = i;
+        // multi-index, map it back (swap) to an input multi-index, and read
+        // through the input's own strides + offset.
+        let mut rem = out_pos;
+        let mut in_off = tensor.offset;
         for d in 0..D {
-            indices[d] = remaining / tensor.strides[d];
+            let id = rem / out_strides[d];
-            remaining %= tensor.strides[d];
+            rem %= out_strides[d];
+            // output dim d corresponds to input dim d, with dim1/dim2 swapped;
+            // the index value along src_dim equals the output index along d
+            let src_dim = if d == dim1 {
+                dim2
+            } else if d == dim2 {
+                dim1
+            } else {
+                d
+            };
+            in_off += id * tensor.strides[src_dim];
         }
-
+        *r = src[in_off];
-        indices.swap(dim1, dim2);
-
-        let mut new_idx = 0;
-        for d in 0..D {
-            new_idx += indices[d] * new_strides[d];
-        }
-
-        result[new_idx] = data[i];
     }
 
     let metal_data = MetalBuffer::from_slice(tensor.device.metal_device(), &result)
@@ -77,7 +91,7 @@ pub fn swap_dims<const D: usize>(
     MetalTensorPrimitive {
         data: Arc::new(metal_data),
         shape: new_shape,
-        strides: new_strides,
+        strides: out_strides,
         device: tensor.device.clone(),
         offset: 0,
     }

crates/core/rtx-backend-metal/src/ops/unary.rs

@@ -131,132 +131,63 @@ pub fn abs<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimi
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }
 
-/// Element-wise sine using Metal GPU kernel.
+/// Element-wise sine.
+///
+/// rtx-metal has no dedicated kernel for this op; computed on host
+/// (same fallback pattern as `reduction::sum_dim`).
 pub fn sin<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimitive<D> {
-    let numel = tensor.numel();
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, numel, MetalBufferUsage::Shared)
+    let result: Vec<f32> = data.iter().map(|x| x.sin()).collect();
-        .expect("Failed to allocate output buffer for sin");
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-
-    #[cfg(target_os = "macos")]
-    {
-        tensor_ops::sin(device, tensor.data(), &mut output)
-            .expect("Failed to execute Metal sin kernel");
-    }
-
-    #[cfg(not(target_os = "macos"))]
-    {
-        let data = tensor.to_vec();
-        let result: Vec<f32> = data.iter().map(|x| x.sin()).collect();
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }
 
-/// Element-wise cosine using Metal GPU kernel.
+/// Element-wise cosine.
+///
+/// rtx-metal has no dedicated kernel for this op; computed on host
+/// (same fallback pattern as `reduction::sum_dim`).
 pub fn cos<const D: usize>(tensor: &MetalTensorPrimitive<D>) -> MetalTensorPrimitive<D> {
-    let numel = tensor.numel();
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, numel, MetalBufferUsage::Shared)
+    let result: Vec<f32> = data.iter().map(|x| x.cos()).collect();
-        .expect("Failed to allocate output buffer for cos");
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-
-    #[cfg(target_os = "macos")]
-    {
-        tensor_ops::cos(device, tensor.data(), &mut output)
-            .expect("Failed to execute Metal cos kernel");
-    }
-
-    #[cfg(not(target_os = "macos"))]
-    {
-        let data = tensor.to_vec();
-        let result: Vec<f32> = data.iter().map(|x| x.cos()).collect();
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }
 
-/// Element-wise power using Metal GPU kernel.
+/// Element-wise power.
+///
+/// rtx-metal has no dedicated kernel for this op; computed on host
+/// (same fallback pattern as `reduction::sum_dim`).
 pub fn pow<const D: usize>(tensor: &MetalTensorPrimitive<D>, exp: f32) -> MetalTensorPrimitive<D> {
-    let numel = tensor.numel();
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, numel, MetalBufferUsage::Shared)
+    let result: Vec<f32> = data.iter().map(|x| x.powf(exp)).collect();
-        .expect("Failed to allocate output buffer for pow");
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-
-    #[cfg(target_os = "macos")]
-    {
-        tensor_ops::pow(device, tensor.data(), &mut output, exp)
-            .expect("Failed to execute Metal pow kernel");
-    }
-
-    #[cfg(not(target_os = "macos"))]
-    {
-        let data = tensor.to_vec();
-        let result: Vec<f32> = data.iter().map(|x| x.powf(exp)).collect();
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }
 
-/// Clamp tensor values to a range using Metal GPU kernel.
+/// Clamp tensor values to a range.
-pub fn clamp<const D: usize>(
+///
-    tensor: &MetalTensorPrimitive<D>,
+/// rtx-metal has no dedicated kernel for this op; computed on host
-    min: f32,
+/// (same fallback pattern as `reduction::sum_dim`).
-    max: f32,
+pub fn clamp<const D: usize>(tensor: &MetalTensorPrimitive<D>, min: f32,
-) -> MetalTensorPrimitive<D> {
+    max: f32) -> MetalTensorPrimitive<D> {
-    let numel = tensor.numel();
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, numel, MetalBufferUsage::Shared)
+    let result: Vec<f32> = data.iter().map(|x| x.clamp(min, max)).collect();
-        .expect("Failed to allocate output buffer for clamp");
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-
-    #[cfg(target_os = "macos")]
-    {
-        tensor_ops::clamp(device, tensor.data(), &mut output, min, max)
-            .expect("Failed to execute Metal clamp kernel");
-    }
-
-    #[cfg(not(target_os = "macos"))]
-    {
-        let data = tensor.to_vec();
-        let result: Vec<f32> = data.iter().map(|x| x.clamp(min, max)).collect();
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }
 
-/// Greater than scalar comparison using Metal GPU kernel.
+/// Greater-than-scalar comparison (1.0 / 0.0 mask).
-pub fn gt_scalar<const D: usize>(
+///
-    tensor: &MetalTensorPrimitive<D>,
+/// rtx-metal has no dedicated kernel for this op; computed on host
-    value: f32,
+/// (same fallback pattern as `reduction::sum_dim`).
-) -> MetalTensorPrimitive<D> {
+pub fn gt_scalar<const D: usize>(tensor: &MetalTensorPrimitive<D>, value: f32) -> MetalTensorPrimitive<D> {
-    let numel = tensor.numel();
     let device = tensor.device.metal_device();
-
+    let data = tensor.to_vec();
-    let mut output = MetalBuffer::new(device, numel, MetalBufferUsage::Shared)
+    let result: Vec<f32> = data.iter().map(|&x| if x > value { 1.0 } else { 0.0 }).collect();
-        .expect("Failed to allocate output buffer for gt_scalar");
+    let output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-
-    #[cfg(target_os = "macos")]
-    {
-        tensor_ops::gt_scalar(device, tensor.data(), &mut output, value)
-            .expect("Failed to execute Metal gt_scalar kernel");
-    }
-
-    #[cfg(not(target_os = "macos"))]
-    {
-        let data = tensor.to_vec();
-        let result: Vec<f32> = data
-            .iter()
-            .map(|&x| if x > value { 1.0 } else { 0.0 })
-            .collect();
-        output = MetalBuffer::from_slice(device, &result).expect("Failed to create buffer");
-    }
-
     MetalTensorPrimitive::new(output, tensor.shape, tensor.device.clone())
 }

crates/core/rtx-backend-metal/tests/backend_parity_tests.rs

@@ -5,7 +5,7 @@
 
 use rtx_backend_metal::MetalDeviceWrapper;
 use rtx_backend_metal::ops::{
-    activation, basic, creation, device as dev_ops, gemm, reduction, unary,
+    activation, basic, creation, device as dev_ops, gemm, reduction, shape, unary,
 };
 
 /// Tolerance for floating-point comparisons.
@@ -608,3 +608,69 @@ fn test_parity_matmul_larger() {
         }
     }
 }
+
+#[test]
+fn test_parity_transpose_2d() {
+    let device = MetalDeviceWrapper::new().expect("Metal device");
+    let data: Vec<f32> = (0..12).map(|x| x as f32).collect(); // [3, 4]
+    let t = creation::from_data(&data, [3, 4], &device);
+    let tt = shape::transpose(&t);
+    // physical transpose: result must be contiguous row-major [4, 3]
+    assert!(tt.is_contiguous(), "transpose must return a contiguous tensor");
+    let got = dev_ops::copy_to_host(&tt);
+    let mut want = vec![0.0f32; 12];
+    for i in 0..3 {
+        for j in 0..4 {
+            want[j * 3 + i] = data[i * 4 + j];
+        }
+    }
+    assert_eq!(got, want);
+    // and transposing back must round-trip
+    let back = dev_ops::copy_to_host(&shape::transpose(&tt));
+    assert_eq!(back, data);
+}
+
+#[test]
+fn test_parity_swap_dims_3d() {
+    let device = MetalDeviceWrapper::new().expect("Metal device");
+    let data: Vec<f32> = (0..24).map(|x| x as f32).collect(); // [2, 3, 4]
+    let t = creation::from_data(&data, [2, 3, 4], &device);
+    let s = shape::swap_dims(&t, 0, 2); // -> [4, 3, 2]
+    assert!(s.is_contiguous());
+    let got = dev_ops::copy_to_host(&s);
+    let mut want = vec![0.0f32; 24];
+    for a in 0..2 {
+        for b in 0..3 {
+            for c in 0..4 {
+                want[c * 6 + b * 2 + a] = data[a * 12 + b * 4 + c];
+            }
+        }
+    }
+    assert_eq!(got, want);
+}
+
+#[test]
+fn test_parity_matmul_after_transpose() {
+    // the exact pattern the autograd matmul backward uses:
+    // grad_a = grad_c @ b^T — a transposed operand must feed MPS correctly.
+    let device = MetalDeviceWrapper::new().expect("Metal device");
+    let a: Vec<f32> = (0..6).map(|x| x as f32 + 1.0).collect(); // [2, 3]
+    let b: Vec<f32> = (0..12).map(|x| (x as f32) * 0.5 - 2.0).collect(); // [4, 3] -> b^T [3, 4]
+    let ta = creation::from_data(&a, [2, 3], &device);
+    let tb = creation::from_data(&b, [4, 3], &device);
+    let c = gemm::matmul(&ta, &shape::transpose(&tb));
+    let got = dev_ops::copy_to_host(&c);
+    let mut want = vec![0.0f32; 8];
+    for i in 0..2 {
+        for j in 0..4 {
+            let mut s = 0.0f32;
+            for k in 0..3 {
+                s += a[i * 3 + k] * b[j * 3 + k];
+            }
+            want[i * 4 + j] = s;
+        }
+    }
+    for (g, w) in got.iter().zip(want.iter()) {
+        assert!((g - w).abs() < 1e-4, "matmul-after-transpose mismatch: {got:?} vs {want:?}");
+    }
+}

crates/core/rtx-backend-metal/tests/index_ops_parity_tests.rs

@@ -0,0 +1,234 @@
+//! Parity tests for Metal `index_select` / `index_add` (sparse one-hot SpMM)
+//! against the CPU reference semantics of the `Backend` trait defaults.
+
+use rtx_backend::Backend;
+use rtx_backend_metal::{MetalBackend, MetalDeviceWrapper};
+
+fn device() -> MetalDeviceWrapper {
+    MetalDeviceWrapper::new().expect("Metal device")
+}
+
+/// CPU reference for index_select (mirrors the Backend trait default).
+fn cpu_index_select(src: &[f32], row_len: usize, indices: &[usize]) -> Vec<f32> {
+    let mut out = Vec::with_capacity(indices.len() * row_len);
+    for &idx in indices {
+        out.extend_from_slice(&src[idx * row_len..(idx + 1) * row_len]);
+    }
+    out
+}
+
+/// CPU reference for index_add (mirrors the Backend trait default).
+fn cpu_index_add(src: &[f32], row_len: usize, indices: &[usize], num_rows: usize) -> Vec<f32> {
+    let mut out = vec![0.0f32; num_rows * row_len];
+    for (i, &idx) in indices.iter().enumerate() {
+        for d in 0..row_len {
+            out[idx * row_len + d] += src[i * row_len + d];
+        }
+    }
+    out
+}
+
+fn assert_close(got: &[f32], want: &[f32]) {
+    assert_eq!(got.len(), want.len(), "length mismatch");
+    for (i, (g, w)) in got.iter().zip(want.iter()).enumerate() {
+        assert!(
+            (g - w).abs() <= 1e-5 + 1e-4 * w.abs(),
+            "mismatch at {i}: got {g}, want {w}"
+        );
+    }
+}
+
+#[test]
+fn test_index_select_d2_basic() {
+    let dev = device();
+    let data: Vec<f32> = (0..20).map(|x| x as f32).collect(); // [5, 4]
+    let t = MetalBackend::from_data(&data, [5, 4], &dev);
+    let indices = [4usize, 0, 2];
+    let out = MetalBackend::index_select(t, &indices);
+    assert_eq!(MetalBackend::shape(&out), [3, 4]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_select(&data, 4, &indices),
+    );
+}
+
+#[test]
+fn test_index_select_duplicate_indices() {
+    let dev = device();
+    let data: Vec<f32> = (0..12).map(|x| x as f32 * 0.5).collect(); // [4, 3]
+    let t = MetalBackend::from_data(&data, [4, 3], &dev);
+    let indices = [1usize, 1, 3, 1, 0, 3];
+    let out = MetalBackend::index_select(t, &indices);
+    assert_eq!(MetalBackend::shape(&out), [6, 3]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_select(&data, 3, &indices),
+    );
+}
+
+#[test]
+fn test_index_select_d1() {
+    let dev = device();
+    let data = [10.0f32, 11.0, 12.0, 13.0, 14.0];
+    let t = MetalBackend::from_data(&data, [5], &dev);
+    let indices = [2usize, 2, 4, 0];
+    let out = MetalBackend::index_select(t, &indices);
+    assert_eq!(MetalBackend::shape(&out), [4]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_select(&data, 1, &indices),
+    );
+}
+
+#[test]
+fn test_index_select_d3() {
+    let dev = device();
+    let data: Vec<f32> = (0..24).map(|x| x as f32).collect(); // [4, 3, 2]
+    let t = MetalBackend::from_data(&data, [4, 3, 2], &dev);
+    let indices = [3usize, 1, 1, 0, 2];
+    let out = MetalBackend::index_select(t, &indices);
+    assert_eq!(MetalBackend::shape(&out), [5, 3, 2]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_select(&data, 6, &indices),
+    );
+}
+
+#[test]
+fn test_index_add_d2_basic() {
+    let dev = device();
+    let data: Vec<f32> = (0..12).map(|x| x as f32).collect(); // [3, 4]
+    let t = MetalBackend::from_data(&data, [3, 4], &dev);
+    let indices = [1usize, 3, 0];
+    let out = MetalBackend::index_add(t, &indices, 5);
+    assert_eq!(MetalBackend::shape(&out), [5, 4]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_add(&data, 4, &indices, 5),
+    );
+}
+
+#[test]
+fn test_index_add_duplicate_indices_accumulate() {
+    let dev = device();
+    let data: Vec<f32> = (0..15).map(|x| x as f32 * 0.25).collect(); // [5, 3]
+    let t = MetalBackend::from_data(&data, [5, 3], &dev);
+    // Rows 0, 2, 4 of input all land on output row 1; rows 1, 3 on row 3.
+    let indices = [1usize, 3, 1, 3, 1];
+    let out = MetalBackend::index_add(t, &indices, 4);
+    assert_eq!(MetalBackend::shape(&out), [4, 3]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_add(&data, 3, &indices, 4),
+    );
+}
+
+#[test]
+fn test_index_add_empty_rows_stay_zero() {
+    let dev = device();
+    let data = [1.0f32, 2.0, 3.0, 4.0]; // [2, 2]
+    let t = MetalBackend::from_data(&data, [2, 2], &dev);
+    // Output rows 0, 2, 4, 5 are never referenced -> must be exactly zero.
+    let indices = [3usize, 1];
+    let out = MetalBackend::index_add(t, &indices, 6);
+    let got = MetalBackend::to_data(&out);
+    assert_close(&got, &cpu_index_add(&data, 2, &indices, 6));
+    for &r in &[0usize, 2, 4, 5] {
+        assert_eq!(got[r * 2], 0.0, "row {r} not zero");
+        assert_eq!(got[r * 2 + 1], 0.0, "row {r} not zero");
+    }
+}
+
+#[test]
+fn test_index_add_d1() {
+    let dev = device();
+    let data = [5.0f32, 7.0, 11.0];
+    let t = MetalBackend::from_data(&data, [3], &dev);
+    let indices = [2usize, 0, 2];
+    let out = MetalBackend::index_add(t, &indices, 4);
+    assert_eq!(MetalBackend::shape(&out), [4]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_add(&data, 1, &indices, 4),
+    );
+}
+
+#[test]
+fn test_index_add_d3() {
+    let dev = device();
+    let data: Vec<f32> = (0..24).map(|x| (x as f32).sin()).collect(); // [4, 2, 3]
+    let t = MetalBackend::from_data(&data, [4, 2, 3], &dev);
+    let indices = [0usize, 2, 0, 1];
+    let out = MetalBackend::index_add(t, &indices, 3);
+    assert_eq!(MetalBackend::shape(&out), [3, 2, 3]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_add(&data, 6, &indices, 3),
+    );
+}
+
+#[test]
+fn test_index_select_large() {
+    let dev = device();
+    let n = 5000usize;
+    let f = 64usize;
+    let data: Vec<f32> = (0..n * f).map(|x| ((x * 2654435761) % 1000) as f32 * 0.001).collect();
+    let t = MetalBackend::from_data(&data, [n, f], &dev);
+    // Pseudo-random gather with repeats, GNN-edge style.
+    let indices: Vec<usize> = (0..3 * n).map(|i| (i * 40503) % n).collect();
+    let out = MetalBackend::index_select(t, &indices);
+    assert_eq!(MetalBackend::shape(&out), [3 * n, f]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_select(&data, f, &indices),
+    );
+}
+
+#[test]
+fn test_index_add_large() {
+    let dev = device();
+    let e = 15000usize;
+    let n = 5000usize;
+    let f = 64usize;
+    let data: Vec<f32> = (0..e * f).map(|x| ((x * 2246822519) % 1000) as f32 * 0.001 - 0.5).collect();
+    let t = MetalBackend::from_data(&data, [e, f], &dev);
+    let indices: Vec<usize> = (0..e).map(|i| (i * 40503) % n).collect();
+    let out = MetalBackend::index_add(t, &indices, n);
+    assert_eq!(MetalBackend::shape(&out), [n, f]);
+    assert_close(
+        &MetalBackend::to_data(&out),
+        &cpu_index_add(&data, f, &indices, n),
+    );
+}
+
+#[test]
+fn test_index_ops_repeated_calls_hit_cache() {
+    let dev = device();
+    let data: Vec<f32> = (0..40).map(|x| x as f32).collect(); // [10, 4]
+    let indices = [9usize, 0, 5, 5, 3];
+    let want = cpu_index_select(&data, 4, &indices);
+    // Same topology, different tensor values, several times over — exercises
+    // the per-thread CSR cache path.
+    for round in 0..5 {
+        let scaled: Vec<f32> = data.iter().map(|x| x * (round as f32 + 1.0)).collect();
+        let t = MetalBackend::from_data(&scaled, [10, 4], &dev);
+        let out = MetalBackend::index_select(t, &indices);
+        let want_scaled: Vec<f32> = want.iter().map(|x| x * (round as f32 + 1.0)).collect();
+        assert_close(&MetalBackend::to_data(&out), &want_scaled);
+    }
+}
+
+#[test]
+fn test_index_select_roundtrip_adjoint() {
+    // index_add(index_select(x, idx), idx, n) counts each row once per
+    // appearance in idx — sanity-check the pair against CPU reference.
+    let dev = device();
+    let data: Vec<f32> = (0..18).map(|x| x as f32 * 0.1).collect(); // [6, 3]
+    let t = MetalBackend::from_data(&data, [6, 3], &dev);
+    let indices = [5usize, 5, 1, 0];
+    let gathered = MetalBackend::index_select(t, &indices);
+    let scattered = MetalBackend::index_add(gathered, &indices, 6);
+    let want_gather = cpu_index_select(&data, 3, &indices);
+    let want = cpu_index_add(&want_gather, 3, &indices, 6);
+    assert_close(&MetalBackend::to_data(&scattered), &want);
+}

crates/core/rtx-backend/tests/backend_parity_tests.rs

@@ -88,6 +88,7 @@ fn test_cpu_baseline_operations() {
 // CUDA Backend Parity Tests
 // ==============================================================================
 
+#[cfg(all(target_os = "linux", target_arch = "x86_64"))]
 mod cuda_parity {
     use super::*;
     use rtx_backend_cpu::{CpuBackend, CpuDevice};
@@ -365,7 +366,8 @@ fn test_backend_feature_detection() {
     // CPU is always available
     println!("  CPU Backend: Available");
 
-    // Check CUDA
+    // Check CUDA (dev-dep only present on x86_64 Linux)
+    #[cfg(all(target_os = "linux", target_arch = "x86_64"))]
     {
         use rtx_backend_cuda::CudaDevice;
         match CudaDevice::new(0) {

crates/specialized/rtx-fsi/src/coupling.rs

@@ -34,6 +34,45 @@
 //! convergence is immediate — which is both why it works and how the
 //! implementation can be tested sharply rather than by "it got there
 //! eventually".
+//!
+//! # IQN-ILS — the vector quasi-Newton coupler
+//!
+//! Aitken's factor is a **scalar**: one relaxation for every interface
+//! degree of freedom. When the coupled map's gain differs across
+//! interface modes — a flag whose tip and root see different added mass,
+//! or a map contaminated by uncorrelated sampling noise — no single
+//! scalar fits, and Aitken grinds or stalls.
+//!
+//! Interface Quasi-Newton with Inverse Least-Squares (IQN-ILS, Degroote,
+//! Bathe & Vierendeels 2009) builds a low-rank secant model of the
+//! interface Jacobian from the residual history instead. Writing one
+//! pass as `x_tilde = pass(x)` with residual `r = x_tilde - x`, each
+//! iteration contributes a column pair `(delta r, delta x_tilde)`; the
+//! update solves the least-squares problem
+//!
+//! ```text
+//! alpha = argmin || V alpha + r ||,     x_next = x_tilde + W alpha
+//! ```
+//!
+//! which is exactly the Newton step for `pass(x) - x = 0` in the
+//! subspace the history spans (verified in the tests against linear maps
+//! with anisotropic gains, where scalar Aitken cannot be exact). Two
+//! properties matter for a noisy embedded-boundary interface:
+//!
+//! - the least-squares projection **filters components of the residual
+//!   that no history column explains** — uncorrelated per-pass sampling
+//!   noise does not steer the update the way it steers a scalar factor;
+//! - columns from the **previous few time steps** can be reused
+//!   ([`IqnIls::with_reuse`]), so even a step that converges in one or
+//!   two passes benefits from a full secant model — the regime a tightly
+//!   time-coupled march actually runs in.
+//!
+//! Near-dependent columns are dropped by a modified Gram–Schmidt filter
+//! whose threshold is **relative to each column's own norm** — an
+//! absolute threshold here would be the same latent scale bug that has
+//! now struck this workspace four times.
+
+use std::collections::VecDeque;
 
 use crate::error::FsiError;
 
@@ -210,6 +249,283 @@ impl Subiterated {
     }
 }
 
+/// How near-dependent a secant column may be to the span of the columns
+/// already accepted before it is dropped, **relative to its own norm**.
+const COLUMN_FILTER: f64 = 1e-8;
+
+/// One secant sample: `(delta residual, delta pass-output)` between two
+/// successive iterations.
+type SecantColumn = (Vec<f64>, Vec<f64>);
+
+/// Interface quasi-Newton driver with inverse least-squares (IQN-ILS).
+///
+/// Keep one instance alive across a time march: with
+/// [`Self::with_reuse`] the secant columns of the last few steps carry
+/// over, and the first pass of a new step already runs against a full
+/// Jacobian model. See the module docs for the method and the tests for
+/// its sharp properties (exact on linear maps, anisotropic gains, scale
+/// invariance, noise stalling at the noise scale instead of diverging).
+#[derive(Debug, Clone, PartialEq)]
+pub struct IqnIls {
+    max_iterations: usize,
+    tolerance: f64,
+    initial_relaxation: f64,
+    steps_retained: usize,
+    /// Newest step first; within a step, newest column first.
+    history: VecDeque<Vec<SecantColumn>>,
+}
+
+impl IqnIls {
+    /// A driver with the given per-step iteration budget and interface
+    /// tolerance. Defaults: first-iteration relaxation 0.5, secant reuse
+    /// over the 2 previous steps.
+    ///
+    /// # Errors
+    /// [`FsiError::InvalidParameter`] for a zero iteration budget or a
+    /// non-positive tolerance.
+    pub fn new(max_iterations: usize, tolerance: f64) -> Result<Self, FsiError> {
+        if max_iterations == 0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "max_iterations",
+                value: 0.0,
+            });
+        }
+        if !tolerance.is_finite() || tolerance <= 0.0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "tolerance",
+                value: tolerance,
+            });
+        }
+        Ok(Self {
+            max_iterations,
+            tolerance,
+            initial_relaxation: 0.5,
+            steps_retained: 2,
+            history: VecDeque::new(),
+        })
+    }
+
+    /// Retain the secant columns of the last `steps` time steps across
+    /// [`Self::solve`] calls (0 = within-step only).
+    #[must_use]
+    pub fn with_reuse(mut self, steps: usize) -> Self {
+        self.steps_retained = steps;
+        self.history.truncate(steps);
+        self
+    }
+
+    /// The relaxation applied when no secant information exists yet
+    /// (the very first pass of the very first step).
+    ///
+    /// # Errors
+    /// [`FsiError::InvalidParameter`] for a non-finite or non-positive
+    /// factor.
+    pub fn with_initial_relaxation(mut self, factor: f64) -> Result<Self, FsiError> {
+        if !factor.is_finite() || factor <= 0.0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "initial relaxation",
+                value: factor,
+            });
+        }
+        self.initial_relaxation = factor;
+        Ok(self)
+    }
+
+    /// Update the interface tolerance for the next [`Self::solve`] call
+    /// (a marching coupler re-budgets per step: the tolerance is
+    /// max(noise floor, a fraction of the step's own increment)).
+    ///
+    /// # Errors
+    /// [`FsiError::InvalidParameter`] for a non-positive tolerance.
+    pub fn set_tolerance(&mut self, tolerance: f64) -> Result<(), FsiError> {
+        if !tolerance.is_finite() || tolerance <= 0.0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "tolerance",
+                value: tolerance,
+            });
+        }
+        self.tolerance = tolerance;
+        Ok(())
+    }
+
+    /// Drive `pass` to an interface fixed point, as
+    /// [`Subiterated::solve`] does, reusing secant history across calls.
+    ///
+    /// # Errors
+    /// The same contract as [`Subiterated::solve`]; on
+    /// [`FsiError::CouplingDiverged`] and
+    /// [`FsiError::CouplingNotConverged`] this step's secant columns are
+    /// still retained — they are genuine samples of the map either way.
+    pub fn solve<F>(&mut self, initial: &[f64], pass: F) -> Result<Converged, FsiError>
+    where
+        F: Fn(&[f64]) -> Vec<f64>,
+    {
+        if initial.is_empty() {
+            return Err(FsiError::EmptyInterface { side: "interface" });
+        }
+
+        let mut x = initial.to_vec();
+        let mut previous: Option<(Vec<f64>, Vec<f64>)> = None; // (r, x_tilde)
+        let mut step_columns: Vec<SecantColumn> = Vec::new();
+        let mut first_norm = None;
+        let mut last_norm = f64::NAN;
+
+        for iteration in 1..=self.max_iterations {
+            let x_tilde = pass(&x);
+            if x_tilde.len() != x.len() {
+                return Err(FsiError::CountMismatch {
+                    field: "coupling pass",
+                    got: x_tilde.len(),
+                    expected: x.len(),
+                });
+            }
+            if let Some(index) = x_tilde.iter().position(|v| !v.is_finite()) {
+                return Err(FsiError::NonFinite {
+                    field: "coupling pass",
+                    index,
+                });
+            }
+
+            let r: Vec<f64> = x_tilde.iter().zip(&x).map(|(new, old)| new - old).collect();
+            let norm = norm_of(&r);
+            let first = *first_norm.get_or_insert(norm);
+            last_norm = norm;
+
+            if let Some((r_prev, xt_prev)) = &previous {
+                step_columns.insert(
+                    0,
+                    (
+                        r.iter().zip(r_prev).map(|(a, b)| a - b).collect(),
+                        x_tilde.iter().zip(xt_prev).map(|(a, b)| a - b).collect(),
+                    ),
+                );
+            }
+
+            if norm <= self.tolerance {
+                self.commit(step_columns);
+                return Ok(Converged {
+                    state: x,
+                    residual: norm,
+                    iterations: iteration,
+                });
+            }
+            if norm > first * DIVERGENCE_FACTOR && first > 0.0 {
+                self.commit(step_columns);
+                return Err(FsiError::CouplingDiverged {
+                    iterations: iteration,
+                    residual: norm,
+                });
+            }
+
+            let columns: Vec<&SecantColumn> = step_columns
+                .iter()
+                .chain(self.history.iter().flatten())
+                .collect();
+            match least_squares_update(&columns, &r) {
+                Some(delta) => {
+                    x = x_tilde.iter().zip(&delta).map(|(a, b)| a + b).collect();
+                }
+                None => {
+                    // No usable secant information yet: one relaxed
+                    // fixed-point step to generate it.
+                    for (value, residual) in x.iter_mut().zip(&r) {
+                        *value += self.initial_relaxation * residual;
+                    }
+                }
+            }
+            previous = Some((r, x_tilde));
+        }
+
+        self.commit(step_columns);
+        Err(FsiError::CouplingNotConverged {
+            iterations: self.max_iterations,
+            residual: last_norm,
+            tolerance: self.tolerance,
+        })
+    }
+
+    /// Retire this step's secant columns into the cross-step history.
+    fn commit(&mut self, step_columns: Vec<SecantColumn>) {
+        if self.steps_retained == 0 || step_columns.is_empty() {
+            return;
+        }
+        self.history.push_front(step_columns);
+        self.history.truncate(self.steps_retained);
+    }
+}
+
+/// The IQN-ILS update `W alpha` with `alpha = argmin || V alpha + r ||`,
+/// via modified Gram–Schmidt with dropping of near-dependent columns
+/// (threshold relative to each column's own norm). `None` when no column
+/// survives — the caller falls back to a relaxed fixed-point step.
+fn least_squares_update(columns: &[&SecantColumn], r: &[f64]) -> Option<Vec<f64>> {
+    if columns.is_empty() {
+        return None;
+    }
+    let n = r.len();
+    // Accepted orthonormal basis q_i, the R entries of each accepted
+    // column, and the index of the original column it came from.
+    let mut basis: Vec<Vec<f64>> = Vec::new();
+    let mut upper: Vec<Vec<f64>> = Vec::new(); // per accepted column: R entries over basis
+    let mut accepted: Vec<usize> = Vec::new();
+    for (index, (v, _)) in columns.iter().enumerate() {
+        debug_assert_eq!(v.len(), n);
+        let original_norm = norm_of(v);
+        if original_norm == 0.0 {
+            continue;
+        }
+        let mut q = v.clone();
+        let mut coefficients = Vec::with_capacity(basis.len());
+        for b in &basis {
+            let dot: f64 = b.iter().zip(&q).map(|(a, c)| a * c).sum();
+            for (qi, bi) in q.iter_mut().zip(b) {
+                *qi -= dot * bi;
+            }
+            coefficients.push(dot);
+        }
+        let remaining = norm_of(&q);
+        if remaining <= COLUMN_FILTER * original_norm {
+            continue;
+        }
+        for value in &mut q {
+            *value /= remaining;
+        }
+        coefficients.push(remaining);
+        basis.push(q);
+        upper.push(coefficients);
+        accepted.push(index);
+        if basis.len() == n {
+            break; // the span is full
+        }
+    }
+    if accepted.is_empty() {
+        return None;
+    }
+    // alpha solves R alpha = Q^T (-r), by back substitution: `upper[j]`
+    // holds column j's entries over basis rows 0..=j.
+    let m = accepted.len();
+    let rhs: Vec<f64> = basis
+        .iter()
+        .map(|q| -q.iter().zip(r).map(|(a, b)| a * b).sum::<f64>())
+        .collect();
+    let mut alpha = vec![0.0; m];
+    for j in (0..m).rev() {
+        let mut sum = rhs[j];
+        for k in j + 1..m {
+            sum -= upper[k][j] * alpha[k];
+        }
+        alpha[j] = sum / upper[j][j];
+    }
+    // W alpha over the accepted columns.
+    let mut delta = vec![0.0; n];
+    for (a, &index) in alpha.iter().zip(&accepted) {
+        for (d, w) in delta.iter_mut().zip(&columns[index].1) {
+            *d += a * w;
+        }
+    }
+    Some(delta)
+}
+
 /// Aitken delta-squared relaxation factor from successive residuals.
 ///
 /// Falls back to the previous factor when the residual barely moved
@@ -423,4 +739,196 @@ mod tests {
             Err(FsiError::NonFinite { .. })
         ));
     }
+
+    // ---- IQN-ILS ----
+
+    /// A linear coupled map with a different gain per interface mode —
+    /// the situation a scalar relaxation factor cannot be exact for.
+    fn anisotropic(gains: &'static [f64]) -> impl Fn(&[f64]) -> Vec<f64> {
+        move |state: &[f64]| state.iter().zip(gains).map(|(x, g)| -g * x + 1.0).collect()
+    }
+
+    #[test]
+    fn iqn_recovers_the_heavy_added_mass_case() {
+        let mut scheme = IqnIls::new(200, 1e-10).expect("valid");
+        let converged = scheme.solve(&[1.0], added_mass(2.5)).expect("converges");
+        assert!(converged.residual < 1e-10);
+    }
+
+    #[test]
+    fn iqn_is_exact_on_a_linear_map_within_dimension_plus_two() {
+        // The least-squares secant model spans the full Jacobian after
+        // `dim` independent columns, so a linear map must converge in at
+        // most dim + 2 passes. Asserting the count pins that the update
+        // is the real IQN-ILS step, not a relaxation that happens to
+        // converge.
+        let gains: &[f64] = &[2.5, -0.8, 3.0, 0.3];
+        let mut scheme = IqnIls::new(200, 1e-12).expect("valid");
+        let converged = scheme
+            .solve(&[1.0, 1.0, 1.0, 1.0], anisotropic(gains))
+            .expect("converges");
+        assert!(
+            converged.iterations <= 6,
+            "expected <= dim + 2 = 6 iterations, took {}",
+            converged.iterations
+        );
+    }
+
+    #[test]
+    fn iqn_beats_aitken_on_anisotropic_gains() {
+        // Mixed attracting/repelling modes: no scalar factor fits both,
+        // so Aitken must grind where the vector secant is exact. This is
+        // the property that makes IQN the standard strong coupler.
+        let gains: &[f64] = &[2.2, -0.9, 1.4, 0.1, 2.9];
+        let initial = [1.0, -1.0, 2.0, 0.5, -0.3];
+        let mut aitken = Subiterated::aitken(500, 1e-10).expect("valid");
+        let mut iqn = IqnIls::new(500, 1e-10).expect("valid");
+        let slow = aitken
+            .solve(&initial, anisotropic(gains))
+            .expect("converges");
+        let fast = iqn.solve(&initial, anisotropic(gains)).expect("converges");
+        assert!(
+            fast.iterations < slow.iterations,
+            "iqn {} vs aitken {}",
+            fast.iterations,
+            slow.iterations
+        );
+    }
+
+    #[test]
+    fn iqn_is_scale_invariant() {
+        // The column filter must be relative to each column's own norm —
+        // the absolute-epsilon species has struck this workspace four
+        // times, once in this very module.
+        let mut scheme = IqnIls::new(200, 1e-20).expect("valid");
+        let converged = scheme
+            .solve(&[1e-9], added_mass(2.5))
+            .expect("IQN must converge regardless of residual scale");
+        assert!(
+            converged.iterations <= 4,
+            "expected the same near-immediate convergence as at scale 1, \
+             took {}",
+            converged.iterations
+        );
+    }
+
+    #[test]
+    fn iqn_reuses_secant_history_across_steps() {
+        // A marching coupler solves the same (linearised) interface
+        // problem step after step. With reuse the second step starts
+        // with a full Jacobian model and must converge in fewer passes
+        // than the first; without reuse it must not.
+        let gains: &[f64] = &[2.5, -0.8, 3.0];
+        let initial = [1.0, 1.0, 1.0];
+
+        let mut with_reuse = IqnIls::new(200, 1e-10).expect("valid").with_reuse(2);
+        let first = with_reuse
+            .solve(&initial, anisotropic(gains))
+            .expect("converges");
+        let second = with_reuse
+            .solve(&initial, anisotropic(gains))
+            .expect("converges");
+        assert!(
+            second.iterations < first.iterations,
+            "reuse should shorten the next step: {} then {}",
+            first.iterations,
+            second.iterations
+        );
+
+        let mut without = IqnIls::new(200, 1e-10).expect("valid").with_reuse(0);
+        let cold_first = without
+            .solve(&initial, anisotropic(gains))
+            .expect("converges");
+        let cold_second = without
+            .solve(&initial, anisotropic(gains))
+            .expect("converges");
+        assert_eq!(
+            cold_first.iterations, cold_second.iterations,
+            "without reuse each step must start cold"
+        );
+    }
+
+    #[test]
+    fn iqn_stalls_at_the_noise_scale_instead_of_diverging() {
+        // The embedded-boundary reality: the pass carries a deterministic
+        // but effectively uncorrelated noise component (mask flips) of a
+        // fixed scale. The coupler must converge to a tolerance ABOVE the
+        // noise scale, and must report (not blow through) one below it.
+        // The noise must vary per PASS, not per state: a continuous
+        // function of the state alone has a genuine fixed point and IQN
+        // legitimately converges onto it to machine precision (the first
+        // draft of this test learned that the hard way). A call counter
+        // models the real thing — successive samplings of the same
+        // geometry never repay the same load once the mask has moved.
+        let noise_scale = 1e-6;
+        let calls = std::cell::Cell::new(0u64);
+        let noisy = |state: &[f64]| -> Vec<f64> {
+            calls.set(calls.get() + 1);
+            let noise = (calls.get() as f64 * 2.399_963).sin() * noise_scale;
+            state.iter().map(|x| -2.5 * x + 1.0 + noise).collect()
+        };
+        let mut above = IqnIls::new(50, 20.0 * noise_scale).expect("valid");
+        let converged = above.solve(&[1.0], noisy).expect("converges above noise");
+        assert!(converged.residual <= 20.0 * noise_scale);
+
+        let mut below = IqnIls::new(50, 1e-12).expect("valid");
+        match below.solve(&[1.0], noisy) {
+            Err(
+                FsiError::CouplingNotConverged { residual, .. }
+                | FsiError::CouplingDiverged { residual, .. },
+            ) => {
+                assert!(
+                    residual < 100.0 * noise_scale,
+                    "stall residual {residual} should sit at the noise scale"
+                );
+            }
+            Ok(converged) => panic!(
+                "cannot genuinely converge below the noise floor \
+                 (residual {})",
+                converged.residual
+            ),
+            Err(other) => panic!("unexpected error species: {other:?}"),
+        }
+    }
+
+    #[test]
+    fn iqn_refusals_match_the_subiterated_contract() {
+        assert!(IqnIls::new(0, 1e-8).is_err());
+        assert!(IqnIls::new(10, 0.0).is_err());
+        assert!(IqnIls::new(10, -1e-8).is_err());
+        assert!(
+            IqnIls::new(10, 1e-8)
+                .expect("valid")
+                .with_initial_relaxation(0.0)
+                .is_err()
+        );
+
+        let mut scheme = IqnIls::new(10, 1e-8).expect("valid");
+        assert!(scheme.solve(&[], added_mass(0.5)).is_err());
+        assert!(matches!(
+            scheme.solve(&[1.0, 2.0], |_: &[f64]| vec![0.0]),
+            Err(FsiError::CountMismatch { .. })
+        ));
+        assert!(matches!(
+            scheme.solve(&[1.0], |_: &[f64]| vec![f64::NAN]),
+            Err(FsiError::NonFinite { .. })
+        ));
+        assert!(scheme.set_tolerance(-1.0).is_err());
+        assert!(scheme.set_tolerance(1e-6).is_ok());
+    }
+
+    #[test]
+    fn iqn_budget_exhaustion_is_reported_not_hidden() {
+        // A pass that ignores its input never generates a secant column
+        // pointing at the fixed point of anything; the budget must be
+        // reported honestly.
+        let mut scheme = IqnIls::new(3, 1e-12)
+            .expect("valid")
+            .with_initial_relaxation(1e-6)
+            .expect("valid");
+        assert!(matches!(
+            scheme.solve(&[1.0], added_mass(0.999)),
+            Err(FsiError::CouplingNotConverged { .. } | FsiError::CouplingDiverged { .. })
+        ));
+    }
 }

crates/specialized/rtx-fsi/src/lib.rs

@@ -38,8 +38,10 @@
 
 mod coupling;
 mod error;
+mod smoothing;
 mod transfer;
 
-pub use coupling::{Converged, Subiterated};
+pub use coupling::{Converged, IqnIls, Subiterated};
 pub use error::FsiError;
+pub use smoothing::smooth_tractions;
 pub use transfer::{FluidFace, WettedSurface};

crates/specialized/rtx-fsi/src/smoothing.rs

@@ -0,0 +1,309 @@
+//! Smoothing sampled tractions along the wetted surface.
+//!
+//! # Why the load needs smoothing at all
+//!
+//! An embedded-boundary fluid samples tractions by reconstructing the
+//! near-wall field from whichever cells are currently fluid. When the
+//! interface moves — even by a vanishing amount — a cell can flip between
+//! fluid and solid and every sample whose stencil contains it jumps by a
+//! finite amount. Mapped through the structure's per-step compliance this
+//! becomes the **interface noise floor**: the smallest displacement
+//! tolerance a partitioned coupling can actually reach (measured on
+//! Turek–Hron FSI2: ~1.3e-4 m per pass at full inflow, and it rides with
+//! the loads). Tighter time coupling than the floor allows is blocked,
+//! because wall-velocity noise is displacement-tolerance / dt.
+//!
+//! The flips are a cell-resolution artefact: the discretisation cannot
+//! represent traction variation below the cell scale in the first place,
+//! so averaging the sampled tractions over a stencil about that scale
+//! removes noise the samples were never entitled to carry.
+//!
+//! # The kernel, and why every factor is continuous
+//!
+//! Each smoothed traction is a normalised weighted average over the
+//! samples within `radius` of it along the surface:
+//!
+//! ```text
+//! t'_i = sum_j k(|s_i - s_j|) c_ij A_j t_j  /  sum_j k(|s_i - s_j|) c_ij A_j
+//! ```
+//!
+//! - `k` is a triangular kernel in **arclength** `s` (cumulative centroid
+//!   distance): samples separated by a gap — e.g. the part of a flag
+//!   buried in its mounting cylinder — sit far apart in arclength and
+//!   never mix.
+//! - `c_ij = max(0, n_i . n_j)^2` keeps averaging from mixing tractions
+//!   across corners: `sigma . n` on the two sides of a corner are loads
+//!   in different directions, and averaging the vectors would manufacture
+//!   a spurious tangential load. The factor is **smooth** in the normals,
+//!   deliberately: a hard angular cutoff would make the smoothed load a
+//!   discontinuous function of the interface geometry, and a coupling
+//!   subiteration bounces on exactly such discontinuities (the
+//!   clamp-don't-drop finding from the spike guard).
+//! - `A_j` weights by face area, so the average is the area-consistent
+//!   one and a constant traction field is reproduced exactly.
+//!
+//! On a uniformly sampled straight stretch the kernel matrix is
+//! symmetric with unit column sums, so the total force over the interior
+//! is conserved exactly; end effects and corners redistribute load only
+//! within a kernel radius, at the scale the sampling could not resolve
+//! anyway.
+
+use nalgebra::Vector3;
+
+use crate::error::FsiError;
+use crate::transfer::FluidFace;
+
+/// Smooth sampled tractions with a triangular moving average of
+/// half-width `radius` in surface arclength, weighted by face area and by
+/// normal similarity (see the module docs for the kernel and its
+/// continuity rationale).
+///
+/// `faces` must be ordered along the surface — arclength is accumulated
+/// from consecutive centroid distances. A `radius` of zero returns the
+/// tractions unchanged.
+///
+/// # Errors
+/// - [`FsiError::CountMismatch`] if `tractions` and `faces` differ in
+///   length.
+/// - [`FsiError::InvalidParameter`] for a negative or non-finite radius.
+/// - [`FsiError::NonFinite`] for a non-finite traction sample.
+pub fn smooth_tractions(
+    faces: &[FluidFace],
+    tractions: &[Vector3<f64>],
+    radius: f64,
+) -> Result<Vec<Vector3<f64>>, FsiError> {
+    if tractions.len() != faces.len() {
+        return Err(FsiError::CountMismatch {
+            field: "tractions",
+            got: tractions.len(),
+            expected: faces.len(),
+        });
+    }
+    if !radius.is_finite() || radius < 0.0 {
+        return Err(FsiError::InvalidParameter {
+            parameter: "smoothing radius",
+            value: radius,
+        });
+    }
+    if let Some(index) = tractions
+        .iter()
+        .position(|t| !t.iter().all(|v| v.is_finite()))
+    {
+        return Err(FsiError::NonFinite {
+            field: "tractions",
+            index,
+        });
+    }
+    if radius == 0.0 || faces.is_empty() {
+        return Ok(tractions.to_vec());
+    }
+
+    // Cumulative arclength along the ordered samples.
+    let mut s = Vec::with_capacity(faces.len());
+    let mut acc = 0.0;
+    s.push(0.0);
+    for pair in faces.windows(2) {
+        acc += (pair[1].centroid - pair[0].centroid).norm();
+        s.push(acc);
+    }
+
+    let mut smoothed = Vec::with_capacity(faces.len());
+    for i in 0..faces.len() {
+        // The window is a contiguous index range because arclength is
+        // monotone in the ordering.
+        let lo = (0..i)
+            .rev()
+            .take_while(|&j| s[i] - s[j] < radius)
+            .last()
+            .unwrap_or(i);
+        let hi = (i + 1..faces.len())
+            .take_while(|&j| s[j] - s[i] < radius)
+            .last()
+            .unwrap_or(i);
+        let mut sum = Vector3::zeros();
+        let mut weight_sum = 0.0;
+        for j in lo..=hi {
+            let kernel = 1.0 - (s[i] - s[j]).abs() / radius;
+            let alignment = faces[i].normal.dot(&faces[j].normal).max(0.0).powi(2);
+            let w = kernel * alignment * faces[j].area;
+            sum += w * tractions[j];
+            weight_sum += w;
+        }
+        // The self term always contributes (kernel 1, alignment 1), so
+        // the denominator cannot vanish for a face with positive area.
+        smoothed.push(sum / weight_sum);
+    }
+    Ok(smoothed)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A straight horizontal stretch of `n` uniformly spaced samples with
+    /// upward normals — the interior of a wall, as the sampler sees it.
+    fn straight_faces(n: usize, spacing: f64) -> Vec<FluidFace> {
+        (0..n)
+            .map(|i| FluidFace {
+                centroid: Vector3::new(i as f64 * spacing, 0.0, 0.0),
+                normal: Vector3::new(0.0, 1.0, 0.0),
+                area: spacing,
+            })
+            .collect()
+    }
+
+    #[test]
+    fn a_constant_field_is_reproduced_exactly() {
+        // The weights are normalised, so any constant must pass through
+        // untouched — including at the ends, where the window truncates.
+        let faces = straight_faces(20, 0.1);
+        let tractions = vec![Vector3::new(3.0, -2.0, 0.0); 20];
+        let smoothed = smooth_tractions(&faces, &tractions, 0.25).unwrap();
+        for t in &smoothed {
+            assert!((t - Vector3::new(3.0, -2.0, 0.0)).norm() < 1e-14);
+        }
+    }
+
+    #[test]
+    fn zero_radius_is_the_identity() {
+        let faces = straight_faces(5, 0.1);
+        let tractions: Vec<_> = (0..5)
+            .map(|i| Vector3::new(i as f64, -(i as f64), 0.0))
+            .collect();
+        let smoothed = smooth_tractions(&faces, &tractions, 0.0).unwrap();
+        assert_eq!(smoothed, tractions);
+    }
+
+    #[test]
+    fn a_single_sample_spike_is_reduced_and_its_force_conserved() {
+        // The mask-flip signature: one sample jumps by a finite amount.
+        // Smoothing must spread it (peak reduced) without losing the
+        // impulse (interior column sums are one on a uniform stretch).
+        let n = 21;
+        let faces = straight_faces(n, 0.1);
+        let mut tractions = vec![Vector3::zeros(); n];
+        tractions[10] = Vector3::new(0.0, 5.0, 0.0);
+        let smoothed = smooth_tractions(&faces, &tractions, 0.25).unwrap();
+
+        let peak = smoothed.iter().map(|t| t.norm()).fold(0.0, f64::max);
+        assert!(
+            peak < 0.6 * 5.0,
+            "spike should spread over the window, peak still {peak}"
+        );
+        let total_before: Vector3<f64> =
+            faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
+        let total_after: Vector3<f64> = faces.iter().zip(&smoothed).map(|(f, t)| t * f.area).sum();
+        assert!(
+            (total_after - total_before).norm() < 1e-12 * total_before.norm(),
+            "interior spike force changed: {} vs {}",
+            total_after.y,
+            total_before.y
+        );
+    }
+
+    #[test]
+    fn tractions_do_not_bleed_across_a_right_angle_corner() {
+        // Two perpendicular stretches meeting at a corner (a flag tip).
+        // sigma.n on the two sides are loads in different directions;
+        // max(0, n_i.n_j)^2 = 0 across the corner, so each side smooths
+        // only among its own.
+        let spacing = 0.1;
+        let mut faces = Vec::new();
+        for i in 0..5 {
+            faces.push(FluidFace {
+                centroid: Vector3::new(i as f64 * spacing, 0.0, 0.0),
+                normal: Vector3::new(0.0, -1.0, 0.0),
+                area: spacing,
+            });
+        }
+        for j in 0..5 {
+            faces.push(FluidFace {
+                centroid: Vector3::new(4.0 * spacing, (j + 1) as f64 * spacing, 0.0),
+                normal: Vector3::new(1.0, 0.0, 0.0),
+                area: spacing,
+            });
+        }
+        let mut tractions = vec![Vector3::new(0.0, -1.0, 0.0); 5];
+        tractions.extend(vec![Vector3::new(2.0, 0.0, 0.0); 5]);
+        let smoothed = smooth_tractions(&faces, &tractions, 0.35).unwrap();
+        for (k, t) in smoothed.iter().enumerate() {
+            if k < 5 {
+                assert!(
+                    (t - Vector3::new(0.0, -1.0, 0.0)).norm() < 1e-14,
+                    "bottom sample {k} contaminated across the corner: {t:?}"
+                );
+            } else {
+                assert!(
+                    (t - Vector3::new(2.0, 0.0, 0.0)).norm() < 1e-14,
+                    "side sample {k} contaminated across the corner: {t:?}"
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn samples_across_an_arclength_gap_do_not_mix() {
+        // A buried stretch (samples skipped) leaves consecutive kept
+        // samples far apart in arclength; the kernel must not reach over.
+        let spacing = 0.1;
+        let mut faces = straight_faces(3, spacing);
+        for i in 0..3 {
+            faces.push(FluidFace {
+                centroid: Vector3::new(10.0 + i as f64 * spacing, 0.0, 0.0),
+                normal: Vector3::new(0.0, 1.0, 0.0),
+                area: spacing,
+            });
+        }
+        let mut tractions = vec![Vector3::new(0.0, 1.0, 0.0); 3];
+        tractions.extend(vec![Vector3::new(0.0, -1.0, 0.0); 3]);
+        let smoothed = smooth_tractions(&faces, &tractions, 0.25).unwrap();
+        for (k, t) in smoothed.iter().enumerate() {
+            let expected = if k < 3 { 1.0 } else { -1.0 };
+            assert!(
+                (t.y - expected).abs() < 1e-14,
+                "sample {k} mixed across the gap: {t:?}"
+            );
+        }
+    }
+
+    #[test]
+    fn mismatched_lengths_are_refused() {
+        let faces = straight_faces(4, 0.1);
+        let tractions = vec![Vector3::zeros(); 3];
+        assert!(matches!(
+            smooth_tractions(&faces, &tractions, 0.1),
+            Err(FsiError::CountMismatch { .. })
+        ));
+    }
+
+    #[test]
+    fn an_invalid_radius_is_refused() {
+        let faces = straight_faces(4, 0.1);
+        let tractions = vec![Vector3::zeros(); 4];
+        assert!(matches!(
+            smooth_tractions(&faces, &tractions, -0.1),
+            Err(FsiError::InvalidParameter { .. })
+        ));
+        assert!(matches!(
+            smooth_tractions(&faces, &tractions, f64::NAN),
+            Err(FsiError::InvalidParameter { .. })
+        ));
+    }
+
+    #[test]
+    fn a_non_finite_traction_is_refused() {
+        let faces = straight_faces(4, 0.1);
+        let mut tractions = vec![Vector3::zeros(); 4];
+        tractions[2].x = f64::NAN;
+        assert!(matches!(
+            smooth_tractions(&faces, &tractions, 0.1),
+            Err(FsiError::NonFinite { .. })
+        ));
+    }
+
+    #[test]
+    fn an_empty_surface_smooths_to_empty() {
+        let smoothed = smooth_tractions(&[], &[], 0.1).unwrap();
+        assert!(smoothed.is_empty());
+    }
+}

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();
+        }
+    }
+}

crates/specialized/rtx-fsi/tests/fsi2_interface_noise.rs

@@ -0,0 +1,416 @@
+//! The FSI2 interface noise floor, measured directly — the lever the
+//! benchmark's resonant cycle waits on.
+//!
+//! # What is being measured
+//!
+//! The partitioned coupling iterates one pass: subcycled fluid march on a
+//! candidate interface geometry, traction sampling on that geometry, one
+//! structure step under the sampled load. The tenth-session study found
+//! that this map is not continuous at small scales: a vanishing interface
+//! change flips embedded-mask cells, the load jumps by a finite amount,
+//! and the jump maps through Newmark's `beta dt_c^2 / m` into an
+//! end-of-step displacement jump. That jump size is the **interface noise
+//! floor**: no coupling tolerance below it is reachable, and the
+//! step-to-step scatter of the accepted interface at the tolerance feeds
+//! the fluid wall-velocity noise of tolerance / dt_c (solver_status.md
+//! §"C2 — FSI2").
+//!
+//! This test measures the floor as the modulus of continuity of the real
+//! pass map, on the real machinery (`fsi2_harness`), at full inflow: a
+//! fine sweep of bending amplitudes (every pass from the same snapshot),
+//! whose maximum successive difference is the mask-flip jump — plus a
+//! decades ladder around zero for the small-scale behaviour.
+//!
+//! # What the 2026-08-21 measurements found (ny = 62, t = 4 s)
+//!
+//! 1. **The floor is NOT in the traction sampling.** Surface smoothing of
+//!    the sampled tractions ([`rtx_fsi::smooth_tractions`]) at radii of
+//!    1–3 cells leaves the flip-scan floor unchanged to 0.2% (3.05e-5 at
+//!    every radius): the flip's load jump is **coherent through the
+//!    fluid field itself** — the mask rebuild shifts the pressure
+//!    solution around the flipped cell and every nearby sample moves
+//!    together — and a surface moving-average preserves exactly such
+//!    coherent shifts. Smoothing therefore stays available but OFF by
+//!    default; reaching for it against this floor is a measured dead
+//!    end. The committed run keeps one smoothed scan alive so this
+//!    attribution stays loud.
+//! 2. **The flip scan does NOT transfer across subcycles** — a
+//!    dt_c^2-scaling hypothesis for the floor (load jump through
+//!    Newmark's `beta dt_c^2 / m`) was tried against a subcycle-2 scan
+//!    and REFUTED in that operationalization: max jump 5.4e-5 vs 3.0e-5
+//!    at subcycle 8, median 15x LARGER. The scan's successive passes
+//!    differ by a fixed geometry increment, so their wall-velocity
+//!    difference is increment / dt_c, and the smooth velocity-response
+//!    trend inflates as dt_c shrinks until it swamps the flip signal.
+//!    A cross-subcycle floor claim needs the OPERATIONAL measurement
+//!    instead:
+//! 3. **The stall measurement** — the real release-step map iterated at
+//!    a deep tolerance with the per-pass residuals traced — is the
+//!    operational floor, and it came out FAR below every march
+//!    tolerance: s8 aitken 3.4e-9 / iqn 1.6e-9 in a 12-pass budget, s2
+//!    aitken 6.5e-10 / iqn 6.3e-10 CONVERGED below 1e-9 in 5–6 passes.
+//!    The flip jumps are events at specific geometries, not a floor
+//!    under every step; a typical step's map is locally smooth. The
+//!    tenth session's subcycle-2 blowup (tolerance held at 2e-4 while
+//!    dt_c shrank — wall-velocity noise tol / dt_c) was a tolerance
+//!    mis-budgeting, not an impassable floor: tighter coupling is open
+//!    at a tolerance the map demonstrably supports (~1e-5 leaves three
+//!    decades of margin). Occasional flip-straddling steps still stall
+//!    at the jump scale — the march's stall-accept handles those.
+//!
+//! Environment knobs: `RTX_NOISE_NY` (default 62), `RTX_NOISE_T` (rigid
+//! march horizon, default 4 s), `RTX_NOISE_SUBCYCLE` (baseline subcycle,
+//! default 8), `RTX_NOISE_SCAN` (flip-scan resolution, default 40
+//! passes), `RTX_NOISE_EPS` (flip-scan amplitude, default 5e-4 m),
+//! `RTX_NOISE_FULL` (nonzero: sweep radii 0–3 h and subcycles {8, 4, 2,
+//! 1} instead of the committed set).
+
+mod fsi2_harness;
+
+use fsi2_harness::{FLAG_X0, FLAG_X1, Fsi2Harness, clamp_left, env_or};
+use rtx_cfd::solvers::incompressible::{EmbeddedPisoSolver, FlowField};
+use rtx_fea::analysis::{AnalysisConfig, ConvergenceCriteria, NonlinearDynamicAnalysis};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::MaterialId;
+
+fn norm(v: &[f64]) -> f64 {
+    v.iter().map(|x| x * x).sum::<f64>().sqrt()
+}
+
+fn sub(a: &[f64], b: &[f64]) -> Vec<f64> {
+    a.iter().zip(b).map(|(x, y)| x - y).collect()
+}
+
+#[test]
+#[allow(clippy::too_many_lines)]
+fn fsi2_interface_noise_floor() {
+    let ny = env_or("RTX_NOISE_NY", 62.0) as usize;
+    let t_probe = env_or("RTX_NOISE_T", 4.0);
+    let base_subcycle = env_or("RTX_NOISE_SUBCYCLE", 8.0) as usize;
+    let scan_passes = env_or("RTX_NOISE_SCAN", 40.0) as usize;
+    let eps_max = env_or("RTX_NOISE_EPS", 5e-4);
+    let full_sweep = env_or("RTX_NOISE_FULL", 0.0) != 0.0;
+    let flag_nx = 35;
+
+    let (mut harness, mut solver, mut field) = Fsi2Harness::build(ny, flag_nx, 0.0);
+    let dt_fluid = harness.dt_fluid;
+
+    // Rigid march to operating loads (the floor rides with the loads —
+    // measuring at startup would understate it by orders of magnitude).
+    let start = std::time::Instant::now();
+    let rigid_steps = (t_probe / dt_fluid).round() as usize;
+    for _ in 0..rigid_steps {
+        futures::executor::block_on(solver.advance(&mut field, dt_fluid)).unwrap();
+    }
+    println!(
+        "  rigid march: {rigid_steps} steps to t = {t_probe:.1} s in {:.0} s wall",
+        start.elapsed().as_secs_f64()
+    );
+
+    let zero_d = vec![0.0; 2 * harness.interface.wetted.len()];
+    let fluid_saved = solver.snapshot();
+    let field_saved = field.clone();
+
+    // A smooth cantilever-bending perturbation pattern, unit tip
+    // amplitude: p_y = ((x - x0)/(x1 - x0))^2, p_x = 0 — the shape a
+    // subiteration increment actually has.
+    let pattern: Vec<f64> = harness
+        .interface
+        .reference
+        .iter()
+        .flat_map(|&(x, _)| {
+            let s = (x - FLAG_X0) / (FLAG_X1 - FLAG_X0);
+            [0.0, s * s]
+        })
+        .collect();
+
+    // The flag stepper is rebuilt per subcycle: the coupled dt (and with
+    // it Newmark's beta dt^2 / m response to a load jump) is exactly
+    // what the scaling measurement varies.
+    let flag_mesh = harness.mesh.clone();
+    let make_analysis = move |subcycle: usize| {
+        let mut db = MaterialDatabase::new();
+        db.add_material(
+            MaterialId(0),
+            LinearElastic::new(fsi2_harness::E_S, fsi2_harness::NU_S)
+                .with_density(fsi2_harness::RHO_S),
+            None,
+        );
+        NonlinearDynamicAnalysis::new(
+            flag_mesh.clone(),
+            db,
+            clamp_left(&flag_mesh),
+            dt_fluid * subcycle as f64,
+            1,
+            AnalysisConfig::default(),
+        )
+        .with_total_lagrangian()
+        .with_convergence_criteria(ConvergenceCriteria {
+            max_iterations: 60,
+            ..ConvergenceCriteria::default()
+        })
+    };
+
+    // One config's floor: pass-map continuity at a given subcycle and
+    // smoothing radius, all passes from the same saved fluid state.
+    let measure = |harness: &mut Fsi2Harness,
+                   solver: &mut EmbeddedPisoSolver,
+                   subcycle: usize,
+                   radius_in_h: f64,
+                   ladder: bool|
+     -> (f64, f64) {
+        harness.smooth_radius = radius_in_h * harness.h;
+        // Reset the shared body geometry: a previous measurement's last
+        // candidate must not leak into this one's initial load sampling
+        // (it did, before this line — a 4.5e-5 phantom first residual).
+        harness.set_geometry(&zero_d, &zero_d);
+        let analysis = make_analysis(subcycle);
+        let mut flag = analysis.stepper().unwrap();
+        let wetted_dofs: Vec<[usize; 2]> = harness
+            .interface
+            .wetted
+            .iter()
+            .map(|&id| {
+                let dofs = flag.node_dofs(id);
+                [dofs[0], dofs[1]]
+            })
+            .collect();
+        let (nodal0, _, _) = harness.sample_load(solver, &field_saved, &zero_d);
+        flag.set_nodal_forces(&nodal0);
+        let flag_state = flag.rest_state().unwrap();
+
+        let mut pass = |d_candidate: &[f64]| -> Vec<f64> {
+            solver.restore(&fluid_saved);
+            let mut trial_field: FlowField = field_saved.clone();
+            harness.advance_subcycled(solver, &mut trial_field, &zero_d, d_candidate, subcycle);
+            let (nodal, _, _) = harness.sample_load(solver, &trial_field, d_candidate);
+            flag.set_nodal_forces(&nodal);
+            let (candidate_state, _) = flag.step(&flag_state).unwrap();
+            let mut d = vec![0.0; 2 * wetted_dofs.len()];
+            for (k, dofs) in wetted_dofs.iter().enumerate() {
+                d[2 * k] = candidate_state.displacement[dofs[0]];
+                d[2 * k + 1] = candidate_state.displacement[dofs[1]];
+            }
+            d
+        };
+
+        let base = pass(&zero_d);
+        if ladder {
+            print!("  subcycle {subcycle} smooth {radius_in_h:.1}h  ladder:");
+            for exp in [-7.0f64, -6.0, -5.0, -4.0] {
+                let eps = 10.0f64.powf(exp);
+                let d: Vec<f64> = pattern.iter().map(|p| eps * p).collect();
+                let response = norm(&sub(&pass(&d), &base));
+                print!("  1e{exp:.0} -> {response:.2e}");
+            }
+            println!();
+        }
+
+        let mut previous = base;
+        let mut max_jump = 0.0f64;
+        let mut jumps = Vec::with_capacity(scan_passes);
+        for k in 1..=scan_passes {
+            let eps = eps_max * k as f64 / scan_passes as f64;
+            let d: Vec<f64> = pattern.iter().map(|p| eps * p).collect();
+            let current = pass(&d);
+            let jump = norm(&sub(&current, &previous));
+            jumps.push(jump);
+            max_jump = max_jump.max(jump);
+            previous = current;
+        }
+        jumps.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        let median_jump = jumps[jumps.len() / 2];
+        println!(
+            "  subcycle {subcycle} smooth {radius_in_h:.1}h  flip scan ({scan_passes} passes to \
+             eps = {eps_max:.1e}): max successive jump {max_jump:.3e}, median {median_jump:.3e}"
+        );
+        (max_jump, median_jump)
+    };
+
+    if full_sweep {
+        for &s in &[base_subcycle, 4, 2, 1] {
+            for r in [0.0, 1.0, 2.0, 3.0] {
+                measure(&mut harness, &mut solver, s, r, r == 0.0);
+            }
+        }
+        println!(
+            "  full sweep total {:.0} s wall",
+            start.elapsed().as_secs_f64()
+        );
+        return;
+    }
+
+    // The stall measurement: iterate the REAL release-step map at an
+    // unreachable tolerance and trace every pass's residual. Where the
+    // subiteration stalls is the operational noise floor — the number a
+    // march's coupling tolerance must sit above — measured per subcycle
+    // and per coupler on the same fluid state.
+    let stall = |harness: &mut Fsi2Harness,
+                 solver: &mut EmbeddedPisoSolver,
+                 subcycle: usize,
+                 coupler: &str|
+     -> (f64, f64, usize) {
+        harness.smooth_radius = 0.0;
+        harness.set_geometry(&zero_d, &zero_d);
+        let analysis = make_analysis(subcycle);
+        let flag = std::cell::RefCell::new(analysis.stepper().unwrap());
+        let wetted_dofs: Vec<[usize; 2]> = harness
+            .interface
+            .wetted
+            .iter()
+            .map(|&id| {
+                let dofs = flag.borrow().node_dofs(id);
+                [dofs[0], dofs[1]]
+            })
+            .collect();
+        let extract = |state: &rtx_fea::analysis::DynamicState| -> Vec<f64> {
+            let mut d = vec![0.0; 2 * wetted_dofs.len()];
+            for (k, dofs) in wetted_dofs.iter().enumerate() {
+                d[2 * k] = state.displacement[dofs[0]];
+                d[2 * k + 1] = state.displacement[dofs[1]];
+            }
+            d
+        };
+        let (nodal0, _, _) = harness.sample_load(solver, &field_saved, &zero_d);
+        flag.borrow_mut().set_nodal_forces(&nodal0);
+        let flag_state = flag.borrow_mut().rest_state().unwrap();
+        // The predictor the march uses: the structure alone under the
+        // committed load.
+        let (predicted, _) = flag.borrow_mut().step(&flag_state).unwrap();
+        let d_pred = extract(&predicted);
+
+        let solver = std::cell::RefCell::new(solver);
+        let trace = std::cell::RefCell::new(Vec::<f64>::new());
+        let harness_ref = &*harness;
+        let pass = |d_candidate: &[f64]| -> Vec<f64> {
+            let mut solver_ref = solver.borrow_mut();
+            solver_ref.restore(&fluid_saved);
+            let mut trial_field: FlowField = field_saved.clone();
+            harness_ref.advance_subcycled(
+                &mut solver_ref,
+                &mut trial_field,
+                &zero_d,
+                d_candidate,
+                subcycle,
+            );
+            let (nodal, _, _) = harness_ref.sample_load(&solver_ref, &trial_field, d_candidate);
+            let mut flag_ref = flag.borrow_mut();
+            flag_ref.set_nodal_forces(&nodal);
+            let (candidate_state, _) = flag_ref.step(&flag_state).unwrap();
+            let d_new = extract(&candidate_state);
+            trace.borrow_mut().push(norm(&sub(&d_new, d_candidate)));
+            d_new
+        };
+
+        // A tolerance at the bottom of what the map could conceivably
+        // support: the point is the trace, not the verdict. (It turned
+        // out to be REACHABLE at subcycle 2 — that reachability is the
+        // finding pinned below.)
+        let budget = 12;
+        let deep = 1e-9;
+        let outcome = if coupler == "iqn" {
+            rtx_fsi::IqnIls::new(budget, deep)
+                .unwrap()
+                .solve(&d_pred, pass)
+        } else {
+            rtx_fsi::Subiterated::aitken(budget, deep)
+                .unwrap()
+                .solve(&d_pred, pass)
+        };
+        drop(outcome); // converged or budget-exhausted — the trace has the data
+        let trace = trace.into_inner();
+        let min = trace.iter().copied().fold(f64::MAX, f64::min);
+        let last = *trace.last().unwrap();
+        println!(
+            "  stall subcycle {subcycle} {coupler}: {} passes, residual first {:.3e} \
+             min {min:.3e} last {last:.3e}",
+            trace.len(),
+            trace.first().unwrap()
+        );
+        (min, last, trace.len())
+    };
+
+    // The committed set: the baseline floor, the smoothed floor (the
+    // attribution guard), the cross-subcycle scan (recorded, unpinned —
+    // see the module docs for why it does not transfer), and the stall
+    // levels per subcycle and coupler.
+    let (floor_base, _) = measure(&mut harness, &mut solver, base_subcycle, 0.0, true);
+    let (floor_smoothed, _) = measure(&mut harness, &mut solver, base_subcycle, 2.0, false);
+    let (floor_tight_scan, _) = measure(&mut harness, &mut solver, base_subcycle / 4, 0.0, true);
+    let stall_8_aitken = stall(&mut harness, &mut solver, base_subcycle, "aitken");
+    let stall_8_iqn = stall(&mut harness, &mut solver, base_subcycle, "iqn");
+    let stall_2_aitken = stall(&mut harness, &mut solver, base_subcycle / 4, "aitken");
+    let stall_2_iqn = stall(&mut harness, &mut solver, base_subcycle / 4, "iqn");
+    println!(
+        "  floors: scan base {floor_base:.3e}, smoothed(2h) {floor_smoothed:.3e}, \
+         subcycle/4 scan {floor_tight_scan:.3e} (trend-contaminated, unpinned); \
+         stalls (min): s8 aitken {:.3e} / iqn {:.3e}, s2 aitken {:.3e} / iqn {:.3e}; \
+         total {:.0} s wall",
+        stall_8_aitken.0,
+        stall_8_iqn.0,
+        stall_2_aitken.0,
+        stall_2_iqn.0,
+        start.elapsed().as_secs_f64()
+    );
+
+    for value in [
+        floor_base,
+        floor_smoothed,
+        floor_tight_scan,
+        stall_8_aitken.0,
+        stall_8_iqn.0,
+        stall_2_aitken.0,
+        stall_2_iqn.0,
+    ] {
+        assert!(value.is_finite() && value > 0.0);
+    }
+
+    // Pins are for the default configuration only.
+    if ny == 62 && base_subcycle == 8 && (t_probe - 4.0).abs() < 1e-9 && scan_passes == 40 {
+        // The unsmoothed floor at these loads, measured 2026-08-21 as
+        // 3.05e-5. The band is generous (the maximum of 40 samples of a
+        // jump process moves between platforms); leaving it is a
+        // material change to the pass map's continuity either way.
+        assert!(
+            (1.0e-5..8.0e-5).contains(&floor_base),
+            "the subcycle-8 noise floor {floor_base:.3e} left its measured \
+             band [1e-5, 8e-5] — the pass map's continuity changed"
+        );
+        // Attribution: smoothing the sampled tractions does NOT move the
+        // floor (measured ratio 1.002) — the flip noise is coherent
+        // through the fluid field. If this ratio ever leaves [0.5, 2],
+        // the noise has moved into the sampling channel and the
+        // smoothing lever is worth revisiting.
+        let attribution = floor_smoothed / floor_base;
+        assert!(
+            (0.5..2.0).contains(&attribution),
+            "smoothing changed the floor by x{attribution:.2} — the noise \
+             channel attribution (coherent-through-the-fluid) no longer holds"
+        );
+        // The stall measurement, 2026-08-21: on a clean release step the
+        // subiteration converges DEEP at both subcycles — s8 aitken
+        // 3.4e-9 / iqn 1.6e-9 (12-pass budget), s2 aitken 6.5e-10 / iqn
+        // 6.3e-10 (converged below 1e-9 in 5-6 passes). The flip-scan
+        // jumps are events at specific geometries, not a floor under
+        // every step: a typical step's map is locally smooth far below
+        // any march tolerance, and the tenth session's subcycle-2 blowup
+        // (tolerance 2e-4 held fixed as dt_c shrank, wall-velocity noise
+        // = tol / dt_c) was a tolerance mis-budgeting, not an
+        // impassable floor. If any of these stalls rises above 1e-7,
+        // the step map's local smoothness is gone and the tight-coupling
+        // tolerance budget must be re-measured before the next ladder.
+        for (label, value) in [
+            ("s8 aitken", stall_8_aitken.0),
+            ("s8 iqn", stall_8_iqn.0),
+            ("s2 aitken", stall_2_aitken.0),
+            ("s2 iqn", stall_2_iqn.0),
+        ] {
+            assert!(
+                value < 1e-7,
+                "{label} stall {value:.3e} rose above 1e-7 — the release \
+                 step's local smoothness is gone; re-measure the \
+                 tight-coupling tolerance budget"
+            );
+        }
+    }
+}

crates/specialized/rtx-fsi/tests/turek_hron_fsi2.rs

@@ -24,6 +24,11 @@
 //! fluid keeps TVD convection (shedding physics; limiter chatter is
 //! harmless in time marching) and the multigrid projection.
 //!
+//! The geometry, load sampling (spike clamp + optional surface
+//! smoothing) and fluid configuration live in `fsi2_harness/`; the
+//! interface-noise-floor probe `fsi2_interface_noise.rs` measures the
+//! same machinery's pass-to-pass continuity.
+//!
 //! # Phases (the validation ladder inside FSI2)
 //!
 //! 1. **Rigid flag** to `t_release`: the ramped inflow over the fixed
@@ -90,42 +95,24 @@
 //! `RTX_FSI2_RTOL` (interface tolerance floor and its
 //! relative-to-increment part), `RTX_FSI2_MAXSUB` (Aitken budget,
 //! default 12), `RTX_FSI2_FLAG_NX` (flag mesh, default 35),
+//! `RTX_FSI2_SMOOTH` (traction smoothing radius in multiples of the cell
+//! size, default 0 = off), `RTX_FSI2_COUPLER` (`aitken` default, or `iqn`
+//! for IQN-ILS with `RTX_FSI2_REUSE` steps of secant history, default 2),
 //! `RTX_FSI2_CSV` (trajectory dump path).
 
+mod fsi2_harness;
+
 use std::cell::RefCell;
 use std::io::Write as _;
-use std::sync::{Arc, RwLock};
 
+use fsi2_harness::{Fsi2Harness, crossing_frequency, env_or, mid_amp};
 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::analysis::{
     AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
 };
-use rtx_fea::assembly::dof_mapping::DofComponent;
-use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
-use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
 use rtx_fea::materials::{LinearElastic, MaterialDatabase};
-use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use rtx_fea::mesh::{MaterialId, NodeId};
-use rtx_fsi::{FluidFace, Subiterated, WettedSurface};
+use rtx_fsi::{IqnIls, Subiterated};
-
-const L: f64 = 2.5;
-const H: f64 = 0.41;
-const RHO_F: f64 = 1000.0;
-const NU_F: f64 = 1e-3;
-const U_MEAN: f64 = 1.0;
-const RHO_S: f64 = 10_000.0;
-const E_S: f64 = 1.4e6;
-const NU_S: f64 = 0.4;
-
-const FLAG_X0: f64 = 0.25;
-const FLAG_X1: f64 = 0.6;
-const FLAG_Y0: f64 = 0.19;
-const FLAG_Y1: f64 = 0.21;
 
 // FEATFLOW level-4, dt 0.001 reference values.
 const REF_UY_MEAN: f64 = 1.30e-3;
@@ -136,221 +123,6 @@ const REF_UX_AMP: f64 = 12.70e-3;
 const REF_DRAG_MEAN: f64 = 215.06;
 const REF_LIFT_AMP: f64 = 237.8;
 
-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.
-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)
-}
-
-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).
-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).
-struct Interface {
-    wetted: Vec<NodeId>,
-    reference: Vec<(f64, f64)>,
-    /// Ordered boundary walk: indices into `wetted` (`usize::MAX` marks
-    /// the fixed anchor vertices inside the cylinder / at the clamp).
-    walk: Vec<(usize, (f64, f64))>,
-}
-
-impl Interface {
-    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`.
-    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).
-    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.
-    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()
-    }
-}
-
-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
-}
-
-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.
-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()))
-}
-
 #[test]
 #[allow(clippy::too_many_lines)]
 fn fsi2_flapping_flag() {
@@ -380,92 +152,16 @@ fn fsi2_flapping_flag() {
     let rtol = env_or("RTX_FSI2_RTOL", 1e-2);
     let max_subiterations_budget = env_or("RTX_FSI2_MAXSUB", 12.0) as usize;
     let flag_nx = env_or("RTX_FSI2_FLAG_NX", 35.0) as usize;
+    let smooth_in_h = env_or("RTX_FSI2_SMOOTH", 0.0);
     let csv_path = std::env::var("RTX_FSI2_CSV").ok();
 
-    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 — the structure and
-    // the transfer need nowhere near the fluid's dt (CSM3 measured 0.23%
-    // frequency error at dt = 5e-3; dt_c here is ~2.6e-3 at ny = 62), and
-    // the per-pass cost is dominated by the structure solve. Within a
-    // pass the interface geometry is interpolated linearly across the
-    // substeps, so the mask still moves less than a cell per fluid step.
-    let dt_fluid = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
     let subcycle = env_or("RTX_FSI2_SUBCYCLE", 8.0) as usize;
+
+    let (harness, mut solver, mut field) = Fsi2Harness::build(ny, flag_nx, smooth_in_h);
+    let dt_fluid = harness.dt_fluid;
     let dt = dt_fluid * subcycle as f64;
-
+    let interface = &harness.interface;
-    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");
-
-    // 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.
     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();
 
     // Phase 1: rigid flag to t_release.
     let start = std::time::Instant::now();
@@ -473,43 +169,10 @@ fn fsi2_flapping_flag() {
     for _ in 0..rigid_steps {
         futures::executor::block_on(solver.advance(&mut field, dt_fluid)).unwrap();
     }
-    // Surface drag and lift on cylinder + flag at the current geometry.
-    let measure_force = |solver: &EmbeddedPisoSolver, field: &FlowField| -> (f64, f64) {
-        let mask = solver.mask().unwrap();
-        let body = solver.body().unwrap();
-        let vertices = 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 * 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, 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 * 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, mu, 0.0, s.x, s.y, s.nx, s.ny,
-            ) {
-                drag += tx * s.ds;
-                lift += ty * s.ds;
-            }
-        }
-        (drag, lift)
-    };
     // The fluid harness check: surface drag on cylinder + flag near the
     // CFD2 value this solver measured on this geometry (ny = 62: ~121;
     // the reference is 136.700 with the boundary layer barely a cell).
-    let (rigid_drag, rigid_lift) = measure_force(&solver, &field);
+    let (rigid_drag, rigid_lift) = harness.measure_force(&solver, &field);
     println!(
         "  rigid phase: {rigid_steps} steps to t = {t_release:.1} s in {:.0} s wall; \
          surface drag {rigid_drag:.1} (CFD2 ref 136.7, this grid measured ~121), \
@@ -521,7 +184,7 @@ fn fsi2_flapping_flag() {
     let mut db = MaterialDatabase::new();
     db.add_material(
         MaterialId(0),
-        LinearElastic::new(E_S, NU_S).with_density(RHO_S),
+        LinearElastic::new(fsi2_harness::E_S, fsi2_harness::NU_S).with_density(fsi2_harness::RHO_S),
         None,
     );
     // A deep Newton budget: a mid-swing subiteration can hand the flag a
@@ -529,9 +192,9 @@ fn fsi2_flapping_flag() {
     // within a period); typical steps converge in 1-2 iterations, and a
     // t = 25.8 s failure at the default budget of 25 is what set this.
     let analysis = NonlinearDynamicAnalysis::new(
-        mesh.clone(),
+        harness.mesh.clone(),
         db,
-        clamp_left(&mesh),
+        fsi2_harness::clamp_left(&harness.mesh),
         dt,
         1,
         AnalysisConfig::default(),
@@ -550,7 +213,7 @@ fn fsi2_flapping_flag() {
             [dofs[0], dofs[1]]
         })
         .collect();
-    let a_dofs = flag.borrow().node_dofs(a_node);
+    let a_dofs = flag.borrow().node_dofs(harness.a_node);
     let extract = |state: &DynamicState| -> Vec<f64> {
         let mut d = vec![0.0; 2 * wetted_dofs.len()];
         for (k, dofs) in wetted_dofs.iter().enumerate() {
@@ -560,88 +223,10 @@ fn fsi2_flapping_flag() {
         d
     };
 
-    // 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).
-    let spiked_total = std::cell::Cell::new(0usize);
-    let sample_load = |solver: &EmbeddedPisoSolver,
-                       field: &FlowField,
-                       d: &[f64]|
-     -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
-        let vertices = 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 * 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, 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;
-                    spiked_total.set(spiked_total.get() + 1);
-                }
-            }
-        }
-        let nodes_now = 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);
-        (
-            interface
-                .wetted
-                .iter()
-                .zip(nodal)
-                .map(|(&id, f)| (id, f))
-                .collect(),
-            conservation,
-            skipped,
-        )
-    };
-
     // Phase 2: release. The flag starts at rest under the current fluid
     // load (consistent initial acceleration — the step response about the
     // steady deflection is the seed perturbation for the instability).
-    let (nodal0, conservation0, _) = sample_load(&solver, &field, &zero_d);
+    let (nodal0, conservation0, _) = harness.sample_load(&solver, &field, &zero_d);
     flag.borrow_mut().set_nodal_forces(&nodal0);
     let mut flag_state = flag.borrow_mut().rest_state().unwrap();
     let mut committed_nodal = nodal0;
@@ -650,6 +235,16 @@ fn fsi2_flapping_flag() {
     let solver = RefCell::new(solver);
     let field = RefCell::new(field);
 
+    // The interface driver: per-step Aitken (the tenth-session default),
+    // or a persistent IQN-ILS whose secant history carries across steps.
+    let coupler_kind = std::env::var("RTX_FSI2_COUPLER").unwrap_or_else(|_| "aitken".into());
+    let reuse = env_or("RTX_FSI2_REUSE", 2.0) as usize;
+    let mut iqn = (coupler_kind == "iqn").then(|| {
+        IqnIls::new(max_subiterations_budget, 1.0)
+            .unwrap()
+            .with_reuse(reuse)
+    });
+
     let coupled_steps = ((t_end - t_release) / dt).round() as usize;
     let mut times = Vec::with_capacity(coupled_steps);
     let mut ux_series = Vec::with_capacity(coupled_steps);
@@ -676,7 +271,7 @@ fn fsi2_flapping_flag() {
         let fluid_saved = solver.borrow().snapshot();
         let field_saved = field.borrow().clone();
         type PassResult = (
-            FlowField,
+            rtx_cfd::solvers::incompressible::FlowField,
             DynamicState,
             Vec<(NodeId, Vector3<f64>)>,
             f64,
@@ -685,36 +280,23 @@ fn fsi2_flapping_flag() {
         let latest: RefCell<Option<PassResult>> = RefCell::new(None);
 
         let pass = |d_candidate: &[f64]| -> Vec<f64> {
-            // Interface velocity of THIS candidate, constant over the step.
-            let ddot: Vec<f64> = d_candidate
-                .iter()
-                .zip(&d_n)
-                .map(|(new, old)| (new - old) / dt)
-                .collect();
             // Subcycled fluid steps from the SAME start-of-step state,
-            // geometry interpolated to each substep's end time.
+            // geometry interpolated to each substep's end time, interface
+            // velocity of THIS candidate constant over the step.
             let mut solver_ref = solver.borrow_mut();
             solver_ref.restore(&fluid_saved);
             let mut trial_field = field_saved.clone();
-            for m in 1..=subcycle {
+            harness.advance_subcycled(
-                let fraction = m as f64 / subcycle as f64;
+                &mut solver_ref,
-                let d_sub: Vec<f64> = d_n
+                &mut trial_field,
-                    .iter()
+                &d_n,
-                    .zip(d_candidate)
+                d_candidate,
-                    .map(|(old, new)| old + fraction * (new - old))
+                subcycle,
-                    .collect();
+            );
-                {
-                    let mut geometry = shared.write().unwrap();
-                    geometry.0 = interface.polygon(&d_sub);
-                    geometry.1 = interface.walk_velocities(&ddot);
-                }
-                futures::executor::block_on(solver_ref.advance(&mut trial_field, dt_fluid))
-                    .unwrap();
-            }
             // Load on the candidate geometry, flag answers from the
             // committed state.
             let (nodal, conservation, skipped) =
-                sample_load(&solver_ref, &trial_field, d_candidate);
+                harness.sample_load(&solver_ref, &trial_field, d_candidate);
             let mut flag_ref = flag.borrow_mut();
             flag_ref.set_nodal_forces(&nodal);
             let (candidate_state, _) = flag_ref.step(&flag_state).unwrap();
@@ -731,8 +313,15 @@ fn fsi2_flapping_flag() {
             .sum::<f64>()
             .sqrt();
         let tol_step = tol_floor.max(rtol * increment);
-        let mut scheme = Subiterated::aitken(max_subiterations_budget, tol_step).unwrap();
+        let outcome = if let Some(iqn) = iqn.as_mut() {
-        match scheme.solve(&d_predicted, pass) {
+            iqn.set_tolerance(tol_step).unwrap();
+            iqn.solve(&d_predicted, pass)
+        } else {
+            Subiterated::aitken(max_subiterations_budget, tol_step)
+                .unwrap()
+                .solve(&d_predicted, pass)
+        };
+        match outcome {
             Ok(converged) => {
                 total_subiterations += converged.iterations;
                 max_subiterations = max_subiterations.max(converged.iterations);
@@ -775,7 +364,7 @@ fn fsi2_flapping_flag() {
         ux_series.push(ux);
         uy_series.push(uy);
         if (step + 1) % 10 == 0 {
-            let (drag, lift) = measure_force(&solver.borrow(), &field.borrow());
+            let (drag, lift) = harness.measure_force(&solver.borrow(), &field.borrow());
             force_times.push(t);
             drag_series.push(drag);
             lift_series.push(lift);
@@ -839,7 +428,7 @@ fn fsi2_flapping_flag() {
          (ref {:.2} ± {:.1}), ux(A) = {:.4} ± {:.4} mm (ref {:.2} ± {:.2}), f = {} Hz \
          (ref {REF_UY_FREQ}); onset amp {:.3e} -> {:.3e} m",
         elapsed,
-        spiked_total.get(),
+        harness.spiked_total.get(),
         t_window.first().unwrap_or(&t_release),
         uy_mid * 1e3,
         uy_amp * 1e3,
@@ -883,7 +472,16 @@ fn fsi2_flapping_flag() {
     // 81.6 mm is not reached by this coupling; the measured state is the
     // wake-forced 3.73 Hz / ±17.3 mm cycle at BOTH grids). If a change
     // moves these numbers, that is a finding either way and must be loud.
-    if ny == 62 && flag_nx == 35 && (t_end - 7.0).abs() < 1e-9 && (t_release - 6.0).abs() < 1e-9 {
+    // The bands pin the UNSMOOTHED, Aitken-coupled sampling (the
+    // defaults): smoothing or the IQN coupler change the load path and
+    // re-pin deliberately.
+    let default_coupling = smooth_in_h == 0.0 && coupler_kind == "aitken";
+    if default_coupling
+        && ny == 62
+        && flag_nx == 35
+        && (t_end - 7.0).abs() < 1e-9
+        && (t_release - 6.0).abs() < 1e-9
+    {
         // The committed default: the release response over [6, 7] s,
         // measured 2026-08-21 as uy mid 3.773 mm, amp 3.792 mm. The band
         // is ±35% for cross-platform floating-point drift in a growing
@@ -898,7 +496,7 @@ fn fsi2_flapping_flag() {
             "uy release-response amp {uy_amp:.4e} outside the measured band \
              [2.4e-3, 5.2e-3]"
         );
-    } else if t_end >= 25.0 {
+    } else if default_coupling && t_end >= 25.0 {
         // Study horizons: the measured attractor of the loosely-coupled
         // (subcycle 8) march — f = 3.729 / 3.728 Hz and uy amp 17.3 mm at
         // ny = 62 / 82 (2026-08-21).