R8-a harness: save/load of the coupled state (continue on the same grid, or extrude the slab onto the full duct)

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

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

@@ -28,13 +28,15 @@
 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::Vector3;
+use nalgebra::{DVector, Vector3};
 use rtx_fea::analysis::{
     AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
 };
@@ -251,7 +253,26 @@ fn fsi2_on_embedded3() {
     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);
-    let mut fl = E3Fluid::build(ny, nz, speed, rest.clone());
+    // `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)"
@@ -268,7 +289,11 @@ fn fsi2_on_embedded3() {
     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 = (t_rigid / dt).round() as usize;
+    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();
@@ -318,7 +343,7 @@ fn fsi2_on_embedded3() {
         tot0[1],
         start.elapsed().as_secs_f64()
     );
-    if rigid_only {
+    if rigid_only && saved.is_none() {
         return;
     }
 
@@ -363,16 +388,34 @@ fn fsi2_on_embedded3() {
         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.width);
             flag.borrow_mut().set_nodal_forces(&nodal0);
-    let mut flag_state = flag.borrow_mut().rest_state().unwrap();
-    let mut committed_nodal = nodal0;
+            let state = flag.borrow_mut().rest_state().unwrap();
             // The fluid's own previous line and centreline (its geometry's history).
-    let mut line_n = Line {
+            let line = Line {
                 t: fl.time(),
                 ..rest.clone()
             };
-    let mut c_fluid_n = zero_c.clone();
+            (state, nodal0, line, zero_c.clone())
+        }
+    };
 
     let coupled_steps = ((t_end - fl.time()) / dt).round() as usize;
     let fl = RefCell::new(fl);
@@ -537,6 +580,34 @@ fn fsi2_on_embedded3() {
             )
             .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);

crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs

@@ -155,7 +155,19 @@ impl E3Fluid {
     /// 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).
-    pub fn build(ny: usize, nz_slab: usize, speed: f64, rest: Line) -> Self {
+    /// `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 {
@@ -249,6 +261,10 @@ impl E3Fluid {
                 }
             }
         }
+        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);

crates/specialized/rtx-fsi/tests/fsi2_embedded3/state.rs

@@ -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
+                    };
+                }
+            }
+        }
+    }
+}