From 5b1621e6ada81e047c838dffda1d7c9414053b4f Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Thu, 17 Sep 2026 16:20:35 -0500 Subject: [PATCH] =?UTF-8?q?embedded3=20item=2011:=20moving=20bodies=20(end?= =?UTF-8?q?-of-step=20mask,=20fresh-cell=20refill,=20space-time=20cut=20ce?= =?UTF-8?q?ll:=20step-averaged=20apertures,=20GCL=20wall=20flux,=20Reynold?= =?UTF-8?q?s-transport=20momentum),=20the=203D=20fresh-cell=20falsifier=20?= =?UTF-8?q?(plate=20/=20circle=20/=20stadium,=20wall=20+=20control-volume?= =?UTF-8?q?=20routes)=20and=20the=20Lipschitz=20sweep;=20ghost=20wall=20re?= =?UTF-8?q?produces=20the=202D=20falsifier=20to=20the=20digit;=20cut=20wal?= =?UTF-8?q?l=205=E2=80=9314=C3=97=20smoother=20on=20the=20circle,=20gates?= =?UTF-8?q?=20not=20met=20(fresh=20cell's=20first=20step);=20wall.rs=20spl?= =?UTF-8?q?it=20(impose.rs)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../incompressible/embedded3/cutwall.rs | 87 ++- .../incompressible/embedded3/impose.rs | 139 +++++ .../solvers/incompressible/embedded3/mod.rs | 1 + .../embedded3/step/cut_predictor.rs | 13 +- .../incompressible/embedded3/step/device.rs | 1 + .../incompressible/embedded3/step/mod.rs | 177 +++++- .../embedded3/step/projection.rs | 28 +- .../solvers/incompressible/embedded3/wall.rs | 203 +++---- .../rtx-cfd/tests/embedded3_embedded_mms.rs | 255 +------- .../rtx-cfd/tests/embedded3_falsifier.rs | 566 ++++++++++++++++++ .../rtx-cfd/tests/embedded3_sphere/mod.rs | 255 ++++++++ .../rtx-cfd/tests/embedded3_wall_lipschitz.rs | 100 ++++ 12 files changed, 1406 insertions(+), 419 deletions(-) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_wall_lipschitz.rs 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 db420a7..4713db4 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -188,9 +188,34 @@ impl Mask { anchor, fluid_cells, cut: Some(cut), + step_apertures: None, + step_open: None, }) } + /// Set the step-averaged apertures and the space-time classification + /// from the previous mask's geometry. + pub fn set_step_apertures(&mut self, old: &Mask) { + let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { + return; + }; + let avg = |a: &[f64], b: &[f64]| -> Vec { + a.iter().zip(b).map(|(x, y)| 0.5 * (x + y)).collect() + }; + let au = avg(&cut.a_u, &old_cut.a_u); + let av = avg(&cut.a_v, &old_cut.a_v); + let aw = avg(&cut.a_w, &old_cut.a_w); + let open = |a: &[f64]| -> Vec { a.iter().map(|&x| x > 0.0).collect() }; + let active = self + .cell_fluid + .iter() + .zip(&old.cell_fluid) + .map(|(&n, &o)| n || o) + .collect(); + self.step_open = Some((open(&au), open(&av), open(&aw), active)); + self.step_apertures = Some((au, av, aw)); + } + pub(super) fn lattice(&self) -> Lattice { Lattice { g: self.grid, @@ -319,6 +344,51 @@ impl Mask { (table, correction) } + /// The moving rigid body's wall fluxes by the discrete geometric + /// conservation law: `(V_c^{n+1} − V_c^n)/dt` per active cell (a dying + /// cell's remaining volume leaves through its step-averaged apertures), + /// the net (the cut geometry's closure defect) redistributed over the + /// wall cells by wall area. + pub fn gcl_flux_table(&self, old: &Mask, dt: f64) -> (Vec, f64) { + let mut table = vec![0.0; self.grid.cells()]; + let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { + return (table, 0.0); + }; + let g = self.grid; + let dv = g.dx * g.dy * g.dz; + let (mut net, mut area) = (0.0, 0.0); + for (idx, entry) in table.iter_mut().enumerate() { + if !self.cell_active(idx) { + continue; + } + *entry = (cut.vol[idx] - old_cut.vol[idx]) * dv / dt; + net += *entry; + let w = cut.wall[idx]; + area += (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); + } + if std::env::var_os("RTX_E3_DEBUG").is_some() { + let dead = (0..table.len()) + .filter(|&i| !self.cell_fluid[i] && old.cell_fluid[i]) + .count(); + let fresh = (0..table.len()) + .filter(|&i| self.cell_fluid[i] && !old.cell_fluid[i]) + .count(); + let (vn, vn1): (f64, f64) = (old_cut.vol.iter().sum(), cut.vol.iter().sum()); + eprintln!( + " gcl: dead {dead} fresh {fresh} net {net:.3e} area {area:.3e} ΣV old {vn:.6} new {vn1:.6} (Δ {:.3e})", + vn1 - vn + ); + } + let correction = if area > 0.0 { net / area } else { 0.0 }; + if correction != 0.0 { + for (idx, w) in cut.wall.iter().enumerate() { + let a = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); + table[idx] -= correction * a; + } + } + (table, correction) + } + /// The volume flux of the surface velocity through a cell's wall into /// the body, `U_b·W_c`, uncorrected. Zero without a cut geometry. pub fn wall_flux(&self, body: &Body, idx: usize, t: f64) -> f64 { @@ -348,14 +418,27 @@ impl Mask { /// shear `Σ_f μ A_w (u_f − U_b)/d_f` over the unknown faces. `None` /// without a cut geometry. pub fn cut_wall_force(&self, body: &Body, f: &Field, mu: f64, t: f64) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_parts(body, f, mu, t)?; + Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]]) + } + + /// The cut-cell load route split into its pressure and shear parts. + pub fn cut_wall_force_parts( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + ) -> Option<([f64; 3], [f64; 3])> { let cut = self.cut.as_ref()?; let g = self.grid; let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let mut pressure = [0.0; 3]; let mut force = [0.0; 3]; for (idx, w) in cut.wall.iter().enumerate() { if self.cell_fluid[idx] { for c in 0..3 { - force[c] += f.p[idx] * w[c]; + pressure[c] += f.p[idx] * w[c]; } } } @@ -395,6 +478,6 @@ impl Mask { } } } - Some(force) + Some((pressure, force)) } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs new file mode 100644 index 0000000..6b0a176 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/impose.rs @@ -0,0 +1,139 @@ +//! The wall's imposition on the velocity field (`impl Mask` continued +//! from `wall.rs`, split for the file-size rule): prescribed faces take +//! the surface velocity, ghost faces their reconstruction from the source +//! field minus the shared flux compatibility correction. + +use super::body::Body; +use super::wall::{FaceKind, Mask}; + +impl Mask { + /// Impose the wall on `(u, v, w)` from the same field. + pub fn impose(&self, body: &Body, u: &mut [f64], v: &mut [f64], w: &mut [f64], t: f64) -> f64 { + let (us, vs, ws) = (u.to_vec(), v.to_vec(), w.to_vec()); + self.impose_from(body, &us, &vs, &ws, u, v, w, t) + } + + /// Solid faces: the surface velocity; ghost faces: the reconstruction + /// from the SOURCE field, minus the shared flux compatibility + /// correction over the flux-carrying ghosts. Returns the correction. + #[allow(clippy::too_many_arguments)] + pub fn impose_from( + &self, + body: &Body, + u_src: &[f64], + v_src: &[f64], + w_src: &[f64], + u: &mut [f64], + v: &mut [f64], + w: &mut [f64], + t: f64, + ) -> f64 { + let g = self.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + for k in 0..nz { + for j in 0..ny { + for i in 1..nx { + let idx = g.uface(k, j, i); + if self.u_kind[idx] == FaceKind::Solid { + u[idx] = body + .surface_velocity( + i as f64 * dx, + (j as f64 + 0.5) * dy, + (k as f64 + 0.5) * dz, + t, + ) + .0; + } + } + } + for j in 1..ny { + for i in 0..nx { + let idx = g.vface(k, j, i); + if self.v_kind[idx] == FaceKind::Solid { + v[idx] = body + .surface_velocity( + (i as f64 + 0.5) * dx, + j as f64 * dy, + (k as f64 + 0.5) * dz, + t, + ) + .1; + } + } + } + } + for k in 0..=nz { + for j in 0..ny { + for i in 0..nx { + let idx = g.wface(k, j, i); + if self.w_kind[idx] == FaceKind::Solid { + w[idx] = body + .surface_velocity( + (i as f64 + 0.5) * dx, + (j as f64 + 0.5) * dy, + k as f64 * dz, + t, + ) + .2; + } + } + } + } + let u_vals: Vec = self + .u_ghosts + .iter() + .map(|gh| gh.reconstruct(u_src)) + .collect(); + let v_vals: Vec = self + .v_ghosts + .iter() + .map(|gh| gh.reconstruct(v_src)) + .collect(); + let w_vals: Vec = self + .w_ghosts + .iter() + .map(|gh| gh.reconstruct(w_src)) + .collect(); + let (au, av, aw) = (dy * dz, dx * dz, dx * dy); + let mut net = 0.0; + let mut area = 0.0; + for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) { + if gh.flux_sign != 0.0 { + net += gh.flux_sign * val * au; + area += au; + } + } + for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) { + if gh.flux_sign != 0.0 { + net += gh.flux_sign * val * av; + area += av; + } + } + for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) { + if gh.flux_sign != 0.0 { + net += gh.flux_sign * val * aw; + area += aw; + } + } + let correction = if area > 0.0 { net / area } else { 0.0 }; + for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) { + u[gh.idx] = val - gh.flux_sign * correction; + } + for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) { + v[gh.idx] = val - gh.flux_sign * correction; + } + for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) { + w[gh.idx] = val - gh.flux_sign * correction; + } + // The periodic seam: the w face at k = nz is the face at k = 0. + for j in 0..ny { + for i in 0..nx { + let (f0, fn_) = (g.wface(0, j, i), g.wface(nz, j, i)); + if self.w_kind[f0] != FaceKind::Fluid && self.w_kind[fn_] == self.w_kind[f0] { + w[fn_] = w[f0]; + } + } + } + correction + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs index d117e76..c83cdb9 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -11,6 +11,7 @@ pub mod cut; pub mod cutwall; pub mod field; pub mod grid; +pub mod impose; pub mod loads; pub mod poisson; pub mod step; diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs index 746be09..50d5176 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs @@ -4,8 +4,9 @@ //! the wall; its faces carry the mass fluxes averaged from the two adjacent //! cells (the 2D face velocities when every aperture is 1), upwind plus the //! TVD correction as the 2D predictor, apertured diffusion, the pressure -//! force `−(p₊ − p₋) α A` (the projection's gradient), the wall's momentum -//! flux `m_w U_b` with `m_w = −Σ m_f` (so a uniform field stays uniform), +//! force `−(p₊ − p₋) α A` (the projection's gradient), the net mass flux +//! times the face's own value (Reynolds transport; a uniform field stays +//! uniform on any wall motion), //! and the implicit wall shear `μ A_w (u − U_b)/d_f`; the time derivative //! carries the inertia floor. @@ -195,8 +196,12 @@ impl Solver { } }; } - // The wall's momentum flux closes the mass balance exactly. - conv -= mass_out * ub; + // Reynolds transport over a volume whose wall moves with the fluid + // on it: `ρV du/dt = −Σ m (u_face − u)` — the net mass flux of the + // control volume (zero for a body at rest, the swept rate + // otherwise) multiplies the face's own value, so a uniform field + // stays uniform on any wall motion. + conv -= mass_out * u0; let p_plus = lat.cell(cell_plus).map_or(0.0, |ci| field.p[ci]); let p_minus = lat.cell(cell_minus).map_or(0.0, |ci| field.p[ci]); let pressure = -(p_plus - p_minus) * cv.alpha * area[c]; 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 fba038b..43e49a5 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 @@ -550,6 +550,7 @@ impl DeviceStep { tm.cg_iterations += cg_iterations as u64; } StepResult { + fresh_cells: 0, converged: final_residual < self.solver.params.tolerance, corrector_steps_performed: total, final_residual, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index c3e5e31..c40b88f 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -99,6 +99,8 @@ pub struct StepResult { pub final_residual: f64, /// CG iterations summed over the step's projections. pub poisson_iterations: usize, + /// Cells that became fluid on this step (a moving body). + pub fresh_cells: usize, } type Vec3Fn = Box (f64, f64, f64) + Send + Sync>; @@ -109,6 +111,8 @@ pub struct Solver { pub(super) momentum_source: Option, boundary_velocity: Option, body: Option, + /// The body moves: the mask is rebuilt at every step's new time. + moving: bool, mask: Option, last_ghost_correction: f64, wall_fluxes: Vec, @@ -137,6 +141,7 @@ impl Solver { momentum_source: None, boundary_velocity: None, body: None, + moving: false, mask: None, last_ghost_correction: 0.0, wall_fluxes: Vec::new(), @@ -164,9 +169,27 @@ impl Solver { /// A static embedded body (the mask is built at initialisation). pub fn set_body(&mut self, body: Body) { self.body = Some(body); + self.moving = false; self.mask = None; } + /// A moving embedded body: the mask is rebuilt at every step's + /// end-of-step geometry (the 2D solver's order — predictor on the old + /// mask, projection on the new one). + pub fn set_moving_body(&mut self, body: Body) { + self.body = Some(body); + self.moving = true; + self.mask = None; + } + + fn build_mask(&self, body: &Body, g: Grid, t: f64) -> Mask { + match self.params.wall_scheme { + WallScheme::GhostBinary => Mask::build(body, g, t, self.params.boundaries), + WallScheme::CutCell => Mask::build_cut(body, g, t, self.params.boundaries), + } + .expect("embedded mask") + } + #[must_use] pub fn body(&self) -> Option<&Body> { self.body.as_ref() @@ -230,24 +253,52 @@ impl Solver { .is_none_or(|m| m.is_fluid_cell(m.grid().cell(k, j, i))) } - // The apertures (1 without a cut geometry). + // The projection's unknowns and equations (space-time on a moving cut + // wall, the fluid predicates otherwise). + #[inline] + pub(super) fn u_is_unknown(&self, k: usize, j: usize, i: usize) -> bool { + self.mask + .as_ref() + .is_none_or(|m| m.u_open(m.grid().uface(k, j, i))) + } + #[inline] + pub(super) fn v_is_unknown(&self, k: usize, j: usize, i: usize) -> bool { + self.mask + .as_ref() + .is_none_or(|m| m.v_open(m.grid().vface(k, j, i))) + } + #[inline] + pub(super) fn w_is_unknown(&self, k: usize, j: usize, i: usize) -> bool { + self.mask + .as_ref() + .is_none_or(|m| m.w_open(m.grid().wface(k, j, i))) + } + #[inline] + pub(super) fn cell_is_active(&self, k: usize, j: usize, i: usize) -> bool { + self.mask + .as_ref() + .is_none_or(|m| m.cell_active(m.grid().cell(k, j, i))) + } + + // The projection's apertures: step-averaged on a moving cut wall + // (1 without a cut geometry). #[inline] pub(super) fn au(&self, k: usize, j: usize, i: usize) -> f64 { self.mask .as_ref() - .map_or(1.0, |m| m.a_u(m.grid().uface(k, j, i))) + .map_or(1.0, |m| m.au_step(m.grid().uface(k, j, i))) } #[inline] pub(super) fn av(&self, k: usize, j: usize, i: usize) -> f64 { self.mask .as_ref() - .map_or(1.0, |m| m.a_v(m.grid().vface(k, j, i))) + .map_or(1.0, |m| m.av_step(m.grid().vface(k, j, i))) } #[inline] pub(super) fn aw(&self, k: usize, j: usize, i: usize) -> f64 { self.mask .as_ref() - .map_or(1.0, |m| m.a_w(m.grid().wface(k, j, i))) + .map_or(1.0, |m| m.aw_step(m.grid().wface(k, j, i))) } /// The surface velocity's compatible flux through the cell's wall at /// the step's new time (cut wall only; the table is rebuilt per step). @@ -402,15 +453,7 @@ impl Solver { let t = self.time; if let Some(body) = &self.body { if self.mask.is_none() { - let mask = match self.params.wall_scheme { - WallScheme::GhostBinary => { - Mask::build(body, field.grid, t, self.params.boundaries) - } - WallScheme::CutCell => { - Mask::build_cut(body, field.grid, t, self.params.boundaries) - } - }; - self.mask = Some(mask.expect("embedded mask")); + self.mask = Some(self.build_mask(body, field.grid, t)); } } self.apply_boundary_normals(field, t); @@ -434,13 +477,50 @@ impl Solver { field.update_old_values(); self.momentum_predictor(field, dt, t_old); self.apply_boundary_normals(field, t_new); + // A moving body: the mask at the end-of-step geometry, the pressure + // of the cells that just became fluid refilled from their + // neighbours (fluid in both masks), the new mask's prescribed and + // ghost values imposed from the previous corrected field. + let mut fresh_cells = 0; + if self.moving { + if let Some(body) = &self.body { + let mut new_mask = self.build_mask(body, field.grid, t_new); + if let Some(old_mask) = &self.mask { + fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); + new_mask.set_step_apertures(old_mask); + } + new_mask.impose_from( + body, + &field.u_old, + &field.v_old, + &field.w_old, + &mut field.u, + &mut field.v, + &mut field.w, + t_new, + ); + if new_mask.cut().is_some() { + let (table, correction) = match &self.mask { + Some(old_mask) => new_mask.gcl_flux_table(old_mask, dt), + None => new_mask.wall_flux_table(body, t_new), + }; + self.wall_fluxes = table; + self.last_ghost_correction = correction; + } + self.mask = Some(new_mask); + } + } field.copy_to_starred(); let mut cut_correction = None; if let (Some(body), Some(mask)) = (&self.body, &self.mask) { if mask.cut().is_some() { - let (table, correction) = mask.wall_flux_table(body, t_new); - self.wall_fluxes = table; - cut_correction = Some(correction); + if self.moving { + cut_correction = Some(self.last_ghost_correction); + } else { + let (table, correction) = mask.wall_flux_table(body, t_new); + self.wall_fluxes = table; + cut_correction = Some(correction); + } } } let mut total = 0; @@ -450,6 +530,12 @@ impl Solver { let sol = self.solve_correction(field, dt, corrector == 0); poisson_iterations += sol.iterations; let mass_residual = self.apply_correction(field, dt); + if std::env::var_os("RTX_E3_DEBUG").is_some() { + eprintln!( + " corrector {corrector}: CG {} it (converged {}), residual {:.3e}, mass {:.3e}", + sol.iterations, sol.converged, sol.residual, mass_residual + ); + } final_residual = mass_residual; total += 1; if mass_residual < self.params.tolerance { @@ -468,6 +554,65 @@ impl Solver { corrector_steps_performed: total, final_residual, poisson_iterations, + fresh_cells, } } } + +/// Refill the pressure of the cells fluid in `new` and not in `old` from +/// their face neighbours fluid in both; returns their count. +fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { + let g = field.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let periodic = new.periodic_z(); + let mut fresh = 0; + let mut refills = Vec::new(); + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let idx = g.cell(k, j, i); + if !(new.is_fluid_cell(idx) && !old.is_fluid_cell(idx)) { + continue; + } + fresh += 1; + let mut sum = 0.0; + let mut count = 0usize; + let mut visit = |nb: usize| { + if new.is_fluid_cell(nb) && old.is_fluid_cell(nb) { + sum += field.p[nb]; + count += 1; + } + }; + if i + 1 < nx { + visit(g.cell(k, j, i + 1)); + } + if i > 0 { + visit(g.cell(k, j, i - 1)); + } + if j + 1 < ny { + visit(g.cell(k, j + 1, i)); + } + if j > 0 { + visit(g.cell(k, j - 1, i)); + } + if k + 1 < nz { + visit(g.cell(k + 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(0, j, i)); + } + if k > 0 { + visit(g.cell(k - 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(nz - 1, j, i)); + } + if count > 0 { + refills.push((idx, sum / count as f64)); + } + } + } + } + for (idx, p) in refills { + field.p[idx] = p; + } + fresh +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs index 7785281..9ab868d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs @@ -29,7 +29,7 @@ impl Solver { for j in 0..ny { for i in 0..nx { let idx = g.cell(k, j, i); - if !self.cell_is_fluid(k, j, i) { + if !self.cell_is_active(k, j, i) { problem.active[idx] = false; continue; } @@ -38,42 +38,42 @@ impl Solver { if b.x1 == outlet { extra += ae_outlet; } - } else if self.u_is_fluid(k, j, i + 1) { + } else if self.u_is_unknown(k, j, i + 1) { problem.ae[idx] = ae_interior * self.au(k, j, i + 1); } if i == 0 { if b.x0 == outlet { extra += ae_outlet; } - } else if self.u_is_fluid(k, j, i) { + } else if self.u_is_unknown(k, j, i) { problem.aw[idx] = ae_interior * self.au(k, j, i); } if j + 1 == ny { if b.y1 == outlet { extra += an_outlet; } - } else if self.v_is_fluid(k, j + 1, i) { + } else if self.v_is_unknown(k, j + 1, i) { problem.an[idx] = an_interior * self.av(k, j + 1, i); } if j == 0 { if b.y0 == outlet { extra += an_outlet; } - } else if self.v_is_fluid(k, j, i) { + } else if self.v_is_unknown(k, j, i) { problem.as_[idx] = an_interior * self.av(k, j, i); } if k + 1 == nz && !periodic { if b.z1 == outlet { extra += at_outlet; } - } else if self.w_is_fluid((k + 1) % nz, j, i) { + } else if self.w_is_unknown((k + 1) % nz, j, i) { problem.at[idx] = at_interior * self.aw((k + 1) % nz, j, i); } if k == 0 && !periodic { if b.z0 == outlet { extra += at_outlet; } - } else if self.w_is_fluid(k, j, i) { + } else if self.w_is_unknown(k, j, i) { problem.ab[idx] = at_interior * self.aw(k, j, i); } problem.extra_diag[idx] = extra; @@ -252,7 +252,7 @@ impl Solver { for j in 0..ny { for i in 0..nx { let idx = g.cell(k, j, i); - if !self.cell_is_fluid(k, j, i) { + if !self.cell_is_active(k, j, i) { field.sp[idx] = 0.0; continue; } @@ -281,7 +281,7 @@ impl Solver { for k in 0..nz { for j in 0..ny { for i in 0..nx { - if self.cell_is_fluid(k, j, i) { + if self.cell_is_active(k, j, i) { let idx = g.cell(k, j, i); p_prime[idx] = field.p_prime[idx]; } @@ -328,7 +328,7 @@ impl Solver { for k in 0..nz { for j in 0..ny { for i in 1..nx { - if self.u_is_fluid(k, j, i) { + if self.u_is_unknown(k, j, i) { let dp_dx = (pp[g.cell(k, j, i)] - pp[g.cell(k, j, i - 1)]) / dx; let f = g.uface(k, j, i); field.u[f] = field.u_star[f] - (dt / rho) * dp_dx; @@ -347,7 +347,7 @@ impl Solver { } for i in 0..nx { for j in 1..ny { - if self.v_is_fluid(k, j, i) { + if self.v_is_unknown(k, j, i) { let dp_dy = (pp[g.cell(k, j, i)] - pp[g.cell(k, j - 1, i)]) / dy; let f = g.vface(k, j, i); field.v[f] = field.v_star[f] - (dt / rho) * dp_dy; @@ -369,7 +369,7 @@ impl Solver { for i in 0..nx { let k_range = if periodic { 0..nz } else { 1..nz }; for k in k_range { - if self.w_is_fluid(k, j, i) { + if self.w_is_unknown(k, j, i) { let below = if k > 0 { k - 1 } else { nz - 1 }; let dp_dz = (pp[g.cell(k, j, i)] - pp[g.cell(below, j, i)]) / dz; let f = g.wface(k, j, i); @@ -394,7 +394,7 @@ impl Solver { for k in 0..nz { for j in 0..ny { for i in 0..nx { - if self.cell_is_fluid(k, j, i) { + if self.cell_is_active(k, j, i) { let idx = g.cell(k, j, i); field.p[idx] += pp[idx]; } @@ -405,7 +405,7 @@ impl Solver { for k in 0..nz { for j in 0..ny { for i in 0..nx { - if !self.cell_is_fluid(k, j, i) { + if !self.cell_is_active(k, j, i) { continue; } let idx = g.cell(k, j, i); diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs index 9ecd473..723f1d4 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -47,7 +47,7 @@ pub(crate) struct StencilNode { #[derive(Debug, Clone)] pub(super) struct Ghost { - idx: usize, + pub(super) idx: usize, x: f64, y: f64, z: f64, @@ -58,7 +58,7 @@ pub(super) struct Ghost { nodes: Vec, /// Outward-from-fluid sign for the compatibility correction (0 when no /// fluid cell is adjacent). - flux_sign: f64, + pub(super) flux_sign: f64, } #[derive(Clone)] @@ -77,6 +77,17 @@ pub struct Mask { /// The cut geometry of the apertured wall (`WallScheme::CutCell`, /// `cutwall.rs`); `None` on the binary ghost wall. pub(super) cut: Option, + /// The step-averaged apertures `½(αⁿ + αⁿ⁺¹)` of a moving cut wall + /// (the space-time continuity: a cell's volume change over the step + /// equals the flux through the apertures it had during it); `None` = + /// the instantaneous ones. + pub(super) step_apertures: Option<(Vec, Vec, Vec)>, + /// The projection's space-time classification on a moving cut wall: + /// a face is an unknown where its step-averaged aperture is positive, + /// a cell has an equation where it holds fluid at either end of the + /// step (a dying cell empties through the apertures it had); `None` = + /// the instantaneous kinds. + pub(super) step_open: Option<(Vec, Vec, Vec, Vec)>, } /// The z lattice position of a query: the lower plane index, the upper @@ -250,7 +261,7 @@ pub(crate) fn linear_fit(pts_w: &[(f64, f64, f64, f64, f64)], at: (f64, f64, f64 } impl Ghost { - fn reconstruct(&self, values: &[f64]) -> f64 { + pub(super) fn reconstruct(&self, values: &[f64]) -> f64 { let mut pts: Vec<(f64, f64, f64, f64, f64)> = self .nodes .iter() @@ -484,9 +495,65 @@ impl Mask { anchor, fluid_cells, cut: None, + step_apertures: None, + step_open: None, }) } + // The projection's unknowns (the instantaneous kinds at rest). + #[inline] + #[must_use] + pub fn u_open(&self, idx: usize) -> bool { + self.step_open + .as_ref() + .map_or(self.u_kind[idx] == FaceKind::Fluid, |o| o.0[idx]) + } + #[inline] + #[must_use] + pub fn v_open(&self, idx: usize) -> bool { + self.step_open + .as_ref() + .map_or(self.v_kind[idx] == FaceKind::Fluid, |o| o.1[idx]) + } + #[inline] + #[must_use] + pub fn w_open(&self, idx: usize) -> bool { + self.step_open + .as_ref() + .map_or(self.w_kind[idx] == FaceKind::Fluid, |o| o.2[idx]) + } + #[inline] + #[must_use] + pub fn cell_active(&self, idx: usize) -> bool { + self.step_open + .as_ref() + .map_or(self.cell_fluid[idx], |o| o.3[idx]) + } + + /// The step-averaged aperture of a u / v / w face (the instantaneous + /// one for a wall at rest). + #[inline] + #[must_use] + pub fn au_step(&self, idx: usize) -> f64 { + self.step_apertures + .as_ref() + .map_or_else(|| self.a_u(idx), |a| a.0[idx]) + } + #[inline] + #[must_use] + pub fn av_step(&self, idx: usize) -> f64 { + self.step_apertures + .as_ref() + .map_or_else(|| self.a_v(idx), |a| a.1[idx]) + } + #[inline] + #[must_use] + pub fn aw_step(&self, idx: usize) -> f64 { + self.step_apertures + .as_ref() + .map_or_else(|| self.a_w(idx), |a| a.2[idx]) + } + /// The cut geometry (apertured wall only). #[must_use] pub fn cut(&self) -> Option<&CutGeometry> { @@ -555,134 +622,4 @@ impl Mask { pub fn periodic_z(&self) -> bool { self.periodic_z } - - /// Impose the wall on `(u, v, w)` from the same field. - pub fn impose(&self, body: &Body, u: &mut [f64], v: &mut [f64], w: &mut [f64], t: f64) -> f64 { - let (us, vs, ws) = (u.to_vec(), v.to_vec(), w.to_vec()); - self.impose_from(body, &us, &vs, &ws, u, v, w, t) - } - - /// Solid faces: the surface velocity; ghost faces: the reconstruction - /// from the SOURCE field, minus the shared flux compatibility - /// correction over the flux-carrying ghosts. Returns the correction. - #[allow(clippy::too_many_arguments)] - pub fn impose_from( - &self, - body: &Body, - u_src: &[f64], - v_src: &[f64], - w_src: &[f64], - u: &mut [f64], - v: &mut [f64], - w: &mut [f64], - t: f64, - ) -> f64 { - let g = self.grid; - let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); - for k in 0..nz { - for j in 0..ny { - for i in 1..nx { - let idx = g.uface(k, j, i); - if self.u_kind[idx] == FaceKind::Solid { - u[idx] = body - .surface_velocity( - i as f64 * dx, - (j as f64 + 0.5) * dy, - (k as f64 + 0.5) * dz, - t, - ) - .0; - } - } - } - for j in 1..ny { - for i in 0..nx { - let idx = g.vface(k, j, i); - if self.v_kind[idx] == FaceKind::Solid { - v[idx] = body - .surface_velocity( - (i as f64 + 0.5) * dx, - j as f64 * dy, - (k as f64 + 0.5) * dz, - t, - ) - .1; - } - } - } - } - for k in 0..=nz { - for j in 0..ny { - for i in 0..nx { - let idx = g.wface(k, j, i); - if self.w_kind[idx] == FaceKind::Solid { - w[idx] = body - .surface_velocity( - (i as f64 + 0.5) * dx, - (j as f64 + 0.5) * dy, - k as f64 * dz, - t, - ) - .2; - } - } - } - } - let u_vals: Vec = self - .u_ghosts - .iter() - .map(|gh| gh.reconstruct(u_src)) - .collect(); - let v_vals: Vec = self - .v_ghosts - .iter() - .map(|gh| gh.reconstruct(v_src)) - .collect(); - let w_vals: Vec = self - .w_ghosts - .iter() - .map(|gh| gh.reconstruct(w_src)) - .collect(); - let (au, av, aw) = (dy * dz, dx * dz, dx * dy); - let mut net = 0.0; - let mut area = 0.0; - for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) { - if gh.flux_sign != 0.0 { - net += gh.flux_sign * val * au; - area += au; - } - } - for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) { - if gh.flux_sign != 0.0 { - net += gh.flux_sign * val * av; - area += av; - } - } - for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) { - if gh.flux_sign != 0.0 { - net += gh.flux_sign * val * aw; - area += aw; - } - } - let correction = if area > 0.0 { net / area } else { 0.0 }; - for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) { - u[gh.idx] = val - gh.flux_sign * correction; - } - for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) { - v[gh.idx] = val - gh.flux_sign * correction; - } - for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) { - w[gh.idx] = val - gh.flux_sign * correction; - } - // The periodic seam: the w face at k = nz is the face at k = 0. - for j in 0..ny { - for i in 0..nx { - let (f0, fn_) = (g.wface(0, j, i), g.wface(nz, j, i)); - if self.w_kind[f0] != FaceKind::Fluid && self.w_kind[fn_] == self.w_kind[f0] { - w[fn_] = w[f0]; - } - } - } - correction - } } diff --git a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs index b21fc95..78d2822 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs @@ -11,255 +11,10 @@ //! are at most the binary wall's at every n, its loads within 10 % at the //! finest rung. -use rtx_cfd::solvers::incompressible::ConvectionScheme; -use rtx_cfd::solvers::incompressible::embedded3::{ - Body, FaceKind, Field, Fluid, Grid, Parameters, Solver, WallScheme, -}; -use std::f64::consts::PI; +mod embedded3_sphere; -const RHO: f64 = 1.0; -const MU: f64 = 0.05; -const C: (f64, f64, f64) = (0.6, 0.45, 0.5); -const R: f64 = 0.2; - -fn u3(x: f64, y: f64, z: f64) -> f64 { - (PI * x).sin() * (PI * y).cos() * (PI * z).cos() -} -fn v3(x: f64, y: f64, z: f64) -> f64 { - (PI * x).cos() * (PI * y).sin() * (PI * z).cos() -} -fn w3(x: f64, y: f64, z: f64) -> f64 { - -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin() -} -fn p3(x: f64, y: f64, z: f64) -> f64 { - (PI * x).sin() * (PI * y).sin() * (PI * z).sin() -} -/// The velocity gradient ∂u_i/∂x_j and the pressure gradient. -fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) { - let (sx, cx) = (PI * x).sin_cos(); - let (sy, cy) = (PI * y).sin_cos(); - let (sz, cz) = (PI * z).sin_cos(); - ( - [ - [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz], - [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz], - [ - 2.0 * PI * sx * cy * sz, - 2.0 * PI * cx * sy * sz, - -2.0 * PI * cx * cy * cz, - ], - ], - [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz], - ) -} -fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let (g, gp) = grads(x, y, z); - let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)]; - let lap = -3.0 * PI * PI; - let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2]; - ( - RHO * conv(0) + gp[0] - MU * lap * u[0], - RHO * conv(1) + gp[1] - MU * lap * u[1], - RHO * conv(2) + gp[2] - MU * lap * u[2], - ) -} -fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) { - let u = if x <= 0.0 || x >= 1.0 { - 0.0 - } else { - u3(x, y, z) - }; - let v = if y <= 0.0 || y >= 1.0 { - 0.0 - } else { - v3(x, y, z) - }; - let w = if z <= 0.0 || z >= 1.0 { - 0.0 - } else { - w3(x, y, z) - }; - (u, v, w) -} - -/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA` -/// over the sphere by a fine Fibonacci quadrature. -fn exact_force_and_flux() -> ([f64; 3], [f64; 3]) { - let n = 200_000; - let golden = PI * (3.0 - 5.0_f64.sqrt()); - let (mut f, mut m) = ([0.0; 3], [0.0; 3]); - let da = 4.0 * PI * R * R / n as f64; - for k in 0..n { - let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64; - let rr = (1.0 - zz * zz).sqrt(); - let th = golden * k as f64; - let nrm = [rr * th.cos(), rr * th.sin(), zz]; - let (x, y, z) = (C.0 + R * nrm[0], C.1 + R * nrm[1], C.2 + R * nrm[2]); - let (g, _) = grads(x, y, z); - let p = p3(x, y, z); - let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)]; - let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2]; - for i in 0..3 { - let mut t = -p * nrm[i]; - for j in 0..3 { - t += MU * (g[i][j] + g[j][i]) * nrm[j]; - } - f[i] += t * da; - m[i] += RHO * u[i] * un * da; - } - } - (f, m) -} - -struct Measurement { - l2_velocity: f64, - max_div: f64, - ghost_correction: f64, - force_surface: [f64; 3], - skipped: usize, - force_cv: [f64; 3], -} - -fn measure(n: usize, scheme: WallScheme) -> Measurement { - let h = 1.0 / n as f64; - let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h); - let mut solver = Solver::new( - Fluid { - density: RHO, - viscosity: MU, - reference_velocity: 1.0, - reference_length: 1.0, - }, - Parameters { - corrector_steps: 2, - tolerance: 1e-8, - convection_scheme: ConvectionScheme::Upwind, - wall_scheme: scheme, - ..Parameters::default() - }, - ); - solver.set_momentum_source(|x, y, z, _t| source3(x, y, z)); - solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z)); - solver.set_body( - Body::sphere(|_t| C, R) - .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))), - ); - let g = Grid::cubic(n, n, n, h); - let mut f = Field::new(g); - solver.initialize(&mut f); - let mut last = solver.advance(&mut f, dt); - for _ in 0..200_000 { - let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone()); - last = solver.advance(&mut f, dt); - let mut change = 0.0_f64; - for (a, b) in - f.u.iter() - .zip(&bu) - .chain(f.v.iter().zip(&bv)) - .chain(f.w.iter().zip(&bw)) - { - change = change.max((a - b).abs()); - } - if change / dt < 1e-6 { - break; - } - } - let mask = solver.mask().expect("mask"); - let (mut sq, mut vol) = (0.0, 0.0); - let dv = h * h * h; - for k in 0..n { - for j in 0..n { - for i in 1..n { - if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid { - let e = f.u[g.uface(k, j, i)] - - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h); - sq += e * e * dv; - vol += dv; - } - } - } - for j in 1..n { - for i in 0..n { - if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid { - let e = f.v[g.vface(k, j, i)] - - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h); - sq += e * e * dv; - vol += dv; - } - } - } - } - for k in 1..n { - for j in 0..n { - for i in 0..n { - if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid { - let e = f.w[g.wface(k, j, i)] - - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h); - sq += e * e * dv; - vol += dv; - } - } - } - } - let body = solver.body().expect("body"); - let t = solver.time(); - // The apertured divergence per unit volume, the porous surface's flux - // through the wall included (the plain divergence on the binary wall). - let mut max_div = 0.0_f64; - let mut at_vol = 1.0; - let mut sum_flux = 0.0; - let (wall_fluxes, _) = mask.wall_flux_table(body, t); - for k in 0..n { - for j in 0..n { - for i in 0..n { - let idx = g.cell(k, j, i); - if mask.is_fluid_cell(idx) { - let flux = (mask.a_u(g.uface(k, j, i + 1)) * f.u[g.uface(k, j, i + 1)] - - mask.a_u(g.uface(k, j, i)) * f.u[g.uface(k, j, i)]) - * h - * h - + (mask.a_v(g.vface(k, j + 1, i)) * f.v[g.vface(k, j + 1, i)] - - mask.a_v(g.vface(k, j, i)) * f.v[g.vface(k, j, i)]) - * h - * h - + (mask.a_w(g.wface(k + 1, j, i)) * f.w[g.wface(k + 1, j, i)] - - mask.a_w(g.wface(k, j, i)) * f.w[g.wface(k, j, i)]) - * h - * h - + wall_fluxes[idx]; - sum_flux += flux.abs(); - if (flux / (h * h * h)).abs() > max_div { - max_div = (flux / (h * h * h)).abs(); - at_vol = mask.vol(idx); - } - } - } - } - } - println!( - " [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}", - last.final_residual - ); - let surface = match scheme { - WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h), - WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce { - f: mask.cut_wall_force(body, &f, MU, t).expect("cut wall"), - samples: 0, - skipped: 0, - }, - }; - let (i0, i1) = (n / 8, n - n / 8); - let src = |x: f64, y: f64, z: f64| source3(x, y, z); - let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1)); - Measurement { - l2_velocity: (sq / vol).sqrt(), - max_div, - ghost_correction: solver.ghost_correction().abs(), - force_surface: surface.f, - skipped: surface.skipped, - force_cv, - } -} +use embedded3_sphere::{C, Measurement, exact_force_and_flux, measure}; +use rtx_cfd::solvers::incompressible::embedded3::WallScheme; fn norm(a: [f64; 3]) -> f64 { (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt() @@ -273,13 +28,13 @@ struct Ladder { } fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder { - let (fe, m) = exact_force_and_flux(); + let (fe, m) = exact_force_and_flux(C); let f_scale = norm(fe); let fcv = [fe[0] - m[0], fe[1] - m[1], fe[2] - m[2]]; println!( " {scheme:?}: exact force {fe:.5?}; momentum flux {m:.5?}; the control-volume route measures {fcv:.5?}" ); - let ms: Vec = resolutions.iter().map(|&n| measure(n, scheme)).collect(); + let ms: Vec = resolutions.iter().map(|&n| measure(n, scheme, C)).collect(); let errors: Vec = ms.iter().map(|x| x.l2_velocity).collect(); let mut se = Vec::new(); let mut ce = Vec::new(); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs new file mode 100644 index 0000000..a1d6255 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs @@ -0,0 +1,566 @@ +//! embedded3 item 11a: the fresh-cell falsifier of the 2D track +//! (`embedded_fresh_cell_falsifier.rs`, omni-cortex +//! `docs/fresh_cell_gcl_campaign.md`) on the 3D solver, per unit span — +//! the rigid Turek–Hron flag (0.35 × 0.02 m) extruded across a periodic +//! slab, oscillating transversely in still fluid at the flag's tip speed +//! (1 m/s peak, 80 mm amplitude) on h = 1/152 at dt = 3.24e-4. Per step: +//! the load per unit span (the ghost wall's traction route over the +//! plate's samples; the cut wall's operator route), a far-field pressure +//! probe, the fluid's kinetic energy, the fresh-cell count. +//! +//! Registered gates (`docs/embedded3_campaign.md` item 11): +//! - GhostBinary reproduces the 2D wall's impulse: energy per flipped +//! column within 30 % of the 2D 0.048 J/m per flipped cell, spike RMS +//! exponent in dt ≈ −1 (published −0.8 for the raw volume source); +//! - CutCell: energy per fresh column ≥ 20× lower, max force spike < 5 % +//! of ½ρU²L, exponent ∈ [−0.3, 0.3]. +//! +//! Default run: dt only, both schemes (minutes on the host); +//! `RTX_E3_FALSIFIER_LADDER=1` runs dt, dt/2, dt/4 and fits the exponent +//! (the gated variant is `#[ignore]`); `RTX_E3_FALSIFIER_NZ` sets the +//! span in cells (default 4); `RTX_E3_FALSIFIER_CSV=` dumps records. + +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; +use std::io::Write as _; + +const RHO: f64 = 1000.0; +const MU: f64 = 1.0; +const N: usize = 152; +const DT_FSI2: f64 = 3.24e-4; +const HX: f64 = 0.175; +const HY: f64 = 0.01; +const AMP: f64 = 0.08; +const U_PEAK: f64 = 1.0; +const CX: f64 = 0.5; +const CY0: f64 = 0.5; +/// The 2D wall's measured energy per flipped cell (J/m at U = 1, h = 1/152). +const ENERGY_2D: f64 = 0.048; + +fn span_cells() -> usize { + std::env::var("RTX_E3_FALSIFIER_NZ") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(4) +} + +fn center_y(t: f64) -> f64 { + CY0 + AMP * (U_PEAK / AMP * t).sin() +} + +fn center_v(t: f64) -> f64 { + U_PEAK * (U_PEAK / AMP * t).cos() +} + +fn plate_sdf(x: f64, y: f64, yc: f64) -> f64 { + let qx = (x - CX).abs() - HX; + let qy = (y - yc).abs() - HY; + let outside = (qx.max(0.0).powi(2) + qy.max(0.0).powi(2)).sqrt(); + outside + qx.max(qy).min(0.0) +} + +const R_CIRCLE: f64 = 0.05; + +/// `RTX_E3_FALSIFIER_BODY=circle`: the 2D falsifier's smooth body (a +/// cylinder of radius 0.05 across the span) instead of the plate. +fn circle_body() -> bool { + std::env::var("RTX_E3_FALSIFIER_BODY").is_ok_and(|v| v == "circle") +} + +/// `RTX_E3_FALSIFIER_BODY=stadium`: the plate with semicircular ends +/// (radius `HY`): the same length and thickness, a smooth interface for +/// the cut geometry's linear interpolant. +fn stadium_body() -> bool { + std::env::var("RTX_E3_FALSIFIER_BODY").is_ok_and(|v| v == "stadium") +} + +fn stadium_sdf(x: f64, y: f64, yc: f64) -> f64 { + let half = HX - HY; + let qx = (x - CX).abs().max(half) - half; + (qx * qx + (y - yc).powi(2)).sqrt() - HY +} + +fn plate(moving: bool) -> Body { + let yc = move |t: f64| if moving { center_y(t) } else { CY0 }; + let vc = move |t: f64| if moving { center_v(t) } else { 0.0 }; + if circle_body() { + return Body::from_sdf(move |x, y, _z, t| { + ((x - CX).powi(2) + (y - yc(t)).powi(2)).sqrt() - R_CIRCLE + }) + .with_surface_velocity(move |_, _, _, t| (0.0, vc(t), 0.0)); + } + if stadium_body() { + return Body::from_sdf(move |x, y, _z, t| stadium_sdf(x, y, yc(t))) + .with_surface_velocity(move |_, _, _, t| (0.0, vc(t), 0.0)); + } + Body::from_sdf(move |x, y, _z, t| plate_sdf(x, y, yc(t))) + .with_surface_velocity(move |_, _, _, t| (0.0, vc(t), 0.0)) +} + +/// The load scale `½ρU²L` of the body (its length across the motion). +fn load_scale() -> f64 { + let l = if circle_body() { + 2.0 * R_CIRCLE + } else { + 2.0 * HX + }; + 0.5 * RHO * U_PEAK * U_PEAK * l +} + +/// Surface samples of the plate at `t`: `(x, y, z, nx, ny, area)` over the +/// four edges at spacing `ds` and `nz` z levels. +fn samples( + t: f64, + moving: bool, + ds: f64, + nz: usize, + dz: f64, +) -> Vec<(f64, f64, f64, f64, f64, f64)> { + let yc = if moving { center_y(t) } else { CY0 }; + let mut out = Vec::new(); + if circle_body() { + let n = ((2.0 * std::f64::consts::PI * R_CIRCLE / ds).ceil() as usize).max(8); + let dth = 2.0 * std::f64::consts::PI / n as f64; + for k in 0..n { + let th = (k as f64 + 0.5) * dth; + let (sn, cs) = th.sin_cos(); + for kz in 0..nz { + out.push(( + CX + R_CIRCLE * cs, + yc + R_CIRCLE * sn, + (kz as f64 + 0.5) * dz, + cs, + sn, + R_CIRCLE * dth * dz, + )); + } + } + return out; + } + if stadium_body() { + let half = HX - HY; + let n_flat = ((2.0 * half / ds).ceil() as usize).max(1); + for k in 0..n_flat { + let x = CX - half + (k as f64 + 0.5) / n_flat as f64 * 2.0 * half; + for kz in 0..nz { + let z = (kz as f64 + 0.5) * dz; + let a = 2.0 * half / n_flat as f64 * dz; + out.push((x, yc + HY, z, 0.0, 1.0, a)); + out.push((x, yc - HY, z, 0.0, -1.0, a)); + } + } + let n_arc = ((std::f64::consts::PI * HY / ds).ceil() as usize).max(4); + for (cx, sign) in [(CX + half, 1.0), (CX - half, -1.0)] { + for k in 0..n_arc { + let th = -std::f64::consts::FRAC_PI_2 + + (k as f64 + 0.5) / n_arc as f64 * std::f64::consts::PI; + let (sn, cs) = th.sin_cos(); + let (nx, ny) = (sign * cs, sn); + for kz in 0..nz { + out.push(( + cx + HY * nx, + yc + HY * ny, + (kz as f64 + 0.5) * dz, + nx, + ny, + std::f64::consts::PI * HY / n_arc as f64 * dz, + )); + } + } + } + return out; + } + let (x0, x1, y0, y1) = (CX - HX, CX + HX, yc - HY, yc + HY); + let mut edge = |ax: f64, ay: f64, bx: f64, by: f64, nx: f64, ny: f64| { + let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); + let n = ((len / ds).ceil() as usize).max(1); + for k in 0..n { + let s = (k as f64 + 0.5) / n as f64; + for kz in 0..nz { + out.push(( + ax + s * (bx - ax), + ay + s * (by - ay), + (kz as f64 + 0.5) * dz, + nx, + ny, + len / n as f64 * dz, + )); + } + } + }; + edge(x0, y0, x1, y0, 0.0, -1.0); + edge(x1, y0, x1, y1, 1.0, 0.0); + edge(x1, y1, x0, y1, 0.0, 1.0); + edge(x0, y1, x0, y0, -1.0, 0.0); + out +} + +struct Record { + t: f64, + /// Load per unit span (the scheme's wall route). + fy: f64, + /// Load per unit span by the control-volume route (a box of whole + /// cells around the body, reading no near-wall value). + fy_cv: f64, + fresh: usize, + skipped: usize, + p_far: f64, + /// Kinetic energy per unit span over the fluid cells. + ke: f64, +} + +struct Run { + records: Vec, + /// The largest kinetic-energy change per step at a step with fresh + /// cells (after the impulsive start) over that step's flipped columns + /// (J/m) — the 2D falsifier's 2.604 J/m over 54 cells = 0.048. + energy_per_flip: f64, + seconds: f64, +} + +fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { + let nz = span_cells(); + let h = 1.0 / N as f64; + let lz = nz as f64 * h; + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 2.0 * HY, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: scheme, + boundaries: Boundaries { + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|_, _, _, _| (0.0, 0.0, 0.0)); + if moving { + solver.set_moving_body(plate(true)); + } else { + solver.set_body(plate(false)); + } + let g = Grid::cubic(N, N, nz, h); + let mut field = Field::new(g); + solver.initialize(&mut field); + let steps = (t_end / dt).round() as usize; + let mut records = Vec::with_capacity(steps); + // The 2D definition: the largest |ΔKE| step's energy over that + // step's flipped columns. + let mut largest_jump = 0.0_f64; + let mut energy_per_flip = 0.0_f64; + let mut ke_prev: Option = None; + let start = std::time::Instant::now(); + let (jp, ip, kp) = ( + (0.92 * N as f64) as usize, + (0.5 * N as f64) as usize, + nz / 2, + ); + for step in 0..steps { + let result = solver.advance(&mut field, dt); + let t = (step + 1) as f64 * dt; + let mask = solver.mask().expect("mask"); + let body = solver.body().expect("body"); + let (mut fy, mut skipped) = (0.0, 0usize); + match scheme { + WallScheme::GhostBinary => { + for (x, y, z, nx, ny, area) in samples(t, moving, 0.5 * h, nz, h) { + match mask.traction_at(body, &field, MU, t, [x, y, z], [nx, ny, 0.0]) { + Some(tr) => fy += tr[1] * area, + None => skipped += 1, + } + } + } + WallScheme::CutCell => { + fy = mask.cut_wall_force(body, &field, MU, t).expect("cut wall")[1]; + } + } + fy /= lz; + let margin = 8; + let fy_cv = mask.control_volume_force( + &field, + dt, + RHO, + MU, + None, + (margin, N - margin, margin, N - margin, 0, nz), + )[1] / lz; + let p_far = field.p[g.cell(kp, jp, ip)]; + let mut ke = 0.0; + for k in 0..nz { + for j in 0..N { + for i in 0..N { + let idx = g.cell(k, j, i); + if mask.is_fluid_cell(idx) { + let uc = 0.5 * (field.u[g.uface(k, j, i)] + field.u[g.uface(k, j, i + 1)]); + let vc = 0.5 * (field.v[g.vface(k, j, i)] + field.v[g.vface(k, j + 1, i)]); + let wc = 0.5 * (field.w[g.wface(k, j, i)] + field.w[g.wface(k + 1, j, i)]); + ke += 0.5 * RHO * (uc * uc + vc * vc + wc * wc) * h * h * h * mask.vol(idx); + } + } + } + } + ke /= lz; + if let Some(prev) = ke_prev { + if step > 30 && result.fresh_cells > 0 && (ke - prev).abs() > largest_jump { + largest_jump = (ke - prev).abs(); + // The plate's event is its row (the 2D divided by the row's + // 54 cells); the circle's is the step's fresh columns. + let columns = if circle_body() { + result.fresh_cells as f64 / nz as f64 + } else { + (2.0 * HX / h).round() + }; + energy_per_flip = largest_jump / columns; + } + } + ke_prev = Some(ke); + records.push(Record { + t, + fy, + fy_cv, + fresh: result.fresh_cells, + skipped, + p_far, + ke, + }); + } + Run { + records, + energy_per_flip, + seconds: start.elapsed().as_secs_f64(), + } +} + +/// Spike series: the load minus its 21-step running median. +fn spikes(f: &[f64]) -> Vec { + let w = 10usize; + (0..f.len()) + .map(|k| { + let lo = k.saturating_sub(w); + let hi = (k + w + 1).min(f.len()); + let mut win: Vec = f[lo..hi].to_vec(); + win.sort_by(|a, b| a.partial_cmp(b).unwrap()); + f[k] - win[win.len() / 2] + }) + .collect() +} + +struct Stats { + rms_force: f64, + rms_spike: f64, + max_spike: f64, + rms_spike_cv: f64, + max_spike_cv: f64, + rms_pfar_spike: f64, + max_pfar_spike: f64, + max_ke_jump: f64, + fresh_total: usize, + skipped_max: usize, +} + +fn stats(records: &[Record], t_lo: f64, t_hi: f64) -> Stats { + let fy: Vec = records.iter().map(|r| r.fy).collect(); + let sp = spikes(&fy); + let fcv: Vec = records.iter().map(|r| r.fy_cv).collect(); + let spc = spikes(&fcv); + let pf: Vec = records.iter().map(|r| r.p_far).collect(); + let spf = spikes(&pf); + let idx: Vec = (0..records.len()) + .filter(|&k| records[k].t >= t_lo && records[k].t <= t_hi) + .collect(); + let rms = |v: &dyn Fn(usize) -> f64| { + (idx.iter().map(|&k| v(k) * v(k)).sum::() / idx.len().max(1) as f64).sqrt() + }; + Stats { + rms_force: rms(&|k| fy[k]), + rms_spike: rms(&|k| sp[k]), + max_spike: idx.iter().map(|&k| sp[k].abs()).fold(0.0, f64::max), + rms_spike_cv: rms(&|k| spc[k]), + max_spike_cv: idx.iter().map(|&k| spc[k].abs()).fold(0.0, f64::max), + rms_pfar_spike: rms(&|k| spf[k]), + max_pfar_spike: idx.iter().map(|&k| spf[k].abs()).fold(0.0, f64::max), + max_ke_jump: idx + .iter() + .filter(|&&k| k > 0) + .map(|&k| (records[k].ke - records[k - 1].ke).abs()) + .fold(0.0, f64::max), + fresh_total: idx.iter().map(|&k| records[k].fresh).sum(), + skipped_max: idx.iter().map(|&k| records[k].skipped).max().unwrap_or(0), + } +} + +fn dump(dir: &str, name: &str, records: &[Record]) { + let path = std::path::Path::new(dir).join(format!("{name}.csv")); + let mut f = std::fs::File::create(path).expect("csv"); + writeln!(f, "t,fy,fy_cv,fresh,skipped,p_far,ke").unwrap(); + for r in records { + writeln!( + f, + "{:.6},{:.6e},{:.6e},{},{},{:.6e},{:.6e}", + r.t, r.fy, r.fy_cv, r.fresh, r.skipped, r.p_far, r.ke + ) + .unwrap(); + } +} + +struct Verdict { + energy_per_flip: f64, + max_spike: f64, + exponent: Option, +} + +fn falsify(scheme: WallScheme, ladder: bool) -> Verdict { + let csv_dir = std::env::var("RTX_E3_FALSIFIER_CSV").ok(); + let period = 2.0 * std::f64::consts::PI * AMP / U_PEAK; + let t_end = 0.3 * period; + let (t_lo, t_hi) = (0.02 * period, 0.28 * period); + let rest = run(scheme, false, DT_FSI2, t_end); + let s0 = stats(&rest.records, t_lo, t_hi); + println!( + " {scheme:?} plate AT REST, dt {DT_FSI2:.2e} ({:.0} s): rms force {:.3e}, rms spike {:.3e}, max spike {:.3e}, fresh {}, skipped max {}", + rest.seconds, s0.rms_force, s0.rms_spike, s0.max_spike, s0.fresh_total, s0.skipped_max + ); + if let Some(d) = &csv_dir { + dump(d, &format!("{scheme:?}_rest"), &rest.records); + } + let dts: Vec = if ladder { + vec![DT_FSI2, DT_FSI2 / 2.0, DT_FSI2 / 4.0] + } else { + vec![DT_FSI2] + }; + let mut points = Vec::new(); + let mut energy = 0.0_f64; + let mut max_spike = 0.0_f64; + for &dt in &dts { + let r = run(scheme, true, dt, t_end); + let s = stats(&r.records, t_lo, t_hi); + println!( + " {scheme:?} plate MOVING, dt {dt:.3e} ({} steps, {:.0} s): rms force {:.3e}, rms spike {:.3e} ({:.1}x rest), max spike {:.3e} N/m ({:.2e} of ½ρU²L), fresh cells {} ({:.2}/step), skipped max {}", + r.records.len(), + r.seconds, + s.rms_force, + s.rms_spike, + s.rms_spike / s0.rms_spike.max(1e-300), + s.max_spike, + s.max_spike / load_scale(), + s.fresh_total, + s.fresh_total as f64 / r.records.len() as f64, + s.skipped_max + ); + println!( + " control-volume route: rms spike {:.3e}, max spike {:.3e} N/m ({:.2e} of ½ρU²L)", + s.rms_spike_cv, + s.max_spike_cv, + s.max_spike_cv / load_scale() + ); + println!( + " far probe p(0.5, 0.92): rms spike {:.3e}, max spike {:.3e}; max |ΔKE| per step {:.3e} J/m; energy per flipped column {:.3e} J/m ({:.2} of the 2D wall's {ENERGY_2D})", + s.rms_pfar_spike, + s.max_pfar_spike, + s.max_ke_jump, + r.energy_per_flip, + r.energy_per_flip / ENERGY_2D + ); + if let Some(d) = &csv_dir { + dump(d, &format!("{scheme:?}_moving_dt{dt:.3e}"), &r.records); + } + assert!(s.rms_force.is_finite() && s.rms_spike.is_finite()); + if dt == DT_FSI2 { + energy = r.energy_per_flip; + max_spike = s.max_spike; + } + points.push((dt, s.rms_spike)); + } + let exponent = (points.len() >= 2).then(|| { + let xs: Vec = points.iter().map(|p| p.0.ln()).collect(); + let ys: Vec = points.iter().map(|p| p.1.ln()).collect(); + let mx = xs.iter().sum::() / xs.len() as f64; + let my = ys.iter().sum::() / ys.len() as f64; + let num: f64 = xs.iter().zip(&ys).map(|(x, y)| (x - mx) * (y - my)).sum(); + let den: f64 = xs.iter().map(|x| (x - mx).powi(2)).sum(); + let e = num / den; + println!( + " {scheme:?}: spike RMS ~ (dt)^{e:.2} across {} time steps", + points.len() + ); + e + }); + Verdict { + energy_per_flip: energy, + max_spike, + exponent, + } +} + +#[test] +fn oscillating_plate_both_walls() { + let ladder = std::env::var("RTX_E3_FALSIFIER_LADDER").is_ok(); + println!( + " body: {}", + if circle_body() { + "circle R 0.05" + } else if stadium_body() { + "stadium 0.35 x 0.02 (semicircular ends)" + } else { + "plate 0.35 x 0.02" + } + ); + let ghost = falsify(WallScheme::GhostBinary, ladder); + let cut = falsify(WallScheme::CutCell, ladder); + println!( + " energy per flipped column: ghost {:.3e}, cut {:.3e} (ratio {:.1}x); max spike: ghost {:.3e}, cut {:.3e} N/m", + ghost.energy_per_flip, + cut.energy_per_flip, + ghost.energy_per_flip / cut.energy_per_flip.max(1e-300), + ghost.max_spike, + cut.max_spike + ); + assert!( + ghost.energy_per_flip > 0.0, + "the binary wall must flip cells" + ); +} + +/// The registered gates on the dt ladder. +#[test] +#[ignore = "item 11's gated ladder (dt, dt/2, dt/4 on both walls; tens of minutes on the host)"] +fn oscillating_plate_gates() { + let ghost = falsify(WallScheme::GhostBinary, true); + let cut = falsify(WallScheme::CutCell, true); + let ratio = ghost.energy_per_flip / cut.energy_per_flip.max(1e-300); + println!( + " GATES: ghost energy per flipped column {:.3e} ({:.2} of 2D), exponent {:.2}; cut energy {:.3e} ({:.1}x lower), max spike {:.3e} N/m ({:.2e} of ½ρU²L), exponent {:.2}", + ghost.energy_per_flip, + ghost.energy_per_flip / ENERGY_2D, + ghost.exponent.unwrap(), + cut.energy_per_flip, + ratio, + cut.max_spike, + cut.max_spike / load_scale(), + cut.exponent.unwrap() + ); + let g2d = ghost.energy_per_flip / ENERGY_2D; + assert!( + (0.7..=1.3).contains(&g2d), + "ghost energy per flip {g2d:.2} of 2D" + ); + assert!(ratio >= 20.0, "cut energy only {ratio:.1}x lower"); + assert!( + cut.max_spike < 0.05 * load_scale(), + "cut max spike {:.3e}", + cut.max_spike + ); + let e = cut.exponent.unwrap(); + assert!((-0.3..=0.3).contains(&e), "cut exponent {e:.2}"); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs new file mode 100644 index 0000000..72177e4 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs @@ -0,0 +1,255 @@ +//! The embedded-sphere manufactured solution shared by the embedded3 +//! wall gates (items 9–11): the fields, the source, the boundary data, +//! the exact surface integrals, and one steady march measured. + +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, FaceKind, Field, Fluid, Grid, Parameters, Solver, WallScheme, +}; +use std::f64::consts::PI; + +pub const RHO: f64 = 1.0; +pub const MU: f64 = 0.05; +/// The sphere's default centre (off-centre so the exact force is not zero by symmetry). +pub const C: (f64, f64, f64) = (0.6, 0.45, 0.5); +pub const R: f64 = 0.2; + +pub fn u3(x: f64, y: f64, z: f64) -> f64 { + (PI * x).sin() * (PI * y).cos() * (PI * z).cos() +} +pub fn v3(x: f64, y: f64, z: f64) -> f64 { + (PI * x).cos() * (PI * y).sin() * (PI * z).cos() +} +pub fn w3(x: f64, y: f64, z: f64) -> f64 { + -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin() +} +pub fn p3(x: f64, y: f64, z: f64) -> f64 { + (PI * x).sin() * (PI * y).sin() * (PI * z).sin() +} +/// The velocity gradient ∂u_i/∂x_j and the pressure gradient. +pub fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) { + let (sx, cx) = (PI * x).sin_cos(); + let (sy, cy) = (PI * y).sin_cos(); + let (sz, cz) = (PI * z).sin_cos(); + ( + [ + [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz], + [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz], + [ + 2.0 * PI * sx * cy * sz, + 2.0 * PI * cx * sy * sz, + -2.0 * PI * cx * cy * cz, + ], + ], + [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz], + ) +} +pub fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let (g, gp) = grads(x, y, z); + let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)]; + let lap = -3.0 * PI * PI; + let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2]; + ( + RHO * conv(0) + gp[0] - MU * lap * u[0], + RHO * conv(1) + gp[1] - MU * lap * u[1], + RHO * conv(2) + gp[2] - MU * lap * u[2], + ) +} +pub fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let u = if x <= 0.0 || x >= 1.0 { + 0.0 + } else { + u3(x, y, z) + }; + let v = if y <= 0.0 || y >= 1.0 { + 0.0 + } else { + v3(x, y, z) + }; + let w = if z <= 0.0 || z >= 1.0 { + 0.0 + } else { + w3(x, y, z) + }; + (u, v, w) +} + +/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA` +/// over the sphere by a fine Fibonacci quadrature. +pub fn exact_force_and_flux(c: (f64, f64, f64)) -> ([f64; 3], [f64; 3]) { + let n = 200_000; + let golden = PI * (3.0 - 5.0_f64.sqrt()); + let (mut f, mut m) = ([0.0; 3], [0.0; 3]); + let da = 4.0 * PI * R * R / n as f64; + for k in 0..n { + let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64; + let rr = (1.0 - zz * zz).sqrt(); + let th = golden * k as f64; + let nrm = [rr * th.cos(), rr * th.sin(), zz]; + let (x, y, z) = (c.0 + R * nrm[0], c.1 + R * nrm[1], c.2 + R * nrm[2]); + let (g, _) = grads(x, y, z); + let p = p3(x, y, z); + let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)]; + let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2]; + for i in 0..3 { + let mut t = -p * nrm[i]; + for j in 0..3 { + t += MU * (g[i][j] + g[j][i]) * nrm[j]; + } + f[i] += t * da; + m[i] += RHO * u[i] * un * da; + } + } + (f, m) +} + +pub struct Measurement { + pub l2_velocity: f64, + pub max_div: f64, + pub ghost_correction: f64, + pub force_surface: [f64; 3], + pub skipped: usize, + pub force_cv: [f64; 3], +} + +/// March the manufactured solution with the sphere at `c` to steady state on grid `n`. +pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement { + let h = 1.0 / n as f64; + let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h); + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + convection_scheme: ConvectionScheme::Upwind, + wall_scheme: scheme, + ..Parameters::default() + }, + ); + solver.set_momentum_source(|x, y, z, _t| source3(x, y, z)); + solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z)); + solver.set_body( + Body::sphere(move |_t| c, R) + .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))), + ); + let g = Grid::cubic(n, n, n, h); + let mut f = Field::new(g); + solver.initialize(&mut f); + let mut last = solver.advance(&mut f, dt); + for _ in 0..200_000 { + let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone()); + last = solver.advance(&mut f, dt); + let mut change = 0.0_f64; + for (a, b) in + f.u.iter() + .zip(&bu) + .chain(f.v.iter().zip(&bv)) + .chain(f.w.iter().zip(&bw)) + { + change = change.max((a - b).abs()); + } + if change / dt < 1e-6 { + break; + } + } + let mask = solver.mask().expect("mask"); + let (mut sq, mut vol) = (0.0, 0.0); + let dv = h * h * h; + for k in 0..n { + for j in 0..n { + for i in 1..n { + if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid { + let e = f.u[g.uface(k, j, i)] + - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h); + sq += e * e * dv; + vol += dv; + } + } + } + for j in 1..n { + for i in 0..n { + if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid { + let e = f.v[g.vface(k, j, i)] + - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h); + sq += e * e * dv; + vol += dv; + } + } + } + } + for k in 1..n { + for j in 0..n { + for i in 0..n { + if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid { + let e = f.w[g.wface(k, j, i)] + - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h); + sq += e * e * dv; + vol += dv; + } + } + } + } + let body = solver.body().expect("body"); + let t = solver.time(); + // The apertured divergence per unit volume, the porous surface's flux + // through the wall included (the plain divergence on the binary wall). + let mut max_div = 0.0_f64; + let mut at_vol = 1.0; + let mut sum_flux = 0.0; + let (wall_fluxes, _) = mask.wall_flux_table(body, t); + for k in 0..n { + for j in 0..n { + for i in 0..n { + let idx = g.cell(k, j, i); + if mask.is_fluid_cell(idx) { + let flux = (mask.a_u(g.uface(k, j, i + 1)) * f.u[g.uface(k, j, i + 1)] + - mask.a_u(g.uface(k, j, i)) * f.u[g.uface(k, j, i)]) + * h + * h + + (mask.a_v(g.vface(k, j + 1, i)) * f.v[g.vface(k, j + 1, i)] + - mask.a_v(g.vface(k, j, i)) * f.v[g.vface(k, j, i)]) + * h + * h + + (mask.a_w(g.wface(k + 1, j, i)) * f.w[g.wface(k + 1, j, i)] + - mask.a_w(g.wface(k, j, i)) * f.w[g.wface(k, j, i)]) + * h + * h + + wall_fluxes[idx]; + sum_flux += flux.abs(); + if (flux / (h * h * h)).abs() > max_div { + max_div = (flux / (h * h * h)).abs(); + at_vol = mask.vol(idx); + } + } + } + } + } + println!( + " [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}", + last.final_residual + ); + let surface = match scheme { + WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h), + WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce { + f: mask.cut_wall_force(body, &f, MU, t).expect("cut wall"), + samples: 0, + skipped: 0, + }, + }; + let (i0, i1) = (n / 8, n - n / 8); + let src = |x: f64, y: f64, z: f64| source3(x, y, z); + let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1)); + Measurement { + l2_velocity: (sq / vol).sqrt(), + max_div, + ghost_correction: solver.ghost_correction().abs(), + force_surface: surface.f, + skipped: surface.skipped, + force_cv, + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_lipschitz.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_lipschitz.rs new file mode 100644 index 0000000..39532c5 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_lipschitz.rs @@ -0,0 +1,100 @@ +//! embedded3 item 11b: the wall's smoothness in the interface position. +//! The manufactured sphere is marched to steady state at `M + 1` centres +//! spaced `h/M` apart across one cell along x; at each the load error +//! `E(δ) = F(δ) − F_exact(δ)` (the exact force moves with the sphere and +//! is subtracted) is measured on the scheme's route. The largest jump of +//! `E` between neighbouring positions, relative to the load, and the +//! Lipschitz quotient `|ΔE| / (Δδ |F|)` are reported for both walls. +//! Registered gate (`docs/embedded3_campaign.md` item 11): the cut wall's +//! largest neighbouring jump < 1 % of the load with a bounded quotient. +//! +//! Default run: 8 positions at n = 24 (about two minutes on the host); +//! the gated `#[ignore]` variant sweeps 40. + +mod embedded3_sphere; + +use embedded3_sphere::{C, exact_force_and_flux, measure}; +use rtx_cfd::solvers::incompressible::embedded3::WallScheme; + +fn norm(a: [f64; 3]) -> f64 { + (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt() +} + +struct Sweep { + /// Largest neighbouring jump of the load error relative to the load. + max_jump: f64, + /// Largest Lipschitz quotient `|ΔE| / (Δδ |F|)` (per unit length). + max_quotient: f64, +} + +fn sweep(n: usize, positions: usize, scheme: WallScheme) -> Sweep { + let h = 1.0 / n as f64; + let step = h / positions as f64; + let mut errors: Vec<[f64; 3]> = Vec::new(); + let mut scale = 0.0_f64; + for m in 0..=positions { + let c = (C.0 + m as f64 * step, C.1, C.2); + let (fe, _) = exact_force_and_flux(c); + let r = measure(n, scheme, c); + let e = [ + r.force_surface[0] - fe[0], + r.force_surface[1] - fe[1], + r.force_surface[2] - fe[2], + ]; + scale = scale.max(norm(fe)); + println!( + " {scheme:?} δ = {:.4} h: F {:.5?} exact {:.5?} error {:.3e} (rel {:.3e})", + m as f64 / positions as f64, + r.force_surface, + fe, + norm(e), + norm(e) / norm(fe) + ); + errors.push(e); + } + let mut max_jump = 0.0_f64; + for w in errors.windows(2) { + let d = norm([w[1][0] - w[0][0], w[1][1] - w[0][1], w[1][2] - w[0][2]]); + max_jump = max_jump.max(d / scale); + } + let max_quotient = max_jump / step; + println!( + " {scheme:?}: largest neighbouring jump {:.3e} of the load (spacing {:.3e} = h/{positions}); Lipschitz quotient {:.3e} per unit length", + max_jump, step, max_quotient + ); + Sweep { + max_jump, + max_quotient, + } +} + +#[test] +fn sphere_load_across_one_cell() { + let ghost = sweep(24, 8, WallScheme::GhostBinary); + let cut = sweep(24, 8, WallScheme::CutCell); + println!( + " jumps: ghost {:.3e}, cut {:.3e} ({:.1}x smaller); quotients: ghost {:.3e}, cut {:.3e}", + ghost.max_jump, + cut.max_jump, + ghost.max_jump / cut.max_jump.max(1e-300), + ghost.max_quotient, + cut.max_quotient + ); + assert!(ghost.max_jump.is_finite() && cut.max_jump.is_finite()); +} + +#[test] +#[ignore = "item 11's gated sweep (40 positions, both walls; tens of minutes on the host)"] +fn sphere_load_lipschitz_gate() { + let ghost = sweep(24, 40, WallScheme::GhostBinary); + let cut = sweep(24, 40, WallScheme::CutCell); + println!( + " GATE: cut largest jump {:.3e} of the load (ghost {:.3e}); cut quotient {:.3e} (ghost {:.3e})", + cut.max_jump, ghost.max_jump, cut.max_quotient, ghost.max_quotient + ); + assert!( + cut.max_jump < 0.01, + "cut-cell jump {:.3e} of the load", + cut.max_jump + ); +}