R8-f: the operator-route wall load on the device (RTX_E3_LOADS_DEVICE=1, default off)
e3_loads.cu evaluates every summand of Mask::cut_wall_force (cells' p W, fluid faces' implicit wall shear) and cut_wall_exchange_parts (diffusive and convective exchange with prescribed neighbours) on the device from the predictor set of the device mask, in the host's expression order with FMA off; the summands are appended as records, downloaded once per pass, sorted into the host's loop order and summed sequentially on the host, so the totals are the host route's to the bit. Records replay into the R8-a load sink (exact zeros dropped; sum-invariant) and into R8-c's WallLoad shape. DeviceStep::cut_wall_loads_device / cut_wall_force_device / check_wall_loads (the G1 instrument, RTX_E3_LOADS_CHECK=1 in the flag test). The flag test and the R8-a harness take the device route under the knob; default paths untouched. cut.rs: E3CutPtrs and DeviceCut::ptrs pub(super). Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
co-authored by
Claude Opus 5.5
Claude Opus 5.5
parent
171da41ed1
commit
2d5ffbabdd
@@ -0,0 +1,237 @@
|
||||
/**
|
||||
* R8-f: the operator-route wall load on the device — `Mask::cut_wall_force`
|
||||
* (the cells' p W and the fluid faces' implicit wall shear, cutwall.rs) and
|
||||
* `Mask::cut_wall_exchange_parts` (the diffusive and convective exchange of
|
||||
* the fluid faces with their prescribed neighbours, exchange.rs), summand
|
||||
* for summand in the host's expression order (FMA contraction off, so each
|
||||
* summand is the host's to the bit). The kernels do not sum: every summand
|
||||
* is appended as a record (key, index, value) to one buffer; the host sorts
|
||||
* the records into its loop order and sums them sequentially, so the totals
|
||||
* are the host's to the bit as well.
|
||||
*
|
||||
* Appended after `e3_step.cu` and `e3_cut.cu` at load (shares E3Params,
|
||||
* E3Ptrs, E3Cut and `cut_face` / `cut_cv` / `face_corr3` / `upwind3`).
|
||||
* The E3Cut pointers are the PREDICTOR set: the instantaneous apertures, the
|
||||
* fluid-kind open flags and the fluid cells of the mask at the solver's time.
|
||||
*/
|
||||
|
||||
/* key = kind << 28 | c << 24 | slot; kind 0 the cell's p W (v[3]), 1 the
|
||||
face's wall shear, 2 an exchange summand (slot = 4 d + 2 side + sub; side
|
||||
0 plus, 1 minus; sub 0 convective, 1 diffusive). */
|
||||
struct LoadRec {
|
||||
unsigned int key;
|
||||
unsigned int idx;
|
||||
double v[3];
|
||||
};
|
||||
|
||||
struct LoadArgs {
|
||||
double mu, rho;
|
||||
int k0, k1; /* the z planes */
|
||||
int shifts; /* the centroid-shift tables are the mask's (face_shifts Some) */
|
||||
int centroid; /* diffusion_centroid */
|
||||
int axis; /* wall_exchange_axis */
|
||||
int conv_off; /* exchange_convection_off */
|
||||
int order; /* the shear closure's order (wall_order) */
|
||||
unsigned int cap; /* the record buffer's capacity */
|
||||
};
|
||||
|
||||
__device__ __forceinline__ void load_push(LoadRec* out, unsigned int* count, unsigned int cap,
|
||||
unsigned int key, unsigned int idx, double v0, double v1, double v2)
|
||||
{
|
||||
unsigned int at = atomicAdd(count, 1u);
|
||||
if (at < cap) {
|
||||
out[at].key = key;
|
||||
out[at].idx = idx;
|
||||
out[at].v[0] = v0;
|
||||
out[at].v[1] = v1;
|
||||
out[at].v[2] = v2;
|
||||
}
|
||||
}
|
||||
|
||||
/* The cells' p W over the planes k0..k1 (thread per cell of those planes):
|
||||
W = −Σ A n from the apertures (cut.rs, `[-sx, -sy, -sz]`); a fluid cell
|
||||
with W = 0 adds exact zeros to every sum and is not recorded. */
|
||||
extern "C" __global__ void e3_loads_cells(E3Params g, E3Ptrs f, E3Cut m, LoadArgs a,
|
||||
LoadRec* __restrict__ out, unsigned int* __restrict__ count)
|
||||
{
|
||||
long long t = (long long)blockIdx.x * blockDim.x + threadIdx.x;
|
||||
long long plane = (long long)g.nx * g.ny;
|
||||
if (t >= plane * (a.k1 - a.k0)) return;
|
||||
long long idx = (long long)a.k0 * plane + t;
|
||||
if (!m.active[idx]) return;
|
||||
int i = (int)(idx % g.nx), j = (int)((idx / g.nx) % g.ny), k = (int)(idx / plane);
|
||||
double ax = g.dy * g.dz, ay = g.dx * g.dz, az = g.dx * g.dy;
|
||||
double sx = (m.a_u[uf3(g, k, j, i + 1)] - m.a_u[uf3(g, k, j, i)]) * ax;
|
||||
double sy = (m.a_v[vf3(g, k, j + 1, i)] - m.a_v[vf3(g, k, j, i)]) * ay;
|
||||
double sz = (m.a_w[wf3(g, k + 1, j, i)] - m.a_w[wf3(g, k, j, i)]) * az;
|
||||
double w0 = -sx, w1 = -sy, w2 = -sz;
|
||||
if (w0 == 0.0 && w1 == 0.0 && w2 == 0.0) return;
|
||||
double p = f.p[idx];
|
||||
load_push(out, count, a.cap, 0u, (unsigned int)idx, p * w0, p * w1, p * w2);
|
||||
}
|
||||
|
||||
__device__ __forceinline__ const double* load_val(const E3Ptrs& f, int c) { return c == 0 ? f.u : (c == 1 ? f.v : f.w); }
|
||||
__device__ __forceinline__ const int* load_open(const E3Cut& m, int c) { return c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w); }
|
||||
|
||||
/* The fluid faces of component c (thread per face of the component's
|
||||
lattice): the implicit wall shear (cutwall.rs) and the exchange with the
|
||||
prescribed neighbours (exchange.rs). */
|
||||
extern "C" __global__ void e3_loads_faces(E3Params g, E3Ptrs f, E3Cut m, LoadArgs a, int c,
|
||||
LoadRec* __restrict__ out, unsigned int* __restrict__ count)
|
||||
{
|
||||
long long t = (long long)blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int ni = c == 0 ? g.nx + 1 : g.nx;
|
||||
int nj = c == 1 ? g.ny + 1 : g.ny;
|
||||
int nk = c == 2 ? g.nz + 1 : g.nz;
|
||||
if (t >= (long long)ni * nj * nk) return;
|
||||
int i = (int)(t % ni), j = (int)((t / ni) % nj), k = (int)(t / ((long long)ni * nj));
|
||||
/* the host's loop ranges */
|
||||
if (c == 0 && (i < 1 || i >= g.nx)) return;
|
||||
if (c == 1 && (j < 1 || j >= g.ny)) return;
|
||||
int ks = a.k0, ke = a.k1;
|
||||
if (c == 2) {
|
||||
int ws = g.periodic_z ? 0 : 1;
|
||||
if (ws > ks) ks = ws;
|
||||
if (g.nz < ke) ke = g.nz;
|
||||
}
|
||||
if (k < ks || k >= ke) return;
|
||||
int fidx = (int)t;
|
||||
const int* open = load_open(m, c);
|
||||
if (!open[fidx]) return;
|
||||
const double* vals = load_val(f, c);
|
||||
double h[3] = { g.dx, g.dy, g.dz };
|
||||
double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy };
|
||||
double alpha, apm[3], app[3], wall[3], distance;
|
||||
cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &distance);
|
||||
double u0 = vals[fidx];
|
||||
unsigned int kc = (unsigned int)c << 24;
|
||||
|
||||
/* The implicit wall shear: mu A_w (c1 (u − U_b) + c2 (u_n − U_b)). */
|
||||
double a_w = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||||
if (a_w != 0.0) {
|
||||
const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
|
||||
double ub = ubt[fidx];
|
||||
/* advancing_factor(0) / d = (1 + 0.5 * 0) / d */
|
||||
double c1 = (1.0 + 0.5 * 0.0) / distance, c2 = 0.0;
|
||||
int nb = -1;
|
||||
if (a.order >= 2) {
|
||||
double n[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w };
|
||||
int d = 0;
|
||||
for (int kk = 1; kk < 3; ++kk) if (fabs(n[kk]) > fabs(n[d])) d = kk;
|
||||
int q[3] = { i, j, k };
|
||||
q[d] -= n[d] > 0.0 ? 1 : -1;
|
||||
int fn = cut_face(g, c, q[0], q[1], q[2]);
|
||||
if (fn >= 0 && cut_ap(m, c)[fn] > 0.0) {
|
||||
double d1 = distance;
|
||||
double d2 = d1 + h[d] * fabs(n[d]);
|
||||
if (!(d2 <= d1)) {
|
||||
c1 = d2 / (d1 * (d2 - d1));
|
||||
c2 = -d1 / (d2 * (d2 - d1));
|
||||
nb = fn;
|
||||
}
|
||||
}
|
||||
}
|
||||
double un = nb >= 0 ? vals[nb] : ub;
|
||||
double v = a.mu * a_w * (c1 * (u0 - ub) + c2 * (un - ub));
|
||||
load_push(out, count, a.cap, (1u << 28) | kc, (unsigned int)fidx,
|
||||
c == 0 ? v : 0.0, c == 1 ? v : 0.0, c == 2 ? v : 0.0);
|
||||
}
|
||||
|
||||
/* The exchange with prescribed neighbours. */
|
||||
const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w);
|
||||
double shift0[3] = { 0.0, 0.0, 0.0 };
|
||||
if (a.shifts) { shift0[0] = sh[3 * fidx]; shift0[1] = sh[3 * fidx + 1]; shift0[2] = sh[3 * fidx + 2]; }
|
||||
int cm[3] = { i, j, k }; cm[c] -= 1; /* cell minus */
|
||||
int cp[3] = { i, j, k }; /* cell plus */
|
||||
const double* apc = cut_ap(m, c);
|
||||
for (int d = 0; d < 3; ++d) {
|
||||
double a_d = area[d];
|
||||
int q[3];
|
||||
q[0] = i; q[1] = j; q[2] = k; q[d] += 1; int f_up1 = cut_face(g, c, q[0], q[1], q[2]);
|
||||
q[0] = i; q[1] = j; q[2] = k; q[d] -= 1; int f_dn1 = cut_face(g, c, q[0], q[1], q[2]);
|
||||
int presc_up = f_up1 >= 0 && !open[f_up1];
|
||||
int presc_dn = f_dn1 >= 0 && !open[f_dn1];
|
||||
if (!presc_up && !presc_dn) continue;
|
||||
q[0] = i; q[1] = j; q[2] = k; q[d] += 2; int f_up2 = cut_face(g, c, q[0], q[1], q[2]);
|
||||
q[0] = i; q[1] = j; q[2] = k; q[d] -= 2; int f_dn2 = cut_face(g, c, q[0], q[1], q[2]);
|
||||
double m_plus, m_minus;
|
||||
if (d == c) {
|
||||
double f_up = (f_up1 >= 0 ? apc[f_up1] : alpha) * (f_up1 >= 0 ? vals[f_up1] : u0);
|
||||
double f_dn = (f_dn1 >= 0 ? apc[f_dn1] : alpha) * (f_dn1 >= 0 ? vals[f_dn1] : u0);
|
||||
double f0 = alpha * u0;
|
||||
m_plus = 0.5 * (f0 + f_up) * a_d;
|
||||
m_minus = 0.5 * (f_dn + f0) * a_d;
|
||||
} else {
|
||||
const double* vd = load_val(f, d);
|
||||
const double* apd = cut_ap(m, d);
|
||||
int q1[3] = { cm[0], cm[1], cm[2] }; q1[d] += 1;
|
||||
int q2[3] = { cp[0], cp[1], cp[2] }; q2[d] += 1;
|
||||
int fa = cut_face(g, d, q1[0], q1[1], q1[2]);
|
||||
int fb = cut_face(g, d, q2[0], q2[1], q2[2]);
|
||||
int fc = cut_face(g, d, cm[0], cm[1], cm[2]);
|
||||
int fd = cut_face(g, d, cp[0], cp[1], cp[2]);
|
||||
double fla = (fa >= 0 ? apd[fa] : 1.0) * (fa >= 0 ? vd[fa] : 0.0);
|
||||
double flb = (fb >= 0 ? apd[fb] : 1.0) * (fb >= 0 ? vd[fb] : 0.0);
|
||||
double flc = (fc >= 0 ? apd[fc] : 1.0) * (fc >= 0 ? vd[fc] : 0.0);
|
||||
double fld = (fd >= 0 ? apd[fd] : 1.0) * (fd >= 0 ? vd[fd] : 0.0);
|
||||
m_plus = 0.5 * (fla + flb) * a_d;
|
||||
m_minus = 0.5 * (flc + fld) * a_d;
|
||||
}
|
||||
/* exchange.rs `solid_spacing`: the centroid spacing in a cross
|
||||
direction without the axis exchange, else `exchange_delta`. */
|
||||
double delta_x = h[d];
|
||||
if (a.axis) {
|
||||
double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||||
if (a_w0 != 0.0) {
|
||||
double n_d = fabs(wall[d]) / a_w0;
|
||||
if (!(n_d < 1e-12)) delta_x = fmin(distance / n_d, h[d]);
|
||||
}
|
||||
}
|
||||
int centroid_spacing = a.centroid && d != c && !a.axis;
|
||||
unsigned int kd = (2u << 28) | kc | (unsigned int)(4 * d);
|
||||
if (presc_up) {
|
||||
double un = vals[f_up1];
|
||||
double up1 = un;
|
||||
double dn1 = f_dn1 >= 0 ? vals[f_dn1] : 0.0;
|
||||
double up2 = f_up2 >= 0 ? vals[f_up2] : 0.0;
|
||||
(void)up1;
|
||||
double delta;
|
||||
if (g.scheme == SCHEME_UPWIND) delta = 0.0;
|
||||
else if (m_plus >= 0.0) delta = face_corr3(g.scheme, f_dn1 >= 0, dn1, u0, un);
|
||||
else delta = face_corr3(g.scheme, f_up2 >= 0, up2, un, u0);
|
||||
double u_face = upwind3(m_plus, u0, un) + delta;
|
||||
if (!a.conv_off) {
|
||||
double v = -a.rho * m_plus * (u_face - u0);
|
||||
double r = -v;
|
||||
load_push(out, count, a.cap, kd | 0u, (unsigned int)fidx,
|
||||
c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
|
||||
}
|
||||
double s = centroid_spacing ? fmin(fmax(h[d] - 1.0 * shift0[d], 0.25 * h[d]), 2.0 * h[d]) : delta_x;
|
||||
double v = a.mu * app[d] * a_d * (un - u0) / s;
|
||||
double r = -v;
|
||||
load_push(out, count, a.cap, kd | 1u, (unsigned int)fidx,
|
||||
c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
|
||||
}
|
||||
if (presc_dn) {
|
||||
double ud = vals[f_dn1];
|
||||
double dn2 = f_dn2 >= 0 ? vals[f_dn2] : 0.0;
|
||||
double up1 = f_up1 >= 0 ? vals[f_up1] : 0.0;
|
||||
double delta;
|
||||
if (g.scheme == SCHEME_UPWIND) delta = 0.0;
|
||||
else if (m_minus >= 0.0) delta = face_corr3(g.scheme, f_dn2 >= 0, dn2, ud, u0);
|
||||
else delta = face_corr3(g.scheme, f_up1 >= 0, up1, u0, ud);
|
||||
double u_face = upwind3(m_minus, ud, u0) + delta;
|
||||
if (!a.conv_off) {
|
||||
double v = a.rho * m_minus * (u_face - u0);
|
||||
double r = -v;
|
||||
load_push(out, count, a.cap, kd | 2u, (unsigned int)fidx,
|
||||
c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
|
||||
}
|
||||
double s = centroid_spacing ? fmin(fmax(h[d] - (-1.0) * shift0[d], 0.25 * h[d]), 2.0 * h[d]) : delta_x;
|
||||
double v = a.mu * apm[d] * a_d * (ud - u0) / s;
|
||||
double r = -v;
|
||||
load_push(out, count, a.cap, kd | 3u, (unsigned int)fidx,
|
||||
c == 0 ? r : 0.0, c == 1 ? r : 0.0, c == 2 ? r : 0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -7,10 +7,12 @@
|
||||
|
||||
mod cut;
|
||||
mod geom;
|
||||
mod loads;
|
||||
mod mask;
|
||||
mod poisson_setup;
|
||||
mod snapshot;
|
||||
|
||||
pub use loads::{loads_device_enabled, DeviceWallLoads, WallLoadRecord, WallLoadsCheck};
|
||||
pub use snapshot::DeviceSnapshot;
|
||||
|
||||
use super::{Side, Solver, StepResult};
|
||||
|
||||
@@ -55,7 +55,7 @@ fn cut_kernels() -> &'static CutKernels {
|
||||
/// `struct E3Cut` in e3_cut.cu: 27 device pointers.
|
||||
#[repr(C)]
|
||||
#[derive(Clone, Copy)]
|
||||
struct E3CutPtrs {
|
||||
pub(super) struct E3CutPtrs {
|
||||
ptrs: [u64; 27],
|
||||
}
|
||||
unsafe impl DeviceRepr for E3CutPtrs {}
|
||||
@@ -751,7 +751,7 @@ impl DeviceCut {
|
||||
}
|
||||
}
|
||||
|
||||
fn ptrs(&self) -> E3CutPtrs {
|
||||
pub(super) fn ptrs(&self) -> E3CutPtrs {
|
||||
let rt = runtime();
|
||||
let s = &rt.stream;
|
||||
let pf = |x: &CudaSlice<f64>| x.device_ptr(s).0;
|
||||
|
||||
@@ -0,0 +1,627 @@
|
||||
//! R8-f: the operator-route wall load on the device (`RTX_E3_LOADS_DEVICE=1`,
|
||||
//! default off). The host route — `Mask::cut_wall_force` (the cells' `p W`
|
||||
//! and the fluid faces' implicit wall shear, `cutwall.rs`) plus
|
||||
//! `Mask::cut_wall_exchange_force` (`exchange.rs`) — loops over every cell
|
||||
//! and face of the grid on the host, after the whole field comes down; in the
|
||||
//! 3D coupled step (1.45 M cells) that was 75 % of the step.
|
||||
//!
|
||||
//! Here the kernels (`e3_loads.cu`) evaluate every summand on the device from
|
||||
//! the predictor set of the device mask (the instantaneous apertures, the
|
||||
//! fluid-kind open flags, the distances, the surface velocities and the
|
||||
//! centroid shifts of the mask at the solver's time — the tables the host
|
||||
//! route reads) and the device fields, in the host's expression order with
|
||||
//! FMA contraction off, so each summand is the host's to the bit. They do
|
||||
//! not sum: each summand is appended as a record to one buffer, which comes
|
||||
//! down once per call; the host sorts the records into its own loop order
|
||||
//! (cells ascending; shear by component then face; exchange by component,
|
||||
//! face, direction, side, convective before diffusive) and sums them
|
||||
//! sequentially — the totals are therefore the host route's TO THE BIT, not
|
||||
//! merely within a tolerance (a device reduction would reorder the sums).
|
||||
//!
|
||||
//! What is not recorded: summands that are exactly zero because their cell
|
||||
//! has no wall (`W = 0`: the host adds `p · 0 = ±0` for every fluid cell).
|
||||
//! Adding `±0` to a sequential sum that starts at `+0` never changes its bits
|
||||
//! (a round-to-nearest sum of non-zero terms is never `−0`, and `x ± 0 = x`
|
||||
//! otherwise), so the totals and every consumer's sums are unchanged; the
|
||||
//! record lists replayed into the sinks drop exact zeros too (R8-c's own
|
||||
//! `cut_wall_loads` already drops all-zero loads).
|
||||
//!
|
||||
//! Supported: the cut-cell wall with the device's own closures (order 1 or
|
||||
//! 2, oblique distance, fine floor, axis exchange, centroid diffusion,
|
||||
//! exchange convection on/off). Not supported (`None`, the caller keeps the
|
||||
//! host route; logged once): the diagnostic load window, the gradient
|
||||
//! weights, the advancing-wall closure, the centroid foot, the S2-7 host
|
||||
//! prototypes, or a device mask not in its predictor phase.
|
||||
|
||||
use super::DeviceStep;
|
||||
use super::cut::Phase;
|
||||
use crate::solvers::incompressible::embedded3::Grid;
|
||||
use crate::solvers::incompressible::embedded3::body::Body;
|
||||
use crate::solvers::incompressible::embedded3::exchange::{in_load_window, to_load_sink};
|
||||
use crate::solvers::incompressible::embedded3::interface::{LoadKind, WallLoad};
|
||||
use crate::solvers::incompressible::embedded3::poisson::device::{cfg, load_module, runtime};
|
||||
use cudarc::driver::{
|
||||
CudaFunction, CudaModule, CudaSlice, DeviceRepr, PushKernelArg, ValidAsZeroBits,
|
||||
};
|
||||
use std::sync::{Arc, OnceLock};
|
||||
|
||||
const LOAD_KERNELS: &str = concat!(
|
||||
include_str!("../../../../../kernels/cuda/e3_step.cu"),
|
||||
include_str!("../../../../../kernels/cuda/e3_cut.cu"),
|
||||
include_str!("../../../../../kernels/cuda/e3_loads.cu")
|
||||
);
|
||||
|
||||
struct LoadKernels {
|
||||
_module: Arc<CudaModule>,
|
||||
cells: CudaFunction,
|
||||
faces: CudaFunction,
|
||||
}
|
||||
|
||||
static LOAD_ONCE: OnceLock<LoadKernels> = OnceLock::new();
|
||||
|
||||
/// The record buffer's capacity (grown when a pass overflows it).
|
||||
static LOADS_CAP: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0);
|
||||
|
||||
fn kernels() -> &'static LoadKernels {
|
||||
LOAD_ONCE.get_or_init(|| {
|
||||
// FMA contraction off: the host's roundings.
|
||||
let module = load_module(LOAD_KERNELS, "e3_loads.cu", true);
|
||||
let f = |name: &str| module.load_function(name).expect(name);
|
||||
LoadKernels {
|
||||
cells: f("e3_loads_cells"),
|
||||
faces: f("e3_loads_faces"),
|
||||
_module: module,
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
/// `RTX_E3_LOADS_DEVICE=1` (read once): callers take the device route.
|
||||
pub fn loads_device_enabled() -> bool {
|
||||
static ON: OnceLock<bool> = OnceLock::new();
|
||||
*ON.get_or_init(|| std::env::var("RTX_E3_LOADS_DEVICE").is_ok_and(|v| v == "1"))
|
||||
}
|
||||
|
||||
fn log_fallback(reason: &str) {
|
||||
static LOGGED: std::sync::Once = std::sync::Once::new();
|
||||
LOGGED.call_once(|| {
|
||||
eprintln!(" R8-f device loads: host route ({reason})");
|
||||
});
|
||||
}
|
||||
|
||||
/// `struct LoadRec` in e3_loads.cu.
|
||||
#[repr(C)]
|
||||
#[derive(Clone, Copy, Debug, Default)]
|
||||
struct LoadRec {
|
||||
key: u32,
|
||||
idx: u32,
|
||||
v: [f64; 3],
|
||||
}
|
||||
unsafe impl DeviceRepr for LoadRec {}
|
||||
unsafe impl ValidAsZeroBits for LoadRec {}
|
||||
|
||||
/// `struct LoadArgs` in e3_loads.cu.
|
||||
#[repr(C)]
|
||||
#[derive(Clone, Copy)]
|
||||
struct LoadArgs {
|
||||
mu: f64,
|
||||
rho: f64,
|
||||
k0: i32,
|
||||
k1: i32,
|
||||
shifts: i32,
|
||||
centroid: i32,
|
||||
axis: i32,
|
||||
conv_off: i32,
|
||||
order: i32,
|
||||
cap: u32,
|
||||
}
|
||||
unsafe impl DeviceRepr for LoadArgs {}
|
||||
unsafe impl ValidAsZeroBits for LoadArgs {}
|
||||
|
||||
/// One summand of the operator route (a force ON the body), in the host's
|
||||
/// loop order.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct WallLoadRecord {
|
||||
/// R8-a's part numbering: 0 the cell's `p W`, 1 the face's wall shear,
|
||||
/// 2 the diffusive and 3 the convective exchange.
|
||||
pub part: usize,
|
||||
/// The face's component (0 for a cell).
|
||||
pub c: usize,
|
||||
/// The cell or face index (in its own lattice).
|
||||
pub index: usize,
|
||||
/// The cell centre or the face position.
|
||||
pub pos: [f64; 3],
|
||||
/// The force; a face summand has only component `c`.
|
||||
pub value: [f64; 3],
|
||||
}
|
||||
|
||||
/// The device route's parts (the host's sums, to the bit) and its records.
|
||||
#[derive(Clone, Debug, Default)]
|
||||
pub struct DeviceWallLoads {
|
||||
pub pressure: [f64; 3],
|
||||
pub shear: [f64; 3],
|
||||
pub diffusive: [f64; 3],
|
||||
pub convective: [f64; 3],
|
||||
pub records: Vec<WallLoadRecord>,
|
||||
/// `Some((k0, k1))` for a plane-restricted route.
|
||||
pub planes: Option<(usize, usize)>,
|
||||
}
|
||||
|
||||
impl DeviceWallLoads {
|
||||
/// The route's total as `Mask::cut_wall_force` forms it.
|
||||
#[must_use]
|
||||
pub fn total(&self) -> [f64; 3] {
|
||||
let (p, s) = (self.pressure, self.shear);
|
||||
let x = [
|
||||
self.diffusive[0] + self.convective[0],
|
||||
self.diffusive[1] + self.convective[1],
|
||||
self.diffusive[2] + self.convective[2],
|
||||
];
|
||||
[p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]
|
||||
}
|
||||
|
||||
/// The summands as the R8-a load sink receives them (position,
|
||||
/// component, part, value), exact zeros dropped.
|
||||
#[must_use]
|
||||
pub fn sink_items(&self) -> Vec<([f64; 3], usize, usize, f64)> {
|
||||
let mut out = Vec::with_capacity(self.records.len() + self.records.len() / 4);
|
||||
for r in &self.records {
|
||||
if r.part == 0 {
|
||||
for c in 0..3 {
|
||||
if r.value[c] != 0.0 {
|
||||
out.push((r.pos, c, 0, r.value[c]));
|
||||
}
|
||||
}
|
||||
} else if r.value[r.c] != 0.0 {
|
||||
out.push((r.pos, r.c, r.part, r.value[r.c]));
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// Hand the summands to the process-wide R8-a load sink (`exchange::
|
||||
/// set_load_sink`), in the host route's order, exact zeros dropped.
|
||||
pub fn replay_load_sink(&self) {
|
||||
for (pos, c, part, v) in self.sink_items() {
|
||||
to_load_sink(pos, c, part, v);
|
||||
}
|
||||
}
|
||||
|
||||
/// R8-c's summands (`Mask::cut_wall_loads`: kind, position, foot, force;
|
||||
/// all-zero loads dropped), the foot by the body's own φ as there.
|
||||
#[must_use]
|
||||
pub fn wall_loads(&self, body: &Body, grid: Grid, t: f64) -> Vec<WallLoad> {
|
||||
let raw: Vec<(LoadKind, [f64; 3], [f64; 3])> = self
|
||||
.records
|
||||
.iter()
|
||||
.filter(|r| r.value != [0.0; 3])
|
||||
.map(|r| {
|
||||
let kind = match r.part {
|
||||
0 => LoadKind::Pressure,
|
||||
1 => LoadKind::Shear,
|
||||
2 => LoadKind::ExchangeDiffusive,
|
||||
_ => LoadKind::ExchangeConvective,
|
||||
};
|
||||
(kind, r.pos, r.value)
|
||||
})
|
||||
.collect();
|
||||
// interface.rs `cut_wall_loads`: two projections onto φ = 0.
|
||||
let eps = 1e-6 * grid.dx.min(grid.dy).min(grid.dz);
|
||||
let foot = |x: [f64; 3]| {
|
||||
let mut q = x;
|
||||
for _ in 0..2 {
|
||||
let s = body.phi(q[0], q[1], q[2], t);
|
||||
let n = body.normal(q[0], q[1], q[2], t, eps);
|
||||
q = [q[0] - s * n.0, q[1] - s * n.1, q[2] - s * n.2];
|
||||
}
|
||||
q
|
||||
};
|
||||
use rayon::prelude::*;
|
||||
raw.par_iter()
|
||||
.map(|&(kind, x, f)| WallLoad {
|
||||
kind,
|
||||
x,
|
||||
foot: foot(x),
|
||||
f,
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// The G1 instrument's verdict for one call ([`DeviceStep::check_wall_loads`]).
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct WallLoadsCheck {
|
||||
/// Host and device totals (per span when plane-restricted).
|
||||
pub host: [f64; 3],
|
||||
pub device: [f64; 3],
|
||||
/// The totals bit-identical.
|
||||
pub totals_identical: bool,
|
||||
/// The host sink's summands (incl. exact zeros), the non-zero ones, and
|
||||
/// the device's replayed ones.
|
||||
pub host_items: usize,
|
||||
pub host_nonzero: usize,
|
||||
pub device_items: usize,
|
||||
/// Non-zero host summands and device summands identical item by item
|
||||
/// (position, component, part, value bits).
|
||||
pub items_identical: bool,
|
||||
/// The sums of the two item lists (component-wise, sequential).
|
||||
pub host_item_sum: [f64; 3],
|
||||
pub device_item_sum: [f64; 3],
|
||||
/// The records summed in the device's (atomic, unsorted) order: the
|
||||
/// relative difference a device-order reduction would have made.
|
||||
pub device_order_rel: f64,
|
||||
/// R8-c's summands (whole body, no planes): identical to the host's
|
||||
/// `cut_wall_loads` (None when not checked).
|
||||
pub r8c_identical: Option<bool>,
|
||||
/// Wall time of the host route (on the downloaded field) and of the
|
||||
/// device route (kernels, download of the records, sort, sums), ms.
|
||||
pub host_ms: f64,
|
||||
pub device_ms: f64,
|
||||
}
|
||||
|
||||
impl DeviceStep {
|
||||
/// Why the device route cannot serve this solver (None = it can).
|
||||
fn loads_unsupported(&self) -> Option<&'static str> {
|
||||
let mask = self.solver.mask()?;
|
||||
if mask.cut.is_none() {
|
||||
return Some("no cut geometry");
|
||||
}
|
||||
let Some(dc) = self.cut.as_ref() else {
|
||||
return Some("no device mask");
|
||||
};
|
||||
if dc.phase != Phase::Predictor {
|
||||
return Some("the device mask is not in its predictor phase");
|
||||
}
|
||||
if !in_load_window(f64::NAN) {
|
||||
return Some("the diagnostic load window is set");
|
||||
}
|
||||
if mask.grad_weights.is_some() {
|
||||
return Some("the gradient weights (S2-5 host prototype)");
|
||||
}
|
||||
if mask.wall_advancing {
|
||||
return Some("the advancing-wall closure");
|
||||
}
|
||||
if mask.wall_foot_centroid {
|
||||
return Some("the centroid foot");
|
||||
}
|
||||
if mask.cv_sides_exact || mask.wall_order2_centroid {
|
||||
return Some("an S2-7 host prototype");
|
||||
}
|
||||
if mask.periodic_z
|
||||
!= (self.solver.params.boundaries.z0
|
||||
== crate::solvers::incompressible::embedded3::step::Side::Periodic)
|
||||
{
|
||||
return Some("the mask's z periodicity differs from the device's");
|
||||
}
|
||||
None
|
||||
}
|
||||
|
||||
/// The operator route's parts and records on the device over the planes
|
||||
/// `planes` (all when `None`), with `mu` the dynamic viscosity (the host
|
||||
/// route's argument). `None` when the device cannot serve it (the caller
|
||||
/// keeps the host route; the reason is logged once).
|
||||
pub fn cut_wall_loads_device(
|
||||
&self,
|
||||
mu: f64,
|
||||
planes: Option<(usize, usize)>,
|
||||
) -> Option<DeviceWallLoads> {
|
||||
if let Some(reason) = self.loads_unsupported() {
|
||||
log_fallback(reason);
|
||||
return None;
|
||||
}
|
||||
self.solver.mask()?;
|
||||
let (recs, _) = self.load_records(mu, planes);
|
||||
let mut records = recs;
|
||||
records.sort_unstable_by_key(|r| (r.key >> 24, r.idx, r.key & 0x00ff_ffff));
|
||||
Some(self.assemble(&records, planes))
|
||||
}
|
||||
|
||||
/// The route's total like `Mask::cut_wall_force` (whole body, `None`
|
||||
/// planes) or `Mask::cut_wall_force_per_span` (the planes, divided by
|
||||
/// their thickness), with the summands handed to the R8-a load sink as
|
||||
/// the host route hands them (exact zeros dropped). `None`: the caller
|
||||
/// keeps the host route.
|
||||
pub fn cut_wall_force_device(
|
||||
&self,
|
||||
mu: f64,
|
||||
planes: Option<(usize, usize)>,
|
||||
) -> Option<[f64; 3]> {
|
||||
let loads = self.cut_wall_loads_device(mu, planes)?;
|
||||
loads.replay_load_sink();
|
||||
let t = loads.total();
|
||||
Some(match planes {
|
||||
None => t,
|
||||
Some((k0, k1)) => {
|
||||
let lz = (k1 - k0) as f64 * self.grid.dz;
|
||||
[t[0] / lz, t[1] / lz, t[2] / lz]
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
/// G1 instrument (`RTX_E3_LOADS_CHECK=1` in the drivers): the device
|
||||
/// route against the host route on the downloaded field — totals, the
|
||||
/// R8-a sink's items and (whole body) R8-c's summands, bit for bit.
|
||||
/// Installs its own capture sink: no load sink may be installed.
|
||||
pub fn check_wall_loads(
|
||||
&self,
|
||||
mu: f64,
|
||||
planes: Option<(usize, usize)>,
|
||||
) -> Option<WallLoadsCheck> {
|
||||
use crate::solvers::incompressible::embedded3::exchange::set_load_sink;
|
||||
use crate::solvers::incompressible::embedded3::field::Field;
|
||||
type Items = Vec<([f64; 3], usize, usize, f64)>;
|
||||
if self.loads_unsupported().is_some() {
|
||||
return None;
|
||||
}
|
||||
let mut field = Field::new(self.grid);
|
||||
self.download(&mut field);
|
||||
let t = self.solver.time();
|
||||
let capture = || {
|
||||
let store = Arc::new(std::sync::Mutex::new(Items::new()));
|
||||
let s = store.clone();
|
||||
let prev = set_load_sink(Some(Box::new(move |pos, c, part, v| {
|
||||
s.lock().expect("capture").push((pos, c, part, v));
|
||||
})));
|
||||
assert!(
|
||||
prev.is_none(),
|
||||
"check_wall_loads: a load sink was already installed"
|
||||
);
|
||||
store
|
||||
};
|
||||
let take = |store: Arc<std::sync::Mutex<Items>>| -> Items {
|
||||
set_load_sink(None);
|
||||
std::mem::take(&mut *store.lock().expect("capture"))
|
||||
};
|
||||
// Host.
|
||||
let (mask, body) = (self.solver.mask()?, self.solver.body()?);
|
||||
let lap = std::time::Instant::now();
|
||||
let store = capture();
|
||||
let host = match planes {
|
||||
None => mask.cut_wall_force(body, &field, mu, t),
|
||||
Some(p) => mask.cut_wall_force_per_span(body, &field, mu, t, p),
|
||||
};
|
||||
let host_items = take(store);
|
||||
let host_ms = lap.elapsed().as_secs_f64() * 1e3;
|
||||
let host = host?;
|
||||
// The device route's own cost (what a caller pays).
|
||||
let lap = std::time::Instant::now();
|
||||
let timed = self.cut_wall_loads_device(mu, planes);
|
||||
let device_ms = lap.elapsed().as_secs_f64() * 1e3;
|
||||
drop(timed);
|
||||
// Device (records in both orders).
|
||||
let (raw, _) = self.load_records(mu, planes);
|
||||
let mut device_order = [0.0f64; 3];
|
||||
for r in &raw {
|
||||
let ll = self.assemble(std::slice::from_ref(r), planes);
|
||||
let v = ll.total();
|
||||
for c in 0..3 {
|
||||
device_order[c] += v[c];
|
||||
}
|
||||
}
|
||||
let mut sorted = raw;
|
||||
sorted.sort_unstable_by_key(|r| (r.key >> 24, r.idx, r.key & 0x00ff_ffff));
|
||||
let loads = self.assemble(&sorted, planes);
|
||||
let store = capture();
|
||||
loads.replay_load_sink();
|
||||
let device_items = take(store);
|
||||
let device = match planes {
|
||||
None => loads.total(),
|
||||
Some((k0, k1)) => {
|
||||
let lz = (k1 - k0) as f64 * self.grid.dz;
|
||||
let v = loads.total();
|
||||
[v[0] / lz, v[1] / lz, v[2] / lz]
|
||||
}
|
||||
};
|
||||
let bits = |v: [f64; 3]| v.map(f64::to_bits);
|
||||
let nonzero: Items = host_items.iter().copied().filter(|x| x.3 != 0.0).collect();
|
||||
let same_item = |a: &([f64; 3], usize, usize, f64), b: &([f64; 3], usize, usize, f64)| {
|
||||
bits(a.0) == bits(b.0) && a.1 == b.1 && a.2 == b.2 && a.3.to_bits() == b.3.to_bits()
|
||||
};
|
||||
let items_identical = nonzero.len() == device_items.len()
|
||||
&& nonzero
|
||||
.iter()
|
||||
.zip(&device_items)
|
||||
.all(|(a, b)| same_item(a, b));
|
||||
let item_sum = |v: &Items| {
|
||||
let mut s = [0.0f64; 3];
|
||||
for x in v {
|
||||
s[x.1] += x.3;
|
||||
}
|
||||
s
|
||||
};
|
||||
let scale = host
|
||||
.iter()
|
||||
.fold(0.0f64, |m, v| m.max(v.abs()))
|
||||
.max(f64::MIN_POSITIVE);
|
||||
let lz = planes.map_or(1.0, |(k0, k1)| (k1 - k0) as f64 * self.grid.dz);
|
||||
let device_order_rel = (0..3)
|
||||
.map(|c| (device_order[c] / lz - host[c]).abs() / scale)
|
||||
.fold(0.0f64, f64::max);
|
||||
let r8c_identical = if planes.is_none() {
|
||||
let mask = self.solver.mask()?;
|
||||
let body = self.solver.body()?;
|
||||
mask.cut_wall_loads(body, &field, mu, t).map(|(h, route)| {
|
||||
let d = loads.wall_loads(body, self.grid, t);
|
||||
bits(route) == bits(loads.total())
|
||||
&& h.len() == d.len()
|
||||
&& h.iter().zip(&d).all(|(a, b)| {
|
||||
a.kind == b.kind
|
||||
&& bits(a.x) == bits(b.x)
|
||||
&& bits(a.foot) == bits(b.foot)
|
||||
&& bits(a.f) == bits(b.f)
|
||||
})
|
||||
})
|
||||
} else {
|
||||
None
|
||||
};
|
||||
Some(WallLoadsCheck {
|
||||
host,
|
||||
device,
|
||||
totals_identical: bits(host) == bits(device),
|
||||
host_items: host_items.len(),
|
||||
host_nonzero: nonzero.len(),
|
||||
device_items: device_items.len(),
|
||||
items_identical,
|
||||
host_item_sum: item_sum(&nonzero),
|
||||
device_item_sum: item_sum(&device_items),
|
||||
device_order_rel,
|
||||
r8c_identical,
|
||||
host_ms,
|
||||
device_ms,
|
||||
})
|
||||
}
|
||||
|
||||
/// The records of one pass (unsorted: the device's append order) and
|
||||
/// the number the kernels produced.
|
||||
fn load_records(&self, mu: f64, planes: Option<(usize, usize)>) -> (Vec<LoadRec>, usize) {
|
||||
let rt = runtime();
|
||||
let k = kernels();
|
||||
let g = self.grid;
|
||||
let mask = self.solver.mask().expect("mask");
|
||||
let (k0, k1) = planes.unwrap_or((0, g.nz));
|
||||
let mut prm = self.params(1.0);
|
||||
// cut_cv reads the order's low bits, the oblique distance (16) and
|
||||
// the fine floor (256): the MASK's settings (the host route's).
|
||||
prm.wall_order = i32::from(mask.wall_order)
|
||||
+ 16 * i32::from(mask.wall_distance_oblique)
|
||||
+ 256 * i32::from(mask.distance_floor_fine);
|
||||
prm.scheme = super::scheme_code(mask.scheme);
|
||||
let mut args = LoadArgs {
|
||||
mu,
|
||||
rho: mask.density,
|
||||
k0: k0 as i32,
|
||||
k1: k1 as i32,
|
||||
shifts: i32::from(mask.face_shifts.is_some()),
|
||||
centroid: i32::from(mask.diffusion_centroid),
|
||||
axis: i32::from(mask.wall_exchange_axis),
|
||||
conv_off: i32::from(mask.exchange_convection_off),
|
||||
order: i32::from(mask.wall_order),
|
||||
cap: 0,
|
||||
};
|
||||
let ptrs = self.ptrs();
|
||||
let cptrs = self.cut.as_ref().expect("cut").ptrs();
|
||||
let counts = [
|
||||
(g.nx + 1) * g.ny * g.nz,
|
||||
g.nx * (g.ny + 1) * g.nz,
|
||||
g.nx * g.ny * (g.nz + 1),
|
||||
];
|
||||
let cells = g.nx * g.ny * (k1 - k0);
|
||||
let mut cap = LOADS_CAP
|
||||
.load(std::sync::atomic::Ordering::Relaxed)
|
||||
.max(1 << 16);
|
||||
loop {
|
||||
args.cap = cap as u32;
|
||||
let mut out = rt
|
||||
.stream
|
||||
.alloc_zeros::<LoadRec>(cap)
|
||||
.expect("alloc records");
|
||||
let mut count = rt.stream.alloc_zeros::<u32>(1).expect("alloc count");
|
||||
unsafe {
|
||||
rt.stream
|
||||
.launch_builder(&k.cells)
|
||||
.arg(&prm)
|
||||
.arg(&ptrs)
|
||||
.arg(&cptrs)
|
||||
.arg(&args)
|
||||
.arg(&mut out)
|
||||
.arg(&mut count)
|
||||
.launch(cfg(cells))
|
||||
.expect("e3_loads_cells");
|
||||
for c in 0..3i32 {
|
||||
rt.stream
|
||||
.launch_builder(&k.faces)
|
||||
.arg(&prm)
|
||||
.arg(&ptrs)
|
||||
.arg(&cptrs)
|
||||
.arg(&args)
|
||||
.arg(&c)
|
||||
.arg(&mut out)
|
||||
.arg(&mut count)
|
||||
.launch(cfg(counts[c as usize]))
|
||||
.expect("e3_loads_faces");
|
||||
}
|
||||
}
|
||||
let n: Vec<u32> = rt.stream.memcpy_dtov(&count).expect("count");
|
||||
let n = n[0] as usize;
|
||||
if n <= cap {
|
||||
let recs: Vec<LoadRec> = if n == 0 {
|
||||
Vec::new()
|
||||
} else {
|
||||
rt.stream.memcpy_dtov(&out.slice(0..n)).expect("records")
|
||||
};
|
||||
LOADS_CAP.store(cap, std::sync::atomic::Ordering::Relaxed);
|
||||
return (recs, n);
|
||||
}
|
||||
cap = n + n / 4;
|
||||
}
|
||||
}
|
||||
|
||||
/// The parts summed in the order of `records` (sequentially, as the host
|
||||
/// loops do) and the records with their positions.
|
||||
fn assemble(&self, records: &[LoadRec], planes: Option<(usize, usize)>) -> DeviceWallLoads {
|
||||
let g = self.grid;
|
||||
let h = [g.dx, g.dy, g.dz];
|
||||
let (nx, ny) = (g.nx, g.ny);
|
||||
let mut out = DeviceWallLoads {
|
||||
planes,
|
||||
records: Vec::with_capacity(records.len()),
|
||||
..DeviceWallLoads::default()
|
||||
};
|
||||
for r in records {
|
||||
let kind = r.key >> 28;
|
||||
let c = ((r.key >> 24) & 0xf) as usize;
|
||||
let slot = (r.key & 0x00ff_ffff) as usize;
|
||||
let idx = r.idx as usize;
|
||||
let (part, pos) = if kind == 0 {
|
||||
let (k, j, i) = (idx / (nx * ny), (idx / nx) % ny, idx % nx);
|
||||
(
|
||||
0,
|
||||
[
|
||||
(i as f64 + 0.5) * g.dx,
|
||||
(j as f64 + 0.5) * g.dy,
|
||||
(k as f64 + 0.5) * g.dz,
|
||||
],
|
||||
)
|
||||
} else {
|
||||
let (ni, nj) = (nx + usize::from(c == 0), ny + usize::from(c == 1));
|
||||
let p = [idx % ni, (idx / ni) % nj, idx / (ni * nj)];
|
||||
let mut x = [0.0; 3];
|
||||
for d in 0..3 {
|
||||
// cutwall.rs `Lattice::face_position`.
|
||||
let off = if d == c { 0.0 } else { 0.5 };
|
||||
x[d] = (p[d] as f64 + off) * h[d];
|
||||
}
|
||||
let part = if kind == 1 {
|
||||
1
|
||||
} else if slot & 1 == 0 {
|
||||
3
|
||||
} else {
|
||||
2
|
||||
};
|
||||
(part, x)
|
||||
};
|
||||
match part {
|
||||
0 => {
|
||||
for cc in 0..3 {
|
||||
out.pressure[cc] += r.v[cc];
|
||||
}
|
||||
}
|
||||
1 => out.shear[c] += r.v[c],
|
||||
2 => out.diffusive[c] += r.v[c],
|
||||
_ => out.convective[c] += r.v[c],
|
||||
}
|
||||
out.records.push(WallLoadRecord {
|
||||
part,
|
||||
c,
|
||||
index: idx,
|
||||
pos,
|
||||
value: r.v,
|
||||
});
|
||||
}
|
||||
// cutwall.rs adds the gradient-weight force (zero: the weights are
|
||||
// off on this route) to the pressure part.
|
||||
for cc in 0..3 {
|
||||
out.pressure[cc] += 0.0;
|
||||
}
|
||||
out
|
||||
}
|
||||
}
|
||||
@@ -650,23 +650,80 @@ fn flag_wake_on_the_device() {
|
||||
let (mut drag_sum, mut lift_min, mut lift_max, mut samples) =
|
||||
(0.0, f64::INFINITY, f64::NEG_INFINITY, 0usize);
|
||||
let mut worst_residual = 0.0_f64;
|
||||
// R8-f: `RTX_E3_LOADS_DEVICE=1` takes the operator route on the device
|
||||
// (the per-span and whole-body totals below); `RTX_E3_LOADS_CHECK=1`
|
||||
// checks the device route against the host on EVERY step (both the
|
||||
// slab's planes and the whole body), bit for bit.
|
||||
let loads_device =
|
||||
rtx_cfd::solvers::incompressible::embedded3::step::device::loads_device_enabled();
|
||||
let loads_check = std::env::var("RTX_E3_LOADS_CHECK").is_ok_and(|v| v == "1");
|
||||
let mut loads_bad = 0usize;
|
||||
for step in 0..steps {
|
||||
let r = device.advance(dt);
|
||||
worst_residual = worst_residual.max(r.final_residual);
|
||||
assert!(r.final_residual.is_finite(), "death at step {step}");
|
||||
let t = device.solver.time();
|
||||
if loads_check {
|
||||
for (name, planes) in [("slab", Some(slab)), ("whole", None)] {
|
||||
let ck = device
|
||||
.check_wall_loads(RHO * NU, planes)
|
||||
.expect("R8-f: the device route is unsupported here");
|
||||
let ok =
|
||||
ck.totals_identical && ck.items_identical && ck.r8c_identical.unwrap_or(true);
|
||||
loads_bad += usize::from(!ok);
|
||||
println!(
|
||||
" loads check step {step} t {t:.5} {name}: {} — totals host ({:+.12e}, {:+.12e}, {:+.12e}) device ({:+.12e}, {:+.12e}, {:+.12e}); items host {} ({} non-zero) device {} {}; item sums x {:+.12e}/{:+.12e} y {:+.12e}/{:+.12e}; R8-c {}; device-order rel {:.2e}; host {:.1} ms device {:.1} ms",
|
||||
if ok { "IDENTICAL" } else { "DIFFERS" },
|
||||
ck.host[0],
|
||||
ck.host[1],
|
||||
ck.host[2],
|
||||
ck.device[0],
|
||||
ck.device[1],
|
||||
ck.device[2],
|
||||
ck.host_items,
|
||||
ck.host_nonzero,
|
||||
ck.device_items,
|
||||
if ck.items_identical {
|
||||
"identical"
|
||||
} else {
|
||||
"DIFFER"
|
||||
},
|
||||
ck.host_item_sum[0],
|
||||
ck.device_item_sum[0],
|
||||
ck.host_item_sum[1],
|
||||
ck.device_item_sum[1],
|
||||
match ck.r8c_identical {
|
||||
Some(true) => "identical",
|
||||
Some(false) => "DIFFERS",
|
||||
None => "n/a",
|
||||
},
|
||||
ck.device_order_rel,
|
||||
ck.host_ms,
|
||||
ck.device_ms
|
||||
);
|
||||
}
|
||||
}
|
||||
let sample = (step + 1) % 10 == 0 || step + 1 == steps;
|
||||
let phase_due = vtk_dir.is_some()
|
||||
&& t >= last_period_start + next_phase as f64 * period / phases as f64
|
||||
&& next_phase < total_phases;
|
||||
if sample || phase_due {
|
||||
let dev = loads_device.then(|| {
|
||||
(
|
||||
device.cut_wall_force_device(RHO * NU, Some(slab)),
|
||||
device.cut_wall_force_device(RHO * NU, None),
|
||||
)
|
||||
});
|
||||
device.download(&mut field);
|
||||
let solver = &device.solver;
|
||||
let mask = solver.mask().expect("mask");
|
||||
let body = solver.body().expect("body");
|
||||
let fs = mask
|
||||
.cut_wall_force_per_span(body, &field, RHO * NU, t, slab)
|
||||
.expect("wall");
|
||||
let fs = match dev {
|
||||
Some((Some(v), _)) => v,
|
||||
_ => mask
|
||||
.cut_wall_force_per_span(body, &field, RHO * NU, t, slab)
|
||||
.expect("wall"),
|
||||
};
|
||||
// The reconstructed wall route on the same slab, per span.
|
||||
let fr = mask
|
||||
.cut_wall_force_reconstructed(body, &field, RHO * NU, t, Some(slab))
|
||||
@@ -675,9 +732,12 @@ fn flag_wake_on_the_device() {
|
||||
[v[0] / lz, v[1] / lz, v[2] / lz]
|
||||
})
|
||||
.expect("reconstructed");
|
||||
let ft = mask
|
||||
.cut_wall_force(body, &field, RHO * NU, t)
|
||||
.expect("wall");
|
||||
let ft = match dev {
|
||||
Some((_, Some(v))) => v,
|
||||
_ => mask
|
||||
.cut_wall_force(body, &field, RHO * NU, t)
|
||||
.expect("wall"),
|
||||
};
|
||||
if sample {
|
||||
samples_seen += 1;
|
||||
}
|
||||
@@ -910,4 +970,8 @@ fn flag_wake_on_the_device() {
|
||||
if let Some(t) = device.timers() {
|
||||
println!(" timers: {t:?}");
|
||||
}
|
||||
if loads_check {
|
||||
println!(" R8-f LOADS CHECK: {loads_bad} calls differ");
|
||||
assert_eq!(loads_bad, 0, "R8-f: the device route differs from the host");
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user