diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs index 56895ea..37816c9 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs @@ -65,6 +65,14 @@ unsafe impl ValidAsZeroBits for E3CutPtrs {} pub(super) struct DeviceCut { /// The surface-velocity tables (kept to share between the phases). ub_host: [Vec; 3], + /// Which phase's apertures / open flags / activity `ptrs` hands out + /// (P1-5: both sets live in the struct and persist across steps). + phase: Phase, + /// The predictor phase's instantaneous apertures, open flags (fluid + /// kind) and activity (fluid cells). + a_pred: [CudaSlice; 3], + open_pred: [CudaSlice; 3], + active_pred: CudaSlice, a: [CudaSlice; 3], d: [CudaSlice; 3], ub: [CudaSlice; 3], @@ -307,8 +315,43 @@ impl DeviceCut { let l_cells = lap.elapsed(); let ub_host = ub; let ub_dev = [up_f(&ub_host[0]), up_f(&ub_host[1]), up_f(&ub_host[2])]; + // The predictor phase's set (instantaneous apertures, fluid kinds, + // fluid cells) — built alongside so that both persist (P1-5). + let open_pred = |c: usize| -> Vec { + (0..counts[c]) + .into_par_iter() + .map(|idx| { + let kind = match c { + 0 => mask.u_kind(idx), + 1 => mask.v_kind(idx), + _ => mask.w_kind(idx), + }; + i32::from(kind == FaceKind::Fluid) + }) + .collect() + }; + let active_pred: Vec = (0..nc).into_par_iter().map(|i| i32::from(mask.is_fluid_cell(i))).collect(); + let (a_pred, open_pred_dev, active_pred_dev) = if projection { + ( + [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)], + [up_i(&open_pred(0)), up_i(&open_pred(1)), up_i(&open_pred(2))], + up_i(&active_pred), + ) + } else { + // A predictor-phase build IS the predictor set; the projection + // set is filled with the same arrays (replaced on the next update). + ( + [up_f(ap_u), up_f(ap_v), up_f(ap_w)], + [up_i(&open[0]), up_i(&open[1]), up_i(&open[2])], + up_i(&active), + ) + }; let built = Self { ub_host, + phase, + a_pred, + open_pred: open_pred_dev, + active_pred: active_pred_dev, a: [up_f(ap_u), up_f(ap_v), up_f(ap_w)], d: [up_f(&cut.d_u), up_f(&cut.d_v), up_f(&cut.d_w)], ub: ub_dev, @@ -347,42 +390,230 @@ impl DeviceCut { /// activity flags and the apertures (end-of-step instead of the step's) /// are recomputed and uploaded. Identical to a fresh `build_with` in /// every table. - pub(super) fn predictor_from(prev: Self, solver: &Solver, g: Grid) -> Option { + pub(super) fn predictor_from(mut prev: Self, solver: &Solver, g: Grid) -> Option { + // P1-5: the predictor set was built or updated with the projection's; + // switching phases is a flag. + let _ = (solver, g); + prev.phase = Phase::Predictor; + Some(prev) + } + + /// P1-5: update the persistent tables in place for the projection phase + /// of the next step of a MOVING body: the per-face tables for the faces + /// of the changed cells and the surface velocities for the imposition + /// band's faces, the per-cell tables for the changed cells (the wall + /// flux for the wall cells too: its compatibility correction is + /// global). `false` when a full build is needed (no changed set, the + /// advancing-wall closure, a size change). + pub(super) fn update(&mut self, solver: &Solver, g: Grid, t: f64) -> bool { + use crate::solvers::incompressible::embedded3::poisson::device_cg::{scatter_f64, scatter_i32, scatter_u32}; use rayon::prelude::*; - let mask = solver.mask()?; - let cut = mask.cut()?; - let rt = runtime(); - let up_f = |v: &[f64]| -> CudaSlice { - rt.stream - .memcpy_stod(if v.is_empty() { &[0.0f64][..] } else { v }) - .expect("upload") - }; - let up_i = |v: &[i32]| -> CudaSlice { rt.stream.memcpy_stod(v).expect("upload") }; + let Some(mask) = solver.mask() else { return false }; + let Some(cut) = mask.cut() else { return false }; + let Some(body) = solver.body() else { return false }; + let Some(changed) = mask.changed_cells() else { return false }; + if mask.wall_advancing || mask.wall_exchange_foot { + return false; + } let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let nxy = nx * ny; + let nc = g.cells(); let counts = [(nx + 1) * ny * nz, nx * (ny + 1) * nz, nx * ny * (nz + 1)]; - let open = |c: usize| -> Vec { - (0..counts[c]) - .into_par_iter() - .map(|idx| { + if self.active.len() != nc || self.a[0].len() != counts[0] { + return false; + } + let lap = Instant::now(); + let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); + // The faces of the changed cells, per component, ascending. + let mut marks: [Vec; 3] = [vec![false; counts[0]], vec![false; counts[1]], vec![false; counts[2]]]; + for &idx in changed { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + marks[0][g.uface(k, j, i)] = true; + marks[0][g.uface(k, j, i + 1)] = true; + marks[1][g.vface(k, j, i)] = true; + marks[1][g.vface(k, j + 1, i)] = true; + marks[2][g.wface(k, j, i)] = true; + marks[2][g.wface(k + 1, j, i)] = true; + } + let touched: [Vec; 3] = [0, 1, 2].map(|c| { + (0..counts[c]).into_par_iter().filter(|&f| marks[c][f]).map(|f| f as u32).collect() + }); + let cells_u32: Vec = changed.iter().map(|&i| i as u32).collect(); + // The imposition band's faces (their surface velocity moves with t). + let band = mask.impose_band().unwrap_or(f64::INFINITY); + let dists: [&[f64]; 3] = [&cut.d_u, &cut.d_v, &cut.d_w]; + let band_faces: [Vec; 3] = [0, 1, 2].map(|c| { + (0..counts[c]).into_par_iter().filter(|&f| dists[c][f].abs() <= band).map(|f| f as u32).collect() + }); + let h = [g.dx, g.dy, g.dz]; + let shifts = mask.face_shift_tables(); + let step = mask.step_apertures(); + let step_open = mask.step_open_flags(); + for c in 0..3 { + let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1)); + let pos = |idx: usize| -> [f64; 3] { + let (k, j, i) = (idx / (nj * ni), (idx / ni) % nj, idx % ni); + [ + (i as f64 + if c == 0 { 0.0 } else { 0.5 }) * h[0], + (j as f64 + if c == 1 { 0.0 } else { 0.5 }) * h[1], + (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2], + ] + }; + let a_inst: &[f64] = match c { 0 => &cut.a_u, 1 => &cut.a_v, _ => &cut.a_w }; + let a_step: &[f64] = match (step, c) { + (Some(a), 0) => &a.0, + (Some(a), 1) => &a.1, + (Some(a), _) => &a.2, + (None, _) => a_inst, + }; + let tf = &touched[c]; + let pick = |v: &[f64]| tf.iter().map(|&f| v[f as usize]).collect::>(); + scatter_f64(tf, &pick(a_step), &mut self.a[c]); + scatter_f64(tf, &pick(a_inst), &mut self.a_pred[c]); + scatter_f64(tf, &pick(dists[c]), &mut self.d[c]); + let open_step: Vec = tf.iter().map(|&f| i32::from(step_open.map_or(a_inst[f as usize] > 0.0, |o| o[c][f as usize]))).collect(); + scatter_i32(tf, &open_step, &mut self.open[c]); + let open_inst: Vec = tf + .iter() + .map(|&f| { let kind = match c { - 0 => mask.u_kind(idx), - 1 => mask.v_kind(idx), - _ => mask.w_kind(idx), + 0 => mask.u_kind(f as usize), + 1 => mask.v_kind(f as usize), + _ => mask.w_kind(f as usize), }; i32::from(kind == FaceKind::Fluid) }) - .collect() - }; - let active: Vec = (0..g.cells()) + .collect(); + scatter_i32(tf, &open_inst, &mut self.open_pred[c]); + if let Some(sh) = shifts { + let idx3: Vec = tf.iter().flat_map(|&f| [3 * f, 3 * f + 1, 3 * f + 2]).collect(); + let val3: Vec = tf.iter().flat_map(|&f| [sh[c][3 * f as usize], sh[c][3 * f as usize + 1], sh[c][3 * f as usize + 2]]).collect(); + scatter_f64(&idx3, &val3, &mut self.shift[c]); + } + // The surface velocity at the foot for the band's faces. + let bf = &band_faces[c]; + let ubv: Vec = bf + .par_iter() + .map(|&f| { + let x = pos(f as usize); + let xf = match (shifts, mask.wall_foot_centroid) { + (Some(sh), true) => [x[0] + sh[c][3 * f as usize], x[1] + sh[c][3 * f as usize + 1], x[2] + sh[c][3 * f as usize + 2]], + _ => x, + }; + mask.surface_velocity_at(body, xf, c, t) + }) + .collect(); + for (&f, &v) in bf.iter().zip(&ubv) { + self.ub_host[c][f as usize] = v; + } + scatter_f64(bf, &ubv, &mut self.ub[c]); + } + let l_faces = lap.elapsed(); + // Per-cell tables for the changed cells. + let active_step: Vec = changed.iter().map(|&i| i32::from(mask.cell_active(i))).collect(); + let active_inst: Vec = changed.iter().map(|&i| i32::from(mask.is_fluid_cell(i))).collect(); + let owner: Vec = changed.iter().map(|&i| mask.master(i).unwrap_or(i) as u32).collect(); + scatter_i32(&cells_u32, &active_step, &mut self.active); + scatter_i32(&cells_u32, &active_inst, &mut self.active_pred); + scatter_u32(&cells_u32, &owner, &mut self.owner); + // The wall flux: the changed cells and every wall cell (the correction). + let wall_flux: &[f64] = solver.wall_fluxes(); + if wall_flux.len() != nc { + return false; + } + let mut wmark = vec![false; nc]; + for &i in changed { + wmark[i] = true; + } + let wcells: Vec = (0..nc) .into_par_iter() - .map(|i| i32::from(mask.is_fluid_cell(i))) + .filter(|&i| wmark[i] || cut.wall[i] != [0.0; 3]) + .map(|i| i as u32) .collect(); - Some(Self { - a: [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)], - open: [up_i(&open(0)), up_i(&open(1)), up_i(&open(2))], - active: up_i(&active), - ..prev - }) + let wvals: Vec = wcells.iter().map(|&i| wall_flux[i as usize]).collect(); + scatter_f64(&wcells, &wvals, &mut self.wall_flux); + // The merged-cell CSR whole (its lengths change). + let mut fold_ptr = vec![0u32; nc + 1]; + let mut merged = 0; + for i in 0..nc { + if let Some(m) = mask.master(i) { + fold_ptr[m + 1] += 1; + merged += 1; + } + } + for i in 0..nc { + fold_ptr[i + 1] += fold_ptr[i]; + } + let mut fold_idx = vec![0u32; merged]; + let mut cursor: Vec = fold_ptr[..nc].to_vec(); + for i in 0..nc { + if let Some(m) = mask.master(i) { + fold_idx[cursor[m] as usize] = i as u32; + cursor[m] += 1; + } + } + let rt = runtime(); + let up_u = |v: &[u32]| -> CudaSlice { + rt.stream + .memcpy_stod(if v.is_empty() { &[0u32][..] } else { v }) + .expect("upload") + }; + self.fold_ptr = up_u(&fold_ptr); + self.fold_idx = up_u(&fold_idx); + self.merged = merged; + self.phase = Phase::Projection; + rt.stream.synchronize().expect("sync"); + if profile { + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + eprintln!( + " table laps (update): faces {:.0} ms ({} / {} / {} touched, {} / {} / {} band), cells {:.0} ms ({} changed)", + ms(l_faces), + touched[0].len(), + touched[1].len(), + touched[2].len(), + band_faces[0].len(), + band_faces[1].len(), + band_faces[2].len(), + ms(lap.elapsed() - l_faces), + changed.len() + ); + } + true + } + + /// `RTX_E3_BAND_CHECK=1`: every device table against a fresh full build. + pub(super) fn check_against(&self, full: &Self) { + let rt = runtime(); + let same_f = |name: &str, a: &CudaSlice, b: &CudaSlice| { + let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); + assert!(x.len() == y.len() && x.iter().zip(&y).all(|(p, q)| p.to_bits() == q.to_bits()), "band tables: {name} differs"); + }; + let same_i = |name: &str, a: &CudaSlice, b: &CudaSlice| { + let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); + assert!(x == y, "band tables: {name} differs"); + }; + let same_u = |name: &str, a: &CudaSlice, b: &CudaSlice| { + let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); + assert!(x == y, "band tables: {name} differs"); + }; + for c in 0..3 { + same_f("a", &self.a[c], &full.a[c]); + same_f("a_pred", &self.a_pred[c], &full.a_pred[c]); + same_f("d", &self.d[c], &full.d[c]); + same_i("open", &self.open[c], &full.open[c]); + same_i("open_pred", &self.open_pred[c], &full.open_pred[c]); + same_f("shift", &self.shift[c], &full.shift[c]); + // ub: the band's faces only (beyond it the full build writes 0, ours keeps stale — never read). + let (x, y) = (rt.stream.memcpy_dtov(&self.ub[c]).expect("dl"), rt.stream.memcpy_dtov(&full.ub[c]).expect("dl")); + let bad = x.iter().zip(&y).filter(|(p, q)| **q != 0.0 && p.to_bits() != q.to_bits()).count(); + assert!(bad == 0, "band tables: ub[{c}] differs on {bad} band faces"); + } + same_f("wall_flux", &self.wall_flux, &full.wall_flux); + same_i("active", &self.active, &full.active); + same_i("active_pred", &self.active_pred, &full.active_pred); + same_u("owner", &self.owner, &full.owner); + same_u("fold_ptr", &self.fold_ptr, &full.fold_ptr); + same_u("fold_idx", &self.fold_idx, &full.fold_idx); } fn ptrs(&self) -> E3CutPtrs { @@ -391,11 +622,15 @@ impl DeviceCut { let pf = |x: &CudaSlice| x.device_ptr(s).0; let pi = |x: &CudaSlice| x.device_ptr(s).0; let pu = |x: &CudaSlice| x.device_ptr(s).0; + let (a, open, active) = match self.phase { + Phase::Projection => (&self.a, &self.open, &self.active), + Phase::Predictor => (&self.a_pred, &self.open_pred, &self.active_pred), + }; E3CutPtrs { ptrs: [ - pf(&self.a[0]), - pf(&self.a[1]), - pf(&self.a[2]), + pf(&a[0]), + pf(&a[1]), + pf(&a[2]), pf(&self.d[0]), pf(&self.d[1]), pf(&self.d[2]), @@ -403,10 +638,10 @@ impl DeviceCut { pf(&self.ub[1]), pf(&self.ub[2]), pf(&self.wall_flux), - pi(&self.open[0]), - pi(&self.open[1]), - pi(&self.open[2]), - pi(&self.active), + pi(&open[0]), + pi(&open[1]), + pi(&open[2]), + pi(active), pu(&self.owner), pu(&self.fold_ptr), pu(&self.fold_idx), @@ -501,7 +736,25 @@ impl DeviceStep { self.upload_after_rebuild(&field); self.host_field = Some(field); let l_up = lap.elapsed(); - self.cut = DeviceCut::build(&self.solver, g, Phase::Projection, t_new); + // P1-5: the persistent tables updated in place when the mask has a + // changed set; a full build otherwise (and under the identity + // reference `RTX_E3_TABLES_REBUILD=1`). + let rebuild_all = std::env::var("RTX_E3_TABLES_REBUILD").is_ok(); + let check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); + self.cut = match self.cut.take() { + Some(mut prev) if !rebuild_all => { + if prev.update(&self.solver, g, t_new) { + if check { + let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new).expect("full"); + prev.check_against(&full); + } + Some(prev) + } else { + DeviceCut::build(&self.solver, g, Phase::Projection, t_new) + } + } + _ => DeviceCut::build(&self.solver, g, Phase::Projection, t_new), + }; let l_tables = lap.elapsed(); if profile { eprintln!( 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 0ec01dd..33c52b1 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs @@ -578,6 +578,12 @@ impl Mask { _ => None, } } + /// The space-time open flags per u / v / w face (a moving cut wall). + #[must_use] + pub fn step_open_flags(&self) -> Option<[&[bool]; 3]> { + self.step_open.as_ref().map(|o| [o.0.as_slice(), o.1.as_slice(), o.2.as_slice()]) + } + /// The cells whose operator rows may have changed since the previous /// step (P1-4), or `None` when every row must be rebuilt. #[must_use]