Merge r8a-span-uniform-fsi (R8-a: the first coupled 3D FSI, span-uniform; default-off). Merge fixes: the load-route observers of R8-a (to_load_sink) and R8-c (sink(LoadKind, …)) both kept on the same summands (both inert by default); the harness's DeviceSdf gains plate: None (R8-c's new field). Merged-tree gate: slab ny 62 default / BAND_CHECK / all-host / uniform plate byte-identical to main's base; device suites 2/2, 3/3; the coupled harness reproduces R8-a's recorded slab march row for row (478 rows)

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
2026-09-25 21:21:40 -05:00
parent fc73929f8e 11fff74c3d
commit 171da41ed1
8 changed files with 1376 additions and 0 deletions
@@ -0,0 +1,685 @@
+//! R8-a: the first coupled 3D FSI — the embedded3 cut-cell fluid (device
+//! path) coupled to the 2D Turek–Hron flag (rtx-fea 35×2 Quad8 SVK, total
+//! Lagrangian, Newmark γ 0.7) under a SPAN-UNIFORM deformation: the
+//! structure's centreline drives the 3D body's polyline each coupled step,
+//! the fluid's operator-route wall load, integrated over the span per unit
+//! width, loads the structure's wetted nodes. Partitioned: per coupled
+//! step the fluid re-runs the same step from a device snapshot for each
+//! Aitken subiteration (the overset harness's pattern; `DeviceStep::
+//! snapshot / restore`).
+//!
+//! The body is the embedded flag test's capsule (a semicircular tip, apex
+//! on A): its 2D counterpart is the overset's SEMICIRCLE line (P5-2 ny 62:
+//! 94.8 mm at 1.914 Hz), not the flat/1.25 mm-corner reference line.
+//!
+//! Knobs `RTX_E3FSI_*`: `NY` (62), `NZ` (4 = the periodic slab; 0 = the
+//! full 0.41 m duct with slip sides), `T_RIGID` (3.0 s of rigid flag),
+//! `T_END` (13.0), `RIGID_ONLY` (1 = stop after the rigid phase: target 1),
+//! `RTOL` (1e-3), `FLOOR` (1.5e-7 on the centreline vector), `MAX_SUBIT`
+//! (12), `STALL_ACCEPT` (5), `GAMMA` (0.7), `SPEED` (3.0 m/s, the band's
+//! surface-speed bound), `CSV` (per-step series), `TRACE` (steps whose
+//! passes are printed).
+//!
+//! `RTX_E3FSI_NY=62 RTX_E3FSI_CSV=<path> RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-fsi \
+//!  --features cuda --test fsi2_embedded3 -- --ignored --nocapture`
+#![cfg(feature = "cuda")]
+
+#[path = "fsi2_embedded3/fluid.rs"]
+mod fluid;
+#[path = "fsi2_harness/mod.rs"]
+mod fsi2_harness;
+#[path = "fsi2_embedded3/state.rs"]
+mod state;
+
+use std::cell::RefCell;
+use std::io::Write as _;
+
+use fluid::{CX, CY, Contribution, E3Fluid, HALF, Line, R_CYL};
+use fsi2_harness::{FSI2, Interface, clamp_left, flag_mesh, median, mid_amp};
+use nalgebra::{DVector, Vector3};
+use rtx_fea::analysis::{
+    AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
+};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::{MaterialId, Mesh, NodeId};
+use rtx_fsi::Subiterated;
+
+pub fn env_f(name: &str, default: f64) -> f64 {
+    std::env::var(name)
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(default)
+}
+
+const X0: f64 = 0.25;
+/// Centreline stations: the element corners at x 0.25, 0.26, …, 0.59 (the
+/// last 10 mm is the capsule's cap).
+const STATIONS: usize = 35;
+
+/// The flag's structure-side bookkeeping: the centreline nodes, the wetted
+/// edges with their reference coordinates, the tip node A.
+struct Flag {
+    centre: Vec<NodeId>,
+    interface: Interface,
+    /// (reference x, wetted index or None for the clamp corner), ascending.
+    bottom: Vec<(f64, Option<usize>)>,
+    top: Vec<(f64, Option<usize>)>,
+    /// (reference y, wetted index), ascending (corners included).
+    tip: Vec<(f64, usize)>,
+    a_node: NodeId,
+}
+
+impl Flag {
+    fn build(mesh: &Mesh) -> Self {
+        let find = |x: f64, y: f64| -> NodeId {
+            *mesh
+                .nodes
+                .iter()
+                .find(|(_, n)| {
+                    let p = n.position();
+                    (p.x - x).abs() < 1e-9 && (p.y - y).abs() < 1e-9
+                })
+                .expect("node")
+                .0
+        };
+        let centre = (0..STATIONS)
+            .map(|k| find(X0 + 0.01 * k as f64, 0.2))
+            .collect();
+        let interface = Interface::build(mesh);
+        let edge = |idx: &[usize]| -> Vec<(f64, Option<usize>)> {
+            let mut v: Vec<(f64, Option<usize>)> = vec![(X0, None)];
+            v.extend(idx.iter().map(|&k| (interface.reference[k].0, Some(k))));
+            v.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
+            v
+        };
+        let bottom = edge(&interface.bottom);
+        let top = edge(&interface.top);
+        let mut tip: Vec<(f64, usize)> = interface
+            .tip
+            .iter()
+            .map(|&k| (interface.reference[k].1, k))
+            .collect();
+        tip.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
+        Self {
+            centre,
+            interface,
+            bottom,
+            top,
+            tip,
+            a_node: find(0.6, 0.2),
+        }
+    }
+}
+
+/// Linear split of `value` at `s` over sorted stations; the share of a
+/// `None` station (the clamp) is dropped.
+fn split<T: Copy>(stations: &[(f64, T)], s: f64, value: f64, mut add: impl FnMut(T, f64)) {
+    let n = stations.len();
+    let s = s.clamp(stations[0].0, stations[n - 1].0);
+    let m = stations.partition_point(|st| st.0 < s).clamp(1, n - 1);
+    let (a, b) = (stations[m - 1], stations[m]);
+    let u = if b.0 > a.0 {
+        (s - a.0) / (b.0 - a.0)
+    } else {
+        0.0
+    };
+    add(a.1, (1.0 - u) * value);
+    add(b.1, u * value);
+}
+
+/// The operator route's contributions onto the wetted nodes (per unit
+/// span): each is placed by its closest point on the fluid's centreline —
+/// the top or bottom edge at the same arc position, or the tip edge by its
+/// lateral offset beyond the last station. Returns the nodal loads, the
+/// flag's (fx, fy) and the cylinder's (fx, fy) per span, and the flag's fy
+/// by part (pressure, shear, diffusive and convective exchange).
+/// `RTX_E3FSI_PARTS` (bit mask, default 15 = all) keeps parts off the
+/// structure (a diagnostic; the reported loads keep every part).
+fn distribute(
+    flag: &Flag,
+    line: &Line,
+    contributions: &[Contribution],
+    width: f64,
+) -> (Vec<(NodeId, Vector3<f64>)>, [f64; 2], [f64; 2], [f64; 4]) {
+    let nw = flag.interface.wetted.len();
+    let mut f = vec![[0.0f64; 2]; nw];
+    let (mut on_flag, mut on_cyl) = ([0.0f64; 2], [0.0f64; 2]);
+    let mut parts_y = [0.0f64; 4];
+    let pts = &line.pts;
+    let mut cum = vec![0.0; pts.len()];
+    for m in 1..pts.len() {
+        cum[m] = cum[m - 1]
+            + ((pts[m][0] - pts[m - 1][0]).powi(2) + (pts[m][1] - pts[m - 1][1]).powi(2)).sqrt();
+    }
+    let parts_on = env_f("RTX_E3FSI_PARTS", 15.0) as usize;
+    for &(pos, c, part, v) in contributions {
+        if c > 1 {
+            continue;
+        }
+        let v = v / width;
+        let (x, y) = (pos[0], pos[1]);
+        let mut best = (f64::INFINITY, 0usize, 0.0f64);
+        for m in 0..pts.len() - 1 {
+            let (a, b) = (pts[m], pts[m + 1]);
+            let (ex, ey) = (b[0] - a[0], b[1] - a[1]);
+            let l2 = ex * ex + ey * ey;
+            let u = (((x - a[0]) * ex + (y - a[1]) * ey) / l2).clamp(0.0, 1.0);
+            let d = ((x - a[0] - u * ex).powi(2) + (y - a[1] - u * ey).powi(2)).sqrt();
+            if d < best.0 {
+                best = (d, m, u);
+            }
+        }
+        let d_cyl = ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL;
+        if d_cyl < best.0 - HALF {
+            on_cyl[c] += v;
+            continue;
+        }
+        on_flag[c] += v;
+        if c == 1 {
+            parts_y[part] += v;
+        }
+        if parts_on & (1 << part) == 0 {
+            continue;
+        }
+        let (_, m, u) = best;
+        let (a, b) = (pts[m], pts[m + 1]);
+        let (ex, ey) = (b[0] - a[0], b[1] - a[1]);
+        let len = (ex * ex + ey * ey).sqrt();
+        let (tx, ty) = (ex / len, ey / len);
+        let (px, py) = (a[0] + u * ex, a[1] + u * ey);
+        // Lateral offset: + on the upper side of the centreline.
+        let eta = tx * (y - py) - ty * (x - px);
+        let along = tx * (x - px) + ty * (y - py);
+        let mut add = |k: usize, w: f64| f[k][c] += w;
+        if m + 2 == pts.len() && u >= 1.0 && along > 0.0 {
+            let yr = 0.2 + eta.clamp(-HALF, HALF);
+            split(&flag.tip, yr, v, &mut add);
+        } else {
+            let xr = X0 + cum[m] + u * len;
+            let edge = if eta >= 0.0 { &flag.top } else { &flag.bottom };
+            split(edge, xr, v, |k, w| {
+                if let Some(k) = k {
+                    add(k, w);
+                }
+            });
+        }
+    }
+    let nodal = flag
+        .interface
+        .wetted
+        .iter()
+        .zip(&f)
+        .map(|(&id, v)| (id, Vector3::new(v[0], v[1], 0.0)))
+        .collect();
+    (nodal, on_flag, on_cyl, parts_y)
+}
+
+/// The fluid's centreline from the structure's centreline displacement `c`
+/// (2 per station) with velocities `(c − c_prev) / dt`.
+fn line_of(t: f64, c: &[f64], c_prev: &[f64], dt: f64) -> Line {
+    let pts = (0..STATIONS)
+        .map(|k| [X0 + 0.01 * k as f64 + c[2 * k], 0.2 + c[2 * k + 1]])
+        .collect();
+    let vel = (0..STATIONS)
+        .map(|k| {
+            [
+                (c[2 * k] - c_prev[2 * k]) / dt,
+                (c[2 * k + 1] - c_prev[2 * k + 1]) / dt,
+            ]
+        })
+        .collect();
+    Line { t, pts, vel }
+}
+
+#[test]
+#[ignore = "R8-a: the coupled FSI2 on embedded3 (GPU; minutes to hours)"]
+fn fsi2_on_embedded3() {
+    let ny = env_f("RTX_E3FSI_NY", 62.0) as usize;
+    let nz = env_f("RTX_E3FSI_NZ", 4.0) as usize;
+    let t_rigid = env_f("RTX_E3FSI_T_RIGID", 3.0);
+    let t_end = env_f("RTX_E3FSI_T_END", 13.0);
+    let rigid_only = env_f("RTX_E3FSI_RIGID_ONLY", 0.0) > 0.5;
+    let rtol = env_f("RTX_E3FSI_RTOL", 1e-3);
+    let floor = env_f("RTX_E3FSI_FLOOR", 1.5e-7);
+    let max_subit = env_f("RTX_E3FSI_MAX_SUBIT", 12.0) as usize;
+    let stall_accept = env_f("RTX_E3FSI_STALL_ACCEPT", 5.0);
+    let gamma = env_f("RTX_E3FSI_GAMMA", 0.7);
+    let speed = env_f("RTX_E3FSI_SPEED", 3.0);
+    let trace = env_f("RTX_E3FSI_TRACE", 0.0) as usize;
+    let csv_path = std::env::var("RTX_E3FSI_CSV").ok();
+    let restore_check = env_f("RTX_E3FSI_RESTORE_CHECK", 0.0) as usize;
+    let case = FSI2;
+
+    let mesh = flag_mesh(35, 2);
+    let flag_geo = Flag::build(&mesh);
+    let zero_c = vec![0.0; 2 * STATIONS];
+    let rest = line_of(0.0, &zero_c, &zero_c, 1.0);
+    // `RTX_E3FSI_LOAD=<dir>`: continue from a saved coupled state (extruded
+    // onto the full duct when the saved nz differs); `RTX_E3FSI_SAVE=<dir>`
+    // saves the state every `RTX_E3FSI_SAVE_EVERY` coupled steps (2000) and at the end.
+    let saved = std::env::var("RTX_E3FSI_LOAD")
+        .ok()
+        .map(|d| state::Saved::load(&d).expect("load the saved state"));
+    let save_dir = std::env::var("RTX_E3FSI_SAVE").ok();
+    let save_every = env_f("RTX_E3FSI_SAVE_EVERY", 2000.0) as usize;
+    let mut fl = E3Fluid::build(ny, nz, speed, rest.clone(), saved.as_ref());
+    if let Some(s) = &saved {
+        println!(
+            "  loaded the coupled state at t {:.4} from {} ({}×{}×{} → nz {})",
+            s.t,
+            std::env::var("RTX_E3FSI_LOAD").unwrap(),
+            s.dims[0],
+            s.dims[1],
+            s.dims[2],
+            fl.grid.nz
+        );
+    }
+    let dt = fl.dt;
+    println!(
+        "  R8-a FSI2 on embedded3: rigid to {t_rigid} s, coupled to {t_end} s; Aitken rtol {rtol:.1e} floor {floor:.1e} max {max_subit} stall accept {stall_accept}; Newmark γ {gamma}; the body's 2D counterpart = the overset SEMICIRCLE line (ny 62: 94.8 mm, 1.914 Hz)"
+    );
+    let mut csv = csv_path.as_ref().map(|p| {
+        let mut f = std::fs::File::create(p).expect("csv");
+        writeln!(
+            f,
+            "t,phase,ux,uy,drag,lift,drag_flag,lift_flag,subit,dres,residual,cg,fresh,sink_dx,sink_dy,lift_p,lift_shear,lift_xdiff,lift_xconv"
+        )
+        .unwrap();
+        f
+    });
+    let start = std::time::Instant::now();
+
+    // Phase 1: the rigid flag (target 1: the rest state vs CFD2 136.7 / 10.53).
+    let rigid_steps = if saved.is_some() {
+        0
+    } else {
+        (t_rigid / dt).round() as usize
+    };
+    let mut last = ([0.0; 3], Vec::new());
+    for step in 0..rigid_steps {
+        let r = fl.step();
+        assert!(r.final_residual.is_finite(), "rigid death at step {step}");
+        if (step + 1) % 50 == 0 || step + 1 == rigid_steps {
+            last = fl.loads();
+            let (tot, contrib) = &last;
+            let (_, on_flag, on_cyl, _) = distribute(&flag_geo, &rest, contrib, fl.load_width);
+            let t = fl.time();
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{t:.6},rigid,0,0,{:.5},{:.5},{:.5},{:.5},0,0,{:.3e},{},{},{:.3e},{:.3e},,,,",
+                    tot[0],
+                    tot[1],
+                    on_flag[0],
+                    on_flag[1],
+                    r.final_residual,
+                    r.poisson_iterations,
+                    r.fresh_cells,
+                    on_flag[0] + on_cyl[0] - tot[0],
+                    on_flag[1] + on_cyl[1] - tot[1]
+                )
+                .unwrap();
+            }
+            if (step + 1) % 500 == 0 || step + 1 == rigid_steps {
+                println!(
+                    "  rigid t {t:.3}: drag/span {:.2} lift/span {:+.2} (flag {:.2} {:+.2}, cylinder {:.2} {:+.2}); residual {:.1e}, CG {}; [{:.0} s]",
+                    tot[0],
+                    tot[1],
+                    on_flag[0],
+                    on_flag[1],
+                    on_cyl[0],
+                    on_cyl[1],
+                    r.final_residual,
+                    r.poisson_iterations,
+                    start.elapsed().as_secs_f64()
+                );
+            }
+        }
+    }
+    let (tot0, contrib0) = last;
+    println!(
+        "  RIGID ny {ny} nz {nz} at t {:.3}: drag/span {:.2} (CFD2 136.7), lift/span {:+.2} (CFD2 10.53); {rigid_steps} steps in {:.0} s",
+        fl.time(),
+        tot0[0],
+        tot0[1],
+        start.elapsed().as_secs_f64()
+    );
+    if rigid_only && saved.is_none() {
+        return;
+    }
+
+    // The structure: FSI2's flag at the coupled step dt.
+    let mut db = MaterialDatabase::new();
+    db.add_material(
+        MaterialId(0),
+        LinearElastic::new(case.e_s, case.nu_s).with_density(case.rho_s),
+        None,
+    );
+    let beta = (gamma + 0.5).powi(2) / 4.0;
+    let analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        db,
+        clamp_left(&mesh),
+        dt,
+        1,
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian()
+    .with_convergence_criteria(ConvergenceCriteria {
+        max_iterations: 60,
+        ..ConvergenceCriteria::default()
+    })
+    .with_newmark_parameters(gamma, beta);
+    let flag = RefCell::new(analysis.stepper().unwrap());
+    let centre_dofs: Vec<[usize; 2]> = flag_geo
+        .centre
+        .iter()
+        .map(|&id| {
+            let d = flag.borrow().node_dofs(id);
+            [d[0], d[1]]
+        })
+        .collect();
+    let a_dofs = flag.borrow().node_dofs(flag_geo.a_node);
+    let extract = |s: &DynamicState| -> Vec<f64> {
+        let mut c = vec![0.0; 2 * STATIONS];
+        for (k, d) in centre_dofs.iter().enumerate() {
+            c[2 * k] = s.displacement[d[0]];
+            c[2 * k + 1] = s.displacement[d[1]];
+        }
+        c
+    };
+
+    let (mut flag_state, mut committed_nodal, mut line_n, mut c_fluid_n) = match &saved {
+        Some(s) => {
+            let nodal: Vec<(NodeId, Vector3<f64>)> = flag_geo
+                .interface
+                .wetted
+                .iter()
+                .enumerate()
+                .map(|(k, &id)| (id, Vector3::new(s.nodal[2 * k], s.nodal[2 * k + 1], 0.0)))
+                .collect();
+            let state = DynamicState {
+                displacement: DVector::from_vec(s.disp.clone()),
+                velocity: DVector::from_vec(s.vel.clone()),
+                acceleration: DVector::from_vec(s.acc.clone()),
+            };
+            (state, nodal, s.line.clone(), s.c_fluid.clone())
+        }
+        None => {
+            let (nodal0, _, _, _) = distribute(&flag_geo, &rest, &contrib0, fl.load_width);
+            flag.borrow_mut().set_nodal_forces(&nodal0);
+            let state = flag.borrow_mut().rest_state().unwrap();
+            // The fluid's own previous line and centreline (its geometry's history).
+            let line = Line {
+                t: fl.time(),
+                ..rest.clone()
+            };
+            (state, nodal0, line, zero_c.clone())
+        }
+    };
+
+    let coupled_steps = ((t_end - fl.time()) / dt).round() as usize;
+    let fl = RefCell::new(fl);
+    let (mut times, mut uy_series, mut ux_series) = (Vec::new(), Vec::new(), Vec::new());
+    let (mut drag_s, mut lift_s) = (Vec::new(), Vec::new());
+    let (mut total_subit, mut max_seen, mut stalled) = (0usize, 0usize, 0usize);
+    let mut death: Option<String> = None;
+    let t_fluid = std::cell::Cell::new(0.0f64);
+    let t_restore = std::cell::Cell::new(0.0f64);
+    let t_loads = std::cell::Cell::new(0.0f64);
+    let t_struct = std::cell::Cell::new(0.0f64);
+    let phase_start = std::time::Instant::now();
+    for step in 0..coupled_steps {
+        let t_old = fl.borrow().time();
+        let t_new = t_old + dt;
+        let predicted = {
+            flag.borrow_mut().set_nodal_forces(&committed_nodal);
+            let (p, _) = flag.borrow_mut().step(&flag_state).unwrap();
+            extract(&p)
+        };
+        let c_struct_n = extract(&flag_state);
+        let snap = fl.borrow_mut().snapshot();
+        let dirty = std::cell::Cell::new(false);
+        type Pass = (
+            DynamicState,
+            Vec<(NodeId, Vector3<f64>)>,
+            Line,
+            Vec<f64>,
+            [f64; 3],
+            [f64; 2],
+            fluid_step::Stats,
+        );
+        let latest: RefCell<Option<Pass>> = RefCell::new(None);
+        let pass = |cand: &[f64]| -> Vec<f64> {
+            let line = line_of(t_new, cand, &c_fluid_n, dt);
+            let mut f = fl.borrow_mut();
+            let tr = std::time::Instant::now();
+            f.set_lines(line_n.clone(), line.clone());
+            if dirty.get() {
+                f.restore(&snap);
+            }
+            dirty.set(true);
+            t_restore.set(t_restore.get() + tr.elapsed().as_secs_f64());
+            let tf = std::time::Instant::now();
+            let r = f.step();
+            t_fluid.set(t_fluid.get() + tf.elapsed().as_secs_f64());
+            if !r.final_residual.is_finite() {
+                return vec![f64::NAN; cand.len()];
+            }
+            let tl = std::time::Instant::now();
+            let (tot, contrib) = f.loads();
+            let (nodal, on_flag, on_cyl, parts_y) =
+                distribute(&flag_geo, &line, &contrib, f.load_width);
+            t_loads.set(t_loads.get() + tl.elapsed().as_secs_f64());
+            let ts = std::time::Instant::now();
+            let mut st = flag.borrow_mut();
+            st.set_nodal_forces(&nodal);
+            let (new_state, _) = st.step(&flag_state).unwrap();
+            t_struct.set(t_struct.get() + ts.elapsed().as_secs_f64());
+            let out = extract(&new_state);
+            if step < trace {
+                let res: f64 = out
+                    .iter()
+                    .zip(cand)
+                    .map(|(a, b)| (a - b).powi(2))
+                    .sum::<f64>()
+                    .sqrt();
+                println!(
+                    "    step {step} pass: |c_new − c_cand| {res:.3e}, tip cand ({:+.4e}, {:+.4e}), load flag ({:+.3}, {:+.3}) cyl ({:+.3}, {:+.3}) total ({:+.3}, {:+.3}), residual {:.1e}, fresh {}",
+                    cand[2 * STATIONS - 2],
+                    cand[2 * STATIONS - 1],
+                    on_flag[0],
+                    on_flag[1],
+                    on_cyl[0],
+                    on_cyl[1],
+                    tot[0],
+                    tot[1],
+                    r.final_residual,
+                    r.fresh_cells
+                );
+            }
+            let stats = fluid_step::Stats {
+                residual: r.final_residual,
+                cg: r.poisson_iterations,
+                fresh: r.fresh_cells,
+                sink: [
+                    on_flag[0] + on_cyl[0] - tot[0],
+                    on_flag[1] + on_cyl[1] - tot[1],
+                ],
+                parts_y,
+            };
+            *latest.borrow_mut() =
+                Some((new_state, nodal, line, cand.to_vec(), tot, on_flag, stats));
+            out
+        };
+        let increment: f64 = predicted
+            .iter()
+            .zip(&c_struct_n)
+            .map(|(a, b)| (a - b).powi(2))
+            .sum::<f64>()
+            .sqrt();
+        let tol = floor.max(rtol * increment);
+        let acceptable = (stall_accept * tol).max(0.1 * increment);
+        let outcome = Subiterated::aitken(max_subit, tol)
+            .unwrap()
+            .solve(&predicted, &pass);
+        // `RTX_E3FSI_RESTORE_CHECK=N`: on the first N steps, the accepted
+        // candidate re-run from the snapshot (the restore's repeatability).
+        if step < restore_check && outcome.is_ok() {
+            let (cand, tot_a, out_a) = {
+                let l = latest.borrow();
+                let l = l.as_ref().expect("a pass ran");
+                (l.3.clone(), l.4, extract(&l.0))
+            };
+            let out_b = pass(&cand);
+            let tot_b = latest.borrow().as_ref().expect("re-run").4;
+            let dout: f64 = out_a
+                .iter()
+                .zip(&out_b)
+                .map(|(a, b)| (a - b).powi(2))
+                .sum::<f64>()
+                .sqrt();
+            println!(
+                "  RESTORE CHECK step {step}: load a ({:+.9e}, {:+.9e}) b ({:+.9e}, {:+.9e}), |Δ structure output| {dout:.3e}",
+                tot_a[0], tot_a[1], tot_b[0], tot_b[1]
+            );
+        }
+        let (iters, dres) = match outcome {
+            Ok(c) => (c.iterations, c.residual),
+            Err(
+                rtx_fsi::FsiError::CouplingNotConverged {
+                    iterations,
+                    residual,
+                    ..
+                }
+                | rtx_fsi::FsiError::CouplingDiverged {
+                    iterations,
+                    residual,
+                },
+            ) if residual < acceptable => {
+                stalled += 1;
+                (iterations, residual)
+            }
+            Err(e) => {
+                println!(
+                    "  R8-a DEATH at coupled step {step} t = {t_new:.4}: {e:?} (increment {increment:.3e}, tol {tol:.3e}, acceptable {acceptable:.3e})"
+                );
+                death = Some(format!("step {step} t {t_new:.4}: {e:?}"));
+                break;
+            }
+        };
+        total_subit += iters;
+        max_seen = max_seen.max(iters);
+        let (new_state, nodal, line, cand, tot, on_flag, stats) =
+            latest.borrow_mut().take().expect("a pass ran");
+        flag_state = new_state;
+        committed_nodal = nodal;
+        line_n = line;
+        c_fluid_n = cand;
+        let ux = flag_state.displacement[a_dofs[0]];
+        let uy = flag_state.displacement[a_dofs[1]];
+        times.push(t_new);
+        ux_series.push(ux);
+        uy_series.push(uy);
+        drag_s.push(tot[0]);
+        lift_s.push(tot[1]);
+        if let Some(f) = csv.as_mut() {
+            writeln!(
+                f,
+                "{t_new:.6},coupled,{ux:.6e},{uy:.6e},{:.5},{:.5},{:.5},{:.5},{iters},{dres:.3e},{:.3e},{},{},{:.3e},{:.3e},{:.4},{:.4},{:.4},{:.4}",
+                tot[0],
+                tot[1],
+                on_flag[0],
+                on_flag[1],
+                stats.residual,
+                stats.cg,
+                stats.fresh,
+                stats.sink[0],
+                stats.sink[1],
+                stats.parts_y[0],
+                stats.parts_y[1],
+                stats.parts_y[2],
+                stats.parts_y[3]
+            )
+            .unwrap();
+        }
+        let at_end = step + 1 == coupled_steps;
+        if let Some(dir) = save_dir
+            .as_ref()
+            .filter(|_| at_end || (step + 1) % save_every == 0)
+        {
+            let f = fl.borrow();
+            let g = f.grid;
+            state::Saved {
+                t: t_new,
+                dims: [g.nx, g.ny, g.nz],
+                h: f.h,
+                u: f.field.u.clone(),
+                v: f.field.v.clone(),
+                w: f.field.w.clone(),
+                p: f.field.p.clone(),
+                disp: flag_state.displacement.as_slice().to_vec(),
+                vel: flag_state.velocity.as_slice().to_vec(),
+                acc: flag_state.acceleration.as_slice().to_vec(),
+                line: line_n.clone(),
+                c_fluid: c_fluid_n.clone(),
+                nodal: committed_nodal
+                    .iter()
+                    .flat_map(|(_, v)| [v.x, v.y])
+                    .collect(),
+            }
+            .save(dir)
+            .expect("save the coupled state");
+        }
+        if (step + 1) % 250 == 0 {
+            let w = &uy_series[uy_series.len().saturating_sub(600)..];
+            let (mid, amp) = mid_amp(w);
+            println!(
+                "  t {t_new:.3} ({} steps): uy(A) {uy:+.4e} ux {ux:+.4e} (last ~1 period mid {mid:+.3e} amp {amp:.3e}), drag {:.1} lift {:+.1}; {:.2} subit/step (max {max_seen}, stalled {stalled}); fluid {:.0} s restore {:.0} s loads {:.0} s structure {:.0} s of {:.0} s",
+                step + 1,
+                tot[0],
+                tot[1],
+                total_subit as f64 / (step + 1) as f64,
+                t_fluid.get(),
+                t_restore.get(),
+                t_loads.get(),
+                t_struct.get(),
+                phase_start.elapsed().as_secs_f64()
+            );
+        }
+    }
+    // Summary: the last two seconds (or what there is).
+    let n = times.len();
+    if n > 0 {
+        let t_last = times[n - 1];
+        let from = times.partition_point(|&t| t < t_last - 2.0);
+        let (mid, amp) = mid_amp(&uy_series[from..]);
+        let (uxm, uxa) = mid_amp(&ux_series[from..]);
+        let f = fsi2_harness::crossing_frequency(&times[from..], &uy_series[from..]);
+        let mut d: Vec<f64> = drag_s[from..].to_vec();
+        let (dm, _) = mid_amp(&d);
+        let (lm, la) = mid_amp(&lift_s[from..]);
+        let dmed = median(&mut d);
+        println!(
+            "  FINAL R8-a ny {ny} nz {nz}: {n} coupled steps to t {t_last:.3}; last 2 s: uy(A) {:.2} ± {:.2} mm, ux(A) {:.2} ± {:.2} mm, f {}, drag mid {dm:.1} (median {dmed:.1}), lift {lm:+.1} ± {la:.1}; {:.2} subit/step (max {max_seen}, stalled {stalled}); death {}; wall {:.0} s",
+            1e3 * mid,
+            1e3 * amp,
+            1e3 * uxm,
+            1e3 * uxa,
+            f.map_or("n/a".into(), |f| format!("{f:.4} Hz")),
+            total_subit as f64 / n as f64,
+            death.as_deref().unwrap_or("none"),
+            start.elapsed().as_secs_f64()
+        );
+    }
+}
+
+mod fluid_step {
+    pub struct Stats {
+        pub residual: f64,
+        pub cg: usize,
+        pub fresh: usize,
+        pub sink: [f64; 2],
+        pub parts_y: [f64; 4],
+    }
+}
@@ -0,0 +1,353 @@
+//! R8-a: the embedded3 3D cut-cell fluid as the fluid side of a
+//! partitioned FSI loop. The body is the embedded flag test's (the circle
+//! wall to wall united with a capsule of half-thickness 10 mm around the
+//! flag's centreline, its apex on the benchmark's tip A), but the
+//! centreline is no longer prescribed: the harness sets it per coupled
+//! step from the 2D structure (span-uniform; the centreline's element-corner
+//! nodes at reference x 0.25 … 0.59, so the capsule's apex sits on A as the
+//! flag test's tip inset puts it), as the pair of lines at the
+//! step's start and end; the body's φ and surface velocity at any time in
+//! between are the linear blend of the two (the solver asks at the step's
+//! two ends only).
+//!
+//! The loads: the operator route (`Mask::cut_wall_force`) with the R8-a
+//! load sink installed — every contribution the route sums, with its
+//! position, is returned for the harness to distribute onto the flag.
+
+use std::cell::RefCell;
+use std::sync::atomic::{AtomicU64, Ordering};
+use std::sync::{Arc, Mutex, RwLock};
+
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::exchange::set_load_sink;
+use rtx_cfd::solvers::incompressible::embedded3::step::device::{DeviceSnapshot, DeviceStep};
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, DeviceSdf, Field, Fluid, Grid, Parameters, Side, Solver, StepResult,
+    WallScheme,
+};
+
+pub const H: f64 = 0.41;
+pub const L: f64 = 2.5;
+pub const CX: f64 = 0.2;
+pub const CY: f64 = 0.2;
+pub const R_CYL: f64 = 0.05;
+/// The capsule's half-thickness (the flag's half-thickness).
+pub const HALF: f64 = 0.01;
+pub const RHO: f64 = 1000.0;
+pub const NU: f64 = 1e-3;
+/// The flag test's CFL velocity (its dt convention, kept for comparability).
+const U_CFL: f64 = 2.25;
+
+/// A centreline at time `t`: points (x, y) and their velocities.
+#[derive(Clone, Debug)]
+pub struct Line {
+    pub t: f64,
+    pub pts: Vec<[f64; 2]>,
+    pub vel: Vec<[f64; 2]>,
+}
+
+/// The step's two lines.
+pub struct Lines {
+    pub a: Line,
+    pub b: Line,
+}
+
+impl Lines {
+    /// The blend at `t` (the start line before `a.t`, the end line after `b.t`).
+    fn at(&self, t: f64) -> Line {
+        let eps = 1e-12 * (1.0 + t.abs());
+        if (t - self.b.t).abs() <= eps || t >= self.b.t || self.b.t <= self.a.t {
+            return self.b.clone();
+        }
+        if t <= self.a.t + eps {
+            return self.a.clone();
+        }
+        let s = (t - self.a.t) / (self.b.t - self.a.t);
+        let mix = |p: &[[f64; 2]], q: &[[f64; 2]]| -> Vec<[f64; 2]> {
+            p.iter()
+                .zip(q)
+                .map(|(p, q)| [p[0] + s * (q[0] - p[0]), p[1] + s * (q[1] - p[1])])
+                .collect()
+        };
+        assert_eq!(
+            self.a.pts.len(),
+            self.b.pts.len(),
+            "blend needs equal point counts"
+        );
+        Line {
+            t,
+            pts: mix(&self.a.pts, &self.b.pts),
+            vel: mix(&self.a.vel, &self.b.vel),
+        }
+    }
+}
+
+/// Bumped on every `set_lines` (the host closures' per-thread cache key).
+static VERSION: AtomicU64 = AtomicU64::new(0);
+
+/// The capsule's signed distance at (x, y) and the centreline velocity at
+/// the closest point (the flag test's `flag_2d_recorded`, cached per thread
+/// and time).
+fn capsule(lines: &RwLock<Lines>, x: f64, y: f64, t: f64) -> (f64, (f64, f64)) {
+    thread_local! {
+        static CACHE: RefCell<(u64, u64, Vec<[f64; 4]>)> = const { RefCell::new((u64::MAX, u64::MAX, Vec::new())) };
+    }
+    let ver = VERSION.load(Ordering::Acquire);
+    CACHE.with(|cell| {
+        let mut c = cell.borrow_mut();
+        if c.0 != t.to_bits() || c.1 != ver {
+            let line = lines.read().expect("lines").at(t);
+            c.2 = line
+                .pts
+                .iter()
+                .zip(&line.vel)
+                .map(|(p, v)| [p[0], p[1], v[0], v[1]])
+                .collect();
+            c.0 = t.to_bits();
+            c.1 = ver;
+        }
+        let pts = &c.2;
+        let mut best = f64::INFINITY;
+        let mut v_best = (0.0, 0.0);
+        for m in 0..pts.len() - 1 {
+            let [ax, ay, avx, avy] = pts[m];
+            let [bx, by, bvx, bvy] = pts[m + 1];
+            let (ex, ey) = (bx - ax, by - ay);
+            let l2 = ex * ex + ey * ey;
+            if l2 == 0.0 {
+                continue;
+            }
+            let u = (((x - ax) * ex + (y - ay) * ey) / l2).clamp(0.0, 1.0);
+            let (px, py) = (ax + u * ex, ay + u * ey);
+            let d = ((x - px).powi(2) + (y - py).powi(2)).sqrt();
+            if d < best {
+                best = d;
+                v_best = (avx + u * (bvx - avx), avy + u * (bvy - avy));
+            }
+        }
+        (best - HALF, v_best)
+    })
+}
+
+fn cylinder(x: f64, y: f64) -> f64 {
+    ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL
+}
+
+/// One load contribution: position, component, part (0 pressure, 1 wall
+/// shear, 2 diffusive exchange, 3 convective exchange), force on the body
+/// (N, over the whole z extent).
+pub type Contribution = ([f64; 3], usize, usize, f64);
+
+pub struct E3Fluid {
+    pub device: DeviceStep,
+    pub field: Field,
+    pub grid: Grid,
+    pub h: f64,
+    pub dt: f64,
+    /// The duct's z extent.
+    pub width: f64,
+    /// The z extent the last `loads` integrated over (per unit span).
+    pub load_width: f64,
+    pub lines: Arc<RwLock<Lines>>,
+    sink: Arc<Mutex<Vec<Contribution>>>,
+}
+
+impl E3Fluid {
+    /// The fluid on the benchmark channel at rung `ny`: `nz_slab > 0` the
+    /// thin slab periodic in z (the flag as a 2D problem), 0 the full 0.41 m
+    /// duct with slip side walls; the 2D inflow (parabolic in y, Ū 1) in
+    /// both. `speed` bounds the flag's surface speed (the narrow band).
+    /// `start`: a saved state (its time, its line, its fields — extruded
+    /// onto every plane when the saved nz differs); the rest flow otherwise.
+    pub fn build(
+        ny: usize,
+        nz_slab: usize,
+        speed: f64,
+        rest: Line,
+        start: Option<&super::state::Saved>,
+    ) -> Self {
+        let rest = match start {
+            Some(s) => s.line.clone(),
+            None => rest,
+        };
+        let h = H / ny as f64;
+        let nx = (L / h).round() as usize;
+        let nz = if nz_slab > 0 {
+            nz_slab
+        } else {
+            (H / h).round() as usize
+        };
+        // The flag test's step (its CFL velocity; `speed` is the band's bound only).
+        let dt = (0.3 * h / U_CFL).min(0.5 * h * h / (6.0 * NU))
+            * super::env_f("RTX_E3FSI_DT_SCALE", 1.0);
+        let boundaries = if nz_slab > 0 {
+            Boundaries {
+                x1: Side::PressureOutlet,
+                z0: Side::Periodic,
+                z1: Side::Periodic,
+                ..Boundaries::default()
+            }
+        } else {
+            Boundaries {
+                x1: Side::PressureOutlet,
+                z0: Side::SlipWall,
+                z1: Side::SlipWall,
+                ..Boundaries::default()
+            }
+        };
+        let mut solver = Solver::new(
+            Fluid {
+                density: RHO,
+                viscosity: RHO * NU,
+                reference_velocity: 1.0,
+                reference_length: 2.0 * R_CYL,
+            },
+            Parameters {
+                corrector_steps: super::env_f("RTX_E3FSI_CORRECTORS", 3.0) as usize,
+                inner_stop_factor: super::env_f("RTX_E3FSI_INNER", 1e-3),
+                tolerance: 1e-8,
+                convection_scheme: ConvectionScheme::TvdVanAlbada,
+                wall_scheme: WallScheme::CutCell,
+                boundaries,
+                max_surface_speed: Some(speed),
+                ..Parameters::default()
+            },
+        );
+        let inflow = |y: f64| 6.0 * y * (H - y) / (H * H);
+        solver.set_boundary_velocity(move |x, y, _z, _t| {
+            if x <= 0.0 {
+                (inflow(y), 0.0, 0.0)
+            } else {
+                (0.0, 0.0, 0.0)
+            }
+        });
+        let lines = Arc::new(RwLock::new(Lines {
+            a: rest.clone(),
+            b: rest,
+        }));
+        VERSION.fetch_add(1, Ordering::AcqRel);
+        let (l1, l2, l3) = (lines.clone(), lines.clone(), lines.clone());
+        let body = Body::from_sdf(move |x, y, _z, t| cylinder(x, y).min(capsule(&l1, x, y, t).0))
+            .with_surface_velocity(move |x, y, _z, t| {
+                let (df, (vx, vy)) = capsule(&l2, x, y, t);
+                if df <= cylinder(x, y) {
+                    (vx, vy, 0.0)
+                } else {
+                    (0.0, 0.0, 0.0)
+                }
+            });
+        let width = nz as f64 * h;
+        let body = body.with_device_sdf(move |t| {
+            let line = l3.read().expect("lines").at(t);
+            DeviceSdf {
+                cyl: [CX, CY, R_CYL],
+                cyl_cut: false,
+                flag_cut: false,
+                zc: 0.5 * width,
+                span: width,
+                r_edge: h,
+                half: HALF,
+                fillet: 0.0,
+                poly: line.pts,
+                vel: line.vel,
+                // R8-c's plate body (merged alongside): the span-uniform harness keeps the polyline.
+                plate: None,
+            }
+        });
+        solver.set_moving_body(body);
+        let g = Grid::cubic(nx, ny, nz, h);
+        let mut field = Field::new(g);
+        for k in 0..nz {
+            for j in 0..ny {
+                let u0 = inflow((j as f64 + 0.5) * h);
+                for i in 0..=nx {
+                    field.u[g.uface(k, j, i)] = u0;
+                }
+            }
+        }
+        if let Some(s) = start {
+            s.fill(&mut field);
+            solver.set_time(s.t);
+        }
+        solver.initialize(&mut field);
+        let mut device = DeviceStep::new(solver, g);
+        device.upload(&field);
+        println!(
+            "  R8-a fluid: ny {ny}, {nx}×{ny}×{nz} = {} cells ({}), h {h:.4e}, dt {dt:.4e}, speed bound {speed} m/s",
+            g.cells(),
+            if nz_slab > 0 {
+                "slab, periodic z"
+            } else {
+                "full duct, slip sides"
+            }
+        );
+        Self {
+            device,
+            field,
+            grid: g,
+            h,
+            dt,
+            width,
+            load_width: width,
+            lines,
+            sink: Arc::new(Mutex::new(Vec::new())),
+        }
+    }
+
+    pub fn time(&self) -> f64 {
+        self.device.solver.time()
+    }
+
+    /// The step's start and end lines (their times are the lines' own).
+    pub fn set_lines(&self, a: Line, b: Line) {
+        *self.lines.write().expect("lines") = Lines { a, b };
+        VERSION.fetch_add(1, Ordering::AcqRel);
+    }
+
+    pub fn step(&mut self) -> StepResult {
+        self.device.advance(self.dt)
+    }
+
+    pub fn snapshot(&mut self) -> DeviceSnapshot {
+        self.device.snapshot()
+    }
+
+    pub fn restore(&mut self, snap: &DeviceSnapshot) {
+        self.device.restore(snap);
+    }
+
+    /// The operator-route force on the whole body per unit span and the
+    /// route's contributions (N over the whole z extent).
+    pub fn loads(&mut self) -> ([f64; 3], Vec<Contribution>) {
+        self.device.download(&mut self.field);
+        let t = self.time();
+        self.sink.lock().expect("sink").clear();
+        let s = self.sink.clone();
+        let previous = set_load_sink(Some(Box::new(move |pos, c, part, v| {
+            s.lock().expect("sink").push((pos, c, part, v));
+        })));
+        assert!(previous.is_none(), "a load sink was already installed");
+        let mask = self.device.solver.mask().expect("mask");
+        let body = self.device.solver.body().expect("body");
+        // `RTX_E3FSI_LOAD_PLANES=n`: the route on the n mid planes only (per
+        // unit span of those planes) — a cost knob for a spanwise-uniform
+        // flow; the whole span by default.
+        let nz = self.grid.nz;
+        let n = (super::env_f("RTX_E3FSI_LOAD_PLANES", 0.0) as usize).min(nz);
+        let f = if n > 0 && n < nz {
+            let k0 = (nz - n) / 2;
+            self.load_width = n as f64 * self.h;
+            mask.cut_wall_force_per_span(body, &self.field, RHO * NU, t, (k0, k0 + n))
+                .expect("cut wall force per span")
+        } else {
+            self.load_width = self.width;
+            let f = mask
+                .cut_wall_force(body, &self.field, RHO * NU, t)
+                .expect("cut wall force");
+            [f[0] / self.width, f[1] / self.width, f[2] / self.width]
+        };
+        set_load_sink(None);
+        let contributions = std::mem::take(&mut *self.sink.lock().expect("sink"));
+        (f, contributions)
+    }
+}
@@ -0,0 +1,162 @@
+//! R8-a: the coupled state on disk — the fluid's fields, the flag's
+//! kinematic state, the fluid's last centreline and the committed load —
+//! so a march can continue on the same grid or be EXTRUDED onto the full
+//! duct (the slab's z-average onto every plane: the 3D solver started on
+//! the 2D problem's own state, rule 16).
+
+use std::io::{Read as _, Write as _};
+use std::path::Path;
+
+use super::fluid::Line;
+use rtx_cfd::solvers::incompressible::embedded3::{Field, Grid};
+
+pub struct Saved {
+    pub t: f64,
+    pub dims: [usize; 3],
+    pub h: f64,
+    pub u: Vec<f64>,
+    pub v: Vec<f64>,
+    pub w: Vec<f64>,
+    pub p: Vec<f64>,
+    pub disp: Vec<f64>,
+    pub vel: Vec<f64>,
+    pub acc: Vec<f64>,
+    pub line: Line,
+    pub c_fluid: Vec<f64>,
+    /// The committed nodal load (fx, fy per wetted node, in wetted order).
+    pub nodal: Vec<f64>,
+}
+
+fn write_vec(dir: &Path, name: &str, v: &[f64]) -> std::io::Result<()> {
+    let mut f = std::fs::File::create(dir.join(format!("{name}.f64")))?;
+    let mut bytes = Vec::with_capacity(8 * v.len());
+    for x in v {
+        bytes.extend_from_slice(&x.to_le_bytes());
+    }
+    f.write_all(&bytes)
+}
+
+fn read_vec(dir: &Path, name: &str) -> std::io::Result<Vec<f64>> {
+    let mut bytes = Vec::new();
+    std::fs::File::open(dir.join(format!("{name}.f64")))?.read_to_end(&mut bytes)?;
+    Ok(bytes
+        .chunks_exact(8)
+        .map(|c| f64::from_le_bytes(c.try_into().unwrap()))
+        .collect())
+}
+
+impl Saved {
+    pub fn save(&self, dir: &str) -> std::io::Result<()> {
+        let d = Path::new(dir);
+        std::fs::create_dir_all(d)?;
+        let flat = |pts: &[[f64; 2]]| pts.iter().flat_map(|p| [p[0], p[1]]).collect::<Vec<_>>();
+        write_vec(d, "u", &self.u)?;
+        write_vec(d, "v", &self.v)?;
+        write_vec(d, "w", &self.w)?;
+        write_vec(d, "p", &self.p)?;
+        write_vec(d, "disp", &self.disp)?;
+        write_vec(d, "vel", &self.vel)?;
+        write_vec(d, "acc", &self.acc)?;
+        write_vec(d, "line_pts", &flat(&self.line.pts))?;
+        write_vec(d, "line_vel", &flat(&self.line.vel))?;
+        write_vec(d, "c_fluid", &self.c_fluid)?;
+        write_vec(d, "nodal", &self.nodal)?;
+        std::fs::write(
+            d.join("meta.txt"),
+            format!(
+                "{:e} {} {} {} {:e} {:e}\n",
+                self.t, self.dims[0], self.dims[1], self.dims[2], self.h, self.line.t
+            ),
+        )
+    }
+
+    pub fn load(dir: &str) -> std::io::Result<Self> {
+        let d = Path::new(dir);
+        let meta = std::fs::read_to_string(d.join("meta.txt"))?;
+        let m: Vec<&str> = meta.split_whitespace().collect();
+        let pairs = |v: Vec<f64>| v.chunks_exact(2).map(|c| [c[0], c[1]]).collect::<Vec<_>>();
+        Ok(Self {
+            t: m[0].parse().unwrap(),
+            dims: [
+                m[1].parse().unwrap(),
+                m[2].parse().unwrap(),
+                m[3].parse().unwrap(),
+            ],
+            h: m[4].parse().unwrap(),
+            u: read_vec(d, "u")?,
+            v: read_vec(d, "v")?,
+            w: read_vec(d, "w")?,
+            p: read_vec(d, "p")?,
+            disp: read_vec(d, "disp")?,
+            vel: read_vec(d, "vel")?,
+            acc: read_vec(d, "acc")?,
+            line: Line {
+                t: m[5].parse().unwrap(),
+                pts: pairs(read_vec(d, "line_pts")?),
+                vel: pairs(read_vec(d, "line_vel")?),
+            },
+            c_fluid: read_vec(d, "c_fluid")?,
+            nodal: read_vec(d, "nodal")?,
+        })
+    }
+
+    /// The saved fields onto `field` (same nx, ny): the saved planes'
+    /// z-average on every plane of the target (w = 0: the 2D problem's
+    /// state); identical planes copy through when nz matches.
+    pub fn fill(&self, field: &mut Field) {
+        let g: Grid = field.grid;
+        let [nx, ny, nzs] = self.dims;
+        assert_eq!(
+            (g.nx, g.ny),
+            (nx, ny),
+            "the saved state's grid differs in x or y"
+        );
+        assert!(
+            (g.dx - self.h).abs() < 1e-12 * self.h,
+            "the saved state's h differs"
+        );
+        let same = g.nz == nzs;
+        let src = Grid::cubic(nx, ny, nzs, self.h);
+        for j in 0..ny {
+            for i in 0..=nx {
+                let mean = (0..nzs).map(|k| self.u[src.uface(k, j, i)]).sum::<f64>() / nzs as f64;
+                for k in 0..g.nz {
+                    field.u[g.uface(k, j, i)] = if same {
+                        self.u[src.uface(k, j, i)]
+                    } else {
+                        mean
+                    };
+                }
+            }
+        }
+        for j in 0..=ny {
+            for i in 0..nx {
+                let mean = (0..nzs).map(|k| self.v[src.vface(k, j, i)]).sum::<f64>() / nzs as f64;
+                for k in 0..g.nz {
+                    field.v[g.vface(k, j, i)] = if same {
+                        self.v[src.vface(k, j, i)]
+                    } else {
+                        mean
+                    };
+                }
+            }
+        }
+        if same {
+            field.w.copy_from_slice(&self.w);
+        } else {
+            field.w.iter_mut().for_each(|w| *w = 0.0);
+        }
+        for j in 0..ny {
+            for i in 0..nx {
+                let mean = (0..nzs).map(|k| self.p[src.cell(k, j, i)]).sum::<f64>() / nzs as f64;
+                for k in 0..g.nz {
+                    field.p[g.cell(k, j, i)] = if same {
+                        self.p[src.cell(k, j, i)]
+                    } else {
+                        mean
+                    };
+                }
+            }
+        }
+    }
+}