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 a2d96a1..caf696a 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -228,8 +228,13 @@ impl Mask { old.and_then(|o| o.cut.as_ref()) .map_or(v, |oc| v.max(oc.vol[idx])) }; + // `RTX_E3_MERGE_FRACTION` overrides the threshold (a study knob). + let threshold = std::env::var("RTX_E3_MERGE_FRACTION") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(MERGE_FRACTION); let small: Vec = (0..n) - .map(|idx| self.cell_active(idx) && frac(idx) < MERGE_FRACTION) + .map(|idx| self.cell_active(idx) && frac(idx) < threshold) .collect(); let mut master = vec![usize::MAX; n]; for idx in (0..n).filter(|&i| small[i]) { @@ -256,17 +261,37 @@ impl Mask { } /// Set the step-averaged apertures and the space-time classification - /// from the previous mask's geometry. + /// from the previous mask's geometry: the trapezoid `½(αⁿ + αⁿ⁺¹)`, + /// or with `inner` intermediate geometries the composite trapezoid + /// over the step (the exact time integral of the space-time cut cell, + /// arXiv 2512.23358, approached as the sub-sampling refines). pub fn set_step_apertures(&mut self, old: &Mask) { + self.set_step_apertures_with(old, &[]); + } + + /// As [`Self::set_step_apertures`] with the apertures of the + /// intermediate geometries `inner` (in time order) inside the step. + pub fn set_step_apertures_with(&mut self, old: &Mask, inner: &[&CutGeometry]) { let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { return; }; - let avg = |a: &[f64], b: &[f64]| -> Vec { - a.iter().zip(b).map(|(x, y)| 0.5 * (x + y)).collect() + let n = inner.len() + 1; + let w_end = 0.5 / n as f64; + let w_in = 1.0 / n as f64; + let avg = |pick: &dyn Fn(&CutGeometry) -> &[f64]| -> Vec { + let a = pick(cut); + let b = pick(old_cut); + let mut out: Vec = a.iter().zip(b).map(|(x, y)| w_end * (x + y)).collect(); + for g in inner { + for (o, v) in out.iter_mut().zip(pick(g)) { + *o += w_in * v; + } + } + out }; - let au = avg(&cut.a_u, &old_cut.a_u); - let av = avg(&cut.a_v, &old_cut.a_v); - let aw = avg(&cut.a_w, &old_cut.a_w); + let au = avg(&|g: &CutGeometry| &g.a_u); + let av = avg(&|g: &CutGeometry| &g.a_v); + let aw = avg(&|g: &CutGeometry| &g.a_w); let open = |a: &[f64]| -> Vec { a.iter().map(|&x| x > 0.0).collect() }; let active = self .cell_fluid 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 50d5176..864c4c7 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 @@ -213,7 +213,13 @@ impl Solver { let a_w = (cv.wall[0] * cv.wall[0] + cv.wall[1] * cv.wall[1] + cv.wall[2] * cv.wall[2]).sqrt(); let shear = mu * a_w / cv.distance; - let v_eff = cv.alpha.max(INERTIA_FLOOR) * h[c] * area[c]; + let fraction = if self.params.momentum_volume_cell_mean { + let vol = |q: [i64; 3]| lat.cell(q).map_or(cv.alpha, |ci| mask.vol(ci)); + 0.5 * (vol(cell_minus) + vol(cell_plus)) + } else { + cv.alpha + }; + let v_eff = fraction.max(INERTIA_FLOOR) * h[c] * area[c]; let inertia = rho * v_eff / dt; (inertia * u0 - conv + diff + pressure + source + shear * ub) / (inertia + shear) } 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 6acf7c4..ffcb5c5 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 @@ -12,6 +12,7 @@ mod projection; use super::Grid; use super::body::Body; +use super::cut::CutGeometry; use super::field::Field; use super::poisson::{PcgCache, Problem, solve_pcg_cached}; use super::wall::{FaceKind, Mask, WallScheme}; @@ -79,6 +80,16 @@ pub struct Parameters { /// (identical results; the far corners keep their sign). `None` = /// every corner every step. pub max_surface_speed: Option, + /// Intermediate geometries per step for the space-time apertures + /// (0 = the trapezoid of the two end points). + pub aperture_substeps: usize, + /// The cut predictor's momentum volume per unit `h A`: the face's own + /// aperture (the registered design, `false`) or the mean fluid + /// fraction of the two cells the face separates (`true`: the control + /// volume's actual fluid content — a fresh face then carries the + /// inertia of the fluid half of its volume instead of a tenth of a + /// cell against O(1) fluxes). + pub momentum_volume_cell_mean: bool, } impl Default for Parameters { @@ -93,6 +104,8 @@ impl Default for Parameters { inner_stop_factor: 1e-2, wall_scheme: WallScheme::GhostBinary, max_surface_speed: None, + aperture_substeps: 0, + momentum_volume_cell_mean: false, } } } @@ -502,7 +515,32 @@ impl Solver { let mut new_mask = self.build_mask(body, field.grid, t_new, dt); if let Some(old_mask) = &self.mask { fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); - new_mask.set_step_apertures(old_mask); + let n_in = self.params.aperture_substeps; + if n_in == 0 { + new_mask.set_step_apertures(old_mask); + } else { + // The intermediate geometries, each within the band of the + // previous one. + let g = field.grid; + let h = g.dx.min(g.dy).min(g.dz); + let t_old = t_new - dt; + let mut inner: Vec = Vec::with_capacity(n_in); + for m in 1..=n_in { + let tm = t_old + dt * m as f64 / (n_in + 1) as f64; + let prev: Option<(&CutGeometry, f64, f64)> = match ( + self.params.max_surface_speed, + inner.last().or(old_mask.cut()), + ) { + (Some(speed), Some(c)) => { + Some((c, 3.0 * h, speed * dt / (n_in + 1) as f64)) + } + _ => None, + }; + inner.push(CutGeometry::build_from(body, g, tm, prev)); + } + let refs: Vec<&CutGeometry> = inner.iter().collect(); + new_mask.set_step_apertures_with(old_mask, &refs); + } } new_mask.impose_from( body, diff --git a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs index a1d6255..8b50a72 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs @@ -241,6 +241,14 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { z1: Side::Periodic, ..Boundaries::default() }, + max_surface_speed: Some(U_PEAK), + // `RTX_E3_SUBSTEPS`: intermediate geometries for the space-time apertures. + aperture_substeps: std::env::var("RTX_E3_SUBSTEPS") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(0), + // `RTX_E3_CELL_MEAN=1`: the momentum volume from the two cells' fractions. + momentum_volume_cell_mean: std::env::var("RTX_E3_CELL_MEAN").is_ok_and(|v| v == "1"), ..Parameters::default() }, ); diff --git a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs index 841f380..fdd4a6b 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_sphere/mod.rs @@ -133,6 +133,7 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement tolerance: 1e-8, convection_scheme: ConvectionScheme::Upwind, wall_scheme: scheme, + momentum_volume_cell_mean: std::env::var("RTX_E3_CELL_MEAN").is_ok_and(|v| v == "1"), ..Parameters::default() }, );