diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs index 0ec824c..0679490 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs @@ -48,42 +48,38 @@ impl Mask { vec![0.0; 3 * sizes[1]], vec![0.0; 3 * sizes[2]], ]; + // P1-3 (f): per face in parallel (each entry from its own geometry). + use rayon::prelude::*; for c in 0..3 { - let (ni, nj, nk) = ( - g.nx + usize::from(c == 0), - g.ny + usize::from(c == 1), - g.nz + usize::from(c == 2), - ); - for k in 0..nk { - for j in 0..nj { - for i in 0..ni { - let p = [i as i64, j as i64, k as i64]; - let Some(f) = lat.face(c, p) else { continue }; - let Some(alpha) = self.aperture(c, p) else { - continue; - }; - if alpha <= 0.0 || alpha >= 1.0 { - continue; - } - // Interior faces only (a control volume needs both cells). - let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c]; - if on_side && !(c == 2 && self.periodic_z) { - continue; - } - let cv = self.cv_geometry(c, p); - let mut n = cv.wall; - n[c] = 0.0; - let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); - if a == 0.0 { - continue; - } - for d in 0..3 { - // `wall` points into the body: the open part lies the other way. - tables[c][3 * f + d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a; - } + let (ni, nj) = (g.nx + usize::from(c == 0), g.ny + usize::from(c == 1)); + tables[c] + .par_chunks_mut(3) + .enumerate() + .for_each(|(f, out)| { + let (k, j, i) = (f / (nj * ni), (f / ni) % nj, f % ni); + let p = [i as i64, j as i64, k as i64]; + debug_assert_eq!(lat.face(c, p), Some(f)); + let Some(alpha) = self.aperture(c, p) else { return }; + if alpha <= 0.0 || alpha >= 1.0 { + return; } - } - } + // Interior faces only (a control volume needs both cells). + let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c]; + if on_side && !(c == 2 && self.periodic_z) { + return; + } + let cv = self.cv_geometry(c, p); + let mut n = cv.wall; + n[c] = 0.0; + let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if a == 0.0 { + return; + } + for d in 0..3 { + // `wall` points into the body: the open part lies the other way. + out[d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a; + } + }); } self.face_shifts = Some(tables); } 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 0b49b42..f916c9c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -128,35 +128,31 @@ impl Mask { let (nx, ny, nz) = (g.nx, g.ny, g.nz); let periodic = b.z0 == Side::Periodic; let allowed = |side: Side| matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall); - let mut cell_fluid = vec![true; g.cells()]; - let mut fluid_cells = 0; - let mut anchor = None; - for k in 0..nz { - for j in 0..ny { - for i in 0..nx { - let idx = g.cell(k, j, i); - let fluid = cut.vol[idx] > 0.0; - cell_fluid[idx] = fluid; - if fluid { - fluid_cells += 1; - if anchor.is_none() { - anchor = Some(idx); - } - } else { - let touches = (i == 0 && !allowed(b.x0)) - || (i + 1 == nx && !allowed(b.x1)) - || (j == 0 && !allowed(b.y0)) - || (j + 1 == ny && !allowed(b.y1)) - || (k == 0 && !allowed(b.z0)) - || (k + 1 == nz && !allowed(b.z1)); - if touches { - return Err(format!( - "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})" - )); - } - } - } + // P1-3 (f): the whole-grid classifications as parallel per-entry + // maps (the same values; the anchor is the smallest fluid index, the + // serial loop's first). + use rayon::prelude::*; + let nxy = nx * ny; + let cell_fluid: Vec = cut.vol.par_iter().map(|&v| v > 0.0).collect(); + let fluid_cells = cell_fluid.par_iter().filter(|&&f| f).count(); + let anchor = (0..g.cells()).into_par_iter().find_first(|&idx| cell_fluid[idx]); + let touching = (0..g.cells()).into_par_iter().find_first(|&idx| { + if cell_fluid[idx] { + return false; } + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + (i == 0 && !allowed(b.x0)) + || (i + 1 == nx && !allowed(b.x1)) + || (j == 0 && !allowed(b.y0)) + || (j + 1 == ny && !allowed(b.y1)) + || (k == 0 && !allowed(b.z0)) + || (k + 1 == nz && !allowed(b.z1)) + }); + if let Some(idx) = touching { + let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); + return Err(format!( + "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})" + )); } let Some(anchor) = anchor else { return Err("embedded body covers the whole domain".into()); @@ -168,32 +164,28 @@ impl Mask { FaceKind::Solid } }; - let mut u_kind = vec![FaceKind::Fluid; g.n_ufaces()]; - let mut v_kind = vec![FaceKind::Fluid; g.n_vfaces()]; - let mut w_kind = vec![FaceKind::Fluid; g.n_wfaces()]; - for k in 0..nz { - for j in 0..ny { - for i in 1..nx { - let f = g.uface(k, j, i); - u_kind[f] = kind(cut.a_u[f]); - } - } - for j in 1..ny { - for i in 0..nx { - let f = g.vface(k, j, i); - v_kind[f] = kind(cut.a_v[f]); - } - } - } - let w_range = if periodic { 0..nz + 1 } else { 1..nz }; - for k in w_range { - for j in 0..ny { - for i in 0..nx { - let f = g.wface(k, j, i); - w_kind[f] = kind(cut.a_w[f]); - } - } - } + // Domain-side faces keep `Fluid` (the serial loops skipped them). + let u_kind: Vec = (0..g.n_ufaces()) + .into_par_iter() + .map(|f| { + let i = f % (nx + 1); + if i == 0 || i == nx { FaceKind::Fluid } else { kind(cut.a_u[f]) } + }) + .collect(); + let v_kind: Vec = (0..g.n_vfaces()) + .into_par_iter() + .map(|f| { + let j = (f / nx) % (ny + 1); + if j == 0 || j == ny { FaceKind::Fluid } else { kind(cut.a_v[f]) } + }) + .collect(); + let w_kind: Vec = (0..g.n_wfaces()) + .into_par_iter() + .map(|f| { + let k = f / nxy; + if !periodic && (k == 0 || k == nz) { FaceKind::Fluid } else { kind(cut.a_w[f]) } + }) + .collect(); let mut mask = Self { grid: g, periodic_z: periodic, @@ -254,9 +246,13 @@ impl Mask { .ok() .and_then(|v| v.parse().ok()) .unwrap_or(MERGE_FRACTION); - let small: Vec = (0..n) - .map(|idx| self.cell_active(idx) && frac(idx) < threshold) - .collect(); + let small: Vec = { + use rayon::prelude::*; + (0..n) + .into_par_iter() + .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]) { let (k, j, i) = g.kji(idx);