diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs index 5ba3a2c..fb14cf0 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/poisson/export.rs @@ -130,13 +130,15 @@ pub fn export_fine(problem: &Problem) -> LevelExport { /// is the identical export without a second clone of the problem and a /// second level build. pub(crate) fn export_fine_from(lv: &Level) -> LevelExport { + use rayon::prelude::*; let (_, coarse_of) = lv.coarsen(); + // Per-entry maps (P1-3 (e)): the same values in the same order. let to_u32 = |v: &[usize]| { - v.iter() + v.par_iter() .map(|&i| if i == usize::MAX { u32::MAX } else { i as u32 }) .collect::>() }; - let cast = |v: &[f64]| v.iter().map(|&x| x as f32).collect::>(); + let cast = |v: &[f64]| v.par_iter().map(|&x| x as f32).collect::>(); LevelExport { nx: lv.problem.nx, ny: lv.problem.ny, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs index 2beb8f2..622990f 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs @@ -151,6 +151,9 @@ pub struct DeviceStep { n_blocks: usize, cg: Option, cg_dt: f64, + /// The host mirror a moving body's rebuild reads and writes, kept + /// between steps (P1-3 (c)). + host_field: Option, /// PERF-3 P2: per corrector index, the basis of the projected initial guess. pub(super) guess: Vec, timers: Option, @@ -231,6 +234,7 @@ impl DeviceStep { n_blocks, cg: None, cg_dt: 0.0, + host_field: None, guess: Vec::new(), timers, initialized: false, @@ -266,6 +270,40 @@ impl DeviceStep { rt.stream.synchronize().expect("sync"); } + /// The arrays a moving body's mask rebuild writes (P1-3 (c)): the + /// velocities (imposed and refilled) and the pressure (refilled) — the + /// correction stays on the device, untouched by the rebuild. + pub fn upload_after_rebuild(&mut self, field: &Field) { + let rt = runtime(); + assert_eq!(field.grid, self.grid); + rt.stream.memcpy_htod(&field.u, &mut self.u).expect("u"); + rt.stream.memcpy_htod(&field.v, &mut self.v).expect("v"); + rt.stream.memcpy_htod(&field.w, &mut self.w).expect("w"); + rt.stream.memcpy_htod(&field.p, &mut self.p).expect("p"); + rt.stream.synchronize().expect("sync"); + } + + /// The arrays a moving body's mask rebuild reads (P1-3 (c)): the + /// corrected velocities, the pressure and the previous velocities. + pub fn download_for_rebuild(&self, field: &mut Field) { + let rt = runtime(); + assert_eq!(field.grid, self.grid); + rt.stream.memcpy_dtoh(&self.u, &mut field.u).expect("u"); + rt.stream.memcpy_dtoh(&self.v, &mut field.v).expect("v"); + rt.stream.memcpy_dtoh(&self.w, &mut field.w).expect("w"); + rt.stream.memcpy_dtoh(&self.p, &mut field.p).expect("p"); + rt.stream + .memcpy_dtoh(&self.u_old, &mut field.u_old) + .expect("u_old"); + rt.stream + .memcpy_dtoh(&self.v_old, &mut field.v_old) + .expect("v_old"); + rt.stream + .memcpy_dtoh(&self.w_old, &mut field.w_old) + .expect("w_old"); + rt.stream.synchronize().expect("sync"); + } + /// The device field into the host mirror. pub fn download(&self, field: &mut Field) { let rt = runtime(); 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 cd85aeb..4ba9126 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 @@ -489,12 +489,17 @@ impl DeviceStep { // `RTX_E3_MOVING_PROFILE=1`: the laps of the moving path (S2-2b). let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); let lap = Instant::now(); - let mut field = Field::new(g); - self.download(&mut field); + // P1-3 (c): the host mirror persists (no 1.1 GB allocation per + // step) and only the arrays the rebuild reads come down — u, v, + // w, p (fresh cells take their neighbours' pressure) and the + // previous velocities (the imposition's source). + let mut field = self.host_field.take().unwrap_or_else(|| Field::new(g)); + self.download_for_rebuild(&mut field); let l_down = lap.elapsed(); fresh_cells = self.solver.rebuild_moving_mask(&mut field, dt, t_new); let l_mask = lap.elapsed(); - self.upload(&field); + 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); let l_tables = lap.elapsed(); 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 index a1f48fe..c58632e 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs @@ -25,12 +25,48 @@ impl Solver { 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 { + // P1-3 (d): every entry depends on its own cell alone — assembled in + // parallel by plane with the serial loop's arithmetic. + use rayon::prelude::*; + let nxy = nx * ny; + // The predicates over the mask alone (the solver holds `Cell`s and + // cannot cross threads; these are its `u_is_unknown` etc. verbatim). + let mask = self.mask.as_ref(); + let cell_is_active = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.cell_active(g.cell(k, j, i))); + let u_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.u_open(g.uface(k, j, i))); + let v_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.v_open(g.vface(k, j, i))); + let w_is_unknown = |k: usize, j: usize, i: usize| mask.is_none_or(|m| m.w_open(g.wface(k, j, i))); + let au = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.au_step(g.uface(k, j, i))); + let av = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.av_step(g.vface(k, j, i))); + let aw = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.aw_step(g.wface(k, j, i))); + let gu = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.grad_weight(0, g.uface(k, j, i))); + let gv = |k: usize, j: usize, i: usize| mask.map_or(1.0, |m| m.grad_weight(1, g.vface(k, j, i))); + let (mut p_active, mut p_ae, mut p_aw, mut p_an, mut p_as, mut p_at, mut p_ab, mut p_extra) = ( + std::mem::take(&mut problem.active), + std::mem::take(&mut problem.ae), + std::mem::take(&mut problem.aw), + std::mem::take(&mut problem.an), + std::mem::take(&mut problem.as_), + std::mem::take(&mut problem.at), + std::mem::take(&mut problem.ab), + std::mem::take(&mut problem.extra_diag), + ); + p_active + .par_chunks_mut(nxy) + .zip(p_ae.par_chunks_mut(nxy)) + .zip(p_aw.par_chunks_mut(nxy)) + .zip(p_an.par_chunks_mut(nxy)) + .zip(p_as.par_chunks_mut(nxy)) + .zip(p_at.par_chunks_mut(nxy)) + .zip(p_ab.par_chunks_mut(nxy)) + .zip(p_extra.par_chunks_mut(nxy)) + .enumerate() + .for_each(|(k, (((((((c_active, c_ae), c_aw), c_an), c_as), c_at), c_ab), c_extra))| { for j in 0..ny { for i in 0..nx { - let idx = g.cell(k, j, i); - if !self.cell_is_active(k, j, i) { - problem.active[idx] = false; + let idx = j * nx + i; + if !cell_is_active(k, j, i) { + c_active[idx] = false; continue; } let mut extra = 0.0; @@ -38,48 +74,56 @@ impl Solver { if b.x1 == outlet { extra += ae_outlet; } - } else if self.u_is_unknown(k, j, i + 1) { - problem.ae[idx] = ae_interior * self.au(k, j, i + 1) * self.gu(k, j, i + 1); + } else if u_is_unknown(k, j, i + 1) { + c_ae[idx] = ae_interior * au(k, j, i + 1) * gu(k, j, i + 1); } if i == 0 { if b.x0 == outlet { extra += ae_outlet; } - } else if self.u_is_unknown(k, j, i) { - problem.aw[idx] = ae_interior * self.au(k, j, i) * self.gu(k, j, i); + } else if u_is_unknown(k, j, i) { + c_aw[idx] = ae_interior * au(k, j, i) * gu(k, j, i); } if j + 1 == ny { if b.y1 == outlet { extra += an_outlet; } - } else if self.v_is_unknown(k, j + 1, i) { - problem.an[idx] = an_interior * self.av(k, j + 1, i) * self.gv(k, j + 1, i); + } else if v_is_unknown(k, j + 1, i) { + c_an[idx] = an_interior * av(k, j + 1, i) * gv(k, j + 1, i); } if j == 0 { if b.y0 == outlet { extra += an_outlet; } - } else if self.v_is_unknown(k, j, i) { - problem.as_[idx] = an_interior * self.av(k, j, i) * self.gv(k, j, i); + } else if v_is_unknown(k, j, i) { + c_as[idx] = an_interior * av(k, j, i) * gv(k, j, i); } if k + 1 == nz && !periodic { if b.z1 == outlet { extra += at_outlet; } - } else if self.w_is_unknown((k + 1) % nz, j, i) { - problem.at[idx] = at_interior * self.aw((k + 1) % nz, j, i); + } else if w_is_unknown((k + 1) % nz, j, i) { + c_at[idx] = at_interior * aw((k + 1) % nz, j, i); } if k == 0 && !periodic { if b.z0 == outlet { extra += at_outlet; } - } else if self.w_is_unknown(k, j, i) { - problem.ab[idx] = at_interior * self.aw(k, j, i); + } else if w_is_unknown(k, j, i) { + c_ab[idx] = at_interior * aw(k, j, i); } - problem.extra_diag[idx] = extra; + c_extra[idx] = extra; } } - } + }); + problem.active = p_active; + problem.ae = p_ae; + problem.aw = p_aw; + problem.an = p_an; + problem.as_ = p_as; + problem.at = p_at; + problem.ab = p_ab; + problem.extra_diag = p_extra; self.merge_small_cells(&mut problem); problem }