diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu index 3aa60ba..afa82cf 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu @@ -275,3 +275,159 @@ extern "C" __global__ void e3_geom_gather( long long s = (long long) w * idx[t]; for (int c = 0; c < w; ++c) out[(long long) w * t + c] = src[s + c]; } + +/* + * R6-2 step 2: the body's surface velocity on the device (the flag test's + * host closure: the centreline's velocity interpolated at the capsule's + * closest point where the span-cut capsule is not farther than the circle, + * zero on the circle, no z component), and its two uses on a moving cut + * mask: + * e3_geom_ub the imposition band's faces: the velocity at the foot + * (`Mask::surface_velocity_at`: the face centre plus its + * centroid shift when the centroid foot is on, φ and its + * unit gradient from the trilinear interpolant of the + * corner φ, the foot x − s n) into the dense table and a + * packed copy for the host mirror; + * e3_geom_impose the band's solid faces (not open at the instant, off the + * domain sides the host loops skip): the velocity at the + * face centre into the field (`Mask::impose_from`, whose + * ghost lists are empty on a cut mask); + * e3_geom_seam the periodic seam: w at k = nz takes k = 0's where both + * are solid. + */ +__device__ double body_velocity(double x, double y, double z, int c, const GeomSdf& s, + const double* __restrict__ poly, const double* __restrict__ vel) +{ + double best = 1.0 / 0.0, vx = 0.0, vy = 0.0; + for (int m = 0; m + 1 < s.npts; ++m) { + double ax = poly[2 * m], ay = poly[2 * m + 1]; + double bx = poly[2 * m + 2], by = poly[2 * m + 3]; + double ex = bx - ax, ey = by - ay; + double l2 = ex * ex + ey * ey; + double u = ((x - ax) * ex + (y - ay) * ey) / l2; + if (u < 0.0) u = 0.0; + if (u > 1.0) u = 1.0; + double px = ax + u * ex, py = ay + u * ey; + double qx = x - px, qy = y - py; + double d = sqrt(qx * qx + qy * qy); + if (d < best) { + best = d; + double avx = vel[2 * m], avy = vel[2 * m + 1]; + double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3]; + vx = avx + u * (bvx - avx); + vy = avy + u * (bvy - avy); + } + } + double df = best - s.half; + if (s.flag_cut) df = span_cut(df, z, s); + double ex0 = x - s.cx, ey0 = y - s.cy; + double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc; + if (s.cyl_cut) dc = span_cut(dc, z, s); + if (df <= dc) return c == 0 ? vx : (c == 1 ? vy : 0.0); + return 0.0; +} + +__device__ __forceinline__ double ub_node(const GeomGrid& g, const double* __restrict__ phi, + long long i, long long j, long long k) +{ + if (i < 0) i = 0; + if (i > g.nx) i = g.nx; + if (j < 0) j = 0; + if (j > g.ny) j = g.ny; + if (k < 0) k = 0; + if (k > g.nz) k = g.nz; + return phi[(k * (g.ny + 1) + j) * (g.nx + 1) + i]; +} + +__device__ __forceinline__ double ub_lerp(double a, double b, double f) { return a + f * (b - a); } + +/* The face centre of face `f` of component `c`. */ +__device__ __forceinline__ void face_centre(const GeomGrid& g, int c, long long f, double* x) +{ + long long ni = g.nx + (c == 0), nj = g.ny + (c == 1); + long long k = f / (nj * ni), j = (f / ni) % nj, i = f % ni; + x[0] = ((double) i + (c == 0 ? 0.0 : 0.5)) * g.dx; + x[1] = ((double) j + (c == 1 ? 0.0 : 0.5)) * g.dy; + x[2] = ((double) k + (c == 2 ? 0.0 : 0.5)) * g.dz; +} + +extern "C" __global__ void e3_geom_ub( + GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel, + const double* __restrict__ phi, int c, const unsigned int* __restrict__ faces, long long nf, + const double* __restrict__ shift, int use_shift, + double* __restrict__ ub, double* __restrict__ packed) +{ + long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (q >= nf) return; + long long f = faces[q]; + double xc[3]; + face_centre(g, c, f, xc); + double x0 = xc[0], x1 = xc[1], x2 = xc[2]; + if (use_shift) { + x0 = x0 + shift[3 * f]; + x1 = x1 + shift[3 * f + 1]; + x2 = x2 + shift[3 * f + 2]; + } + /* the interpolant's value and unit gradient */ + double gx = x0 / g.dx, gy = x1 / g.dy, gz = x2 / g.dz; + long long i0 = (long long) floor(gx), j0 = (long long) floor(gy), k0 = (long long) floor(gz); + double fx = gx - (double) i0, fy = gy - (double) j0, fz = gz - (double) k0; + double c000 = ub_node(g, phi, i0, j0, k0), c100 = ub_node(g, phi, i0 + 1, j0, k0); + double c010 = ub_node(g, phi, i0, j0 + 1, k0), c110 = ub_node(g, phi, i0 + 1, j0 + 1, k0); + double c001 = ub_node(g, phi, i0, j0, k0 + 1), c101 = ub_node(g, phi, i0 + 1, j0, k0 + 1); + double c011 = ub_node(g, phi, i0, j0 + 1, k0 + 1), c111 = ub_node(g, phi, i0 + 1, j0 + 1, k0 + 1); + double c00 = ub_lerp(c000, c100, fx); + double c10 = ub_lerp(c010, c110, fx); + double c01 = ub_lerp(c001, c101, fx); + double c11 = ub_lerp(c011, c111, fx); + double c0 = ub_lerp(c00, c10, fy); + double c1 = ub_lerp(c01, c11, fy); + double sd = ub_lerp(c0, c1, fz); + double dx0 = ub_lerp(c100 - c000, c110 - c010, fy); + double dx1 = ub_lerp(c101 - c001, c111 - c011, fy); + double px = ub_lerp(dx0, dx1, fz) / g.dx; + double dy0 = ub_lerp(c010 - c000, c110 - c100, fx); + double dy1 = ub_lerp(c011 - c001, c111 - c101, fx); + double py = ub_lerp(dy0, dy1, fz) / g.dy; + double dz0 = ub_lerp(c001 - c000, c101 - c100, fx); + double dz1 = ub_lerp(c011 - c010, c111 - c110, fx); + double pz = ub_lerp(dz0, dz1, fy) / g.dz; + double norm = sqrt(px * px + py * py + pz * pz); + double n0 = 1.0, n1 = 0.0, n2 = 0.0; + if (norm > 0.0) { + n0 = px / norm; + n1 = py / norm; + n2 = pz / norm; + } + double v = body_velocity(x0 - sd * n0, x1 - sd * n1, x2 - sd * n2, c, s, poly, vel); + ub[f] = v; + packed[q] = v; +} + +extern "C" __global__ void e3_geom_impose( + GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel, + int c, const unsigned int* __restrict__ faces, long long nf, + const int* __restrict__ open, double* __restrict__ field) +{ + long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (q >= nf) return; + long long f = faces[q]; + if (open[f] != 0) return; + long long ni = g.nx + (c == 0), nj = g.ny + (c == 1); + long long j = (f / ni) % nj, i = f % ni; + /* the host loops: u over i in 1..nx, v over j in 1..ny, every w */ + if (c == 0 && (i == 0 || i == g.nx)) return; + if (c == 1 && (j == 0 || j == g.ny)) return; + double xc[3]; + face_centre(g, c, f, xc); + field[f] = body_velocity(xc[0], xc[1], xc[2], c, s, poly, vel); +} + +extern "C" __global__ void e3_geom_seam(GeomGrid g, const int* __restrict__ open_w, double* __restrict__ w) +{ + long long n = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long nxy = (long long) g.nx * g.ny; + if (n >= nxy) return; + long long fn = (long long) g.nz * nxy + n; + if (open_w[n] == 0 && open_w[fn] == 0) w[fn] = w[n]; +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs index 3a8fae7..2fe83ec 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs @@ -33,6 +33,13 @@ pub struct DeviceSdf { pub fillet: f64, /// The centreline polyline's points (x, y) at `t`. pub poly: Vec<[f64; 2]>, + /// R6-2 step 2: the centreline's velocity (vx, vy) per point at `t`, the + /// device form of the surface velocity (`e3_geom.cu` `e3_geom_ub`): the + /// velocity interpolated at the capsule's closest point where the + /// capsule (span-cut) is not farther than the circle, zero on the circle, + /// no z component — the flag test's host closure. Empty: the surface + /// velocity stays on the host. + pub vel: Vec<[f64; 2]>, } /// One surface quadrature point: position, unit normal out of the solid, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs index 4af568d..4739409 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs @@ -48,6 +48,9 @@ pub struct MaskUpdate { pub correction: f64, /// The merged cells. pub merged: usize, + /// R6-2 step 2 (`RTX_E3_UB_DEVICE=1`): the device imposes the wall on + /// the uploaded field (the host rebuild skips `impose_from`). + pub impose_on_device: bool, } /// The instantaneous arrays of a retired mask, kept for the mask two steps later. 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 9d30146..4d1d0be 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 @@ -675,8 +675,11 @@ impl DeviceCut { (&mut self.a_pred, &mut self.d) } - /// `RTX_E3_BAND_CHECK=1`: every device table against a fresh full build. - pub(super) fn check_against(&self, full: &Self) { + /// `RTX_E3_BAND_CHECK=1`: every device table against a fresh full build + /// (the surface velocity on every face within `band` of the wall, the + /// imposition band — zeros included; beyond it the full build writes 0 + /// and ours keeps stale values, never read). + pub(super) fn check_against(&self, full: &Self, band: f64) { let rt = runtime(); let same_f = |name: &str, a: &CudaSlice, b: &CudaSlice| { let (x, y) = ( @@ -709,16 +712,22 @@ impl DeviceCut { 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). + // ub: the band's faces (|d| ≤ band, by the full build's distances). let (x, y) = ( rt.stream.memcpy_dtov(&self.ub[c]).expect("dl"), rt.stream.memcpy_dtov(&full.ub[c]).expect("dl"), ); + let d: Vec = rt.stream.memcpy_dtov(&full.d[c]).expect("dl"); let bad = x .iter() .zip(&y) - .filter(|(p, q)| **q != 0.0 && p.to_bits() != q.to_bits()) + .zip(&d) + .filter(|((p, q), dd)| dd.abs() <= band && p.to_bits() != q.to_bits()) .count(); + assert!( + x.len() == y.len() && y.len() == d.len(), + "band tables: ub[{c}] lengths differ" + ); assert!(bad == 0, "band tables: ub[{c}] differs on {bad} band faces"); } same_f("wall_flux", &self.wall_flux, &full.wall_flux); @@ -774,6 +783,51 @@ impl DeviceCut { } impl DeviceStep { + /// R6-2 step 2 under `RTX_E3_BAND_CHECK=1`: the device-imposed velocities + /// against the host's `impose_from` on the uploaded (un-imposed) mirror, + /// bit for bit on every face. + fn check_device_impose(&mut self, t: f64) { + let rt = runtime(); + let field = self.host_field.as_ref().expect("the host mirror"); + let mask = self.solver.mask().expect("mask"); + let body = self.solver.body().expect("body"); + let (mut u, mut v, mut w) = (field.u.clone(), field.v.clone(), field.w.clone()); + mask.impose_from( + body, + &field.u_old, + &field.v_old, + &field.w_old, + &mut u, + &mut v, + &mut w, + t, + ); + let mut bad = [0usize; 3]; + for (c, (host, dev)) in [(&u, &self.u), (&v, &self.v), (&w, &self.w)] + .into_iter() + .enumerate() + { + let d: Vec = rt.stream.memcpy_dtov(dev).expect("dl"); + bad[c] = host + .iter() + .zip(&d) + .filter(|(a, b)| a.to_bits() != b.to_bits()) + .count() + + host.len().abs_diff(d.len()); + } + if bad == [0; 3] { + eprintln!(" impose check t {t:.6}: device imposition IDENTICAL to the host (u, v, w every face)"); + } else { + eprintln!( + " impose check t {t:.6}: DIFFERS — {} / {} / {} faces", + bad[0], bad[1], bad[2] + ); + if std::env::var("RTX_E3_GEOM_CHECK_SOFT").is_err() { + panic!("impose check: the device imposition differs from the host"); + } + } + } + /// One step on the cut-cell mask (the host `advance` with the cut /// predictor, the apertured merged continuity and the owner-read /// corrections). @@ -858,6 +912,7 @@ impl DeviceStep { None => true, }; let mut device_mask: Option = None; + let mut ub_body: Option = None; if geom_on { if let Some(dc) = self.cut.as_mut() { let fresh_geom = self.geom.is_none(); @@ -900,6 +955,11 @@ impl DeviceStep { upd.merged ); } + // R6-2 step 2 (`RTX_E3_UB_DEVICE=1`): the band's surface + // velocity and the imposition on the device. + ub_body = super::geom::UbBody::new(&self.solver, t_new); + let mut upd = upd; + upd.impose_on_device = ub_body.is_some(); self.solver.pending_mask = Some(upd); } } @@ -923,17 +983,38 @@ impl DeviceStep { // 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(); + let ub_band = self + .solver + .mask() + .and_then(|m| m.impose_band()) + .unwrap_or(f64::INFINITY); self.cut = match self.cut.take() { // R6-2: the device classification wrote the tables; the // surface velocities remain. Some(mut prev) if device_mask.is_some() => { let dm = self.dmask.as_mut().expect("device mask"); let merged = device_mask.expect("merged"); - assert!(prev.update_after_device_mask(&self.solver, g, t_new, dm, merged)); + let device_side = match (self.geom.as_ref(), ub_body.as_ref()) { + (Some(gm), Some(ubb)) => { + Some((gm, ubb, [&mut self.u, &mut self.v, &mut self.w])) + } + _ => None, + }; + assert!(prev.update_after_device_mask( + &self.solver, + g, + t_new, + dm, + device_side, + merged + )); + if check && ub_body.is_some() { + self.check_device_impose(t_new); + } if check { let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new) .expect("full"); - prev.check_against(&full); + prev.check_against(&full, ub_band); eprintln!(" table check t {t_new:.6}: device tables IDENTICAL to the full build"); } Some(prev) @@ -943,7 +1024,7 @@ impl DeviceStep { if check { let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new) .expect("full"); - prev.check_against(&full); + prev.check_against(&full, ub_band); } Some(prev) } else { diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs index 5f85929..7747112 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs @@ -12,6 +12,7 @@ //! down. use super::cut::DeviceCut; +use crate::solvers::incompressible::embedded3::body::DeviceSdf; use crate::solvers::incompressible::embedded3::cut::{CutGeometry, GeomPool}; use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime}; use crate::solvers::incompressible::embedded3::step::Solver; @@ -32,6 +33,9 @@ struct GeomKernels { faces: CudaFunction, cells: CudaFunction, gather: CudaFunction, + ub: CudaFunction, + impose: CudaFunction, + seam: CudaFunction, } static GEOM_ONCE: OnceLock = OnceLock::new(); @@ -46,6 +50,9 @@ fn kernels() -> &'static GeomKernels { faces: f("e3_geom_faces"), cells: f("e3_geom_cells"), gather: f("e3_geom_gather"), + ub: f("e3_geom_ub"), + impose: f("e3_geom_impose"), + seam: f("e3_geom_seam"), _module: module, } }) @@ -116,6 +123,80 @@ pub(super) fn log_fallback(reason: &str) { }); } +fn geom_grid(g: Grid) -> GeomGrid { + GeomGrid { + nx: g.nx as i32, + ny: g.ny as i32, + nz: g.nz as i32, + dx: g.dx, + dy: g.dy, + dz: g.dz, + } +} + +fn geom_sdf(sdf: &DeviceSdf) -> GeomSdf { + GeomSdf { + cx: sdf.cyl[0], + cy: sdf.cyl[1], + rc: sdf.cyl[2], + zc: sdf.zc, + span: sdf.span, + r_edge: sdf.r_edge, + half: sdf.half, + fillet: sdf.fillet, + cyl_cut: i32::from(sdf.cyl_cut), + flag_cut: i32::from(sdf.flag_cut), + npts: sdf.poly.len() as i32, + } +} + +/// Interleaved (x, y) pairs on the device (one dummy entry when empty). +fn upload_pairs(v: &[[f64; 2]]) -> CudaSlice { + let flat: Vec = v.iter().flat_map(|p| [p[0], p[1]]).collect(); + runtime() + .stream + .memcpy_stod(if flat.is_empty() { + &[0.0f64][..] + } else { + &flat + }) + .expect("upload pairs") +} + +/// R6-2 step 2 (`RTX_E3_UB_DEVICE=1`, default off): the band's surface +/// velocity at the foot and the end-of-rebuild imposition on the device. +pub(crate) fn ub_device_enabled() -> bool { + std::env::var("RTX_E3_UB_DEVICE").is_ok_and(|v| v == "1") +} + +/// R6-2 step 2: the body's surface velocity in the device's form at one +/// time (the polyline and its velocities uploaded once per step). +pub(super) struct UbBody { + gs: GeomSdf, + poly: CudaSlice, + vel: CudaSlice, +} + +impl UbBody { + /// `None` unless `RTX_E3_UB_DEVICE=1` (default off until its + /// two-period gate) and the body has a device form of its surface + /// velocity (`DeviceSdf::vel`, one per polyline point). + pub(super) fn new(solver: &Solver, t: f64) -> Option { + if !ub_device_enabled() { + return None; + } + let sdf = solver.body()?.device_sdf(t)?; + if sdf.vel.is_empty() || sdf.vel.len() != sdf.poly.len() { + return None; + } + Some(Self { + gs: geom_sdf(&sdf), + poly: upload_pairs(&sdf.poly), + vel: upload_pairs(&sdf.vel), + }) + } +} + /// The persistent device geometry: corner φ and bound, the evaluated /// corners, the cell volumes and wall vectors (the face tables live in /// `DeviceCut`). @@ -147,6 +228,108 @@ fn par_copy(src: &[T]) -> Vec { } impl DeviceGeom { + /// R6-2 step 2: the imposition on the `nf` band faces of component `c` + /// listed in `faces` — the solid ones (`open` = the instantaneous kind + /// is fluid) take the body's surface velocity at the face centre — and, + /// for `c = 2` on a periodic mask, the seam. + #[allow(clippy::too_many_arguments)] + pub(super) fn impose( + body: &UbBody, + g: Grid, + c: usize, + faces: &CudaSlice, + nf: usize, + open: &CudaSlice, + field: &mut CudaSlice, + ) { + if nf == 0 { + return; + } + let rt = runtime(); + let k = kernels(); + let gg = geom_grid(g); + let ci = c as i32; + let nf64 = nf as i64; + unsafe { + rt.stream + .launch_builder(&k.impose) + .arg(&gg) + .arg(&body.gs) + .arg(&body.poly) + .arg(&body.vel) + .arg(&ci) + .arg(faces) + .arg(&nf64) + .arg(open) + .arg(field) + .launch(cfg(nf)) + .expect("e3_geom_impose"); + } + } + + /// The periodic seam after the imposition (`Mask::impose_from`'s last loop). + pub(super) fn seam(g: Grid, open_w: &CudaSlice, w: &mut CudaSlice) { + let rt = runtime(); + let k = kernels(); + let gg = geom_grid(g); + unsafe { + rt.stream + .launch_builder(&k.seam) + .arg(&gg) + .arg(open_w) + .arg(w) + .launch(cfg(g.nx * g.ny)) + .expect("e3_geom_seam"); + } + } + + /// R6-2 step 2: the surface velocity at the foot for the `nf` faces of + /// component `c` listed in `faces` (`Mask::surface_velocity_at` on the + /// device, from this state's corner φ — the mirror's): written into the + /// dense table `ub` and packed into `packed` (the host mirror's copy). + #[allow(clippy::too_many_arguments)] + pub(super) fn surface_velocity( + &self, + body: &UbBody, + g: Grid, + c: usize, + faces: &CudaSlice, + nf: usize, + shift: Option<&CudaSlice>, + ub: &mut CudaSlice, + packed: &mut CudaSlice, + ) { + if nf == 0 { + return; + } + let rt = runtime(); + let k = kernels(); + let gg = geom_grid(g); + let ci = c as i32; + let nf64 = nf as i64; + let use_shift = i32::from(shift.is_some()); + // Without the shift the kernel never reads it: any f64 buffer stands in. + let shift = shift.unwrap_or(&self.vol); + unsafe { + rt.stream + .launch_builder(&k.ub) + .arg(&gg) + .arg(&body.gs) + .arg(&body.poly) + .arg(&body.vel) + .arg(&self.phi) + .arg(&ci) + .arg(faces) + .arg(&nf64) + .arg(shift) + .arg(&use_shift) + .arg(ub) + .arg(packed) + .launch(cfg(nf)) + .expect("e3_geom_ub"); + } + } + pub(super) fn new(g: Grid) -> Self { let rt = runtime(); let nn = (g.nx + 1) * (g.ny + 1) * (g.nz + 1); @@ -215,36 +398,9 @@ impl DeviceGeom { full = true; } } - let gg = GeomGrid { - nx: nx as i32, - ny: ny as i32, - nz: nz as i32, - dx: g.dx, - dy: g.dy, - dz: g.dz, - }; - let gs = GeomSdf { - cx: sdf.cyl[0], - cy: sdf.cyl[1], - rc: sdf.cyl[2], - zc: sdf.zc, - span: sdf.span, - r_edge: sdf.r_edge, - half: sdf.half, - fillet: sdf.fillet, - cyl_cut: i32::from(sdf.cyl_cut), - flag_cut: i32::from(sdf.flag_cut), - npts: sdf.poly.len() as i32, - }; - let poly: Vec = sdf.poly.iter().flat_map(|p| [p[0], p[1]]).collect(); - let d_poly = rt - .stream - .memcpy_stod(if poly.is_empty() { - &[0.0f64][..] - } else { - &poly - }) - .expect("poly"); + let gg = geom_grid(g); + let gs = geom_sdf(&sdf); + let d_poly = upload_pairs(&sdf.poly); let (has_prev, band, motion) = match prev { Some((_, band, motion)) => (1i32, band, motion), None => (0i32, 0.0, 0.0), diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/mask.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/mask.rs index 5a19349..a7cf7ad 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/mask.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/mask.rs @@ -10,7 +10,7 @@ //! kernels overwrite them. use super::cut::{DeviceCut, Phase}; -use super::geom::DeviceGeom; +use super::geom::{DeviceGeom, UbBody}; use crate::solvers::incompressible::embedded3::cutwall::MERGE_FRACTION; use crate::solvers::incompressible::embedded3::maskupdate::MaskUpdate; use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime}; @@ -157,6 +157,8 @@ pub(super) struct DeviceMask { out_face_a: [Option>; 3], out_face_shift: [Option>; 3], out_wall_a: Option>, + /// Step 2: the band's surface velocities, packed (the host mirror's copy). + out_ub: Option>, /// The laps of the last run (ms): lists, cells + merging, faces, GCL, fold. pub(super) laps: [f64; 6], } @@ -190,6 +192,7 @@ impl DeviceMask { out_face_a: [None, None, None], out_face_shift: [None, None, None], out_wall_a: None, + out_ub: None, laps: [0.0; 6], } } @@ -597,31 +600,31 @@ impl DeviceMask { face_shift, correction, merged, + impose_on_device: false, } } - /// The imposition band's faces per component (|d| <= band), ascending. - pub(super) fn band_faces(&mut self, dc: &DeviceCut, g: Grid, band: f64) -> [Vec; 3] { + /// The imposition band's faces of component `c` (|d| ≤ `band`), + /// ascending, compacted into `list_tmp`; their count. + fn band_list(&mut self, dc: &DeviceCut, g: Grid, band: f64, c: usize) -> usize { let rt = runtime(); let k = kernels(); - let counts = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()]; - [0usize, 1, 2].map(|c| { - let n64 = counts[c] as i64; - unsafe { - rt.stream - .launch_builder(&k.band_flags) - .arg(&n64) - .arg(&dc.d[c]) - .arg(&band) - .arg(&mut self.flags) - .launch(cfg(counts[c])) - .expect("e3_mask_band_flags"); - } - let mut tmp = self.list_tmp.take(); - let n = self.compact(counts[c], &mut tmp); - self.list_tmp = tmp; - head(self.list_tmp.as_ref().expect("band"), n) - }) + let n = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()][c]; + let n64 = n as i64; + unsafe { + rt.stream + .launch_builder(&k.band_flags) + .arg(&n64) + .arg(&dc.d[c]) + .arg(&band) + .arg(&mut self.flags) + .launch(cfg(n)) + .expect("e3_mask_band_flags"); + } + let mut tmp = self.list_tmp.take(); + let nb = self.compact(n, &mut tmp); + self.list_tmp = tmp; + nb } } @@ -629,13 +632,19 @@ impl DeviceCut { /// R6-2: finish the projection tables after the device classification /// (`DeviceMask::run` wrote every table but the surface velocities): the /// surface velocity at the foot for the imposition band's faces (listed - /// on the device), the merged count, the phase. + /// on the device) — on the device from `geom`'s corner φ when the body + /// has a device form of its surface velocity (step 2; the host mirror + /// by a packed download), on the host otherwise — the merged count, the + /// phase. `RTX_E3_BAND_CHECK=1`: the device values against the host's + /// `surface_velocity_at`, bit for bit on every band face. + #[allow(clippy::too_many_arguments)] pub(super) fn update_after_device_mask( &mut self, solver: &Solver, g: Grid, t: f64, dm: &mut DeviceMask, + device_side: Option<(&DeviceGeom, &UbBody, [&mut CudaSlice; 3])>, merged: usize, ) -> bool { use crate::solvers::incompressible::embedded3::poisson::device_cg::scatter_f64; @@ -645,12 +654,24 @@ impl DeviceCut { }; let lap = Instant::now(); let band = mask.impose_band().unwrap_or(f64::INFINITY); - let band_faces = dm.band_faces(self, g, band); - let l_band = lap.elapsed(); let h = [g.dx, g.dy, g.dz]; let shifts = mask.face_shift_tables(); + let centroid_foot = shifts.is_some() && mask.wall_foot_centroid; let (nx, ny) = (g.nx, g.ny); - for (c, bf) in band_faces.iter().enumerate() { + let check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); + let on_device = device_side.is_some(); + let (geom, ubody, mut fields) = match device_side { + Some((gm, ubb, f)) => (Some(gm), Some(ubb), Some(f)), + None => (None, None, None), + }; + let mut l_band = std::time::Duration::ZERO; + let mut nbf = [0usize; 3]; + let mut checked = 0usize; + for c in 0..3 { + let lb = Instant::now(); + let nb = dm.band_list(self, g, band, c); + l_band += lb.elapsed(); + nbf[c] = nb; 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); @@ -660,25 +681,73 @@ impl DeviceCut { (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2], ] }; - 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; + let host_ub = |bf: &[u32]| -> 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() + }; + let list = dm.list_tmp.as_ref().expect("band list"); + let bf = head(list, nb); + match (geom, ubody) { + (Some(gm), Some(ubb)) => { + ensure(&mut dm.out_ub, nb); + let packed = dm.out_ub.as_mut().expect("packed ub"); + let shift = centroid_foot.then_some(&self.shift[c]); + gm.surface_velocity(ubb, g, c, list, nb, shift, &mut self.ub[c], packed); + let ubv = head(dm.out_ub.as_ref().expect("packed ub"), nb); + for (&f, &v) in bf.iter().zip(&ubv) { + self.ub_host[c][f as usize] = v; + } + if check { + let reference = host_ub(&bf); + let bad = reference + .iter() + .zip(&ubv) + .filter(|(a, b)| a.to_bits() != b.to_bits()) + .count(); + if bad > 0 { + eprintln!( + " ub check t {t:.6}: DIFFERS — component {c}: {bad} of {nb} band faces" + ); + if std::env::var("RTX_E3_GEOM_CHECK_SOFT").is_err() { + panic!( + "ub check: the device surface velocity differs from the host" + ); + } + } + checked += nb; + } + // The imposition on the uploaded field: the same band list. + let field = &mut *fields.as_mut().expect("fields")[c]; + DeviceGeom::impose(ubb, g, c, list, nb, &self.open_pred[c], field); + if c == 2 { + DeviceGeom::seam(g, &self.open_pred[2], field); + } + } + _ => { + let ubv = host_ub(&bf); + for (&f, &v) in bf.iter().zip(&ubv) { + self.ub_host[c][f as usize] = v; + } + scatter_f64(&bf, &ubv, &mut self.ub[c]); + } } - scatter_f64(bf, &ubv, &mut self.ub[c]); + } + if check && on_device { + eprintln!( + " ub check t {t:.6}: device surface velocity IDENTICAL to the host on {checked} band faces" + ); } self.merged = merged; self.phase = Phase::Projection; @@ -686,12 +755,13 @@ impl DeviceCut { runtime().stream.synchronize().expect("sync"); if std::env::var("RTX_E3_MOVING_PROFILE").is_ok() { eprintln!( - " table laps (device mask): band list {:.0} ms, ub {:.0} ms ({} / {} / {} band faces)", + " table laps (device mask): band list {:.0} ms, ub {:.0} ms ({} / {} / {} band faces, {})", l_band.as_secs_f64() * 1e3, (lap.elapsed() - l_band).as_secs_f64() * 1e3, - band_faces[0].len(), - band_faces[1].len(), - band_faces[2].len() + nbf[0], + nbf[1], + nbf[2], + if on_device { "device, imposed" } else { "host" } ); } true 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 1c96e5e..8876d4a 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 @@ -152,17 +152,25 @@ impl Solver { let fresh_cells = refill_fresh_in(&old_mask, &new_mask, &mut field.p, &upd.changed); let l_step = lap.elapsed(); let sub = std::time::Instant::now(); - let body = self.body.as_ref().expect("body"); - new_mask.impose_from( - body, - &field.u_old, - &field.v_old, - &field.w_old, - &mut field.u, - &mut field.v, - &mut field.w, - t_new, + // R6-2 step 2: the device imposes after the upload (a cut mask has + // no ghosts: the imposition is the band's solid faces and the seam). + assert!( + !upd.impose_on_device || new_mask.ghost_faces() == 0, + "R6-2 step 2: the device imposition has no ghost reconstruction" ); + if !upd.impose_on_device { + let body = self.body.as_ref().expect("body"); + new_mask.impose_from( + body, + &field.u_old, + &field.v_old, + &field.w_old, + &mut field.u, + &mut field.v, + &mut field.w, + t_new, + ); + } let s_impose = sub.elapsed(); self.last_ghost_correction = upd.correction; if let Some((r, p_ref, fresh_ref, table, correction)) = reference { diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index f45a2bd..58dfc15 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -415,6 +415,14 @@ fn flag_wake_on_the_device() { .collect(), None => analytic_polyline(t).iter().map(|p| [p.0, p.1]).collect(), }, + // R6-2 step 2: the centreline's velocity per point (the analytic mode is transverse). + vel: match recorded() { + Some(rec) => recorded_polyline(rec, t) + .iter() + .map(|p| [p.2, p.3]) + .collect(), + None => analytic_polyline(t).iter().map(|p| [0.0, p.2]).collect(), + }, }); solver.set_moving_body(body); let g = Grid::cubic(nx, ny_grid, nz, h);