diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/e3_mask.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_mask.cu new file mode 100644 index 0000000..fab1e3a --- /dev/null +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/e3_mask.cu @@ -0,0 +1,428 @@ +/** + * embedded3 R6-2: the moving mask's classification on the device, behind + * `RTX_E3_MASK_DEVICE=1` (with `RTX_E3_GEOM_DEVICE=1`). From the device + * geometry of the step (R6-1: the end-of-step apertures in DeviceCut's + * predictor set, the volumes and wall vectors in DeviceGeom) and the + * previous step's (snapshots taken before the geometry kernels overwrite + * them), the host `Mask` expressions in fp64 with FMA contraction off: + * + * e3_mask_changed the changed set (touched by either build, dilated by + * one; `set_step_apertures_with`); + * e3_mask_face_flags the faces of the changed cells per component; + * e3_mask_cells per changed cell: fluid (vol > 0), the space-time + * activity, the GCL entry (`gcl_flux_table`); + * e3_mask_merge per changed cell: the virtual-merging master + * (`compute_merging` with the old mask); + * e3_mask_faces per face of a changed cell: the kind, the step + * aperture and open flag (`set_step_apertures`), the + * open-part centroid shift (`compute_face_shifts`); + * e3_mask_wall_flags / e3_mask_wall_area / e3_mask_wall_flux + * the GCL table's wall cells, areas and the table + * with the compatibility correction; + * e3_mask_slave_flags / e3_mask_fold_count / e3_mask_fold_fill + * the merged cells' CSR (slaves ascending per master); + * e3_mask_band_flags the imposition band's faces (|d| <= band); + * e3_scan_* block scans: stream compaction (ascending) and the + * exclusive scan of counts. + * Every per-entry value depends on its own inputs alone; the two serial + * sums of the GCL table (in ascending order) are the host's. + */ + +typedef unsigned int u32; +typedef unsigned char u8; + +struct MaskGrid { + int nx, ny, nz, periodic; + double dx, dy, dz; +}; + +#define SCAN_BLOCK 1024 + +/* Exclusive block scan of v over a 1024-thread block; *total = the block sum. */ +__device__ u32 block_excl_scan(u32 v, u32* total) +{ + __shared__ u32 warp_sums[32]; + int lane = threadIdx.x & 31, wid = threadIdx.x >> 5; + u32 x = v; + for (int o = 1; o < 32; o <<= 1) { + u32 y = __shfl_up_sync(0xffffffffu, x, o); + if (lane >= o) x += y; + } + if (lane == 31) warp_sums[wid] = x; + __syncthreads(); + if (wid == 0) { + u32 s = warp_sums[lane]; + for (int o = 1; o < 32; o <<= 1) { + u32 y = __shfl_up_sync(0xffffffffu, s, o); + if (lane >= o) s += y; + } + warp_sums[lane] = s; + } + __syncthreads(); + u32 incl = x + (wid > 0 ? warp_sums[wid - 1] : 0u); + *total = warp_sums[31]; + __syncthreads(); + return incl - v; +} + +/* Per block of 1024: the number of flags set. */ +extern "C" __global__ void e3_scan_count_flags(long long n, const u8* __restrict__ flags, u32* __restrict__ sums) +{ + long long i = (long long) blockIdx.x * SCAN_BLOCK + threadIdx.x; + u32 total; + block_excl_scan(i < n ? (u32) (flags[i] != 0) : 0u, &total); + if (threadIdx.x == 0) sums[blockIdx.x] = total; +} + +/* Per block of 1024: the sum of the values. */ +extern "C" __global__ void e3_scan_count_vals(long long n, const u32* __restrict__ vals, u32* __restrict__ sums) +{ + long long i = (long long) blockIdx.x * SCAN_BLOCK + threadIdx.x; + u32 total; + block_excl_scan(i < n ? vals[i] : 0u, &total); + if (threadIdx.x == 0) sums[blockIdx.x] = total; +} + +/* One block: the block sums scanned exclusively in place; sums[nb] = the total. */ +extern "C" __global__ void e3_scan_top(int nb, u32* __restrict__ sums) +{ + __shared__ u32 carry; + if (threadIdx.x == 0) carry = 0u; + __syncthreads(); + for (int base = 0; base < nb; base += SCAN_BLOCK) { + int i = base + threadIdx.x; + u32 v = i < nb ? sums[i] : 0u; + u32 total; + u32 ex = block_excl_scan(v, &total); + u32 c = carry; + if (i < nb) sums[i] = c + ex; + __syncthreads(); + if (threadIdx.x == 0) carry = c + total; + __syncthreads(); + } + if (threadIdx.x == 0) sums[nb] = carry; +} + +/* Stream compaction: the flagged indices, ascending, at their scanned offsets. */ +extern "C" __global__ void e3_scan_compact(long long n, const u8* __restrict__ flags, const u32* __restrict__ offsets, u32* __restrict__ out) +{ + long long i = (long long) blockIdx.x * SCAN_BLOCK + threadIdx.x; + u32 v = i < n ? (u32) (flags[i] != 0) : 0u; + u32 total; + u32 ex = block_excl_scan(v, &total); + if (v) out[offsets[blockIdx.x] + ex] = (u32) i; +} + +/* The exclusive scan of the values. */ +extern "C" __global__ void e3_scan_values(long long n, const u32* __restrict__ vals, const u32* __restrict__ offsets, u32* __restrict__ out) +{ + long long i = (long long) blockIdx.x * SCAN_BLOCK + threadIdx.x; + u32 total; + u32 ex = block_excl_scan(i < n ? vals[i] : 0u, &total); + if (i < n) out[i] = offsets[blockIdx.x] + ex; +} + +/* The changed set: touched by either build, dilated by one (periodic z wraps). */ +extern "C" __global__ void e3_mask_changed(MaskGrid g, const u8* __restrict__ t_new, const u8* __restrict__ t_old, u8* __restrict__ flag) +{ + long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long nx = g.nx, ny = g.ny, nz = g.nz, nxy = nx * ny; + if (idx >= nxy * nz) return; + long long k = idx / nxy, j = (idx % nxy) / nx, i = idx % nx; +#define T(q) (t_new[q] | t_old[q]) + u8 f = T(idx) + || (i + 1 < nx && T(idx + 1)) || (i > 0 && T(idx - 1)) + || (j + 1 < ny && T(idx + nx)) || (j > 0 && T(idx - nx)) + || (k + 1 < nz && T(idx + nxy)) || (k > 0 && T(idx - nxy)) + || (g.periodic && nz > 1 && k + 1 == nz && T(idx - (nz - 1) * nxy)) + || (g.periodic && nz > 1 && k == 0 && T(idx + (nz - 1) * nxy)); +#undef T + flag[idx] = f; +} + +/* The faces of the changed cells, component c (a face is its two cells' face). */ +extern "C" __global__ void e3_mask_face_flags(MaskGrid g, int c, const u8* __restrict__ changed, u8* __restrict__ flag) +{ + long long f = (long long) blockIdx.x * blockDim.x + threadIdx.x; + long long ni = g.nx + (c == 0), nj = g.ny + (c == 1), nk = g.nz + (c == 2); + if (f >= ni * nj * nk) return; + long long k = f / (nj * ni), j = (f / ni) % nj, i = f % ni; + long long nx = g.nx, nxy = (long long) g.nx * g.ny; + long long cell = (k * g.ny + j) * nx + i; + long long own = c == 0 ? i : (c == 1 ? j : k); + long long lim = c == 0 ? g.nx : (c == 1 ? g.ny : g.nz); + long long step = c == 0 ? 1 : (c == 1 ? nx : nxy); + u8 m = 0; + if (own < lim) m |= changed[cell]; + if (own > 0) m |= changed[cell - step]; + flag[f] = m; +} + +/* Per changed cell: fluid at the new geometry, the space-time activity (fluid + * at either end), the GCL entry `(V^{n+1} − V^n) dv / dt` of an active cell; + * the device tables written (activity: step and predictor sets). Out flags: + * bit 0 fluid, bit 1 active. */ +extern "C" __global__ void e3_mask_cells( + int n, const u32* __restrict__ list, const double* __restrict__ vol, const double* __restrict__ vol_prev, + double dv, double dt, int* __restrict__ active, int* __restrict__ active_pred, + u8* __restrict__ out_flags, double* __restrict__ out_entry) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + u32 idx = list[t]; + int fn = vol[idx] > 0.0; + int fo = vol_prev[idx] > 0.0; + int act = fn || fo; + active[idx] = act; + active_pred[idx] = fn; + out_flags[t] = (u8) (fn | (act << 1)); + out_entry[t] = act ? (vol[idx] - vol_prev[idx]) * dv / dt : 0.0; +} + +/* lat.cell: the cell at (i, j, k) with the periodic wrap in z; -1 outside. */ +__device__ __forceinline__ long long lat_cell(const MaskGrid& g, long long i, long long j, long long k) +{ + if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return -1; + if (g.periodic) { + k %= g.nz; + if (k < 0) k += g.nz; + } else if (k < 0 || k >= g.nz) { + return -1; + } + return (k * g.ny + j) * g.nx + i; +} + +__device__ __forceinline__ int is_small(long long q, const double* vol, const double* vol_prev, const int* active, double threshold) +{ + return active[q] && fmax(vol[q], vol_prev[q]) < threshold; +} + +/* Per changed cell: the master (the active, not small face neighbour of + * largest new volume, the first on ties) of a small cell. */ +extern "C" __global__ void e3_mask_merge( + MaskGrid g, int n, const u32* __restrict__ list, const double* __restrict__ vol, + const double* __restrict__ vol_prev, const int* __restrict__ active, double threshold, + u32* __restrict__ owner, u32* __restrict__ out_master) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + long long idx = list[t]; + long long nxy = (long long) g.nx * g.ny; + long long k = idx / nxy, j = (idx % nxy) / g.nx, i = idx % g.nx; + u32 master = 0xffffffffu; + if (is_small(idx, vol, vol_prev, active, threshold)) { + int have = 0; + double best = 0.0; + long long bi = 0; + for (int d = 0; d < 3; ++d) { + for (int side = -1; side <= 1; side += 2) { + long long q = lat_cell(g, i + (d == 0) * side, j + (d == 1) * side, k + (d == 2) * side); + if (q < 0) continue; + double v = vol[q]; + if (active[q] && !is_small(q, vol, vol_prev, active, threshold) && (!have || v > best)) { + have = 1; + best = v; + bi = q; + } + } + } + if (have) master = (u32) bi; + } + owner[idx] = master == 0xffffffffu ? (u32) idx : master; + out_master[t] = master; +} + +/* lat.face: the index of the face of component c at (i, j, k), -1 outside + * (periodic z wraps by nz whatever the component). */ +__device__ __forceinline__ long long lat_face(const MaskGrid& g, int c, long long i, long long j, long long k) +{ + long long ni = g.nx + (c == 0), nj = g.ny + (c == 1), nk = g.nz + (c == 2); + if (i < 0 || i >= ni || j < 0 || j >= nj) return -1; + if (g.periodic) { + k %= g.nz; + if (k < 0) k += g.nz; + } else if (k < 0 || k >= nk) { + return -1; + } + if (c == 0) return (k * g.ny + j) * (g.nx + 1) + i; + if (c == 1) return (k * (g.ny + 1) + j) * g.nx + i; + return (k * g.ny + j) * g.nx + i; +} + +struct Apertures { + const double* a[3]; +}; + +/* aperture(c, p).unwrap_or(dflt) */ +__device__ __forceinline__ double ap_or(const MaskGrid& g, const Apertures& A, int c, const long long p[3], double dflt) +{ + long long f = lat_face(g, c, p[0], p[1], p[2]); + return f < 0 ? dflt : A.a[c][f]; +} + +/* Per face of a changed cell, component c: the kind (bit 0: fluid), the step + * aperture ½(αⁿ⁺¹ + αⁿ) and its open flag (bit 1), the open-part centroid + * shift (`compute_face_shifts`, with `has_shift`). */ +extern "C" __global__ void e3_mask_faces( + MaskGrid g, int c, int n, const u32* __restrict__ list, Apertures A, const double* __restrict__ a_prev, + double* __restrict__ a_step, int* __restrict__ open, int* __restrict__ open_pred, int has_shift, + double* __restrict__ shift, u8* __restrict__ out_flags, double* __restrict__ out_a, double* __restrict__ out_shift) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + long long f = list[t]; + 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; + double a_new = A.a[c][f]; + int side; + if (c == 0) side = i == 0 || i == g.nx; + else if (c == 1) side = j == 0 || j == g.ny; + else side = !g.periodic && (k == 0 || k == g.nz); + int fluid = side || a_new > 0.0; + double as = 0.5 * (a_new + a_prev[f]); + int op = as > 0.0; + a_step[f] = as; + open[f] = op; + open_pred[f] = fluid; + out_flags[t] = (u8) (fluid | (op << 1)); + out_a[t] = as; + if (!has_shift) return; + double out[3] = {0.0, 0.0, 0.0}; + long long p[3] = {i, j, k}; + double h[3] = {g.dx, g.dy, g.dz}; + long long nlim[3] = {g.nx, g.ny, g.nz}; + long long fa = lat_face(g, c, p[0], p[1], p[2]); + /* `self.aperture(c, p)` is Some for a face inside the grid */ + double alpha = A.a[c][fa]; + int on_side = p[c] == 0 || p[c] == nlim[c]; + if (alpha > 0.0 && alpha < 1.0 && !(on_side && !(c == 2 && g.periodic))) { + /* cv_geometry(c, p): the control volume's side apertures and wall */ + double area[3] = {g.dy * g.dz, g.dx * g.dz, g.dx * g.dy}; + long long cm[3] = {p[0], p[1], p[2]}; + cm[c] -= 1; + double ap[3][2]; + for (int d = 0; d < 3; ++d) { + if (d == c) { + long long q[3] = {p[0], p[1], p[2]}; + q[c] -= 1; + double am = ap_or(g, A, c, q, alpha); + q[c] += 2; + double apl = ap_or(g, A, c, q, alpha); + ap[d][0] = 0.5 * (am + alpha); + ap[d][1] = 0.5 * (alpha + apl); + } else { + long long cmd[3] = {cm[0], cm[1], cm[2]}; + long long cpd[3] = {p[0], p[1], p[2]}; + cmd[d] += 1; + cpd[d] += 1; + double minus = 0.5 * (ap_or(g, A, d, cm, 1.0) + ap_or(g, A, d, p, 1.0)); + double plus = 0.5 * (ap_or(g, A, d, cmd, 1.0) + ap_or(g, A, d, cpd, 1.0)); + ap[d][0] = minus; + ap[d][1] = plus; + } + } + double nw[3]; + for (int d = 0; d < 3; ++d) nw[d] = -(ap[d][1] - ap[d][0]) * area[d]; + nw[c] = 0.0; + double a = sqrt(nw[0] * nw[0] + nw[1] * nw[1] + nw[2] * nw[2]); + if (a != 0.0) { + for (int d = 0; d < 3; ++d) out[d] = -0.5 * h[d] * (1.0 - alpha) * nw[d] / a; + } + } + for (int d = 0; d < 3; ++d) { + shift[3 * f + d] = out[d]; + out_shift[3 * t + d] = out[d]; + } +} + +/* The GCL table's wall cells: active with a nonzero wall vector. */ +extern "C" __global__ void e3_mask_wall_flags(long long nc, const int* __restrict__ active, const double* __restrict__ wall, u8* __restrict__ flag) +{ + long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= nc) return; + const double* w = wall + 3 * idx; + flag[idx] = active[idx] && (w[0] != 0.0 || w[1] != 0.0 || w[2] != 0.0); +} + +/* The wall areas of the listed cells. */ +extern "C" __global__ void e3_mask_wall_area(int n, const u32* __restrict__ list, const double* __restrict__ wall, double* __restrict__ out) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + const double* w = wall + 3 * (long long) list[t]; + out[t] = sqrt(w[0] * w[0] + w[1] * w[1] + w[2] * w[2]); +} + +/* The GCL wall-flux table: the entries of the changed active cells, less the + * correction per unit wall area on every cell with a wall. */ +extern "C" __global__ void e3_mask_wall_flux( + long long nc, const u8* __restrict__ changed, const int* __restrict__ active, const double* __restrict__ vol, + const double* __restrict__ vol_prev, const double* __restrict__ wall, double dv, double dt, double correction, + double* __restrict__ table) +{ + long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= nc) return; + double v = 0.0; + if (changed[idx] && active[idx]) v = (vol[idx] - vol_prev[idx]) * dv / dt; + if (correction != 0.0) { + const double* w = wall + 3 * idx; + if (w[0] != 0.0 || w[1] != 0.0 || w[2] != 0.0) { + double a = sqrt(w[0] * w[0] + w[1] * w[1] + w[2] * w[2]); + v -= correction * a; + } + } + table[idx] = v; +} + +/* The merged cells: owner != self. */ +extern "C" __global__ void e3_mask_slave_flags(long long nc, const u32* __restrict__ owner, u8* __restrict__ flag) +{ + long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= nc) return; + flag[idx] = owner[idx] != (u32) idx; +} + +/* The slave count per master (counts zeroed by the caller). */ +extern "C" __global__ void e3_mask_fold_count(int n, const u32* __restrict__ slaves, const u32* __restrict__ owner, u32* __restrict__ counts) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + atomicAdd(&counts[owner[slaves[t]]], 1u); +} + +/* The CSR's slave lists, ascending per master: a slave's rank = its master's + * slaves (face neighbours of the master, distinct) of smaller index. */ +extern "C" __global__ void e3_mask_fold_fill( + MaskGrid g, int n, const u32* __restrict__ slaves, const u32* __restrict__ owner, const u32* __restrict__ fold_ptr, + u32* __restrict__ fold_idx) +{ + int t = blockIdx.x * blockDim.x + threadIdx.x; + if (t >= n) return; + long long s = slaves[t]; + long long m = owner[s]; + long long nxy = (long long) g.nx * g.ny; + long long k = m / nxy, j = (m % nxy) / g.nx, i = m % g.nx; + long long seen[6]; + int ns = 0; + u32 rank = 0; + for (int d = 0; d < 3; ++d) { + for (int side = -1; side <= 1; side += 2) { + long long q = lat_cell(g, i + (d == 0) * side, j + (d == 1) * side, k + (d == 2) * side); + if (q < 0 || q == m) continue; + int dup = 0; + for (int r = 0; r < ns; ++r) dup |= seen[r] == q; + if (dup) continue; + seen[ns++] = q; + if (q < s && owner[q] == (u32) m) ++rank; + } + } + fold_idx[fold_ptr[m] + rank] = (u32) s; +} + +/* The imposition band's faces: |d| <= band. */ +extern "C" __global__ void e3_mask_band_flags(long long n, const double* __restrict__ d, double band, u8* __restrict__ flag) +{ + long long f = (long long) blockIdx.x * blockDim.x + threadIdx.x; + if (f >= n) return; + flag[f] = fabs(d[f]) <= band; +} 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 1eb2124..76ef796 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs @@ -13,13 +13,13 @@ //! the computed one as the aperture closes (the shear coefficient grows as //! `1/A_f`), which is what makes the wall smooth in the interface position. -use super::Grid; use super::body::Body; use super::cut::CutGeometry; use super::exchange::in_load_window; use super::field::Field; use super::step::{Boundaries, Side}; use super::wall::{FaceKind, Mask}; +use super::Grid; /// The inertia floor: the momentum volume's fraction in the time /// derivative is at least this. @@ -32,6 +32,12 @@ pub(super) const DISTANCE_FLOOR_FINE: f64 = 0.01; /// step) stays below this shares its pressure unknown with a neighbour. pub(super) const MERGE_FRACTION: f64 = 0.1; +/// `RTX_E3_MOVING_PROFILE` (read once). +pub(super) fn profile_on() -> bool { + static ON: std::sync::OnceLock = std::sync::OnceLock::new(); + *ON.get_or_init(|| std::env::var("RTX_E3_MOVING_PROFILE").is_ok()) +} + /// Lattice addressing of faces and cells with the periodic wrap in z as /// data: a face of component `c` at `p = [i, j, k]` (its own coordinate is /// the face index, the others the cell's), a cell at `[i, j, k]`. @@ -138,6 +144,7 @@ impl Mask { /// 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 lap = std::time::Instant::now(); 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); @@ -148,7 +155,9 @@ impl Mask { 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 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; @@ -182,24 +191,69 @@ impl Mask { .into_par_iter() .map(|f| { let i = f % (nx + 1); - if i == 0 || i == nx { FaceKind::Fluid } else { kind(cut.a_u[f]) } + 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]) } + 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]) } + if !periodic && (k == 0 || k == nz) { + FaceKind::Fluid + } else { + kind(cut.a_w[f]) + } }) .collect(); - let mut mask = Self { + let mut mask = Self::from_parts( + g, + periodic, + cell_fluid, + [u_kind, v_kind, w_kind], + anchor, + fluid_cells, + cut, + ); + let l_class = lap.elapsed(); + mask.compute_merging(None); + if profile_on() { + eprintln!( + " from_cut laps: classification {:.0} ms, merging {:.0} ms", + l_class.as_secs_f64() * 1e3, + (lap.elapsed() - l_class).as_secs_f64() * 1e3 + ); + } + Ok(mask) + } + + /// The mask struct of a cut classification (every closure flag at its + /// default; the caller sets them). + pub(super) fn from_parts( + g: Grid, + periodic: bool, + cell_fluid: Vec, + kinds: [Vec; 3], + anchor: usize, + fluid_cells: usize, + cut: CutGeometry, + ) -> Self { + let [u_kind, v_kind, w_kind] = kinds; + Self { grid: g, periodic_z: periodic, cell_fluid, @@ -234,9 +288,7 @@ impl Mask { grad_weights: None, diffusion_centroid: false, face_shifts: None, - }; - mask.compute_merging(None); - Ok(mask) + } } /// The virtual merging map: a small cell (fraction < `MERGE_FRACTION` @@ -306,6 +358,7 @@ impl Mask { let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { return; }; + let lap = std::time::Instant::now(); let n = inner.len() + 1; let w_end = 0.5 / n as f64; let w_in = 1.0 / n as f64; @@ -336,12 +389,17 @@ impl Mask { || (periodic && nz > 1 && k == 0 && t(idx + (nz - 1) * nxy)) }) .collect(); + let l_changed = lap.elapsed(); // P1-5 (j): with the trapezoid alone and a previous step's arrays, a // face of no changed cell keeps its step aperture (its corners were // untouched in both builds, so αⁿ⁻¹ = αⁿ = αⁿ⁺¹): the old mask's // arrays move over and only the changed cells' faces and the // changed cells' activity are recomputed, with the same expressions. - let prev = if inner.is_empty() { old.step_apertures.take().zip(old.step_open.take()) } else { None }; + let prev = if inner.is_empty() { + old.step_apertures.take().zip(old.step_open.take()) + } else { + None + }; if let Some(((mut au, mut av, mut aw), (mut ou, mut ov, mut ow, mut active))) = prev { let mix = |a: &[f64], b: &[f64], f: usize| w_end * (a[f] + b[f]); for &idx in &changed { @@ -387,7 +445,16 @@ impl Mask { self.step_open = Some((open(&au), open(&av), open(&aw), active)); self.step_apertures = Some((au, av, aw)); } + let l_apert = lap.elapsed(); self.compute_merging(Some(old)); + if profile_on() { + eprintln!( + " step-aperture laps: changed set {:.0} ms, apertures {:.0} ms, merging {:.0} ms", + l_changed.as_secs_f64() * 1e3, + (l_apert - l_changed).as_secs_f64() * 1e3, + (lap.elapsed() - l_apert).as_secs_f64() * 1e3 + ); + } self.changed_cells = Some(changed); } @@ -449,7 +516,11 @@ impl Mask { } // S2-7: the sides' own apertures instead of the whole-face averages. if self.cv_sides_exact { - if let Some(exact) = self.cut.as_ref().and_then(|cut| self.exact_cv_sides(cut, c, p)) { + if let Some(exact) = self + .cut + .as_ref() + .and_then(|cut| self.exact_cv_sides(cut, c, p)) + { ap = exact; } } @@ -715,6 +786,7 @@ impl Mask { let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else { return (table, 0.0); }; + let lap = std::time::Instant::now(); let g = self.grid; let dv = g.dx * g.dy * g.dz; let (mut net, mut area) = (0.0, 0.0); @@ -759,6 +831,7 @@ impl Mask { vn1 - vn ); } + let l_sums = lap.elapsed(); let correction = if area > 0.0 { net / area } else { 0.0 }; if correction != 0.0 { // Every cell with a wall (the others subtract exactly 0). @@ -770,6 +843,13 @@ impl Mask { } } } + if profile_on() { + eprintln!( + " gcl laps: lists + sums {:.0} ms, correction {:.0} ms", + l_sums.as_secs_f64() * 1e3, + (lap.elapsed() - l_sums).as_secs_f64() * 1e3 + ); + } (table, correction) } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs new file mode 100644 index 0000000..3c2b828 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/maskupdate.rs @@ -0,0 +1,373 @@ +//! R6-2 (`RTX_E3_MASK_DEVICE=1`): the host mirror of a moving mask whose +//! classification ran on the device. The device computes, on the changed set +//! of the step (the cells touched by either build, dilated by one) and the +//! faces of those cells, exactly the values the host rebuild computes over +//! the whole grid (`Mask::from_cut`, `compute_face_shifts`, +//! `set_step_apertures`, `compute_merging`, `gcl_flux_table`) and hands them +//! over compactly; the host mask is the previous one's arrays with those +//! entries replaced. Outside the changed set every one of these arrays is +//! unchanged over two steps (their inputs are the apertures and volumes of +//! the touched cells' neighbourhoods), so the instantaneous arrays come from +//! the mask retired one step earlier and the step arrays from the old mask. + +use super::cut::CutGeometry; +use super::step::{Boundaries, Side}; +use super::wall::{FaceKind, Mask}; +use super::Grid; +use rayon::prelude::*; + +/// `RTX_E3_MASK_DEVICE=1` (with `RTX_E3_GEOM_DEVICE=1`, which it requires). +#[must_use] +pub fn mask_device_enabled() -> bool { + let on = |k: &str| std::env::var(k).is_ok_and(|v| v == "1"); + on("RTX_E3_MASK_DEVICE") && on("RTX_E3_GEOM_DEVICE") +} + +/// The device classification of one moving step, compact. +#[derive(Debug, Default)] +pub struct MaskUpdate { + /// The changed cells, ascending. + pub changed: Vec, + /// Per changed cell: bit 0 fluid at the new geometry, bit 1 active over the step. + pub cell_flags: Vec, + /// Per changed cell: the merging master (`u32::MAX` = its own row). + pub master: Vec, + /// The faces of the changed cells per component, ascending. + pub faces: [Vec; 3], + /// Per face: bit 0 kind fluid, bit 1 step-open. + pub face_flags: [Vec; 3], + /// Per face: the step-averaged aperture. + pub face_a: [Vec; 3], + /// Per face: the centroid shift (three interleaved; empty without the centroid diffusion). + pub face_shift: [Vec; 3], + /// The GCL table's compatibility correction (the table itself stays on the device). + pub correction: f64, + /// The merged cells. + pub merged: usize, +} + +/// The instantaneous arrays of a retired mask, kept for the mask two steps later. +#[derive(Default)] +pub struct MaskArrays { + pub(super) generation: u64, + pub(super) cell_fluid: Vec, + pub(super) kinds: [Vec; 3], + pub(super) face_shifts: Option<[Vec; 3]>, +} + +/// A raw pointer that crosses threads for the disjoint writes of a scatter. +struct Shared(*mut T); +impl Clone for Shared { + fn clone(&self) -> Self { + *self + } +} +impl Copy for Shared {} +unsafe impl Send for Shared {} +unsafe impl Sync for Shared {} +impl Shared { + /// SAFETY: the caller writes each in-range index from one thread only. + unsafe fn write(self, i: usize, v: T) { + unsafe { *self.0.add(i) = v }; + } +} + +/// An index of a compact list. +trait Ix: Copy + Sync { + fn ix(self) -> usize; +} +impl Ix for u32 { + fn ix(self) -> usize { + self as usize + } +} +impl Ix for usize { + fn ix(self) -> usize { + self + } +} + +/// `dst[idx[t]] = f(t)` in parallel; `idx` holds distinct in-range indices. +fn scatter(dst: &mut [T], idx: &[I], f: impl Fn(usize) -> T + Sync) { + let n = dst.len(); + assert!( + idx.iter().all(|&i| i.ix() < n), + "scatter index out of range" + ); + let p = Shared(dst.as_mut_ptr()); + idx.par_iter() + .enumerate() + .with_min_len(4096) + .for_each(|(t, &i)| { + // SAFETY: the indices are distinct (ascending compaction output) and in range. + unsafe { p.write(i.ix(), f(t)) }; + }); +} + +fn allowed(side: Side) -> bool { + matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall) +} + +impl Mask { + /// The mask at `cut` from the device classification `upd`, the old + /// mask (generation n − 1: its step arrays and merging map move over) + /// and `pool` (generation n − 2's instantaneous arrays; the old mask's + /// are cloned without it). Classification flags (scheme, closures) are + /// the caller's, as for [`Mask::from_cut`]. + pub(crate) fn from_update( + cut: CutGeometry, + g: Grid, + b: Boundaries, + upd: &MaskUpdate, + old: &mut Mask, + pool: Option, + ) -> Result { + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let nxy = nx * ny; + let periodic = b.z0 == Side::Periodic; + let (mut cell_fluid, kinds, mut face_shifts) = match pool { + Some(p) => (p.cell_fluid, p.kinds, p.face_shifts), + None => ( + old.cell_fluid.clone(), + [old.u_kind.clone(), old.v_kind.clone(), old.w_kind.clone()], + old.face_shifts.clone(), + ), + }; + let [mut u_kind, mut v_kind, mut w_kind] = kinds; + let ch = &upd.changed; + scatter(&mut cell_fluid, ch, |t| upd.cell_flags[t] & 1 != 0); + // The touching check: only a changed cell can newly reach a side. + let touches = |idx: usize| { + 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 ch.iter().any(|&idx| !cell_fluid[idx] && touches(idx)) { + let first = (0..g.cells()) + .into_par_iter() + .find_first(|&idx| !cell_fluid[idx] && touches(idx)) + .expect("touching cell"); + let (k, j, i) = (first / nxy, (first % nxy) / nx, first % nx); + return Err(format!( + "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})" + )); + } + // The fluid count and the anchor (the smallest fluid index). + let delta: i64 = ch + .iter() + .map(|&idx| i64::from(cell_fluid[idx]) - i64::from(old.cell_fluid[idx])) + .sum(); + let fluid_cells = (old.fluid_cells as i64 + delta) as usize; + let anchor = if cell_fluid[old.anchor] { + match ch.iter().find(|&&idx| cell_fluid[idx]) { + Some(&idx) if idx < old.anchor => Some(idx), + _ => Some(old.anchor), + } + } else { + (0..g.cells()) + .into_par_iter() + .find_first(|&idx| cell_fluid[idx]) + }; + let Some(anchor) = anchor else { + return Err("embedded body covers the whole domain".into()); + }; + // The faces: kinds, the step arrays (moved from the old mask), the shifts. + let (Some((mut au, mut av, mut aw)), Some((mut ou, mut ov, mut ow, mut active))) = + (old.step_apertures.take(), old.step_open.take()) + else { + return Err("R6-2: the old mask has no step arrays".into()); + }; + let kind_of = |bits: u8| { + if bits & 1 != 0 { + FaceKind::Fluid + } else { + FaceKind::Solid + } + }; + { + let kinds: [&mut Vec; 3] = [&mut u_kind, &mut v_kind, &mut w_kind]; + let steps: [&mut Vec; 3] = [&mut au, &mut av, &mut aw]; + let opens: [&mut Vec; 3] = [&mut ou, &mut ov, &mut ow]; + for (c, ((kd, st), op)) in kinds.into_iter().zip(steps).zip(opens).enumerate() { + let (list, flags, a) = (&upd.faces[c], &upd.face_flags[c], &upd.face_a[c]); + scatter(kd, list.as_slice(), |t| kind_of(flags[t])); + scatter(st, list.as_slice(), |t| a[t]); + scatter(op, list.as_slice(), |t| flags[t] & 2 != 0); + } + } + scatter(&mut active, ch, |t| upd.cell_flags[t] & 2 != 0); + if let Some(tables) = face_shifts.as_mut() { + for c in 0..3 { + let (list, sh) = (&upd.faces[c], &upd.face_shift[c]); + assert_eq!(sh.len(), 3 * list.len(), "R6-2: face shifts missing"); + let tab = &mut tables[c]; + let n = tab.len(); + assert!( + list.iter().all(|&f| 3 * f as usize + 2 < n), + "shift index out of range" + ); + let p = Shared(tab.as_mut_ptr()); + list.par_iter() + .enumerate() + .with_min_len(4096) + .for_each(|(t, &f)| { + let f = f as usize; + // SAFETY: distinct faces, three entries each, in range. + unsafe { + p.write(3 * f, sh[3 * t]); + p.write(3 * f + 1, sh[3 * t + 1]); + p.write(3 * f + 2, sh[3 * t + 2]); + } + }); + } + } + let mut merge_master = std::mem::take(&mut old.merge_master); + if merge_master.len() != g.cells() { + return Err("R6-2: the old mask has no merging map".into()); + } + scatter(&mut merge_master, ch, |t| match upd.master[t] { + u32::MAX => usize::MAX, + m => m as usize, + }); + let mut mask = Self::from_parts( + g, + periodic, + cell_fluid, + [u_kind, v_kind, w_kind], + anchor, + fluid_cells, + cut, + ); + mask.step_apertures = Some((au, av, aw)); + mask.step_open = Some((ou, ov, ow, active)); + mask.merge_master = merge_master; + mask.face_shifts = face_shifts; + mask.changed_cells = Some(upd.changed.clone()); + Ok(mask) + } + + /// The instantaneous arrays of this (retiring) mask for the pool. + pub(crate) fn into_pool(&mut self, generation: u64, cell_fluid: Vec) -> MaskArrays { + MaskArrays { + generation, + cell_fluid, + kinds: [ + std::mem::take(&mut self.u_kind), + std::mem::take(&mut self.v_kind), + std::mem::take(&mut self.w_kind), + ], + face_shifts: self.face_shifts.take(), + } + } + + /// `RTX_E3_BAND_CHECK=1`: this mask's classification arrays against the + /// host rebuild's `reference`, bit for bit; the differences named. + pub(crate) fn compare_classification(&self, reference: &Mask) -> Vec { + let mut report = Vec::new(); + let mut eq = |name: &str, ok: bool| { + if !ok { + report.push(format!("{name} differs")); + } + }; + let bits = |a: &[f64], b: &[f64]| { + a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()) + }; + eq("cell_fluid", self.cell_fluid == reference.cell_fluid); + eq("u_kind", self.u_kind == reference.u_kind); + eq("v_kind", self.v_kind == reference.v_kind); + eq("w_kind", self.w_kind == reference.w_kind); + eq("anchor", self.anchor == reference.anchor); + eq("fluid_cells", self.fluid_cells == reference.fluid_cells); + eq("merge_master", self.merge_master == reference.merge_master); + eq( + "changed_cells", + self.changed_cells == reference.changed_cells, + ); + match (&self.step_apertures, &reference.step_apertures) { + (Some(a), Some(r)) => { + eq("step a_u", bits(&a.0, &r.0)); + eq("step a_v", bits(&a.1, &r.1)); + eq("step a_w", bits(&a.2, &r.2)); + } + (a, r) => eq("step apertures presence", a.is_some() == r.is_some()), + } + match (&self.step_open, &reference.step_open) { + (Some(a), Some(r)) => { + eq("step open u", a.0 == r.0); + eq("step open v", a.1 == r.1); + eq("step open w", a.2 == r.2); + eq("step active", a.3 == r.3); + } + (a, r) => eq("step open presence", a.is_some() == r.is_some()), + } + match (&self.face_shifts, &reference.face_shifts) { + (Some(a), Some(r)) => { + for c in 0..3 { + eq(&format!("face shifts[{c}]"), bits(&a[c], &r[c])); + } + } + (a, r) => eq("face shifts presence", a.is_some() == r.is_some()), + } + report + } +} + +/// Refill the pressure of the cells fluid in `new` and not in `old` from +/// their face neighbours fluid in both; returns their count — the host +/// `refill_fresh_cells` over the changed cells (a fresh cell is touched), +/// in the same ascending order. +pub(crate) fn refill_fresh_in(old: &Mask, new: &Mask, p: &mut [f64], changed: &[usize]) -> usize { + let g = new.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let periodic = new.periodic_z; + let mut fresh = 0; + let mut refills = Vec::new(); + for &idx in changed { + if !(new.cell_fluid[idx] && !old.cell_fluid[idx]) { + continue; + } + fresh += 1; + let (k, j, i) = g.kji(idx); + let mut sum = 0.0; + let mut count = 0usize; + let mut visit = |nb: usize| { + if new.cell_fluid[nb] && old.cell_fluid[nb] { + sum += p[nb]; + count += 1; + } + }; + if i + 1 < nx { + visit(g.cell(k, j, i + 1)); + } + if i > 0 { + visit(g.cell(k, j, i - 1)); + } + if j + 1 < ny { + visit(g.cell(k, j + 1, i)); + } + if j > 0 { + visit(g.cell(k, j - 1, i)); + } + if k + 1 < nz { + visit(g.cell(k + 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(0, j, i)); + } + if k > 0 { + visit(g.cell(k - 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(nz - 1, j, i)); + } + if count > 0 { + refills.push((idx, sum / count as f64)); + } + } + for (idx, v) in refills { + p[idx] = v; + } + fresh +} 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 e04393d..5f55ba2 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs @@ -16,6 +16,7 @@ pub mod field; pub mod grid; pub mod impose; pub mod loads; +pub mod maskupdate; pub mod poisson; pub mod reconstruct; pub mod step; 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 acc1359..34d508f 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 @@ -7,12 +7,13 @@ mod cut; mod geom; +mod mask; use super::{Side, Solver, StepResult}; -use crate::solvers::incompressible::embedded3::Grid; use crate::solvers::incompressible::embedded3::field::Field; use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime}; use crate::solvers::incompressible::embedded3::poisson::device_cg::DeviceCg; +use crate::solvers::incompressible::embedded3::Grid; use crate::solvers::incompressible::poisson::MultigridParameters; use crate::solvers::incompressible::simple::ConvectionScheme; use cudarc::driver::{ @@ -165,6 +166,8 @@ pub struct DeviceStep { steps_since_hierarchy: usize, /// R6-1: the persistent device cut geometry (`RTX_E3_GEOM_DEVICE=1`). geom: Option, + /// R6-2: the device classification (`RTX_E3_MASK_DEVICE=1`). + dmask: Option, } impl DeviceStep { @@ -244,6 +247,7 @@ impl DeviceStep { cut, steps_since_hierarchy: 0, geom: None, + dmask: 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 2c1687c..167f21e 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 @@ -5,13 +5,13 @@ //! kernels (`e3_cut.cu`, appended to `e3_step.cu` at load). use super::{DeviceStep, E3Params, E3Ptrs, StepResult}; -use crate::solvers::incompressible::embedded3::Grid; use crate::solvers::incompressible::embedded3::field::Field; -use crate::solvers::incompressible::embedded3::poisson::GuessBasis; use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime}; use crate::solvers::incompressible::embedded3::poisson::device_cg::DeviceCg; +use crate::solvers::incompressible::embedded3::poisson::GuessBasis; use crate::solvers::incompressible::embedded3::step::Solver; use crate::solvers::incompressible::embedded3::wall::FaceKind; +use crate::solvers::incompressible::embedded3::Grid; use crate::solvers::incompressible::poisson::MultigridParameters; use cudarc::driver::{ CudaFunction, CudaModule, CudaSlice, DevicePtr, DeviceRepr, PushKernelArg, ValidAsZeroBits, @@ -64,27 +64,27 @@ unsafe impl ValidAsZeroBits for E3CutPtrs {} /// The static cut-cell mask on the device. pub(super) struct DeviceCut { /// The surface-velocity tables (kept to share between the phases). - ub_host: [Vec; 3], + pub(super) ub_host: [Vec; 3], /// Which phase's apertures / open flags / activity `ptrs` hands out /// (P1-5: both sets live in the struct and persist across steps). - phase: Phase, + pub(super) phase: Phase, /// The predictor phase's instantaneous apertures, open flags (fluid /// kind) and activity (fluid cells). - a_pred: [CudaSlice; 3], - open_pred: [CudaSlice; 3], - active_pred: CudaSlice, - a: [CudaSlice; 3], - d: [CudaSlice; 3], - ub: [CudaSlice; 3], - wall_flux: CudaSlice, - open: [CudaSlice; 3], - active: CudaSlice, - owner: CudaSlice, - fold_ptr: CudaSlice, - fold_idx: CudaSlice, + pub(super) a_pred: [CudaSlice; 3], + pub(super) open_pred: [CudaSlice; 3], + pub(super) active_pred: CudaSlice, + pub(super) a: [CudaSlice; 3], + pub(super) d: [CudaSlice; 3], + pub(super) ub: [CudaSlice; 3], + pub(super) wall_flux: CudaSlice, + pub(super) open: [CudaSlice; 3], + pub(super) active: CudaSlice, + pub(super) owner: CudaSlice, + pub(super) fold_ptr: CudaSlice, + pub(super) fold_idx: CudaSlice, cell_flux: CudaSlice, /// The open-part centroid shifts per face (S2-5; one dummy entry when off). - shift: [CudaSlice; 3], + pub(super) shift: [CudaSlice; 3], /// The wall's normal velocity into the fluid per face (A3-i; one dummy entry when off). vn: [CudaSlice; 3], /// The wall velocity at the two axis feet along each direction per face, @@ -234,7 +234,8 @@ impl DeviceCut { let delta = mask.exchange_delta(&cv, d); let mut xf = x; xf[d] += sign * delta; - out[2 * d + side] = mask.surface_velocity_at(body, xf, c, t); + out[2 * d + side] = + mask.surface_velocity_at(body, xf, c, t); } } } @@ -333,11 +334,18 @@ impl DeviceCut { }) .collect() }; - let active_pred: Vec = (0..nc).into_par_iter().map(|i| i32::from(mask.is_fluid_cell(i))).collect(); + let active_pred: Vec = (0..nc) + .into_par_iter() + .map(|i| i32::from(mask.is_fluid_cell(i))) + .collect(); let (a_pred, open_pred_dev, active_pred_dev) = if projection { ( [up_f(&cut.a_u), up_f(&cut.a_v), up_f(&cut.a_w)], - [up_i(&open_pred(0)), up_i(&open_pred(1)), up_i(&open_pred(2))], + [ + up_i(&open_pred(0)), + up_i(&open_pred(1)), + up_i(&open_pred(2)), + ], up_i(&active_pred), ) } else { @@ -410,12 +418,20 @@ impl DeviceCut { /// global). `false` when a full build is needed (no changed set, the /// advancing-wall closure, a size change). pub(super) fn update(&mut self, solver: &Solver, g: Grid, t: f64) -> bool { - use crate::solvers::incompressible::embedded3::poisson::device_cg::{scatter_f64, scatter_i32, scatter_u32}; + use crate::solvers::incompressible::embedded3::poisson::device_cg::{ + scatter_f64, scatter_i32, scatter_u32, + }; use rayon::prelude::*; - let Some(mask) = solver.mask() else { return false }; + let Some(mask) = solver.mask() else { + return false; + }; let Some(cut) = mask.cut() else { return false }; - let Some(body) = solver.body() else { return false }; - let Some(changed) = mask.changed_cells() else { return false }; + let Some(body) = solver.body() else { + return false; + }; + let Some(changed) = mask.changed_cells() else { + return false; + }; if mask.wall_advancing || mask.wall_exchange_foot { return false; } @@ -429,7 +445,11 @@ impl DeviceCut { let lap = Instant::now(); let profile = std::env::var("RTX_E3_MOVING_PROFILE").is_ok(); // The faces of the changed cells, per component, ascending. - let mut marks: [Vec; 3] = [vec![false; counts[0]], vec![false; counts[1]], vec![false; counts[2]]]; + let mut marks: [Vec; 3] = [ + vec![false; counts[0]], + vec![false; counts[1]], + vec![false; counts[2]], + ]; for &idx in changed { let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx); marks[0][g.uface(k, j, i)] = true; @@ -440,19 +460,30 @@ impl DeviceCut { marks[2][g.wface(k + 1, j, i)] = true; } let touched: [Vec; 3] = [0, 1, 2].map(|c| { - (0..counts[c]).into_par_iter().filter(|&f| marks[c][f]).map(|f| f as u32).collect() + (0..counts[c]) + .into_par_iter() + .filter(|&f| marks[c][f]) + .map(|f| f as u32) + .collect() }); let cells_u32: Vec = changed.iter().map(|&i| i as u32).collect(); + let l_touched = lap.elapsed(); // The imposition band's faces (their surface velocity moves with t). let band = mask.impose_band().unwrap_or(f64::INFINITY); let dists: [&[f64]; 3] = [&cut.d_u, &cut.d_v, &cut.d_w]; let band_faces: [Vec; 3] = [0, 1, 2].map(|c| { - (0..counts[c]).into_par_iter().filter(|&f| dists[c][f].abs() <= band).map(|f| f as u32).collect() + (0..counts[c]) + .into_par_iter() + .filter(|&f| dists[c][f].abs() <= band) + .map(|f| f as u32) + .collect() }); let h = [g.dx, g.dy, g.dz]; let shifts = mask.face_shift_tables(); let step = mask.step_apertures(); let step_open = mask.step_open_flags(); + let l_band = lap.elapsed(); + let mut l_ub = std::time::Duration::ZERO; for c in 0..3 { let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1)); let pos = |idx: usize| -> [f64; 3] { @@ -463,7 +494,11 @@ impl DeviceCut { (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2], ] }; - let a_inst: &[f64] = match c { 0 => &cut.a_u, 1 => &cut.a_v, _ => &cut.a_w }; + let a_inst: &[f64] = match c { + 0 => &cut.a_u, + 1 => &cut.a_v, + _ => &cut.a_w, + }; let a_step: &[f64] = match (step, c) { (Some(a), 0) => &a.0, (Some(a), 1) => &a.1, @@ -477,7 +512,12 @@ impl DeviceCut { 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(); + let open_step: Vec = tf + .iter() + .map(|&f| { + i32::from(step_open.map_or(a_inst[f as usize] > 0.0, |o| o[c][f as usize])) + }) + .collect(); scatter_i32(tf, &open_step, &mut self.open[c]); let open_inst: Vec = tf .iter() @@ -492,18 +532,35 @@ impl DeviceCut { .collect(); scatter_i32(tf, &open_inst, &mut self.open_pred[c]); if let Some(sh) = shifts { - let idx3: Vec = tf.iter().flat_map(|&f| [3 * f, 3 * f + 1, 3 * f + 2]).collect(); - let val3: Vec = tf.iter().flat_map(|&f| [sh[c][3 * f as usize], sh[c][3 * f as usize + 1], sh[c][3 * f as usize + 2]]).collect(); + let idx3: Vec = tf + .iter() + .flat_map(|&f| [3 * f, 3 * f + 1, 3 * f + 2]) + .collect(); + let val3: Vec = tf + .iter() + .flat_map(|&f| { + [ + sh[c][3 * f as usize], + sh[c][3 * f as usize + 1], + sh[c][3 * f as usize + 2], + ] + }) + .collect(); scatter_f64(&idx3, &val3, &mut self.shift[c]); } // The surface velocity at the foot for the band's faces. + let lu = Instant::now(); let bf = &band_faces[c]; let ubv: Vec = bf .par_iter() .map(|&f| { let x = pos(f as usize); let xf = match (shifts, mask.wall_foot_centroid) { - (Some(sh), true) => [x[0] + sh[c][3 * f as usize], x[1] + sh[c][3 * f as usize + 1], x[2] + sh[c][3 * f as usize + 2]], + (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) @@ -513,15 +570,26 @@ impl DeviceCut { self.ub_host[c][f as usize] = v; } scatter_f64(bf, &ubv, &mut self.ub[c]); + l_ub += lu.elapsed(); } let l_faces = lap.elapsed(); // Per-cell tables for the changed cells. - let active_step: Vec = changed.iter().map(|&i| i32::from(mask.cell_active(i))).collect(); - let active_inst: Vec = changed.iter().map(|&i| i32::from(mask.is_fluid_cell(i))).collect(); - let owner: Vec = changed.iter().map(|&i| mask.master(i).unwrap_or(i) as u32).collect(); + let active_step: Vec = changed + .iter() + .map(|&i| i32::from(mask.cell_active(i))) + .collect(); + let active_inst: Vec = changed + .iter() + .map(|&i| i32::from(mask.is_fluid_cell(i))) + .collect(); + let owner: Vec = changed + .iter() + .map(|&i| mask.master(i).unwrap_or(i) as u32) + .collect(); scatter_i32(&cells_u32, &active_step, &mut self.active); scatter_i32(&cells_u32, &active_inst, &mut self.active_pred); scatter_u32(&cells_u32, &owner, &mut self.owner); + let l_owner = lap.elapsed(); // The wall flux: the changed cells and every wall cell (the correction). let wall_flux: &[f64] = solver.wall_fluxes(); if wall_flux.len() != nc { @@ -538,6 +606,7 @@ impl DeviceCut { .collect(); let wvals: Vec = wcells.iter().map(|&i| wall_flux[i as usize]).collect(); scatter_f64(&wcells, &wvals, &mut self.wall_flux); + let l_wflux = lap.elapsed(); // The merged-cell CSR whole (its lengths change). let mut fold_ptr = vec![0u32; nc + 1]; let mut merged = 0; @@ -584,6 +653,16 @@ impl DeviceCut { ms(lap.elapsed() - l_faces), changed.len() ); + eprintln!( + " update laps: touched lists {:.0} ms, band scan {:.0} ms, face scatters {:.0} ms, ub {:.0} ms, cell scatters {:.0} ms, wall flux {:.0} ms, fold CSR {:.0} ms", + ms(l_touched), + ms(l_band - l_touched), + ms(l_faces - l_band - l_ub), + ms(l_ub), + ms(l_owner - l_faces), + ms(l_wflux - l_owner), + ms(lap.elapsed() - l_wflux) + ); } true } @@ -600,15 +679,27 @@ impl DeviceCut { pub(super) fn check_against(&self, full: &Self) { let rt = runtime(); let same_f = |name: &str, a: &CudaSlice, b: &CudaSlice| { - let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); - assert!(x.len() == y.len() && x.iter().zip(&y).all(|(p, q)| p.to_bits() == q.to_bits()), "band tables: {name} differs"); + let (x, y) = ( + rt.stream.memcpy_dtov(a).expect("dl"), + rt.stream.memcpy_dtov(b).expect("dl"), + ); + assert!( + x.len() == y.len() && x.iter().zip(&y).all(|(p, q)| p.to_bits() == q.to_bits()), + "band tables: {name} differs" + ); }; let same_i = |name: &str, a: &CudaSlice, b: &CudaSlice| { - let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); + let (x, y) = ( + rt.stream.memcpy_dtov(a).expect("dl"), + rt.stream.memcpy_dtov(b).expect("dl"), + ); assert!(x == y, "band tables: {name} differs"); }; let same_u = |name: &str, a: &CudaSlice, b: &CudaSlice| { - let (x, y) = (rt.stream.memcpy_dtov(a).expect("dl"), rt.stream.memcpy_dtov(b).expect("dl")); + let (x, y) = ( + rt.stream.memcpy_dtov(a).expect("dl"), + rt.stream.memcpy_dtov(b).expect("dl"), + ); assert!(x == y, "band tables: {name} differs"); }; for c in 0..3 { @@ -619,8 +710,15 @@ impl DeviceCut { same_i("open_pred", &self.open_pred[c], &full.open_pred[c]); same_f("shift", &self.shift[c], &full.shift[c]); // ub: the band's faces only (beyond it the full build writes 0, ours keeps stale — never read). - let (x, y) = (rt.stream.memcpy_dtov(&self.ub[c]).expect("dl"), rt.stream.memcpy_dtov(&full.ub[c]).expect("dl")); - let bad = x.iter().zip(&y).filter(|(p, q)| **q != 0.0 && p.to_bits() != q.to_bits()).count(); + let (x, y) = ( + rt.stream.memcpy_dtov(&self.ub[c]).expect("dl"), + rt.stream.memcpy_dtov(&full.ub[c]).expect("dl"), + ); + let bad = x + .iter() + .zip(&y) + .filter(|(p, q)| **q != 0.0 && p.to_bits() != q.to_bits()) + .count(); assert!(bad == 0, "band tables: ub[{c}] differs on {bad} band faces"); } same_f("wall_flux", &self.wall_flux, &full.wall_flux); @@ -749,17 +847,52 @@ impl DeviceStep { // 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. + let mut device_mask: Option = None; if super::geom::enabled() { if let Some(dc) = self.cut.as_mut() { + let fresh_geom = self.geom.is_none(); let geom = self .geom .get_or_insert_with(|| super::geom::DeviceGeom::new(g)); + // R6-2: the previous step's apertures, volumes and + // evaluated cells, before the geometry kernels overwrite them. + let dm = if super::mask::enabled() && !fresh_geom { + let dm = self + .dmask + .get_or_insert_with(|| super::mask::DeviceMask::new(g)); + dm.snapshot(dc, geom); + Some(dm) + } else { + None + }; 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 let Some(dm) = dm { + if dm.ready(&self.solver, geom) { + let lm = Instant::now(); + let upd = dm.run(&self.solver, g, dt, dc, geom); + device_mask = Some(upd.merged); + if profile { + let l = dm.laps; + eprintln!( + " moving laps: device mask {:.0} ms (lists {:.0}, cells + merging {:.0}, faces {:.0}, GCL {:.0}, fold {:.0}; {} changed, {} merged)", + lm.elapsed().as_secs_f64() * 1e3, + l[0], + l[1], + l[2], + l[3], + l[4], + upd.changed.len(), + upd.merged + ); + } + self.solver.pending_mask = Some(upd); + } + } } } if profile { @@ -781,10 +914,25 @@ impl DeviceStep { let rebuild_all = std::env::var("RTX_E3_TABLES_REBUILD").is_ok(); let check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); self.cut = match self.cut.take() { + // 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)); + if check { + let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new) + .expect("full"); + prev.check_against(&full); + eprintln!(" table check t {t_new:.6}: device tables IDENTICAL to the full build"); + } + Some(prev) + } Some(mut prev) if !rebuild_all => { if prev.update(&self.solver, g, t_new) { if check { - let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new).expect("full"); + let full = DeviceCut::build(&self.solver, g, Phase::Projection, t_new) + .expect("full"); prev.check_against(&full); } Some(prev) @@ -829,7 +977,10 @@ impl DeviceStep { smoother: self.solver.params.poisson_smoother, ..MultigridParameters::default() }; - let changed: Option> = self.solver.mask().and_then(|m| m.changed_cells().map(|c| c.to_vec())); + let changed: Option> = self + .solver + .mask() + .and_then(|m| m.changed_cells().map(|c| c.to_vec())); cg.refresh_with(problem, ¶ms, changed.as_deref()); } rt.stream 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 3eb1b5b..39fcdd6 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 @@ -10,11 +10,11 @@ //! 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 crate::solvers::incompressible::embedded3::Grid; use cudarc::driver::{ CudaFunction, CudaModule, CudaSlice, DeviceRepr, PushKernelArg, ValidAsZeroBits, }; @@ -106,6 +106,8 @@ pub(super) struct DeviceGeom { generation: u64, chain_start: u64, history: std::collections::VecDeque<(u64, [Vec; 3], Vec)>, + /// R6-2: the last build was a full pass (the device state was re-synced). + pub(super) last_full: bool, } /// Generations of band lists kept (a recycled array older than this is copied instead). @@ -131,9 +133,16 @@ impl DeviceGeom { generation: 0, chain_start: 1, history: std::collections::VecDeque::new(), + last_full: true, } } + /// R6-2: the cell volumes, wall vectors (3 interleaved) and evaluated + /// cells of the last build. + pub(super) fn cell_tables(&self) -> (&CudaSlice, &CudaSlice, &CudaSlice) { + (&self.vol, &self.wall, &self.touched_cell) + } + /// 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 @@ -212,6 +221,7 @@ impl DeviceGeom { None => (0i32, 0.0, 0.0), }; let full_i = i32::from(full); + self.last_full = full; unsafe { rt.stream .launch_builder(&k.phi) 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 new file mode 100644 index 0000000..57221b2 --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/mask.rs @@ -0,0 +1,699 @@ +//! R6-2: the moving mask's classification on the device (`e3_mask.cu`), +//! behind `RTX_E3_MASK_DEVICE=1` (which needs `RTX_E3_GEOM_DEVICE=1`). After +//! the device geometry of the step (R6-1), the changed set, the faces of the +//! changed cells, the kinds, the step apertures and open flags, the space-time +//! activity, the merging masters, the centroid shifts, the GCL wall-flux +//! table and the merged cells' CSR are computed on the device and written +//! straight into `DeviceCut`'s persistent tables; the host receives the same +//! values compactly (`MaskUpdate`) for its mirror. The previous step's +//! apertures, volumes and evaluated cells are snapshotted before the geometry +//! kernels overwrite them. + +use super::cut::{DeviceCut, Phase}; +use super::geom::DeviceGeom; +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}; +use crate::solvers::incompressible::embedded3::step::Solver; +use crate::solvers::incompressible::embedded3::wall::WallScheme; +use crate::solvers::incompressible::embedded3::Grid; +use cudarc::driver::{ + CudaFunction, CudaModule, CudaSlice, DevicePtr, DeviceRepr, LaunchConfig, PushKernelArg, + ValidAsZeroBits, +}; +use std::sync::{Arc, OnceLock}; +use std::time::Instant; + +const MASK_KERNELS: &str = include_str!("../../../../../kernels/cuda/e3_mask.cu"); +const SCAN_BLOCK: usize = 1024; + +struct MaskKernels { + _module: Arc, + count_flags: CudaFunction, + count_vals: CudaFunction, + top: CudaFunction, + compact: CudaFunction, + values: CudaFunction, + changed: CudaFunction, + face_flags: CudaFunction, + cells: CudaFunction, + merge: CudaFunction, + faces: CudaFunction, + wall_flags: CudaFunction, + wall_area: CudaFunction, + wall_flux: CudaFunction, + slave_flags: CudaFunction, + fold_count: CudaFunction, + fold_fill: CudaFunction, + band_flags: CudaFunction, +} + +static MASK_ONCE: OnceLock = OnceLock::new(); + +fn kernels() -> &'static MaskKernels { + MASK_ONCE.get_or_init(|| { + // FMA contraction off: the host's roundings. + let module = load_module(MASK_KERNELS, "e3_mask.cu", true); + let f = |name: &str| module.load_function(name).expect(name); + MaskKernels { + count_flags: f("e3_scan_count_flags"), + count_vals: f("e3_scan_count_vals"), + top: f("e3_scan_top"), + compact: f("e3_scan_compact"), + values: f("e3_scan_values"), + changed: f("e3_mask_changed"), + face_flags: f("e3_mask_face_flags"), + cells: f("e3_mask_cells"), + merge: f("e3_mask_merge"), + faces: f("e3_mask_faces"), + wall_flags: f("e3_mask_wall_flags"), + wall_area: f("e3_mask_wall_area"), + wall_flux: f("e3_mask_wall_flux"), + slave_flags: f("e3_mask_slave_flags"), + fold_count: f("e3_mask_fold_count"), + fold_fill: f("e3_mask_fold_fill"), + band_flags: f("e3_mask_band_flags"), + _module: module, + } + }) +} + +/// `struct MaskGrid` in e3_mask.cu. +#[repr(C)] +#[derive(Clone, Copy)] +struct MaskGrid { + nx: i32, + ny: i32, + nz: i32, + periodic: i32, + dx: f64, + dy: f64, + dz: f64, +} +unsafe impl DeviceRepr for MaskGrid {} +unsafe impl ValidAsZeroBits for MaskGrid {} + +/// `struct Apertures` in e3_mask.cu: the three instantaneous aperture tables. +#[repr(C)] +#[derive(Clone, Copy)] +struct Apertures { + a: [u64; 3], +} +unsafe impl DeviceRepr for Apertures {} +unsafe impl ValidAsZeroBits for Apertures {} + +/// `RTX_E3_MASK_DEVICE=1` with `RTX_E3_GEOM_DEVICE=1`. +pub(super) fn enabled() -> bool { + crate::solvers::incompressible::embedded3::maskupdate::mask_device_enabled() +} + +fn scan_cfg(n: usize) -> LaunchConfig { + LaunchConfig { + grid_dim: (n.div_ceil(SCAN_BLOCK).max(1) as u32, 1, 1), + block_dim: (SCAN_BLOCK as u32, 1, 1), + shared_mem_bytes: 0, + } +} + +/// A device buffer grown on demand (never shrunk). +fn ensure(buf: &mut Option>, n: usize) { + let n = n.max(1); + if buf.as_ref().is_none_or(|b| b.len() < n) { + let cap = n + n / 4; + *buf = Some(runtime().stream.alloc_zeros::(cap).expect("alloc")); + } +} + +/// Download the first `n` entries. +fn head(buf: &CudaSlice, n: usize) -> Vec { + if n == 0 { + return Vec::new(); + } + runtime() + .stream + .memcpy_dtov(&buf.slice(0..n)) + .expect("download") +} + +/// The persistent state of the device classification. +pub(super) struct DeviceMask { + /// The previous step's instantaneous apertures, volumes and evaluated cells. + a_prev: [CudaSlice; 3], + vol_prev: CudaSlice, + touched_prev: CudaSlice, + /// The snapshot was taken this step (before the geometry kernels). + snapped: bool, + changed_flag: CudaSlice, + flags: CudaSlice, + sums: CudaSlice, + counts: CudaSlice, + changed_list: Option>, + face_list: [Option>; 3], + list_tmp: Option>, + out_cell_flags: Option>, + out_entry: Option>, + out_master: Option>, + out_face_flags: [Option>; 3], + out_face_a: [Option>; 3], + out_face_shift: [Option>; 3], + out_wall_a: Option>, + /// The laps of the last run (ms): lists, cells + merging, faces, GCL, fold. + pub(super) laps: [f64; 6], +} + +impl DeviceMask { + pub(super) fn new(g: Grid) -> Self { + let rt = runtime(); + let nc = g.cells(); + let counts = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()]; + let nmax = nc.max(counts[0]).max(counts[1]).max(counts[2]) + 1; + let z8 = |n: usize| rt.stream.alloc_zeros::(n).expect("alloc"); + Self { + a_prev: counts.map(|n| rt.stream.alloc_zeros::(n).expect("alloc")), + vol_prev: rt.stream.alloc_zeros::(nc).expect("alloc"), + touched_prev: z8(nc), + snapped: false, + changed_flag: z8(nc), + flags: z8(nmax), + sums: rt + .stream + .alloc_zeros::(nmax.div_ceil(SCAN_BLOCK) + 1) + .expect("alloc"), + counts: rt.stream.alloc_zeros::(nc).expect("alloc"), + changed_list: None, + face_list: [None, None, None], + list_tmp: None, + out_cell_flags: None, + out_entry: None, + out_master: None, + out_face_flags: [None, None, None], + out_face_a: [None, None, None], + out_face_shift: [None, None, None], + out_wall_a: None, + laps: [0.0; 6], + } + } + + /// Before the geometry kernels: the current instantaneous apertures, + /// volumes and evaluated cells become the previous step's. + pub(super) fn snapshot(&mut self, dc: &DeviceCut, geom: &DeviceGeom) { + let rt = runtime(); + let (vol, _, touched) = geom.cell_tables(); + for c in 0..3 { + rt.stream + .memcpy_dtod(&dc.a_pred[c], &mut self.a_prev[c]) + .expect("a_prev"); + } + rt.stream + .memcpy_dtod(vol, &mut self.vol_prev) + .expect("vol_prev"); + rt.stream + .memcpy_dtod(touched, &mut self.touched_prev) + .expect("touched_prev"); + self.snapped = true; + } + + /// Whether the classification can run on the device this step: the + /// snapshot is of the geometry the host holds as the previous one (no + /// re-sync), the old mask is a moving one (step arrays, merging map), + /// and no closure the device does not carry is on. + pub(super) fn ready(&self, solver: &Solver, geom: &DeviceGeom) -> bool { + let p = &solver.params; + let Some(m) = solver.mask() else { return false }; + self.snapped + && !geom.last_full + && p.wall_scheme == WallScheme::CutCell + && p.aperture_substeps == 0 + && !p.wall_advancing + && !p.wall_exchange_foot + && !p.pressure_centroid + && std::env::var("RTX_E3_TABLES_REBUILD").is_err() + && m.step_apertures().is_some() + && m.step_open_flags().is_some() + && m.changed_cells().is_some() + && m.face_shift_tables().is_some() == p.diffusion_centroid + } + + /// Compact the first `n` flags of `self.flags` into `out` (ascending); the count. + fn compact(&mut self, n: usize, out: &mut Option>) -> usize { + let rt = runtime(); + let k = kernels(); + let nb = n.div_ceil(SCAN_BLOCK).max(1); + let n64 = n as i64; + let nb32 = nb as i32; + unsafe { + rt.stream + .launch_builder(&k.count_flags) + .arg(&n64) + .arg(&self.flags) + .arg(&mut self.sums) + .launch(scan_cfg(n)) + .expect("e3_scan_count_flags"); + rt.stream + .launch_builder(&k.top) + .arg(&nb32) + .arg(&mut self.sums) + .launch(LaunchConfig { + grid_dim: (1, 1, 1), + block_dim: (SCAN_BLOCK as u32, 1, 1), + shared_mem_bytes: 0, + }) + .expect("e3_scan_top"); + } + let total = rt + .stream + .memcpy_dtov(&self.sums.slice(nb..nb + 1)) + .expect("count")[0] as usize; + ensure(out, total); + let dst = out.as_mut().expect("list"); + unsafe { + rt.stream + .launch_builder(&k.compact) + .arg(&n64) + .arg(&self.flags) + .arg(&self.sums) + .arg(dst) + .launch(scan_cfg(n)) + .expect("e3_scan_compact"); + } + total + } + + /// The device classification of the step (the geometry just built, the + /// snapshot the previous one): `dc`'s projection tables written, the + /// compact values returned for the host mirror. + pub(super) fn run( + &mut self, + solver: &Solver, + g: Grid, + dt: f64, + dc: &mut DeviceCut, + geom: &DeviceGeom, + ) -> MaskUpdate { + let rt = runtime(); + let k = kernels(); + let lap = Instant::now(); + self.snapped = false; + let nc = g.cells(); + let counts = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()]; + let mg = MaskGrid { + nx: g.nx as i32, + ny: g.ny as i32, + nz: g.nz as i32, + periodic: i32::from(solver.params.boundaries.periodic_z()), + dx: g.dx, + dy: g.dy, + dz: g.dz, + }; + let (vol, wall, touched) = geom.cell_tables(); + // The changed set. + unsafe { + rt.stream + .launch_builder(&k.changed) + .arg(&mg) + .arg(touched) + .arg(&self.touched_prev) + .arg(&mut self.changed_flag) + .launch(cfg(nc)) + .expect("e3_mask_changed"); + } + rt.stream + .memcpy_dtod(&self.changed_flag, &mut self.flags.slice_mut(0..nc)) + .expect("flags"); + let mut changed_list = self.changed_list.take(); + let n_ch = self.compact(nc, &mut changed_list); + self.changed_list = changed_list; + // The faces of the changed cells. + let mut n_faces = [0usize; 3]; + for c in 0..3 { + let ci = c as i32; + unsafe { + rt.stream + .launch_builder(&k.face_flags) + .arg(&mg) + .arg(&ci) + .arg(&self.changed_flag) + .arg(&mut self.flags) + .launch(cfg(counts[c])) + .expect("e3_mask_face_flags"); + } + let mut fl = self.face_list[c].take(); + n_faces[c] = self.compact(counts[c], &mut fl); + self.face_list[c] = fl; + } + let ch_list = self.changed_list.as_ref().expect("changed"); + let changed: Vec = head(ch_list, n_ch) + .into_iter() + .map(|i| i as usize) + .collect(); + let faces: [Vec; 3] = + [0, 1, 2].map(|c| head(self.face_list[c].as_ref().expect("faces"), n_faces[c])); + let l_lists = lap.elapsed(); + // Per changed cell: classification, activity, GCL entry; then the masters. + ensure(&mut self.out_cell_flags, n_ch); + ensure(&mut self.out_entry, n_ch); + ensure(&mut self.out_master, n_ch); + let dv = g.dx * g.dy * g.dz; + let n_ch32 = n_ch as i32; + let threshold = std::env::var("RTX_E3_MERGE_FRACTION") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(MERGE_FRACTION); + if n_ch > 0 { + unsafe { + rt.stream + .launch_builder(&k.cells) + .arg(&n_ch32) + .arg(ch_list) + .arg(vol) + .arg(&self.vol_prev) + .arg(&dv) + .arg(&dt) + .arg(&mut dc.active) + .arg(&mut dc.active_pred) + .arg(self.out_cell_flags.as_mut().expect("out")) + .arg(self.out_entry.as_mut().expect("out")) + .launch(cfg(n_ch)) + .expect("e3_mask_cells"); + rt.stream + .launch_builder(&k.merge) + .arg(&mg) + .arg(&n_ch32) + .arg(ch_list) + .arg(vol) + .arg(&self.vol_prev) + .arg(&dc.active) + .arg(&threshold) + .arg(&mut dc.owner) + .arg(self.out_master.as_mut().expect("out")) + .launch(cfg(n_ch)) + .expect("e3_mask_merge"); + } + } + let cell_flags = head(self.out_cell_flags.as_ref().expect("out"), n_ch); + let entries = head(self.out_entry.as_ref().expect("out"), n_ch); + let master = head(self.out_master.as_ref().expect("out"), n_ch); + let l_cells = lap.elapsed(); + // Per face of a changed cell: kind, step aperture and open flag, shift. + let s = &rt.stream; + let ap = Apertures { + a: [0, 1, 2].map(|c| dc.a_pred[c].device_ptr(s).0), + }; + let has_shift = + i32::from(solver.params.diffusion_centroid && dc.shift[0].len() == 3 * counts[0]); + let mut face_flags: [Vec; 3] = Default::default(); + let mut face_a: [Vec; 3] = Default::default(); + let mut face_shift: [Vec; 3] = Default::default(); + for c in 0..3 { + let n = n_faces[c]; + ensure(&mut self.out_face_flags[c], n); + ensure(&mut self.out_face_a[c], n); + ensure(&mut self.out_face_shift[c], 3 * n); + if n > 0 { + let (ci, n32) = (c as i32, n as i32); + unsafe { + rt.stream + .launch_builder(&k.faces) + .arg(&mg) + .arg(&ci) + .arg(&n32) + .arg(self.face_list[c].as_ref().expect("faces")) + .arg(&ap) + .arg(&self.a_prev[c]) + .arg(&mut dc.a[c]) + .arg(&mut dc.open[c]) + .arg(&mut dc.open_pred[c]) + .arg(&has_shift) + .arg(&mut dc.shift[c]) + .arg(self.out_face_flags[c].as_mut().expect("out")) + .arg(self.out_face_a[c].as_mut().expect("out")) + .arg(self.out_face_shift[c].as_mut().expect("out")) + .launch(cfg(n)) + .expect("e3_mask_faces"); + } + } + face_flags[c] = head(self.out_face_flags[c].as_ref().expect("out"), n); + face_a[c] = head(self.out_face_a[c].as_ref().expect("out"), n); + if has_shift != 0 { + face_shift[c] = head(self.out_face_shift[c].as_ref().expect("out"), 3 * n); + } + } + let l_faces = lap.elapsed(); + // The GCL table: the entries' and the wall areas' sums in ascending + // order on the host (the serial loops' order), the table on the device. + let nc64 = nc as i64; + unsafe { + rt.stream + .launch_builder(&k.wall_flags) + .arg(&nc64) + .arg(&dc.active) + .arg(wall) + .arg(&mut self.flags) + .launch(cfg(nc)) + .expect("e3_mask_wall_flags"); + } + let mut tmp = self.list_tmp.take(); + let n_wall = self.compact(nc, &mut tmp); + self.list_tmp = tmp; + ensure(&mut self.out_wall_a, n_wall); + if n_wall > 0 { + let nw32 = n_wall as i32; + unsafe { + rt.stream + .launch_builder(&k.wall_area) + .arg(&nw32) + .arg(self.list_tmp.as_ref().expect("list")) + .arg(wall) + .arg(self.out_wall_a.as_mut().expect("out")) + .launch(cfg(n_wall)) + .expect("e3_mask_wall_area"); + } + } + let areas = head(self.out_wall_a.as_ref().expect("out"), n_wall); + let mut net = 0.0; + for (t, &e) in entries.iter().enumerate() { + if cell_flags[t] & 2 != 0 { + net += e; + } + } + let mut area = 0.0; + for &a in &areas { + area += a; + } + let correction = if area > 0.0 { net / area } else { 0.0 }; + unsafe { + rt.stream + .launch_builder(&k.wall_flux) + .arg(&nc64) + .arg(&self.changed_flag) + .arg(&dc.active) + .arg(vol) + .arg(&self.vol_prev) + .arg(wall) + .arg(&dv) + .arg(&dt) + .arg(&correction) + .arg(&mut dc.wall_flux) + .launch(cfg(nc)) + .expect("e3_mask_wall_flux"); + } + let l_gcl = lap.elapsed(); + // The merged cells' CSR. + unsafe { + rt.stream + .launch_builder(&k.slave_flags) + .arg(&nc64) + .arg(&dc.owner) + .arg(&mut self.flags) + .launch(cfg(nc)) + .expect("e3_mask_slave_flags"); + } + let mut tmp = self.list_tmp.take(); + let merged = self.compact(nc, &mut tmp); + self.list_tmp = tmp; + let slaves = self.list_tmp.as_ref().expect("slaves"); + rt.stream.memset_zeros(&mut self.counts).expect("counts"); + if dc.fold_ptr.len() != nc + 1 { + dc.fold_ptr = rt.stream.alloc_zeros::(nc + 1).expect("alloc"); + } + if dc.fold_idx.len() != merged.max(1) { + dc.fold_idx = rt.stream.alloc_zeros::(merged.max(1)).expect("alloc"); + } + let m32 = merged as i32; + let nb = nc.div_ceil(SCAN_BLOCK).max(1); + let nb32 = nb as i32; + unsafe { + if merged > 0 { + rt.stream + .launch_builder(&k.fold_count) + .arg(&m32) + .arg(slaves) + .arg(&dc.owner) + .arg(&mut self.counts) + .launch(cfg(merged)) + .expect("e3_mask_fold_count"); + } + rt.stream + .launch_builder(&k.count_vals) + .arg(&nc64) + .arg(&self.counts) + .arg(&mut self.sums) + .launch(scan_cfg(nc)) + .expect("e3_scan_count_vals"); + rt.stream + .launch_builder(&k.top) + .arg(&nb32) + .arg(&mut self.sums) + .launch(LaunchConfig { + grid_dim: (1, 1, 1), + block_dim: (SCAN_BLOCK as u32, 1, 1), + shared_mem_bytes: 0, + }) + .expect("e3_scan_top"); + rt.stream + .launch_builder(&k.values) + .arg(&nc64) + .arg(&self.counts) + .arg(&self.sums) + .arg(&mut dc.fold_ptr) + .launch(scan_cfg(nc)) + .expect("e3_scan_values"); + } + rt.stream + .memcpy_htod(&[merged as u32][..], &mut dc.fold_ptr.slice_mut(nc..nc + 1)) + .expect("fold_ptr end"); + if merged > 0 { + unsafe { + rt.stream + .launch_builder(&k.fold_fill) + .arg(&mg) + .arg(&m32) + .arg(slaves) + .arg(&dc.owner) + .arg(&dc.fold_ptr) + .arg(&mut dc.fold_idx) + .launch(cfg(merged)) + .expect("e3_mask_fold_fill"); + } + } + rt.stream.synchronize().expect("sync"); + let l_fold = lap.elapsed(); + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + self.laps = [ + ms(l_lists), + ms(l_cells - l_lists), + ms(l_faces - l_cells), + ms(l_gcl - l_faces), + ms(l_fold - l_gcl), + ms(l_fold), + ]; + MaskUpdate { + changed, + cell_flags, + master, + faces, + face_flags, + face_a, + face_shift, + correction, + merged, + } + } + + /// 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] { + 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) + }) + } +} + +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. + pub(super) fn update_after_device_mask( + &mut self, + solver: &Solver, + g: Grid, + t: f64, + dm: &mut DeviceMask, + merged: usize, + ) -> bool { + use crate::solvers::incompressible::embedded3::poisson::device_cg::scatter_f64; + use rayon::prelude::*; + let (Some(mask), Some(body)) = (solver.mask(), solver.body()) else { + return false; + }; + 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 (nx, ny) = (g.nx, g.ny); + for (c, bf) in band_faces.iter().enumerate() { + let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1)); + let pos = |idx: usize| -> [f64; 3] { + let (k, j, i) = (idx / (nj * ni), (idx / ni) % nj, idx % ni); + [ + (i as f64 + if c == 0 { 0.0 } else { 0.5 }) * h[0], + (j as f64 + if c == 1 { 0.0 } else { 0.5 }) * h[1], + (k as f64 + if c == 2 { 0.0 } else { 0.5 }) * h[2], + ] + }; + let ubv: Vec = bf + .par_iter() + .map(|&f| { + let x = pos(f as usize); + let xf = match (shifts, mask.wall_foot_centroid) { + (Some(sh), true) => [ + x[0] + sh[c][3 * f as usize], + x[1] + sh[c][3 * f as usize + 1], + x[2] + sh[c][3 * f as usize + 2], + ], + _ => x, + }; + mask.surface_velocity_at(body, xf, c, t) + }) + .collect(); + for (&f, &v) in bf.iter().zip(&ubv) { + self.ub_host[c][f as usize] = v; + } + scatter_f64(bf, &ubv, &mut self.ub[c]); + } + self.merged = merged; + self.phase = Phase::Projection; + self.geom_written = false; + 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)", + 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() + ); + } + true + } +} 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 ca230af..47bfbfa 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 @@ -11,12 +11,12 @@ mod moving; mod predictor; mod projection; -use super::Grid; use super::body::Body; use super::cut::CutGeometry; use super::field::Field; -use super::poisson::{PcgCache, Problem, solve_pcg_cached}; +use super::poisson::{solve_pcg_cached, PcgCache, Problem}; use super::wall::{FaceKind, Mask, WallScheme}; +use super::Grid; use crate::solvers::incompressible::poisson::{ MgPrecision, MgSmoother, MultigridParameters, PoissonSolution, }; @@ -304,6 +304,14 @@ pub struct Solver { /// and the generation of `apertures_old`. pub(super) geom_pool: super::cut::GeomPool, pub(super) apertures_old_gen: u64, + /// R6-2: the next moving rebuild's classification, computed on the + /// device (`RTX_E3_MASK_DEVICE=1`); `rebuild_moving_mask` assembles the + /// mask from it instead of classifying the whole grid. + pub(super) pending_mask: Option, + /// R6-2: the instantaneous arrays of the mask retired at the previous + /// rebuild, and the rebuild counter (the mask generation). + pub(super) mask_pool: Option, + pub(super) mask_gen: u64, } impl Solver { @@ -340,6 +348,9 @@ impl Solver { pending_cut: None, geom_pool: super::cut::GeomPool::default(), apertures_old_gen: 0, + pending_mask: None, + mask_pool: None, + mask_gen: 0, } } @@ -403,23 +414,7 @@ impl Solver { } .map(|mut m| { let lap_closures = std::time::Instant::now(); - m.scheme = self.params.convection_scheme; - m.density = self.fluid.density; - m.wall_order = self.params.wall_order; - m.wall_distance_oblique = self.params.wall_distance_oblique; - m.wall_exchange_axis = self.params.wall_exchange_axis; - m.diffusion_transverse = - self.params.diffusion_transverse && self.params.diffusion_centroid; - m.distance_floor_fine = self.params.distance_floor_fine; - m.wall_advancing = self.params.wall_advancing; - m.exchange_convection_off = self.params.exchange_convection_off; - m.cv_sides_exact = self.params.cv_sides_exact; - m.wall_order2_centroid = self.params.wall_order2_centroid; - m.wall_exchange_foot = self.params.wall_exchange_foot; - m.conv_sides_exact = self.params.conv_sides_exact; - m.wall_flux_true_normal = self.params.wall_flux_true_normal; - m.wall_foot_centroid = self.params.wall_foot_centroid; - m.diffusion_centroid = self.params.diffusion_centroid; + self.configure_mask(&mut m); if self.params.diffusion_centroid { m.compute_face_shifts(); } @@ -437,6 +432,26 @@ impl Solver { .expect("embedded mask") } + /// The solver's closure flags on a freshly classified mask. + pub(super) fn configure_mask(&self, m: &mut Mask) { + m.scheme = self.params.convection_scheme; + m.density = self.fluid.density; + m.wall_order = self.params.wall_order; + m.wall_distance_oblique = self.params.wall_distance_oblique; + m.wall_exchange_axis = self.params.wall_exchange_axis; + m.diffusion_transverse = self.params.diffusion_transverse && self.params.diffusion_centroid; + m.distance_floor_fine = self.params.distance_floor_fine; + m.wall_advancing = self.params.wall_advancing; + m.exchange_convection_off = self.params.exchange_convection_off; + m.cv_sides_exact = self.params.cv_sides_exact; + m.wall_order2_centroid = self.params.wall_order2_centroid; + m.wall_exchange_foot = self.params.wall_exchange_foot; + m.conv_sides_exact = self.params.conv_sides_exact; + m.wall_flux_true_normal = self.params.wall_flux_true_normal; + m.wall_foot_centroid = self.params.wall_foot_centroid; + m.diffusion_centroid = self.params.diffusion_centroid; + } + /// R5 instrument: record the cut predictor's momentum terms per face /// from the next step on (host path only; a no-op on the device path). pub fn enable_term_probe(&self) { 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 7f5565e..1c96e5e 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,6 +14,9 @@ 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 { + if let Some(upd) = self.pending_mask.take() { + return self.rebuild_from_update(field, dt, t_new, upd); + } let pending = self.pending_cut.take(); let Some(body) = &self.body else { return 0; @@ -22,8 +25,11 @@ impl Solver { let lap = std::time::Instant::now(); let mut new_mask = self.build_mask(body, field.grid, t_new, dt, pending); let l_build = lap.elapsed(); + let mut l_refill = std::time::Duration::ZERO; if let Some(old_mask) = self.mask.as_mut() { + let lr = std::time::Instant::now(); fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); + l_refill = lr.elapsed(); let n_in = self.params.aperture_substeps; if n_in == 0 { new_mask.set_step_apertures(old_mask); @@ -73,64 +79,205 @@ impl Solver { self.last_ghost_correction = correction; } let s_gcl = sub.elapsed(); - // P1-5 (a): the old mask is replaced below — its arrays MOVE into - // the previous-step records instead of being cloned (1 GB per step - // at ny 124); the values are the same. - 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() { - 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() - .zip(fluid.par_iter()) - .for_each(|(v, &f)| { - if !f { - *v = 0.0; - } - }); - 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()]; - } + if let Some(old) = self.mask.take() { + self.retire_mask(old, field.grid); } + self.mask_gen += 1; if std::env::var("RTX_E3_MOVING_PROFILE").is_ok() { let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; eprintln!( - " mask laps: build_mask {:.0} ms, refill + step apertures + merging {:.0} ms, impose + GCL + volumes {:.0} ms (impose {:.0}, GCL table {:.0}, old apertures + volumes {:.0})", + " mask laps: build_mask {:.0} ms, refill + step apertures + merging {:.0} ms, impose + GCL + volumes {:.0} ms (impose {:.0}, GCL table {:.0}, old apertures + volumes {:.0}); refill {:.0} ms", ms(l_build), ms(l_step - l_build), ms(lap.elapsed() - l_step), ms(s_impose), ms(s_gcl - s_impose), - ms(sub.elapsed() - s_gcl) + ms(sub.elapsed() - s_gcl), + ms(l_refill) ); } self.mask = Some(new_mask); fresh_cells } + + /// R6-2: the moving rebuild from the device classification `upd` (the + /// geometry's mirror in `pending_cut`): the mask assembled from the old + /// one's arrays and the compact device values, the fresh cells refilled + /// over the changed set, the imposition as on the host path; the GCL + /// table stays on the device (its correction comes back). With + /// `RTX_E3_BAND_CHECK=1` the host path runs on copies and every array is + /// compared bit for bit. + fn rebuild_from_update( + &mut self, + field: &mut Field, + dt: f64, + t_new: f64, + upd: super::super::maskupdate::MaskUpdate, + ) -> usize { + use super::super::maskupdate::refill_fresh_in; + let g = field.grid; + let b = self.params.boundaries; + let cut = self + .pending_cut + .take() + .expect("R6-2: the device geometry's mirror"); + let mut old_mask = self.mask.take().expect("R6-2: the previous mask"); + let lap = std::time::Instant::now(); + let check = std::env::var("RTX_E3_BAND_CHECK").is_ok(); + // The reference: the host rebuild on copies. + let reference = check.then(|| { + let mut r = Mask::from_cut(cut.clone(), g, b).expect("reference mask"); + self.configure_mask(&mut r); + if self.params.diffusion_centroid { + r.compute_face_shifts(); + } + let mut p_ref = field.p.clone(); + let fresh_ref = refill_fresh_cells_p(&old_mask, &r, g, &mut p_ref); + let mut old_copy = old_mask.clone(); + r.set_step_apertures(&mut old_copy); + let (table, correction) = r.gcl_flux_table(&old_mask, dt); + (r, p_ref, fresh_ref, table, correction) + }); + let l_ref = lap.elapsed(); + let generation = self.mask_gen + 1; + let pool = self + .mask_pool + .take() + .filter(|p| p.generation + 2 == generation); + let pooled = pool.is_some(); + let mut new_mask = + Mask::from_update(cut, g, b, &upd, &mut old_mask, pool).expect("embedded mask"); + self.configure_mask(&mut new_mask); + let l_build = lap.elapsed(); + 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, + ); + let s_impose = sub.elapsed(); + self.last_ghost_correction = upd.correction; + if let Some((r, p_ref, fresh_ref, table, correction)) = reference { + let mut report = new_mask.compare_classification(&r); + if correction.to_bits() != upd.correction.to_bits() { + report.push(format!( + "GCL correction differs ({correction:e} vs {:e})", + upd.correction + )); + } + if fresh_ref != fresh_cells + || p_ref + .iter() + .zip(&field.p) + .any(|(a, b)| a.to_bits() != b.to_bits()) + { + report.push("fresh-cell refill differs".into()); + } + // The device tables' full reference build reads the host table. + self.wall_fluxes = table; + if report.is_empty() { + eprintln!( + " mask check t {t_new:.6}: device classification IDENTICAL to the host rebuild ({} changed cells, {} / {} / {} faces, pool {})", + upd.changed.len(), + upd.faces[0].len(), + upd.faces[1].len(), + upd.faces[2].len(), + pooled + ); + } else { + eprintln!( + " mask check t {t_new:.6}: DIFFERS — {}", + report.join("; ") + ); + if std::env::var("RTX_E3_GEOM_CHECK_SOFT").is_err() { + panic!("mask check: device classification differs from the host rebuild"); + } + } + } + let s_gcl = sub.elapsed(); + self.retire_mask(old_mask, g); + self.mask_gen = generation; + if std::env::var("RTX_E3_MOVING_PROFILE").is_ok() { + let ms = |d: std::time::Duration| d.as_secs_f64() * 1e3; + eprintln!( + " mask laps: build_mask {:.0} ms, refill + step apertures + merging {:.0} ms, impose + GCL + volumes {:.0} ms (impose {:.0}, GCL table {:.0}, old apertures + volumes {:.0}); refill {:.0} ms; device-mask assembly (pool {}), check reference {:.0} ms", + ms(l_build - l_ref), + ms(l_step - l_build), + ms(lap.elapsed() - l_step), + ms(s_impose), + ms(s_gcl - s_impose), + ms(sub.elapsed() - s_gcl), + ms(l_step - l_build), + pooled, + ms(l_ref) + ); + } + self.mask = Some(new_mask); + fresh_cells + } + /// Retire the previous mask: its apertures and volumes MOVE into the + /// previous-step records (P1-5 (a): no 1 GB clone per step at ny 124), + /// its geometry arrays to the pool. + fn retire_mask(&mut self, mut old: Mask, g: super::super::Grid) { + // P1-5 (a): the old mask is replaced below — its arrays MOVE into + // the previous-step records instead of being cloned (1 GB per step + // at ny 124); the values are the same. + let fluid = std::mem::take(&mut old.cell_fluid); + let keep = super::super::maskupdate::mask_device_enabled(); + if let Some(c) = old.cut.as_mut() { + 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() + .zip(fluid.par_iter()) + .for_each(|(v, &f)| { + if !f { + *v = 0.0; + } + }); + 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; g.cells()]; + } + // R6-2: the instantaneous arrays wait for the mask two steps on. + if keep { + self.mask_pool = Some(old.into_pool(self.mask_gen, fluid)); + } + } } /// Refill the pressure of the cells fluid in `new` and not in `old` from /// their face neighbours fluid in both; returns their count. fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { let g = field.grid; + refill_fresh_cells_p(old, new, g, &mut field.p) +} + +fn refill_fresh_cells_p(old: &Mask, new: &Mask, g: super::super::Grid, p: &mut [f64]) -> usize { let (nx, ny, nz) = (g.nx, g.ny, g.nz); let periodic = new.periodic_z(); let mut fresh = 0; @@ -147,7 +294,7 @@ fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { let mut count = 0usize; let mut visit = |nb: usize| { if new.is_fluid_cell(nb) && old.is_fluid_cell(nb) { - sum += field.p[nb]; + sum += p[nb]; count += 1; } }; @@ -179,8 +326,8 @@ fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { } } } - for (idx, p) in refills { - field.p[idx] = p; + for (idx, v) in refills { + p[idx] = v; } fresh }