diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs index ab517a6..5afd88b 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs +++ b/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=`: continue from a saved coupled state (extruded + // onto the full duct when the saved nz differs); `RTX_E3FSI_SAVE=` + // 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 (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; - // The fluid's own previous line and centreline (its geometry's history). - let mut line_n = Line { - t: fl.time(), - ..rest.clone() + let (mut flag_state, mut committed_nodal, mut line_n, mut c_fluid_n) = match &saved { + Some(s) => { + let nodal: Vec<(NodeId, Vector3)> = 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 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 mut c_fluid_n = 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); diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs index 26efa05..dee7d79 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs +++ b/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); diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/state.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/state.rs new file mode 100644 index 0000000..a15307f --- /dev/null +++ b/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, + pub v: Vec, + pub w: Vec, + pub p: Vec, + pub disp: Vec, + pub vel: Vec, + pub acc: Vec, + pub line: Line, + pub c_fluid: Vec, + /// The committed nodal load (fx, fy per wetted node, in wetted order). + pub nodal: Vec, +} + +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> { + 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::>(); + 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 { + 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| v.chunks_exact(2).map(|c| [c[0], c[1]]).collect::>(); + 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::() / 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::() / 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::() / 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 + }; + } + } + } + } +}