From f6add276c080fe9b813f5345be9902c13342794c Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 18 Sep 2026 03:59:57 -0500 Subject: [PATCH] =?UTF-8?q?embedded3:=20cut=20predictor=20convection=20car?= =?UTF-8?q?ries=20=CF=81=20(host=20+=20e3=5Fcut.cu;=20density-scaling=20pi?= =?UTF-8?q?n);=20operator=20load=20route=20includes=20the=20wall=20exchang?= =?UTF-8?q?e=20(exchange.rs);=20reconstructed=5Fparts,=20probe=20aperture?= =?UTF-8?q?=20floor=20knob;=20dfg=5Fsplit=20diagnostic=20test?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/src/kernels/cuda/e3_cut.cu | 3 +- .../incompressible/embedded3/cutwall.rs | 30 +++- .../incompressible/embedded3/exchange.rs | 143 ++++++++++++++++++ .../solvers/incompressible/embedded3/loads.rs | 22 ++- .../solvers/incompressible/embedded3/mod.rs | 1 + .../embedded3/step/cut_predictor.rs | 4 +- .../incompressible/embedded3/step/mod.rs | 5 + .../solvers/incompressible/embedded3/wall.rs | 7 + .../tests/embedded3_density_scaling.rs | 101 +++++++++++++ .../rtx-cfd/tests/embedded3_dfg_2z.rs | 17 +++ .../rtx-cfd/tests/embedded3_dfg_split.rs | 130 ++++++++++++++++ 11 files changed, 454 insertions(+), 9 deletions(-) create mode 100644 crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_density_scaling.rs create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_dfg_split.rs diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu index c1698e8..be389ba 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu @@ -176,7 +176,8 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu if (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d]; else if (sides[d][0] == SIDE_VELOCITY) diff -= mu * g_minus * a_d * (u0 - beyond_m) / (0.5 * h[d]); } - conv -= mass_out * u0; + /* The mass fluxes above are volume fluxes: the momentum flux carries rho. */ + conv = rho * (conv - mass_out * u0); int cpi = cut_cell(g, cp[0], cp[1], cp[2]); int cmi = cut_cell(g, cm[0], cm[1], cm[2]); double p_plus = cpi >= 0 ? f.p[cpi] : 0.0; 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 e37a1cb..12fcf8b 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -207,6 +207,8 @@ impl Mask { step_apertures: None, step_open: None, merge_master: Vec::new(), + scheme: crate::solvers::incompressible::ConvectionScheme::Upwind, + density: 1.0, }; mask.compute_merging(None); Ok(mask) @@ -564,7 +566,8 @@ impl Mask { /// 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]]) + let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, None)?; + Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]) } /// The reconstructed wall route (S2-1 remedy): on every wall polygon @@ -585,12 +588,26 @@ impl Mask { t: f64, planes: Option<(usize, usize)>, ) -> Option<[f64; 3]> { + let (p, s) = self.cut_wall_force_reconstructed_parts(body, f, mu, t, planes)?; + Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]]) + } + + /// The reconstructed route split into its pressure and shear parts. + pub fn cut_wall_force_reconstructed_parts( + &self, + body: &Body, + f: &Field, + mu: f64, + t: f64, + planes: Option<(usize, usize)>, + ) -> Option<([f64; 3], [f64; 3])> { let cut = self.cut.as_ref()?; let g = self.grid; let (k0, k1) = planes.unwrap_or((0, g.nz)); let h = g.dx.min(g.dy).min(g.dz); let (d1, d2) = (h, 2.0 * h); let mut force = [0.0; 3]; + let mut shear = [0.0; 3]; for (idx, w) in cut.wall.iter().enumerate() { let area = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt(); if area == 0.0 || !self.cell_fluid[idx] { @@ -645,10 +662,10 @@ impl Mask { let dn = wall_gradient(t1[c] - ts[c], t2[c] - ts[c]); // Traction on the body = −(fluid stress on the fluid side): // the shear the fluid exerts on the wall along +t. - force[c] += mu * dn * area; + shear[c] += mu * dn * area; } } - Some(force) + Some((force, shear)) } /// The cut-cell load route restricted to the cells (and faces) of the @@ -663,8 +680,13 @@ impl Mask { (k0, k1): (usize, usize), ) -> Option<[f64; 3]> { let (p, s) = self.cut_wall_force_parts_in(body, f, mu, t, Some((k0, k1)))?; + let x = self.cut_wall_exchange_force(body, f, mu, self.density, t, Some((k0, k1)))?; let lz = (k1 - k0) as f64 * self.grid.dz; - Some([(p[0] + s[0]) / lz, (p[1] + s[1]) / lz, (p[2] + s[2]) / lz]) + Some([ + (p[0] + s[0] + x[0]) / lz, + (p[1] + s[1] + x[1]) / lz, + (p[2] + s[2] + x[2]) / lz, + ]) } /// The cut-cell load route split into its pressure and shear parts. diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs new file mode 100644 index 0000000..05599ed --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs @@ -0,0 +1,143 @@ +//! The wall-exchange part of the operator load route (S2-1 remedy, second +//! form). A fluid face control volume next to a prescribed (ghost or +//! solid) face still exchanges momentum with it: the diffusive flux +//! `μ g A (U_b − u_f)/h` through the half-aperture face at the cell +//! centre and the convective flux `ρ m (u_face − u_f)` on the same face. +//! Both are forces the wall exerts on the fluid that the closure polygon's +//! shear `μ A_w (u_f − U_b)/d_f` does not carry, so the operator route +//! read short of the box route by exactly this exchange (a conservation +//! gap that did not shrink with h). Summed here with the predictor's own +//! flux formulas, the operator route closes the discrete momentum balance. +use super::body::Body; +use super::field::Field; +use super::wall::{FaceKind, Mask}; +use crate::solvers::incompressible::ConvectionScheme; + +impl Mask { + /// The momentum the fluid's face control volumes exchange with the + /// prescribed faces beside them, as a force on the body (the negative + /// of the force on the fluid), over the z planes `planes` (all when + /// `None`). `None` without a cut geometry. + pub fn cut_wall_exchange_force( + &self, + body: &Body, + f: &Field, + mu: f64, + rho: f64, + t: f64, + planes: Option<(usize, usize)>, + ) -> Option<[f64; 3]> { + let _ = body; + let _ = t; + self.cut.as_ref()?; + let g = self.grid; + let lat = self.lattice(); + let h = [g.dx, g.dy, g.dz]; + let area = [g.dy * g.dz, g.dx * g.dz, g.dx * g.dy]; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let (k0, k1) = planes.unwrap_or((0, nz)); + let w_range = if self.periodic_z { 0..nz } else { 1..nz }; + let vals: [&[f64]; 3] = [&f.u, &f.v, &f.w]; + let scheme = self.scheme; + let kind = |cc: usize, idx: usize| match cc { + 0 => self.u_kind[idx], + 1 => self.v_kind[idx], + _ => self.w_kind[idx], + }; + let val = |cc: usize, q: [i64; 3]| lat.face(cc, q).map(|i| vals[cc][i]); + let ap = |cc: usize, q: [i64; 3]| self.aperture(cc, q); + let e = |d: usize| { + let mut v = [0i64; 3]; + v[d] = 1; + v + }; + let add = + |a: [i64; 3], b: [i64; 3], s: i64| [a[0] + s * b[0], a[1] + s * b[1], a[2] + s * b[2]]; + let upwind = |m: f64, up: f64, dn: f64| if m >= 0.0 { up } else { dn }; + let mut force = [0.0; 3]; + for c in 0..3 { + let (ir, jr, kr) = match c { + 0 => (1..nx, 0..ny, k0..k1), + 1 => (0..nx, 1..ny, k0..k1), + _ => (0..nx, 0..ny, w_range.start.max(k0)..w_range.end.min(k1)), + }; + let ec = e(c); + for k in kr { + for j in jr.clone() { + for i in ir.clone() { + let p = [i as i64, j as i64, k as i64]; + let idx = lat.face(c, p).expect("face"); + if kind(c, idx) != FaceKind::Fluid { + continue; + } + let cv = self.cv_geometry(c, p); + let u0 = vals[c][idx]; + let cell_minus = add(p, ec, -1); + let cell_plus = p; + for d in 0..3 { + let ed = e(d); + let a_d = area[d]; + let up1 = val(c, add(p, ed, 1)); + let up2 = val(c, add(p, ed, 2)); + let dn1 = val(c, add(p, ed, -1)); + let dn2 = val(c, add(p, ed, -2)); + let (m_plus, m_minus) = if d == c { + let f_up = + ap(c, add(p, ec, 1)).unwrap_or(cv.alpha) * up1.unwrap_or(u0); + let f_dn = + ap(c, add(p, ec, -1)).unwrap_or(cv.alpha) * dn1.unwrap_or(u0); + let f0 = cv.alpha * u0; + (0.5 * (f0 + f_up) * a_d, 0.5 * (f_dn + f0) * a_d) + } else { + let flux = |q: [i64; 3]| { + ap(d, q).unwrap_or(1.0) * val(d, q).unwrap_or(0.0) + }; + ( + 0.5 * (flux(add(cell_minus, ed, 1)) + + flux(add(cell_plus, ed, 1))) + * a_d, + 0.5 * (flux(cell_minus) + flux(cell_plus)) * a_d, + ) + }; + // Plus side: a prescribed neighbour face. + if let Some(fp) = lat.face(c, add(p, ed, 1)) { + if kind(c, fp) != FaceKind::Fluid { + let un = vals[c][fp]; + let delta = if scheme == ConvectionScheme::Upwind { + 0.0 + } else if m_plus >= 0.0 { + scheme.face_correction(dn1, u0, un) + } else { + scheme.face_correction(up2, un, u0) + }; + let u_face = upwind(m_plus, u0, un) + delta; + let on_fluid = -rho * m_plus * (u_face - u0) + + mu * cv.ap[d][1] * a_d * (un - u0) / h[d]; + force[c] -= on_fluid; + } + } + // Minus side. + if let Some(fm) = lat.face(c, add(p, ed, -1)) { + if kind(c, fm) != FaceKind::Fluid { + let ud = vals[c][fm]; + let delta = if scheme == ConvectionScheme::Upwind { + 0.0 + } else if m_minus >= 0.0 { + scheme.face_correction(dn2, ud, u0) + } else { + scheme.face_correction(up1, u0, ud) + }; + let u_face = upwind(m_minus, ud, u0) + delta; + let on_fluid = rho * m_minus * (u_face - u0) + + mu * cv.ap[d][0] * a_d * (ud - u0) / h[d]; + force[c] -= on_fluid; + } + } + } + } + } + } + } + Some(force) + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs index b75a3b2..72b186f 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs @@ -13,6 +13,16 @@ use super::body::Body; use super::field::Field; use super::wall::{FaceKind, Mask, linear_fit, stencil_nodes, z_planes}; +/// Minimum face aperture for a velocity node to enter a probe's +/// interpolation (`RTX_E3_PROBE_MIN_APERTURE`, default 0.5; 0 keeps every +/// fluid face, the reading before S2-1's fix). +fn probe_min_aperture() -> f64 { + std::env::var("RTX_E3_PROBE_MIN_APERTURE") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(0.5) +} + #[derive(Debug, Clone, Copy, PartialEq)] pub struct SurfaceForce { pub f: [f64; 3], @@ -154,10 +164,16 @@ impl Mask { let foot_c = (foot.0, foot.1, if zq == z { foot.2 } else { zq }); let nodes = stencil_nodes((x, y, zq), c, g, self.periodic_z(), |_| None); let values: &[f64] = [&f.u, &f.v, &f.w][c]; + // A cut face's velocity lives on the open part of the face, not + // at its centre: a node whose centre lies in the body would put + // the wall profile's value h/2 too deep and bias the wall + // gradient at O(1). Nodes below the minimum aperture are dropped + // and the point's own wall intercept takes their place. + let min_aperture = probe_min_aperture(); let fluid = |idx: usize| match c { - 0 => self.u_kind(idx) == FaceKind::Fluid, - 1 => self.v_kind(idx) == FaceKind::Fluid, - _ => self.w_kind(idx) == FaceKind::Fluid, + 0 => self.u_kind(idx) == FaceKind::Fluid && self.a_u(idx) >= min_aperture, + 1 => self.v_kind(idx) == FaceKind::Fluid && self.a_v(idx) >= min_aperture, + _ => self.w_kind(idx) == FaceKind::Fluid && self.a_w(idx) >= min_aperture, }; if nodes.iter().all(|n| fluid(n.idx)) { out[c] = nodes.iter().map(|n| n.weight * values[n.idx]).sum(); 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 c46facf..d060461 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -9,6 +9,7 @@ pub mod body; pub mod cut; pub mod cutwall; +pub mod exchange; pub mod export_vtk; pub mod field; pub mod grid; 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 864c4c7..7798d10 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 @@ -201,7 +201,9 @@ impl Solver { // 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; + // The mass fluxes above are volume fluxes: the momentum flux carries ρ + // (inertia, diffusion and the pressure are dynamic). + let conv = rho * (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/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index ffcb5c5..3aea502 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 @@ -213,6 +213,11 @@ impl Solver { Mask::build_cut_from(body, g, t, self.params.boundaries, prev) } } + .map(|mut m| { + m.scheme = self.params.convection_scheme; + m.density = self.fluid.density; + m + }) .expect("embedded mask") } 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 a3adfe5..2cd3881 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -92,6 +92,11 @@ pub struct Mask { /// (`usize::MAX` = its own row) — a small cell shares its pressure /// unknown with its largest active face neighbour in the projection. pub(super) merge_master: Vec, + /// The predictor's convection scheme (the exchange route replicates + /// its limited fluxes on the faces next to prescribed ones). + pub(super) scheme: crate::solvers::incompressible::ConvectionScheme, + /// The fluid's density (the exchange route's convective flux). + pub(super) density: f64, } /// The z lattice position of a query: the lower plane index, the upper @@ -502,6 +507,8 @@ impl Mask { step_apertures: None, step_open: None, merge_master: Vec::new(), + scheme: crate::solvers::incompressible::ConvectionScheme::Upwind, + density: 1.0, }) } diff --git a/crates/specialized/rtx-cfd/tests/embedded3_density_scaling.rs b/crates/specialized/rtx-cfd/tests/embedded3_density_scaling.rs new file mode 100644 index 0000000..84c160b --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_density_scaling.rs @@ -0,0 +1,101 @@ +//! Density-scaling pin for the cut-cell predictor (host): the same flow at +//! `ρ` and `1000 ρ` with `μ` scaled alike is the same velocity field and a +//! pressure scaled by 1000 — every term of the momentum equation carries +//! `ρ` (the convection term used to be a bare volume flux times velocity, +//! which starved every ρ = 1000 cut-cell run of convection). +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; + +fn run(rho: f64, moving: bool) -> Field { + let n = 16; + let h = 1.0 / n as f64; + let g = Grid::cubic(2 * n, n, n, h); + let nu = 1e-2; + let mut solver = Solver::new( + Fluid { + density: rho, + viscosity: rho * nu, + reference_velocity: 1.0, + reference_length: 0.3, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-11, + convection_scheme: ConvectionScheme::TvdVanAlbada, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + }, + max_surface_speed: if moving { Some(0.5) } else { None }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|x, _, _, _| { + if x <= 0.0 { + (1.0, 0.0, 0.0) + } else { + (0.0, 0.0, 0.0) + } + }); + let xc = move |t: f64| 0.7 + if moving { 0.1 * (3.0 * t).sin() } else { 0.0 }; + let body = Body::from_sdf(move |x, y, z, t| { + ((x - xc(t)).powi(2) + (y - 0.5_f64).powi(2) + (z - 0.5_f64).powi(2)).sqrt() - 0.15 + }) + .with_surface_velocity(move |_, _, _, t| { + (if moving { 0.3 * (3.0 * t).cos() } else { 0.0 }, 0.0, 0.0) + }); + if moving { + solver.set_moving_body(body); + } else { + solver.set_body(body); + } + let mut field = Field::new(g); + for k in 0..n { + for j in 0..n { + for i in 0..=2 * n { + field.u[g.uface(k, j, i)] = 1.0; + } + } + } + solver.initialize(&mut field); + let dt = 0.2 * h; + for _ in 0..40 { + solver.advance(&mut field, dt); + } + field +} + +fn compare(moving: bool) { + let a = run(1.0, moving); + let b = run(1000.0, moving); + let max = |x: &[f64], y: &[f64], s: f64| { + x.iter() + .zip(y) + .map(|(p, q)| (p - q / s).abs()) + .fold(0.0, f64::max) + }; + let du = max(&a.u, &b.u, 1.0) + .max(max(&a.v, &b.v, 1.0)) + .max(max(&a.w, &b.w, 1.0)); + let dp = max(&a.p, &b.p, 1000.0); + let pscale = a.p.iter().fold(0.0f64, |m, p| m.max(p.abs())); + println!(" moving {moving}: max |Δu| {du:.3e}, max |Δp/1000| {dp:.3e} (p scale {pscale:.3e})"); + assert!(du < 1e-9, "velocity is not density-invariant: {du:.3e}"); + assert!( + dp < 1e-9 * pscale.max(1.0), + "pressure does not scale with density: {dp:.3e}" + ); +} + +#[test] +fn cut_cell_flow_is_density_invariant_at_rest() { + compare(false); +} + +#[test] +fn cut_cell_flow_is_density_invariant_moving() { + compare(true); +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs index 9480924..a8ddc29 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_2z.rs @@ -201,6 +201,23 @@ fn dfg_3d_2z_on_the_device() { println!(" instant written to {}", path.display()); } let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples"); + { + // S2-1 diagnosis: each wall route split into its pressure and shear parts. + let body = solver.body().expect("body"); + let (po, so) = mask + .cut_wall_force_parts(body, &field, RHO * NU, solver.time()) + .expect("parts"); + let (pr, sr) = mask + .cut_wall_force_reconstructed_parts(body, &field, RHO * NU, solver.time(), None) + .expect("parts"); + println!( + " SPLIT ny {ny}: operator c_D pressure {:.4} + shear {:.4}; reconstructed pressure {:.4} + shear {:.4}", + coef * po[0], + coef * so[0], + coef * pr[0], + coef * sr[0] + ); + } println!( " FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; reconstructed {:.4}, {:.2e} from CV) c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {:.5}) Δp {dp:.4} — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s", ((cd - cd_cv) / cd).abs(), diff --git a/crates/specialized/rtx-cfd/tests/embedded3_dfg_split.rs b/crates/specialized/rtx-cfd/tests/embedded3_dfg_split.rs new file mode 100644 index 0000000..309a6ab --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_dfg_split.rs @@ -0,0 +1,130 @@ +//! S2-1 diagnosis (host): DFG 3D-2Z at a coarse rung with every load +//! route split into pressure and shear parts — which part of the wall +//! routes departs from the box route. `RTX_E3_DFG_NY` (default 31). +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; + +const H: f64 = 0.41; +const L: f64 = 2.5; +const D: f64 = 0.1; +const CX: f64 = 0.5; +const CY: f64 = 0.2; +const U_M: f64 = 0.45; +const U_BAR: f64 = 4.0 / 9.0 * U_M; +const RHO: f64 = 1.0; +const NU: f64 = 1e-3; + +fn inflow(y: f64, z: f64) -> f64 { + 16.0 * U_M * y * z * (H - y) * (H - z) / (H * H * H * H) +} + +#[test] +#[ignore = "host DFG at ny 31 with the routes split (about half an hour)"] +fn dfg_routes_split_on_the_host() { + let ny: usize = std::env::var("RTX_E3_DFG_NY") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(31); + let h = H / ny as f64; + let nx = (L / h).round() as usize; + let nz = ny; + let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU)); + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: RHO * NU, + reference_velocity: U_BAR, + reference_length: D, + }, + Parameters { + corrector_steps: 2, + tolerance: 1e-8, + convection_scheme: ConvectionScheme::TvdVanAlbada, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + x1: Side::PressureOutlet, + ..Boundaries::default() + }, + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|x, y, z, _t| { + if x <= 0.0 { + (inflow(y, z), 0.0, 0.0) + } else { + (0.0, 0.0, 0.0) + } + }); + solver.set_body(Body::extruded( + rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D), + H, + )); + let g = Grid::cubic(nx, ny, nz, h); + let mut field = Field::new(g); + for k in 0..nz { + for j in 0..ny { + let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h); + for i in 0..=nx { + field.u[g.uface(k, j, i)] = u0; + } + } + } + solver.initialize(&mut field); + let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H); + let steps = (8.0 / dt).ceil() as usize; + let start = std::time::Instant::now(); + let mut last = (0.0, 0.0); + for step in 0..steps { + let r = solver.advance(&mut field, dt); + if (step + 1) % (steps / 20).max(1) == 0 || step + 1 == steps { + let t = solver.time(); + let mask = solver.mask().unwrap(); + let body = solver.body().unwrap(); + let mu = RHO * NU; + let (po, so) = mask.cut_wall_force_parts(body, &field, mu, t).unwrap(); + let ex = mask + .cut_wall_exchange_force(body, &field, mu, RHO, t, None) + .unwrap(); + let (pr, sr) = mask + .cut_wall_force_reconstructed_parts(body, &field, mu, t, None) + .unwrap(); + let margin = 3.0 * D; + let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2); + let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2); + let bx = ( + ci(CX - margin), + ci(CX + margin), + cj(CY - 0.15), + cj(CY + 0.15), + 0, + nz, + ); + let fcv = mask.control_volume_force_with_walls(&field, dt, RHO, mu, None, bx, true); + let cd = |f: [f64; 3]| coef * f[0]; + println!( + " t {t:7.3}: c_D operator {:.4} (p {:.4} + s {:.4} + exchange {:.4}) | reconstructed {:.4} (p {:.4} + s {:.4}) | box {:.4}; residual {:.1e} [{:.0} s]", + cd(po) + cd(so) + cd(ex), + cd(po), + cd(so), + cd(ex), + cd(pr) + cd(sr), + cd(pr), + cd(sr), + cd(fcv), + r.final_residual, + start.elapsed().as_secs_f64() + ); + let now = (cd(po) + cd(so) + cd(ex), cd(fcv)); + if (now.0 - last.0).abs() < 1e-4 * now.0.abs() + && (now.1 - last.1).abs() < 1e-4 * now.1.abs() + && t > 2.0 + { + println!(" settled"); + break; + } + last = now; + } + } +}