#![allow(dead_code)]
use crate::math::{Fix128, Vec3Fix};
#[cfg(not(feature = "std"))]
use alloc::collections::BTreeMap;
#[cfg(not(feature = "std"))]
use alloc::vec;
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;
#[cfg(feature = "std")]
use std::collections::BTreeMap;
#[derive(Clone, Copy, Debug, PartialEq, Eq, Default)]
#[non_exhaustive]
pub enum FaceBc {
#[default]
Fluid,
Wall {
velocity: Vec3Fix,
},
SlipWall,
Inflow {
normal_velocity: Fix128,
},
Outflow,
}
impl FaceBc {
#[must_use]
pub fn is_wall(self) -> bool {
matches!(self, Self::Wall { .. } | Self::SlipWall)
}
#[must_use]
pub fn is_inflow(self) -> bool {
matches!(self, Self::Inflow { .. })
}
#[must_use]
pub fn blocks_pressure(self) -> bool {
self.is_wall() || self.is_inflow()
}
#[must_use]
pub fn no_slip_velocity(self) -> Option<Vec3Fix> {
match self {
Self::Wall { velocity } => Some(velocity),
_ => None,
}
}
}
fn read_face_bc(map: &BTreeMap<usize, FaceBc>, solid: &[bool], ix: usize) -> FaceBc {
if let Some(bc) = map.get(&ix) {
return *bc;
}
if solid[ix] {
FaceBc::Wall {
velocity: Vec3Fix::ZERO,
}
} else {
FaceBc::Fluid
}
}
fn write_face_bc(map: &mut BTreeMap<usize, FaceBc>, ix: usize, bc: FaceBc) {
let default_wall = FaceBc::Wall {
velocity: Vec3Fix::ZERO,
};
if bc == FaceBc::Fluid || bc == default_wall {
map.remove(&ix);
} else {
map.insert(ix, bc);
}
}
fn wall_or_fluid(solid: bool) -> FaceBc {
if solid {
FaceBc::Wall {
velocity: Vec3Fix::ZERO,
}
} else {
FaceBc::Fluid
}
}
fn wall_between(a: Option<FaceBc>, b: Option<FaceBc>) -> Option<Vec3Fix> {
let side = |bc: Option<FaceBc>| match bc {
None => Some(None),
Some(bc) => bc.no_slip_velocity().map(Some),
};
match (side(a)?, side(b)?) {
(Some(p), Some(q)) => Some(Vec3Fix::new(
(p.x + q.x).half(),
(p.y + q.y).half(),
(p.z + q.z).half(),
)),
(Some(p), None) | (None, Some(p)) => Some(p),
(None, None) => None,
}
}
#[derive(Clone, Debug)]
pub struct MacGrid {
pub nx: usize,
pub ny: usize,
pub nz: usize,
pub dx: Fix128,
pub u: Vec<Fix128>,
pub v: Vec<Fix128>,
pub w: Vec<Fix128>,
pub pressure: Vec<Fix128>,
pub u_solid: Vec<bool>,
pub v_solid: Vec<bool>,
pub w_solid: Vec<bool>,
u_face_bc: BTreeMap<usize, FaceBc>,
v_face_bc: BTreeMap<usize, FaceBc>,
w_face_bc: BTreeMap<usize, FaceBc>,
}
impl MacGrid {
#[must_use]
pub fn new(nx: usize, ny: usize, nz: usize, dx: Fix128) -> Self {
Self {
nx,
ny,
nz,
dx,
u: vec![Fix128::ZERO; (nx + 1) * ny * nz],
v: vec![Fix128::ZERO; nx * (ny + 1) * nz],
w: vec![Fix128::ZERO; nx * ny * (nz + 1)],
pressure: vec![Fix128::ZERO; nx * ny * nz],
u_solid: vec![false; (nx + 1) * ny * nz],
v_solid: vec![false; nx * (ny + 1) * nz],
w_solid: vec![false; nx * ny * (nz + 1)],
u_face_bc: BTreeMap::new(),
v_face_bc: BTreeMap::new(),
w_face_bc: BTreeMap::new(),
}
}
#[inline]
#[must_use]
pub(crate) fn idx_u(&self, i: usize, j: usize, k: usize) -> usize {
i + (self.nx + 1) * (j + self.ny * k)
}
#[inline]
#[must_use]
pub(crate) fn idx_v(&self, i: usize, j: usize, k: usize) -> usize {
i + self.nx * (j + (self.ny + 1) * k)
}
#[inline]
#[must_use]
pub(crate) fn idx_w(&self, i: usize, j: usize, k: usize) -> usize {
i + self.nx * (j + self.ny * k)
}
#[inline]
#[must_use]
fn idx_c(&self, i: usize, j: usize, k: usize) -> usize {
i + self.nx * (j + self.ny * k)
}
#[must_use]
pub fn u(&self, i: usize, j: usize, k: usize) -> Fix128 {
if i > self.nx || j >= self.ny || k >= self.nz {
return Fix128::ZERO;
}
self.u[self.idx_u(i, j, k)]
}
#[must_use]
pub fn v(&self, i: usize, j: usize, k: usize) -> Fix128 {
if i >= self.nx || j > self.ny || k >= self.nz {
return Fix128::ZERO;
}
self.v[self.idx_v(i, j, k)]
}
#[must_use]
pub fn w(&self, i: usize, j: usize, k: usize) -> Fix128 {
if i >= self.nx || j >= self.ny || k > self.nz {
return Fix128::ZERO;
}
self.w[self.idx_w(i, j, k)]
}
#[must_use]
pub fn pressure(&self, i: usize, j: usize, k: usize) -> Fix128 {
if i >= self.nx || j >= self.ny || k >= self.nz {
return Fix128::ZERO;
}
self.pressure[self.idx_c(i, j, k)]
}
#[must_use]
pub fn cell_velocity(&self, i: usize, j: usize, k: usize) -> (Fix128, Fix128, Fix128) {
let u_c = (self.u(i, j, k) + self.u(i + 1, j, k)).half();
let v_c = (self.v(i, j, k) + self.v(i, j + 1, k)).half();
let w_c = (self.w(i, j, k) + self.w(i, j, k + 1)).half();
(u_c, v_c, w_c)
}
#[must_use]
pub fn is_u_solid(&self, i: usize, j: usize, k: usize) -> bool {
if i > self.nx || j >= self.ny || k >= self.nz {
return false;
}
self.u_solid[self.idx_u(i, j, k)]
}
#[must_use]
pub fn is_v_solid(&self, i: usize, j: usize, k: usize) -> bool {
if i >= self.nx || j > self.ny || k >= self.nz {
return false;
}
self.v_solid[self.idx_v(i, j, k)]
}
#[must_use]
pub fn is_w_solid(&self, i: usize, j: usize, k: usize) -> bool {
if i >= self.nx || j >= self.ny || k > self.nz {
return false;
}
self.w_solid[self.idx_w(i, j, k)]
}
pub fn set_u_solid(&mut self, i: usize, j: usize, k: usize, solid: bool) {
self.set_u_bc(i, j, k, wall_or_fluid(solid));
}
pub fn set_v_solid(&mut self, i: usize, j: usize, k: usize, solid: bool) {
self.set_v_bc(i, j, k, wall_or_fluid(solid));
}
pub fn set_w_solid(&mut self, i: usize, j: usize, k: usize, solid: bool) {
self.set_w_bc(i, j, k, wall_or_fluid(solid));
}
#[must_use]
pub fn u_bc(&self, i: usize, j: usize, k: usize) -> FaceBc {
if i > self.nx || j >= self.ny || k >= self.nz {
return FaceBc::Fluid;
}
read_face_bc(&self.u_face_bc, &self.u_solid, self.idx_u(i, j, k))
}
#[must_use]
pub fn v_bc(&self, i: usize, j: usize, k: usize) -> FaceBc {
if i >= self.nx || j > self.ny || k >= self.nz {
return FaceBc::Fluid;
}
read_face_bc(&self.v_face_bc, &self.v_solid, self.idx_v(i, j, k))
}
#[must_use]
pub fn w_bc(&self, i: usize, j: usize, k: usize) -> FaceBc {
if i >= self.nx || j >= self.ny || k > self.nz {
return FaceBc::Fluid;
}
read_face_bc(&self.w_face_bc, &self.w_solid, self.idx_w(i, j, k))
}
pub fn set_u_bc(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
if i > self.nx || j >= self.ny || k >= self.nz {
return;
}
let ix = self.idx_u(i, j, k);
self.u_solid[ix] = bc.is_wall();
write_face_bc(&mut self.u_face_bc, ix, bc);
}
pub fn set_v_bc(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
if i >= self.nx || j > self.ny || k >= self.nz {
return;
}
let ix = self.idx_v(i, j, k);
self.v_solid[ix] = bc.is_wall();
write_face_bc(&mut self.v_face_bc, ix, bc);
}
pub fn set_w_bc(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
if i >= self.nx || j >= self.ny || k > self.nz {
return;
}
let ix = self.idx_w(i, j, k);
self.w_solid[ix] = bc.is_wall();
write_face_bc(&mut self.w_face_bc, ix, bc);
}
pub fn set_closed_box_walls(&mut self) {
let wall = FaceBc::Wall {
velocity: Vec3Fix::ZERO,
};
let along = |n: usize| if n == 1 { FaceBc::SlipWall } else { wall };
let (bc_x, bc_y, bc_z) = (along(self.nx), along(self.ny), along(self.nz));
for k in 0..self.nz {
for j in 0..self.ny {
self.set_u_bc(0, j, k, bc_x);
self.set_u_bc(self.nx, j, k, bc_x);
}
}
for k in 0..self.nz {
for i in 0..self.nx {
self.set_v_bc(i, 0, k, bc_y);
self.set_v_bc(i, self.ny, k, bc_y);
}
}
for j in 0..self.ny {
for i in 0..self.nx {
self.set_w_bc(i, j, 0, bc_z);
self.set_w_bc(i, j, self.nz, bc_z);
}
}
}
pub fn enforce_solid_faces(&mut self) {
for (val, solid) in self.u.iter_mut().zip(self.u_solid.iter()) {
if *solid {
*val = Fix128::ZERO;
}
}
for (val, solid) in self.v.iter_mut().zip(self.v_solid.iter()) {
if *solid {
*val = Fix128::ZERO;
}
}
for (val, solid) in self.w.iter_mut().zip(self.w_solid.iter()) {
if *solid {
*val = Fix128::ZERO;
}
}
}
pub fn enforce_face_boundaries(&mut self) {
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..=self.nx {
let ix = self.idx_u(i, j, k);
match self.u_bc(i, j, k) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => self.u[ix] = normal_velocity,
_ => self.u[ix] = Fix128::ZERO,
}
}
}
}
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..=self.nx {
if matches!(self.u_bc(i, j, k), FaceBc::Outflow) {
let inner = if i > 0 { i - 1 } else { i + 1 };
let ix = self.idx_u(i, j, k);
self.u[ix] = self.u(inner, j, k);
}
}
}
}
for k in 0..self.nz {
for j in 0..=self.ny {
for i in 0..self.nx {
let ix = self.idx_v(i, j, k);
match self.v_bc(i, j, k) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => self.v[ix] = normal_velocity,
_ => self.v[ix] = Fix128::ZERO,
}
}
}
}
for k in 0..self.nz {
for j in 0..=self.ny {
for i in 0..self.nx {
if matches!(self.v_bc(i, j, k), FaceBc::Outflow) {
let inner = if j > 0 { j - 1 } else { j + 1 };
let ix = self.idx_v(i, j, k);
self.v[ix] = self.v(i, inner, k);
}
}
}
}
for k in 0..=self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
let ix = self.idx_w(i, j, k);
match self.w_bc(i, j, k) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => self.w[ix] = normal_velocity,
_ => self.w[ix] = Fix128::ZERO,
}
}
}
}
for k in 0..=self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
if matches!(self.w_bc(i, j, k), FaceBc::Outflow) {
let inner = if k > 0 { k - 1 } else { k + 1 };
let ix = self.idx_w(i, j, k);
self.w[ix] = self.w(i, j, inner);
}
}
}
}
}
#[inline]
#[must_use]
pub(crate) fn u_blocks_pressure(&self, i: usize, j: usize, k: usize) -> bool {
self.is_u_solid(i, j, k)
|| (!self.u_face_bc.is_empty() && self.u_bc(i, j, k).blocks_pressure())
}
#[inline]
#[must_use]
pub(crate) fn v_blocks_pressure(&self, i: usize, j: usize, k: usize) -> bool {
self.is_v_solid(i, j, k)
|| (!self.v_face_bc.is_empty() && self.v_bc(i, j, k).blocks_pressure())
}
#[inline]
#[must_use]
pub(crate) fn w_blocks_pressure(&self, i: usize, j: usize, k: usize) -> bool {
self.is_w_solid(i, j, k)
|| (!self.w_face_bc.is_empty() && self.w_bc(i, j, k).blocks_pressure())
}
#[inline]
#[must_use]
fn u_is_inflow(&self, i: usize, j: usize, k: usize) -> bool {
!self.u_face_bc.is_empty() && self.u_bc(i, j, k).is_inflow()
}
#[inline]
#[must_use]
fn v_is_inflow(&self, i: usize, j: usize, k: usize) -> bool {
!self.v_face_bc.is_empty() && self.v_bc(i, j, k).is_inflow()
}
#[inline]
#[must_use]
fn w_is_inflow(&self, i: usize, j: usize, k: usize) -> bool {
!self.w_face_bc.is_empty() && self.w_bc(i, j, k).is_inflow()
}
#[must_use]
pub fn u_wall_across_y(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let jf = if forward { j + 1 } else { j };
wall_between(
(i > 0).then(|| self.v_bc(i - 1, jf, k)),
(i < self.nx).then(|| self.v_bc(i, jf, k)),
)
}
#[must_use]
pub fn u_wall_across_z(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let kf = if forward { k + 1 } else { k };
wall_between(
(i > 0).then(|| self.w_bc(i - 1, j, kf)),
(i < self.nx).then(|| self.w_bc(i, j, kf)),
)
}
#[must_use]
pub fn v_wall_across_x(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let if_ = if forward { i + 1 } else { i };
wall_between(
(j > 0).then(|| self.u_bc(if_, j - 1, k)),
(j < self.ny).then(|| self.u_bc(if_, j, k)),
)
}
#[must_use]
pub fn v_wall_across_z(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let kf = if forward { k + 1 } else { k };
wall_between(
(j > 0).then(|| self.w_bc(i, j - 1, kf)),
(j < self.ny).then(|| self.w_bc(i, j, kf)),
)
}
#[must_use]
pub fn w_wall_across_x(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let if_ = if forward { i + 1 } else { i };
wall_between(
(k > 0).then(|| self.u_bc(if_, j, k - 1)),
(k < self.nz).then(|| self.u_bc(if_, j, k)),
)
}
#[must_use]
pub fn w_wall_across_y(&self, i: usize, j: usize, k: usize, forward: bool) -> Option<Vec3Fix> {
let jf = if forward { j + 1 } else { j };
wall_between(
(k > 0).then(|| self.v_bc(i, jf, k - 1)),
(k < self.nz).then(|| self.v_bc(i, jf, k)),
)
}
#[must_use]
pub fn divergence(&self, i: usize, j: usize, k: usize) -> Fix128 {
if self.dx.is_zero() {
return Fix128::ZERO;
}
let du = self.u(i + 1, j, k) - self.u(i, j, k);
let dv = self.v(i, j + 1, k) - self.v(i, j, k);
let dw = self.w(i, j, k + 1) - self.w(i, j, k);
(du + dv + dw) / self.dx
}
}
struct PoissonMask {
nx: usize,
ny: usize,
nz: usize,
open: Vec<[bool; 6]>,
}
impl PoissonMask {
fn from_grid(grid: &MacGrid) -> Self {
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
let mut open = vec![[true; 6]; nx * ny * nz];
for k in 0..nz {
for j in 0..ny {
for i in 0..nx {
open[i + nx * (j + ny * k)] = [
!grid.u_blocks_pressure(i, j, k),
!grid.u_blocks_pressure(i + 1, j, k),
!grid.v_blocks_pressure(i, j, k),
!grid.v_blocks_pressure(i, j + 1, k),
!grid.w_blocks_pressure(i, j, k),
!grid.w_blocks_pressure(i, j, k + 1),
];
}
}
}
Self { nx, ny, nz, open }
}
#[inline]
fn degree(&self, c: usize) -> i64 {
self.open[c].iter().filter(|&&o| o).count() as i64
}
#[inline]
fn neighbour_sum(&self, p: &[Fix128], i: usize, j: usize, k: usize) -> Fix128 {
let c = i + self.nx * (j + self.ny * k);
let o = self.open[c];
let mut acc = Fix128::ZERO;
if o[0] && i > 0 {
acc = acc + p[c - 1];
}
if o[1] && i + 1 < self.nx {
acc = acc + p[c + 1];
}
if o[2] && j > 0 {
acc = acc + p[c - self.nx];
}
if o[3] && j + 1 < self.ny {
acc = acc + p[c + self.nx];
}
if o[4] && k > 0 {
acc = acc + p[c - self.nx * self.ny];
}
if o[5] && k + 1 < self.nz {
acc = acc + p[c + self.nx * self.ny];
}
acc
}
fn apply(&self, p: &[Fix128], out: &mut [Fix128]) {
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
let c = i + self.nx * (j + self.ny * k);
out[c] =
p[c] * Fix128::from_int(-self.degree(c)) + self.neighbour_sum(p, i, j, k);
}
}
}
}
}
fn inverse_degrees(mask: &PoissonMask, n: usize) -> Vec<Fix128> {
let mut inv = vec![Fix128::ZERO; n];
for (c, slot) in inv.iter_mut().enumerate() {
let deg = mask.degree(c);
if deg > 0 {
*slot = Fix128::from_ratio(1, deg);
}
}
inv
}
fn subtract_pressure_gradient(grid: &mut MacGrid, coeff: Fix128) {
for k in 0..grid.nz {
for j in 0..grid.ny {
for i in 0..=grid.nx {
let ix = grid.idx_u(i, j, k);
if grid.u_solid[ix] {
grid.u[ix] = Fix128::ZERO;
continue;
}
if grid.u_is_inflow(i, j, k) {
continue;
}
let hi = if i < grid.nx {
grid.pressure(i, j, k)
} else {
Fix128::ZERO
};
let lo = if i > 0 {
grid.pressure(i - 1, j, k)
} else {
Fix128::ZERO
};
grid.u[ix] = grid.u[ix] - coeff * (hi - lo);
}
}
}
for k in 0..grid.nz {
for j in 0..=grid.ny {
for i in 0..grid.nx {
let ix = grid.idx_v(i, j, k);
if grid.v_solid[ix] {
grid.v[ix] = Fix128::ZERO;
continue;
}
if grid.v_is_inflow(i, j, k) {
continue;
}
let hi = if j < grid.ny {
grid.pressure(i, j, k)
} else {
Fix128::ZERO
};
let lo = if j > 0 {
grid.pressure(i, j - 1, k)
} else {
Fix128::ZERO
};
grid.v[ix] = grid.v[ix] - coeff * (hi - lo);
}
}
}
for k in 0..=grid.nz {
for j in 0..grid.ny {
for i in 0..grid.nx {
let ix = grid.idx_w(i, j, k);
if grid.w_solid[ix] {
grid.w[ix] = Fix128::ZERO;
continue;
}
if grid.w_is_inflow(i, j, k) {
continue;
}
let hi = if k < grid.nz {
grid.pressure(i, j, k)
} else {
Fix128::ZERO
};
let lo = if k > 0 {
grid.pressure(i, j, k - 1)
} else {
Fix128::ZERO
};
grid.w[ix] = grid.w[ix] - coeff * (hi - lo);
}
}
}
}
fn poisson_rhs(grid: &MacGrid, scale: Fix128) -> Vec<Fix128> {
let mut rhs = vec![Fix128::ZERO; grid.nx * grid.ny * grid.nz];
for k in 0..grid.nz {
for j in 0..grid.ny {
for i in 0..grid.nx {
rhs[i + grid.nx * (j + grid.ny * k)] = grid.divergence(i, j, k) * scale;
}
}
}
rhs
}
pub fn project_pressure(grid: &mut MacGrid, dt_s: Fix128, density_kg_m3: Fix128, iterations: u32) {
project_pressure_red_black_gs(grid, dt_s, density_kg_m3, iterations);
}
mod multigrid_decomposed;
#[cfg(feature = "std")]
pub(crate) use multigrid_decomposed::project_pressure_multigrid_decomposed;
#[cfg(feature = "std")]
mod rank_threads;
#[cfg(feature = "std")]
pub use rank_threads::{
project_pressure_distributed, project_pressure_distributed_with_report,
DistributedProjectionReport,
};
const MG_PRE_SMOOTH: u32 = 1;
const MG_POST_SMOOTH: u32 = 1;
const MG_CORRECTION_SCALE_NUM: i64 = 2;
const MG_CORRECTION_SCALE_DEN: i64 = 1;
const MG_COARSE_VISITS: u32 = 2;
struct MgLevel {
nx: usize,
ny: usize,
nz: usize,
cond: Vec<[i64; 6]>,
}
impl MgLevel {
fn from_mask(mask: &PoissonMask) -> Self {
Self {
nx: mask.nx,
ny: mask.ny,
nz: mask.nz,
cond: mask.open.iter().map(|o| o.map(i64::from)).collect(),
}
}
fn cells(&self) -> usize {
self.nx * self.ny * self.nz
}
fn degree(&self, c: usize) -> i64 {
self.cond[c].iter().sum()
}
fn coarsen(&self) -> Self {
let fx = if self.nx > 1 { 2 } else { 1 };
let fy = if self.ny > 1 { 2 } else { 1 };
let fz = if self.nz > 1 { 2 } else { 1 };
let (cx, cy, cz) = (self.nx / fx, self.ny / fy, self.nz / fz);
let mut cond = vec![[0i64; 6]; cx * cy * cz];
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
let c = i + self.nx * (j + self.ny * k);
let (ci, cj, ck) = (i / fx, j / fy, k / fz);
let target = &mut cond[ci + cx * (cj + cy * ck)];
let f = self.cond[c];
if i % fx == 0 {
target[0] += f[0];
}
if i % fx == fx - 1 {
target[1] += f[1];
}
if j % fy == 0 {
target[2] += f[2];
}
if j % fy == fy - 1 {
target[3] += f[3];
}
if k % fz == 0 {
target[4] += f[4];
}
if k % fz == fz - 1 {
target[5] += f[5];
}
}
}
}
Self {
nx: cx,
ny: cy,
nz: cz,
cond,
}
}
#[inline]
fn neighbour_sum(&self, p: &[Fix128], i: usize, j: usize, k: usize) -> Fix128 {
let c = i + self.nx * (j + self.ny * k);
let f = self.cond[c];
let mut acc = Fix128::ZERO;
if f[0] != 0 && i > 0 {
acc = acc + p[c - 1] * Fix128::from_int(f[0]);
}
if f[1] != 0 && i + 1 < self.nx {
acc = acc + p[c + 1] * Fix128::from_int(f[1]);
}
if f[2] != 0 && j > 0 {
acc = acc + p[c - self.nx] * Fix128::from_int(f[2]);
}
if f[3] != 0 && j + 1 < self.ny {
acc = acc + p[c + self.nx] * Fix128::from_int(f[3]);
}
if f[4] != 0 && k > 0 {
acc = acc + p[c - self.nx * self.ny] * Fix128::from_int(f[4]);
}
if f[5] != 0 && k + 1 < self.nz {
acc = acc + p[c + self.nx * self.ny] * Fix128::from_int(f[5]);
}
acc
}
fn inverse_degrees(&self) -> Vec<Fix128> {
(0..self.cells())
.map(|c| {
let deg = self.degree(c);
if deg > 0 {
Fix128::from_ratio(1, deg)
} else {
Fix128::ZERO
}
})
.collect()
}
fn smooth(&self, inv_deg: &[Fix128], p: &mut [Fix128], rhs: &[Fix128], iterations: u32) {
for _ in 0..iterations {
for colour in 0..2usize {
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
if (i + j + k) % 2 != colour {
continue;
}
let c = i + self.nx * (j + self.ny * k);
let nb = self.neighbour_sum(p, i, j, k);
p[c] = (nb - rhs[c]) * inv_deg[c];
}
}
}
}
}
}
fn residual(&self, p: &[Fix128], rhs: &[Fix128], res: &mut [Fix128]) {
for k in 0..self.nz {
for j in 0..self.ny {
for i in 0..self.nx {
let c = i + self.nx * (j + self.ny * k);
let ap =
p[c] * Fix128::from_int(-self.degree(c)) + self.neighbour_sum(p, i, j, k);
res[c] = rhs[c] - ap;
}
}
}
}
}
fn mg_vcycle(
levels: &[MgLevel],
invs: &[Vec<Fix128>],
ps: &mut [Vec<Fix128>],
rhss: &mut [Vec<Fix128>],
) {
let fine = &levels[0];
let (Some((p, p_rest)), Some((rhs, rhs_rest))) = (ps.split_first_mut(), rhss.split_first_mut())
else {
return;
};
if levels.len() == 1 {
for x in p.iter_mut() {
*x = Fix128::ZERO;
}
fine.smooth(&invs[0], p, rhs, 1);
return;
}
let coarse = &levels[1];
fine.smooth(&invs[0], p, rhs, MG_PRE_SMOOTH);
let mut res = vec![Fix128::ZERO; fine.cells()];
fine.residual(p, rhs, &mut res);
let (fx, fy, fz) = (
fine.nx / coarse.nx,
fine.ny / coarse.ny,
fine.nz / coarse.nz,
);
let rc = &mut rhs_rest[0];
for x in rc.iter_mut() {
*x = Fix128::ZERO;
}
for k in 0..fine.nz {
for j in 0..fine.ny {
for i in 0..fine.nx {
let ci = i / fx + coarse.nx * (j / fy + coarse.ny * (k / fz));
rc[ci] = rc[ci] + res[i + fine.nx * (j + fine.ny * k)];
}
}
}
for x in p_rest[0].iter_mut() {
*x = Fix128::ZERO;
}
for _ in 0..MG_COARSE_VISITS {
mg_vcycle(&levels[1..], &invs[1..], p_rest, rhs_rest);
}
let scale = Fix128::from_ratio(MG_CORRECTION_SCALE_NUM, MG_CORRECTION_SCALE_DEN);
let e = &p_rest[0];
for k in 0..fine.nz {
for j in 0..fine.ny {
for i in 0..fine.nx {
let ci = i / fx + coarse.nx * (j / fy + coarse.ny * (k / fz));
let c = i + fine.nx * (j + fine.ny * k);
p[c] = p[c] + e[ci] * scale;
}
}
}
fine.smooth(&invs[0], p, rhs, MG_POST_SMOOTH);
}
pub fn project_pressure_multigrid(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
cycles: u32,
) {
let pow2 = |n: usize| n.is_power_of_two();
if cycles == 0
|| !(pow2(grid.nx) && pow2(grid.ny) && pow2(grid.nz))
|| grid.dx.is_zero()
|| density_kg_m3.is_zero()
|| dt_s.is_zero()
{
return;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let rhs0 = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let mut levels = vec![MgLevel::from_mask(&mask)];
while let Some(next) = levels
.last()
.filter(|l| l.cells() > 1)
.map(MgLevel::coarsen)
{
levels.push(next);
}
let invs: Vec<Vec<Fix128>> = levels.iter().map(MgLevel::inverse_degrees).collect();
let mut ps: Vec<Vec<Fix128>> = levels
.iter()
.map(|l| vec![Fix128::ZERO; l.cells()])
.collect();
let mut rhss: Vec<Vec<Fix128>> = levels
.iter()
.map(|l| vec![Fix128::ZERO; l.cells()])
.collect();
ps[0].copy_from_slice(&grid.pressure);
rhss[0].copy_from_slice(&rhs0);
for _ in 0..cycles {
mg_vcycle(&levels, &invs, &mut ps, &mut rhss);
}
grid.pressure.copy_from_slice(&ps[0]);
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
}
pub(crate) fn project_pressure_red_black_gs(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
) {
if grid.dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() {
return;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let n = grid.nx * grid.ny * grid.nz;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
for _ in 0..iterations {
for colour in 0..2u32 {
for k in 0..grid.nz {
for j in 0..grid.ny {
for i in 0..grid.nx {
if ((i + j + k) as u32 % 2) != colour {
continue;
}
let idx = i + grid.nx * (j + grid.ny * k);
let neighbours = mask.neighbour_sum(&grid.pressure, i, j, k);
grid.pressure[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
}
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
}
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub(crate) enum HaloSchedule {
EverySweep,
EveryIteration,
}
fn slab_bounds(nz: usize, ranks: usize, rank: usize) -> (usize, usize) {
(rank * nz / ranks, (rank + 1) * nz / ranks)
}
fn slab_owner(bounds: &[(usize, usize)], k: usize) -> Option<usize> {
bounds.iter().position(|&(k0, k1)| k0 <= k && k < k1)
}
fn poison_beyond_halo(
buf: &mut [Fix128],
nz: usize,
plane: usize,
(k0, k1): (usize, usize),
sentinel: Fix128,
) {
let lo = k0.saturating_sub(1);
let hi = (k1 + 1).min(nz);
for k in 0..nz {
if k >= lo && k < hi {
continue;
}
let base = k * plane;
for slot in &mut buf[base..base + plane] {
*slot = sentinel;
}
}
}
pub(crate) trait RankTransport {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128];
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize);
}
pub(crate) struct LocalTransport {
plane: usize,
slabs: Vec<Vec<Fix128>>,
}
impl LocalTransport {
pub(crate) fn new(ranks: usize, nz: usize, plane: usize) -> Self {
Self {
plane,
slabs: vec![vec![Fix128::ZERO; nz * plane]; ranks],
}
}
}
impl RankTransport for LocalTransport {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
&mut self.slabs[rank]
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
let base = layer * self.plane;
let incoming = self.slabs[src][base..base + self.plane].to_vec();
self.slabs[dst][base..base + self.plane].copy_from_slice(&incoming);
}
}
fn exchange_slab_halos<T: RankTransport>(transport: &mut T, bounds: &[(usize, usize)], nz: usize) {
for (r, &(k0, k1)) in bounds.iter().enumerate() {
if k0 == k1 {
continue; }
for layer in [k0.checked_sub(1), (k1 < nz).then_some(k1)]
.into_iter()
.flatten()
{
let Some(src) = slab_owner(bounds, layer) else {
continue;
};
transport.deliver_layer(src, r, layer);
}
}
}
pub(crate) fn project_pressure_decomposed(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
ranks: usize,
schedule: HaloSchedule,
) {
let mut transport = LocalTransport::new(ranks, grid.nz, grid.nx * grid.ny);
project_pressure_decomposed_over(
grid,
dt_s,
density_kg_m3,
iterations,
ranks,
schedule,
&mut transport,
);
}
pub(crate) fn project_pressure_banded(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
ranks: usize,
) {
if grid.dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() || ranks == 0 {
return;
}
let nz = grid.nz;
let (nx, ny) = (grid.nx, grid.ny);
let plane = cell_plane(nx, ny);
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let mut faces: Vec<SlabFaces> = bounds
.iter()
.map(|&b| SlabFaces::from_grid(grid, b))
.collect();
let conditions: Vec<SlabFaceConditions> = bounds
.iter()
.map(|&b| SlabFaceConditions::from_grid(grid, b))
.collect();
enforce_slab_face_boundaries_over(&mut faces, &conditions);
let mut transport = LocalSlabTransport::from_field(&bounds, nz, plane, 1, &grid.pressure);
project_pressure_slab_local_over(
&mut faces,
dt_s,
density_kg_m3,
iterations,
HaloSchedule::EverySweep,
&mut transport,
);
for (r, rank_faces) in faces.iter().enumerate() {
write_band_back(grid, rank_faces, transport.slab(r));
}
}
fn write_band_back(grid: &mut MacGrid, rank_faces: &SlabFaces, band: &SlabStorage) {
let (nx, ny) = (grid.nx, grid.ny);
let plane = cell_plane(nx, ny);
let (k0, k1) = rank_faces.owned();
for k in k0..k1 {
let layer = band
.layer(k)
.expect("a rank's own layer is resident in its band");
grid.pressure[k * plane..(k + 1) * plane].copy_from_slice(layer);
let (u, _) = rank_faces.u_layer(k);
for j in 0..ny {
for i in 0..=nx {
let ix = grid.idx_u(i, j, k);
grid.u[ix] = u[i + (nx + 1) * j];
}
}
let (v, _) = rank_faces.v_layer(k);
for j in 0..=ny {
for i in 0..nx {
let ix = grid.idx_v(i, j, k);
grid.v[ix] = v[i + nx * j];
}
}
}
for k in rank_faces.w_written() {
let (w, _) = rank_faces.w_layer(k);
for j in 0..ny {
for i in 0..nx {
let ix = grid.idx_w(i, j, k);
grid.w[ix] = w[i + nx * j];
}
}
}
}
fn project_pressure_decomposed_over<T: RankTransport>(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
ranks: usize,
schedule: HaloSchedule,
transport: &mut T,
) {
if grid.dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() || ranks == 0 {
return;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
let n = nx * ny * nz;
let plane = nx * ny;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
let sentinel = Fix128::from_int(1_000_000);
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
for (r, &b) in bounds.iter().enumerate() {
let buf = transport.slab_mut(r);
assert_eq!(
buf.len(),
n,
"transport handed rank {r} a slab of {} cells instead of the {n} this stage needs",
buf.len(),
);
buf.copy_from_slice(&grid.pressure);
poison_beyond_halo(buf, nz, plane, b, sentinel);
}
for _ in 0..iterations {
for colour in 0..2u32 {
for (r, &(k0, k1)) in bounds.iter().enumerate() {
let buf = transport.slab_mut(r);
for k in k0..k1 {
for j in 0..ny {
for i in 0..nx {
if ((i + j + k) as u32 % 2) != colour {
continue;
}
let idx = i + nx * (j + ny * k);
let neighbours = mask.neighbour_sum(buf, i, j, k);
buf[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
}
if schedule == HaloSchedule::EverySweep {
exchange_slab_halos(transport, &bounds, nz);
}
}
if schedule == HaloSchedule::EveryIteration {
exchange_slab_halos(transport, &bounds, nz);
}
}
for (r, &(k0, k1)) in bounds.iter().enumerate() {
let buf = transport.slab_mut(r);
for k in k0..k1 {
let base = k * plane;
grid.pressure[base..base + plane].copy_from_slice(&buf[base..base + plane]);
}
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
}
fn gather_slabs_to_root<T: RankTransport>(transport: &mut T, bounds: &[(usize, usize)]) {
for (r, &(k0, k1)) in bounds.iter().enumerate().skip(1) {
for layer in k0..k1 {
transport.deliver_layer(r, 0, layer);
}
}
}
#[allow(clippy::too_many_arguments)]
pub(crate) fn project_pressure_decomposed_on_rank<T: RankTransport>(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
ranks: usize,
schedule: HaloSchedule,
my_rank: usize,
transport: &mut T,
) {
if grid.dx.is_zero()
|| density_kg_m3.is_zero()
|| dt_s.is_zero()
|| ranks == 0
|| my_rank >= ranks
{
return;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
let n = nx * ny * nz;
let plane = nx * ny;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
let sentinel = Fix128::from_int(1_000_000);
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let (k0, k1) = bounds[my_rank];
{
let buf = transport.slab_mut(my_rank);
assert_eq!(
buf.len(),
n,
"transport handed rank {my_rank} a slab of {} cells instead of the {n} this stage needs",
buf.len(),
);
buf.copy_from_slice(&grid.pressure);
poison_beyond_halo(buf, nz, plane, (k0, k1), sentinel);
}
for _ in 0..iterations {
for colour in 0..2u32 {
{
let buf = transport.slab_mut(my_rank);
for k in k0..k1 {
for j in 0..ny {
for i in 0..nx {
if ((i + j + k) as u32 % 2) != colour {
continue;
}
let idx = i + nx * (j + ny * k);
let neighbours = mask.neighbour_sum(buf, i, j, k);
buf[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
}
if schedule == HaloSchedule::EverySweep {
exchange_slab_halos(transport, &bounds, nz);
}
}
if schedule == HaloSchedule::EveryIteration {
exchange_slab_halos(transport, &bounds, nz);
}
}
gather_slabs_to_root(transport, &bounds);
if my_rank == 0 {
{
let buf = transport.slab_mut(0);
grid.pressure.copy_from_slice(buf);
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
}
}
#[cfg(feature = "std")]
const WIRE_BYTES_PER_CELL: usize = 16;
#[cfg(feature = "std")]
pub(crate) struct SocketTransport<S> {
my_rank: usize,
plane: usize,
buf: Vec<Fix128>,
links: Vec<Option<S>>,
wire: Vec<u8>,
}
#[cfg(feature = "std")]
impl<S> SocketTransport<S> {
pub(crate) fn new(my_rank: usize, cells: usize, plane: usize, links: Vec<Option<S>>) -> Self {
Self {
my_rank,
plane,
buf: vec![Fix128::ZERO; cells],
links,
wire: vec![0u8; plane * WIRE_BYTES_PER_CELL],
}
}
}
#[cfg(feature = "std")]
impl<S: std::io::Read + std::io::Write> RankTransport for SocketTransport<S> {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
assert_eq!(
rank, self.my_rank,
"rank {} was asked for rank {rank}'s slab, which is in another address \
space: the driver is not rank-local",
self.my_rank,
);
&mut self.buf
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
if src == dst {
return;
}
let base = layer * self.plane;
if self.my_rank == src {
for (cell, chunk) in self.buf[base..base + self.plane]
.iter()
.zip(self.wire.chunks_exact_mut(WIRE_BYTES_PER_CELL))
{
chunk[..8].copy_from_slice(&cell.hi.to_le_bytes());
chunk[8..].copy_from_slice(&cell.lo.to_le_bytes());
}
let link = self.links[dst]
.as_mut()
.expect("no stream to the rank this delivery is addressed to");
link.write_all(&self.wire)
.expect("send a halo layer to the peer rank");
link.flush().expect("flush a halo layer to the peer rank");
} else if self.my_rank == dst {
let link = self.links[src]
.as_mut()
.expect("no stream from the rank this delivery comes from");
link.read_exact(&mut self.wire)
.expect("receive a halo layer from the peer rank");
for (cell, chunk) in self.buf[base..base + self.plane]
.iter_mut()
.zip(self.wire.chunks_exact(WIRE_BYTES_PER_CELL))
{
let hi =
i64::from_le_bytes(chunk[..8].try_into().expect("8 bytes of a 16-byte cell"));
let lo =
u64::from_le_bytes(chunk[8..].try_into().expect("8 bytes of a 16-byte cell"));
*cell = Fix128::from_raw(hi, lo);
}
}
}
}
pub(crate) const HALO_MISSING_BELOW: &str =
"slab halo too narrow: an open -z face of an owned cell reads the layer below the slab, \
which is not resident";
pub(crate) const HALO_MISSING_ABOVE: &str =
"slab halo too narrow: an open +z face of an owned cell reads the layer above the slab, \
which is not resident";
pub(crate) const DELIVERY_SOURCE_LACKS_LAYER: &str =
"slab delivery: the sending rank does not hold the layer it was asked to send";
pub(crate) const DELIVERY_DESTINATION_LACKS_LAYER: &str =
"slab delivery: the receiving rank has no room for the layer it was asked to receive, \
so its halo is narrower than the exchange schedule";
#[inline]
fn cell_plane(nx: usize, ny: usize) -> usize {
nx * ny
}
#[inline]
fn u_plane(nx: usize, ny: usize) -> usize {
(nx + 1) * ny
}
#[inline]
fn v_plane(nx: usize, ny: usize) -> usize {
nx * (ny + 1)
}
pub(crate) struct SlabStorage {
plane: usize,
lo: usize,
hi: usize,
cells: Vec<Fix128>,
}
impl SlabStorage {
pub(crate) fn for_slab(plane: usize, nz: usize, (k0, k1): (usize, usize), halo: usize) -> Self {
let (lo, hi) = if k0 == k1 {
(k0, k0)
} else {
(k0.saturating_sub(halo), (k1 + halo).min(nz))
};
Self {
plane,
lo,
hi,
cells: vec![Fix128::ZERO; (hi - lo) * plane],
}
}
pub(crate) fn resident(&self) -> (usize, usize) {
(self.lo, self.hi)
}
pub(crate) fn layer(&self, k: usize) -> Option<&[Fix128]> {
let row = k.checked_sub(self.lo)?;
let base = row.checked_mul(self.plane)?;
self.cells.get(base..base + self.plane)
}
pub(crate) fn layer_mut(&mut self, k: usize) -> Option<&mut [Fix128]> {
let row = k.checked_sub(self.lo)?;
let base = row.checked_mul(self.plane)?;
self.cells.get_mut(base..base + self.plane)
}
pub(crate) fn sweep_window(&mut self, k: usize) -> SweepWindow<'_> {
let row = k.checked_sub(self.lo).filter(|r| self.lo + r < self.hi);
let row = row.unwrap_or_else(|| {
panic!(
"a sweep asked for layer {k}, which rank's slab does not hold (resident {}..{})",
self.lo, self.hi,
)
});
let plane = self.plane;
let (lower, rest) = self.cells.split_at_mut(row * plane);
let (centre, upper) = rest.split_at_mut(plane);
SweepWindow {
below: (row > 0).then(|| &lower[(row - 1) * plane..row * plane]),
centre,
above: (!upper.is_empty()).then(|| &upper[..plane]),
}
}
pub(crate) fn bytes(&self) -> SlabBytes {
SlabBytes {
pressure: self.cells.capacity() * size_of::<Fix128>(),
..SlabBytes::ZERO
}
}
}
pub(crate) struct SweepWindow<'a> {
below: Option<&'a [Fix128]>,
centre: &'a mut [Fix128],
above: Option<&'a [Fix128]>,
}
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub(crate) struct FaceFlags {
solid: bool,
blocks_pressure: bool,
inflow: bool,
}
impl FaceFlags {
fn of_u(grid: &MacGrid, i: usize, j: usize, k: usize) -> Self {
Self {
solid: grid.u_solid[grid.idx_u(i, j, k)],
blocks_pressure: grid.u_blocks_pressure(i, j, k),
inflow: grid.u_is_inflow(i, j, k),
}
}
fn of_v(grid: &MacGrid, i: usize, j: usize, k: usize) -> Self {
Self {
solid: grid.v_solid[grid.idx_v(i, j, k)],
blocks_pressure: grid.v_blocks_pressure(i, j, k),
inflow: grid.v_is_inflow(i, j, k),
}
}
fn of_w(grid: &MacGrid, i: usize, j: usize, k: usize) -> Self {
Self {
solid: grid.w_solid[grid.idx_w(i, j, k)],
blocks_pressure: grid.w_blocks_pressure(i, j, k),
inflow: grid.w_is_inflow(i, j, k),
}
}
}
pub(crate) struct SlabFaces {
nx: usize,
ny: usize,
nz: usize,
dx: Fix128,
k0: usize,
k1: usize,
u: Vec<Fix128>,
u_flags: Vec<FaceFlags>,
v: Vec<Fix128>,
v_flags: Vec<FaceFlags>,
w: Vec<Fix128>,
w_flags: Vec<FaceFlags>,
}
impl SlabFaces {
pub(crate) fn new(
nx: usize,
ny: usize,
nz: usize,
dx: Fix128,
(k0, k1): (usize, usize),
) -> Self {
let owned = k1 - k0;
let w_layers = if owned == 0 { 0 } else { owned + 1 };
Self {
nx,
ny,
nz,
dx,
k0,
k1,
u: vec![Fix128::ZERO; owned * u_plane(nx, ny)],
u_flags: vec![FaceFlags::default(); owned * u_plane(nx, ny)],
v: vec![Fix128::ZERO; owned * v_plane(nx, ny)],
v_flags: vec![FaceFlags::default(); owned * v_plane(nx, ny)],
w: vec![Fix128::ZERO; w_layers * cell_plane(nx, ny)],
w_flags: vec![FaceFlags::default(); w_layers * cell_plane(nx, ny)],
}
}
pub(crate) fn from_grid(grid: &MacGrid, (k0, k1): (usize, usize)) -> Self {
let mut faces = Self::new(grid.nx, grid.ny, grid.nz, grid.dx, (k0, k1));
if k0 == k1 {
return faces;
}
let (nx, ny) = (grid.nx, grid.ny);
for k in k0..k1 {
let row = (k - k0) * u_plane(nx, ny);
for j in 0..ny {
for i in 0..=nx {
let at = row + i + (nx + 1) * j;
faces.u[at] = grid.u[grid.idx_u(i, j, k)];
faces.u_flags[at] = FaceFlags::of_u(grid, i, j, k);
}
}
let row = (k - k0) * v_plane(nx, ny);
for j in 0..=ny {
for i in 0..nx {
let at = row + i + nx * j;
faces.v[at] = grid.v[grid.idx_v(i, j, k)];
faces.v_flags[at] = FaceFlags::of_v(grid, i, j, k);
}
}
}
for k in k0..=k1 {
let row = (k - k0) * cell_plane(nx, ny);
for j in 0..ny {
for i in 0..nx {
let at = row + i + nx * j;
faces.w[at] = grid.w[grid.idx_w(i, j, k)];
faces.w_flags[at] = FaceFlags::of_w(grid, i, j, k);
}
}
}
faces
}
pub(crate) fn owned(&self) -> (usize, usize) {
(self.k0, self.k1)
}
fn owned_row(&self, k: usize, plane: usize) -> usize {
assert!(
k >= self.k0 && k < self.k1,
"layer {k} is not owned by this slab (owns {}..{})",
self.k0,
self.k1,
);
(k - self.k0) * plane
}
fn w_row(&self, k: usize) -> usize {
assert!(
k >= self.k0 && k <= self.k1 && self.k0 != self.k1,
"Z-face layer {k} is not held by this slab (holds {}..={})",
self.k0,
self.k1,
);
(k - self.k0) * cell_plane(self.nx, self.ny)
}
fn u_layer(&self, k: usize) -> (&[Fix128], &[FaceFlags]) {
let plane = u_plane(self.nx, self.ny);
let row = self.owned_row(k, plane);
(&self.u[row..row + plane], &self.u_flags[row..row + plane])
}
fn v_layer(&self, k: usize) -> (&[Fix128], &[FaceFlags]) {
let plane = v_plane(self.nx, self.ny);
let row = self.owned_row(k, plane);
(&self.v[row..row + plane], &self.v_flags[row..row + plane])
}
fn w_layer(&self, k: usize) -> (&[Fix128], &[FaceFlags]) {
let plane = cell_plane(self.nx, self.ny);
let row = self.w_row(k);
(&self.w[row..row + plane], &self.w_flags[row..row + plane])
}
pub(crate) fn u_layer_mut(&mut self, k: usize) -> (&mut [Fix128], &mut [FaceFlags]) {
let plane = u_plane(self.nx, self.ny);
let row = self.owned_row(k, plane);
(
&mut self.u[row..row + plane],
&mut self.u_flags[row..row + plane],
)
}
pub(crate) fn v_layer_mut(&mut self, k: usize) -> (&mut [Fix128], &mut [FaceFlags]) {
let plane = v_plane(self.nx, self.ny);
let row = self.owned_row(k, plane);
(
&mut self.v[row..row + plane],
&mut self.v_flags[row..row + plane],
)
}
pub(crate) fn w_layer_mut(&mut self, k: usize) -> (&mut [Fix128], &mut [FaceFlags]) {
let plane = cell_plane(self.nx, self.ny);
let row = self.w_row(k);
(
&mut self.w[row..row + plane],
&mut self.w_flags[row..row + plane],
)
}
fn w_written(&self) -> core::ops::Range<usize> {
if self.k0 == self.k1 {
return self.k0..self.k0;
}
let end = if self.k1 == self.nz {
self.nz + 1
} else {
self.k1
};
self.k0..end
}
pub(crate) fn bytes(&self) -> SlabBytes {
let values = self.u.capacity() + self.v.capacity() + self.w.capacity();
let flags = self.u_flags.capacity() + self.v_flags.capacity() + self.w_flags.capacity();
SlabBytes {
faces: values * size_of::<Fix128>(),
face_flags: flags * size_of::<FaceFlags>(),
..SlabBytes::ZERO
}
}
}
pub(crate) struct SlabStencil {
k0: usize,
k1: usize,
plane: usize,
open: Vec<[bool; 6]>,
inv_deg: Vec<Fix128>,
rhs: Vec<Fix128>,
}
impl SlabStencil {
pub(crate) fn build(faces: &SlabFaces, scale: Fix128) -> Self {
assert!(
!faces.dx.is_zero(),
"a slab stencil needs a non-zero cell spacing",
);
let (nx, ny) = (faces.nx, faces.ny);
let (k0, k1) = faces.owned();
let plane = cell_plane(nx, ny);
let cells = (k1 - k0) * plane;
let mut open = Vec::with_capacity(cells);
let mut inv_deg = Vec::with_capacity(cells);
let mut rhs = Vec::with_capacity(cells);
let u_row = nx + 1;
for k in k0..k1 {
let (u, u_flags) = faces.u_layer(k);
let (v, v_flags) = faces.v_layer(k);
let (w_lo, w_lo_flags) = faces.w_layer(k);
let (w_hi, w_hi_flags) = faces.w_layer(k + 1);
for j in 0..ny {
for i in 0..nx {
let o = [
!u_flags[i + u_row * j].blocks_pressure,
!u_flags[i + 1 + u_row * j].blocks_pressure,
!v_flags[i + nx * j].blocks_pressure,
!v_flags[i + nx * (j + 1)].blocks_pressure,
!w_lo_flags[i + nx * j].blocks_pressure,
!w_hi_flags[i + nx * j].blocks_pressure,
];
let degree = o.iter().filter(|&&face| face).count() as i64;
inv_deg.push(if degree > 0 {
Fix128::from_ratio(1, degree)
} else {
Fix128::ZERO
});
let du = u[i + 1 + u_row * j] - u[i + u_row * j];
let dv = v[i + nx * (j + 1)] - v[i + nx * j];
let dw = w_hi[i + nx * j] - w_lo[i + nx * j];
rhs.push((du + dv + dw) / faces.dx * scale);
open.push(o);
}
}
}
Self {
k0,
k1,
plane,
open,
inv_deg,
rhs,
}
}
fn row(&self, k: usize) -> usize {
assert!(
k >= self.k0 && k < self.k1,
"layer {k} is not owned by this stencil (owns {}..{})",
self.k0,
self.k1,
);
(k - self.k0) * self.plane
}
pub(crate) fn bytes(&self) -> SlabBytes {
SlabBytes {
open: self.open.capacity() * size_of::<[bool; 6]>(),
inverse_degrees: self.inv_deg.capacity() * size_of::<Fix128>(),
rhs: self.rhs.capacity() * size_of::<Fix128>(),
..SlabBytes::ZERO
}
}
}
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub(crate) struct SlabBytes {
pub(crate) pressure: usize,
pub(crate) faces: usize,
pub(crate) face_flags: usize,
pub(crate) open: usize,
pub(crate) inverse_degrees: usize,
pub(crate) rhs: usize,
}
impl SlabBytes {
pub(crate) const ZERO: Self = Self {
pressure: 0,
faces: 0,
face_flags: 0,
open: 0,
inverse_degrees: 0,
rhs: 0,
};
pub(crate) fn add(self, other: Self) -> Self {
Self {
pressure: self.pressure + other.pressure,
faces: self.faces + other.faces,
face_flags: self.face_flags + other.face_flags,
open: self.open + other.open,
inverse_degrees: self.inverse_degrees + other.inverse_degrees,
rhs: self.rhs + other.rhs,
}
}
pub(crate) fn total(self) -> usize {
self.pressure + self.faces + self.face_flags + self.open + self.inverse_degrees + self.rhs
}
}
pub(crate) trait SlabTransport {
fn slab_mut(&mut self, rank: usize) -> &mut SlabStorage;
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize);
}
pub(crate) struct LocalSlabTransport {
slabs: Vec<SlabStorage>,
staging: Vec<Fix128>,
}
impl LocalSlabTransport {
pub(crate) fn from_field(
bounds: &[(usize, usize)],
nz: usize,
plane: usize,
halo: usize,
field: &[Fix128],
) -> Self {
assert_eq!(
field.len(),
nz * plane,
"the field a slab decomposition starts from must cover the whole domain",
);
let mut slabs = Vec::with_capacity(bounds.len());
for &b in bounds {
let mut slab = SlabStorage::for_slab(plane, nz, b, halo);
let (lo, hi) = slab.resident();
for k in lo..hi {
let layer = slab
.layer_mut(k)
.expect("a layer inside the band this slab just reported");
layer.copy_from_slice(&field[k * plane..(k + 1) * plane]);
}
slabs.push(slab);
}
Self {
slabs,
staging: vec![Fix128::ZERO; plane],
}
}
pub(crate) fn slab(&self, rank: usize) -> &SlabStorage {
&self.slabs[rank]
}
}
impl SlabTransport for LocalSlabTransport {
fn slab_mut(&mut self, rank: usize) -> &mut SlabStorage {
&mut self.slabs[rank]
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
if src == dst {
return;
}
let from = self.slabs[src]
.layer(layer)
.expect(DELIVERY_SOURCE_LACKS_LAYER);
self.staging.copy_from_slice(from);
let into = self.slabs[dst]
.layer_mut(layer)
.expect(DELIVERY_DESTINATION_LACKS_LAYER);
into.copy_from_slice(&self.staging);
}
}
fn exchange_slab_halos_local<T: SlabTransport>(
transport: &mut T,
bounds: &[(usize, usize)],
nz: usize,
) {
for (r, &(k0, k1)) in bounds.iter().enumerate() {
if k0 == k1 {
continue; }
for layer in [k0.checked_sub(1), (k1 < nz).then_some(k1)]
.into_iter()
.flatten()
{
let Some(src) = slab_owner(bounds, layer) else {
continue;
};
transport.deliver_layer(src, r, layer);
}
}
}
fn sweep_slab_colour(
storage: &mut SlabStorage,
faces: &SlabFaces,
stencil: &SlabStencil,
colour: u32,
) {
let (nx, ny, nz) = (faces.nx, faces.ny, faces.nz);
let (k0, k1) = faces.owned();
for k in k0..k1 {
let row = stencil.row(k);
let SweepWindow {
below,
centre,
above,
} = storage.sweep_window(k);
for j in 0..ny {
for i in 0..nx {
if ((i + j + k) as u32 % 2) != colour {
continue;
}
let here = i + nx * j;
let cell = row + here;
let o = stencil.open[cell];
let mut acc = Fix128::ZERO;
if o[0] && i > 0 {
acc = acc + centre[here - 1];
}
if o[1] && i + 1 < nx {
acc = acc + centre[here + 1];
}
if o[2] && j > 0 {
acc = acc + centre[here - nx];
}
if o[3] && j + 1 < ny {
acc = acc + centre[here + nx];
}
if o[4] && k > 0 {
acc = acc + below.expect(HALO_MISSING_BELOW)[here];
}
if o[5] && k + 1 < nz {
acc = acc + above.expect(HALO_MISSING_ABOVE)[here];
}
centre[here] = (acc - stencil.rhs[cell]) * stencil.inv_deg[cell];
}
}
}
}
fn subtract_slab_pressure_gradient(faces: &mut SlabFaces, storage: &SlabStorage, coeff: Fix128) {
let (nx, ny, nz) = (faces.nx, faces.ny, faces.nz);
let (k0, k1) = faces.owned();
let u_row = nx + 1;
for k in k0..k1 {
let here = storage
.layer(k)
.expect("a rank's own layer, which its band holds by construction");
let (u, u_flags) = faces.u_layer_mut(k);
for j in 0..ny {
for i in 0..=nx {
let at = i + u_row * j;
if u_flags[at].solid {
u[at] = Fix128::ZERO;
continue;
}
if u_flags[at].inflow {
continue;
}
let hi = if i < nx {
here[i + nx * j]
} else {
Fix128::ZERO
};
let lo = if i > 0 {
here[i - 1 + nx * j]
} else {
Fix128::ZERO
};
u[at] = u[at] - coeff * (hi - lo);
}
}
let (v, v_flags) = faces.v_layer_mut(k);
for j in 0..=ny {
for i in 0..nx {
let at = i + nx * j;
if v_flags[at].solid {
v[at] = Fix128::ZERO;
continue;
}
if v_flags[at].inflow {
continue;
}
let hi = if j < ny {
here[i + nx * j]
} else {
Fix128::ZERO
};
let lo = if j > 0 {
here[i + nx * (j - 1)]
} else {
Fix128::ZERO
};
v[at] = v[at] - coeff * (hi - lo);
}
}
}
for k in faces.w_written() {
let above = if k < nz {
Some(
storage
.layer(k)
.expect("a rank's own layer, which its band holds by construction"),
)
} else {
None
};
let below = if k > 0 {
Some(storage.layer(k - 1).expect(HALO_MISSING_BELOW))
} else {
None
};
let (w, w_flags) = faces.w_layer_mut(k);
for j in 0..ny {
for i in 0..nx {
let at = i + nx * j;
if w_flags[at].solid {
w[at] = Fix128::ZERO;
continue;
}
if w_flags[at].inflow {
continue;
}
let hi = above.map_or(Fix128::ZERO, |layer| layer[at]);
let lo = below.map_or(Fix128::ZERO, |layer| layer[at]);
w[at] = w[at] - coeff * (hi - lo);
}
}
}
}
pub(crate) fn project_pressure_slab_local_over<T: SlabTransport>(
faces: &mut [SlabFaces],
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
schedule: HaloSchedule,
transport: &mut T,
) {
let Some(first) = faces.first() else {
return;
};
let (nx, ny, nz, dx) = (first.nx, first.ny, first.nz, first.dx);
if dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() {
return;
}
let ranks = faces.len();
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
for (r, rank_faces) in faces.iter().enumerate() {
assert_eq!(
rank_faces.owned(),
bounds[r],
"rank {r} was handed the faces of layers {:?} but the decomposition gives it {:?}",
rank_faces.owned(),
bounds[r],
);
assert_eq!(
(rank_faces.nx, rank_faces.ny, rank_faces.nz),
(nx, ny, nz),
"rank {r}'s faces describe a different grid from rank 0's",
);
}
let scale = density_kg_m3 * dx * dx / dt_s;
let stencils: Vec<SlabStencil> = faces
.iter()
.map(|rank_faces| SlabStencil::build(rank_faces, scale))
.collect();
for _ in 0..iterations {
for colour in 0..2u32 {
for r in 0..ranks {
sweep_slab_colour(transport.slab_mut(r), &faces[r], &stencils[r], colour);
}
if schedule == HaloSchedule::EverySweep {
exchange_slab_halos_local(transport, &bounds, nz);
}
}
if schedule == HaloSchedule::EveryIteration {
exchange_slab_halos_local(transport, &bounds, nz);
}
}
let inv_dx = Fix128::ONE / dx;
let coeff = dt_s / density_kg_m3 * inv_dx;
for (r, rank_faces) in faces.iter_mut().enumerate() {
subtract_slab_pressure_gradient(rank_faces, transport.slab_mut(r), coeff);
}
}
pub(crate) const SLAB_FACE_HALO_MISSING: &str =
"slab face conditions: the first Z-face of this rank is an outflow and copies the face \
one layer below the slab, which the rank below has to hand over";
pub(crate) struct SlabFaceConditions {
nx: usize,
ny: usize,
k0: usize,
k1: usize,
u: Vec<(usize, FaceBc)>,
v: Vec<(usize, FaceBc)>,
w: Vec<(usize, FaceBc)>,
}
impl SlabFaceConditions {
pub(crate) fn new(nx: usize, ny: usize, (k0, k1): (usize, usize)) -> Self {
Self {
nx,
ny,
k0,
k1,
u: Vec::new(),
v: Vec::new(),
w: Vec::new(),
}
}
pub(crate) fn from_grid(grid: &MacGrid, (k0, k1): (usize, usize)) -> Self {
let (nx, ny) = (grid.nx, grid.ny);
let mut out = Self::new(nx, ny, (k0, k1));
if k0 == k1 {
return out;
}
let (up, vp, cp) = (u_plane(nx, ny), v_plane(nx, ny), cell_plane(nx, ny));
for (&ix, &bc) in grid.u_face_bc.range(k0 * up..k1 * up) {
let at = ix % up;
out.set_u(at % (nx + 1), at / (nx + 1), ix / up, bc);
}
for (&ix, &bc) in grid.v_face_bc.range(k0 * vp..k1 * vp) {
let at = ix % vp;
out.set_v(at % nx, at / nx, ix / vp, bc);
}
for (&ix, &bc) in grid.w_face_bc.range(k0 * cp..(k1 + 1) * cp) {
let at = ix % cp;
out.set_w(at % nx, at / nx, ix / cp, bc);
}
out
}
pub(crate) fn set_u(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
assert!(i <= self.nx && j < self.ny, "({i}, {j}) is not an X-face");
let at = self.owned_base(k, u_plane(self.nx, self.ny)) + i + (self.nx + 1) * j;
Self::write(&mut self.u, at, bc);
}
pub(crate) fn set_v(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
assert!(i < self.nx && j <= self.ny, "({i}, {j}) is not a Y-face");
let at = self.owned_base(k, v_plane(self.nx, self.ny)) + i + self.nx * j;
Self::write(&mut self.v, at, bc);
}
pub(crate) fn set_w(&mut self, i: usize, j: usize, k: usize, bc: FaceBc) {
assert!(i < self.nx && j < self.ny, "({i}, {j}) is not a Z-face");
assert!(
k >= self.k0 && k <= self.k1 && self.k0 != self.k1,
"Z-face layer {k} is not held by this slab (holds {}..={})",
self.k0,
self.k1,
);
let at = (k - self.k0) * cell_plane(self.nx, self.ny) + i + self.nx * j;
Self::write(&mut self.w, at, bc);
}
fn owned_base(&self, k: usize, plane: usize) -> usize {
assert!(
k >= self.k0 && k < self.k1,
"layer {k} is not owned by this slab (owns {}..{})",
self.k0,
self.k1,
);
(k - self.k0) * plane
}
fn write(entries: &mut Vec<(usize, FaceBc)>, at: usize, bc: FaceBc) {
let default_wall = FaceBc::Wall {
velocity: Vec3Fix::ZERO,
};
let slot = entries.partition_point(|&(other, _)| other < at);
let present = entries.get(slot).is_some_and(|&(other, _)| other == at);
if bc == FaceBc::Fluid || bc == default_wall {
if present {
entries.remove(slot);
}
} else if present {
entries[slot].1 = bc;
} else {
entries.insert(slot, (at, bc));
}
}
}
fn layer_conditions(entries: &[(usize, FaceBc)], base: usize, plane: usize) -> &[(usize, FaceBc)] {
let lo = entries.partition_point(|&(at, _)| at < base);
let hi = entries.partition_point(|&(at, _)| at < base + plane);
&entries[lo..hi]
}
#[inline]
fn condition_at(entries: &[(usize, FaceBc)], next: &mut usize, at: usize, solid: bool) -> FaceBc {
while *next < entries.len() && entries[*next].0 < at {
*next += 1;
}
match entries.get(*next) {
Some(&(other, bc)) if other == at => bc,
_ => wall_or_fluid(solid),
}
}
pub(crate) fn enforce_slab_face_boundaries(
faces: &mut SlabFaces,
cond: &SlabFaceConditions,
below_w: Option<&[Fix128]>,
) {
let (k0, k1) = faces.owned();
assert_eq!(
(cond.k0, cond.k1),
(k0, k1),
"the conditions describe layers {:?} and the faces {:?}",
(cond.k0, cond.k1),
(k0, k1),
);
assert_eq!(
(cond.nx, cond.ny),
(faces.nx, faces.ny),
"the conditions describe a different grid from the faces",
);
if k0 == k1 {
return; }
let (nx, ny) = (faces.nx, faces.ny);
let u_row = nx + 1;
let up = u_plane(nx, ny);
for k in k0..k1 {
let base = (k - k0) * up;
let entries = layer_conditions(&cond.u, base, up);
let (u, flags) = faces.u_layer_mut(k);
let mut next = 0usize;
for j in 0..ny {
for i in 0..=nx {
let at = i + u_row * j;
match condition_at(entries, &mut next, base + at, flags[at].solid) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => u[at] = normal_velocity,
_ => u[at] = Fix128::ZERO,
}
}
}
for &(entry, bc) in entries {
if bc == FaceBc::Outflow {
let at = entry - base;
let (i, j) = (at % u_row, at / u_row);
let inner = if i > 0 { i - 1 } else { i + 1 };
u[at] = if inner <= nx {
u[inner + u_row * j]
} else {
Fix128::ZERO
};
}
}
}
let vp = v_plane(nx, ny);
for k in k0..k1 {
let base = (k - k0) * vp;
let entries = layer_conditions(&cond.v, base, vp);
let (v, flags) = faces.v_layer_mut(k);
let mut next = 0usize;
for j in 0..=ny {
for i in 0..nx {
let at = i + nx * j;
match condition_at(entries, &mut next, base + at, flags[at].solid) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => v[at] = normal_velocity,
_ => v[at] = Fix128::ZERO,
}
}
}
for &(entry, bc) in entries {
if bc == FaceBc::Outflow {
let at = entry - base;
let (i, j) = (at % nx, at / nx);
let inner = if j > 0 { j - 1 } else { j + 1 };
v[at] = if inner <= ny {
v[i + nx * inner]
} else {
Fix128::ZERO
};
}
}
}
let cp = cell_plane(nx, ny);
for k in k0..=k1 {
let base = (k - k0) * cp;
let entries = layer_conditions(&cond.w, base, cp);
let (w, flags) = faces.w_layer_mut(k);
let mut next = 0usize;
for j in 0..ny {
for i in 0..nx {
let at = i + nx * j;
match condition_at(entries, &mut next, base + at, flags[at].solid) {
FaceBc::Fluid | FaceBc::Outflow => {}
FaceBc::Inflow { normal_velocity } => w[at] = normal_velocity,
_ => w[at] = Fix128::ZERO,
}
}
}
}
let mut source = Vec::new();
for k in k0..=k1 {
let base = (k - k0) * cp;
let entries = layer_conditions(&cond.w, base, cp);
if !entries.iter().any(|&(_, bc)| bc == FaceBc::Outflow) {
continue;
}
if source.is_empty() {
source = vec![Fix128::ZERO; cp];
}
let inner = if k > 0 { k - 1 } else { k + 1 };
if inner < k0 {
source.copy_from_slice(below_w.expect(SLAB_FACE_HALO_MISSING));
} else {
assert!(
inner <= k1,
"the face below Z-face layer {k} is layer {inner}, which is above the \
band {k0}..={k1}: only a rank owning nothing could reach this, and one \
has already returned",
);
source.copy_from_slice(faces.w_layer(inner).0);
}
let (w, _) = faces.w_layer_mut(k);
for &(entry, bc) in entries {
if bc == FaceBc::Outflow {
let at = entry - base;
w[at] = source[at];
}
}
}
}
pub(crate) trait PlaneChannel {
fn send_plane(&mut self, dst: usize, layer: usize, plane: &[Fix128]);
fn recv_plane(&mut self, src: usize, layer: usize, plane: &mut [Fix128]);
}
pub(crate) struct LocalPlaneChannel {
staged: Vec<Fix128>,
layer: Option<usize>,
}
impl LocalPlaneChannel {
pub(crate) fn new(plane: usize) -> Self {
Self {
staged: vec![Fix128::ZERO; plane],
layer: None,
}
}
pub(crate) fn is_drained(&self) -> bool {
self.layer.is_none()
}
}
impl PlaneChannel for LocalPlaneChannel {
fn send_plane(&mut self, _dst: usize, layer: usize, plane: &[Fix128]) {
assert!(
self.layer.is_none(),
"layer {layer} was offered while layer {:?} was still in flight",
self.layer,
);
self.staged.copy_from_slice(plane);
self.layer = Some(layer);
}
fn recv_plane(&mut self, _src: usize, layer: usize, plane: &mut [Fix128]) {
let staged = self
.layer
.take()
.expect("a rank asked for a plane that no rank had offered");
assert_eq!(
staged, layer,
"a rank asked for layer {layer} and layer {staged} was in flight",
);
plane.copy_from_slice(&self.staged);
}
}
pub(crate) fn enforce_slab_face_boundaries_on_rank<C: PlaneChannel>(
faces: &mut SlabFaces,
cond: &SlabFaceConditions,
ranks: usize,
my_rank: usize,
channel: &mut C,
) {
if ranks == 0 || my_rank >= ranks {
return;
}
let nz = faces.nz;
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
enforce_slab_face_boundaries_in_chain(faces, cond, &bounds, my_rank, channel);
}
fn enforce_slab_face_boundaries_in_chain<C: PlaneChannel>(
faces: &mut SlabFaces,
cond: &SlabFaceConditions,
bounds: &[(usize, usize)],
my_rank: usize,
channel: &mut C,
) {
let (k0, k1) = bounds[my_rank];
assert_eq!(
faces.owned(),
(k0, k1),
"rank {my_rank} was handed the faces of layers {:?} but the decomposition gives it {:?}",
faces.owned(),
(k0, k1),
);
if k0 == k1 {
return; }
let plane = cell_plane(faces.nx, faces.ny);
let mut below = Vec::new();
if k0 > 0 {
let src = slab_owner(bounds, k0 - 1)
.expect("the layer below an owned layer is owned by some rank");
below = vec![Fix128::ZERO; plane];
channel.recv_plane(src, k0 - 1, &mut below);
}
let below_w = if k0 > 0 { Some(&below[..]) } else { None };
enforce_slab_face_boundaries(faces, cond, below_w);
if let Some(dst) = slab_owner(bounds, k1) {
let (w, _) = faces.w_layer(k1 - 1);
channel.send_plane(dst, k1 - 1, w);
}
}
pub(crate) fn enforce_slab_face_boundaries_over(
faces: &mut [SlabFaces],
cond: &[SlabFaceConditions],
) {
let Some(first) = faces.first() else {
return;
};
let plane = cell_plane(first.nx, first.ny);
let ranks = faces.len();
assert_eq!(
cond.len(),
ranks,
"{} sets of conditions for {ranks} ranks",
cond.len(),
);
let mut channel = LocalPlaneChannel::new(plane);
for (r, rank_faces) in faces.iter_mut().enumerate() {
enforce_slab_face_boundaries_on_rank(rank_faces, &cond[r], ranks, r, &mut channel);
}
assert!(
channel.is_drained(),
"a Z-face layer was passed on to a rank that never read it",
);
}
#[cfg(feature = "std")]
fn encode_plane(values: &[Fix128], wire: &mut [u8]) {
for (cell, chunk) in values
.iter()
.zip(wire.chunks_exact_mut(WIRE_BYTES_PER_CELL))
{
chunk[..8].copy_from_slice(&cell.hi.to_le_bytes());
chunk[8..].copy_from_slice(&cell.lo.to_le_bytes());
}
}
#[cfg(feature = "std")]
fn decode_plane(wire: &[u8], values: &mut [Fix128]) {
for (cell, chunk) in values
.iter_mut()
.zip(wire.chunks_exact(WIRE_BYTES_PER_CELL))
{
let hi = i64::from_le_bytes(chunk[..8].try_into().expect("8 bytes of a 16-byte cell"));
let lo = u64::from_le_bytes(chunk[8..].try_into().expect("8 bytes of a 16-byte cell"));
*cell = Fix128::from_raw(hi, lo);
}
}
#[cfg(feature = "std")]
pub(crate) struct SlabSocketTransport<S> {
my_rank: usize,
plane: usize,
slab: SlabStorage,
links: Vec<Option<S>>,
wire: Vec<u8>,
}
#[cfg(feature = "std")]
impl<S> SlabSocketTransport<S> {
pub(crate) fn new(
my_rank: usize,
plane: usize,
slab: SlabStorage,
links: Vec<Option<S>>,
) -> Self {
Self {
my_rank,
plane,
slab,
links,
wire: vec![0u8; plane * WIRE_BYTES_PER_CELL],
}
}
pub(crate) fn slab(&self) -> &SlabStorage {
&self.slab
}
}
#[cfg(feature = "std")]
impl<S: std::io::Read + std::io::Write> SlabTransport for SlabSocketTransport<S> {
fn slab_mut(&mut self, rank: usize) -> &mut SlabStorage {
assert_eq!(
rank, self.my_rank,
"rank {} was asked for rank {rank}'s band, which is in another address space: \
the driver is not rank-local",
self.my_rank,
);
&mut self.slab
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
if src == dst {
return;
}
if self.my_rank == src {
let from = self.slab.layer(layer).expect(DELIVERY_SOURCE_LACKS_LAYER);
encode_plane(from, &mut self.wire);
let link = self.links[dst]
.as_mut()
.expect("no stream to the rank this delivery is addressed to");
link.write_all(&self.wire)
.expect("send a halo layer to the peer rank");
link.flush().expect("flush a halo layer to the peer rank");
} else if self.my_rank == dst {
let link = self.links[src]
.as_mut()
.expect("no stream from the rank this delivery comes from");
link.read_exact(&mut self.wire)
.expect("receive a halo layer from the peer rank");
let into = self
.slab
.layer_mut(layer)
.expect(DELIVERY_DESTINATION_LACKS_LAYER);
decode_plane(&self.wire, into);
}
}
}
#[cfg(feature = "std")]
impl<S: std::io::Read + std::io::Write> PlaneChannel for SlabSocketTransport<S> {
fn send_plane(&mut self, dst: usize, layer: usize, plane: &[Fix128]) {
assert_eq!(
plane.len(),
self.plane,
"layer {layer} is {} values and a plane on this wire is {}",
plane.len(),
self.plane,
);
encode_plane(plane, &mut self.wire);
let link = self.links[dst]
.as_mut()
.expect("no stream to the rank this plane is addressed to");
link.write_all(&self.wire)
.expect("send a Z-face layer to the peer rank");
link.flush().expect("flush a Z-face layer to the peer rank");
}
fn recv_plane(&mut self, src: usize, layer: usize, plane: &mut [Fix128]) {
assert_eq!(
plane.len(),
self.plane,
"layer {layer} is {} values and a plane on this wire is {}",
plane.len(),
self.plane,
);
let link = self.links[src]
.as_mut()
.expect("no stream from the rank this plane comes from");
link.read_exact(&mut self.wire)
.expect("receive a Z-face layer from the peer rank");
decode_plane(&self.wire, plane);
}
}
#[allow(clippy::too_many_arguments)]
pub(crate) fn project_pressure_slab_local_on_rank<T: SlabTransport>(
faces: &mut SlabFaces,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
ranks: usize,
schedule: HaloSchedule,
my_rank: usize,
transport: &mut T,
) -> SlabBytes {
if faces.dx.is_zero()
|| density_kg_m3.is_zero()
|| dt_s.is_zero()
|| ranks == 0
|| my_rank >= ranks
{
return SlabBytes::ZERO;
}
let nz = faces.nz;
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
assert_eq!(
faces.owned(),
bounds[my_rank],
"rank {my_rank} was handed the faces of layers {:?} but the decomposition gives it {:?}",
faces.owned(),
bounds[my_rank],
);
let scale = density_kg_m3 * faces.dx * faces.dx / dt_s;
let stencil = SlabStencil::build(faces, scale);
for _ in 0..iterations {
for colour in 0..2u32 {
sweep_slab_colour(transport.slab_mut(my_rank), faces, &stencil, colour);
if schedule == HaloSchedule::EverySweep {
exchange_slab_halos_local(transport, &bounds, nz);
}
}
if schedule == HaloSchedule::EveryIteration {
exchange_slab_halos_local(transport, &bounds, nz);
}
}
let inv_dx = Fix128::ONE / faces.dx;
let coeff = dt_s / density_kg_m3 * inv_dx;
subtract_slab_pressure_gradient(faces, transport.slab_mut(my_rank), coeff);
stencil.bytes()
}
pub(crate) fn project_pressure_jacobi(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
iterations: u32,
) {
if grid.dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() {
return;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let n = grid.nx * grid.ny * grid.nz;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
let mut p_new = grid.pressure.clone();
for _ in 0..iterations {
for k in 0..grid.nz {
for j in 0..grid.ny {
for i in 0..grid.nx {
let idx = i + grid.nx * (j + grid.ny * k);
let neighbours = mask.neighbour_sum(&grid.pressure, i, j, k);
p_new[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
grid.pressure.clone_from(&p_new);
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
}
pub(crate) fn project_pressure_bicgstab(
grid: &mut MacGrid,
dt_s: Fix128,
density_kg_m3: Fix128,
max_iterations: u32,
tolerance: Fix128,
) -> BicgstabStats {
let default_stats = BicgstabStats {
iterations: 0,
final_residual: Fix128::ZERO,
converged: true,
};
if grid.dx.is_zero() || density_kg_m3.is_zero() || dt_s.is_zero() {
return default_stats;
}
grid.enforce_face_boundaries();
let scale = density_kg_m3 * grid.dx * grid.dx / dt_s;
let n = grid.nx * grid.ny * grid.nz;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let mut diag = vec![Fix128::ZERO; n];
for (c, slot) in diag.iter_mut().enumerate() {
*slot = Fix128::from_int(-mask.degree(c)); }
let mut x = grid.pressure.clone();
let mut r = vec![Fix128::ZERO; n];
mask.apply(&x, &mut r);
for i in 0..n {
r[i] = rhs[i] - r[i];
}
let r_hat = r.clone();
let mut p_vec = r.clone();
let mut rho_prev = dot(&r_hat, &r);
let mut alpha = Fix128::ONE;
let mut omega = Fix128::ONE;
let mut v_vec = vec![Fix128::ZERO; n];
let mut y_vec = vec![Fix128::ZERO; n];
let mut z_vec = vec![Fix128::ZERO; n];
let mut s_vec = vec![Fix128::ZERO; n];
let mut t_vec = vec![Fix128::ZERO; n];
let mut iterations = 0_u32;
let mut residual = linf_norm(&r);
let mut converged = residual < tolerance;
while iterations < max_iterations && !converged {
iterations += 1;
let rho = dot(&r_hat, &r);
if rho.is_zero() {
break;
}
if iterations > 1 {
let beta = (rho / rho_prev) * (alpha / omega);
for i in 0..n {
p_vec[i] = r[i] + beta * (p_vec[i] - omega * v_vec[i]);
}
}
for i in 0..n {
y_vec[i] = if diag[i].is_zero() {
p_vec[i]
} else {
p_vec[i] / diag[i]
};
}
mask.apply(&y_vec, &mut v_vec);
let denom = dot(&r_hat, &v_vec);
if denom.is_zero() {
break;
}
alpha = rho / denom;
for i in 0..n {
s_vec[i] = r[i] - alpha * v_vec[i];
}
let s_norm = linf_norm(&s_vec);
if s_norm < tolerance {
for i in 0..n {
x[i] = x[i] + alpha * y_vec[i];
}
residual = s_norm;
converged = true;
break;
}
for i in 0..n {
z_vec[i] = if diag[i].is_zero() {
s_vec[i]
} else {
s_vec[i] / diag[i]
};
}
mask.apply(&z_vec, &mut t_vec);
let tt = dot(&t_vec, &t_vec);
if tt.is_zero() {
for i in 0..n {
x[i] = x[i] + alpha * y_vec[i];
}
residual = s_norm;
converged = residual < tolerance;
break;
}
omega = dot(&t_vec, &s_vec) / tt;
for i in 0..n {
x[i] = x[i] + alpha * y_vec[i] + omega * z_vec[i];
}
for i in 0..n {
r[i] = s_vec[i] - omega * t_vec[i];
}
residual = linf_norm(&r);
converged = residual < tolerance;
rho_prev = rho;
}
grid.pressure = x;
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density_kg_m3 * inv_dx);
BicgstabStats {
iterations,
final_residual: residual,
converged,
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub struct BicgstabStats {
pub iterations: u32,
pub final_residual: Fix128,
pub converged: bool,
}
fn dot(a: &[Fix128], b: &[Fix128]) -> Fix128 {
let mut acc = Fix128::ZERO;
for i in 0..a.len() {
acc = acc + a[i] * b[i];
}
acc
}
fn linf_norm(v: &[Fix128]) -> Fix128 {
let mut best = Fix128::ZERO;
for &x in v {
let ax = x.abs();
if ax > best {
best = ax;
}
}
best
}
fn split(p_over_dx: Fix128, axis_offset: Fix128) -> (usize, Fix128) {
let shifted = p_over_dx - axis_offset;
if shifted.hi < 0 {
(0, Fix128::ZERO)
} else {
(
shifted.hi as usize,
Fix128 {
hi: 0,
lo: shifted.lo,
},
)
}
}
pub(crate) fn sample_u_trilinear(grid: &MacGrid, pos_m: Vec3Fix) -> Fix128 {
if grid.dx.is_zero() {
return Fix128::ZERO;
}
let inv_dx = Fix128::ONE / grid.dx;
let fx = pos_m.x * inv_dx;
let fy = pos_m.y * inv_dx;
let fz = pos_m.z * inv_dx;
let (i, u) = split(fx, Fix128::ZERO);
let (j, v) = split(fy, Fix128::from_ratio(1, 2));
let (k, w) = split(fz, Fix128::from_ratio(1, 2));
let i0 = i.min(grid.nx);
let i1 = (i + 1).min(grid.nx);
let j0 = j.min(grid.ny - 1);
let j1 = (j + 1).min(grid.ny - 1);
let k0 = k.min(grid.nz - 1);
let k1 = (k + 1).min(grid.nz - 1);
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let c00 = grid.u(i0, j0, k0) * om_u + grid.u(i1, j0, k0) * u;
let c10 = grid.u(i0, j1, k0) * om_u + grid.u(i1, j1, k0) * u;
let c01 = grid.u(i0, j0, k1) * om_u + grid.u(i1, j0, k1) * u;
let c11 = grid.u(i0, j1, k1) * om_u + grid.u(i1, j1, k1) * u;
let c0 = c00 * om_v + c10 * v;
let c1 = c01 * om_v + c11 * v;
c0 * om_w + c1 * w
}
pub(crate) fn sample_v_trilinear(grid: &MacGrid, pos_m: Vec3Fix) -> Fix128 {
if grid.dx.is_zero() {
return Fix128::ZERO;
}
let inv_dx = Fix128::ONE / grid.dx;
let fx = pos_m.x * inv_dx;
let fy = pos_m.y * inv_dx;
let fz = pos_m.z * inv_dx;
let (i, u) = split(fx, Fix128::from_ratio(1, 2));
let (j, v) = split(fy, Fix128::ZERO);
let (k, w) = split(fz, Fix128::from_ratio(1, 2));
let i0 = i.min(grid.nx - 1);
let i1 = (i + 1).min(grid.nx - 1);
let j0 = j.min(grid.ny);
let j1 = (j + 1).min(grid.ny);
let k0 = k.min(grid.nz - 1);
let k1 = (k + 1).min(grid.nz - 1);
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let c00 = grid.v(i0, j0, k0) * om_v + grid.v(i0, j1, k0) * v;
let c10 = grid.v(i1, j0, k0) * om_v + grid.v(i1, j1, k0) * v;
let c01 = grid.v(i0, j0, k1) * om_v + grid.v(i0, j1, k1) * v;
let c11 = grid.v(i1, j0, k1) * om_v + grid.v(i1, j1, k1) * v;
let c0 = c00 * om_u + c10 * u;
let c1 = c01 * om_u + c11 * u;
c0 * om_w + c1 * w
}
pub(crate) fn sample_w_trilinear(grid: &MacGrid, pos_m: Vec3Fix) -> Fix128 {
if grid.dx.is_zero() {
return Fix128::ZERO;
}
let inv_dx = Fix128::ONE / grid.dx;
let fx = pos_m.x * inv_dx;
let fy = pos_m.y * inv_dx;
let fz = pos_m.z * inv_dx;
let (i, u) = split(fx, Fix128::from_ratio(1, 2));
let (j, v) = split(fy, Fix128::from_ratio(1, 2));
let (k, w) = split(fz, Fix128::ZERO);
let i0 = i.min(grid.nx - 1);
let i1 = (i + 1).min(grid.nx - 1);
let j0 = j.min(grid.ny - 1);
let j1 = (j + 1).min(grid.ny - 1);
let k0 = k.min(grid.nz);
let k1 = (k + 1).min(grid.nz);
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let c00 = grid.w(i0, j0, k0) * om_w + grid.w(i0, j0, k1) * w;
let c10 = grid.w(i1, j0, k0) * om_w + grid.w(i1, j0, k1) * w;
let c01 = grid.w(i0, j1, k0) * om_w + grid.w(i0, j1, k1) * w;
let c11 = grid.w(i1, j1, k0) * om_w + grid.w(i1, j1, k1) * w;
let c0 = c00 * om_u + c10 * u;
let c1 = c01 * om_u + c11 * u;
c0 * om_v + c1 * v
}
pub(crate) fn sample_u_range(grid: &MacGrid, pos_m: Vec3Fix) -> (Fix128, Fix128) {
if grid.dx.is_zero() {
return (Fix128::ZERO, Fix128::ZERO);
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, _) = split(pos_m.x * inv_dx, Fix128::ZERO);
let (j, _) = split(pos_m.y * inv_dx, Fix128::from_ratio(1, 2));
let (k, _) = split(pos_m.z * inv_dx, Fix128::from_ratio(1, 2));
let i0 = i.min(grid.nx);
let i1 = (i + 1).min(grid.nx);
let j0 = j.min(grid.ny - 1);
let j1 = (j + 1).min(grid.ny - 1);
let k0 = k.min(grid.nz - 1);
let k1 = (k + 1).min(grid.nz - 1);
corner_range(&[
grid.u(i0, j0, k0),
grid.u(i1, j0, k0),
grid.u(i0, j1, k0),
grid.u(i1, j1, k0),
grid.u(i0, j0, k1),
grid.u(i1, j0, k1),
grid.u(i0, j1, k1),
grid.u(i1, j1, k1),
])
}
pub(crate) fn sample_v_range(grid: &MacGrid, pos_m: Vec3Fix) -> (Fix128, Fix128) {
if grid.dx.is_zero() {
return (Fix128::ZERO, Fix128::ZERO);
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, _) = split(pos_m.x * inv_dx, Fix128::from_ratio(1, 2));
let (j, _) = split(pos_m.y * inv_dx, Fix128::ZERO);
let (k, _) = split(pos_m.z * inv_dx, Fix128::from_ratio(1, 2));
let i0 = i.min(grid.nx - 1);
let i1 = (i + 1).min(grid.nx - 1);
let j0 = j.min(grid.ny);
let j1 = (j + 1).min(grid.ny);
let k0 = k.min(grid.nz - 1);
let k1 = (k + 1).min(grid.nz - 1);
corner_range(&[
grid.v(i0, j0, k0),
grid.v(i1, j0, k0),
grid.v(i0, j1, k0),
grid.v(i1, j1, k0),
grid.v(i0, j0, k1),
grid.v(i1, j0, k1),
grid.v(i0, j1, k1),
grid.v(i1, j1, k1),
])
}
pub(crate) fn sample_w_range(grid: &MacGrid, pos_m: Vec3Fix) -> (Fix128, Fix128) {
if grid.dx.is_zero() {
return (Fix128::ZERO, Fix128::ZERO);
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, _) = split(pos_m.x * inv_dx, Fix128::from_ratio(1, 2));
let (j, _) = split(pos_m.y * inv_dx, Fix128::from_ratio(1, 2));
let (k, _) = split(pos_m.z * inv_dx, Fix128::ZERO);
let i0 = i.min(grid.nx - 1);
let i1 = (i + 1).min(grid.nx - 1);
let j0 = j.min(grid.ny - 1);
let j1 = (j + 1).min(grid.ny - 1);
let k0 = k.min(grid.nz);
let k1 = (k + 1).min(grid.nz);
corner_range(&[
grid.w(i0, j0, k0),
grid.w(i1, j0, k0),
grid.w(i0, j1, k0),
grid.w(i1, j1, k0),
grid.w(i0, j0, k1),
grid.w(i1, j0, k1),
grid.w(i0, j1, k1),
grid.w(i1, j1, k1),
])
}
fn corner_range(corners: &[Fix128; 8]) -> (Fix128, Fix128) {
let mut lo = corners[0];
let mut hi = corners[0];
for &c in &corners[1..] {
if c < lo {
lo = c;
}
if c > hi {
hi = c;
}
}
(lo, hi)
}
#[must_use]
pub fn g2p_velocity(grid: &MacGrid, pos_m: Vec3Fix) -> Vec3Fix {
if grid.dx.is_zero() {
return Vec3Fix::default();
}
Vec3Fix::new(
sample_u_trilinear(grid, pos_m),
sample_v_trilinear(grid, pos_m),
sample_w_trilinear(grid, pos_m),
)
}
fn deposit_u_trilinear<S: FaceSink>(sink: &mut S, pos_m: Vec3Fix, vx: Fix128) {
let grid = sink.shape();
if grid.dx.is_zero() {
return;
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, u) = split(pos_m.x * inv_dx, Fix128::ZERO);
let (j, v) = split(pos_m.y * inv_dx, Fix128::from_ratio(1, 2));
let (k, w) = split(pos_m.z * inv_dx, Fix128::from_ratio(1, 2));
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let corners = [
(i, j, k, om_u * om_v * om_w),
(i + 1, j, k, u * om_v * om_w),
(i, j + 1, k, om_u * v * om_w),
(i + 1, j + 1, k, u * v * om_w),
(i, j, k + 1, om_u * om_v * w),
(i + 1, j, k + 1, u * om_v * w),
(i, j + 1, k + 1, om_u * v * w),
(i + 1, j + 1, k + 1, u * v * w),
];
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
for (ci, cj, ck, weight) in corners {
if ci <= nx && cj < ny && ck < nz {
let ix = sink.shape().idx_u(ci, cj, ck);
sink.add_u(ix, weight, vx);
}
}
}
fn deposit_v_trilinear<S: FaceSink>(sink: &mut S, pos_m: Vec3Fix, vy: Fix128) {
let grid = sink.shape();
if grid.dx.is_zero() {
return;
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, u) = split(pos_m.x * inv_dx, Fix128::from_ratio(1, 2));
let (j, v) = split(pos_m.y * inv_dx, Fix128::ZERO);
let (k, w) = split(pos_m.z * inv_dx, Fix128::from_ratio(1, 2));
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let corners = [
(i, j, k, om_u * om_v * om_w),
(i + 1, j, k, u * om_v * om_w),
(i, j + 1, k, om_u * v * om_w),
(i + 1, j + 1, k, u * v * om_w),
(i, j, k + 1, om_u * om_v * w),
(i + 1, j, k + 1, u * om_v * w),
(i, j + 1, k + 1, om_u * v * w),
(i + 1, j + 1, k + 1, u * v * w),
];
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
for (ci, cj, ck, weight) in corners {
if ci < nx && cj <= ny && ck < nz {
let ix = sink.shape().idx_v(ci, cj, ck);
sink.add_v(ix, weight, vy);
}
}
}
fn deposit_w_trilinear<S: FaceSink>(sink: &mut S, pos_m: Vec3Fix, vz: Fix128) {
let grid = sink.shape();
if grid.dx.is_zero() {
return;
}
let inv_dx = Fix128::ONE / grid.dx;
let (i, u) = split(pos_m.x * inv_dx, Fix128::from_ratio(1, 2));
let (j, v) = split(pos_m.y * inv_dx, Fix128::from_ratio(1, 2));
let (k, w) = split(pos_m.z * inv_dx, Fix128::ZERO);
let om_u = Fix128::ONE - u;
let om_v = Fix128::ONE - v;
let om_w = Fix128::ONE - w;
let corners = [
(i, j, k, om_u * om_v * om_w),
(i + 1, j, k, u * om_v * om_w),
(i, j + 1, k, om_u * v * om_w),
(i + 1, j + 1, k, u * v * om_w),
(i, j, k + 1, om_u * om_v * w),
(i + 1, j, k + 1, u * om_v * w),
(i, j + 1, k + 1, om_u * v * w),
(i + 1, j + 1, k + 1, u * v * w),
];
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
for (ci, cj, ck, weight) in corners {
if ci < nx && cj < ny && ck <= nz {
let ix = sink.shape().idx_w(ci, cj, ck);
sink.add_w(ix, weight, vz);
}
}
}
pub(crate) fn p2g_trilinear<S: FaceSink>(sink: &mut S, pos_m: Vec3Fix, vel_m_per_s: Vec3Fix) {
deposit_u_trilinear(sink, pos_m, vel_m_per_s.x);
deposit_v_trilinear(sink, pos_m, vel_m_per_s.y);
deposit_w_trilinear(sink, pos_m, vel_m_per_s.z);
}
pub(crate) fn p2g_nearest<S: FaceSink>(sink: &mut S, pos_m: Vec3Fix, vel_m_per_s: Vec3Fix) {
let shape = sink.shape();
if shape.dx.is_zero() {
return;
}
let inv_dx = Fix128::ONE / shape.dx;
let (Some(ix), Some(iy), Some(iz)) = (
nearest_cell(pos_m.x, inv_dx, shape.nx),
nearest_cell(pos_m.y, inv_dx, shape.ny),
nearest_cell(pos_m.z, inv_dx, shape.nz),
) else {
return;
};
let half = Fix128::from_ratio(1, 2);
let (u_lo, u_hi) = (shape.idx_u(ix, iy, iz), shape.idx_u(ix + 1, iy, iz));
let (v_lo, v_hi) = (shape.idx_v(ix, iy, iz), shape.idx_v(ix, iy + 1, iz));
let (w_lo, w_hi) = (shape.idx_w(ix, iy, iz), shape.idx_w(ix, iy, iz + 1));
sink.add_u(u_lo, half, vel_m_per_s.x);
sink.add_u(u_hi, half, vel_m_per_s.x);
sink.add_v(v_lo, half, vel_m_per_s.y);
sink.add_v(v_hi, half, vel_m_per_s.y);
sink.add_w(w_lo, half, vel_m_per_s.z);
sink.add_w(w_hi, half, vel_m_per_s.z);
}
fn nearest_cell(p_m: Fix128, inv_dx: Fix128, n: usize) -> Option<usize> {
if p_m < Fix128::ZERO || n == 0 {
return None;
}
let c = p_m.checked_mul(inv_dx)?;
if c < Fix128::ZERO {
return None;
}
let i = usize::try_from(c.hi).ok()?;
if i < n {
Some(i)
} else if i == n && c.lo == 0 {
Some(n - 1)
} else {
None
}
}
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
#[non_exhaustive]
pub enum ParticleScatter {
#[default]
Trilinear,
Nearest,
}
pub(crate) trait FaceSink {
fn shape(&self) -> &MacGrid;
fn add_u(&mut self, ix: usize, weight: Fix128, v: Fix128);
fn add_v(&mut self, ix: usize, weight: Fix128, v: Fix128);
fn add_w(&mut self, ix: usize, weight: Fix128, v: Fix128);
}
impl FaceSink for MacGrid {
#[inline]
fn shape(&self) -> &MacGrid {
self
}
#[inline]
fn add_u(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.u[ix] = self.u[ix] + weight * v;
}
#[inline]
fn add_v(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.v[ix] = self.v[ix] + weight * v;
}
#[inline]
fn add_w(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.w[ix] = self.w[ix] + weight * v;
}
}
#[inline]
fn fix_raw(x: Fix128) -> i128 {
(i128::from(x.hi) << 64) | i128::from(x.lo)
}
#[inline]
fn fix_from_raw(r: i128) -> Fix128 {
Fix128 {
hi: (r >> 64) as i64,
lo: r as u64,
}
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct I256 {
hi: i128,
lo: u128,
}
impl I256 {
const ZERO: Self = Self { hi: 0, lo: 0 };
#[inline]
fn wrapping_add(self, o: Self) -> Self {
let (lo, carry) = self.lo.overflowing_add(o.lo);
Self {
hi: self.hi.wrapping_add(o.hi).wrapping_add(i128::from(carry)),
lo,
}
}
#[inline]
fn wrapping_neg(self) -> Self {
let lo = (!self.lo).wrapping_add(1);
Self {
hi: (!self.hi).wrapping_add(i128::from(lo == 0)),
lo,
}
}
fn from_product(a: i128, b: i128) -> Self {
let (hi, lo) = umul_128(a.unsigned_abs(), b.unsigned_abs());
let mag = Self { hi: hi as i128, lo };
if (a < 0) != (b < 0) {
mag.wrapping_neg()
} else {
mag
}
}
}
fn umul_128(a: u128, b: u128) -> (u128, u128) {
const M: u128 = u64::MAX as u128;
let (a0, a1) = (a & M, a >> 64);
let (b0, b1) = (b & M, b >> 64);
let ll = a0 * b0;
let lh = a0 * b1;
let hl = a1 * b0;
let hh = a1 * b1;
let mid = (ll >> 64) + (lh & M) + (hl & M);
let lo = ((mid & M) << 64) | (ll & M);
let hi = hh + (lh >> 64) + (hl >> 64) + (mid >> 64);
(hi, lo)
}
fn udiv_256_by_128(n_hi: u128, n_lo: u128, d: u128) -> (u128, u128, u128) {
debug_assert!(d != 0);
let top = if n_hi != 0 {
255 - n_hi.leading_zeros()
} else if n_lo != 0 {
127 - n_lo.leading_zeros()
} else {
return (0, 0, 0);
};
let (mut q_hi, mut q_lo, mut rem) = (0u128, 0u128, 0u128);
for i in (0..=top).rev() {
let bit = if i >= 128 {
(n_hi >> (i - 128)) & 1
} else {
(n_lo >> i) & 1
};
let carry = rem >> 127;
rem = (rem << 1) | bit;
if carry == 1 || rem >= d {
rem = rem.wrapping_sub(d);
if i >= 128 {
q_hi |= 1u128 << (i - 128);
} else {
q_lo |= 1u128 << i;
}
}
}
(q_hi, q_lo, rem)
}
fn div_trunc_i256_by_u128(n: I256, d: u128) -> i128 {
let negative = n.hi < 0;
let mag = if negative { n.wrapping_neg() } else { n };
let (_, q, _) = udiv_256_by_128(mag.hi as u128, mag.lo, d);
if negative {
(q as i128).wrapping_neg()
} else {
q as i128
}
}
#[derive(Clone, Copy)]
struct Moments {
p: I256,
d: u128,
}
impl Moments {
const ZERO: Self = Self {
p: I256::ZERO,
d: 0,
};
#[inline]
fn add(&mut self, weight: Fix128, v: Fix128) {
let w = fix_raw(weight);
self.p = self.p.wrapping_add(I256::from_product(w, fix_raw(v)));
self.d = self.d.wrapping_add(w as u128);
}
#[inline]
fn mean(&self) -> Option<Fix128> {
if self.d == 0 {
None
} else {
Some(fix_from_raw(div_trunc_i256_by_u128(self.p, self.d)))
}
}
}
struct WideFaces {
shape: MacGrid,
u: Vec<Moments>,
v: Vec<Moments>,
w: Vec<Moments>,
}
impl WideFaces {
fn new(grid: &MacGrid) -> Self {
Self {
shape: MacGrid::new(grid.nx, grid.ny, grid.nz, grid.dx),
u: vec![Moments::ZERO; grid.u.len()],
v: vec![Moments::ZERO; grid.v.len()],
w: vec![Moments::ZERO; grid.w.len()],
}
}
}
impl FaceSink for WideFaces {
#[inline]
fn shape(&self) -> &MacGrid {
&self.shape
}
#[inline]
fn add_u(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.u[ix].add(weight, v);
}
#[inline]
fn add_v(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.v[ix].add(weight, v);
}
#[inline]
fn add_w(&mut self, ix: usize, weight: Fix128, v: Fix128) {
self.w[ix].add(weight, v);
}
}
pub fn p2g_normalized(grid: &mut MacGrid, particles: &[(Vec3Fix, Vec3Fix)]) {
p2g_normalized_with(grid, particles, ParticleScatter::Trilinear);
}
pub fn p2g_normalized_with(
grid: &mut MacGrid,
particles: &[(Vec3Fix, Vec3Fix)],
scatter: ParticleScatter,
) {
if grid.dx.is_zero() {
return;
}
let mut acc = WideFaces::new(grid);
for &(pos, vel) in particles {
if pos.x < Fix128::ZERO || pos.y < Fix128::ZERO || pos.z < Fix128::ZERO {
continue;
}
match scatter {
ParticleScatter::Trilinear => p2g_trilinear(&mut acc, pos, vel),
ParticleScatter::Nearest => p2g_nearest(&mut acc, pos, vel),
}
}
let resolve = |dst: &mut [Fix128], src: &[Moments]| {
for (out, m) in dst.iter_mut().zip(src) {
if let Some(mean) = m.mean() {
*out = mean;
}
}
};
resolve(&mut grid.u, &acc.u);
resolve(&mut grid.v, &acc.v);
resolve(&mut grid.w, &acc.w);
}
#[cfg(test)]
#[allow(clippy::disallowed_methods)]
mod tests {
use super::*;
#[cfg(not(feature = "std"))]
use alloc::format;
#[test]
fn mac_grid_dimensions() {
let g = MacGrid::new(4, 3, 2, Fix128::ONE);
assert_eq!(g.u.len(), 5 * 3 * 2);
assert_eq!(g.v.len(), 4 * 4 * 2);
assert_eq!(g.w.len(), 4 * 3 * 3);
assert_eq!(g.pressure.len(), 24);
}
#[test]
fn mac_grid_zero_initialised() {
let g = MacGrid::new(3, 3, 3, Fix128::ONE);
assert_eq!(g.u(0, 0, 0), Fix128::ZERO);
assert_eq!(g.v(1, 1, 1), Fix128::ZERO);
assert_eq!(g.pressure(0, 0, 0), Fix128::ZERO);
}
#[test]
fn mac_grid_out_of_range_returns_zero() {
let g = MacGrid::new(3, 3, 3, Fix128::ONE);
assert_eq!(g.u(10, 0, 0), Fix128::ZERO);
assert_eq!(g.pressure(5, 0, 0), Fix128::ZERO);
}
#[test]
fn cell_velocity_averages_faces() {
let mut g = MacGrid::new(2, 1, 1, Fix128::ONE);
let i0 = g.idx_u(0, 0, 0);
let i1 = g.idx_u(1, 0, 0);
g.u[i0] = Fix128::from_int(2);
g.u[i1] = Fix128::from_int(4);
let (uc, _, _) = g.cell_velocity(0, 0, 0);
assert_eq!(uc, Fix128::from_int(3));
}
#[test]
fn divergence_zero_for_still_grid() {
let g = MacGrid::new(3, 3, 3, Fix128::ONE);
assert_eq!(g.divergence(1, 1, 1), Fix128::ZERO);
}
#[test]
fn divergence_positive_for_expanding_flow() {
let mut g = MacGrid::new(3, 3, 3, Fix128::ONE);
let ix = g.idx_u(2, 1, 1);
g.u[ix] = Fix128::ONE;
let d = g.divergence(1, 1, 1);
assert!(d > Fix128::ZERO);
}
#[test]
fn project_zero_flow_stays_zero() {
let mut g = MacGrid::new(3, 3, 3, Fix128::ONE);
project_pressure(
&mut g,
Fix128::from_ratio(1, 60),
Fix128::from_int(1000),
10,
);
assert_eq!(g.divergence(1, 1, 1), Fix128::ZERO);
}
#[test]
fn project_reduces_divergence() {
let mut g = MacGrid::new(4, 4, 4, Fix128::ONE);
for k in 0..4 {
for j in 0..4 {
for i in 0..=4 {
let ix = g.idx_u(i, j, k);
g.u[ix] = Fix128::from_int(i as i64);
}
}
}
let div_before = g.divergence(2, 2, 2);
project_pressure(
&mut g,
Fix128::from_ratio(1, 60),
Fix128::from_int(1000),
50,
);
let div_after = g.divergence(2, 2, 2);
assert!(div_after.abs() < div_before.abs());
}
#[test]
fn p2g_scatters_x_velocity() {
let mut g = MacGrid::new(4, 4, 4, Fix128::ONE);
let pos = Vec3Fix::new(
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
);
let vel = Vec3Fix::new(Fix128::from_int(4), Fix128::ZERO, Fix128::ZERO);
p2g_nearest(&mut g, pos, vel);
assert_eq!(g.u(1, 1, 1), Fix128::from_int(2));
assert_eq!(g.u(2, 1, 1), Fix128::from_int(2));
}
#[test]
fn g2p_reads_cell_average_velocity() {
let mut g = MacGrid::new(4, 4, 4, Fix128::ONE);
let i0 = g.idx_u(1, 1, 1);
let i1 = g.idx_u(2, 1, 1);
g.u[i0] = Fix128::from_int(2);
g.u[i1] = Fix128::from_int(4);
let vel = g2p_velocity(
&g,
Vec3Fix::new(
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
),
);
assert_eq!(vel.x, Fix128::from_int(3));
}
#[test]
fn p2g_out_of_range_ignored() {
let mut g = MacGrid::new(2, 2, 2, Fix128::ONE);
let pos = Vec3Fix::new(
Fix128::from_int(100),
Fix128::from_int(100),
Fix128::from_int(100),
);
let vel = Vec3Fix::new(Fix128::ONE, Fix128::ONE, Fix128::ONE);
p2g_nearest(&mut g, pos, vel);
assert_eq!(g.u(0, 0, 0), Fix128::ZERO);
}
mod wide_p2g {
use super::*;
struct Rng(u64);
impl Rng {
fn next(&mut self) -> u64 {
self.0 ^= self.0 << 13;
self.0 ^= self.0 >> 7;
self.0 ^= self.0 << 17;
self.0
}
fn u128(&mut self) -> u128 {
(u128::from(self.next()) << 64) | u128::from(self.next())
}
fn u128_any_len(&mut self) -> u128 {
let bits = (self.next() % 129) as u32;
if bits == 0 {
0
} else {
self.u128() >> (128 - bits)
}
}
}
const M128: u128 = u128::MAX;
fn ref_umul(a: u128, b: u128) -> (u128, u128) {
let la: Vec<u128> = (0..4).map(|i| (a >> (32 * i)) & 0xffff_ffff).collect();
let lb: Vec<u128> = (0..4).map(|i| (b >> (32 * i)) & 0xffff_ffff).collect();
let mut col = [0u128; 9];
for i in 0..4 {
for j in 0..4 {
col[i + j] += la[i] * lb[j];
}
}
let mut carry = 0u128;
let mut limbs = [0u128; 8];
for k in 0..8 {
let t = col[k] + carry;
limbs[k] = t & 0xffff_ffff;
carry = t >> 32;
}
assert_eq!(carry, 0);
let lo = (0..4).fold(0u128, |acc, k| acc | (limbs[k] << (32 * k)));
let hi = (0..4).fold(0u128, |acc, k| acc | (limbs[k + 4] << (32 * k)));
(hi, lo)
}
fn ref_udiv(n_hi: u128, n_lo: u128, d: u128) -> (u128, u128, u128) {
let mut bits: Vec<bool> = Vec::new();
for i in (0..128).rev() {
bits.push((n_hi >> i) & 1 == 1);
}
for i in (0..128).rev() {
bits.push((n_lo >> i) & 1 == 1);
}
let mut rem: u128 = 0;
let mut qbits: Vec<bool> = Vec::new();
for b in bits {
let top = rem >> 127 == 1;
rem = (rem << 1) | u128::from(b);
if top || rem >= d {
rem = rem.wrapping_sub(d);
qbits.push(true);
} else {
qbits.push(false);
}
}
let to = |s: &[bool]| s.iter().fold(0u128, |a, &b| (a << 1) | u128::from(b));
(to(&qbits[..128]), to(&qbits[128..]), rem)
}
fn mul_q_d(q_hi: u128, q_lo: u128, d: u128) -> (u128, u128) {
let (h0, l0) = umul_128(q_lo, d);
let (h1, l1) = umul_128(q_hi, d);
assert_eq!(h1, 0, "q·d overflows 256 bits");
(h0.checked_add(l1).expect("q·d overflows 256 bits"), l0)
}
fn add_256(a: (u128, u128), b: u128) -> (u128, u128) {
let (lo, c) = a.1.overflowing_add(b);
(a.0.checked_add(u128::from(c)).expect("overflow"), lo)
}
fn i256_from_i128(x: i128) -> I256 {
I256 {
hi: if x < 0 { -1 } else { 0 },
lo: x as u128,
}
}
fn neg_mag(h: u128, l: u128) -> I256 {
I256 {
hi: h as i128,
lo: l,
}
.wrapping_neg()
}
#[test]
fn umul_matches_the_limb_schoolbook_reference() {
let mut r = Rng(0x1357_9bdf_2468_ace0);
assert_eq!(umul_128(M128, M128), (M128 - 1, 1));
assert_eq!(umul_128(0, M128), (0, 0));
assert_eq!(umul_128(1, M128), (0, M128));
assert_eq!(umul_128(1 << 127, 2), (1, 0));
for _ in 0..20_000 {
let (a, b) = (r.u128_any_len(), r.u128_any_len());
assert_eq!(umul_128(a, b), ref_umul(a, b), "{a} * {b}");
}
}
#[test]
fn small_division_matches_the_builtin_u128_operators() {
let mut r = Rng(0x0bad_cafe_f00d_1234);
for _ in 0..20_000 {
let n = r.u128_any_len();
let d = r.u128_any_len().max(1);
assert_eq!(udiv_256_by_128(0, n, d), (0, n / d, n % d), "{n} / {d}");
}
assert_eq!(udiv_256_by_128(0, 0, 1), (0, 0, 0));
assert_eq!(udiv_256_by_128(0, 7, 1), (0, 7, 0));
assert_eq!(udiv_256_by_128(0, M128, 1), (0, M128, 0));
assert_eq!(udiv_256_by_128(0, M128, M128), (0, 1, 0));
assert_eq!(udiv_256_by_128(0, M128 - 1, M128), (0, 0, M128 - 1));
}
#[test]
fn wide_division_satisfies_the_division_identity_and_matches_the_reference() {
let mut r = Rng(0x5eed_5eed_5eed_5eed);
let mut big_divisor = 0;
for round in 0..30_000 {
let n_hi = match round % 4 {
0 => 0,
1 => r.u128_any_len(),
_ => r.u128(),
};
let n_lo = r.u128_any_len();
let d = if round % 3 == 0 {
(1u128 << 127) | r.u128_any_len()
} else {
r.u128_any_len().max(1)
};
if d >> 127 == 1 {
big_divisor += 1;
}
let (qh, ql, rem) = udiv_256_by_128(n_hi, n_lo, d);
assert!(rem < d, "remainder {rem} >= divisor {d}");
let prod = mul_q_d(qh, ql, d);
assert_eq!(add_256(prod, rem), (n_hi, n_lo), "q·d + r != n");
assert_eq!((qh, ql, rem), ref_udiv(n_hi, n_lo, d));
}
assert!(big_divisor > 5_000);
let (qh, ql, rem) = udiv_256_by_128(M128, M128, 1);
assert_eq!((qh, ql, rem), (M128, M128, 0));
let (qh, ql, rem) = udiv_256_by_128(M128, M128, M128);
assert_eq!((qh, ql, rem), (1, 1, 0));
}
#[test]
fn an_exactly_divisible_dividend_leaves_a_zero_remainder() {
let mut r = Rng(0xfeed_beef_0001_0002);
for _ in 0..5_000 {
let d = r.u128_any_len().max(1);
let q = r.u128_any_len();
let (h, l) = umul_128(q, d);
assert_eq!(udiv_256_by_128(h, l, d), (0, q, 0), "{q} * {d}");
}
}
#[test]
fn i256_product_sum_and_negation() {
let mut r = Rng(0x7777_1111_9999_3333);
for _ in 0..10_000 {
let sh = r.next() % 60;
let a = (r.u128() as i128) >> (65 + sh);
let b = (r.u128() as i128) >> (65 + sh);
assert_eq!(I256::from_product(a, b), i256_from_i128(a * b), "{a} * {b}");
}
assert_eq!(
I256::from_product(i128::MIN, i128::MIN),
I256 {
hi: 1 << 126,
lo: 0
}
);
let p = I256::from_product(i128::MIN, i128::MAX);
assert_eq!(p, neg_mag((1 << 126) - 1, 1 << 127));
let a = I256 { hi: 0, lo: M128 };
let one = I256 { hi: 0, lo: 1 };
assert_eq!(a.wrapping_add(one), I256 { hi: 1, lo: 0 });
assert_eq!(I256 { hi: 1, lo: 0 }.wrapping_neg(), I256 { hi: -1, lo: 0 });
assert_eq!(I256::ZERO.wrapping_neg(), I256::ZERO);
let min = I256 {
hi: i128::MIN,
lo: 0,
};
assert_eq!(min.wrapping_neg(), min);
let x = I256 {
hi: 12345,
lo: 6789,
};
assert_eq!(x.wrapping_add(x.wrapping_neg()), I256::ZERO);
}
fn build(q: u128, d: u128, r: u128, neg: bool) -> I256 {
let (h, l) = umul_128(q, d);
let (h, l) = add_256((h, l), r);
if neg {
neg_mag(h, l)
} else {
I256 {
hi: h as i128,
lo: l,
}
}
}
#[test]
fn signed_division_truncates_toward_zero_at_the_i128_boundaries() {
let max = i128::MAX as u128; let min_mag = 1u128 << 127; for &d in &[
1u128,
2,
3,
1 << 64,
(1 << 64) + 7,
1 << 127,
(1 << 127) + 5,
M128,
] {
for r in [0, d - 1, d / 2] {
assert_eq!(
div_trunc_i256_by_u128(build(max, d, r, false), d),
i128::MAX,
"+max d={d} r={r}"
);
assert_eq!(
div_trunc_i256_by_u128(build(max, d, r, true), d),
-i128::MAX,
"-max d={d} r={r}"
);
assert_eq!(
div_trunc_i256_by_u128(build(min_mag, d, r, true), d),
i128::MIN,
"-2^127 d={d} r={r}"
);
}
}
assert_eq!(div_trunc_i256_by_u128(i256_from_i128(5), 7), 0);
assert_eq!(div_trunc_i256_by_u128(i256_from_i128(-5), 7), 0);
assert_eq!(div_trunc_i256_by_u128(I256::ZERO, 1), 0);
let mut r = Rng(0x2468_ace0_1357_9bdf);
for _ in 0..20_000 {
let n = (r.u128() as i128) >> (r.next() % 120);
let d = r.u128_any_len().max(1);
if d > i128::MAX as u128 {
continue;
}
assert_eq!(
div_trunc_i256_by_u128(i256_from_i128(n), d),
n / d as i128,
"{n} / {d}"
);
}
assert_eq!(
div_trunc_i256_by_u128(i256_from_i128(i128::MAX), 1),
i128::MAX
);
assert_eq!(
div_trunc_i256_by_u128(i256_from_i128(i128::MIN), 1),
i128::MIN
);
}
#[test]
fn a_numerator_beyond_2_pow_191_divides_back_to_the_velocity() {
for &v in &[
i128::MAX,
i128::MIN,
i128::MIN + 1,
1i128 << 100,
-(1i128 << 126) - 3,
] {
for n in 1u64..=20 {
let w = 1u128 << 64;
let mut p = I256::ZERO;
for _ in 0..n {
p = p.wrapping_add(I256::from_product(w as i128, v));
}
let d = u128::from(n) * w;
assert_eq!(div_trunc_i256_by_u128(p, d), v, "v={v} n={n}");
}
}
}
#[test]
fn the_most_negative_numerator_is_handled_as_a_magnitude_of_2_pow_255() {
let min = I256 {
hi: i128::MIN,
lo: 0,
};
assert_eq!(div_trunc_i256_by_u128(min, M128), i128::MIN);
assert_eq!(
udiv_256_by_128(1u128 << 127, 0, M128),
(0, 1u128 << 127, 1u128 << 127)
);
}
}
#[test]
fn p2g_trilinear_deposits_to_8_corners() {
let mut g = MacGrid::new(4, 4, 4, Fix128::ONE);
let pos = Vec3Fix::new(
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
);
let vel = Vec3Fix::new(Fix128::from_int(4), Fix128::ZERO, Fix128::ZERO);
p2g_trilinear(&mut g, pos, vel);
assert_eq!(g.u(1, 1, 1), Fix128::from_int(2));
assert_eq!(g.u(2, 1, 1), Fix128::from_int(2));
}
#[test]
fn g2p_trilinear_between_faces_interpolates() {
let mut g = MacGrid::new(4, 4, 4, Fix128::ONE);
let i0 = g.idx_u(1, 1, 1);
let i1 = g.idx_u(2, 1, 1);
g.u[i0] = Fix128::from_int(10);
g.u[i1] = Fix128::from_int(20);
let pos = Vec3Fix::new(
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
Fix128::from_ratio(15, 10),
);
let v = g2p_velocity(&g, pos);
assert_eq!(v.x, Fix128::from_int(15));
}
#[test]
fn split_negative_position_clamps_to_zero() {
let (base, _) = split(Fix128::from_int(-3), Fix128::ZERO);
assert_eq!(base, 0);
}
fn seed_divergent_flow(nx: usize) -> MacGrid {
let mut g = MacGrid::new(nx, nx, nx, Fix128::ONE);
for i in 0..=nx {
for j in 0..nx {
for k in 0..nx {
let ix = g.idx_u(i, j, k);
if ix < g.u.len() {
g.u[ix] = Fix128::from_int(i as i64);
}
}
}
}
g
}
#[test]
fn bicgstab_reduces_divergence_below_jacobi_iterations() {
let mut g = seed_divergent_flow(4);
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let div_before = g.divergence(2, 2, 2).abs();
let stats = project_pressure_bicgstab(&mut g, dt, rho, 15, Fix128::from_ratio(1, 10_000));
let div_after = g.divergence(2, 2, 2).abs();
assert!(div_after < div_before);
assert!(stats.iterations <= 15);
}
#[test]
fn bicgstab_matches_jacobi_within_tolerance() {
let mut g_bicg = seed_divergent_flow(4);
let mut g_jac = g_bicg.clone();
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let _ = project_pressure_bicgstab(&mut g_bicg, dt, rho, 30, Fix128::from_ratio(1, 100_000));
project_pressure_jacobi(&mut g_jac, dt, rho, 200);
let div_bicg = g_bicg.divergence(2, 2, 2).abs();
let div_jac = g_jac.divergence(2, 2, 2).abs();
assert!(div_bicg < Fix128::from_ratio(1, 10));
assert!(div_jac < Fix128::from_ratio(1, 10));
}
#[derive(Clone, Copy)]
enum VisitOrder {
Natural,
Reverse,
Strided,
}
fn colour_cells(nx: usize, ny: usize, nz: usize, colour: u32, order: VisitOrder) -> Vec<usize> {
let mut v = Vec::new();
for k in 0..nz {
for j in 0..ny {
for i in 0..nx {
if ((i + j + k) as u32 % 2) == colour {
v.push(i + nx * (j + ny * k));
}
}
}
}
match order {
VisitOrder::Natural => v,
VisitOrder::Reverse => {
v.reverse();
v
}
VisitOrder::Strided => {
let n = v.len();
let mut out = Vec::with_capacity(n);
let mut seen = vec![false; n];
let mut idx = 0usize;
for _ in 0..n {
while seen[idx] {
idx = (idx + 1) % n;
}
seen[idx] = true;
out.push(v[idx]);
idx = (idx + 7) % n;
}
out
}
}
}
fn reference_red_black(
grid: &mut MacGrid,
dt_s: Fix128,
density: Fix128,
iterations: u32,
order: VisitOrder,
) {
grid.enforce_face_boundaries();
let scale = density * grid.dx * grid.dx / dt_s;
let n = grid.nx * grid.ny * grid.nz;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
let (nx, ny, nz) = (grid.nx, grid.ny, grid.nz);
let plan = [
colour_cells(nx, ny, nz, 0, order),
colour_cells(nx, ny, nz, 1, order),
];
for _ in 0..iterations {
for cells in &plan {
for &idx in cells {
let i = idx % nx;
let j = (idx / nx) % ny;
let k = idx / (nx * ny);
let neighbours = mask.neighbour_sum(&grid.pressure, i, j, k);
grid.pressure[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density * inv_dx);
}
fn reference_plain_gs(grid: &mut MacGrid, dt_s: Fix128, density: Fix128, iterations: u32) {
grid.enforce_face_boundaries();
let scale = density * grid.dx * grid.dx / dt_s;
let n = grid.nx * grid.ny * grid.nz;
let rhs = poisson_rhs(grid, scale);
let mask = PoissonMask::from_grid(grid);
let inv_deg = inverse_degrees(&mask, n);
for _ in 0..iterations {
for k in 0..grid.nz {
for j in 0..grid.ny {
for i in 0..grid.nx {
let idx = i + grid.nx * (j + grid.ny * k);
let neighbours = mask.neighbour_sum(&grid.pressure, i, j, k);
grid.pressure[idx] = (neighbours - rhs[idx]) * inv_deg[idx];
}
}
}
}
let inv_dx = Fix128::ONE / grid.dx;
subtract_pressure_gradient(grid, dt_s / density * inv_dx);
}
fn grids_are_bit_equal(a: &MacGrid, b: &MacGrid) -> bool {
a.pressure == b.pressure && a.u == b.u && a.v == b.v && a.w == b.w
}
#[test]
fn red_black_sweep_is_independent_of_visit_order() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_divergent_flow(5);
let mut solver_grid = base.clone();
project_pressure_red_black_gs(&mut solver_grid, dt, rho, 8);
for order in [
VisitOrder::Natural,
VisitOrder::Reverse,
VisitOrder::Strided,
] {
let mut reference = base.clone();
reference_red_black(&mut reference, dt, rho, 8, order);
assert!(
grids_are_bit_equal(&solver_grid, &reference),
"red-black result changed with the visit order: the sweep is not \
order-independent, so the parallel path cannot be bit-identical",
);
}
}
#[test]
fn slab_decomposition_reproduces_the_monolithic_pressure_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
for &(n, ranks) in &[
(8usize, 1usize),
(8, 2),
(8, 4),
(8, 8),
(7, 2),
(7, 3),
(7, 4),
(5, 4),
(3, 4),
] {
let base = seed_divergent_flow(n);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let mut split = base.clone();
project_pressure_decomposed(&mut split, dt, rho, 6, ranks, HaloSchedule::EverySweep);
assert!(
grids_are_bit_equal(&monolithic, &split),
"{n}³ grid split across {ranks} slabs did not reproduce the \
monolithic solve: either one halo layer is not enough for the \
7-point stencil, or a rank read past its halo",
);
}
}
#[test]
fn a_halo_exchanged_once_per_iteration_diverges_from_the_monolithic_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_divergent_flow(8);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let mut stale = base.clone();
project_pressure_decomposed(&mut stale, dt, rho, 6, 4, HaloSchedule::EveryIteration);
assert!(
!grids_are_bit_equal(&monolithic, &stale),
"a halo one sweep out of date still reproduced the monolithic solve, \
so the bit-equality test above is not actually testing the exchange",
);
}
struct StagedTransport {
cells: usize,
plane: usize,
flat: Vec<Fix128>,
mailbox: Vec<Fix128>,
}
impl StagedTransport {
fn new(ranks: usize, cells: usize, plane: usize) -> Self {
Self {
cells,
plane,
flat: vec![Fix128::ZERO; ranks * cells],
mailbox: vec![Fix128::ZERO; plane],
}
}
}
impl RankTransport for StagedTransport {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
let base = rank * self.cells;
&mut self.flat[base..base + self.cells]
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
let off = layer * self.plane;
let from = src * self.cells + off;
self.mailbox
.copy_from_slice(&self.flat[from..from + self.plane]);
let to = dst * self.cells + off;
self.flat[to..to + self.plane].copy_from_slice(&self.mailbox);
}
}
struct UndeliveredTransport(LocalTransport);
impl RankTransport for UndeliveredTransport {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
self.0.slab_mut(rank)
}
fn deliver_layer(&mut self, _src: usize, _dst: usize, _layer: usize) {}
}
#[test]
fn a_second_rank_transport_reproduces_the_in_process_one_bit_for_bit() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
for &(n, ranks) in &[(8usize, 2usize), (8, 4), (8, 8), (7, 3), (5, 4), (3, 4)] {
let base = seed_divergent_flow(n);
let mut via_local = base.clone();
project_pressure_decomposed(
&mut via_local,
dt,
rho,
6,
ranks,
HaloSchedule::EverySweep,
);
let mut via_staged = base.clone();
let mut staged = StagedTransport::new(ranks, n * n * n, n * n);
project_pressure_decomposed_over(
&mut via_staged,
dt,
rho,
6,
ranks,
HaloSchedule::EverySweep,
&mut staged,
);
assert!(
grids_are_bit_equal(&via_local, &via_staged),
"{n}³ grid over {ranks} slabs: the staged transport disagreed with \
the in-process one, so the solve depends on how a layer is \
carried and not only on which layer arrives when",
);
}
}
#[test]
fn a_transport_that_never_delivers_does_not_reproduce_the_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_divergent_flow(8);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let mut undelivered = base.clone();
let mut transport = UndeliveredTransport(LocalTransport::new(4, 8, 8 * 8));
project_pressure_decomposed_over(
&mut undelivered,
dt,
rho,
6,
4,
HaloSchedule::EverySweep,
&mut transport,
);
assert!(
!grids_are_bit_equal(&monolithic, &undelivered),
"slabs that never received a halo still reproduced the monolithic \
solve, so the transport tests above are not measuring the delivery",
);
}
#[cfg(feature = "std")]
const XPROC_ROLE: &str = "ALICE_PHYSICS_XPROC_ROLE";
#[cfg(feature = "std")]
const XPROC_PORT: &str = "ALICE_PHYSICS_XPROC_PORT";
#[cfg(feature = "std")]
const XPROC_SIZE: &str = "ALICE_PHYSICS_XPROC_SIZE";
#[cfg(feature = "std")]
const XPROC_FAULT: &str = "ALICE_PHYSICS_XPROC_FAULT";
#[cfg(feature = "std")]
const XPROC_CHILD: &str = "eulerian_grid::tests::cross_process_rank1_child";
#[cfg(feature = "std")]
const XPROC_ITERATIONS: u32 = 6;
#[cfg(feature = "std")]
fn xproc_dt() -> Fix128 {
Fix128::from_ratio(1, 100)
}
#[cfg(feature = "std")]
fn xproc_rho() -> Fix128 {
Fix128::from_int(1000)
}
#[cfg(feature = "std")]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum CrossFault {
None,
ShiftOneLayer,
DropFirstDelivery,
ChunkedWire,
}
#[cfg(feature = "std")]
impl CrossFault {
fn as_str(self) -> &'static str {
match self {
Self::None => "none",
Self::ShiftOneLayer => "shift",
Self::DropFirstDelivery => "drop",
Self::ChunkedWire => "chunked",
}
}
fn parse(name: &str) -> Self {
match name {
"none" => Self::None,
"shift" => Self::ShiftOneLayer,
"drop" => Self::DropFirstDelivery,
"chunked" => Self::ChunkedWire,
other => panic!("unknown cross-process fault `{other}`"),
}
}
}
#[cfg(feature = "std")]
struct FaultyCrossing<S> {
inner: SocketTransport<S>,
fault: CrossFault,
seen: usize,
}
#[cfg(feature = "std")]
impl<S: std::io::Read + std::io::Write> RankTransport for FaultyCrossing<S> {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
self.inner.slab_mut(rank)
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
let index = self.seen;
self.seen += 1;
match self.fault {
CrossFault::None => self.inner.deliver_layer(src, dst, layer),
CrossFault::ShiftOneLayer => {
self.inner.deliver_layer(src, dst, layer.saturating_sub(1));
}
CrossFault::DropFirstDelivery => {
if index > 0 {
self.inner.deliver_layer(src, dst, layer);
}
}
CrossFault::ChunkedWire => unreachable!(
"the chunked wire is a stream wrapper, not a delivery fault; \
`solve_as_rank` routes it around `FaultyCrossing`"
),
}
}
}
#[cfg(feature = "std")]
struct Chunked<S> {
inner: S,
read_calls: usize,
read_bytes: usize,
write_calls: usize,
write_bytes: usize,
}
#[cfg(feature = "std")]
impl<S: std::io::Read> std::io::Read for Chunked<S> {
fn read(&mut self, buf: &mut [u8]) -> std::io::Result<usize> {
if buf.is_empty() {
return Ok(0);
}
let n = self.inner.read(&mut buf[..1])?;
if n > 0 {
self.read_calls += 1;
self.read_bytes += n;
}
Ok(n)
}
}
#[cfg(feature = "std")]
impl<S: std::io::Write> std::io::Write for Chunked<S> {
fn write(&mut self, buf: &[u8]) -> std::io::Result<usize> {
if buf.is_empty() {
return Ok(0);
}
let n = self.inner.write(&buf[..1])?;
if n > 0 {
self.write_calls += 1;
self.write_bytes += n;
}
Ok(n)
}
fn flush(&mut self) -> std::io::Result<()> {
self.inner.flush()
}
}
#[cfg(feature = "std")]
fn solve_as_rank(
grid: &mut MacGrid,
n: usize,
my_rank: usize,
fault: CrossFault,
link: std::net::TcpStream,
) {
use std::time::Duration;
link.set_read_timeout(Some(Duration::from_secs(30)))
.expect("set a read timeout on the peer link");
link.set_write_timeout(Some(Duration::from_secs(30)))
.expect("set a write timeout on the peer link");
if fault == CrossFault::ChunkedWire {
let mut links: Vec<Option<Chunked<std::net::TcpStream>>> = vec![None, None];
links[1 - my_rank] = Some(Chunked {
inner: link,
read_calls: 0,
read_bytes: 0,
write_calls: 0,
write_bytes: 0,
});
let mut transport = SocketTransport::new(my_rank, n * n * n, n * n, links);
project_pressure_decomposed_on_rank(
grid,
xproc_dt(),
xproc_rho(),
XPROC_ITERATIONS,
2,
HaloSchedule::EverySweep,
my_rank,
&mut transport,
);
let wire = transport.links[1 - my_rank]
.as_ref()
.expect("the peer link is still held by the transport");
let one_layer = n * n * WIRE_BYTES_PER_CELL;
assert!(
wire.read_calls == wire.read_bytes
&& wire.write_calls == wire.write_bytes
&& wire.read_bytes >= one_layer
&& wire.write_bytes >= one_layer,
"rank {my_rank}: the chunked wire saw {} reads moving {} bytes and {} \
writes moving {} bytes; the gate needs calls == bytes in both \
directions and at least one layer ({one_layer} bytes) each way",
wire.read_calls,
wire.read_bytes,
wire.write_calls,
wire.write_bytes
);
return;
}
let mut links: Vec<Option<std::net::TcpStream>> = vec![None, None];
links[1 - my_rank] = Some(link);
let socket = SocketTransport::new(my_rank, n * n * n, n * n, links);
if fault == CrossFault::None {
let mut transport = socket;
project_pressure_decomposed_on_rank(
grid,
xproc_dt(),
xproc_rho(),
XPROC_ITERATIONS,
2,
HaloSchedule::EverySweep,
my_rank,
&mut transport,
);
} else {
let mut transport = FaultyCrossing {
inner: socket,
fault,
seen: 0,
};
project_pressure_decomposed_on_rank(
grid,
xproc_dt(),
xproc_rho(),
XPROC_ITERATIONS,
2,
HaloSchedule::EverySweep,
my_rank,
&mut transport,
);
}
}
#[cfg(feature = "std")]
#[test]
fn cross_process_rank1_child() {
let Ok(role) = std::env::var(XPROC_ROLE) else {
return;
};
assert_eq!(role, "rank1", "unknown cross-process role `{role}`");
let n: usize = std::env::var(XPROC_SIZE)
.expect("grid size from the parent")
.parse()
.expect("grid size is a number");
let port: u16 = std::env::var(XPROC_PORT)
.expect("loopback port from the parent")
.parse()
.expect("loopback port is a number");
let fault =
CrossFault::parse(&std::env::var(XPROC_FAULT).expect("fault mode from the parent"));
let link = std::net::TcpStream::connect(("127.0.0.1", port))
.expect("connect to rank 0 on the loopback port it is listening on");
let mut grid = seed_divergent_flow(n);
solve_as_rank(&mut grid, n, 1, fault, link);
assert!(
grid.pressure.iter().all(|p| *p == Fix128::ZERO),
"rank 1 assembled a pressure field into its own grid: the gather is \
supposed to land on rank 0 alone",
);
}
#[cfg(feature = "std")]
fn pressure_field_from_two_processes(n: usize, fault: CrossFault) -> MacGrid {
use std::io::ErrorKind;
use std::net::TcpListener;
use std::process::{Command, Stdio};
use std::time::{Duration, Instant};
let listener =
TcpListener::bind(("127.0.0.1", 0)).expect("bind a loopback port for rank 1");
let port = listener
.local_addr()
.expect("the bound loopback address")
.port();
let exe = std::env::current_exe().expect("path of this test binary");
let mut child = Command::new(exe)
.args(["--exact", XPROC_CHILD, "--test-threads=1"])
.env(XPROC_ROLE, "rank1")
.env(XPROC_PORT, port.to_string())
.env(XPROC_SIZE, n.to_string())
.env(XPROC_FAULT, fault.as_str())
.stdin(Stdio::null())
.stdout(Stdio::null())
.stderr(Stdio::null())
.spawn()
.expect("re-execute this test binary as rank 1");
listener
.set_nonblocking(true)
.expect("poll for rank 1 instead of blocking on it");
let deadline = Instant::now() + Duration::from_secs(60);
let link = loop {
match listener.accept() {
Ok((stream, _)) => break stream,
Err(e) if e.kind() == ErrorKind::WouldBlock => {}
Err(e) => panic!("accepting rank 1 failed: {e}"),
}
if let Some(status) = child.try_wait().expect("poll rank 1") {
panic!(
"rank 1 exited ({status}) without connecting: does \
`{XPROC_CHILD}` still name a test?"
);
}
assert!(
Instant::now() < deadline,
"rank 1 did not connect within 60 s",
);
std::thread::sleep(Duration::from_millis(2));
};
link.set_nonblocking(false)
.expect("switch the accepted stream back to blocking");
let mut grid = seed_divergent_flow(n);
solve_as_rank(&mut grid, n, 0, fault, link);
let status = child.wait().expect("wait for rank 1 to exit");
assert!(status.success(), "rank 1 exited with {status}");
grid
}
#[cfg(feature = "std")]
#[test]
fn two_processes_reproduce_the_single_process_pressure_solve() {
for &n in &[8usize, 7] {
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(
&mut monolithic,
xproc_dt(),
xproc_rho(),
XPROC_ITERATIONS,
);
let distributed = pressure_field_from_two_processes(n, CrossFault::None);
assert_eq!(
monolithic.pressure, distributed.pressure,
"{n}³ over two processes produced a different pressure field \
than the single-process solve",
);
assert!(
grids_are_bit_equal(&monolithic, &distributed),
"{n}³ over two processes matched on pressure but not on the \
projected velocities",
);
}
}
#[cfg(feature = "std")]
#[test]
fn a_cross_process_halo_shifted_by_one_layer_does_not_reproduce_the_solve() {
let n = 8usize;
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(&mut monolithic, xproc_dt(), xproc_rho(), XPROC_ITERATIONS);
let shifted = pressure_field_from_two_processes(n, CrossFault::ShiftOneLayer);
assert_ne!(
monolithic.pressure, shifted.pressure,
"halos delivered one layer off still reproduced the single-process \
solve, so the two-process test is not measuring what crosses the \
boundary",
);
}
#[cfg(feature = "std")]
#[test]
fn a_cross_process_delivery_dropped_once_does_not_reproduce_the_solve() {
let n = 8usize;
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(&mut monolithic, xproc_dt(), xproc_rho(), XPROC_ITERATIONS);
let dropped = pressure_field_from_two_processes(n, CrossFault::DropFirstDelivery);
assert_ne!(
monolithic.pressure, dropped.pressure,
"a run missing one halo delivery still reproduced the single-process \
solve, so the two-process test is not measuring the delivery",
);
}
#[cfg(feature = "std")]
#[test]
fn a_cross_process_wire_split_into_single_bytes_reproduces_the_solve() {
for &n in &[8usize, 7] {
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(
&mut monolithic,
xproc_dt(),
xproc_rho(),
XPROC_ITERATIONS,
);
let chunked = pressure_field_from_two_processes(n, CrossFault::ChunkedWire);
assert_eq!(
monolithic.pressure, chunked.pressure,
"{n}³ over a wire split into single bytes produced a different \
pressure field than the single-process solve",
);
assert!(
grids_are_bit_equal(&monolithic, &chunked),
"{n}³ over a single-byte wire matched on pressure but not on the \
projected velocities",
);
}
}
#[test]
fn the_rank_local_driver_on_a_single_rank_matches_the_monolithic_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let n = 8usize;
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let mut alone = seed_divergent_flow(n);
let mut transport = LocalTransport::new(1, n, n * n);
project_pressure_decomposed_on_rank(
&mut alone,
dt,
rho,
6,
1,
HaloSchedule::EverySweep,
0,
&mut transport,
);
assert!(
grids_are_bit_equal(&monolithic, &alone),
"the rank-local driver on a single rank did not reproduce the \
monolithic solve",
);
}
struct DeliveryRecorder {
calls: Vec<(usize, usize, usize)>,
}
impl RankTransport for DeliveryRecorder {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
unreachable!("a schedule walk asked the recorder for rank {rank}'s slab")
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
self.calls.push((src, dst, layer));
}
}
#[allow(clippy::type_complexity)]
fn recorded_schedule(
nz: usize,
ranks: usize,
) -> (Vec<(usize, usize, usize)>, Vec<(usize, usize, usize)>) {
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let mut halo = DeliveryRecorder { calls: Vec::new() };
exchange_slab_halos(&mut halo, &bounds, nz);
let mut gather = DeliveryRecorder { calls: Vec::new() };
gather_slabs_to_root(&mut gather, &bounds);
(halo.calls, gather.calls)
}
#[test]
fn only_more_than_two_ranks_give_a_rank_deliveries_it_is_not_party_to() {
for ranks in 1usize..=6 {
for nz in [ranks, 2 * ranks + 1, 3 * ranks] {
let (halo, gather) = recorded_schedule(nz, ranks);
let bystanding = |r: usize| {
halo.iter()
.chain(gather.iter())
.filter(|&&(src, dst, _)| src != r && dst != r)
.count()
};
if ranks <= 2 {
for r in 0..ranks {
assert_eq!(
bystanding(r),
0,
"{nz} layers over {ranks} ranks gave rank {r} a delivery it \
is not party to, so the two-rank harness does exercise the \
no-op branch after all",
);
}
} else {
assert!(
bystanding(0) > 0,
"{nz} layers over {ranks} ranks left rank 0 party to every \
delivery, so running it as rank 0 of a multi-process solve \
would not reach the no-op branch",
);
}
}
}
}
#[test]
fn the_halo_exchange_couples_only_neighbouring_slabs() {
for ranks in 1usize..=6 {
for nz in [1, ranks, 2 * ranks + 1, 3 * ranks] {
let bounds: Vec<(usize, usize)> =
(0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let (halo, gather) = recorded_schedule(nz, ranks);
for &(src, dst, _) in &halo {
assert_ne!(src, dst, "a halo delivery had one rank at both ends");
let (lo, hi) = (src.min(dst), src.max(dst));
for (r, &(k0, k1)) in bounds.iter().enumerate().take(hi).skip(lo + 1) {
assert_eq!(
k0, k1,
"the halo delivery {src} -> {dst} reached over rank {r}, \
which owns layers {k0}..{k1}",
);
}
}
for &(src, dst, layer) in &gather {
assert_eq!(dst, 0, "a gather delivery landed on rank {dst}, not root");
assert_ne!(src, 0, "the root was asked to gather from itself");
assert_eq!(
slab_owner(&bounds, layer),
Some(src),
"rank {src} gathered layer {layer}, which it does not own",
);
}
}
}
let (halo, gather) = recorded_schedule(9, 3);
assert!(
halo.iter()
.all(|&(src, dst, _)| (src, dst) != (0, 2) && (src, dst) != (2, 0)),
"a halo crossed directly between the first and last of three slabs",
);
assert!(
gather.iter().any(|&(src, dst, _)| (src, dst) == (2, 0)),
"the gather never crossed from the last of three slabs to the root",
);
}
#[cfg(feature = "std")]
const XPROC3_RANK: &str = "ALICE_PHYSICS_XPROC3_RANK";
#[cfg(feature = "std")]
const XPROC3_PORT: &str = "ALICE_PHYSICS_XPROC3_PORT";
#[cfg(feature = "std")]
const XPROC3_SIZE: &str = "ALICE_PHYSICS_XPROC3_SIZE";
#[cfg(feature = "std")]
const XPROC3_FAULT: &str = "ALICE_PHYSICS_XPROC3_FAULT";
#[cfg(feature = "std")]
const XPROC3_CHILD: &str = "eulerian_grid::tests::cross_process_three_rank_child";
#[cfg(feature = "std")]
const XPROC3_ITERATIONS: u32 = 6;
#[cfg(feature = "std")]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum BystanderMode {
Skip,
Deliver,
}
#[cfg(feature = "std")]
impl BystanderMode {
fn as_str(self) -> &'static str {
match self {
Self::Skip => "skip",
Self::Deliver => "deliver",
}
}
fn parse(name: &str) -> Self {
match name {
"skip" => Self::Skip,
"deliver" => Self::Deliver,
other => panic!("unknown bystander mode `{other}`"),
}
}
}
#[cfg(feature = "std")]
struct BystanderDelivers<S> {
inner: SocketTransport<S>,
}
#[cfg(feature = "std")]
impl<S: std::io::Read + std::io::Write> RankTransport for BystanderDelivers<S> {
fn slab_mut(&mut self, rank: usize) -> &mut [Fix128] {
self.inner.slab_mut(rank)
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
let mine = self.inner.my_rank;
if mine == src || mine == dst {
self.inner.deliver_layer(src, dst, layer);
} else {
self.inner.deliver_layer(mine, dst, layer);
}
}
}
#[cfg(feature = "std")]
fn bound_peer_link(link: &std::net::TcpStream) {
use std::time::Duration;
link.set_read_timeout(Some(Duration::from_secs(30)))
.expect("set a read timeout on a peer link");
link.set_write_timeout(Some(Duration::from_secs(30)))
.expect("set a write timeout on a peer link");
}
#[cfg(feature = "std")]
fn reap_child_ranks(kids: &mut [std::process::Child]) {
use std::time::{Duration, Instant};
let deadline = Instant::now() + Duration::from_secs(60);
for kid in kids {
loop {
match kid.try_wait() {
Ok(Some(_)) | Err(_) => break,
Ok(None) if Instant::now() >= deadline => {
let _ = kid.kill();
let _ = kid.wait();
break;
}
Ok(None) => std::thread::sleep(Duration::from_millis(2)),
}
}
}
}
#[cfg(feature = "std")]
fn accept_before_deadline(listener: &std::net::TcpListener, who: &str) -> std::net::TcpStream {
use std::io::ErrorKind;
use std::time::{Duration, Instant};
listener
.set_nonblocking(true)
.expect("poll for a peer rank instead of blocking on it");
let deadline = Instant::now() + Duration::from_secs(60);
loop {
match listener.accept() {
Ok((stream, _)) => {
stream
.set_nonblocking(false)
.expect("switch an accepted stream back to blocking");
bound_peer_link(&stream);
return stream;
}
Err(e) if e.kind() == ErrorKind::WouldBlock => {}
Err(e) => panic!("accepting {who} failed: {e}"),
}
assert!(
Instant::now() < deadline,
"{who} did not connect within 60 s"
);
std::thread::sleep(Duration::from_millis(2));
}
}
#[cfg(feature = "std")]
struct SpawnedRanks {
kids: Vec<std::process::Child>,
}
#[cfg(feature = "std")]
impl Drop for SpawnedRanks {
fn drop(&mut self) {
reap_child_ranks(&mut self.kids);
}
}
#[cfg(feature = "std")]
impl SpawnedRanks {
fn accept(&mut self, listener: &std::net::TcpListener, rank: usize) -> std::net::TcpStream {
use std::io::ErrorKind;
use std::time::{Duration, Instant};
listener
.set_nonblocking(true)
.expect("poll for a child rank instead of blocking on it");
let deadline = Instant::now() + Duration::from_secs(60);
let link = loop {
match listener.accept() {
Ok((stream, _)) => break stream,
Err(e) if e.kind() == ErrorKind::WouldBlock => {}
Err(e) => panic!("accepting rank {rank} failed: {e}"),
}
for (slot, kid) in self.kids.iter_mut().enumerate() {
if let Some(status) = kid.try_wait().expect("poll a child rank") {
panic!(
"rank {} exited ({status}) before rank {rank} connected: does \
`{XPROC3_CHILD}` still name a test?",
slot + 1,
);
}
}
assert!(
Instant::now() < deadline,
"rank {rank} did not connect within 60 s",
);
std::thread::sleep(Duration::from_millis(2));
};
link.set_nonblocking(false)
.expect("switch the accepted stream back to blocking");
bound_peer_link(&link);
link
}
fn join(&mut self) {
for (slot, kid) in self.kids.iter_mut().enumerate() {
let status = kid.wait().expect("wait for a child rank to exit");
assert!(status.success(), "rank {} exited with {status}", slot + 1);
}
}
}
#[cfg(feature = "std")]
fn solve_as_rank_of_three(
grid: &mut MacGrid,
n: usize,
my_rank: usize,
mode: BystanderMode,
links: Vec<Option<std::net::TcpStream>>,
) {
let socket = SocketTransport::new(my_rank, n * n * n, n * n, links);
match mode {
BystanderMode::Skip => {
let mut transport = socket;
project_pressure_decomposed_on_rank(
grid,
xproc_dt(),
xproc_rho(),
XPROC3_ITERATIONS,
3,
HaloSchedule::EverySweep,
my_rank,
&mut transport,
);
}
BystanderMode::Deliver => {
let mut transport = BystanderDelivers { inner: socket };
project_pressure_decomposed_on_rank(
grid,
xproc_dt(),
xproc_rho(),
XPROC3_ITERATIONS,
3,
HaloSchedule::EverySweep,
my_rank,
&mut transport,
);
}
}
}
#[cfg(feature = "std")]
#[test]
fn cross_process_three_rank_child() {
use std::io::{Read, Write};
use std::net::{TcpListener, TcpStream};
let Ok(rank) = std::env::var(XPROC3_RANK) else {
return;
};
let my_rank: usize = rank.parse().expect("the rank is a number");
assert!(
my_rank == 1 || my_rank == 2,
"rank {my_rank} is not a child rank of the three-rank harness",
);
let n: usize = std::env::var(XPROC3_SIZE)
.expect("grid size from rank 0")
.parse()
.expect("grid size is a number");
let port: u16 = std::env::var(XPROC3_PORT)
.expect("loopback port from rank 0")
.parse()
.expect("loopback port is a number");
let mode =
BystanderMode::parse(&std::env::var(XPROC3_FAULT).expect("bystander mode from rank 0"));
let mut root = TcpStream::connect(("127.0.0.1", port))
.expect("connect to rank 0 on the loopback port it is listening on");
bound_peer_link(&root);
let mut links: Vec<Option<TcpStream>> = vec![None, None, None];
if my_rank == 1 {
let peer_listener =
TcpListener::bind(("127.0.0.1", 0)).expect("bind a loopback port for rank 2");
let peer_port = peer_listener
.local_addr()
.expect("the bound loopback address")
.port();
root.write_all(&peer_port.to_le_bytes())
.expect("tell rank 0 where rank 2 should dial");
root.flush().expect("flush the peer port to rank 0");
links[2] = Some(accept_before_deadline(&peer_listener, "rank 2"));
links[0] = Some(root);
} else {
let mut peer_port = [0u8; 2];
root.read_exact(&mut peer_port)
.expect("learn from rank 0 where rank 1 is listening");
let peer = TcpStream::connect(("127.0.0.1", u16::from_le_bytes(peer_port)))
.expect("connect to rank 1 on the port it is listening on");
bound_peer_link(&peer);
links[1] = Some(peer);
links[0] = Some(root);
}
let mut grid = seed_divergent_flow(n);
solve_as_rank_of_three(&mut grid, n, my_rank, mode, links);
assert!(
grid.pressure.iter().all(|p| *p == Fix128::ZERO),
"rank {my_rank} assembled a pressure field into its own grid: the \
gather is supposed to land on rank 0 alone",
);
}
#[cfg(feature = "std")]
fn pressure_field_from_three_processes(n: usize, mode: BystanderMode) -> MacGrid {
use std::io::{Read, Write};
use std::net::TcpListener;
use std::process::{Command, Stdio};
let listeners: Vec<TcpListener> = (1usize..=2)
.map(|rank| {
TcpListener::bind(("127.0.0.1", 0))
.unwrap_or_else(|e| panic!("bind a loopback port for rank {rank}: {e}"))
})
.collect();
let exe = std::env::current_exe().expect("path of this test binary");
let mut ranks = SpawnedRanks {
kids: listeners
.iter()
.enumerate()
.map(|(slot, listener)| {
let port = listener
.local_addr()
.expect("the bound loopback address")
.port();
Command::new(&exe)
.args(["--exact", XPROC3_CHILD, "--test-threads=1"])
.env(XPROC3_RANK, (slot + 1).to_string())
.env(XPROC3_PORT, port.to_string())
.env(XPROC3_SIZE, n.to_string())
.env(XPROC3_FAULT, mode.as_str())
.stdin(Stdio::null())
.stdout(Stdio::null())
.stderr(Stdio::null())
.spawn()
.unwrap_or_else(|e| {
panic!("re-execute this test binary as rank {}: {e}", slot + 1)
})
})
.collect(),
};
let mut to_rank1 = ranks.accept(&listeners[0], 1);
let mut to_rank2 = ranks.accept(&listeners[1], 2);
let mut peer_port = [0u8; 2];
to_rank1
.read_exact(&mut peer_port)
.expect("learn from rank 1 where rank 2 should dial");
to_rank2
.write_all(&peer_port)
.expect("tell rank 2 where rank 1 is listening");
to_rank2.flush().expect("flush the peer port to rank 2");
let mut grid = seed_divergent_flow(n);
solve_as_rank_of_three(
&mut grid,
n,
0,
mode,
vec![None, Some(to_rank1), Some(to_rank2)],
);
ranks.join();
grid
}
#[cfg(feature = "std")]
#[test]
fn three_processes_reproduce_the_single_process_pressure_solve() {
for &n in &[9usize, 8] {
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(
&mut monolithic,
xproc_dt(),
xproc_rho(),
XPROC3_ITERATIONS,
);
let distributed = pressure_field_from_three_processes(n, BystanderMode::Skip);
assert_eq!(
monolithic.pressure, distributed.pressure,
"{n}³ over three processes produced a different pressure field \
than the single-process solve",
);
assert!(
grids_are_bit_equal(&monolithic, &distributed),
"{n}³ over three processes matched on pressure but not on the \
projected velocities",
);
}
}
#[cfg(feature = "std")]
#[test]
fn a_rank_that_performs_deliveries_it_is_not_party_to_does_not_reproduce_the_solve() {
let n = 6usize;
let mut monolithic = seed_divergent_flow(n);
project_pressure_red_black_gs(&mut monolithic, xproc_dt(), xproc_rho(), XPROC3_ITERATIONS);
let skipped = pressure_field_from_three_processes(n, BystanderMode::Skip);
assert_eq!(
monolithic.pressure, skipped.pressure,
"{n}³ over three processes did not reproduce the single-process \
solve, so the comparison below is not isolating the mode",
);
let meddled = std::panic::catch_unwind(|| {
pressure_field_from_three_processes(n, BystanderMode::Deliver)
});
match meddled {
Ok(field) => assert_ne!(
monolithic.pressure, field.pressure,
"ranks that performed the deliveries they are not party to still \
reproduced the single-process solve, so the three-process test \
is not measuring the bystander branch",
),
Err(_) => {
}
}
}
#[test]
fn slab_bounds_partition_every_layer_exactly_once() {
for nz in 1usize..=16 {
for ranks in 1usize..=20 {
let bounds: Vec<(usize, usize)> =
(0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
for k in 0..nz {
let owners = bounds.iter().filter(|&&(k0, k1)| k0 <= k && k < k1).count();
assert_eq!(
owners, 1,
"layer {k} of {nz} has {owners} owners across {ranks} ranks",
);
}
let owned: usize = bounds.iter().map(|&(k0, k1)| k1 - k0).sum();
assert_eq!(owned, nz, "{ranks} ranks over {nz} layers own {owned}");
}
}
}
#[test]
fn colour_blind_gauss_seidel_drifts_from_red_black() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_divergent_flow(5);
let mut red_black = base.clone();
project_pressure_red_black_gs(&mut red_black, dt, rho, 8);
let mut colour_blind = base.clone();
reference_plain_gs(&mut colour_blind, dt, rho, 8);
assert!(
!grids_are_bit_equal(&red_black, &colour_blind),
"a colour-blind sweep agreed with the red-black one, so the \
order-independence test has no discriminating power",
);
}
#[test]
fn bicgstab_stats_reports_convergence_status() {
let mut g = seed_divergent_flow(4);
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let stats =
project_pressure_bicgstab(&mut g, dt, rho, 50, Fix128::from_ratio(1, 1_000_000));
assert!(stats.iterations >= 1);
assert!(stats.iterations <= 50);
}
#[test]
fn bicgstab_no_op_on_zero_dt() {
let mut g = seed_divergent_flow(3);
let stats = project_pressure_bicgstab(
&mut g,
Fix128::ZERO,
Fix128::from_int(1000),
10,
Fix128::from_ratio(1, 10_000),
);
assert_eq!(stats.iterations, 0);
assert!(stats.converged);
}
#[test]
fn poisson_operator_diagonal_matches_the_preconditioner() {
let nx = 3;
let ny = 3;
let nz = 3;
let mut grid = MacGrid::new(nx, ny, nz, Fix128::from_ratio(1, 8));
grid.set_u_solid(1, 1, 1, true); let mask = PoissonMask::from_grid(&grid);
let n = nx * ny * nz;
for c in 0..n {
let mut e = vec![Fix128::ZERO; n];
e[c] = Fix128::ONE;
let mut out = vec![Fix128::ZERO; n];
mask.apply(&e, &mut out);
assert_eq!(
out[c],
Fix128::from_int(-mask.degree(c)),
"diagonal of cell {c} disagrees with degree()"
);
}
}
#[test]
fn poisson_operator_is_symmetric() {
let nx = 3;
let ny = 3;
let nz = 2;
let mut grid = MacGrid::new(nx, ny, nz, Fix128::from_ratio(1, 8));
grid.set_closed_box_walls();
grid.set_v_solid(1, 1, 0, true);
let mask = PoissonMask::from_grid(&grid);
let n = nx * ny * nz;
let mut columns = Vec::with_capacity(n);
for c in 0..n {
let mut e = vec![Fix128::ZERO; n];
e[c] = Fix128::ONE;
let mut out = vec![Fix128::ZERO; n];
mask.apply(&e, &mut out);
columns.push(out);
}
for (a, col_a) in columns.iter().enumerate() {
for (b, col_b) in columns.iter().enumerate() {
assert_eq!(col_a[b], col_b[a], "A[{b},{a}] != A[{a},{b}]");
}
}
}
#[test]
fn closed_box_degree_drops_to_the_open_face_count() {
let mut grid = MacGrid::new(2, 2, 1, Fix128::from_ratio(1, 8));
grid.set_closed_box_walls();
let mask = PoissonMask::from_grid(&grid);
for c in 0..4 {
assert_eq!(mask.degree(c), 2, "corner cell of a sealed 2x2x1 slab");
}
let mut open_slab = MacGrid::new(2, 2, 1, Fix128::from_ratio(1, 8));
open_slab.set_w_solid(0, 0, 0, true);
open_slab.set_w_solid(0, 0, 1, true);
let open_mask = PoissonMask::from_grid(&open_slab);
assert_eq!(
open_mask.degree(0),
4,
"only the z faces are walled, the four lateral faces stay open"
);
}
#[test]
fn walls_are_held_at_zero_by_the_projection() {
let mut grid = MacGrid::new(4, 4, 4, Fix128::from_ratio(1, 8));
grid.u.fill(Fix128::ONE);
grid.set_closed_box_walls();
project_pressure(&mut grid, Fix128::from_ratio(1, 100), Fix128::ONE, 50);
for j in 0..4 {
for k in 0..4 {
assert_eq!(grid.u(0, j, k), Fix128::ZERO);
assert_eq!(grid.u(4, j, k), Fix128::ZERO);
}
}
}
const SLAB_CASES: [(usize, usize); 9] = [
(8, 1),
(8, 2),
(8, 4),
(8, 8),
(7, 2),
(7, 3),
(7, 4),
(5, 4),
(3, 4),
];
fn seed_slab_case(n: usize, walls: bool, initial_pressure: bool) -> MacGrid {
let mut grid = seed_divergent_flow(n);
if walls {
grid.set_closed_box_walls();
}
if initial_pressure {
for (c, slot) in grid.pressure.iter_mut().enumerate() {
*slot = Fix128::from_ratio((c % 5) as i64 - 2, 7);
}
}
grid
}
const UNWRITTEN: Fix128 = Fix128::from_int(1_000_000);
fn assemble_slabs(
template: &MacGrid,
faces: &[SlabFaces],
transport: &LocalSlabTransport,
) -> MacGrid {
let mut out = template.clone();
out.pressure.fill(UNWRITTEN);
out.u.fill(UNWRITTEN);
out.v.fill(UNWRITTEN);
out.w.fill(UNWRITTEN);
let (nx, ny) = (out.nx, out.ny);
let plane = cell_plane(nx, ny);
for (r, rank_faces) in faces.iter().enumerate() {
let (k0, k1) = rank_faces.owned();
for k in k0..k1 {
let layer = transport
.slab(r)
.layer(k)
.expect("a rank's own layer is resident in its band");
out.pressure[k * plane..(k + 1) * plane].copy_from_slice(layer);
let (u, _) = rank_faces.u_layer(k);
for j in 0..ny {
for i in 0..=nx {
let ix = out.idx_u(i, j, k);
out.u[ix] = u[i + (nx + 1) * j];
}
}
let (v, _) = rank_faces.v_layer(k);
for j in 0..=ny {
for i in 0..nx {
let ix = out.idx_v(i, j, k);
out.v[ix] = v[i + nx * j];
}
}
}
for k in rank_faces.w_written() {
let (w, _) = rank_faces.w_layer(k);
for j in 0..ny {
for i in 0..nx {
let ix = out.idx_w(i, j, k);
out.w[ix] = w[i + nx * j];
}
}
}
}
out
}
fn solve_slab_local(
base: &MacGrid,
dt_s: Fix128,
density: Fix128,
iterations: u32,
ranks: usize,
halo: usize,
schedule: HaloSchedule,
) -> MacGrid {
let mut enforced = base.clone();
enforced.enforce_face_boundaries();
let nz = enforced.nz;
let plane = cell_plane(enforced.nx, enforced.ny);
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let mut faces: Vec<SlabFaces> = bounds
.iter()
.map(|&b| SlabFaces::from_grid(&enforced, b))
.collect();
let mut transport =
LocalSlabTransport::from_field(&bounds, nz, plane, halo, &enforced.pressure);
project_pressure_slab_local_over(
&mut faces,
dt_s,
density,
iterations,
schedule,
&mut transport,
);
assemble_slabs(base, &faces, &transport)
}
#[test]
fn the_banded_projection_imposes_the_face_conditions_itself() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
for &(n, ranks) in &SLAB_CASES {
let base = seed_slab_scene(n, SlabScene::ClosedBox);
assert!(
base.u.iter().any(|v| !v.is_zero()),
"the seed must carry velocity on the faces the walls will zero"
);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let mut banded = base.clone();
project_pressure_banded(&mut banded, dt, rho, 6, ranks);
assert!(
grids_are_bit_equal(&monolithic, &banded),
"{n}³ over {ranks} ranks: the banded projection on an unenforced sealed box \
differs from the monolithic solve",
);
}
}
#[test]
fn slab_local_storage_reproduces_the_monolithic_pressure_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
for &(n, ranks) in &SLAB_CASES {
for &(walls, initial_pressure) in
&[(false, false), (true, false), (true, true), (false, true)]
{
let base = seed_slab_case(n, walls, initial_pressure);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let assembled =
solve_slab_local(&base, dt, rho, 6, ranks, 1, HaloSchedule::EverySweep);
assert!(
grids_are_bit_equal(&monolithic, &assembled),
"{n}³ over {ranks} slab-local ranks (walls {walls}, initial pressure \
{initial_pressure}) did not reproduce the monolithic solve: either a \
localised array is indexed wrongly, a face is written by the wrong rank, \
or the halo arrives with the wrong contents",
);
}
}
}
#[test]
fn a_stale_slab_halo_does_not_reproduce_the_monolithic_solve() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_slab_case(8, true, true);
let mut monolithic = base.clone();
project_pressure_red_black_gs(&mut monolithic, dt, rho, 6);
let stale = solve_slab_local(&base, dt, rho, 6, 4, 1, HaloSchedule::EveryIteration);
assert!(
!grids_are_bit_equal(&monolithic, &stale),
"a halo one sweep out of date still reproduced the monolithic solve, so the \
bit-equality test above is not actually testing the exchange",
);
}
#[test]
#[should_panic(expected = "slab halo too narrow")]
fn a_halo_narrower_than_the_stencil_aborts_instead_of_returning_a_number() {
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let base = seed_slab_case(8, false, false);
let _ = solve_slab_local(&base, dt, rho, 6, 2, 0, HaloSchedule::EverySweep);
}
#[test]
#[should_panic(expected = "no room for the layer it was asked to receive")]
fn a_slab_delivery_into_a_rank_without_room_aborts() {
let nz = 8usize;
let plane = 4usize;
let bounds: Vec<(usize, usize)> = (0..2).map(|r| slab_bounds(nz, 2, r)).collect();
let field = vec![Fix128::ZERO; nz * plane];
let mut transport = LocalSlabTransport::from_field(&bounds, nz, plane, 0, &field);
exchange_slab_halos_local(&mut transport, &bounds, nz);
}
struct SlabDeliveryRecorder {
calls: Vec<(usize, usize, usize)>,
}
impl SlabTransport for SlabDeliveryRecorder {
fn slab_mut(&mut self, rank: usize) -> &mut SlabStorage {
unreachable!("a schedule walk asked the recorder for rank {rank}'s slab")
}
fn deliver_layer(&mut self, src: usize, dst: usize, layer: usize) {
self.calls.push((src, dst, layer));
}
}
#[test]
fn both_exchange_walks_deliver_the_same_layers_in_the_same_order() {
for ranks in 1usize..=6 {
for nz in [ranks, 2 * ranks + 1, 3 * ranks, 1] {
let bounds: Vec<(usize, usize)> =
(0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let (full_length, _) = recorded_schedule(nz, ranks);
let mut local = SlabDeliveryRecorder { calls: Vec::new() };
exchange_slab_halos_local(&mut local, &bounds, nz);
assert_eq!(
local.calls, full_length,
"{nz} layers over {ranks} ranks: the slab-local exchange walks a \
different sequence from the full-length one",
);
}
}
}
#[test]
fn every_face_layer_is_updated_by_exactly_one_rank() {
for nz in 1usize..=16 {
for ranks in 1usize..=20 {
let mut cells = vec![0usize; nz];
let mut z_faces = vec![0usize; nz + 1];
for r in 0..ranks {
let b = slab_bounds(nz, ranks, r);
let faces = SlabFaces::new(2, 2, nz, Fix128::ONE, b);
for owned in &mut cells[b.0..b.1] {
*owned += 1;
}
for k in faces.w_written() {
z_faces[k] += 1;
}
}
assert!(
cells.iter().all(|&c| c == 1),
"{nz} layers over {ranks} ranks: cell layers swept {cells:?} times",
);
assert!(
z_faces.iter().all(|&c| c == 1),
"{nz} layers over {ranks} ranks: Z-face layers updated {z_faces:?} times",
);
}
}
}
#[test]
fn the_measured_working_set_is_the_one_the_layer_counts_imply() {
let n = 6usize;
let ranks = 3usize;
let grid = seed_slab_case(n, true, false);
let plane = cell_plane(n, n);
let scale = Fix128::ONE;
for r in 0..ranks {
let b = slab_bounds(n, ranks, r);
let faces = SlabFaces::from_grid(&grid, b);
let storage = SlabStorage::for_slab(plane, n, b, 1);
let stencil = SlabStencil::build(&faces, scale);
let bytes = storage.bytes().add(faces.bytes()).add(stencil.bytes());
let (lo, hi) = storage.resident();
let owned = b.1 - b.0;
assert_eq!(
bytes.pressure,
(hi - lo) * plane * size_of::<Fix128>(),
"rank {r}: the pressure band is not the resident layers",
);
assert_eq!(
bytes.faces,
(owned * (u_plane(n, n) + v_plane(n, n)) + (owned + 1) * plane)
* size_of::<Fix128>(),
"rank {r}: the face arrays are not the owned layers plus one Z-face layer",
);
assert_eq!(
bytes.open,
owned * plane * size_of::<[bool; 6]>(),
"rank {r}: the mask is not the owned layers",
);
assert_eq!(bytes.rhs, owned * plane * size_of::<Fix128>());
assert_eq!(bytes.inverse_degrees, owned * plane * size_of::<Fix128>());
assert!(
bytes.total() > 0,
"rank {r} reported an empty working set, which would make every \
comparison against it vacuous",
);
}
}
#[test]
fn the_slab_local_working_set_is_a_fraction_of_the_full_length_one() {
let n = 64usize;
let ranks = 8usize;
let grid = seed_slab_case(n, true, false);
let cells = n * n * n;
let plane = cell_plane(n, n);
let mask = PoissonMask::from_grid(&grid);
let inv_deg = inverse_degrees(&mask, cells);
let rhs = poisson_rhs(&grid, Fix128::ONE);
let full_length_transport = LocalTransport::new(1, n, plane);
let full_length = (grid.u.capacity()
+ grid.v.capacity()
+ grid.w.capacity()
+ grid.pressure.capacity()
+ inv_deg.capacity()
+ rhs.capacity()
+ full_length_transport.slabs[0].capacity())
* size_of::<Fix128>()
+ mask.open.capacity() * size_of::<[bool; 6]>();
let b = slab_bounds(n, ranks, ranks / 2);
let faces = SlabFaces::from_grid(&grid, b);
let storage = SlabStorage::for_slab(plane, n, b, 1);
let stencil = SlabStencil::build(&faces, Fix128::ONE);
let slab_local = storage
.bytes()
.add(faces.bytes())
.add(stencil.bytes())
.total();
assert!(
slab_local * 4 < full_length,
"{n}³ over {ranks} ranks: slab-local storage is {slab_local} B per rank against \
{full_length} B for the full-length path, which is less than the 4x the \
decomposition is for",
);
}
#[cfg(feature = "std")]
#[test]
#[ignore = "allocates over 1 GiB and runs for tens of seconds: the 1e8-cell residency measurement"]
fn one_rank_of_a_hundred_million_cells_fits_in_slab_local_storage() {
let n = 464usize; let ranks = 8usize;
let rank = 4usize; let dx = Fix128::from_ratio(1, 100);
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let plane = cell_plane(n, n);
let b = slab_bounds(n, ranks, rank);
let (k0, k1) = b;
let t0 = std::time::Instant::now();
let mut faces = SlabFaces::new(n, n, n, dx, b);
for k in k0..k1 {
let base = k * u_plane(n, n);
let (u, _) = faces.u_layer_mut(k);
for (i, slot) in u.iter_mut().enumerate() {
*slot = probe_u(base + i);
}
}
let mut storage = SlabStorage::for_slab(plane, n, b, 1);
let seeded = t0.elapsed();
let t1 = std::time::Instant::now();
let stencil = SlabStencil::build(&faces, rho * dx * dx / dt);
let stencil_built = t1.elapsed();
let bytes = storage.bytes().add(faces.bytes()).add(stencil.bytes());
let t2 = std::time::Instant::now();
for colour in 0..2u32 {
sweep_slab_colour(&mut storage, &faces, &stencil, colour);
}
let swept = t2.elapsed();
let owned = (k1 - k0) * plane;
let gib = |b: usize| b as f64 / (1024.0 * 1024.0 * 1024.0);
println!(
"rank {rank} of {ranks}, {n}³ = {} cells, owns layers {k0}..{k1} ({owned} cells)",
n * n * n,
);
println!(
" pressure band {:.3} GiB / faces {:.3} GiB / face conditions {:.3} GiB",
gib(bytes.pressure),
gib(bytes.faces),
gib(bytes.face_flags),
);
println!(
" mask {:.3} GiB / inverse degrees {:.3} GiB / right-hand side {:.3} GiB",
gib(bytes.open),
gib(bytes.inverse_degrees),
gib(bytes.rhs),
);
println!(
" total per rank {:.3} GiB, times {ranks} ranks = {:.3} GiB",
gib(bytes.total()),
gib(bytes.total() * ranks),
);
println!(
" seed {seeded:?} / stencil {stencil_built:?} / one iteration (both colours) \
{swept:?} = {:.1} ns per owned cell",
swept.as_secs_f64() * 1.0e9 / owned as f64,
);
assert!(
bytes.total() < 2 * 1024 * 1024 * 1024,
"one rank of {n}³ over {ranks} ranks needs {:.3} GiB, which is not the \
gigabyte-class residency the decomposition is for",
gib(bytes.total()),
);
}
fn probe_u(global_face: usize) -> Fix128 {
Fix128::from_f64(((global_face % 7) as f64 - 3.0) * 0.1)
}
#[cfg(feature = "std")]
#[test]
#[ignore = "allocates about 2 GiB twice over: the multi-rank timing measurement"]
fn the_whole_slab_decomposition_costs_what_the_monolithic_solve_does_per_cell() {
let n = 256usize;
let ranks = 8usize;
let dx = Fix128::from_ratio(1, 100);
let dt = Fix128::from_ratio(1, 100);
let rho = Fix128::from_int(1000);
let plane = cell_plane(n, n);
let cells = n * n * n;
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(n, ranks, r)).collect();
let gib = |b: usize| b as f64 / (1024.0 * 1024.0 * 1024.0);
let monolithic = {
let mut grid = MacGrid::new(n, n, n, dx);
for (i, slot) in grid.u.iter_mut().enumerate() {
*slot = probe_u(i);
}
let t = std::time::Instant::now();
project_pressure_red_black_gs(&mut grid, dt, rho, 1);
t.elapsed()
};
let t0 = std::time::Instant::now();
let mut faces: Vec<SlabFaces> = bounds
.iter()
.map(|&b| {
let mut f = SlabFaces::new(n, n, n, dx, b);
for k in b.0..b.1 {
let base = k * u_plane(n, n);
let (u, _) = f.u_layer_mut(k);
for (i, slot) in u.iter_mut().enumerate() {
*slot = probe_u(base + i);
}
}
f
})
.collect();
let mut transport = LocalSlabTransport {
slabs: bounds
.iter()
.map(|&b| SlabStorage::for_slab(plane, n, b, 1))
.collect(),
staging: vec![Fix128::ZERO; plane],
};
let seeded = t0.elapsed();
let t1 = std::time::Instant::now();
project_pressure_slab_local_over(
&mut faces,
dt,
rho,
1,
HaloSchedule::EverySweep,
&mut transport,
);
let solved = t1.elapsed();
let per_rank: Vec<usize> = faces
.iter()
.enumerate()
.map(|(r, f)| {
transport
.slab(r)
.bytes()
.add(f.bytes())
.add(SlabStencil::build(f, rho * dx * dx / dt).bytes())
.total()
})
.collect();
let total: usize = per_rank.iter().sum();
println!(
"{n}³ = {cells} cells, one red-black iteration from the same field:\n \
monolithic {monolithic:?} = {:.1} ns/cell\n \
{ranks} slab-local {solved:?} = {:.1} ns/cell (seed {seeded:?})\n \
working set {:.3} GiB over {ranks} ranks, largest rank {:.3} GiB",
monolithic.as_secs_f64() * 1.0e9 / cells as f64,
solved.as_secs_f64() * 1.0e9 / cells as f64,
gib(total),
gib(per_rank.iter().copied().max().unwrap_or(0)),
);
assert!(
solved.as_secs_f64() < 4.0 * monolithic.as_secs_f64(),
"the slab-local decomposition took {solved:?} against the monolithic solve's \
{monolithic:?}, so localising the storage has moved the cost into the index \
arithmetic instead of leaving it in the memory system",
);
}
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum SlabScene {
Open,
ClosedBox,
ZOutflow,
}
impl SlabScene {
fn as_str(self) -> &'static str {
match self {
Self::Open => "open",
Self::ClosedBox => "closed",
Self::ZOutflow => "zoutflow",
}
}
fn parse(name: &str) -> Self {
match name {
"open" => Self::Open,
"closed" => Self::ClosedBox,
"zoutflow" => Self::ZOutflow,
other => panic!("unknown slab scene `{other}`"),
}
}
}
fn seed_slab_scene(n: usize, scene: SlabScene) -> MacGrid {
let mut grid = seed_divergent_flow(n);
let varied_pressure = |grid: &mut MacGrid| {
for (c, slot) in grid.pressure.iter_mut().enumerate() {
*slot = Fix128::from_ratio((c % 5) as i64 - 2, 7);
}
};
match scene {
SlabScene::Open => {}
SlabScene::ClosedBox => {
grid.set_closed_box_walls();
varied_pressure(&mut grid);
}
SlabScene::ZOutflow => {
for k in 0..=n {
for j in 0..n {
for i in 0..n {
let ix = grid.idx_w(i, j, k);
grid.w[ix] =
Fix128::from_ratio(((k * n + i + 2 * j) % 11) as i64 - 5, 4);
}
}
}
varied_pressure(&mut grid);
for j in 0..n {
for i in 0..n {
grid.set_w_bc(i, j, 0, FaceBc::SlipWall);
grid.set_w_bc(i, j, n, FaceBc::SlipWall);
for k in (1..n).filter(|k| k % 2 == 1) {
grid.set_w_bc(i, j, k, FaceBc::Outflow);
}
}
}
for k in 0..n {
for j in 0..n {
grid.set_u_bc(
0,
j,
k,
FaceBc::Inflow {
normal_velocity: Fix128::from_ratio(3, 4),
},
);
grid.set_u_bc(n, j, k, FaceBc::Outflow);
}
for i in 0..n {
grid.set_v_bc(
i,
0,
k,
FaceBc::Wall {
velocity: Vec3Fix::ZERO,
},
);
grid.set_v_bc(
i,
n,
k,
FaceBc::Wall {
velocity: Vec3Fix::ZERO,
},
);
}
}
}
}
grid
}
fn z_outflow_crossings(nz: usize, ranks: usize) -> usize {
(0..ranks)
.map(|r| slab_bounds(nz, ranks, r))
.filter(|&(k0, k1)| k0 != k1 && k0 > 0 && k0 < nz && k0 % 2 == 1)
.count()
}
fn band_for_rank(
plane: usize,
nz: usize,
b: (usize, usize),
halo: usize,
field: &[Fix128],
) -> SlabStorage {
let mut slab = SlabStorage::for_slab(plane, nz, b, halo);
let (lo, hi) = slab.resident();
for k in lo..hi {
slab.layer_mut(k)
.expect("a layer inside the band this slab just reported")
.copy_from_slice(&field[k * plane..(k + 1) * plane]);
}
slab
}
fn assemble_slab_faces(template: &MacGrid, faces: &[SlabFaces]) -> MacGrid {
let mut out = template.clone();
out.u.fill(UNWRITTEN);
out.v.fill(UNWRITTEN);
out.w.fill(UNWRITTEN);
let (nx, ny) = (out.nx, out.ny);
for rank_faces in faces {
let (k0, k1) = rank_faces.owned();
for k in k0..k1 {
let (u, _) = rank_faces.u_layer(k);
for j in 0..ny {
for i in 0..=nx {
let ix = out.idx_u(i, j, k);
out.u[ix] = u[i + (nx + 1) * j];
}
}
let (v, _) = rank_faces.v_layer(k);
for j in 0..=ny {
for i in 0..nx {
let ix = out.idx_v(i, j, k);
out.v[ix] = v[i + nx * j];
}
}
}
for k in rank_faces.w_written() {
let (w, _) = rank_faces.w_layer(k);
for j in 0..ny {
for i in 0..nx {
let ix = out.idx_w(i, j, k);
out.w[ix] = w[i + nx * j];
}
}
}
}
out
}
fn enforce_in_slabs(base: &MacGrid, ranks: usize) -> Vec<SlabFaces> {
let nz = base.nz;
let bounds: Vec<(usize, usize)> = (0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let mut faces: Vec<SlabFaces> = bounds
.iter()
.map(|&b| SlabFaces::from_grid(base, b))
.collect();
let cond: Vec<SlabFaceConditions> = bounds
.iter()
.map(|&b| SlabFaceConditions::from_grid(base, b))
.collect();
enforce_slab_face_boundaries_over(&mut faces, &cond);
faces
}
#[test]
fn distributed_face_enforcement_matches_the_monolithic_one() {
let mut crossings = 0usize;
for &(n, ranks) in &SLAB_CASES {
crossings += z_outflow_crossings(n, ranks);
for scene in [SlabScene::Open, SlabScene::ClosedBox, SlabScene::ZOutflow] {
let base = seed_slab_scene(n, scene);
let mut monolithic = base.clone();
monolithic.enforce_face_boundaries();
let faces = enforce_in_slabs(&base, ranks);
let assembled = assemble_slab_faces(&base, &faces);
assert_eq!(
assembled.u, monolithic.u,
"{n}³ over {ranks} ranks ({scene:?}): the X-face conditions the ranks \
imposed differ from the monolithic ones",
);
assert_eq!(
assembled.v, monolithic.v,
"{n}³ over {ranks} ranks ({scene:?}): the Y-face conditions the ranks \
imposed differ from the monolithic ones",
);
assert_eq!(
assembled.w, monolithic.w,
"{n}³ over {ranks} ranks ({scene:?}): the Z-face conditions the ranks \
imposed differ from the monolithic ones — the layer crossing a slab \
boundary is the one an outflow face reads",
);
}
}
assert!(
crossings > 0,
"no rank of any case had a first Z-face that reads across a slab boundary, so \
the crossing this stage adds was never exercised",
);
}
#[test]
fn both_ranks_sharing_a_z_face_layer_enforce_it_to_the_same_value() {
let mut shared = 0usize;
for &(n, ranks) in &SLAB_CASES {
let base = seed_slab_scene(n, SlabScene::ZOutflow);
let mut monolithic = base.clone();
monolithic.enforce_face_boundaries();
for (r, rank_faces) in enforce_in_slabs(&base, ranks).iter().enumerate() {
let (k0, k1) = rank_faces.owned();
if k0 == k1 || k1 == n {
continue; }
let (w, _) = rank_faces.w_layer(k1);
for j in 0..n {
for i in 0..n {
assert_eq!(
w[i + n * j],
monolithic.w[monolithic.idx_w(i, j, k1)],
"{n}³ over {ranks} ranks: rank {r} recomputed the conditions on \
Z-face layer {k1}, which rank above it writes, and got a \
different value at ({i}, {j})",
);
}
}
shared += 1;
}
}
assert!(
shared > 0,
"no rank held a Z-face layer it does not write, so the recomputation this test \
exists to check never happened",
);
}
#[test]
fn a_z_outflow_on_the_domain_boundary_alone_needs_no_face_halo() {
let mut ran = 0usize;
let mut outflow_faces = 0usize;
for &(n, ranks) in &SLAB_CASES {
let mut base = seed_divergent_flow(n);
for j in 0..n {
for i in 0..n {
base.set_w_bc(i, j, 0, FaceBc::Outflow);
base.set_w_bc(i, j, n, FaceBc::Outflow);
}
}
let mut monolithic = base.clone();
monolithic.enforce_face_boundaries();
let bounds: Vec<(usize, usize)> =
(0..ranks).map(|r| slab_bounds(n, ranks, r)).collect();
let mut faces: Vec<SlabFaces> = bounds
.iter()
.map(|&b| SlabFaces::from_grid(&base, b))
.collect();
for (r, rank_faces) in faces.iter_mut().enumerate() {
let cond = SlabFaceConditions::from_grid(&base, bounds[r]);
outflow_faces += cond
.w
.iter()
.filter(|&&(_, bc)| bc == FaceBc::Outflow)
.count();
enforce_slab_face_boundaries(rank_faces, &cond, None);
ran += 1;
}
let assembled = assemble_slab_faces(&base, &faces);
assert_eq!(
assembled.w, monolithic.w,
"{n}³ over {ranks} ranks: a Z-outflow on the domain boundary alone did not \
reproduce the monolithic enforcement without a halo",
);
}
let expected_ranks: usize = SLAB_CASES.iter().map(|&(_, ranks)| ranks).sum();
assert_eq!(
ran, expected_ranks,
"only {ran} of {expected_ranks} ranks ran, so the ones that did not are not \
evidence of anything",
);
assert!(
outflow_faces > 0,
"the ranks held no outflow Z-face between them, so running without a halo \
proves nothing",
);
#[cfg(feature = "std")]
println!(
"boundary-only Z-outflow: {ran} ranks over {} cases imposed {outflow_faces} \
outflow Z-faces with nothing handed over",
SLAB_CASES.len(),
);
}
#[test]
#[should_panic(expected = "which the rank below has to hand over")]
fn an_interior_z_outflow_is_what_makes_the_face_halo_load_bearing() {
let n = 8usize;
let base = seed_slab_scene(n, SlabScene::ZOutflow);
let b = slab_bounds(n, 8, 1);
assert_eq!(
b,
(1, 2),
"the rank this test aims at owns a different layer"
);
let mut faces = SlabFaces::from_grid(&base, b);
let cond = SlabFaceConditions::from_grid(&base, b);
enforce_slab_face_boundaries(&mut faces, &cond, None);
}
fn seed_low_end_outflow(n: usize) -> MacGrid {
let mut grid = seed_divergent_flow(n);
for (c, slot) in grid.u.iter_mut().enumerate() {
*slot = Fix128::from_ratio(((c * 7) % 13) as i64 - 6, 4);
}
for (c, slot) in grid.v.iter_mut().enumerate() {
*slot = Fix128::from_ratio(((c * 5) % 11) as i64 - 5, 4);
}
for (c, slot) in grid.w.iter_mut().enumerate() {
*slot = Fix128::from_ratio(((c * 3) % 17) as i64 - 8, 4);
}
let rest = FaceBc::Wall {
velocity: Vec3Fix::ZERO,
};
for k in 0..n {
for j in 0..n {
grid.set_u_bc(0, j, k, FaceBc::Outflow);
let inward = if j % 2 == 0 {
FaceBc::Inflow {
normal_velocity: Fix128::from_ratio(3, 4),
}
} else {
rest
};
grid.set_u_bc(1, j, k, inward);
grid.set_u_bc(n, j, k, FaceBc::Outflow);
}
for i in 0..n {
grid.set_v_bc(i, 0, k, FaceBc::Outflow);
let inward = if i % 2 == 0 {
FaceBc::Wall {
velocity: Vec3Fix::new(
Fix128::from_ratio(1, 2),
Fix128::ZERO,
Fix128::ZERO,
),
}
} else {
FaceBc::SlipWall
};
grid.set_v_bc(i, 1, k, inward);
grid.set_v_bc(i, n, k, FaceBc::Outflow);
}
}
for j in 0..n {
for i in 0..n {
grid.set_w_bc(i, j, 0, FaceBc::Outflow);
let inward = if (i + j) % 2 == 0 {
FaceBc::Inflow {
normal_velocity: Fix128::from_ratio(-5, 8),
}
} else {
rest
};
grid.set_w_bc(i, j, 1, inward);
grid.set_w_bc(i, j, n, FaceBc::Outflow);
}
}
grid
}
#[test]
fn low_end_outflow_next_to_inflow_or_wall_matches_the_monolithic_enforcement() {
let mut exercised = 0usize;
#[cfg(feature = "std")]
let mut rows = Vec::new();
let mut differing = 0usize;
for &(n, ranks) in &SLAB_CASES {
if n < 2 {
continue;
}
let base = seed_low_end_outflow(n);
let mut monolithic = base.clone();
monolithic.enforce_face_boundaries();
for k in 0..n {
for j in 0..n {
if base.u(1, j, k) != monolithic.u(1, j, k) {
exercised += 1;
}
}
}
let assembled = assemble_slab_faces(&base, &enforce_in_slabs(&base, ranks));
let count =
|a: &[Fix128], b: &[Fix128]| a.iter().zip(b).filter(|(x, y)| x != y).count();
let (du, dv, dw) = (
count(&assembled.u, &monolithic.u),
count(&assembled.v, &monolithic.v),
count(&assembled.w, &monolithic.w),
);
differing += du + dv + dw;
#[cfg(feature = "std")]
rows.push(format!("| {n}³ | {ranks} | {du} | {dv} | {dw} |"));
}
#[cfg(feature = "std")]
println!(
"faces differing from the monolithic enforcement:\n| grid | ranks | X | Y | Z |\n|---|---|---|---|---|\n{}",
rows.join("\n"),
);
assert!(
exercised > 0,
"no low-end inward neighbour changed under enforcement, so the order is not tested",
);
assert_eq!(
differing, 0,
"the slab-local enforcement disagrees with the monolithic one next to a low-end \
outflow face",
);
}
#[test]
fn conditions_named_face_by_face_match_the_ones_read_from_a_grid() {
let n = 6usize;
let ranks = 3usize;
let base = seed_slab_scene(n, SlabScene::ZOutflow);
for r in 0..ranks {
let b = slab_bounds(n, ranks, r);
let from_grid = SlabFaceConditions::from_grid(&base, b);
let mut by_hand = SlabFaceConditions::new(n, n, b);
if b.0 != b.1 {
for k in b.0..b.1 {
for j in 0..n {
for i in 0..=n {
by_hand.set_u(i, j, k, base.u_bc(i, j, k));
}
}
for j in 0..=n {
for i in 0..n {
by_hand.set_v(i, j, k, base.v_bc(i, j, k));
}
}
}
for k in b.0..=b.1 {
for j in 0..n {
for i in 0..n {
by_hand.set_w(i, j, k, base.w_bc(i, j, k));
}
}
}
}
assert_eq!(
from_grid.u, by_hand.u,
"rank {r}: the X-face conditions differ"
);
assert_eq!(
from_grid.v, by_hand.v,
"rank {r}: the Y-face conditions differ"
);
assert_eq!(
from_grid.w, by_hand.w,
"rank {r}: the Z-face conditions differ"
);
assert!(
!from_grid.w.is_empty() || b.0 == b.1,
"rank {r} read no Z-face conditions out of a scene that sets them on every \
layer, so the comparison above is vacuous",
);
}
}
struct PlaneRecorder {
plane: usize,
sent: Vec<(usize, usize)>,
received: Vec<(usize, usize)>,
}
impl PlaneChannel for PlaneRecorder {
fn send_plane(&mut self, dst: usize, layer: usize, plane: &[Fix128]) {
assert_eq!(
plane.len(),
self.plane,
"a plane of the wrong size was offered"
);
self.sent.push((dst, layer));
}
fn recv_plane(&mut self, src: usize, layer: usize, plane: &mut [Fix128]) {
assert_eq!(
plane.len(),
self.plane,
"a plane of the wrong size was asked for"
);
self.received.push((src, layer));
}
}
#[test]
fn the_face_condition_pipeline_carries_one_plane_per_slab_boundary() {
for ranks in 1usize..=6 {
for nz in [ranks, 2 * ranks + 1, 3 * ranks, 1] {
let bounds: Vec<(usize, usize)> =
(0..ranks).map(|r| slab_bounds(nz, ranks, r)).collect();
let busy: Vec<usize> = (0..ranks).filter(|&r| bounds[r].0 != bounds[r].1).collect();
let expected: Vec<(usize, usize, usize)> = busy
.windows(2)
.map(|pair| (pair[0], pair[1], bounds[pair[0]].1 - 1))
.collect();
let plane = cell_plane(2, 2);
let mut recorder = PlaneRecorder {
plane,
sent: Vec::new(),
received: Vec::new(),
};
for (r, &b) in bounds.iter().enumerate() {
let mut faces = SlabFaces::new(2, 2, nz, Fix128::ONE, b);
let cond = SlabFaceConditions::new(2, 2, b);
enforce_slab_face_boundaries_on_rank(
&mut faces,
&cond,
ranks,
r,
&mut recorder,
);
}
assert_eq!(
recorder.sent.len(),
recorder.received.len(),
"{nz} layers over {ranks} ranks: {} planes were offered and {} asked for",
recorder.sent.len(),
recorder.received.len(),
);
let walked: Vec<(usize, usize, usize)> = recorder
.received
.iter()
.zip(&recorder.sent)
.map(|(&(src, layer), &(dst, sent_layer))| {
assert_eq!(
layer, sent_layer,
"{nz} layers over {ranks} ranks: layer {sent_layer} was offered \
and layer {layer} asked for",
);
(src, dst, layer)
})
.collect();
assert_eq!(
walked, expected,
"{nz} layers over {ranks} ranks: the pipeline carried a different \
sequence of planes than the partition implies",
);
}
}
}
#[cfg(feature = "std")]
#[test]
fn the_two_wire_encodings_agree() {
use std::io::Cursor;
let plane = 5usize;
let values: Vec<Fix128> = (0..plane)
.map(|c| Fix128::from_ratio(c as i64 * 7 - 11, 3))
.collect();
let mut by_hand = vec![0u8; plane * WIRE_BYTES_PER_CELL];
encode_plane(&values, &mut by_hand);
let mut transport: SocketTransport<Cursor<Vec<u8>>> = SocketTransport::new(
0,
plane * 2,
plane,
vec![None, Some(Cursor::new(Vec::new()))],
);
transport.slab_mut(0)[plane..].copy_from_slice(&values);
transport.deliver_layer(0, 1, 1);
let on_the_wire = transport.links[1]
.as_ref()
.expect("the link the delivery went out on")
.get_ref()
.clone();
assert_eq!(
on_the_wire, by_hand,
"the slab transport's encoding differs from the full-length transport's, so the \
two cannot be peers on one wire",
);
assert!(
!by_hand.iter().all(|&b| b == 0),
"the encoding produced only zero bytes, so the comparison above is vacuous",
);
let mut back = vec![Fix128::ZERO; plane];
decode_plane(&by_hand, &mut back);
assert_eq!(back, values, "decoding is not the inverse of encoding");
}
#[cfg(feature = "std")]
const XSLAB_RANK: &str = "ALICE_PHYSICS_XSLAB_RANK";
#[cfg(feature = "std")]
const XSLAB_RANKS: &str = "ALICE_PHYSICS_XSLAB_RANKS";
#[cfg(feature = "std")]
const XSLAB_BROKER: &str = "ALICE_PHYSICS_XSLAB_BROKER";
#[cfg(feature = "std")]
const XSLAB_SIZE: &str = "ALICE_PHYSICS_XSLAB_SIZE";
#[cfg(feature = "std")]
const XSLAB_SCENE: &str = "ALICE_PHYSICS_XSLAB_SCENE";
#[cfg(feature = "std")]
const XSLAB_FAULT: &str = "ALICE_PHYSICS_XSLAB_FAULT";
#[cfg(feature = "std")]
const XSLAB_KIND: &str = "ALICE_PHYSICS_XSLAB_KIND";
#[cfg(feature = "std")]
const XSLAB_ITERS: &str = "ALICE_PHYSICS_XSLAB_ITERS";
#[cfg(feature = "std")]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum SlabRunKind {
MonolithicReference,
Banded,
}
#[cfg(feature = "std")]
impl SlabRunKind {
fn as_str(self) -> &'static str {
match self {
Self::MonolithicReference => "monolithic-reference",
Self::Banded => "banded",
}
}
fn parse(name: &str) -> Self {
match name {
"monolithic-reference" => Self::MonolithicReference,
"banded" => Self::Banded,
other => panic!("unknown slab run kind `{other}`"),
}
}
}
#[cfg(feature = "std")]
const XSLAB_CHILD: &str = "eulerian_grid::tests::cross_process_slab_rank_child";
#[cfg(feature = "std")]
const XSLAB_ITERATIONS: u32 = 6;
#[cfg(feature = "std")]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum SlabCrossFault {
None,
NoPressureHalo,
BlankFaceHalo,
}
#[cfg(feature = "std")]
impl SlabCrossFault {
fn as_str(self) -> &'static str {
match self {
Self::None => "none",
Self::NoPressureHalo => "nohalo",
Self::BlankFaceHalo => "blankface",
}
}
fn parse(name: &str) -> Self {
match name {
"none" => Self::None,
"nohalo" => Self::NoPressureHalo,
"blankface" => Self::BlankFaceHalo,
other => panic!("unknown slab cross-process fault `{other}`"),
}
}
}
#[cfg(feature = "std")]
struct DroppedHalo<T> {
inner: T,
}
#[cfg(feature = "std")]
impl<T: SlabTransport> SlabTransport for DroppedHalo<T> {
fn slab_mut(&mut self, rank: usize) -> &mut SlabStorage {
self.inner.slab_mut(rank)
}
fn deliver_layer(&mut self, _src: usize, _dst: usize, _layer: usize) {}
}
#[cfg(feature = "std")]
struct BlankedPlane<C> {
inner: C,
scratch: Vec<Fix128>,
}
#[cfg(feature = "std")]
impl<C: PlaneChannel> PlaneChannel for BlankedPlane<C> {
fn send_plane(&mut self, dst: usize, layer: usize, plane: &[Fix128]) {
self.inner.send_plane(dst, layer, plane);
}
fn recv_plane(&mut self, src: usize, layer: usize, plane: &mut [Fix128]) {
self.inner.recv_plane(src, layer, &mut self.scratch);
plane.fill(Fix128::ZERO);
}
}
#[cfg(feature = "std")]
fn compare_band(reference: &MacGrid, faces: &SlabFaces, slab: &SlabStorage) -> [u32; 4] {
let (nx, ny) = (reference.nx, reference.ny);
let plane = cell_plane(nx, ny);
let (k0, k1) = faces.owned();
let mut bad_cells = 0u32;
let mut bad_faces = 0u32;
let mut first_bad = u32::MAX;
let mut compared = 0u32;
let mut note = |layer: usize| {
first_bad = first_bad.min(layer as u32);
};
for k in k0..k1 {
let mine = slab.layer(k).expect("a rank's own layer is resident");
let theirs = &reference.pressure[k * plane..(k + 1) * plane];
for (a, b) in mine.iter().zip(theirs) {
compared += 1;
if a != b {
bad_cells += 1;
note(k);
}
}
let (u, _) = faces.u_layer(k);
for j in 0..ny {
for i in 0..=nx {
if u[i + (nx + 1) * j] != reference.u[reference.idx_u(i, j, k)] {
bad_faces += 1;
note(k);
}
}
}
let (v, _) = faces.v_layer(k);
for j in 0..=ny {
for i in 0..nx {
if v[i + nx * j] != reference.v[reference.idx_v(i, j, k)] {
bad_faces += 1;
note(k);
}
}
}
}
for k in faces.w_written() {
let (w, _) = faces.w_layer(k);
for j in 0..ny {
for i in 0..nx {
if w[i + nx * j] != reference.w[reference.idx_w(i, j, k)] {
bad_faces += 1;
note(k);
}
}
}
}
[bad_cells, bad_faces, first_bad, compared]
}
#[cfg(feature = "std")]
fn solve_slab_as_rank(
n: usize,
ranks: usize,
my_rank: usize,
scene: SlabScene,
fault: SlabCrossFault,
links: Vec<Option<std::net::TcpStream>>,
) -> [u32; 4] {
let (dt, rho) = (xproc_dt(), xproc_rho());
let base = seed_slab_scene(n, scene);
let b = slab_bounds(n, ranks, my_rank);
let plane = cell_plane(n, n);
let mut faces = SlabFaces::from_grid(&base, b);
let cond = SlabFaceConditions::from_grid(&base, b);
let band = band_for_rank(plane, n, b, 1, &base.pressure);
let mut socket = SlabSocketTransport::new(my_rank, plane, band, links);
if fault == SlabCrossFault::BlankFaceHalo {
let mut channel = BlankedPlane {
inner: socket,
scratch: vec![Fix128::ZERO; plane],
};
enforce_slab_face_boundaries_on_rank(&mut faces, &cond, ranks, my_rank, &mut channel);
socket = channel.inner;
} else {
enforce_slab_face_boundaries_on_rank(&mut faces, &cond, ranks, my_rank, &mut socket);
}
if fault == SlabCrossFault::NoPressureHalo {
let mut dropped = DroppedHalo { inner: socket };
project_pressure_slab_local_on_rank(
&mut faces,
dt,
rho,
XSLAB_ITERATIONS,
ranks,
HaloSchedule::EverySweep,
my_rank,
&mut dropped,
);
socket = dropped.inner;
} else {
project_pressure_slab_local_on_rank(
&mut faces,
dt,
rho,
XSLAB_ITERATIONS,
ranks,
HaloSchedule::EverySweep,
my_rank,
&mut socket,
);
}
let mut reference = base.clone();
project_pressure_red_black_gs(&mut reference, dt, rho, XSLAB_ITERATIONS);
compare_band(&reference, &faces, socket.slab())
}
#[cfg(feature = "std")]
fn join_slab_mesh(
my_rank: usize,
ranks: usize,
ports: &[u16],
listener: &std::net::TcpListener,
) -> Vec<Option<std::net::TcpStream>> {
use std::io::{Read, Write};
use std::net::TcpStream;
let mut links: Vec<Option<TcpStream>> = (0..ranks).map(|_| None).collect();
for lower in 0..my_rank {
let mut link = TcpStream::connect(("127.0.0.1", ports[lower]))
.unwrap_or_else(|e| panic!("rank {my_rank} dialling rank {lower}: {e}"));
bound_peer_link(&link);
link.write_all(&(my_rank as u16).to_le_bytes())
.expect("announce which rank is dialling");
link.flush().expect("flush the rank announcement");
links[lower] = Some(link);
}
for _ in my_rank + 1..ranks {
let mut link = accept_before_deadline(listener, "a higher rank");
let mut who = [0u8; 2];
link.read_exact(&mut who)
.expect("learn which rank dialled in");
let peer = usize::from(u16::from_le_bytes(who));
assert!(
peer > my_rank && peer < ranks && links[peer].is_none(),
"rank {my_rank} was dialled by rank {peer}, which is not a higher rank it \
is still waiting for",
);
links[peer] = Some(link);
}
links
}
#[cfg(feature = "std")]
#[test]
fn cross_process_slab_rank_child() {
use std::io::{Read, Write};
use std::net::{TcpListener, TcpStream};
let Ok(rank) = std::env::var(XSLAB_RANK) else {
return;
};
let my_rank: usize = rank.parse().expect("the rank is a number");
let ranks: usize = std::env::var(XSLAB_RANKS)
.expect("rank count from rank 0")
.parse()
.expect("the rank count is a number");
assert!(
my_rank > 0 && my_rank < ranks,
"rank {my_rank} is not a child rank of a {ranks}-rank harness",
);
let n: usize = std::env::var(XSLAB_SIZE)
.expect("grid size from rank 0")
.parse()
.expect("grid size is a number");
let broker: u16 = std::env::var(XSLAB_BROKER)
.expect("broker port from rank 0")
.parse()
.expect("broker port is a number");
let scene = SlabScene::parse(&std::env::var(XSLAB_SCENE).expect("scene from rank 0"));
let fault = SlabCrossFault::parse(&std::env::var(XSLAB_FAULT).expect("fault from rank 0"));
let listener = TcpListener::bind(("127.0.0.1", 0)).expect("bind a loopback port");
let my_port = listener
.local_addr()
.expect("the bound loopback address")
.port();
let mut to_root = TcpStream::connect(("127.0.0.1", broker))
.expect("connect to rank 0 on the port it brokers on");
bound_peer_link(&to_root);
let mut hello = [0u8; 4];
hello[..2].copy_from_slice(&(my_rank as u16).to_le_bytes());
hello[2..].copy_from_slice(&my_port.to_le_bytes());
to_root.write_all(&hello).expect("announce rank and port");
to_root.flush().expect("flush the announcement");
let mut table = vec![0u8; 2 * ranks];
to_root
.read_exact(&mut table)
.expect("learn the port of every rank from rank 0");
let ports: Vec<u16> = table
.chunks_exact(2)
.map(|pair| u16::from_le_bytes([pair[0], pair[1]]))
.collect();
let links = join_slab_mesh(my_rank, ranks, &ports, &listener);
let report: [u64; 4] =
match SlabRunKind::parse(&std::env::var(XSLAB_KIND).expect("run kind from rank 0")) {
SlabRunKind::MonolithicReference => {
let counts = solve_slab_as_rank(n, ranks, my_rank, scene, fault, links);
let mut wide = [0u64; 4];
for (slot, value) in wide.iter_mut().zip(counts) {
*slot = u64::from(value);
}
wide
}
SlabRunKind::Banded => {
let iterations: u32 = std::env::var(XSLAB_ITERS)
.expect("iteration count from rank 0")
.parse()
.expect("the iteration count is a number");
solve_slab_banded_as_rank(n, ranks, my_rank, scene, fault, iterations, links)
}
};
let mut bytes = [0u8; 32];
for (slot, value) in bytes.chunks_exact_mut(8).zip(report) {
slot.copy_from_slice(&value.to_le_bytes());
}
to_root.write_all(&bytes).expect("report to rank 0");
to_root.flush().expect("flush the report");
}
#[cfg(feature = "std")]
struct SlabRanks {
kids: Vec<std::process::Child>,
}
#[cfg(feature = "std")]
impl Drop for SlabRanks {
fn drop(&mut self) {
reap_child_ranks(&mut self.kids);
}
}
#[cfg(feature = "std")]
impl SlabRanks {
fn accept(&mut self, listener: &std::net::TcpListener) -> std::net::TcpStream {
use std::io::ErrorKind;
use std::time::{Duration, Instant};
listener
.set_nonblocking(true)
.expect("poll for a child rank instead of blocking on it");
let deadline = Instant::now() + Duration::from_secs(60);
let link = loop {
match listener.accept() {
Ok((stream, _)) => break stream,
Err(e) if e.kind() == ErrorKind::WouldBlock => {}
Err(e) => panic!("accepting a child rank failed: {e}"),
}
for (slot, kid) in self.kids.iter_mut().enumerate() {
if let Some(status) = kid.try_wait().expect("poll a child rank") {
panic!(
"rank {} exited ({status}) before connecting: does \
`{XSLAB_CHILD}` still name a test?",
slot + 1,
);
}
}
assert!(
Instant::now() < deadline,
"a child rank did not connect within 60 s",
);
std::thread::sleep(Duration::from_millis(2));
};
link.set_nonblocking(false)
.expect("switch the accepted stream back to blocking");
bound_peer_link(&link);
link
}
fn join(&mut self) {
for (slot, kid) in self.kids.iter_mut().enumerate() {
let status = kid.wait().expect("wait for a child rank to exit");
assert!(status.success(), "rank {} exited with {status}", slot + 1);
}
}
}
#[cfg(feature = "std")]
fn open_slab_mesh(
ranks: usize,
env: &[(&'static str, String)],
) -> (
SlabRanks,
Vec<Option<std::net::TcpStream>>,
Vec<Option<std::net::TcpStream>>,
) {
use std::io::{Read, Write};
use std::net::TcpListener;
use std::process::{Command, Stdio};
assert!(ranks >= 2, "a cross-process case needs at least two ranks");
let broker = TcpListener::bind(("127.0.0.1", 0)).expect("bind the broker port");
let broker_port = broker
.local_addr()
.expect("the bound loopback address")
.port();
let mine = TcpListener::bind(("127.0.0.1", 0)).expect("bind rank 0's data port");
let my_port = mine.local_addr().expect("the bound address").port();
let exe = std::env::current_exe().expect("path of this test binary");
let mut kids = SlabRanks {
kids: (1..ranks)
.map(|rank| {
let mut cmd = Command::new(&exe);
cmd.args(["--exact", XSLAB_CHILD, "--test-threads=1"])
.env(XSLAB_RANK, rank.to_string())
.env(XSLAB_RANKS, ranks.to_string())
.env(XSLAB_BROKER, broker_port.to_string())
.stdin(Stdio::null())
.stdout(Stdio::null())
.stderr(Stdio::null());
for &(name, ref value) in env {
cmd.env(name, value);
}
cmd.spawn().unwrap_or_else(|e| {
panic!("re-execute this test binary as rank {rank}: {e}")
})
})
.collect(),
};
let mut brokered: Vec<Option<std::net::TcpStream>> = (0..ranks).map(|_| None).collect();
let mut ports = vec![0u16; ranks];
ports[0] = my_port;
for _ in 1..ranks {
let mut link = kids.accept(&broker);
let mut hello = [0u8; 4];
link.read_exact(&mut hello)
.expect("a child rank's announcement");
let rank = usize::from(u16::from_le_bytes([hello[0], hello[1]]));
assert!(
rank > 0 && rank < ranks && brokered[rank].is_none(),
"rank {rank} announced itself twice or is not a child rank",
);
ports[rank] = u16::from_le_bytes([hello[2], hello[3]]);
brokered[rank] = Some(link);
}
let table: Vec<u8> = ports.iter().flat_map(|p| p.to_le_bytes()).collect();
for link in brokered.iter_mut().flatten() {
link.write_all(&table).expect("hand a child the port table");
link.flush().expect("flush the port table");
}
let links = join_slab_mesh(0, ranks, &ports, &mine);
(kids, brokered, links)
}
#[cfg(feature = "std")]
fn collect_slab_reports(brokered: &mut [Option<std::net::TcpStream>]) -> Vec<[u64; 4]> {
use std::io::Read;
let mut reports = vec![[0u64; 4]; brokered.len()];
for (rank, link) in brokered.iter_mut().enumerate() {
let Some(link) = link else { continue };
let mut bytes = [0u8; 32];
link.read_exact(&mut bytes)
.unwrap_or_else(|e| panic!("rank {rank}'s report: {e}"));
for (slot, chunk) in reports[rank].iter_mut().zip(bytes.chunks_exact(8)) {
*slot = u64::from_le_bytes(chunk.try_into().expect("8 bytes of a report"));
}
}
reports
}
#[cfg(feature = "std")]
fn slab_reports_across_processes(
n: usize,
ranks: usize,
scene: SlabScene,
fault: SlabCrossFault,
) -> Vec<[u32; 4]> {
let (mut kids, mut brokered, links) = open_slab_mesh(
ranks,
&[
(XSLAB_KIND, SlabRunKind::MonolithicReference.as_str().into()),
(XSLAB_SIZE, n.to_string()),
(XSLAB_SCENE, scene.as_str().into()),
(XSLAB_FAULT, fault.as_str().into()),
],
);
let mine = solve_slab_as_rank(n, ranks, 0, scene, fault, links);
kids.join();
let mut reports = collect_slab_reports(&mut brokered);
let narrow = |wide: [u64; 4]| -> [u32; 4] {
let mut out = [0u32; 4];
for (slot, value) in out.iter_mut().zip(wide) {
*slot = u32::try_from(value).expect("a report field of a comparison fits in u32");
}
out
};
let mut out: Vec<[u32; 4]> = reports.drain(..).map(narrow).collect();
out[0] = mine;
out
}
#[cfg(feature = "std")]
fn assert_slab_processes_agree(n: usize, ranks: usize, scene: SlabScene) {
let reports = slab_reports_across_processes(n, ranks, scene, SlabCrossFault::None);
for (rank, report) in reports.iter().enumerate() {
assert_eq!(
*report,
[0, 0, u32::MAX, report[3]],
"{n}³ over {ranks} processes ({scene:?}): rank {rank} found {} pressure \
cells and {} face values differing from the single-process solve, first at \
layer {}",
report[0],
report[1],
report[2],
);
}
let compared: u32 = reports.iter().map(|report| report[3]).sum();
assert_eq!(
compared as usize,
n * n * n,
"{n}³ over {ranks} processes ({scene:?}): the ranks compared {compared} cells \
between them instead of the whole domain, so the agreement above is partly \
vacuous",
);
}
#[cfg(feature = "std")]
#[test]
fn three_processes_reproduce_the_monolithic_slab_local_solve() {
for &(n, scene) in &[
(9usize, SlabScene::ClosedBox),
(9, SlabScene::ZOutflow),
(8, SlabScene::ZOutflow),
] {
if scene == SlabScene::ZOutflow {
assert!(
z_outflow_crossings(n, 3) > 0,
"{n}³ over 3 ranks has no Z-face layer crossing a slab boundary, so this \
case does not exercise the face-condition pipeline",
);
}
assert_slab_processes_agree(n, 3, scene);
}
}
#[cfg(feature = "std")]
#[test]
fn four_processes_one_of_which_owns_nothing_reproduce_the_monolithic_slab_local_solve() {
let bounds: Vec<(usize, usize)> = (0..4).map(|r| slab_bounds(3, 4, r)).collect();
assert!(
bounds.iter().any(|&(k0, k1)| k0 == k1),
"this case was chosen because a rank owns nothing and none does: {bounds:?}",
);
assert!(
z_outflow_crossings(3, 4) > 0,
"3³ over 4 ranks has no Z-face layer crossing a slab boundary",
);
assert_slab_processes_agree(3, 4, SlabScene::ZOutflow);
}
#[cfg(feature = "std")]
#[test]
fn a_cross_process_slab_run_without_the_pressure_halo_does_not_reproduce_the_solve() {
let reports = slab_reports_across_processes(
8,
3,
SlabScene::ZOutflow,
SlabCrossFault::NoPressureHalo,
);
assert!(
reports.iter().any(|r| r[0] > 0),
"three processes that exchanged no halo at all still reproduced the \
single-process pressure field, so the cross-process oracle is not measuring the \
exchange: {reports:?}",
);
}
#[cfg(feature = "std")]
#[test]
fn a_cross_process_slab_run_without_the_z_face_layer_does_not_reproduce_the_solve() {
assert!(
z_outflow_crossings(8, 3) > 0,
"8³ over 3 ranks has no Z-face layer crossing a slab boundary, so blanking it \
could not change anything",
);
let reports =
slab_reports_across_processes(8, 3, SlabScene::ZOutflow, SlabCrossFault::BlankFaceHalo);
assert!(
reports.iter().any(|r| r[0] > 0 || r[1] > 0),
"three processes that imposed their face conditions without the layer below \
still reproduced the single-process solve, so the face-condition pipeline is \
not measured: {reports:?}",
);
}
#[cfg(feature = "std")]
#[test]
#[should_panic(expected = "the driver is not rank-local")]
fn a_slab_transport_refuses_to_hand_out_another_ranks_band() {
use std::io::Cursor;
let n = 4usize;
let plane = cell_plane(n, n);
let b = slab_bounds(n, 2, 0);
let band = SlabStorage::for_slab(plane, n, b, 1);
let links: Vec<Option<Cursor<Vec<u8>>>> = vec![None, Some(Cursor::new(Vec::new()))];
let mut transport = SlabSocketTransport::new(0, plane, band, links);
let _ = transport.slab_mut(1);
}
const FOLD_BASIS: u64 = 0xcbf2_9ce4_8422_2325;
const FOLD_PRIME: u64 = 0x0000_0100_0000_01b3;
fn fold_word(acc: u64, word: u64) -> u64 {
let mut out = acc;
for byte in word.to_le_bytes() {
out = (out ^ u64::from(byte)).wrapping_mul(FOLD_PRIME);
}
out
}
fn fold_fix(acc: u64, value: Fix128) -> u64 {
fold_word(fold_word(acc, value.hi as u64), value.lo)
}
fn fold_band(faces: &SlabFaces, slab: &SlabStorage) -> [u64; 4] {
let (k0, k1) = faces.owned();
let mut cells = FOLD_BASIS;
let mut values = FOLD_BASIS;
let mut cell_count = 0u64;
let mut face_count = 0u64;
for k in k0..k1 {
let mine = slab.layer(k).expect("a rank's own layer is resident");
for &p in mine {
cells = fold_fix(cells, p);
cell_count += 1;
}
let (u, _) = faces.u_layer(k);
for &x in u {
values = fold_fix(values, x);
face_count += 1;
}
let (v, _) = faces.v_layer(k);
for &x in v {
values = fold_fix(values, x);
face_count += 1;
}
}
for k in faces.w_written() {
let (w, _) = faces.w_layer(k);
for &x in w {
values = fold_fix(values, x);
face_count += 1;
}
}
[cells, values, cell_count, face_count]
}
#[test]
fn the_fold_separates_values_that_differ_in_their_lowest_bit() {
let base = Fix128 { hi: 3, lo: 1 << 40 };
for (what, other) in [
(
"the lowest fractional bit",
Fix128 {
hi: base.hi,
lo: base.lo | 1,
},
),
(
"the highest fractional bit",
Fix128 {
hi: base.hi,
lo: base.lo | 1 << 63,
},
),
(
"the integer word",
Fix128 {
hi: base.hi + 1,
lo: base.lo,
},
),
] {
assert_ne!(
fold_fix(FOLD_BASIS, base),
fold_fix(FOLD_BASIS, other),
"two values differing in {what} ({base:?} and {other:?}) fold to the same \
number, so a distributed run that differs there would be reported as \
bit-exact",
);
}
}
fn scene_u_bc(n: usize, scene: SlabScene, i: usize) -> FaceBc {
match scene {
SlabScene::Open => FaceBc::Fluid,
SlabScene::ClosedBox if i == 0 || i == n => box_wall(n),
SlabScene::ClosedBox => FaceBc::Fluid,
SlabScene::ZOutflow if i == 0 => FaceBc::Inflow {
normal_velocity: Fix128::from_ratio(3, 4),
},
SlabScene::ZOutflow if i == n => FaceBc::Outflow,
SlabScene::ZOutflow => FaceBc::Fluid,
}
}
fn scene_v_bc(n: usize, scene: SlabScene, j: usize) -> FaceBc {
match scene {
SlabScene::Open => FaceBc::Fluid,
SlabScene::ClosedBox if j == 0 || j == n => box_wall(n),
SlabScene::ZOutflow if j == 0 || j == n => FaceBc::Wall {
velocity: Vec3Fix::ZERO,
},
SlabScene::ClosedBox | SlabScene::ZOutflow => FaceBc::Fluid,
}
}
fn scene_w_bc(n: usize, scene: SlabScene, k: usize) -> FaceBc {
match scene {
SlabScene::Open => FaceBc::Fluid,
SlabScene::ClosedBox if k == 0 || k == n => box_wall(n),
SlabScene::ClosedBox => FaceBc::Fluid,
SlabScene::ZOutflow if k == 0 || k == n => FaceBc::SlipWall,
SlabScene::ZOutflow if k % 2 == 1 => FaceBc::Outflow,
SlabScene::ZOutflow => FaceBc::Fluid,
}
}
fn box_wall(n: usize) -> FaceBc {
if n == 1 {
FaceBc::SlipWall
} else {
FaceBc::Wall {
velocity: Vec3Fix::ZERO,
}
}
}
fn scene_flags(bc: FaceBc) -> FaceFlags {
FaceFlags {
solid: bc.is_wall(),
blocks_pressure: bc.blocks_pressure(),
inflow: bc.is_inflow(),
}
}
fn scene_u_velocity(n: usize, global_face: usize) -> Fix128 {
Fix128::from_int((global_face % (n + 1)) as i64)
}
fn scene_w_velocity(n: usize, scene: SlabScene, i: usize, j: usize, k: usize) -> Fix128 {
match scene {
SlabScene::Open | SlabScene::ClosedBox => Fix128::ZERO,
SlabScene::ZOutflow => Fix128::from_ratio(((k * n + i + 2 * j) % 11) as i64 - 5, 4),
}
}
fn scene_pressure(scene: SlabScene, global_cell: usize) -> Fix128 {
match scene {
SlabScene::Open => Fix128::ZERO,
SlabScene::ClosedBox | SlabScene::ZOutflow => {
Fix128::from_ratio((global_cell % 5) as i64 - 2, 7)
}
}
}
fn band_for_scene(
n: usize,
scene: SlabScene,
b: (usize, usize),
) -> (SlabFaces, SlabFaceConditions, SlabStorage) {
let (k0, k1) = b;
let plane = cell_plane(n, n);
let mut faces = SlabFaces::new(n, n, n, Fix128::ONE, b);
let mut cond = SlabFaceConditions::new(n, n, b);
let mut band = SlabStorage::for_slab(plane, n, b, 1);
let (lo, hi) = band.resident();
for k in lo..hi {
let layer = band
.layer_mut(k)
.expect("a layer inside the band this slab just reported");
for (c, slot) in layer.iter_mut().enumerate() {
*slot = scene_pressure(scene, k * plane + c);
}
}
if k0 == k1 {
return (faces, cond, band);
}
for k in k0..k1 {
let row = k * u_plane(n, n);
let (u, flags) = faces.u_layer_mut(k);
for j in 0..n {
for i in 0..=n {
let at = i + (n + 1) * j;
u[at] = scene_u_velocity(n, row + at);
flags[at] = scene_flags(scene_u_bc(n, scene, i));
}
}
for j in 0..n {
for i in 0..=n {
let bc = scene_u_bc(n, scene, i);
if bc != FaceBc::Fluid {
cond.set_u(i, j, k, bc);
}
}
}
let (v, flags) = faces.v_layer_mut(k);
for j in 0..=n {
for i in 0..n {
flags[i + n * j] = scene_flags(scene_v_bc(n, scene, j));
}
}
debug_assert!(
v.iter().all(|&x| x == Fix128::ZERO),
"no scene seeds a Y-face velocity, so the band starts at zero like the grid",
);
for j in 0..=n {
for i in 0..n {
let bc = scene_v_bc(n, scene, j);
if bc != FaceBc::Fluid {
cond.set_v(i, j, k, bc);
}
}
}
}
for k in k0..=k1 {
let (w, flags) = faces.w_layer_mut(k);
for j in 0..n {
for i in 0..n {
let at = i + n * j;
w[at] = scene_w_velocity(n, scene, i, j, k);
flags[at] = scene_flags(scene_w_bc(n, scene, k));
}
}
for j in 0..n {
for i in 0..n {
let bc = scene_w_bc(n, scene, k);
if bc != FaceBc::Fluid {
cond.set_w(i, j, k, bc);
}
}
}
}
(faces, cond, band)
}
#[cfg(feature = "std")]
fn solve_slab_banded_as_rank(
n: usize,
ranks: usize,
my_rank: usize,
scene: SlabScene,
fault: SlabCrossFault,
iterations: u32,
links: Vec<Option<std::net::TcpStream>>,
) -> [u64; 4] {
let (dt, rho) = (xproc_dt(), xproc_rho());
let b = slab_bounds(n, ranks, my_rank);
let plane = cell_plane(n, n);
let (mut faces, cond, band) = band_for_scene(n, scene, b);
let mut socket = SlabSocketTransport::new(my_rank, plane, band, links);
if fault == SlabCrossFault::BlankFaceHalo {
let mut channel = BlankedPlane {
inner: socket,
scratch: vec![Fix128::ZERO; plane],
};
enforce_slab_face_boundaries_on_rank(&mut faces, &cond, ranks, my_rank, &mut channel);
socket = channel.inner;
} else {
enforce_slab_face_boundaries_on_rank(&mut faces, &cond, ranks, my_rank, &mut socket);
}
if fault == SlabCrossFault::NoPressureHalo {
let mut dropped = DroppedHalo { inner: socket };
project_pressure_slab_local_on_rank(
&mut faces,
dt,
rho,
iterations,
ranks,
HaloSchedule::EverySweep,
my_rank,
&mut dropped,
);
socket = dropped.inner;
} else {
project_pressure_slab_local_on_rank(
&mut faces,
dt,
rho,
iterations,
ranks,
HaloSchedule::EverySweep,
my_rank,
&mut socket,
);
}
fold_band(&faces, socket.slab())
}
#[cfg(feature = "std")]
fn banded_slab_folds_across_processes(
n: usize,
ranks: usize,
scene: SlabScene,
fault: SlabCrossFault,
iterations: u32,
) -> Vec<[u64; 4]> {
let (mut kids, mut brokered, links) = open_slab_mesh(
ranks,
&[
(XSLAB_KIND, SlabRunKind::Banded.as_str().into()),
(XSLAB_SIZE, n.to_string()),
(XSLAB_SCENE, scene.as_str().into()),
(XSLAB_FAULT, fault.as_str().into()),
(XSLAB_ITERS, iterations.to_string()),
],
);
let mine = solve_slab_banded_as_rank(n, ranks, 0, scene, fault, iterations, links);
kids.join();
let mut folds = collect_slab_reports(&mut brokered);
folds[0] = mine;
folds
}
#[cfg(feature = "std")]
fn fold_reference_band(reference: &MacGrid, n: usize, ranks: usize, rank: usize) -> [u64; 4] {
let b = slab_bounds(n, ranks, rank);
let faces = SlabFaces::from_grid(reference, b);
let band = band_for_rank(cell_plane(n, n), n, b, 1, &reference.pressure);
fold_band(&faces, &band)
}
#[cfg(feature = "std")]
fn assert_banded_slab_processes_agree(
n: usize,
ranks: usize,
scene: SlabScene,
iterations: u32,
) -> (Vec<[u64; 4]>, [u64; 4]) {
let folds =
banded_slab_folds_across_processes(n, ranks, scene, SlabCrossFault::None, iterations);
let mut reference = seed_slab_scene(n, scene);
project_pressure_red_black_gs(&mut reference, xproc_dt(), xproc_rho(), iterations);
let mut cells = 0u64;
let mut faces = 0u64;
for (rank, got) in folds.iter().enumerate() {
let want = fold_reference_band(&reference, n, ranks, rank);
assert_eq!(
*got, want,
"{n}³ over {ranks} processes ({scene:?}, {iterations} iterations): rank {rank} \
folded its own layers to {got:?} against {want:?} for the single-process \
solve, so the distributed run is not bit-exact over them",
);
cells += got[2];
faces += got[3];
}
assert_eq!(
cells,
(n * n * n) as u64,
"{n}³ over {ranks} processes ({scene:?}): the ranks folded {cells} pressure cells \
between them instead of the whole domain, so the agreement above is partly vacuous",
);
let want_faces =
(n * u_plane(n, n) + n * v_plane(n, n) + (n + 1) * cell_plane(n, n)) as u64;
assert_eq!(
faces, want_faces,
"{n}³ over {ranks} processes ({scene:?}): the ranks folded {faces} face values \
between them instead of the {want_faces} the domain has, so the velocity half of \
the agreement above is over a subset",
);
let whole = fold_reference_band(&reference, n, 1, 0);
assert_eq!(
whole[2],
(n * n * n) as u64,
"the whole-domain fold covered {} cells instead of {}",
whole[2],
n * n * n,
);
(folds, whole)
}
#[test]
fn the_band_local_scene_matches_the_one_read_from_a_grid() {
let mut sparse_entries = 0usize;
let mut solid_faces = 0usize;
let mut nonzero_pressure = 0usize;
let mut inflow_faces = 0usize;
let mut cases = 0usize;
for &(n, ranks) in SLAB_CASES.iter().chain(&[(1, 2), (1, 3), (2, 2), (2, 5)]) {
for scene in [SlabScene::Open, SlabScene::ClosedBox, SlabScene::ZOutflow] {
let grid = seed_slab_scene(n, scene);
let plane = cell_plane(n, n);
for rank in 0..ranks {
cases += 1;
let b = slab_bounds(n, ranks, rank);
let (faces, cond, band) = band_for_scene(n, scene, b);
let want_faces = SlabFaces::from_grid(&grid, b);
let want_cond = SlabFaceConditions::from_grid(&grid, b);
let want_band = band_for_rank(plane, n, b, 1, &grid.pressure);
let where_ = format!("{n}³ over {ranks} ranks ({scene:?}), rank {rank}");
assert_eq!(faces.u, want_faces.u, "{where_}: X-face velocities");
assert_eq!(faces.v, want_faces.v, "{where_}: Y-face velocities");
assert_eq!(faces.w, want_faces.w, "{where_}: Z-face velocities");
assert_eq!(faces.u_flags, want_faces.u_flags, "{where_}: X-face flags");
assert_eq!(faces.v_flags, want_faces.v_flags, "{where_}: Y-face flags");
assert_eq!(faces.w_flags, want_faces.w_flags, "{where_}: Z-face flags");
assert_eq!(cond.u, want_cond.u, "{where_}: X-face conditions");
assert_eq!(cond.v, want_cond.v, "{where_}: Y-face conditions");
assert_eq!(cond.w, want_cond.w, "{where_}: Z-face conditions");
let (lo, hi) = band.resident();
assert_eq!(
(lo, hi),
want_band.resident(),
"{where_}: resident pressure layers",
);
for k in lo..hi {
assert_eq!(
band.layer(k),
want_band.layer(k),
"{where_}: initial pressure of layer {k}",
);
nonzero_pressure += band
.layer(k)
.expect("a resident layer")
.iter()
.filter(|&&p| p != Fix128::ZERO)
.count();
}
sparse_entries += cond.u.len() + cond.v.len() + cond.w.len();
solid_faces += faces
.u_flags
.iter()
.chain(&faces.v_flags)
.chain(&faces.w_flags)
.filter(|f| f.blocks_pressure)
.count();
inflow_faces += faces.u_flags.iter().filter(|f| f.inflow).count();
}
}
}
assert!(cases > 0, "no case ran");
assert!(
sparse_entries > 0,
"every case compared an empty condition list, so the sparse half of the \
construction is untested",
);
assert!(
solid_faces > 0,
"no case produced a face that blocks pressure, so the flags are all defaults",
);
assert!(
inflow_faces > 0,
"no case produced an inflow face, so the one predicate the sparse map carries and \
the dense flag cannot is untested",
);
assert!(
nonzero_pressure > 0,
"every case started from zero pressure, so a band that forgot to seed its halo \
would compare equal",
);
}
#[cfg(feature = "std")]
#[test]
fn the_banded_cross_process_slab_path_reproduces_the_monolithic_solve() {
let cases = [
(9usize, 3usize, SlabScene::ClosedBox, XSLAB_ITERATIONS),
(8, 3, SlabScene::ZOutflow, 1),
(3, 4, SlabScene::ZOutflow, XSLAB_ITERATIONS + 3),
(6, 2, SlabScene::Open, XSLAB_ITERATIONS),
];
assert!(
cases
.iter()
.any(|&(.., iterations)| iterations != XSLAB_ITERATIONS),
"every case runs {XSLAB_ITERATIONS} iterations, so a rank that ignored the count it \
was handed would pass this test",
);
for &(n, ranks, scene, iterations) in &cases {
if scene == SlabScene::ZOutflow {
assert!(
z_outflow_crossings(n, ranks) > 0,
"{n}³ over {ranks} ranks has no Z-face layer crossing a slab boundary, so \
this case does not exercise the face-condition pipeline",
);
}
assert_banded_slab_processes_agree(n, ranks, scene, iterations);
}
}
#[cfg(feature = "std")]
#[test]
fn a_banded_cross_process_slab_run_without_the_pressure_halo_does_not_reproduce_the_solve() {
let (n, ranks) = (8usize, 3usize);
let folds = banded_slab_folds_across_processes(
n,
ranks,
SlabScene::ZOutflow,
SlabCrossFault::NoPressureHalo,
XSLAB_ITERATIONS,
);
let mut reference = seed_slab_scene(n, SlabScene::ZOutflow);
project_pressure_red_black_gs(&mut reference, xproc_dt(), xproc_rho(), XSLAB_ITERATIONS);
let differing = folds
.iter()
.enumerate()
.filter(|&(rank, got)| *got != fold_reference_band(&reference, n, ranks, rank))
.count();
assert!(
differing > 0,
"{n}³ over {ranks} processes that exchanged no pressure halo at all still folded \
to the single-process solve, so the banded oracle is not measuring the exchange: \
{folds:?}",
);
}
#[cfg(feature = "std")]
#[test]
fn a_banded_cross_process_slab_run_without_the_z_face_layer_does_not_reproduce_the_solve() {
let (n, ranks) = (8usize, 3usize);
assert!(
z_outflow_crossings(n, ranks) > 0,
"{n}³ over {ranks} ranks has no Z-face layer crossing a slab boundary, so blanking \
it could not change anything",
);
let folds = banded_slab_folds_across_processes(
n,
ranks,
SlabScene::ZOutflow,
SlabCrossFault::BlankFaceHalo,
XSLAB_ITERATIONS,
);
let mut reference = seed_slab_scene(n, SlabScene::ZOutflow);
project_pressure_red_black_gs(&mut reference, xproc_dt(), xproc_rho(), XSLAB_ITERATIONS);
let differing = folds
.iter()
.enumerate()
.filter(|&(rank, got)| *got != fold_reference_band(&reference, n, ranks, rank))
.count();
assert!(
differing > 0,
"{n}³ over {ranks} processes that imposed their face conditions without the layer \
below still folded to the single-process solve, so the face-condition pipeline is \
not measured: {folds:?}",
);
}
#[cfg(feature = "std")]
const XSLAB_BIG_SIZE: &str = "ALICE_PHYSICS_XSLAB_BIG_SIZE";
#[cfg(feature = "std")]
const XSLAB_BIG_RANKS: &str = "ALICE_PHYSICS_XSLAB_BIG_RANKS";
#[cfg(feature = "std")]
const XSLAB_BIG_SCENE: &str = "ALICE_PHYSICS_XSLAB_BIG_SCENE";
#[cfg(feature = "std")]
const XSLAB_BIG_ITERATIONS: &str = "ALICE_PHYSICS_XSLAB_BIG_ITERATIONS";
#[cfg(feature = "std")]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
struct BigSlabRun {
n: usize,
ranks: usize,
scene: SlabScene,
iterations: u32,
}
#[cfg(feature = "std")]
fn parse_big_slab_run(read: &dyn Fn(&str) -> Option<String>) -> Result<BigSlabRun, String> {
let number = |name: &str| -> Result<usize, String> {
let raw = read(name).ok_or_else(|| {
format!(
"{name} is not set: the large cross-process run takes its size from the \
environment and has no default, because a default would make a mistyped \
variable name look like a successful run at the size you asked for",
)
})?;
let value: usize = raw
.trim()
.parse()
.map_err(|_| format!("{name} is `{raw}`, which is not a number"))?;
Ok(value)
};
let n = number(XSLAB_BIG_SIZE)?;
if n == 0 {
return Err(format!(
"{XSLAB_BIG_SIZE} is 0, so there is no domain to solve"
));
}
let ranks = number(XSLAB_BIG_RANKS)?;
if ranks < 2 {
return Err(format!(
"{XSLAB_BIG_RANKS} is {ranks}: a cross-process run needs at least two \
processes, and the point of this one is more than two",
));
}
let scene_name = read(XSLAB_BIG_SCENE).ok_or_else(|| {
format!("{XSLAB_BIG_SCENE} is not set: name `open`, `closed` or `zoutflow`")
})?;
let scene = match scene_name.trim() {
"open" => SlabScene::Open,
"closed" => SlabScene::ClosedBox,
"zoutflow" => SlabScene::ZOutflow,
other => {
return Err(format!(
"{XSLAB_BIG_SCENE} is `{other}`, which is not `open`, `closed` or `zoutflow`",
));
}
};
let iterations = u32::try_from(number(XSLAB_BIG_ITERATIONS)?)
.map_err(|_| format!("{XSLAB_BIG_ITERATIONS} does not fit in a u32"))?;
if iterations == 0 {
return Err(format!(
"{XSLAB_BIG_ITERATIONS} is 0, so neither side would solve anything and the two \
would agree about the field they started from",
));
}
Ok(BigSlabRun {
n,
ranks,
scene,
iterations,
})
}
#[cfg(feature = "std")]
#[test]
fn a_mistyped_large_run_request_is_refused_rather_than_defaulted() {
let set = |pairs: &[(&str, &str)]| -> Vec<(String, String)> {
pairs
.iter()
.map(|&(k, v)| (k.to_owned(), v.to_owned()))
.collect()
};
let reader = |pairs: Vec<(String, String)>| {
move |name: &str| -> Option<String> {
pairs
.iter()
.find(|(k, _)| k == name)
.map(|(_, v)| v.clone())
}
};
let good = [
(XSLAB_BIG_SIZE, "512"),
(XSLAB_BIG_RANKS, "8"),
(XSLAB_BIG_SCENE, "closed"),
(XSLAB_BIG_ITERATIONS, "6"),
];
let with = |name: &str, value: &str| -> Vec<(String, String)> {
let mut pairs = set(&good);
for entry in &mut pairs {
if entry.0 == name {
entry.1 = value.to_owned();
}
}
pairs
};
let without = |name: &str| -> Vec<(String, String)> {
set(&good).into_iter().filter(|(k, _)| k != name).collect()
};
assert_eq!(
parse_big_slab_run(&reader(set(&good))),
Ok(BigSlabRun {
n: 512,
ranks: 8,
scene: SlabScene::ClosedBox,
iterations: 6,
}),
"a complete request is accepted exactly as given",
);
for (what, pairs, must_name) in [
("no size at all", without(XSLAB_BIG_SIZE), XSLAB_BIG_SIZE),
("a size of zero", with(XSLAB_BIG_SIZE, "0"), XSLAB_BIG_SIZE),
(
"a size that is not a number",
with(XSLAB_BIG_SIZE, "512³"),
XSLAB_BIG_SIZE,
),
("no rank count", without(XSLAB_BIG_RANKS), XSLAB_BIG_RANKS),
("a single rank", with(XSLAB_BIG_RANKS, "1"), XSLAB_BIG_RANKS),
("no scene", without(XSLAB_BIG_SCENE), XSLAB_BIG_SCENE),
(
"an unknown scene",
with(XSLAB_BIG_SCENE, "ClosedBox"),
XSLAB_BIG_SCENE,
),
(
"no iteration count",
without(XSLAB_BIG_ITERATIONS),
XSLAB_BIG_ITERATIONS,
),
(
"zero iterations",
with(XSLAB_BIG_ITERATIONS, "0"),
XSLAB_BIG_ITERATIONS,
),
] {
let got = parse_big_slab_run(&reader(pairs));
let message = match got {
Err(message) => message,
Ok(run) => panic!("{what} was accepted as {run:?} instead of being refused"),
};
assert!(
message.contains(must_name),
"{what} was refused with `{message}`, which does not name {must_name}, so the \
message does not say which variable to fix",
);
}
}
#[cfg(feature = "std")]
#[test]
#[ignore = "measurement: takes its size from the environment and allocates by the \
gigabyte; the path it runs is pinned at test-suite sizes by \
the_banded_cross_process_slab_path_reproduces_the_monolithic_solve"]
fn a_large_banded_cross_process_slab_run_reproduces_the_monolithic_solve() {
let run = match parse_big_slab_run(&|name| std::env::var(name).ok()) {
Ok(run) => run,
Err(why) => panic!("{why}"),
};
let cells = run.n * run.n * run.n;
println!(
"{n}³ = {cells} cells over {ranks} processes ({scene:?}), {iterations} red-black \
iterations, bit-exact against the single-process solve",
n = run.n,
ranks = run.ranks,
scene = run.scene,
iterations = run.iterations,
);
let started = std::time::Instant::now();
let (folds, whole) =
assert_banded_slab_processes_agree(run.n, run.ranks, run.scene, run.iterations);
println!(
" agreed in {:?} ({:.1} ns per cell per iteration, distributed and single-process \
runs together)",
started.elapsed(),
started.elapsed().as_secs_f64() * 1.0e9 / (cells * run.iterations as usize) as f64,
);
println!(
" WHOLE-DOMAIN FOLD n={} pressure={:#018x} faces={:#018x} cells={} face_values={}",
run.n, whole[0], whole[1], whole[2], whole[3]
);
for (rank, fold) in folds.iter().enumerate() {
println!(
" band fold rank={rank}/{} pressure={:#018x} faces={:#018x} cells={} \
face_values={}",
run.ranks, fold[0], fold[1], fold[2], fold[3]
);
}
}
}