diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs index ba192a2..4760802 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs @@ -37,6 +37,8 @@ pub struct CutGeometry { /// narrow band) — a cell whose corners are all untouched has unchanged /// apertures, volume and activity. pub touched: Vec, + /// Per cell: any of its eight corners touched (computed once per build). + pub touched_cell: Vec, } impl CutGeometry { @@ -191,6 +193,21 @@ impl CutGeometry { let sz = (a_w[g.wface(k + 1, j, i)] - a_w[g.wface(k, j, i)]) * az; *w = [-sx, -sy, -sz]; }); + let touched_cell: Vec = (0..g.cells()) + .into_par_iter() + .map(|idx| { + let (k, j, i) = (idx / (ny * nx), (idx / nx) % ny, idx % nx); + let mut t = false; + for dk in 0..2 { + for dj in 0..2 { + for di in 0..2 { + t |= touched[Self::node(g, k + dk, j + dj, i + di)]; + } + } + } + t + }) + .collect(); Self { grid, phi, @@ -204,12 +221,16 @@ impl CutGeometry { d_w, bound, touched, + touched_cell, } } - /// The cells with a touched corner (P1-4). + /// The cells with a touched corner (P1-4; the stored map). #[must_use] pub fn touched_cells(&self) -> Vec { + if !self.touched_cell.is_empty() { + return self.touched_cell.clone(); + } use rayon::prelude::*; let g = self.grid; let (nx, ny) = (g.nx, g.ny); 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 38d627e..ad7b385 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -283,73 +283,99 @@ impl Mask { /// 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) { + pub fn set_step_apertures(&mut self, old: &mut 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]) { + pub fn set_step_apertures_with(&mut self, old: &mut Mask, inner: &[&CutGeometry]) { let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { return; }; 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(&|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 - .iter() - .zip(&old.cell_fluid) - .map(|(&n, &o)| n || o) - .collect(); - self.step_open = Some((open(&au), open(&av), open(&aw), active)); - self.step_apertures = Some((au, av, aw)); + use rayon::prelude::*; + let g = self.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let nxy = nx * ny; // P1-4: the changed set — cells touched by either build (this step's // apertures average the two geometries), dilated by one. - let t_new = cut.touched_cells(); - let t_old = old_cut.touched_cells(); - self.compute_merging(Some(old)); - { - use rayon::prelude::*; - let g = self.grid; - let (nx, ny, nz) = (g.nx, g.ny, g.nz); - let nxy = nx * ny; - let periodic = self.periodic_z; - let changed: Vec = (0..g.cells()) - .into_par_iter() - .filter(|&idx| { - let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); - let t = |q: usize| t_new[q] || t_old[q]; - if t(idx) { - return true; + let t_new = &cut.touched_cell; + let t_old = &old_cut.touched_cell; + let periodic = self.periodic_z; + let changed: Vec = (0..g.cells()) + .into_par_iter() + .filter(|&idx| { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + let t = |q: usize| t_new[q] || t_old[q]; + if t(idx) { + return true; + } + (i + 1 < nx && t(idx + 1)) + || (i > 0 && t(idx - 1)) + || (j + 1 < ny && t(idx + nx)) + || (j > 0 && t(idx - nx)) + || (k + 1 < nz && t(idx + nxy)) + || (k > 0 && t(idx - nxy)) + || (periodic && nz > 1 && k + 1 == nz && t(idx - (nz - 1) * nxy)) + || (periodic && nz > 1 && k == 0 && t(idx + (nz - 1) * nxy)) + }) + .collect(); + // P1-5 (j): with the trapezoid alone and a previous step's arrays, a + // face of no changed cell keeps its step aperture (its corners were + // untouched in both builds, so αⁿ⁻¹ = αⁿ = αⁿ⁺¹): the old mask's + // arrays move over and only the changed cells' faces and the + // changed cells' activity are recomputed, with the same expressions. + let prev = if inner.is_empty() { old.step_apertures.take().zip(old.step_open.take()) } else { None }; + if let Some(((mut au, mut av, mut aw), (mut ou, mut ov, mut ow, mut active))) = prev { + let mix = |a: &[f64], b: &[f64], f: usize| w_end * (a[f] + b[f]); + for &idx in &changed { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + for f in [g.uface(k, j, i), g.uface(k, j, i + 1)] { + au[f] = mix(&cut.a_u, &old_cut.a_u, f); + ou[f] = au[f] > 0.0; + } + for f in [g.vface(k, j, i), g.vface(k, j + 1, i)] { + av[f] = mix(&cut.a_v, &old_cut.a_v, f); + ov[f] = av[f] > 0.0; + } + for f in [g.wface(k, j, i), g.wface(k + 1, j, i)] { + aw[f] = mix(&cut.a_w, &old_cut.a_w, f); + ow[f] = aw[f] > 0.0; + } + active[idx] = self.cell_fluid[idx] || old.cell_fluid[idx]; + } + self.step_open = Some((ou, ov, ow, active)); + self.step_apertures = Some((au, av, aw)); + } else { + 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 gi in inner { + for (o, v) in out.iter_mut().zip(pick(gi)) { + *o += w_in * v; } - (i + 1 < nx && t(idx + 1)) - || (i > 0 && t(idx - 1)) - || (j + 1 < ny && t(idx + nx)) - || (j > 0 && t(idx - nx)) - || (k + 1 < nz && t(idx + nxy)) - || (k > 0 && t(idx - nxy)) - || (periodic && nz > 1 && k + 1 == nz && t(idx - (nz - 1) * nxy)) - || (periodic && nz > 1 && k == 0 && t(idx + (nz - 1) * nxy)) - }) + } + out + }; + let au = avg(&|gg: &CutGeometry| &gg.a_u); + let av = avg(&|gg: &CutGeometry| &gg.a_v); + let aw = avg(&|gg: &CutGeometry| &gg.a_w); + let open = |a: &[f64]| -> Vec { a.iter().map(|&x| x > 0.0).collect() }; + let active = self + .cell_fluid + .iter() + .zip(&old.cell_fluid) + .map(|(&n, &o)| n || o) .collect(); - self.changed_cells = Some(changed); + self.step_open = Some((open(&au), open(&av), open(&aw), active)); + self.step_apertures = Some((au, av, aw)); } + self.compute_merging(Some(old)); + self.changed_cells = Some(changed); } pub(super) fn lattice(&self) -> Lattice { diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs index 43931c5..7b6e4ac 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs @@ -21,7 +21,7 @@ impl Solver { let lap = std::time::Instant::now(); let mut new_mask = self.build_mask(body, field.grid, t_new, dt); let l_build = lap.elapsed(); - if let Some(old_mask) = &self.mask { + if let Some(old_mask) = self.mask.as_mut() { fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); let n_in = self.params.aperture_substeps; if n_in == 0 {