diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs index 36bedb3..a760f84 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -942,6 +942,16 @@ impl Mask { for c in 0..3 { v[c] = f.p[idx] * w[c]; pressure[c] += v[c]; + super::exchange::to_load_sink( + [ + (i as f64 + 0.5) * g.dx, + (j as f64 + 0.5) * g.dy, + (k as f64 + 0.5) * g.dz, + ], + c, + 0, + v[c], + ); } let x = [ (i as f64 + 0.5) * g.dx, @@ -1000,6 +1010,7 @@ impl Mask { let mut fv = [0.0; 3]; fv[c] = v; sink(LoadKind::Shear, x, fv); + super::exchange::to_load_sink(x, c, 1, v); } } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs index c00f147..d6c543d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs @@ -30,6 +30,36 @@ pub(super) fn in_load_window(x: f64) -> bool { .is_none_or(|(x0, x1)| x >= x0 && x < x1) } +/// R8-a: a sink for the operator load route's CONTRIBUTIONS — each term +/// the route sums (part 0 the cell's `p W`, 1 the face's wall shear, 2 and +/// 3 the face's diffusive and convective exchange) is handed to it as +/// (position, component, part, force on the body) while the totals are +/// summed as before. +/// A coupled harness distributes them onto its structure. Process-wide; +/// `None` (the default) hands nothing and the routes' sums are unchanged. +/// The gradient-weight term (the host prototype `pressure_centroid`) is +/// not handed. +pub type LoadSink = Box; + +static LOAD_SINK: std::sync::Mutex> = std::sync::Mutex::new(None); +static LOAD_SINK_ON: std::sync::atomic::AtomicBool = std::sync::atomic::AtomicBool::new(false); + +/// Install (or clear, with `None`) the load sink; returns the previous one. +pub fn set_load_sink(sink: Option) -> Option { + let mut guard = LOAD_SINK.lock().expect("load sink"); + LOAD_SINK_ON.store(sink.is_some(), std::sync::atomic::Ordering::SeqCst); + std::mem::replace(&mut *guard, sink) +} + +#[inline] +pub(super) fn to_load_sink(pos: [f64; 3], c: usize, part: usize, value: f64) { + if LOAD_SINK_ON.load(std::sync::atomic::Ordering::Relaxed) { + if let Some(f) = LOAD_SINK.lock().expect("load sink").as_mut() { + f(pos, c, part, value); + } + } +} + impl Mask { /// The cut-cell load route: the force on the body from the operators /// themselves — `Σ_c p_c W_c` over the cells plus the implicit wall @@ -223,10 +253,12 @@ impl Mask { let v = -rho * m_plus * (u_face - u0); convective[c] -= v; sink(LoadKind::ExchangeConvective, xf, comp(c, -v)); + to_load_sink(xf, c, 3, -v); } let v = mu * cv.ap[d][1] * a_d * (un - u0) / solid_spacing(1.0); force[c] -= v; sink(LoadKind::ExchangeDiffusive, xf, comp(c, -v)); + to_load_sink(xf, c, 2, -v); } } // Minus side. @@ -245,11 +277,13 @@ impl Mask { let v = rho * m_minus * (u_face - u0); convective[c] -= v; sink(LoadKind::ExchangeConvective, xf, comp(c, -v)); + to_load_sink(xf, c, 3, -v); } let v = mu * cv.ap[d][0] * a_d * (ud - u0) / solid_spacing(-1.0); force[c] -= v; sink(LoadKind::ExchangeDiffusive, xf, comp(c, -v)); + to_load_sink(xf, c, 2, -v); } } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs index 2b62e30..999412c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs @@ -9,6 +9,9 @@ mod cut; mod geom; mod mask; mod poisson_setup; +mod snapshot; + +pub use snapshot::DeviceSnapshot; use super::{Side, Solver, StepResult}; use crate::solvers::incompressible::embedded3::field::Field; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/snapshot.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/snapshot.rs new file mode 100644 index 0000000..7137a2e --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/snapshot.rs @@ -0,0 +1,124 @@ +//! R8-a: the step's snapshot and restore for a partitioned FSI loop — the +//! coupled harness re-runs the same fluid step once per subiteration with +//! a new candidate body position, so the device fields and the moving +//! body's host state at the start of the step must come back exactly. +//! +//! Nothing here runs unless a caller asks for it: every existing path is +//! untouched (byte-identical by construction). +//! +//! What a restore rebuilds instead of copying (the persistent R6 state +//! whose band history would otherwise describe the rejected pass): the +//! device geometry and classification (`geom`, `dmask`; the next step +//! re-syncs them from the restored mask, as the first moving step of a run +//! does), the predictor tables (`DeviceCut::build` at the restored time, +//! the identity reference of `predictor_from`), the multigrid hierarchy +//! and the projected-guess basis. The restored step is the same step up to +//! the CG's initial guess (the Poisson solves' own tolerance). + +use super::DeviceStep; +use super::cut::{DeviceCut, Phase}; +use crate::solvers::incompressible::embedded3::poisson::device::runtime; +use crate::solvers::incompressible::embedded3::wall::Mask; +use cudarc::driver::CudaSlice; + +/// The start-of-step state of a [`DeviceStep`] (device copies of the +/// fields plus the moving body's host records). +pub struct DeviceSnapshot { + time: f64, + fields: Vec>, + mask: Option, + vol_old: Vec, + apertures_old: Option<[Vec; 3]>, + wall_fluxes: Vec, + last_ghost_correction: f64, + mask_gen: u64, +} + +impl DeviceSnapshot { + /// The solver time the snapshot was taken at. + #[must_use] + pub fn time(&self) -> f64 { + self.time + } +} + +impl DeviceStep { + fn field_slots(&mut self) -> [&mut CudaSlice; 8] { + [ + &mut self.u, + &mut self.v, + &mut self.w, + &mut self.p, + &mut self.p_prime, + &mut self.u_old, + &mut self.v_old, + &mut self.w_old, + ] + } + + /// R8-a: the start-of-step state (device fields copied on the device, + /// the host mask cloned). Take it between steps. + pub fn snapshot(&mut self) -> DeviceSnapshot { + let rt = runtime(); + let mut fields = Vec::with_capacity(8); + for src in self.field_slots() { + let mut dst = rt.stream.alloc_zeros::(src.len()).expect("alloc"); + rt.stream.memcpy_dtod(&*src, &mut dst).expect("snapshot"); + fields.push(dst); + } + rt.stream.synchronize().expect("sync"); + let s = &self.solver; + let mask = s.mask.clone().map(|mut m| { + // A clone is a host generation: its arrays must not return to + // the device geometry's recycling pool under a live generation. + if let Some(c) = m.cut.as_mut() { + c.generation = 0; + } + m + }); + DeviceSnapshot { + time: s.time, + fields, + mask, + vol_old: s.vol_old.clone(), + apertures_old: s.apertures_old.clone(), + wall_fluxes: s.wall_fluxes.clone(), + last_ghost_correction: s.last_ghost_correction, + mask_gen: s.mask_gen, + } + } + + /// R8-a: back to `snap` (taken on this stepper). The body's functions + /// must answer for the snapshot's time as they did when it was taken + /// (the predictor tables are rebuilt from them). + pub fn restore(&mut self, snap: &DeviceSnapshot) { + let rt = runtime(); + for (dst, src) in self.field_slots().into_iter().zip(&snap.fields) { + rt.stream.memcpy_dtod(src, dst).expect("restore"); + } + rt.stream.synchronize().expect("sync"); + { + let s = &mut self.solver; + s.time = snap.time; + s.mask = snap.mask.clone(); + s.vol_old = snap.vol_old.clone(); + s.apertures_old = snap.apertures_old.clone(); + s.apertures_old_gen = 0; + s.wall_fluxes = snap.wall_fluxes.clone(); + s.last_ghost_correction = snap.last_ghost_correction; + s.mask_gen = snap.mask_gen; + s.pending_cut = None; + s.pending_mask = None; + s.mask_pool = None; + s.geom_pool = Default::default(); + } + self.geom = None; + self.dmask = None; + self.cg = None; + self.steps_since_hierarchy = 0; + self.guess.clear(); + if self.cut.is_some() { + self.cut = DeviceCut::build(&self.solver, self.grid, Phase::Predictor, snap.time); + } + } +} diff --git a/crates/specialized/rtx-fsi/Cargo.toml b/crates/specialized/rtx-fsi/Cargo.toml index 1313d60..7ffe3ab 100644 --- a/crates/specialized/rtx-fsi/Cargo.toml +++ b/crates/specialized/rtx-fsi/Cargo.toml @@ -21,5 +21,9 @@ futures = { workspace = true } rtx-cfd = { workspace = true } rtx-fea = { workspace = true } +[features] +# R8-a: the coupled FSI on the embedded3 device fluid (`tests/fsi2_embedded3.rs`). +cuda = ["rtx-cfd/cuda"] + [lints] workspace = true diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs new file mode 100644 index 0000000..99cf106 --- /dev/null +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs @@ -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= 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, + interface: Interface, + /// (reference x, wetted index or None for the clamp corner), ascending. + bottom: Vec<(f64, Option)>, + top: Vec<(f64, Option)>, + /// (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)> { + let mut v: Vec<(f64, Option)> = 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(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; 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=`: 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)" + ); + 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 { + 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)> = 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 = 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)>, + Line, + Vec, + [f64; 3], + [f64; 2], + fluid_step::Stats, + ); + let latest: RefCell> = RefCell::new(None); + let pass = |cand: &[f64]| -> Vec { + 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::() + .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::() + .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::() + .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(×[from..], &uy_series[from..]); + let mut d: Vec = 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], + } +} diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs new file mode 100644 index 0000000..4ba4ef8 --- /dev/null +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3/fluid.rs @@ -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, 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>, + sink: Arc>>, +} + +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) { + 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) + } +} 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 + }; + } + } + } + } +}