diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu new file mode 100644 index 0000000..3aa60ba --- /dev/null +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu @@ -0,0 +1,277 @@ +/** + * embedded3 R6-1: the body's φ and the cut geometry on the device — the + * host `cut.rs::CutGeometry::build_from` expression for expression in fp64 + * (compiled with FMA contraction off; `/` and `sqrt` are IEEE round-to- + * nearest in double on the device, as on the host), and the flag test's + * φ (`tests/embedded3_flag_wake.rs`: the circle, the capsule around the + * centreline polyline, the span cut with rounded edges, the fillet union). + * + * Kernels: + * e3_geom_phi per corner: the narrow band's keep-or-evaluate (bound − + * motion > band keeps φ), φ from the polyline otherwise; + * e3_geom_faces per face of one component: the aperture by the two + * Kuhn triangles and the face-centre φ, where a corner was + * evaluated (every face on a full pass); + * e3_geom_cells per cell: the fluid volume by the six Kuhn tetrahedra + * and the wall vector by closure of the face apertures; + * e3_geom_gather compact copies for the host mirror. + */ + +struct GeomSdf { + double cx, cy, rc; /* the circle */ + double zc, span, r_edge;/* the span cut */ + double half; /* capsule half-thickness */ + double fillet; /* root fillet radius (0: min) */ + int cyl_cut, flag_cut; /* cut to the span */ + int npts; /* polyline points (x, y interleaved) */ +}; + +struct GeomGrid { + int nx, ny, nz; + double dx, dy, dz; +}; + +__device__ __forceinline__ double rs_min(double a, double b) { return b < a ? b : a; } +__device__ __forceinline__ double rs_max(double a, double b) { return b > a ? b : a; } + +/* The span cut with rounded edges (flag_3d / cylinder_3d). */ +__device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s) +{ + double r = s.r_edge; + double q1 = d2 + r; + double q2 = fabs(z - s.zc) - 0.5 * s.span + r; + double m1 = rs_max(q1, 0.0), m2 = rs_max(q2, 0.0); + double outside = sqrt(m1 * m1 + m2 * m2); + return outside + rs_min(rs_max(q1, q2), 0.0) - r; +} + +__device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, const double* __restrict__ poly) +{ + /* the circle */ + 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); + /* the capsule: distance to the polyline */ + double best = 1.0 / 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 df = best - s.half; + if (s.flag_cut) df = span_cut(df, z, s); + double r = s.fillet; + if (r > 0.0 && dc < r && df < r) { + double a = r - dc, b = r - df; + return r - sqrt(a * a + b * b); + } + return rs_min(dc, df); +} + +/* Per corner. `has_prev`: keep φ where the decayed bound stays beyond the band. */ +extern "C" __global__ void e3_geom_phi( + GeomGrid g, GeomSdf s, const double* __restrict__ poly, + int has_prev, double band, double motion, + double* __restrict__ phi, double* __restrict__ bound, unsigned char* __restrict__ touched) +{ + long long n = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long nn = (long long) (g.nx + 1) * (g.ny + 1) * (g.nz + 1); + if (n >= nn) return; + long long nxy = (long long) (g.nx + 1) * (g.ny + 1); + int k = (int) (n / nxy); + int j = (int) ((n / (g.nx + 1)) % (g.ny + 1)); + int i = (int) (n % (g.nx + 1)); + if (has_prev) { + double b = bound[n] - motion; + if (b > band) { + bound[n] = b; + touched[n] = 0; + return; + } + } + double v = geom_phi_at((double) i * g.dx, (double) j * g.dy, (double) k * g.dz, s, poly); + phi[n] = v; + bound[n] = fabs(v); + touched[n] = 1; +} + +__device__ __forceinline__ long long gnode(const GeomGrid& g, int k, int j, int i) +{ + return ((long long) k * (g.ny + 1) + j) * (g.nx + 1) + i; +} + +/* cut.rs tri_area_fraction */ +__device__ double tri_frac(double p0, double p1, double p2) +{ + double v[3] = {p0, p1, p2}; + int pos = (p0 >= 0.0) + (p1 >= 0.0) + (p2 >= 0.0); + if (pos == 0) return 0.0; + if (pos == 3) return 1.0; + if (pos == 1) { + int a = v[0] >= 0.0 ? 0 : (v[1] >= 0.0 ? 1 : 2); + int b = (a + 1) % 3, c = (a + 2) % 3; + return (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c])); + } + int a = v[0] < 0.0 ? 0 : (v[1] < 0.0 ? 1 : 2); + int b = (a + 1) % 3, c = (a + 2) % 3; + return 1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c])); +} + +__device__ __forceinline__ double quad_frac(double q00, double q10, double q01, double q11) +{ + return 0.5 * (tri_frac(q00, q10, q11) + tri_frac(q00, q11, q01)); +} + +/* Per face of component c (0 u, 1 v, 2 w): aperture and face-centre φ. */ +extern "C" __global__ void e3_geom_faces( + GeomGrid g, int c, int full, const double* __restrict__ phi, const unsigned char* __restrict__ touched, + double* __restrict__ a, double* __restrict__ d) +{ + long long f = (long long) blockIdx.x * blockDim.x + threadIdx.x; + int ni = g.nx + (c == 0), nj = g.ny + (c == 1), nk = g.nz + (c == 2); + if (f >= (long long) ni * nj * nk) return; + int k = (int) (f / ((long long) nj * ni)); + int j = (int) ((f / ni) % nj); + int i = (int) (f % ni); + long long n00, n10, n01, n11; + if (c == 0) { /* x-face: (j, k) */ + n00 = gnode(g, k, j, i); n10 = gnode(g, k, j + 1, i); + n01 = gnode(g, k + 1, j, i); n11 = gnode(g, k + 1, j + 1, i); + } else if (c == 1) { /* y-face: (i, k) */ + n00 = gnode(g, k, j, i); n10 = gnode(g, k, j, i + 1); + n01 = gnode(g, k + 1, j, i); n11 = gnode(g, k + 1, j, i + 1); + } else { /* z-face: (i, j) */ + n00 = gnode(g, k, j, i); n10 = gnode(g, k, j, i + 1); + n01 = gnode(g, k, j + 1, i); n11 = gnode(g, k, j + 1, i + 1); + } + if (!full && !(touched[n00] | touched[n10] | touched[n01] | touched[n11])) return; + double q00 = phi[n00], q10 = phi[n10], q01 = phi[n01], q11 = phi[n11]; + a[f] = quad_frac(q00, q10, q01, q11); + d[f] = 0.25 * (q00 + q10 + q01 + q11); +} + +__device__ __forceinline__ double det3(const double* a, const double* b, const double* c) +{ + return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + + a[2] * (b[0] * c[1] - b[1] * c[0]); +} + +__device__ double tet_vol(const double* p0, const double* p1, const double* p2, const double* p3) +{ + double e1[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; + double e2[3] = {p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; + double e3[3] = {p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2]}; + return fabs(det3(e1, e2, e3)) / 6.0; +} + +__device__ __forceinline__ void lerp3(const double* a, const double* b, double t, double* out) +{ + out[0] = a[0] + t * (b[0] - a[0]); + out[1] = a[1] + t * (b[1] - a[1]); + out[2] = a[2] + t * (b[2] - a[2]); +} + +/* cut.rs tet_fluid_volume */ +__device__ double tet_fluid(const double pts[4][3], const double* v) +{ + double total = tet_vol(pts[0], pts[1], pts[2], pts[3]); + int pos[4], neg[4], np = 0, nn = 0; + for (int q = 0; q < 4; ++q) { + if (v[q] >= 0.0) pos[np++] = q; + } + for (int q = 0; q < 4; ++q) { + if (v[q] < 0.0) neg[nn++] = q; + } + if (np == 0) return 0.0; + if (np == 4) return total; + if (np == 1 || np == 3) { + int a = np == 1 ? pos[0] : neg[0]; + const int* o = np == 1 ? neg : pos; + double pb[3], pc[3], pd[3]; + lerp3(pts[a], pts[o[0]], v[a] / (v[a] - v[o[0]]), pb); + lerp3(pts[a], pts[o[1]], v[a] / (v[a] - v[o[1]]), pc); + lerp3(pts[a], pts[o[2]], v[a] / (v[a] - v[o[2]]), pd); + double t = tet_vol(pts[a], pb, pc, pd); + return np == 1 ? t : total - t; + } + int a = pos[0], b = pos[1], c = neg[0], d = neg[1]; + double pac[3], pad[3], pbc[3], pbd[3]; + lerp3(pts[a], pts[c], v[a] / (v[a] - v[c]), pac); + lerp3(pts[a], pts[d], v[a] / (v[a] - v[d]), pad); + lerp3(pts[b], pts[c], v[b] / (v[b] - v[c]), pbc); + lerp3(pts[b], pts[d], v[b] / (v[b] - v[d]), pbd); + return tet_vol(pts[a], pts[b], pbc, pbd) + tet_vol(pts[a], pac, pbc, pbd) + + tet_vol(pts[a], pac, pad, pbd); +} + +/* The Kuhn split around (0,0,0)–(1,1,1), corners as (x, y, z). */ +__constant__ int KUHN[6][4][3] = { + {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {1, 1, 1}}, + {{0, 0, 0}, {1, 0, 0}, {1, 0, 1}, {1, 1, 1}}, + {{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 1, 1}}, + {{0, 0, 0}, {0, 1, 0}, {0, 1, 1}, {1, 1, 1}}, + {{0, 0, 0}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}}, + {{0, 0, 0}, {0, 0, 1}, {0, 1, 1}, {1, 1, 1}}, +}; + +/* Per cell: fluid volume fraction and the wall vector (3 interleaved). */ +extern "C" __global__ void e3_geom_cells( + GeomGrid g, int full, const double* __restrict__ phi, const unsigned char* __restrict__ touched, + const double* __restrict__ a_u, const double* __restrict__ a_v, const double* __restrict__ a_w, + double* __restrict__ vol, double* __restrict__ wall, unsigned char* __restrict__ touched_cell) +{ + long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long nxy = (long long) g.nx * g.ny; + if (idx >= nxy * g.nz) return; + int k = (int) (idx / nxy); + int j = (int) ((idx / g.nx) % g.ny); + int i = (int) (idx % g.nx); + unsigned char t = 0; + for (int dk = 0; dk < 2; ++dk) + for (int dj = 0; dj < 2; ++dj) + for (int di = 0; di < 2; ++di) + t |= touched[gnode(g, k + dk, j + dj, i + di)]; + touched_cell[idx] = t; + if (!full && !t) return; + double fluid = 0.0; + for (int q = 0; q < 6; ++q) { + double pts[4][3]; + double vals[4]; + for (int m = 0; m < 4; ++m) { + pts[m][0] = (double) KUHN[q][m][0]; + pts[m][1] = (double) KUHN[q][m][1]; + pts[m][2] = (double) KUHN[q][m][2]; + vals[m] = phi[gnode(g, k + KUHN[q][m][2], j + KUHN[q][m][1], i + KUHN[q][m][0])]; + } + fluid += tet_fluid(pts, vals); + } + vol[idx] = fluid; + double ax = g.dy * g.dz, ay = g.dx * g.dz, az = g.dx * g.dy; + long long fu0 = ((long long) k * g.ny + j) * (g.nx + 1) + i; + long long fv0 = ((long long) k * (g.ny + 1) + j) * g.nx + i; + long long fw0 = ((long long) k * g.ny + j) * g.nx + i; + double sx = (a_u[fu0 + 1] - a_u[fu0]) * ax; + double sy = (a_v[fv0 + g.nx] - a_v[fv0]) * ay; + double sz = (a_w[fw0 + nxy] - a_w[fw0]) * az; + wall[3 * idx] = -sx; + wall[3 * idx + 1] = -sy; + wall[3 * idx + 2] = -sz; +} + +/* out[w t + c] = src[w idx[t] + c], c < w (compact copies for the host mirror). */ +extern "C" __global__ void e3_geom_gather( + int n, int w, const unsigned int* __restrict__ idx, const double* __restrict__ src, double* __restrict__ out) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + long long s = (long long) w * idx[t]; + for (int c = 0; c < w; ++c) out[(long long) w * t + c] = src[s + c]; +} 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 2d48bb7..2662ac5 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs @@ -6,6 +6,34 @@ type SdfFn = Box f64 + Send + Sync>; type VelFn = Box (f64, f64, f64) + Send + Sync>; type SamplerFn = Box Vec + Send + Sync>; +type DeviceSdfFn = Box DeviceSdf + Send + Sync>; + +/// R6-1: a body's φ in a form a device kernel evaluates (`e3_geom.cu` +/// `e3_geom_phi`), the same arithmetic as the flag test's host closure: +/// a circle in (x, y) (optionally cut to the flag's span with rounded +/// edges) united with a capsule of half-thickness `half` around the +/// centreline polyline `poly` at `t` (optionally cut to the span), by a +/// concave fillet of radius `fillet` (the plain `min` at 0). The host +/// passes the polyline per step (a small upload). +#[derive(Debug, Clone)] +pub struct DeviceSdf { + /// The circle: centre (x, y) and radius. + pub cyl: [f64; 3], + /// The circle cut to the span (`RTX_E3_FLAG_CYL_SPAN=flag`). + pub cyl_cut: bool, + /// The capsule cut to the span (a finite-span flag). + pub flag_cut: bool, + /// The span's centre in z, the span and the edge radius. + pub zc: f64, + pub span: f64, + pub r_edge: f64, + /// The capsule's half-thickness. + pub half: f64, + /// The root fillet radius (0 = the plain union). + pub fillet: f64, + /// The centreline polyline's points (x, y) at `t`. + pub poly: Vec<[f64; 2]>, +} /// One surface quadrature point: position, unit normal out of the solid, /// the area it represents. @@ -24,6 +52,7 @@ pub struct Body { phi: SdfFn, velocity: Option, sampler: Option, + device_sdf: Option, } impl Body { @@ -35,6 +64,7 @@ impl Body { phi: Box::new(phi), velocity: None, sampler: None, + device_sdf: None, } } @@ -123,6 +153,24 @@ impl Body { self } + /// R6-1: attach the device form of φ (`RTX_E3_GEOM_DEVICE=1` evaluates + /// it on the band instead of the host closure; the caller guarantees + /// the two are the same arithmetic). + #[must_use] + pub fn with_device_sdf(mut self, f: F) -> Self + where + F: Fn(f64) -> DeviceSdf + Send + Sync + 'static, + { + self.device_sdf = Some(Box::new(f)); + self + } + + /// The device form of φ at `t`, when the body has one. + #[must_use] + pub fn device_sdf(&self, t: f64) -> Option { + self.device_sdf.as_ref().map(|f| f(t)) + } + #[inline] #[must_use] pub fn phi(&self, x: f64, y: f64, z: f64, t: f64) -> f64 { 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 4760802..4b33e1e 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs @@ -39,6 +39,72 @@ pub struct CutGeometry { pub touched: Vec, /// Per cell: any of its eight corners touched (computed once per build). pub touched_cell: Vec, + /// R6-1: the device geometry's generation of this build (0 = built on + /// the host); a device generation's arrays are recycled (`GeomPool`). + pub generation: u64, +} + +/// R6-1: the arrays of retired device-built geometries, kept for the next +/// device build's host mirror to overwrite in place of fresh allocations +/// (the generation says which band lists bring an array up to date). +#[derive(Debug, Default)] +pub struct GeomPool { + f64s: Vec<(u64, u8, Vec)>, + walls: Vec<(u64, Vec<[f64; 3]>)>, + flags: Vec>, +} + +impl GeomPool { + /// Array kinds: apertures u/v/w, face distances u/v/w, φ, bound. + pub const A: [u8; 3] = [0, 1, 2]; + pub const D: [u8; 3] = [3, 4, 5]; + pub const PHI: u8 = 7; + pub const BOUND: u8 = 8; + + /// Keep `v` (generation `generation` of kind `kind`); host-built + /// generations (0) are dropped. + pub fn put(&mut self, generation: u64, kind: u8, v: Vec) { + if generation == 0 { + return; + } + self.f64s.retain(|e| e.1 != kind); + self.f64s.push((generation, kind, v)); + } + + #[must_use] + pub fn take(&mut self, kind: u8) -> Option<(u64, Vec)> { + let at = self.f64s.iter().position(|e| e.1 == kind)?; + let (g, _, v) = self.f64s.swap_remove(at); + Some((g, v)) + } + + #[must_use] + pub fn take_wall(&mut self) -> Option<(u64, Vec<[f64; 3]>)> { + self.walls.pop() + } + + #[must_use] + pub fn take_flags(&mut self) -> Option> { + self.flags.pop() + } + + /// A retired geometry's arrays that are not kept elsewhere. + pub fn put_cut(&mut self, c: CutGeometry) { + let g = c.generation; + if g == 0 { + return; + } + self.put(g, Self::D[0], c.d_u); + self.put(g, Self::D[1], c.d_v); + self.put(g, Self::D[2], c.d_w); + self.put(g, Self::PHI, c.phi); + self.put(g, Self::BOUND, c.bound); + self.walls.clear(); + self.walls.push((g, c.wall)); + self.flags.truncate(2); + self.flags.push(c.touched); + self.flags.push(c.touched_cell); + } } impl CutGeometry { @@ -222,6 +288,7 @@ impl CutGeometry { bound, touched, touched_cell, + generation: 0, } } 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 ad7b385..1eb2124 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -124,7 +124,20 @@ impl Mask { b: Boundaries, prev: Option<(&CutGeometry, f64, f64)>, ) -> Result { + let lap = std::time::Instant::now(); let cut = CutGeometry::build_from(body, g, t, prev); + if std::env::var("RTX_E3_MOVING_PROFILE").is_ok() { + eprintln!( + " mask laps: cut geometry (host, within build_mask) {:.0} ms", + lap.elapsed().as_secs_f64() * 1e3 + ); + } + Self::from_cut(cut, g, b) + } + + /// R6-1: the mask of a cut geometry built elsewhere (the device's + /// `DeviceGeom` mirror); `build_cut_from` is this after `CutGeometry::build_from`. + pub fn from_cut(cut: CutGeometry, g: Grid, b: Boundaries) -> Result { 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); diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs index 8171147..e04393d 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -21,7 +21,7 @@ pub mod reconstruct; pub mod step; pub mod wall; -pub use body::{Body, SurfaceSample}; +pub use body::{Body, DeviceSdf, SurfaceSample}; pub use cut::CutGeometry; pub use export_vtk::write_vtk; pub use field::Field; 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 622990f..acc1359 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 @@ -6,6 +6,7 @@ //! off so the predictors are the host's arithmetic to the bit. mod cut; +mod geom; use super::{Side, Solver, StepResult}; use crate::solvers::incompressible::embedded3::Grid; @@ -162,6 +163,8 @@ pub struct DeviceStep { cut: Option, /// Steps since the multigrid hierarchy was last rebuilt (moving bodies). steps_since_hierarchy: usize, + /// R6-1: the persistent device cut geometry (`RTX_E3_GEOM_DEVICE=1`). + geom: Option, } impl DeviceStep { @@ -240,6 +243,7 @@ impl DeviceStep { initialized: false, cut, steps_since_hierarchy: 0, + geom: None, } } 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 37816c9..2c1687c 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 @@ -91,6 +91,9 @@ pub(super) struct DeviceCut { /// six interleaved (S2-7b `axisfoot`; one dummy entry when off). foot: [CudaSlice; 3], pub(super) merged: usize, + /// R6-1: this step's predictor apertures and distances were written by + /// the device geometry (`update` skips their scatters). + pub(super) geom_written: bool, } /// Which phase the tables serve: the predictor reads the instantaneous @@ -369,6 +372,7 @@ impl DeviceCut { vn: [up_f(&vn[0]), up_f(&vn[1]), up_f(&vn[2])], foot: [up_f(&foot[0]), up_f(&foot[1]), up_f(&foot[2])], merged, + geom_written: false, }; if profile { let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; @@ -469,8 +473,10 @@ impl DeviceCut { 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]); + if !self.geom_written { + 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 @@ -562,6 +568,7 @@ impl DeviceCut { self.fold_idx = up_u(&fold_idx); self.merged = merged; self.phase = Phase::Projection; + self.geom_written = false; rt.stream.synchronize().expect("sync"); if profile { let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; @@ -581,6 +588,14 @@ impl DeviceCut { true } + /// R6-1: the predictor apertures and the face distances, for the device + /// geometry to write in place. + pub(super) fn geometry_targets( + &mut self, + ) -> (&mut [CudaSlice; 3], &mut [CudaSlice; 3]) { + (&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) { let rt = runtime(); @@ -731,6 +746,30 @@ impl DeviceStep { let mut field = self.host_field.take().unwrap_or_else(|| Field::new(g)); self.download_for_rebuild(&mut field); let l_down = lap.elapsed(); + // R6-1 (`RTX_E3_GEOM_DEVICE=1`): the end-of-step cut geometry on + // the device, its face tables written into the persistent + // predictor set; the host classifies the mirror. + if super::geom::enabled() { + if let Some(dc) = self.cut.as_mut() { + let geom = self + .geom + .get_or_insert_with(|| super::geom::DeviceGeom::new(g)); + let mut pool = std::mem::take(&mut self.solver.geom_pool); + let built = geom.build(&self.solver, &mut pool, g, t_new, dt, dc); + self.solver.geom_pool = pool; + if let Some(mirror) = built { + self.solver.pending_cut = Some(mirror); + dc.geom_written = true; + } + } + if profile { + eprintln!( + " moving laps: device geometry {:.0} ms", + (lap.elapsed() - l_down).as_secs_f64() * 1e3 + ); + } + } + let l_geom = lap.elapsed(); fresh_cells = self.solver.rebuild_moving_mask(&mut field, dt, t_new); let l_mask = lap.elapsed(); self.upload_after_rebuild(&field); @@ -760,7 +799,7 @@ impl DeviceStep { eprintln!( " moving laps: download {:.0} ms, host mask {:.0} ms, upload {:.0} ms, projection tables {:.0} ms", l_down.as_secs_f64() * 1e3, - (l_mask - l_down).as_secs_f64() * 1e3, + (l_mask - l_geom).as_secs_f64() * 1e3, (l_up - l_mask).as_secs_f64() * 1e3, (l_tables - l_up).as_secs_f64() * 1e3 ); 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 new file mode 100644 index 0000000..3eb1b5b --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs @@ -0,0 +1,624 @@ +//! R6-1: the body's φ and the cut geometry on the device (`e3_geom.cu`), +//! behind `RTX_E3_GEOM_DEVICE=1`. The corner φ (the narrow band's +//! keep-or-evaluate, φ from the body's `DeviceSdf` polyline), the face +//! apertures and face-centre φ (written straight into the persistent +//! `DeviceCut` predictor apertures and distances), the cell volumes and +//! wall vectors — `CutGeometry::build_from` in fp64 with the host's +//! operation order. The host mirror (the `CutGeometry` the rest of the +//! moving rebuild classifies) is the previous geometry with the band's +//! entries replaced by compact device copies: only what changed comes +//! down. + +use super::cut::DeviceCut; +use crate::solvers::incompressible::embedded3::Grid; +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; +use crate::solvers::incompressible::embedded3::wall::WallScheme; +use cudarc::driver::{ + CudaFunction, CudaModule, CudaSlice, DeviceRepr, PushKernelArg, ValidAsZeroBits, +}; +use rayon::prelude::*; +use std::sync::{Arc, OnceLock}; +use std::time::Instant; + +const GEOM_KERNELS: &str = include_str!("../../../../../kernels/cuda/e3_geom.cu"); + +struct GeomKernels { + _module: Arc, + phi: CudaFunction, + faces: CudaFunction, + cells: CudaFunction, + gather: CudaFunction, +} + +static GEOM_ONCE: OnceLock = OnceLock::new(); + +fn kernels() -> &'static GeomKernels { + GEOM_ONCE.get_or_init(|| { + // FMA contraction off: the host's roundings. + let module = load_module(GEOM_KERNELS, "e3_geom.cu", true); + let f = |name: &str| module.load_function(name).expect(name); + GeomKernels { + phi: f("e3_geom_phi"), + faces: f("e3_geom_faces"), + cells: f("e3_geom_cells"), + gather: f("e3_geom_gather"), + _module: module, + } + }) +} + +/// `struct GeomSdf` in e3_geom.cu. +#[repr(C)] +#[derive(Clone, Copy)] +struct GeomSdf { + cx: f64, + cy: f64, + rc: f64, + zc: f64, + span: f64, + r_edge: f64, + half: f64, + fillet: f64, + cyl_cut: i32, + flag_cut: i32, + npts: i32, +} +unsafe impl DeviceRepr for GeomSdf {} +unsafe impl ValidAsZeroBits for GeomSdf {} + +/// `struct GeomGrid` in e3_geom.cu. +#[repr(C)] +#[derive(Clone, Copy)] +struct GeomGrid { + nx: i32, + ny: i32, + nz: i32, + dx: f64, + dy: f64, + dz: f64, +} +unsafe impl DeviceRepr for GeomGrid {} +unsafe impl ValidAsZeroBits for GeomGrid {} + +/// `RTX_E3_GEOM_DEVICE=1`. +pub(super) fn enabled() -> bool { + std::env::var("RTX_E3_GEOM_DEVICE").is_ok_and(|v| v == "1") +} + +/// The persistent device geometry: corner φ and bound, the evaluated +/// corners, the cell volumes and wall vectors (the face tables live in +/// `DeviceCut`). +pub(super) struct DeviceGeom { + phi: CudaSlice, + bound: CudaSlice, + touched: CudaSlice, + vol: CudaSlice, + wall: CudaSlice, + touched_cell: CudaSlice, + /// The mirror this state produced last (its φ buffer's address): the + /// host's previous geometry is ours only while it is still that one. + synced: usize, + /// The generation of the last build, the first of the current chain + /// (the last re-sync), and the band lists (faces per component, cells) + /// of the recent generations: what changed from one to the next. + generation: u64, + chain_start: u64, + history: std::collections::VecDeque<(u64, [Vec; 3], Vec)>, +} + +/// Generations of band lists kept (a recycled array older than this is copied instead). +const HISTORY: usize = 4; + +fn par_copy(src: &[T]) -> Vec { + src.par_iter().copied().collect() +} + +impl DeviceGeom { + pub(super) fn new(g: Grid) -> Self { + let rt = runtime(); + let nn = (g.nx + 1) * (g.ny + 1) * (g.nz + 1); + let nc = g.cells(); + Self { + phi: rt.stream.alloc_zeros::(nn).expect("alloc"), + bound: rt.stream.alloc_zeros::(nn).expect("alloc"), + touched: rt.stream.alloc_zeros::(nn).expect("alloc"), + vol: rt.stream.alloc_zeros::(nc).expect("alloc"), + wall: rt.stream.alloc_zeros::(3 * nc).expect("alloc"), + touched_cell: rt.stream.alloc_zeros::(nc).expect("alloc"), + synced: 0, + generation: 0, + chain_start: 1, + history: std::collections::VecDeque::new(), + } + } + + /// The cut geometry of the solver's body at `t` (the moving rebuild's + /// end-of-step geometry, the band of the current mask moved by the + /// step `dt`) on the device; the face tables written into `dc`'s + /// predictor apertures and distances. Returns the host mirror, or + /// `None` when the body has no device form (the host builds it). + pub(super) fn build( + &mut self, + solver: &Solver, + pool: &mut GeomPool, + g: Grid, + t: f64, + dt: f64, + dc: &mut DeviceCut, + ) -> Option { + if solver.params.wall_scheme != WallScheme::CutCell || solver.params.aperture_substeps != 0 + { + return None; + } + let body = solver.body()?; + let sdf = body.device_sdf(t)?; + let rt = runtime(); + let k = kernels(); + let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); + let lap = Instant::now(); + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let nn = (nx + 1) * (ny + 1) * (nz + 1); + let nc = g.cells(); + let counts = [(nx + 1) * ny * nz, nx * (ny + 1) * nz, nx * ny * (nz + 1)]; + let prev = solver.band_prev(g, dt); + // The device state is the previous geometry's only if the host's + // previous geometry is the mirror we produced; otherwise re-sync. + let mut full = prev.is_none(); + if let Some((p, _, _)) = prev { + if p.phi.as_ptr() as usize != self.synced || p.phi.len() != nn { + rt.stream + .memcpy_htod(&p.phi, &mut self.phi) + .expect("sync phi"); + rt.stream + .memcpy_htod(&p.bound, &mut self.bound) + .expect("sync bound"); + 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 (has_prev, band, motion) = match prev { + Some((_, band, motion)) => (1i32, band, motion), + None => (0i32, 0.0, 0.0), + }; + let full_i = i32::from(full); + unsafe { + rt.stream + .launch_builder(&k.phi) + .arg(&gg) + .arg(&gs) + .arg(&d_poly) + .arg(&has_prev) + .arg(&band) + .arg(&motion) + .arg(&mut self.phi) + .arg(&mut self.bound) + .arg(&mut self.touched) + .launch(cfg(nn)) + .expect("e3_geom_phi"); + } + let (a_dev, d_dev) = dc.geometry_targets(); + for c in 0..3 { + let ci = c as i32; + unsafe { + rt.stream + .launch_builder(&k.faces) + .arg(&gg) + .arg(&ci) + .arg(&full_i) + .arg(&self.phi) + .arg(&self.touched) + .arg(&mut a_dev[c]) + .arg(&mut d_dev[c]) + .launch(cfg(counts[c])) + .expect("e3_geom_faces"); + } + } + unsafe { + rt.stream + .launch_builder(&k.cells) + .arg(&gg) + .arg(&full_i) + .arg(&self.phi) + .arg(&self.touched) + .arg(&a_dev[0]) + .arg(&a_dev[1]) + .arg(&a_dev[2]) + .arg(&mut self.vol) + .arg(&mut self.wall) + .arg(&mut self.touched_cell) + .launch(cfg(nc)) + .expect("e3_geom_cells"); + } + let l_launch = lap.elapsed(); + let generation_new = self.generation + 1; + if full { + self.history.clear(); + self.chain_start = generation_new; + } + // The host mirror, while the kernels run: the corner pass (the same + // keep-or-evaluate arithmetic) into recycled buffers. + let recycle = |pool: &mut GeomPool, kind: u8, n: usize| -> Vec { + match pool.take(kind) { + Some((_, v)) if v.len() == n => v, + _ => vec![0.0; n], + } + }; + let mut phi = recycle(pool, GeomPool::PHI, nn); + let mut bound = recycle(pool, GeomPool::BOUND, nn); + let mut touched = match pool.take_flags() { + Some(v) if v.len() == nn => v, + _ => vec![true; nn], + }; + phi.par_iter_mut() + .zip(bound.par_iter_mut()) + .zip(touched.par_iter_mut()) + .enumerate() + .for_each(|(n, ((phi_n, bound_n), touched_n))| { + if let Some((p, band, motion)) = prev { + let b = p.bound[n] - motion; + if b > band { + *phi_n = p.phi[n]; + *bound_n = b; + *touched_n = false; + return; + } + } + *phi_n = 0.0; + *bound_n = 0.0; + *touched_n = true; + }); + let node_list: Vec = (0..nn) + .into_par_iter() + .filter(|&n| touched[n]) + .map(|n| n as u32) + .collect(); + // The evaluated cells from the device (its corner flags are the + // host's: the same test on the same bounds). + let tc_dev: Vec = rt + .stream + .memcpy_dtov(&self.touched_cell) + .expect("touched_cell"); + let mut touched_cell = match pool.take_flags() { + Some(v) if v.len() == nc => v, + _ => vec![false; nc], + }; + touched_cell + .par_iter_mut() + .zip(tc_dev.par_iter()) + .for_each(|(t, &d)| *t = d != 0); + let cell_list: Vec = (0..nc) + .into_par_iter() + .filter(|&i| touched_cell[i]) + .map(|i| i as u32) + .collect(); + // The faces of the evaluated cells (every face with an evaluated + // corner is one), per component, ascending. + let face_lists: [Vec; 3] = [0usize, 1, 2].map(|c| { + let mut v: Vec = cell_list + .par_iter() + .flat_map_iter(|&idx| { + let idx = idx as usize; + let (k, j, i) = (idx / (ny * nx), (idx / nx) % ny, idx % nx); + let pair = match c { + 0 => [g.uface(k, j, i), g.uface(k, j, i + 1)], + 1 => [g.vface(k, j, i), g.vface(k, j + 1, i)], + _ => [g.wface(k, j, i), g.wface(k + 1, j, i)], + }; + pair.into_iter().map(|f| f as u32) + }) + .collect(); + v.par_sort_unstable(); + v.dedup(); + v + }); + self.history + .push_back((generation_new, face_lists.clone(), cell_list.clone())); + while self.history.len() > HISTORY { + self.history.pop_front(); + } + let l_lists = lap.elapsed(); + // Each table: a recycled array of generation g0 brought up to date by + // the union of the band lists since g0, or a copy of the previous one. + let history = &self.history; + let chain_start = self.chain_start; + type Band = (u64, [Vec; 3], Vec); + let since = |g0: u64, pick: &dyn Fn(&Band) -> &[u32]| -> Option> { + if g0 < chain_start + || g0 >= generation_new + || !(g0 + 1..=generation_new).all(|q| history.iter().any(|h| h.0 == q)) + { + return None; + } + let mut v: Vec = history + .iter() + .filter(|h| h.0 > g0) + .flat_map(|h| pick(h).iter().copied()) + .collect(); + v.par_sort_unstable(); + v.dedup(); + Some(v) + }; + let prev_face: [[&[f64]; 3]; 2] = match prev { + Some((p, _, _)) => [[&p.a_u, &p.a_v, &p.a_w], [&p.d_u, &p.d_v, &p.d_w]], + None => [[&[], &[], &[]], [&[], &[], &[]]], + }; + let mut recycled = 0usize; + let mut face_tab: Vec<(Vec, Vec)> = Vec::with_capacity(6); + for (t, kinds) in [GeomPool::A, GeomPool::D].iter().enumerate() { + for c in 0..3 { + let got = pool.take(kinds[c]).filter(|(_, v)| v.len() == counts[c]); + let entry = match (prev, got) { + (Some(_), Some((g0, v))) => match since(g0, &|h| &h.1[c]) { + Some(idx) => { + recycled += 1; + (v, idx) + } + None => (par_copy(prev_face[t][c]), face_lists[c].clone()), + }, + (Some(_), None) => (par_copy(prev_face[t][c]), face_lists[c].clone()), + (None, _) => (vec![0.0; counts[c]], face_lists[c].clone()), + }; + face_tab.push(entry); + } + } + let (mut wall, wall_idx) = match (prev, pool.take_wall().filter(|(_, v)| v.len() == nc)) { + (Some(_), Some((g0, v))) => match since(g0, &|h| &h.2) { + Some(idx) => { + recycled += 1; + (v, idx) + } + None => (par_copy(&prev.unwrap().0.wall), cell_list.clone()), + }, + (Some((p, _, _)), None) => (par_copy(&p.wall), cell_list.clone()), + (None, _) => (vec![[0.0; 3]; nc], cell_list.clone()), + }; + let mut vol = match prev { + Some((p, _, _)) => par_copy(&p.vol), + None => vec![0.0; nc], + }; + let l_copy = lap.elapsed(); + // The band's entries: gathered on the device into one buffer. + let total = node_list.len() + + face_tab.iter().map(|e| e.1.len()).sum::() + + cell_list.len() + + 3 * wall_idx.len(); + let mut packed = rt.stream.alloc_zeros::(total.max(1)).expect("alloc"); + let mut off = 0usize; + let mut gather = + |list: &[u32], w: usize, src: &CudaSlice, packed: &mut CudaSlice| { + if list.is_empty() { + return; + } + let d_idx = rt.stream.memcpy_stod(list).expect("idx"); + let n = list.len() as i32; + let wi = w as i32; + let mut dst = packed.slice_mut(off..off + w * list.len()); + unsafe { + rt.stream + .launch_builder(&k.gather) + .arg(&n) + .arg(&wi) + .arg(&d_idx) + .arg(src) + .arg(&mut dst) + .launch(cfg(list.len())) + .expect("e3_geom_gather"); + } + off += w * list.len(); + }; + gather(&node_list, 1, &self.phi, &mut packed); + for (e, entry) in face_tab.iter().enumerate() { + let (t, c) = (e / 3, e % 3); + let src = if t == 0 { &a_dev[c] } else { &d_dev[c] }; + gather(&entry.1, 1, src, &mut packed); + } + gather(&cell_list, 1, &self.vol, &mut packed); + gather(&wall_idx, 3, &self.wall, &mut packed); + let host: Vec = rt.stream.memcpy_dtov(&packed).expect("download"); + let l_down = lap.elapsed(); + let mut off = 0usize; + let mut take = |n: usize| { + let s = &host[off..off + n]; + off += n; + s + }; + let vals = take(node_list.len()); + for (&n, &v) in node_list.iter().zip(vals) { + phi[n as usize] = v; + bound[n as usize] = v.abs(); + } + for entry in face_tab.iter_mut() { + let vals = take(entry.1.len()); + for (&f, &v) in entry.1.iter().zip(vals) { + entry.0[f as usize] = v; + } + } + let vv = take(cell_list.len()); + for (&i, &v) in cell_list.iter().zip(vv) { + vol[i as usize] = v; + } + let vw = take(3 * wall_idx.len()); + for (m, &i) in wall_idx.iter().enumerate() { + wall[i as usize] = [vw[3 * m], vw[3 * m + 1], vw[3 * m + 2]]; + } + let mut tabs = face_tab.into_iter().map(|e| e.0); + let mut next = || tabs.next().expect("table"); + let (a_u, a_v, a_w, d_u, d_v, d_w) = (next(), next(), next(), next(), next(), next()); + let mirror = CutGeometry { + grid: g, + phi, + a_u, + a_v, + a_w, + vol, + wall, + d_u, + d_v, + d_w, + bound, + touched, + touched_cell, + generation: generation_new, + }; + self.generation = generation_new; + self.synced = mirror.phi.as_ptr() as usize; + if profile { + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + eprintln!( + " geom laps (device): launch {:.0} ms, host corner pass + lists {:.0} ms, copy previous {:.0} ms, gather + download {:.0} ms, scatter {:.0} ms ({} corners / {} {} {} faces / {} cells; {} of 7 tables recycled{})", + ms(l_launch), + ms(l_lists - l_launch), + ms(l_copy - l_lists), + ms(l_down - l_copy), + ms(lap.elapsed() - l_down), + node_list.len(), + face_lists[0].len(), + face_lists[1].len(), + face_lists[2].len(), + cell_list.len(), + recycled, + if full { ", full pass" } else { "" } + ); + } + if std::env::var("RTX_E3_BAND_CHECK").is_ok() { + self.check(solver, g, t, dt, &mirror, dc); + } + Some(mirror) + } + + /// `RTX_E3_BAND_CHECK=1` with the device geometry: the mirror and the + /// device tables against the host's `CutGeometry::build_from`, bit for bit. + fn check( + &self, + solver: &Solver, + g: Grid, + t: f64, + dt: f64, + mirror: &CutGeometry, + dc: &mut DeviceCut, + ) { + let rt = runtime(); + let body = solver.body().expect("body"); + let reference = CutGeometry::build_from(body, g, t, solver.band_prev(g, dt)); + let mut report: Vec = Vec::new(); + let cmp = |report: &mut Vec, name: &str, x: &[f64], y: &[f64]| { + assert_eq!(x.len(), y.len(), "geom check: {name} length"); + let mut bad = 0usize; + let mut max_abs = 0.0f64; + for (p, q) in x.iter().zip(y) { + if p.to_bits() != q.to_bits() { + bad += 1; + max_abs = max_abs.max((p - q).abs()); + } + } + if bad > 0 { + report.push(format!("{name}: {bad} differ (max |Δ| {max_abs:.3e})")); + } + }; + let flat = + |w: &[[f64; 3]]| -> Vec { w.iter().flat_map(|v| v.iter().copied()).collect() }; + // The mirror. + cmp(&mut report, "mirror phi", &mirror.phi, &reference.phi); + cmp(&mut report, "mirror bound", &mirror.bound, &reference.bound); + cmp(&mut report, "mirror a_u", &mirror.a_u, &reference.a_u); + cmp(&mut report, "mirror a_v", &mirror.a_v, &reference.a_v); + cmp(&mut report, "mirror a_w", &mirror.a_w, &reference.a_w); + cmp(&mut report, "mirror d_u", &mirror.d_u, &reference.d_u); + cmp(&mut report, "mirror d_v", &mirror.d_v, &reference.d_v); + cmp(&mut report, "mirror d_w", &mirror.d_w, &reference.d_w); + cmp(&mut report, "mirror vol", &mirror.vol, &reference.vol); + cmp( + &mut report, + "mirror wall", + &flat(&mirror.wall), + &flat(&reference.wall), + ); + let flags_ok = + mirror.touched == reference.touched && mirror.touched_cell == reference.touched_cell; + if !flags_ok { + report.push("mirror touched flags differ".into()); + } + // The device tables, whole. + let dl = |s: &CudaSlice| -> Vec { rt.stream.memcpy_dtov(s).expect("dl") }; + cmp(&mut report, "device phi", &dl(&self.phi), &reference.phi); + cmp( + &mut report, + "device bound", + &dl(&self.bound), + &reference.bound, + ); + cmp(&mut report, "device vol", &dl(&self.vol), &reference.vol); + cmp( + &mut report, + "device wall", + &dl(&self.wall), + &flat(&reference.wall), + ); + let (a_dev, d_dev) = dc.geometry_targets(); + let refs_a: [&[f64]; 3] = [&reference.a_u, &reference.a_v, &reference.a_w]; + let refs_d: [&[f64]; 3] = [&reference.d_u, &reference.d_v, &reference.d_w]; + for c in 0..3 { + cmp( + &mut report, + &format!("device a[{c}]"), + &dl(&a_dev[c]), + refs_a[c], + ); + cmp( + &mut report, + &format!("device d[{c}]"), + &dl(&d_dev[c]), + refs_d[c], + ); + } + let tc: Vec = rt.stream.memcpy_dtov(&self.touched_cell).expect("dl"); + if tc + .iter() + .zip(&reference.touched_cell) + .any(|(&p, &q)| (p != 0) != q) + { + report.push("device touched_cell differs".into()); + } + if report.is_empty() { + eprintln!( + " geom check t {t:.6}: device geometry IDENTICAL to the host build (mirror + device tables)" + ); + } else { + eprintln!(" geom check t {t:.6}: DIFFERS — {}", report.join("; ")); + if std::env::var("RTX_E3_GEOM_CHECK_SOFT").is_err() { + panic!("geom check: device geometry differs from the host build"); + } + } + } +} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index cb755e8..ca230af 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -296,6 +296,14 @@ pub struct Solver { initialized: bool, /// `(setup ns, iterate ns, solves, CG iterations)` summed. poisson_profile: (u64, u64, u64, u64), + /// R6-1: the next moving rebuild's cut geometry, built on the device + /// (`RTX_E3_GEOM_DEVICE=1`); `build_mask` classifies it instead of + /// evaluating φ on the host. + pub(super) pending_cut: Option, + /// R6-1: retired device-built geometry arrays (empty on the host path) + /// and the generation of `apertures_old`. + pub(super) geom_pool: super::cut::GeomPool, + pub(super) apertures_old_gen: u64, } impl Solver { @@ -329,6 +337,9 @@ impl Solver { time: 0.0, initialized: false, poisson_profile: (0, 0, 0, 0), + pending_cut: None, + geom_pool: super::cut::GeomPool::default(), + apertures_old_gen: 0, } } @@ -362,17 +373,33 @@ impl Solver { self.mask = None; } - fn build_mask(&self, body: &Body, g: Grid, t: f64, dt: f64) -> Mask { + /// The narrow band a moving cut rebuild re-evaluates φ in: the current + /// geometry, the band half-width `3h` and the body's largest motion + /// over `dt` (`None` without a speed bound or a cut mask). + pub(super) fn band_prev(&self, g: Grid, dt: f64) -> Option<(&CutGeometry, f64, f64)> { + let h = g.dx.min(g.dy).min(g.dz); + match (self.params.max_surface_speed, &self.mask) { + (Some(speed), Some(m)) => m.cut().map(|c| (c, 3.0 * h, speed * dt)), + _ => None, + } + } + + fn build_mask( + &self, + body: &Body, + g: Grid, + t: f64, + dt: f64, + pending: Option, + ) -> Mask { match self.params.wall_scheme { WallScheme::GhostBinary => Mask::build(body, g, t, self.params.boundaries), - WallScheme::CutCell => { - let h = g.dx.min(g.dy).min(g.dz); - let prev = match (self.params.max_surface_speed, &self.mask) { - (Some(speed), Some(m)) => m.cut().map(|c| (c, 3.0 * h, speed * dt)), - _ => None, - }; - Mask::build_cut_from(body, g, t, self.params.boundaries, prev) - } + WallScheme::CutCell => match pending { + Some(cut) => Mask::from_cut(cut, g, self.params.boundaries), + None => { + Mask::build_cut_from(body, g, t, self.params.boundaries, self.band_prev(g, dt)) + } + }, } .map(|mut m| { let lap_closures = std::time::Instant::now(); @@ -729,7 +756,7 @@ impl Solver { let t = self.time; if let Some(body) = &self.body { if self.mask.is_none() { - self.mask = Some(self.build_mask(body, field.grid, t, 0.0)); + self.mask = Some(self.build_mask(body, field.grid, t, 0.0, None)); } } self.apply_boundary_normals(field, t); 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 7b6e4ac..7f5565e 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 @@ -14,12 +14,13 @@ impl Solver { /// step-averaged apertures and the GCL wall-flux table (cut wall). /// Returns the fresh-cell count. `field` holds the predicted field. pub fn rebuild_moving_mask(&mut self, field: &mut Field, dt: f64, t_new: f64) -> usize { + let pending = self.pending_cut.take(); let Some(body) = &self.body else { return 0; }; let mut fresh_cells = 0; let lap = std::time::Instant::now(); - let mut new_mask = self.build_mask(body, field.grid, t_new, dt); + let mut new_mask = self.build_mask(body, field.grid, t_new, dt, pending); let l_build = lap.elapsed(); if let Some(old_mask) = self.mask.as_mut() { fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); @@ -78,11 +79,19 @@ impl Solver { if let Some(mut old) = self.mask.take() { let fluid = std::mem::take(&mut old.cell_fluid); if let Some(c) = old.cut.as_mut() { - self.apertures_old = Some([ + let retired = self.apertures_old.replace([ std::mem::take(&mut c.a_u), std::mem::take(&mut c.a_v), std::mem::take(&mut c.a_w), ]); + // R6-1: a device-built generation's arrays go back to the + // pool (the host path's generation is 0: dropped as before). + if let Some(r) = retired { + for (kind, v) in super::super::cut::GeomPool::A.into_iter().zip(r) { + self.geom_pool.put(self.apertures_old_gen, kind, v); + } + } + self.apertures_old_gen = c.generation; let mut vol = std::mem::take(&mut c.vol); use rayon::prelude::*; vol.par_iter_mut() @@ -93,6 +102,9 @@ impl Solver { } }); self.vol_old = vol; + if let Some(c) = old.cut.take() { + self.geom_pool.put_cut(c); + } } else { self.apertures_old = None; self.vol_old = vec![0.0; field.grid.cells()]; diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index a8dba59..36d9dba 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -26,7 +26,8 @@ use embedded3_flag_kinematics::{Recorded, recorded}; use rtx_cfd::solvers::incompressible::ConvectionScheme; use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep; use rtx_cfd::solvers::incompressible::embedded3::{ - Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, write_vtk, + Body, Boundaries, DeviceSdf, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, + write_vtk, }; use std::io::Write as _; @@ -141,8 +142,7 @@ fn flag_2d_recorded(rec: &Recorded, x: f64, y: f64, t: f64) -> (f64, (f64, f64)) POLY.with(|cell| { let mut c = cell.borrow_mut(); if c.0.to_bits() != t.to_bits() { - c.1 = rec.at(t, body_cy()); - inset_last(&mut c.1, tip_inset()); + c.1 = recorded_polyline(rec, t); c.0 = t; } let pts = &c.1; @@ -165,8 +165,36 @@ fn flag_2d_recorded(rec: &Recorded, x: f64, y: f64, t: f64) -> (f64, (f64, f64)) }) } +/// The recorded centreline at `t` with the tip inset. +fn recorded_polyline(rec: &Recorded, t: f64) -> Vec<(f64, f64, f64, f64)> { + let mut pts = rec.at(t, body_cy()); + inset_last(&mut pts, tip_inset()); + pts +} + +/// The analytic centreline's segments. +const N: usize = 40; + +/// The analytic centreline at `t` (x, y, transverse velocity), the tip inset. +fn analytic_polyline(t: f64) -> [(f64, f64, f64); N + 1] { + let mut p = [(0.0, 0.0, 0.0); N + 1]; + for (m, q) in p.iter_mut().enumerate() { + let s = m as f64 / N as f64; + let (d, v) = deflection(s, t); + *q = (FLAG_X0 + s * FLAG_LEN, body_cy() + d, v); + } + let inset = tip_inset(); + if inset > 0.0 { + let (ax, ay, _) = p[N - 1]; + let (bx, by, bv) = p[N]; + let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); + let f = (1.0 - inset / len).max(0.0); + p[N] = (ax + f * (bx - ax), ay + f * (by - ay), bv); + } + p +} + fn flag_2d_analytic(x: f64, y: f64, t: f64) -> (f64, f64) { - const N: usize = 40; // The centreline polyline at `t`, once per thread and time (PERF-3 // P1-2): the solver asks for the surface velocity at ~10⁶ faces per // step and each call rebuilt the 41 points (four hyperbolic / trigonometric @@ -178,19 +206,7 @@ fn flag_2d_analytic(x: f64, y: f64, t: f64) -> (f64, f64) { let pts = POLYLINE.with(|cell| { let mut c = cell.borrow_mut(); if c.0.to_bits() != t.to_bits() { - for (m, p) in c.1.iter_mut().enumerate() { - let s = m as f64 / N as f64; - let (d, v) = deflection(s, t); - *p = (FLAG_X0 + s * FLAG_LEN, body_cy() + d, v); - } - let inset = tip_inset(); - if inset > 0.0 { - let (ax, ay, _) = c.1[N - 1]; - let (bx, by, bv) = c.1[N]; - let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); - let f = (1.0 - inset / len).max(0.0); - c.1[N] = (ax + f * (bx - ax), ay + f * (by - ay), bv); - } + c.1 = analytic_polyline(t); c.0 = t; } c.1 @@ -380,6 +396,26 @@ fn flag_wake_on_the_device() { (0.0, 0.0, 0.0) } }); + // R6-1: the same φ in the device's form (`RTX_E3_GEOM_DEVICE=1`): the + // circle, the capsule around the step's centreline, the span cuts. + let body = body.with_device_sdf(move |t| DeviceSdf { + cyl: [CX, cy, R_CYL], + cyl_cut: !(flag_span() >= duct_depth() + || !std::env::var("RTX_E3_FLAG_CYL_SPAN").is_ok_and(|v| v == "flag")), + flag_cut: flag_span() < duct_depth(), + zc: 0.5 * duct_depth(), + span: flag_span(), + r_edge, + half: FLAG_HALF, + fillet: r_fillet, + poly: match recorded() { + Some(rec) => recorded_polyline(rec, t) + .iter() + .map(|p| [p.0, p.1]) + .collect(), + None => analytic_polyline(t).iter().map(|p| [p.0, p.1]).collect(), + }, + }); solver.set_moving_body(body); let g = Grid::cubic(nx, ny_grid, nz, h); let mut field = Field::new(g);