diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/field.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/field.rs new file mode 100644 index 0000000..63b639e --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/field.rs @@ -0,0 +1,98 @@ +//! The 3D staggered field: u on `(nx + 1) × ny × nz` faces, v on +//! `nx × (ny + 1) × nz`, w on `nx × ny × (nz + 1)`, p on the cells; flat +//! storage (see [`super::Grid`] for the index conventions). + +use super::Grid; + +#[derive(Debug, Clone)] +pub struct Field { + pub grid: Grid, + pub u: Vec, + pub v: Vec, + pub w: Vec, + pub p: Vec, + pub u_old: Vec, + pub v_old: Vec, + pub w_old: Vec, + pub u_star: Vec, + pub v_star: Vec, + pub w_star: Vec, + pub p_prime: Vec, + /// The continuity source of the projection. + pub sp: Vec, +} + +impl Field { + #[must_use] + pub fn new(grid: Grid) -> Self { + let nu = (grid.nx + 1) * grid.ny * grid.nz; + let nv = grid.nx * (grid.ny + 1) * grid.nz; + let nw = grid.nx * grid.ny * (grid.nz + 1); + let nc = grid.cells(); + Self { + grid, + u: vec![0.0; nu], + v: vec![0.0; nv], + w: vec![0.0; nw], + p: vec![0.0; nc], + u_old: vec![0.0; nu], + v_old: vec![0.0; nv], + w_old: vec![0.0; nw], + u_star: vec![0.0; nu], + v_star: vec![0.0; nv], + w_star: vec![0.0; nw], + p_prime: vec![0.0; nc], + sp: vec![0.0; nc], + } + } + + pub fn update_old_values(&mut self) { + self.u_old.copy_from_slice(&self.u); + self.v_old.copy_from_slice(&self.v); + self.w_old.copy_from_slice(&self.w); + } + + pub fn copy_to_starred(&mut self) { + self.u_star.copy_from_slice(&self.u); + self.v_star.copy_from_slice(&self.v); + self.w_star.copy_from_slice(&self.w); + } + + /// `max |∇·u|` over the cells (per unit volume). + #[must_use] + pub fn max_divergence(&self) -> f64 { + let g = self.grid; + let mut worst = 0.0_f64; + for k in 0..g.nz { + for j in 0..g.ny { + for i in 0..g.nx { + let div = (self.u[g.uface(k, j, i + 1)] - self.u[g.uface(k, j, i)]) / g.dx + + (self.v[g.vface(k, j + 1, i)] - self.v[g.vface(k, j, i)]) / g.dy + + (self.w[g.wface(k + 1, j, i)] - self.w[g.wface(k, j, i)]) / g.dz; + worst = worst.max(div.abs()); + } + } + } + worst + } + + /// Kinetic energy `½ Σ (u² + v² + w²) dV` over the cells (face values + /// averaged to the cell). + #[must_use] + pub fn kinetic_energy(&self, rho: f64) -> f64 { + let g = self.grid; + let dv = g.dx * g.dy * g.dz; + let mut e = 0.0; + for k in 0..g.nz { + for j in 0..g.ny { + for i in 0..g.nx { + let uc = 0.5 * (self.u[g.uface(k, j, i)] + self.u[g.uface(k, j, i + 1)]); + let vc = 0.5 * (self.v[g.vface(k, j, i)] + self.v[g.vface(k, j + 1, i)]); + let wc = 0.5 * (self.w[g.wface(k, j, i)] + self.w[g.wface(k + 1, j, i)]); + e += 0.5 * rho * (uc * uc + vc * vc + wc * wc) * dv; + } + } + } + e + } +} 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 34002e2..2a6036e 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -6,7 +6,11 @@ //! Layout: cells `(k, j, i)` row-major, `cell = (k·ny + j)·nx + i`; //! u faces on `(nx + 1)·ny·nz`, v on `nx·(ny + 1)·nz`, w on `nx·ny·(nz + 1)`. +pub mod field; pub mod grid; pub mod poisson; +pub mod step; +pub use field::Field; pub use grid::Grid; +pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult}; 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 new file mode 100644 index 0000000..6cd36be --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -0,0 +1,363 @@ +//! The host reference of the PISO step: the 2D embedded solver's predictor +//! and projection transcribed expression for expression, the z terms +//! APPENDED after the 2D expression so that at `nz = 1` (z sides slip, +//! `dz = 1`) every number is the 2D solver's. The fluid predicates are the +//! wall's hooks (item 9). + +mod predictor; +mod projection; + +use super::Grid; +use super::field::Field; +use super::poisson::{PcgCache, Problem, solve_pcg_cached}; +use crate::solvers::incompressible::poisson::{ + MgPrecision, MgSmoother, MultigridParameters, PoissonSolution, +}; +use crate::solvers::incompressible::simple::ConvectionScheme; + +/// Boundary type of one domain side (the 2D `SideBoundary` semantics, plus +/// `Periodic`, allowed on the z pair only in Stage 1). +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum Side { + #[default] + Velocity, + SlipWall, + PressureOutlet, + Periodic, +} + +/// The six sides: `x0` at `x = 0`, `x1` at `x = nx dx`, and so on. +#[derive(Debug, Clone, Copy, Default)] +pub struct Boundaries { + pub x0: Side, + pub x1: Side, + pub y0: Side, + pub y1: Side, + pub z0: Side, + pub z1: Side, +} + +impl Boundaries { + fn any_outlet(self) -> bool { + [self.x0, self.x1, self.y0, self.y1, self.z0, self.z1].contains(&Side::PressureOutlet) + } + + pub(super) fn periodic_z(self) -> bool { + self.z0 == Side::Periodic + } +} + +/// The fluid and the reference scales (the 2D `CfdConfig` fields used). +#[derive(Debug, Clone, Copy)] +pub struct Fluid { + pub density: f64, + pub viscosity: f64, + pub reference_velocity: f64, + pub reference_length: f64, +} + +#[derive(Debug, Clone)] +pub struct Parameters { + pub corrector_steps: usize, + /// Tolerance on the normalised mass imbalance after correction. + pub tolerance: f64, + pub boundaries: Boundaries, + pub poisson_smoother: MgSmoother, + pub poisson_precision: MgPrecision, + pub convection_scheme: ConvectionScheme, + /// Relative part of the pressure solve's inner stop (the 2D 1e-2). + pub inner_stop_factor: f64, +} + +impl Default for Parameters { + fn default() -> Self { + Self { + corrector_steps: 2, + tolerance: 1e-6, + boundaries: Boundaries::default(), + poisson_smoother: MgSmoother::Lexicographic, + poisson_precision: MgPrecision::F64, + convection_scheme: ConvectionScheme::Upwind, + inner_stop_factor: 1e-2, + } + } +} + +/// One step's outcome. +#[derive(Debug, Clone, Copy)] +pub struct StepResult { + pub converged: bool, + pub corrector_steps_performed: usize, + pub final_residual: f64, + /// CG iterations summed over the step's projections. + pub poisson_iterations: usize, +} + +type Vec3Fn = Box (f64, f64, f64) + Send + Sync>; + +pub struct Solver { + pub fluid: Fluid, + pub params: Parameters, + pub(super) momentum_source: Option, + boundary_velocity: Option, + pcg_cache: PcgCache, + time: f64, + initialized: bool, + /// `(setup ns, iterate ns, solves, CG iterations)` summed. + poisson_profile: (u64, u64, u64, u64), +} + +impl Solver { + #[must_use] + pub fn new(fluid: Fluid, params: Parameters) -> Self { + let b = params.boundaries; + assert_eq!( + b.z0 == Side::Periodic, + b.z1 == Side::Periodic, + "periodic z needs both z sides periodic" + ); + for side in [b.x0, b.x1, b.y0, b.y1] { + assert_ne!(side, Side::Periodic, "Stage 1: periodic only in z"); + } + Self { + fluid, + params, + momentum_source: None, + boundary_velocity: None, + pcg_cache: PcgCache::default(), + time: 0.0, + initialized: false, + poisson_profile: (0, 0, 0, 0), + } + } + + pub fn set_momentum_source(&mut self, f: F) + where + F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static, + { + self.momentum_source = Some(Box::new(f)); + } + + pub fn set_boundary_velocity(&mut self, f: F) + where + F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static, + { + self.boundary_velocity = Some(Box::new(f)); + } + + #[must_use] + pub fn time(&self) -> f64 { + self.time + } + + pub fn set_time(&mut self, t: f64) { + self.time = t; + } + + #[must_use] + pub fn poisson_profile(&self) -> (u64, u64, u64, u64) { + self.poisson_profile + } + + pub(super) fn boundary(&self, x: f64, y: f64, z: f64, t: f64) -> (f64, f64, f64) { + self.boundary_velocity + .as_ref() + .map_or((0.0, 0.0, 0.0), |f| f(x, y, z, t)) + } + + // The fluid predicates: everything is fluid until the wall arrives. + #[inline] + fn u_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool { + true + } + #[inline] + fn v_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool { + true + } + #[inline] + fn w_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool { + true + } + #[inline] + fn cell_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool { + true + } + + pub(super) fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 { + if face_velocity >= 0.0 { + upstream + } else { + downstream + } + } + + /// Stamp the prescribed normal velocities on the six sides at `t`. + pub fn apply_boundary_normals(&self, field: &mut Field, t: f64) { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let b = self.params.boundaries; + let outlet = Side::PressureOutlet; + for k in 0..nz { + let z = (k as f64 + 0.5) * dz; + for j in 0..ny { + let y = (j as f64 + 0.5) * dy; + if b.x0 != outlet { + field.u[g.uface(k, j, 0)] = self.boundary(0.0, y, z, t).0; + } + if b.x1 != outlet { + field.u[g.uface(k, j, nx)] = self.boundary(nx as f64 * dx, y, z, t).0; + } + } + for i in 0..nx { + let x = (i as f64 + 0.5) * dx; + if b.y0 != outlet { + field.v[g.vface(k, 0, i)] = self.boundary(x, 0.0, z, t).1; + } + if b.y1 != outlet { + field.v[g.vface(k, ny, i)] = self.boundary(x, ny as f64 * dy, z, t).1; + } + } + } + if !b.periodic_z() { + for j in 0..ny { + let y = (j as f64 + 0.5) * dy; + for i in 0..nx { + let x = (i as f64 + 0.5) * dx; + if b.z0 != outlet { + field.w[g.wface(0, j, i)] = self.boundary(x, y, 0.0, t).2; + } + if b.z1 != outlet { + field.w[g.wface(nz, j, i)] = self.boundary(x, y, nz as f64 * dz, t).2; + } + } + } + } + } + + /// The explicit predictor on the fluid faces; outlet faces zero-gradient. + fn momentum_predictor(&self, field: &mut Field, dt: f64, t_old: f64) { + let g = field.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let b = self.params.boundaries; + let periodic = b.periodic_z(); + for k in 0..nz { + for j in 0..ny { + for i in 1..nx { + if !self.u_is_fluid(k, j, i) { + continue; + } + let rhs = self.u_rhs(field, k, j, i, t_old); + let f = g.uface(k, j, i); + field.u[f] = field.u_old[f] + dt * rhs; + } + } + } + for k in 0..nz { + for j in 1..ny { + for i in 0..nx { + if !self.v_is_fluid(k, j, i) { + continue; + } + let rhs = self.v_rhs(field, k, j, i, t_old); + let f = g.vface(k, j, i); + field.v[f] = field.v_old[f] + dt * rhs; + } + } + } + let k_range = if periodic { 0..nz } else { 1..nz }; + for k in k_range { + for j in 0..ny { + for i in 0..nx { + if !self.w_is_fluid(k, j, i) { + continue; + } + let rhs = self.w_rhs(field, k, j, i, t_old); + let f = g.wface(k, j, i); + field.w[f] = field.w_old[f] + dt * rhs; + } + } + } + if periodic { + for j in 0..ny { + for i in 0..nx { + field.w[g.wface(nz, j, i)] = field.w[g.wface(0, j, i)]; + } + } + } + let outlet = Side::PressureOutlet; + for k in 0..nz { + for j in 0..ny { + if b.x0 == outlet { + field.u[g.uface(k, j, 0)] = field.u[g.uface(k, j, 1)]; + } + if b.x1 == outlet { + field.u[g.uface(k, j, nx)] = field.u[g.uface(k, j, nx - 1)]; + } + } + for i in 0..nx { + if b.y0 == outlet { + field.v[g.vface(k, 0, i)] = field.v[g.vface(k, 1, i)]; + } + if b.y1 == outlet { + field.v[g.vface(k, ny, i)] = field.v[g.vface(k, ny - 1, i)]; + } + } + } + for j in 0..ny { + for i in 0..nx { + if b.z0 == outlet { + field.w[g.wface(0, j, i)] = field.w[g.wface(1, j, i)]; + } + if b.z1 == outlet { + field.w[g.wface(nz, j, i)] = field.w[g.wface(nz - 1, j, i)]; + } + } + } + field.copy_to_starred(); + } + + /// Stamp the `t = time` boundary data (lazily called by the first step). + pub fn initialize(&mut self, field: &mut Field) { + let t = self.time; + self.apply_boundary_normals(field, t); + self.initialized = true; + } + + /// One step of `dt`: predictor, correctors, clock. + pub fn advance(&mut self, field: &mut Field, dt: f64) -> StepResult { + assert!( + dt > 0.0 && dt.is_finite(), + "time step must be positive and finite, got {dt}" + ); + if !self.initialized { + self.initialize(field); + } + let t_old = self.time; + let t_new = t_old + dt; + field.update_old_values(); + self.momentum_predictor(field, dt, t_old); + self.apply_boundary_normals(field, t_new); + field.copy_to_starred(); + let mut total = 0; + let mut final_residual = f64::INFINITY; + let mut poisson_iterations = 0; + for corrector in 0..self.params.corrector_steps.max(1) { + let sol = self.solve_correction(field, dt, corrector == 0); + poisson_iterations += sol.iterations; + let mass_residual = self.apply_correction(field, dt); + final_residual = mass_residual; + total += 1; + if mass_residual < self.params.tolerance { + break; + } + field.copy_to_starred(); + } + self.time = t_new; + StepResult { + converged: final_residual < self.params.tolerance, + corrector_steps_performed: total, + final_residual, + poisson_iterations, + } + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/predictor.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/predictor.rs new file mode 100644 index 0000000..0ab87e5 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/predictor.rs @@ -0,0 +1,679 @@ +//! The three explicit predictors of the 3D PISO step (`Solver`): the +//! 2D `u_rhs`/`v_rhs` expression for expression with the z terms appended, +//! and `w_rhs` as the v pattern turned along z. Split from `piso_host.rs` +//! for the file-size rule; `impl Solver` continues here. + +use super::{Side, Solver}; +use crate::solvers::incompressible::embedded3::field::Field; +use crate::solvers::incompressible::simple::ConvectionScheme; + +impl Solver { + /// Plane above `k` (wrapping when periodic). + #[inline] + fn k_up(&self, k: usize, nz: usize) -> Option { + if k + 1 < nz { + Some(k + 1) + } else if self.params.boundaries.periodic_z() { + Some(0) + } else { + None + } + } + + #[inline] + fn k_down(&self, k: usize, nz: usize) -> Option { + if k > 0 { + Some(k - 1) + } else if self.params.boundaries.periodic_z() { + Some(nz - 1) + } else { + None + } + } + + /// The predictor's right-hand side on the u face `(k, j, i)`, `i = 1..nx`: + /// the 2D `u_rhs` expression for expression, then `− conv_z + diff_z`. + #[allow(clippy::too_many_lines)] + pub(crate) fn u_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let rho = self.fluid.density; + let nu = self.fluid.viscosity / rho; + let b = self.params.boundaries; + let velocity = Side::Velocity; + let uo = &field.u_old; + let vo = &field.v_old; + let wo = &field.w_old; + let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii); + let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii); + let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii); + let zc = (k as f64 + 0.5) * dz; + let u_p = uo[uf(k, j, i)]; + + let ue_face = 0.5 * (uo[uf(k, j, i)] + uo[uf(k, j, i + 1)]); + let uw_face = 0.5 * (uo[uf(k, j, i - 1)] + uo[uf(k, j, i)]); + + let south_is_wall = j == 0; + let north_is_wall = j + 1 == ny; + + let vn_face = 0.5 * (vo[vf(k, j + 1, i - 1)] + vo[vf(k, j + 1, i)]); + let vs_face = 0.5 * (vo[vf(k, j, i - 1)] + vo[vf(k, j, i)]); + + let beyond_north = if b.y1 == velocity { + self.boundary(i as f64 * dx, ny as f64 * dy, zc, t_old).0 + } else { + u_p + }; + let beyond_south = if b.y0 == velocity { + self.boundary(i as f64 * dx, 0.0, zc, t_old).0 + } else { + u_p + }; + + let conv_x = (ue_face * Self::upwind(ue_face, uo[uf(k, j, i)], uo[uf(k, j, i + 1)]) + - uw_face * Self::upwind(uw_face, uo[uf(k, j, i - 1)], uo[uf(k, j, i)])) + / dx; + let conv_y = (vn_face + * if north_is_wall { + Self::upwind(vn_face, u_p, beyond_north) + } else { + Self::upwind(vn_face, uo[uf(k, j, i)], uo[uf(k, j + 1, i)]) + } + - vs_face + * if south_is_wall { + Self::upwind(vs_face, beyond_south, u_p) + } else { + Self::upwind(vs_face, uo[uf(k, j - 1, i)], uo[uf(k, j, i)]) + }) + / dy; + + let scheme = self.params.convection_scheme; + let mut conv_x = conv_x; + let mut conv_y = conv_y; + if scheme != ConvectionScheme::Upwind { + let delta_e = if ue_face >= 0.0 { + scheme.face_correction( + Some(uo[uf(k, j, i - 1)]), + uo[uf(k, j, i)], + uo[uf(k, j, i + 1)], + ) + } else { + let far = (i + 2 <= nx).then(|| uo[uf(k, j, i + 2)]); + scheme.face_correction(far, uo[uf(k, j, i + 1)], uo[uf(k, j, i)]) + }; + let delta_w = if uw_face >= 0.0 { + let far = (i >= 2).then(|| uo[uf(k, j, i - 2)]); + scheme.face_correction(far, uo[uf(k, j, i - 1)], uo[uf(k, j, i)]) + } else { + scheme.face_correction( + Some(uo[uf(k, j, i + 1)]), + uo[uf(k, j, i)], + uo[uf(k, j, i - 1)], + ) + }; + let delta_n = if north_is_wall { + 0.0 + } else if vn_face >= 0.0 { + let far = (j >= 1).then(|| uo[uf(k, j - 1, i)]); + scheme.face_correction(far, uo[uf(k, j, i)], uo[uf(k, j + 1, i)]) + } else { + let far = (j + 2 < ny).then(|| uo[uf(k, j + 2, i)]); + scheme.face_correction(far, uo[uf(k, j + 1, i)], uo[uf(k, j, i)]) + }; + let delta_s = if south_is_wall { + 0.0 + } else if vs_face >= 0.0 { + let far = (j >= 2).then(|| uo[uf(k, j - 2, i)]); + scheme.face_correction(far, uo[uf(k, j - 1, i)], uo[uf(k, j, i)]) + } else { + let far = (j + 1 < ny).then(|| uo[uf(k, j + 1, i)]); + scheme.face_correction(far, uo[uf(k, j, i)], uo[uf(k, j - 1, i)]) + }; + conv_x += (ue_face * delta_e - uw_face * delta_w) / dx; + conv_y += (vn_face * delta_n - vs_face * delta_s) / dy; + } + + let diff_x = nu * (uo[uf(k, j, i + 1)] - 2.0 * u_p + uo[uf(k, j, i - 1)]) / (dx * dx); + + let flux_north = if north_is_wall { + if b.y1 == velocity { + let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, zc, t_old).0; + nu * (u_wall - u_p) / (0.5 * dy) + } else { + 0.0 + } + } else { + nu * (uo[uf(k, j + 1, i)] - u_p) / dy + }; + let flux_south = if south_is_wall { + if b.y0 == velocity { + let u_wall = self.boundary(i as f64 * dx, 0.0, zc, t_old).0; + nu * (u_p - u_wall) / (0.5 * dy) + } else { + 0.0 + } + } else { + nu * (u_p - uo[uf(k, j - 1, i)]) / dy + }; + let diff_y = (flux_north - flux_south) / dy; + + let pressure_gradient = + -(field.p[g.cell(k, j, i)] - field.p[g.cell(k, j, i - 1)]) / (rho * dx); + + let body_force = self.momentum_source.as_ref().map_or(0.0, |f| { + f(i as f64 * dx, (j as f64 + 0.5) * dy, zc, t_old).0 / rho + }); + + let rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force; + + // --- the z terms, the y pattern turned along k --- + let ku = self.k_up(k, nz); + let kd = self.k_down(k, nz); + let top_is_wall = ku.is_none(); + let bottom_is_wall = kd.is_none(); + // The w faces above/below the u face: on top of the cells west and + // east of it (face k + 1 of cell k is face index k + 1; periodic: + // the face at k = nz equals the face at 0). + let wt_face = 0.5 * (wo[wf(k + 1, j, i - 1)] + wo[wf(k + 1, j, i)]); + let wb_face = 0.5 * (wo[wf(k, j, i - 1)] + wo[wf(k, j, i)]); + let beyond_top = if b.z1 == velocity { + self.boundary(i as f64 * dx, (j as f64 + 0.5) * dy, nz as f64 * dz, t_old) + .0 + } else { + u_p + }; + let beyond_bottom = if b.z0 == velocity { + self.boundary(i as f64 * dx, (j as f64 + 0.5) * dy, 0.0, t_old) + .0 + } else { + u_p + }; + let u_up = ku.map(|kk| uo[uf(kk, j, i)]); + let u_dn = kd.map(|kk| uo[uf(kk, j, i)]); + let mut conv_z = (wt_face + * match u_up { + Some(un) => Self::upwind(wt_face, u_p, un), + None => Self::upwind(wt_face, u_p, beyond_top), + } + - wb_face + * match u_dn { + Some(ud) => Self::upwind(wb_face, ud, u_p), + None => Self::upwind(wb_face, beyond_bottom, u_p), + }) + / dz; + if scheme != ConvectionScheme::Upwind { + let far_up2 = ku + .and_then(|kk| self.k_up(kk, nz)) + .map(|kk| uo[uf(kk, j, i)]); + let far_dn2 = kd + .and_then(|kk| self.k_down(kk, nz)) + .map(|kk| uo[uf(kk, j, i)]); + let delta_t = if top_is_wall { + 0.0 + } else if wt_face >= 0.0 { + scheme.face_correction(u_dn, u_p, u_up.unwrap_or(u_p)) + } else { + scheme.face_correction(far_up2, u_up.unwrap_or(u_p), u_p) + }; + let delta_b = if bottom_is_wall { + 0.0 + } else if wb_face >= 0.0 { + scheme.face_correction(far_dn2, u_dn.unwrap_or(u_p), u_p) + } else { + scheme.face_correction(u_up, u_p, u_dn.unwrap_or(u_p)) + }; + conv_z += (wt_face * delta_t - wb_face * delta_b) / dz; + } + let flux_top = match u_up { + Some(un) => nu * (un - u_p) / dz, + None => { + if b.z1 == velocity { + nu * (beyond_top - u_p) / (0.5 * dz) + } else { + 0.0 + } + } + }; + let flux_bottom = match u_dn { + Some(ud) => nu * (u_p - ud) / dz, + None => { + if b.z0 == velocity { + nu * (u_p - beyond_bottom) / (0.5 * dz) + } else { + 0.0 + } + } + }; + let diff_z = (flux_top - flux_bottom) / dz; + + rhs_2d - conv_z + diff_z + } + + /// The v face `(k, j, i)`, `j = 1..ny`: the 2D `v_rhs` then the z terms. + #[allow(clippy::too_many_lines)] + pub(crate) fn v_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let rho = self.fluid.density; + let nu = self.fluid.viscosity / rho; + let b = self.params.boundaries; + let velocity = Side::Velocity; + let uo = &field.u_old; + let vo = &field.v_old; + let wo = &field.w_old; + let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii); + let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii); + let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii); + let zc = (k as f64 + 0.5) * dz; + let v_p = vo[vf(k, j, i)]; + + let vn_face = 0.5 * (vo[vf(k, j, i)] + vo[vf(k, j + 1, i)]); + let vs_face = 0.5 * (vo[vf(k, j - 1, i)] + vo[vf(k, j, i)]); + + let west_is_wall = i == 0; + let east_is_wall = i + 1 == nx; + + let ue_face = 0.5 * (uo[uf(k, j - 1, i + 1)] + uo[uf(k, j, i + 1)]); + let uw_face = 0.5 * (uo[uf(k, j - 1, i)] + uo[uf(k, j, i)]); + let beyond_east = if b.x1 == velocity { + self.boundary(nx as f64 * dx, j as f64 * dy, zc, t_old).1 + } else { + v_p + }; + let beyond_west = if b.x0 == velocity { + self.boundary(0.0, j as f64 * dy, zc, t_old).1 + } else { + v_p + }; + + let conv_y = (vn_face * Self::upwind(vn_face, vo[vf(k, j, i)], vo[vf(k, j + 1, i)]) + - vs_face * Self::upwind(vs_face, vo[vf(k, j - 1, i)], vo[vf(k, j, i)])) + / dy; + let conv_x = (ue_face + * if east_is_wall { + Self::upwind(ue_face, v_p, beyond_east) + } else { + Self::upwind(ue_face, vo[vf(k, j, i)], vo[vf(k, j, i + 1)]) + } + - uw_face + * if west_is_wall { + Self::upwind(uw_face, beyond_west, v_p) + } else { + Self::upwind(uw_face, vo[vf(k, j, i - 1)], vo[vf(k, j, i)]) + }) + / dx; + + let scheme = self.params.convection_scheme; + let mut conv_x = conv_x; + let mut conv_y = conv_y; + if scheme != ConvectionScheme::Upwind { + let delta_n = if vn_face >= 0.0 { + scheme.face_correction( + Some(vo[vf(k, j - 1, i)]), + vo[vf(k, j, i)], + vo[vf(k, j + 1, i)], + ) + } else { + let far = (j + 2 <= ny).then(|| vo[vf(k, j + 2, i)]); + scheme.face_correction(far, vo[vf(k, j + 1, i)], vo[vf(k, j, i)]) + }; + let delta_s = if vs_face >= 0.0 { + let far = (j >= 2).then(|| vo[vf(k, j - 2, i)]); + scheme.face_correction(far, vo[vf(k, j - 1, i)], vo[vf(k, j, i)]) + } else { + scheme.face_correction( + Some(vo[vf(k, j + 1, i)]), + vo[vf(k, j, i)], + vo[vf(k, j - 1, i)], + ) + }; + let delta_e = if east_is_wall { + 0.0 + } else if ue_face >= 0.0 { + let far = (i >= 1).then(|| vo[vf(k, j, i - 1)]); + scheme.face_correction(far, vo[vf(k, j, i)], vo[vf(k, j, i + 1)]) + } else { + let far = (i + 2 < nx).then(|| vo[vf(k, j, i + 2)]); + scheme.face_correction(far, vo[vf(k, j, i + 1)], vo[vf(k, j, i)]) + }; + let delta_w = if west_is_wall { + 0.0 + } else if uw_face >= 0.0 { + let far = (i >= 2).then(|| vo[vf(k, j, i - 2)]); + scheme.face_correction(far, vo[vf(k, j, i - 1)], vo[vf(k, j, i)]) + } else { + let far = (i + 1 < nx).then(|| vo[vf(k, j, i + 1)]); + scheme.face_correction(far, vo[vf(k, j, i)], vo[vf(k, j, i - 1)]) + }; + conv_y += (vn_face * delta_n - vs_face * delta_s) / dy; + conv_x += (ue_face * delta_e - uw_face * delta_w) / dx; + } + + let diff_y = nu * (vo[vf(k, j + 1, i)] - 2.0 * v_p + vo[vf(k, j - 1, i)]) / (dy * dy); + + let flux_east = if east_is_wall { + if b.x1 == velocity { + let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, zc, t_old).1; + nu * (v_wall - v_p) / (0.5 * dx) + } else { + 0.0 + } + } else { + nu * (vo[vf(k, j, i + 1)] - v_p) / dx + }; + let flux_west = if west_is_wall { + if b.x0 == velocity { + let v_wall = self.boundary(0.0, j as f64 * dy, zc, t_old).1; + nu * (v_p - v_wall) / (0.5 * dx) + } else { + 0.0 + } + } else { + nu * (v_p - vo[vf(k, j, i - 1)]) / dx + }; + let diff_x = (flux_east - flux_west) / dx; + + let pressure_gradient = + -(field.p[g.cell(k, j, i)] - field.p[g.cell(k, j - 1, i)]) / (rho * dy); + + let body_force = self.momentum_source.as_ref().map_or(0.0, |f| { + f((i as f64 + 0.5) * dx, j as f64 * dy, zc, t_old).1 / rho + }); + + let rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force; + + // --- z terms --- + let ku = self.k_up(k, nz); + let kd = self.k_down(k, nz); + let top_is_wall = ku.is_none(); + let bottom_is_wall = kd.is_none(); + let wt_face = 0.5 * (wo[wf(k + 1, j - 1, i)] + wo[wf(k + 1, j, i)]); + let wb_face = 0.5 * (wo[wf(k, j - 1, i)] + wo[wf(k, j, i)]); + let beyond_top = if b.z1 == velocity { + self.boundary((i as f64 + 0.5) * dx, j as f64 * dy, nz as f64 * dz, t_old) + .1 + } else { + v_p + }; + let beyond_bottom = if b.z0 == velocity { + self.boundary((i as f64 + 0.5) * dx, j as f64 * dy, 0.0, t_old) + .1 + } else { + v_p + }; + let v_up = ku.map(|kk| vo[vf(kk, j, i)]); + let v_dn = kd.map(|kk| vo[vf(kk, j, i)]); + let mut conv_z = (wt_face + * match v_up { + Some(vn) => Self::upwind(wt_face, v_p, vn), + None => Self::upwind(wt_face, v_p, beyond_top), + } + - wb_face + * match v_dn { + Some(vd) => Self::upwind(wb_face, vd, v_p), + None => Self::upwind(wb_face, beyond_bottom, v_p), + }) + / dz; + if scheme != ConvectionScheme::Upwind { + let far_up2 = ku + .and_then(|kk| self.k_up(kk, nz)) + .map(|kk| vo[vf(kk, j, i)]); + let far_dn2 = kd + .and_then(|kk| self.k_down(kk, nz)) + .map(|kk| vo[vf(kk, j, i)]); + let delta_t = if top_is_wall { + 0.0 + } else if wt_face >= 0.0 { + scheme.face_correction(v_dn, v_p, v_up.unwrap_or(v_p)) + } else { + scheme.face_correction(far_up2, v_up.unwrap_or(v_p), v_p) + }; + let delta_b = if bottom_is_wall { + 0.0 + } else if wb_face >= 0.0 { + scheme.face_correction(far_dn2, v_dn.unwrap_or(v_p), v_p) + } else { + scheme.face_correction(v_up, v_p, v_dn.unwrap_or(v_p)) + }; + conv_z += (wt_face * delta_t - wb_face * delta_b) / dz; + } + let flux_top = match v_up { + Some(vn) => nu * (vn - v_p) / dz, + None => { + if b.z1 == velocity { + nu * (beyond_top - v_p) / (0.5 * dz) + } else { + 0.0 + } + } + }; + let flux_bottom = match v_dn { + Some(vd) => nu * (v_p - vd) / dz, + None => { + if b.z0 == velocity { + nu * (v_p - beyond_bottom) / (0.5 * dz) + } else { + 0.0 + } + } + }; + let diff_z = (flux_top - flux_bottom) / dz; + + rhs_2d - conv_z + diff_z + } + + /// The w face `(k, j, i)` between cells `k − 1` (wrapping when periodic) + /// and `k`: the v pattern with z as its own direction and x, y transverse. + #[allow(clippy::too_many_lines)] + pub(crate) fn w_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let rho = self.fluid.density; + let nu = self.fluid.viscosity / rho; + let b = self.params.boundaries; + let velocity = Side::Velocity; + let periodic = b.periodic_z(); + let uo = &field.u_old; + let vo = &field.v_old; + let wo = &field.w_old; + let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii); + let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii); + let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii); + // The cells below and above this face (k = 0 only when periodic). + let k_below = if k > 0 { k - 1 } else { nz - 1 }; + let k_above = k % nz; + // Own-direction neighbours: the faces k − 1 and k + 1 (wrapping). + let w_dn_idx = if k > 0 { + wf(k - 1, j, i) + } else { + wf(nz - 1, j, i) + }; + let w_up_idx = if k + 1 == nz && periodic { + wf(0, j, i) + } else { + wf(k + 1, j, i) + }; + let zf = k as f64 * dz; + let w_p = wo[wf(k, j, i)]; + + let wt_face = 0.5 * (wo[wf(k, j, i)] + wo[w_up_idx]); + let wb_face = 0.5 * (wo[w_dn_idx] + wo[wf(k, j, i)]); + + let west_is_wall = i == 0; + let east_is_wall = i + 1 == nx; + let south_is_wall = j == 0; + let north_is_wall = j + 1 == ny; + + let ue_face = 0.5 * (uo[uf(k_below, j, i + 1)] + uo[uf(k_above, j, i + 1)]); + let uw_face = 0.5 * (uo[uf(k_below, j, i)] + uo[uf(k_above, j, i)]); + let vn_face = 0.5 * (vo[vf(k_below, j + 1, i)] + vo[vf(k_above, j + 1, i)]); + let vs_face = 0.5 * (vo[vf(k_below, j, i)] + vo[vf(k_above, j, i)]); + let beyond_east = if b.x1 == velocity { + self.boundary(nx as f64 * dx, (j as f64 + 0.5) * dy, zf, t_old) + .2 + } else { + w_p + }; + let beyond_west = if b.x0 == velocity { + self.boundary(0.0, (j as f64 + 0.5) * dy, zf, t_old).2 + } else { + w_p + }; + let beyond_north = if b.y1 == velocity { + self.boundary((i as f64 + 0.5) * dx, ny as f64 * dy, zf, t_old) + .2 + } else { + w_p + }; + let beyond_south = if b.y0 == velocity { + self.boundary((i as f64 + 0.5) * dx, 0.0, zf, t_old).2 + } else { + w_p + }; + + let conv_z = (wt_face * Self::upwind(wt_face, w_p, wo[w_up_idx]) + - wb_face * Self::upwind(wb_face, wo[w_dn_idx], w_p)) + / dz; + let conv_x = (ue_face + * if east_is_wall { + Self::upwind(ue_face, w_p, beyond_east) + } else { + Self::upwind(ue_face, w_p, wo[wf(k, j, i + 1)]) + } + - uw_face + * if west_is_wall { + Self::upwind(uw_face, beyond_west, w_p) + } else { + Self::upwind(uw_face, wo[wf(k, j, i - 1)], w_p) + }) + / dx; + let conv_y = (vn_face + * if north_is_wall { + Self::upwind(vn_face, w_p, beyond_north) + } else { + Self::upwind(vn_face, w_p, wo[wf(k, j + 1, i)]) + } + - vs_face + * if south_is_wall { + Self::upwind(vs_face, beyond_south, w_p) + } else { + Self::upwind(vs_face, wo[wf(k, j - 1, i)], w_p) + }) + / dy; + + let scheme = self.params.convection_scheme; + let mut conv_x = conv_x; + let mut conv_y = conv_y; + let mut conv_z = conv_z; + if scheme != ConvectionScheme::Upwind { + // Own direction: far nodes two faces away (wrapping when periodic). + let far_up2 = if periodic { + Some(wo[wf((k + 2) % nz, j, i)]) + } else { + (k + 2 <= nz).then(|| wo[wf(k + 2, j, i)]) + }; + let far_dn2 = if periodic { + Some(wo[wf((k + nz - 2) % nz, j, i)]) + } else { + (k >= 2).then(|| wo[wf(k - 2, j, i)]) + }; + let delta_t = if wt_face >= 0.0 { + scheme.face_correction(Some(wo[w_dn_idx]), w_p, wo[w_up_idx]) + } else { + scheme.face_correction(far_up2, wo[w_up_idx], w_p) + }; + let delta_b = if wb_face >= 0.0 { + scheme.face_correction(far_dn2, wo[w_dn_idx], w_p) + } else { + scheme.face_correction(Some(wo[w_up_idx]), w_p, wo[w_dn_idx]) + }; + let delta_e = if east_is_wall { + 0.0 + } else if ue_face >= 0.0 { + let far = (i >= 1).then(|| wo[wf(k, j, i - 1)]); + scheme.face_correction(far, w_p, wo[wf(k, j, i + 1)]) + } else { + let far = (i + 2 < nx).then(|| wo[wf(k, j, i + 2)]); + scheme.face_correction(far, wo[wf(k, j, i + 1)], w_p) + }; + let delta_w = if west_is_wall { + 0.0 + } else if uw_face >= 0.0 { + let far = (i >= 2).then(|| wo[wf(k, j, i - 2)]); + scheme.face_correction(far, wo[wf(k, j, i - 1)], w_p) + } else { + let far = (i + 1 < nx).then(|| wo[wf(k, j, i + 1)]); + scheme.face_correction(far, w_p, wo[wf(k, j, i - 1)]) + }; + let delta_n = if north_is_wall { + 0.0 + } else if vn_face >= 0.0 { + let far = (j >= 1).then(|| wo[wf(k, j - 1, i)]); + scheme.face_correction(far, w_p, wo[wf(k, j + 1, i)]) + } else { + let far = (j + 2 < ny).then(|| wo[wf(k, j + 2, i)]); + scheme.face_correction(far, wo[wf(k, j + 1, i)], w_p) + }; + let delta_s = if south_is_wall { + 0.0 + } else if vs_face >= 0.0 { + let far = (j >= 2).then(|| wo[wf(k, j - 2, i)]); + scheme.face_correction(far, wo[wf(k, j - 1, i)], w_p) + } else { + let far = (j + 1 < ny).then(|| wo[wf(k, j + 1, i)]); + scheme.face_correction(far, w_p, wo[wf(k, j - 1, i)]) + }; + conv_z += (wt_face * delta_t - wb_face * delta_b) / dz; + conv_x += (ue_face * delta_e - uw_face * delta_w) / dx; + conv_y += (vn_face * delta_n - vs_face * delta_s) / dy; + } + + let diff_z = nu * (wo[w_up_idx] - 2.0 * w_p + wo[w_dn_idx]) / (dz * dz); + let flux_east = if east_is_wall { + if b.x1 == velocity { + nu * (beyond_east - w_p) / (0.5 * dx) + } else { + 0.0 + } + } else { + nu * (wo[wf(k, j, i + 1)] - w_p) / dx + }; + let flux_west = if west_is_wall { + if b.x0 == velocity { + nu * (w_p - beyond_west) / (0.5 * dx) + } else { + 0.0 + } + } else { + nu * (w_p - wo[wf(k, j, i - 1)]) / dx + }; + let diff_x = (flux_east - flux_west) / dx; + let flux_north = if north_is_wall { + if b.y1 == velocity { + nu * (beyond_north - w_p) / (0.5 * dy) + } else { + 0.0 + } + } else { + nu * (wo[wf(k, j + 1, i)] - w_p) / dy + }; + let flux_south = if south_is_wall { + if b.y0 == velocity { + nu * (w_p - beyond_south) / (0.5 * dy) + } else { + 0.0 + } + } else { + nu * (w_p - wo[wf(k, j - 1, i)]) / dy + }; + let diff_y = (flux_north - flux_south) / dy; + + let pressure_gradient = + -(field.p[g.cell(k_above, j, i)] - field.p[g.cell(k_below, j, i)]) / (rho * dz); + let body_force = self.momentum_source.as_ref().map_or(0.0, |f| { + f((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy, zf, t_old).2 / rho + }); + + -conv_x - conv_y - conv_z + diff_x + diff_y + diff_z + pressure_gradient + body_force + } +} 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 new file mode 100644 index 0000000..f9bf472 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs @@ -0,0 +1,418 @@ +//! The projection: the pressure-correction operator, the boundary and +//! source tables, the solve and apply halves of one corrector. + +use super::{Boundaries, Side, Solver}; +use crate::solvers::incompressible::embedded3::Grid; +use crate::solvers::incompressible::embedded3::field::Field; +use crate::solvers::incompressible::embedded3::poisson::{Problem, solve_pcg_cached}; +use crate::solvers::incompressible::poisson::{MultigridParameters, PoissonSolution}; + +impl Solver { + /// The pressure-correction OPERATOR (the 2D `poisson_problem` with the z + /// faces): coefficient `dt · A / δ` across every fluid interior face, + /// zero across prescribed ones, the outlet's Dirichlet half a cell out + /// as `extra_diag`; the right-hand side left zero. + pub(crate) fn poisson_operator(&self, g: Grid, dt: f64) -> Problem { + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let b = self.params.boundaries; + let outlet = Side::PressureOutlet; + let periodic = b.periodic_z(); + let ae_interior = dt * (dy * dz) / dx; + let an_interior = dt * (dx * dz) / dy; + let at_interior = dt * (dx * dy) / dz; + let ae_outlet = dt * (dy * dz) / (0.5 * dx); + let an_outlet = dt * (dx * dz) / (0.5 * dy); + let at_outlet = dt * (dx * dy) / (0.5 * dz); + let mut problem = Problem::new(nx, ny, nz); + problem.periodic_z = periodic; + for k in 0..nz { + 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) { + problem.active[idx] = false; + continue; + } + let mut extra = 0.0; + if i + 1 == nx { + if b.x1 == outlet { + extra += ae_outlet; + } + } else if self.u_is_fluid(k, j, i + 1) { + problem.ae[idx] = ae_interior; + } + if i == 0 { + if b.x0 == outlet { + extra += ae_outlet; + } + } else if self.u_is_fluid(k, j, i) { + problem.aw[idx] = ae_interior; + } + if j + 1 == ny { + if b.y1 == outlet { + extra += an_outlet; + } + } else if self.v_is_fluid(k, j + 1, i) { + problem.an[idx] = an_interior; + } + if j == 0 { + if b.y0 == outlet { + extra += an_outlet; + } + } else if self.v_is_fluid(k, j, i) { + problem.as_[idx] = an_interior; + } + if k + 1 == nz && !periodic { + if b.z1 == outlet { + extra += at_outlet; + } + } else if self.w_is_fluid((k + 1) % nz, j, i) { + problem.at[idx] = at_interior; + } + if k == 0 && !periodic { + if b.z0 == outlet { + extra += at_outlet; + } + } else if self.w_is_fluid(k, j, i) { + problem.ab[idx] = at_interior; + } + problem.extra_diag[idx] = extra; + } + } + } + problem + } + + /// The operator with `field.sp` as the right-hand side. + pub(crate) fn poisson_problem(&self, field: &Field, dt: f64) -> Problem { + let mut problem = self.poisson_operator(field.grid, dt); + problem.rhs.copy_from_slice(&field.sp); + problem + } + + /// The anchor cell of a pure-Neumann projection (the first fluid cell, + /// the 2D `(1, 1)` at `k = 0`), or `None` with an outlet. + pub(crate) fn anchor_cell(&self, g: Grid) -> Option { + (!self.params.boundaries.any_outlet()).then_some(g.cell(0, 1, 1)) + } + + /// The inner stop of a projection from the source scale (the 2D rule). + pub(crate) fn inner_stop(&self, g: Grid, source_scale: f64) -> f64 { + (self.params.inner_stop_factor * source_scale) + .max(0.1 * self.params.tolerance * self.reference_flux(g)) + + 1e-14 + } + + /// The boundary function on the six sides at every point a predictor or + /// the stamping reads (`[side][component]`, see the device driver): + /// x sides: u at `(k, j)`, v at `(k, j = 0..=ny)`, w at `(k = 0..=nz, j)`; + /// y sides: u at `(k, i = 0..=nx)`, v at `(k, i)`, w at `(k = 0..=nz, i)`; + /// z sides: u at `(j, i = 0..=nx)`, v at `(j = 0..=ny, i)`, w at `(j, i)`. + #[allow(clippy::type_complexity)] + pub fn boundary_tables(&self, g: Grid, t: f64) -> [[Vec; 3]; 6] { + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let xc = |i: usize| (i as f64 + 0.5) * dx; + let yc = |j: usize| (j as f64 + 0.5) * dy; + let zc = |k: usize| (k as f64 + 0.5) * dz; + let mut out: [[Vec; 3]; 6] = Default::default(); + for (s, x) in [(0usize, 0.0), (1, nx as f64 * dx)] { + let mut u = vec![0.0; ny * nz]; + let mut v = vec![0.0; (ny + 1) * nz]; + let mut w = vec![0.0; ny * (nz + 1)]; + for k in 0..nz { + for j in 0..ny { + u[k * ny + j] = self.boundary(x, yc(j), zc(k), t).0; + } + for j in 0..=ny { + v[k * (ny + 1) + j] = self.boundary(x, j as f64 * dy, zc(k), t).1; + } + } + for k in 0..=nz { + for j in 0..ny { + w[k * ny + j] = self.boundary(x, yc(j), k as f64 * dz, t).2; + } + } + out[s] = [u, v, w]; + } + for (s, y) in [(2usize, 0.0), (3, ny as f64 * dy)] { + let mut u = vec![0.0; (nx + 1) * nz]; + let mut v = vec![0.0; nx * nz]; + let mut w = vec![0.0; nx * (nz + 1)]; + for k in 0..nz { + for i in 0..=nx { + u[k * (nx + 1) + i] = self.boundary(i as f64 * dx, y, zc(k), t).0; + } + for i in 0..nx { + v[k * nx + i] = self.boundary(xc(i), y, zc(k), t).1; + } + } + for k in 0..=nz { + for i in 0..nx { + w[k * nx + i] = self.boundary(xc(i), y, k as f64 * dz, t).2; + } + } + out[s] = [u, v, w]; + } + for (s, z) in [(4usize, 0.0), (5, nz as f64 * dz)] { + let mut u = vec![0.0; (nx + 1) * ny]; + let mut v = vec![0.0; nx * (ny + 1)]; + let mut w = vec![0.0; nx * ny]; + for j in 0..ny { + for i in 0..=nx { + u[j * (nx + 1) + i] = self.boundary(i as f64 * dx, yc(j), z, t).0; + } + for i in 0..nx { + w[j * nx + i] = self.boundary(xc(i), yc(j), z, t).2; + } + } + for j in 0..=ny { + for i in 0..nx { + v[j * nx + i] = self.boundary(xc(i), j as f64 * dy, z, t).1; + } + } + out[s] = [u, v, w]; + } + out + } + + /// The momentum source at every u, v, w face at time `t`. + pub fn source_tables(&self, g: Grid, t: f64) -> (Vec, Vec, Vec) { + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let mut su = vec![0.0; (nx + 1) * ny * nz]; + let mut sv = vec![0.0; nx * (ny + 1) * nz]; + let mut sw = vec![0.0; nx * ny * (nz + 1)]; + if let Some(f) = self.momentum_source.as_ref() { + for k in 0..nz { + for j in 0..ny { + for i in 0..=nx { + su[g.uface(k, j, i)] = f( + i as f64 * dx, + (j as f64 + 0.5) * dy, + (k as f64 + 0.5) * dz, + t, + ) + .0; + } + } + } + for k in 0..nz { + for j in 0..=ny { + for i in 0..nx { + sv[g.vface(k, j, i)] = f( + (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 { + sw[g.wface(k, j, i)] = f( + (i as f64 + 0.5) * dx, + (j as f64 + 0.5) * dy, + k as f64 * dz, + t, + ) + .2; + } + } + } + } + (su, sv, sw) + } + + pub(crate) fn reference_flux(&self, g: Grid) -> f64 { + self.fluid.density + * self.fluid.reference_velocity + * self.fluid.reference_length + * (g.nz as f64 * g.dz) + } + + /// The solve half of a projection: `sp` from `u*`, then `p'`. + pub(crate) fn solve_correction( + &mut self, + field: &mut Field, + dt: f64, + warm_start: bool, + ) -> PoissonSolution { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let rho = self.fluid.density; + let mut source_scale = 0.0; + for k in 0..nz { + 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) { + field.sp[idx] = 0.0; + continue; + } + let divergence_flux = rho + * ((field.u_star[g.uface(k, j, i + 1)] - field.u_star[g.uface(k, j, i)]) + * (dy * dz) + + (field.v_star[g.vface(k, j + 1, i)] + - field.v_star[g.vface(k, j, i)]) + * (dx * dz) + + (field.w_star[g.wface(k + 1, j, i)] + - field.w_star[g.wface(k, j, i)]) + * (dx * dy)); + field.sp[idx] = -divergence_flux; + source_scale += divergence_flux.abs(); + } + } + } + let inner_stop = self.inner_stop(g, source_scale); + let problem = self.poisson_problem(field, dt); + let mut p_prime = vec![0.0; g.cells()]; + if warm_start { + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + if self.cell_is_fluid(k, j, i) { + let idx = g.cell(k, j, i); + p_prime[idx] = field.p_prime[idx]; + } + } + } + } + } + let anchor_cell = self.anchor_cell(g); + let params = MultigridParameters { + precision: self.params.poisson_precision, + smoother: self.params.poisson_smoother, + ..MultigridParameters::default() + }; + let solution = solve_pcg_cached( + &problem, + &mut p_prime, + ¶ms, + inner_stop, + anchor_cell, + &mut self.pcg_cache, + ); + let (s0, i0, c0, k0) = self.poisson_profile; + self.poisson_profile = ( + s0 + solution.setup_ns, + i0 + solution.iterate_ns, + c0 + 1, + k0 + solution.iterations as u64, + ); + field.p_prime.copy_from_slice(&p_prime); + solution + } + + /// The apply half: correct the fluid faces from `u*`, add `p'` to `p`, + /// return the normalised mass imbalance. + pub(crate) fn apply_correction(&self, field: &mut Field, dt: f64) -> f64 { + let g = field.grid; + let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); + let rho = self.fluid.density; + let b = self.params.boundaries; + let outlet = Side::PressureOutlet; + let periodic = b.periodic_z(); + let pp = &field.p_prime; + for k in 0..nz { + for j in 0..ny { + for i in 1..nx { + if self.u_is_fluid(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; + } + } + if b.x0 == outlet { + let dp_dx = (pp[g.cell(k, j, 0)] - 0.0) / (0.5 * dx); + let f = g.uface(k, j, 0); + field.u[f] = field.u_star[f] - (dt / rho) * dp_dx; + } + if b.x1 == outlet { + let dp_dx = (0.0 - pp[g.cell(k, j, nx - 1)]) / (0.5 * dx); + let f = g.uface(k, j, nx); + field.u[f] = field.u_star[f] - (dt / rho) * dp_dx; + } + } + for i in 0..nx { + for j in 1..ny { + if self.v_is_fluid(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; + } + } + if b.y0 == outlet { + let dp_dy = (pp[g.cell(k, 0, i)] - 0.0) / (0.5 * dy); + let f = g.vface(k, 0, i); + field.v[f] = field.v_star[f] - (dt / rho) * dp_dy; + } + if b.y1 == outlet { + let dp_dy = (0.0 - pp[g.cell(k, ny - 1, i)]) / (0.5 * dy); + let f = g.vface(k, ny, i); + field.v[f] = field.v_star[f] - (dt / rho) * dp_dy; + } + } + } + for j in 0..ny { + 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) { + 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); + field.w[f] = field.w_star[f] - (dt / rho) * dp_dz; + } + } + if periodic { + field.w[g.wface(nz, j, i)] = field.w[g.wface(0, j, i)]; + } + if b.z0 == outlet { + let dp_dz = (pp[g.cell(0, j, i)] - 0.0) / (0.5 * dz); + let f = g.wface(0, j, i); + field.w[f] = field.w_star[f] - (dt / rho) * dp_dz; + } + if b.z1 == outlet { + let dp_dz = (0.0 - pp[g.cell(nz - 1, j, i)]) / (0.5 * dz); + let f = g.wface(nz, j, i); + field.w[f] = field.w_star[f] - (dt / rho) * dp_dz; + } + } + } + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + if self.cell_is_fluid(k, j, i) { + let idx = g.cell(k, j, i); + field.p[idx] += pp[idx]; + } + } + } + } + let mut mass_imbalance = 0.0; + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + if !self.cell_is_fluid(k, j, i) { + continue; + } + let divergence_flux = rho + * ((field.u[g.uface(k, j, i + 1)] - field.u[g.uface(k, j, i)]) * (dy * dz) + + (field.v[g.vface(k, j + 1, i)] - field.v[g.vface(k, j, i)]) + * (dx * dz) + + (field.w[g.wface(k + 1, j, i)] - field.w[g.wface(k, j, i)]) + * (dx * dy)); + mass_imbalance += divergence_flux.abs(); + } + } + } + let reference_flux = self.reference_flux(g); + if reference_flux > 0.0 { + mass_imbalance / reference_flux + } else { + mass_imbalance + } + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_beltrami.rs b/crates/specialized/rtx-cfd/tests/embedded3_beltrami.rs new file mode 100644 index 0000000..4ebfae5 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_beltrami.rs @@ -0,0 +1,269 @@ +//! embedded3 gate 5: the Ethier–Steinman (Beltrami) exact unsteady +//! solution on the unit cube with time-dependent Dirichlet data — the +//! transient machinery with no source term. (1) The L2 velocity error at +//! `T` falls under `dt ~ h²` refinement at order ≥ 0.75; (2) the kinetic +//! energy decay follows the closed form within the discretisation error. + +use rtx_cfd::solvers::incompressible::embedded3::{Field, Fluid, Grid, Parameters, Solver}; +use std::f64::consts::PI; + +const RHO: f64 = 1.0; +const NU: f64 = 0.02; +const A: f64 = PI / 4.0; +const D: f64 = PI / 2.0; +const T_END: f64 = 0.25; + +fn exact(x: f64, y: f64, z: f64, t: f64) -> (f64, f64, f64) { + let decay = (-D * D * NU * t).exp(); + let u = -A + * ((A * x).exp() * (A * y + D * z).sin() + (A * z).exp() * (A * x + D * y).cos()) + * decay; + let v = -A + * ((A * y).exp() * (A * z + D * x).sin() + (A * x).exp() * (A * y + D * z).cos()) + * decay; + let w = -A + * ((A * z).exp() * (A * x + D * y).sin() + (A * y).exp() * (A * z + D * x).cos()) + * decay; + (u, v, w) +} + +/// The boundary data as FACE AVERAGES (3 × 3 Gauss over the face): the +/// exact field has a non-zero normal velocity on the closed box, and its +/// face-centre samples leave an O(h²) net inflow a pure-Neumann projection +/// can only spread uniformly; the face-averaged fluxes of a divergence-free +/// field sum to zero to quadrature accuracy, so the box is compatible. +fn face_averaged(h: f64) -> impl Fn(f64, f64, f64, f64) -> (f64, f64, f64) { + const G: [f64; 3] = [-0.774_596_669_241_483_4, 0.0, 0.774_596_669_241_483_4]; + const W: [f64; 3] = [5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0]; + move |x: f64, y: f64, z: f64, t: f64| { + let on_x = x <= 0.0 || x >= 1.0; + let on_y = y <= 0.0 || y >= 1.0; + let on_z = z <= 0.0 || z >= 1.0; + if !(on_x || on_y || on_z) { + return exact(x, y, z, t); + } + let (mut u, mut v, mut w) = (0.0, 0.0, 0.0); + for (a, wa) in G.iter().zip(&W) { + for (b, wb) in G.iter().zip(&W) { + let (xx, yy, zz) = if on_x { + (x, y + 0.5 * h * a, z + 0.5 * h * b) + } else if on_y { + (x + 0.5 * h * a, y, z + 0.5 * h * b) + } else { + (x + 0.5 * h * a, y + 0.5 * h * b, z) + }; + let e = exact(xx, yy, zz, t); + u += 0.25 * wa * wb * e.0; + v += 0.25 * wa * wb * e.1; + w += 0.25 * wa * wb * e.2; + } + } + (u, v, w) + } +} + +struct Measurement { + l2: f64, + /// `max |div − mean(div)|`, with the mean reported (the residual + /// incompatibility of the face-averaged data: quadrature level). + max_div: f64, + mean_div: f64, + energy_ratio: f64, + steps: usize, +} + +fn measure(n: usize) -> Measurement { + let h = 1.0 / n as f64; + // dt ~ h²: the diffusion limit with a margin. + let dt = 0.25 * h * h / (4.0 * NU); + let steps = (T_END / dt).ceil() as usize; + let dt = T_END / steps as f64; + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: NU * RHO, + reference_velocity: 1.0, + reference_length: 1.0, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(face_averaged(h)); + let g = Grid { + nx: n, + ny: n, + nz: n, + dx: h, + dy: h, + dz: h, + }; + let mut f = Field::new(g); + for k in 0..n { + for j in 0..n { + for i in 0..=n { + f.u[g.uface(k, j, i)] = exact( + i as f64 * h, + (j as f64 + 0.5) * h, + (k as f64 + 0.5) * h, + 0.0, + ) + .0; + } + } + } + for k in 0..n { + for j in 0..=n { + for i in 0..n { + f.v[g.vface(k, j, i)] = exact( + (i as f64 + 0.5) * h, + j as f64 * h, + (k as f64 + 0.5) * h, + 0.0, + ) + .1; + } + } + } + for k in 0..=n { + for j in 0..n { + for i in 0..n { + f.w[g.wface(k, j, i)] = exact( + (i as f64 + 0.5) * h, + (j as f64 + 0.5) * h, + k as f64 * h, + 0.0, + ) + .2; + } + } + } + let e0 = f.kinetic_energy(RHO); + for _ in 0..steps { + solver.advance(&mut f, dt); + } + let e1 = f.kinetic_energy(RHO); + 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 { + let e = f.u[g.uface(k, j, i)] + - exact( + i as f64 * h, + (j as f64 + 0.5) * h, + (k as f64 + 0.5) * h, + T_END, + ) + .0; + sq += e * e * dv; + vol += dv; + } + } + } + for k in 0..n { + for j in 1..n { + for i in 0..n { + let e = f.v[g.vface(k, j, i)] + - exact( + (i as f64 + 0.5) * h, + j as f64 * h, + (k as f64 + 0.5) * h, + T_END, + ) + .1; + sq += e * e * dv; + vol += dv; + } + } + } + for k in 1..n { + for j in 0..n { + for i in 0..n { + let e = f.w[g.wface(k, j, i)] + - exact( + (i as f64 + 0.5) * h, + (j as f64 + 0.5) * h, + k as f64 * h, + T_END, + ) + .2; + sq += e * e * dv; + vol += dv; + } + } + } + let mut divs = Vec::with_capacity(n * n * n); + for k in 0..n { + for j in 0..n { + for i in 0..n { + divs.push( + (f.u[g.uface(k, j, i + 1)] - f.u[g.uface(k, j, i)]) / h + + (f.v[g.vface(k, j + 1, i)] - f.v[g.vface(k, j, i)]) / h + + (f.w[g.wface(k + 1, j, i)] - f.w[g.wface(k, j, i)]) / h, + ); + } + } + } + let mean_div = divs.iter().sum::() / divs.len() as f64; + let max_div = divs + .iter() + .fold(0.0_f64, |m, d| m.max((d - mean_div).abs())); + Measurement { + l2: (sq / vol).sqrt(), + max_div, + mean_div, + energy_ratio: e1 / e0, + steps, + } +} + +#[test] +fn beltrami_error_falls_under_space_time_refinement() { + let resolutions = [8usize, 16, 32]; + let ms: Vec = resolutions.iter().map(|&n| measure(n)).collect(); + let exact_ratio = (-2.0 * D * D * NU * T_END).exp(); + let errors: Vec = ms.iter().map(|m| m.l2).collect(); + for (i, &n) in resolutions.iter().enumerate() { + let rate = if i == 0 { + " -".to_string() + } else { + format!("{:5.2}", (errors[i - 1] / errors[i]).log2()) + }; + println!( + " n = {n:2} ({:5} steps) L2 = {:.6e} order {rate} E(T)/E(0) = {:.5} (exact {exact_ratio:.5}, error {:.2e}) max |div − mean| {:.2e} (mean {:.2e})", + ms[i].steps, + ms[i].l2, + ms[i].energy_ratio, + (ms[i].energy_ratio - exact_ratio).abs(), + ms[i].max_div, + ms[i].mean_div + ); + } + assert!( + errors.windows(2).all(|w| w[1] < w[0]), + "errors not monotone: {errors:?}" + ); + for w in errors.windows(2) { + let rate = (w[0] / w[1]).log2(); + assert!( + rate >= 0.75, + "observed order {rate:.3} below 0.75; errors {errors:?}" + ); + } + for m in &ms { + assert!(m.max_div < 1e-6, "max |div − mean| {:.3e}", m.max_div); + } + // The energy decay: the discretisation error at each rung bounds it. + let mut e_err: Vec = ms + .iter() + .map(|m| (m.energy_ratio - exact_ratio).abs()) + .collect(); + assert!( + e_err.windows(2).all(|w| w[1] < w[0]), + "energy error not falling: {e_err:?}" + ); + e_err.clear(); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_mms.rs b/crates/specialized/rtx-cfd/tests/embedded3_mms.rs new file mode 100644 index 0000000..52e922d --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_mms.rs @@ -0,0 +1,347 @@ +//! embedded3 gate 4: (a) z-invariance identity — the 2D manufactured +//! problem (`embedded_mms.rs`, no body) on the 2D embedded solver with the +//! multigrid Poisson vs the 3D solver at nz = 1 (dz = 1, z slip): the same +//! values over 200 steps; (b) a three-dimensional manufactured solution +//! marched to steady state on `n³` cubes: errors monotone under +//! refinement, orders in [0.75, 2.3], TVD's error below upwind's, and the +//! field divergence-free on every cell. + +use rtx_cfd::CfdConfig; +use rtx_cfd::solvers::incompressible::embedded3::{ + Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, +}; +use rtx_cfd::solvers::incompressible::{ + ConvectionScheme, EmbeddedParameters, EmbeddedPisoSolver, FlowField, PoissonSolverKind, +}; +use std::f64::consts::PI; + +const RHO: f64 = 1.0; +const MU: f64 = 0.05; + +// ---- the 2D manufactured problem (embedded_mms.rs) ---- +fn u2(x: f64, y: f64) -> f64 { + (PI * x).sin() * (PI * y).cos() +} +fn v2(x: f64, y: f64) -> f64 { + -(PI * x).cos() * (PI * y).sin() +} +fn source2(x: f64, y: f64) -> (f64, f64) { + let fx = RHO * 0.5 * PI * (2.0 * PI * x).sin() + + 2.0 * PI * PI * MU * u2(x, y) + + PI * (PI * x).cos() * (PI * y).sin(); + let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin() + + 2.0 * PI * PI * MU * v2(x, y) + + PI * (PI * x).sin() * (PI * y).cos(); + (fx, fy) +} +fn boundary2(x: f64, y: f64) -> (f64, f64) { + let u = if x <= 0.0 || x >= 1.0 { 0.0 } else { u2(x, y) }; + let v = if y <= 0.0 || y >= 1.0 { 0.0 } else { v2(x, y) }; + (u, v) +} + +fn fluid() -> Fluid { + Fluid { + density: RHO, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + } +} + +fn time_step(n: usize) -> f64 { + let h = 1.0 / n as f64; + 0.4 * (h * h / (4.0 * MU / RHO)).min(h) +} + +#[tokio::test] +async fn nz_one_reproduces_the_two_d_embedded_mms_march() { + let n = 16; + let h = 1.0 / n as f64; + let dt = time_step(n); + let config = CfdConfig::new() + .with_density(RHO) + .with_viscosity(MU) + .with_reference_velocity(1.0) + .with_reference_length(1.0); + let mut two = EmbeddedPisoSolver::new( + config, + EmbeddedParameters { + corrector_steps: 2, + tolerance: 1e-8, + poisson_solver: PoissonSolverKind::Multigrid, + ..EmbeddedParameters::default() + }, + ) + .expect("2D"); + two.set_momentum_source(|x, y, _| source2(x, y)); + two.set_boundary_velocity(|x, y, _| boundary2(x, y)); + let mut three = Solver::new( + fluid(), + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + boundaries: Boundaries { + z0: Side::SlipWall, + z1: Side::SlipWall, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + three.set_momentum_source(|x, y, _z, _t| { + let (fx, fy) = source2(x, y); + (fx, fy, 0.0) + }); + three.set_boundary_velocity(|x, y, _z, _t| { + let (u, v) = boundary2(x, y); + (u, v, 0.0) + }); + let mut a = FlowField::new(n, n, h, h).expect("field"); + let g = Grid { + nx: n, + ny: n, + nz: 1, + dx: h, + dy: h, + dz: 1.0, + }; + let mut b = Field::new(g); + for j in 0..n { + let y = (j as f64 + 0.5) * h; + a.u[(j, 0)] = boundary2(0.0, y).0; + a.u[(j, n)] = boundary2(1.0, y).0; + b.u[g.uface(0, j, 0)] = boundary2(0.0, y).0; + b.u[g.uface(0, j, n)] = boundary2(1.0, y).0; + } + for i in 0..n { + let x = (i as f64 + 0.5) * h; + a.v[(0, i)] = boundary2(x, 0.0).1; + a.v[(n, i)] = boundary2(x, 1.0).1; + b.v[g.vface(0, 0, i)] = boundary2(x, 0.0).1; + b.v[g.vface(0, n, i)] = boundary2(x, 1.0).1; + } + for step in 0..200 { + two.advance(&mut a, dt).await.expect("2D step"); + three.advance(&mut b, dt); + let mut worst = 0.0_f64; + for j in 0..n { + for i in 0..=n { + worst = worst.max((a.u[(j, i)] - b.u[g.uface(0, j, i)]).abs()); + } + } + for j in 0..=n { + for i in 0..n { + worst = worst.max((a.v[(j, i)] - b.v[g.vface(0, j, i)]).abs()); + } + } + for j in 0..n { + for i in 0..n { + worst = worst.max((a.p[(j, i)] - b.p[g.cell(0, j, i)]).abs()); + } + } + assert!( + worst == 0.0, + "step {step}: 3D departs from the 2D embedded MMS march by {worst:.3e}" + ); + } + println!(" 200 steps of the manufactured problem value-identical to the 2D embedded solver"); +} + +// ---- the 3D manufactured solution ---- +// u = sin πx cos πy cos πz, v = cos πx sin πy cos πz, w = −2 cos πx cos πy sin πz +// (divergence-free), p = sin πx sin πy sin πz; source = ρ(u·∇)u + ∇p − μ∇²u. +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() +} +fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let (sx, cx) = (PI * x).sin_cos(); + let (sy, cy) = (PI * y).sin_cos(); + let (sz, cz) = (PI * z).sin_cos(); + let (u, v, w) = (u3(x, y, z), v3(x, y, z), w3(x, y, z)); + // Gradients. + let (ux, uy, uz) = (PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz); + let (vx, vy, vz) = (-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz); + let (wx, wy, wz) = ( + 2.0 * PI * sx * cy * sz, + 2.0 * PI * cx * sy * sz, + -2.0 * PI * cx * cy * cz, + ); + let (px, py, pz) = (PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz); + // ∇²(product of three π-trig functions) = −3π² (itself). + let lap = -3.0 * PI * PI; + let fx = RHO * (u * ux + v * uy + w * uz) + px - MU * lap * u; + let fy = RHO * (u * vx + v * vy + w * vz) + py - MU * lap * v; + let fz = RHO * (u * wx + v * wy + w * wz) + pz - MU * lap * w; + (fx, fy, fz) +} +/// The exact field on the cube's boundary with the normal components +/// snapped to their analytic zero. +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) +} + +struct Measurement { + l2_velocity: f64, + max_div: f64, + steps: usize, +} + +fn measure3(n: usize, scheme: ConvectionScheme) -> Measurement { + let h = 1.0 / n as f64; + let dt = time_step(n); + let mut solver = Solver::new( + fluid(), + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + convection_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)); + let g = Grid { + nx: n, + ny: n, + nz: n, + dx: h, + dy: h, + dz: h, + }; + let mut f = Field::new(g); + solver.initialize(&mut f); + let mut steps = 0; + for step in 0..200_000 { + let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone()); + solver.advance(&mut f, dt); + steps = step + 1; + 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 (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 { + 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 k in 0..n { + for j in 1..n { + for i in 0..n { + 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 { + 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; + } + } + } + Measurement { + l2_velocity: (sq / vol).sqrt(), + max_div: f.max_divergence(), + steps, + } +} + +fn ladder(resolutions: &[usize], scheme: ConvectionScheme) -> Vec { + let ms: Vec = resolutions.iter().map(|&n| measure3(n, scheme)).collect(); + let errors: Vec = ms.iter().map(|m| m.l2_velocity).collect(); + for (i, &n) in resolutions.iter().enumerate() { + let rate = if i == 0 { + " -".to_string() + } else { + format!("{:5.2}", (errors[i - 1] / errors[i]).log2()) + }; + println!( + " {scheme:?} n = {n:3} ({:5} steps) L2 velocity {:.6e} order {rate} max |div| {:.3e}", + ms[i].steps, errors[i], ms[i].max_div + ); + } + assert!( + errors.windows(2).all(|w| w[1] < w[0]), + "{scheme:?}: errors not monotone {errors:?}" + ); + for w in errors.windows(2) { + let rate = (w[0] / w[1]).log2(); + assert!( + rate > 0.75 && rate < 2.3, + "{scheme:?}: observed order {rate:.3} outside [0.75, 2.3]; errors {errors:?}" + ); + } + for m in &ms { + assert!(m.max_div < 1e-5, "{scheme:?}: max |div| {:.3e}", m.max_div); + } + ms +} + +#[test] +fn embedded3_mms_orders() { + let resolutions = [12usize, 24]; + let up = ladder(&resolutions, ConvectionScheme::Upwind); + let tvd = ladder(&resolutions, ConvectionScheme::TvdVanAlbada); + let tvd_rate = (tvd[0].l2_velocity / tvd[1].l2_velocity).log2(); + assert!(tvd_rate > 1.1, "TVD order {tvd_rate:.3} not above 1.1"); + for (a, b) in up.iter().zip(&tvd) { + assert!( + b.l2_velocity < a.l2_velocity, + "TVD error not below upwind's" + ); + } +} + +#[test] +#[ignore = "the three-rung ladder to n = 48 (minutes on the host)"] +fn embedded3_mms_orders_three_rungs() { + let resolutions = [12usize, 24, 48]; + ladder(&resolutions, ConvectionScheme::Upwind); + ladder(&resolutions, ConvectionScheme::TvdVanAlbada); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_poiseuille.rs b/crates/specialized/rtx-cfd/tests/embedded3_poiseuille.rs new file mode 100644 index 0000000..ae47023 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_poiseuille.rs @@ -0,0 +1,271 @@ +//! embedded3 gate 6: plane Poiseuille flow on the 3D solver. (a) At +//! `nz = 1` (`dz = 1`, z sides slip) the 3D host step reproduces the 2D +//! embedded solver (no body, multigrid Poisson) to the value over 200 +//! steps; (b) the 2D `poiseuille.rs` gates on the 3D field at nz = 1 and on +//! a periodic-z extrusion: `|u − û| < 1e-7`, `max |v|, |w| < 1e-7`, +//! `p spread < 1e-6`, with û the discrete channel profile. + +use rtx_cfd::CfdConfig; +use rtx_cfd::solvers::incompressible::embedded3::{ + Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, +}; +use rtx_cfd::solvers::incompressible::{ + EmbeddedParameters, EmbeddedPisoSolver, FlowField, PoissonSolverKind, +}; + +const MU: f64 = 0.1; +const G: f64 = 0.8; + +fn discrete_profile(n: usize) -> Vec { + let h = 1.0 / n as f64; + let rhs_value = -G * h * h / MU; + let mut diag = vec![-2.0; n]; + diag[0] = -3.0; + diag[n - 1] = -3.0; + let mut rhs = vec![rhs_value; n]; + let upper = vec![1.0; n]; + for j in 1..n { + let factor = 1.0 / diag[j - 1]; + diag[j] -= factor * upper[j - 1]; + rhs[j] -= factor * rhs[j - 1]; + } + let mut u = vec![0.0; n]; + u[n - 1] = rhs[n - 1] / diag[n - 1]; + for j in (0..n - 1).rev() { + u[j] = (rhs[j] - upper[j] * u[j + 1]) / diag[j]; + } + u +} + +fn fluid() -> Fluid { + Fluid { + density: 1.0, + viscosity: MU, + reference_velocity: 1.0, + reference_length: 1.0, + } +} + +/// The inlet/outlet carry the discrete profile (the embedded solvers +/// re-stamp every Velocity side from the boundary function each step); +/// the walls are no-slip. +fn profile_boundary(n: usize) -> impl Fn(f64, f64) -> f64 + Clone { + let u_hat = discrete_profile(n); + let h = 1.0 / n as f64; + move |x: f64, y: f64| { + if x <= 0.0 || x >= 1.0 { + let j = ((y / h - 0.5).round().max(0.0) as usize).min(n - 1); + u_hat[j] + } else { + 0.0 + } + } +} + +fn solver3(n: usize, nz_periodic: bool) -> Solver { + let z = if nz_periodic { + Side::Periodic + } else { + Side::SlipWall + }; + let mut s = Solver::new( + fluid(), + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + boundaries: Boundaries { + z0: z, + z1: z, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + s.set_momentum_source(|_x, _y, _z, _t| (G, 0.0, 0.0)); + let ub = profile_boundary(n); + s.set_boundary_velocity(move |x, y, _z, _t| (ub(x, y), 0.0, 0.0)); + s +} + +fn field3(n: usize, nz: usize, dz: f64) -> Field { + let h = 1.0 / n as f64; + let g = Grid { + nx: n, + ny: n, + nz, + dx: h, + dy: h, + dz, + }; + let mut f = Field::new(g); + let u_hat = discrete_profile(n); + for k in 0..nz { + for (j, &uj) in u_hat.iter().enumerate() { + f.u[g.uface(k, j, 0)] = uj; + f.u[g.uface(k, j, n)] = uj; + } + } + f +} + +fn field2(n: usize) -> FlowField { + let h = 1.0 / n as f64; + let mut f = FlowField::new(n, n, h, h).expect("field"); + let u_hat = discrete_profile(n); + for (j, &uj) in u_hat.iter().enumerate() { + f.u[(j, 0)] = uj; + f.u[(j, n)] = uj; + } + f +} + +fn dt_for(n: usize) -> f64 { + let h = 1.0 / n as f64; + 0.4 * (h * h / (4.0 * MU)).min(h) +} + +/// (a) the value identity at nz = 1. +#[tokio::test] +async fn nz_one_is_the_two_d_embedded_solver() { + let n = 16; + let dt = dt_for(n); + let config = CfdConfig::new() + .with_density(1.0) + .with_viscosity(MU) + .with_reference_velocity(1.0) + .with_reference_length(1.0); + let mut two = EmbeddedPisoSolver::new( + config, + EmbeddedParameters { + corrector_steps: 2, + tolerance: 1e-8, + poisson_solver: PoissonSolverKind::Multigrid, + ..EmbeddedParameters::default() + }, + ) + .expect("2D solver"); + two.set_momentum_source(|_x, _y, _t| (G, 0.0)); + let ub = profile_boundary(n); + two.set_boundary_velocity(move |x, y, _t| (ub(x, y), 0.0)); + let mut three = solver3(n, false); + let mut a = field2(n); + let mut b = field3(n, 1, 1.0); + let g = b.grid; + let mut signed_zero = 0usize; + for step in 0..200 { + two.advance(&mut a, dt).await.expect("2D step"); + three.advance(&mut b, dt); + let mut worst = 0.0_f64; + for j in 0..n { + for i in 0..=n { + let (x, y) = (a.u[(j, i)], b.u[g.uface(0, j, i)]); + if x != y { + worst = worst.max((x - y).abs()); + } else if x.to_bits() != y.to_bits() { + signed_zero += 1; + } + } + } + for j in 0..=n { + for i in 0..n { + let (x, y) = (a.v[(j, i)], b.v[g.vface(0, j, i)]); + if x != y { + worst = worst.max((x - y).abs()); + } else if x.to_bits() != y.to_bits() { + signed_zero += 1; + } + } + } + for j in 0..n { + for i in 0..n { + let (x, y) = (a.p[(j, i)], b.p[g.cell(0, j, i)]); + if x != y { + worst = worst.max((x - y).abs()); + } + } + } + assert!( + worst == 0.0, + "step {step}: the 3D field departs from the 2D embedded solver by {worst:.3e}" + ); + } + println!( + " 200 steps value-identical to the 2D embedded solver ({signed_zero} ±0 sign differences)" + ); +} + +struct Measurement { + max_u_vs_discrete: f64, + max_v: f64, + max_w: f64, + p_spread: f64, + steps: usize, +} + +fn measure(n: usize, nz: usize, dz: f64, periodic: bool) -> Measurement { + let dt = dt_for(n); + let mut solver = solver3(n, periodic); + let mut f = field3(n, nz, dz); + let g = f.grid; + let u_hat = discrete_profile(n); + let mut steps = 0; + for step in 0..200_000 { + let before = f.u.clone(); + solver.advance(&mut f, dt); + steps = step + 1; + let change = + f.u.iter() + .zip(&before) + .fold(0.0_f64, |m, (a, b)| m.max((a - b).abs())) + / dt; + if change < 1e-8 { + break; + } + } + let mut max_u_vs_discrete = 0.0_f64; + for k in 0..nz { + for (j, &uj) in u_hat.iter().enumerate() { + for i in 1..n { + max_u_vs_discrete = max_u_vs_discrete.max((f.u[g.uface(k, j, i)] - uj).abs()); + } + } + } + let max_v = f.v.iter().fold(0.0_f64, |m, v| m.max(v.abs())); + let max_w = f.w.iter().fold(0.0_f64, |m, v| m.max(v.abs())); + let (mut p_min, mut p_max) = (f64::INFINITY, f64::NEG_INFINITY); + for &p in &f.p { + p_min = p_min.min(p); + p_max = p_max.max(p); + } + Measurement { + max_u_vs_discrete, + max_v, + max_w, + p_spread: p_max - p_min, + steps, + } +} + +/// (b) the 2D gates on the 3D field. +#[test] +fn poiseuille_is_the_discrete_profile_in_three_d() { + for (n, nz, dz, periodic) in [ + (16usize, 1usize, 1.0, false), + (16, 4, 1.0 / 16.0, true), + (32, 1, 1.0, false), + ] { + let m = measure(n, nz, dz, periodic); + println!( + " n {n} nz {nz} periodic {periodic}: {} steps; |u − û| {:.3e}, max |v| {:.3e}, max |w| {:.3e}, p spread {:.3e}", + m.steps, m.max_u_vs_discrete, m.max_v, m.max_w, m.p_spread + ); + assert!( + m.max_u_vs_discrete < 1e-7, + "u departs from the discrete profile by {:.3e}", + m.max_u_vs_discrete + ); + assert!(m.max_v < 1e-7, "spurious transverse flow {:.3e}", m.max_v); + assert!(m.max_w < 1e-7, "spurious spanwise flow {:.3e}", m.max_w); + assert!(m.p_spread < 1e-6, "spurious pressure {:.3e}", m.p_spread); + } +}