diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a73fd88..805a587 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -53,11 +53,10 @@ jobs: - name: Build run: cmake --build build -j"$(nproc)" - # Exercises MPI decomposition, halo exchange and 20 full time steps. - # A crash, an abort or a non-zero exit fails the job. It does NOT check - # that the answer is right -- that is what the regression test adds next. - - name: Smoke test (32^3, 2 ranks) - working-directory: build - run: | - printf '1.0 0.1 1600\n0.05 20 1\n32 32 32\n' > input.in - mpirun -n 2 --oversubscribe ./imexlbm + # Runs a 32^3 Taylor-Green case on 1, 2 and 4 ranks and asserts the two + # exact invariants of the discrete dynamics (mass, momentum) plus + # agreement between rank counts. Subsumes a smoke test: a crash, an abort + # or a non-zero exit fails it too. See the header of the script for what + # each check is for and why kinetic-energy monotonicity is NOT asserted. + - name: Regression test + run: ./tests/regression.sh build/imexlbm diff --git a/src/lbm.cpp b/src/lbm.cpp index 409af00..1a19557 100644 --- a/src/lbm.cpp +++ b/src/lbm.cpp @@ -239,6 +239,72 @@ void LBM::ComputeMacroscopic() LBM_PROF_END(P_MACRO); }; +// --------------------------------------------------------------------------- +// Globally reduced diagnostics. +// +// In this pressure-based formulation sum_i f_i = 3p and sum_i f_i e_i = u, and +// the equilibrium reproduces both moments exactly. With sum_i t_i = 1, +// sum_i t_i e_i = 0 and sum_i t_i e_ia e_ib = delta_ab / 3: +// +// sum_i feq_i = 3p + 4.5*(u^2/3) - 1.5*u^2 = 3p = sum_i f_i +// sum_i feq_i e_i = 3 u_b delta_ab / 3 = u_a = sum_i f_i e_i +// +// so BGK collision leaves both untouched, and streaming on a fully periodic +// domain is a permutation of f. Mass and momentum are therefore invariants of +// the discrete dynamics up to floating-point rounding. A relative drift much +// larger than 1e-12 is a bug -- in streaming, in the halo exchange, or a race -- +// not "numerical error". +// +// ke is NOT conserved: dissipating it is what the Taylor-Green vortex is for. +// It is reported so a test can check that it decays, and that its trajectory +// does not depend on how the domain was split across ranks. +// +// Cost is one 5-double MPI_Allreduce per call. The caller already issues an +// MPI_Barrier at the same cadence, so the marginal cost is in the noise even +// at full machine scale. +// --------------------------------------------------------------------------- +void LBM::Conserved(double &mass, double &mom_x, double &mom_y, double &mom_z, + double &ke) +{ + double lmass = 0.0, lmx = 0.0, lmy = 0.0, lmz = 0.0, lke = 0.0; + + Kokkos::parallel_reduce( + "conserved", + mdrange_policy3({l_s[0], l_s[1], l_s[2]}, {l_e[0], l_e[1], l_e[2]}), + KOKKOS_CLASS_LAMBDA(const int i, const int j, const int k, + double &m, double &mx, double &my, double &mz, + double &e_kin) { + double pl = 0.0, ul = 0.0, vl = 0.0, wl = 0.0; + for (int ii = 0; ii < Q27; ++ii) + { + const double fv = f(ii, i, j, k); + pl += fv; + ul += fv * e(ii, 0); + vl += fv * e(ii, 1); + wl += fv * e(ii, 2); + } + m += pl; + mx += ul; + my += vl; + mz += wl; + e_kin += 0.5 * (ul * ul + vl * vl + wl * wl); + }, + Kokkos::Sum(lmass), Kokkos::Sum(lmx), + Kokkos::Sum(lmy), Kokkos::Sum(lmz), + Kokkos::Sum(lke)); + Kokkos::fence(); + + double local[5] = {lmass, lmx, lmy, lmz, lke}; + double total[5]; + MPI_Allreduce(local, total, 5, MPI_DOUBLE, MPI_SUM, comm); + + mass = total[0]; + mom_x = total[1]; + mom_y = total[2]; + mom_z = total[3]; + ke = total[4]; +}; + void LBM::MPIoutput(int n) { // MPI_IO diff --git a/src/lbm.hpp b/src/lbm.hpp index 85307df..4f2cb3a 100644 --- a/src/lbm.hpp +++ b/src/lbm.hpp @@ -197,6 +197,11 @@ struct LBM void MPIoutput(int n); void Output(int n); + // Globally reduced conserved quantities; see the definition in lbm.cpp for + // why mass and momentum are exact invariants of the discrete dynamics. + void Conserved(double &mass, double &mom_x, double &mom_y, double &mom_z, + double &ke); + void pack_f(buffer_f ff); void unpack_f(buffer_f ff); diff --git a/src/main.cpp b/src/main.cpp index 28fa46f..d83308e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -79,7 +79,23 @@ if (l1.me == 0) { printf("step %6d | interval %8.4f s | total %8.4f s | %10.2f MLUPS\n", it, dt_int, now - start, mlups); } - t_last = now; + + // Collective: every rank must call it. Mass and momentum are + // exact invariants of the discrete dynamics (see LBM::Conserved), + // so a drift here is a bug, not rounding. Deliberately printed on + // its own line with no "MLUPS" token, because tools/parse_scaling.sh + // selects lines with awk '/MLUPS/'. + double mass, px, py, pz, ke; + l1.Conserved(mass, px, py, pz, ke); + if (l1.me == 0) + { + printf("cons step %6d | mass %.15e | px %.15e | py %.15e | pz %.15e | ke %.15e\n", + it, mass, px, py, pz, ke); + } + + // Re-taken after the diagnostic so its cost is not charged to + // the next interval's MLUPS. + t_last = MPI_Wtime(); // l1.MPIoutput(it / s1.inter); } } diff --git a/tests/regression.sh b/tests/regression.sh new file mode 100755 index 0000000..75ea358 --- /dev/null +++ b/tests/regression.sh @@ -0,0 +1,133 @@ +#!/usr/bin/env bash +# +# tests/regression.sh -- physical regression test for IMEXLBM. +# +# usage: tests/regression.sh [path-to-imexlbm] (default ./imexlbm) +# +# What is checked, and why these are the right things to check: +# +# 1. Mass conservation. sum_i f_i is an exact invariant of the discrete +# dynamics -- BGK collision reproduces it and streaming on a periodic +# domain is a permutation (see the comment on LBM::Conserved). Any drift +# beyond rounding is a bug in streaming, in the halo exchange, or a race. +# +# 2. Momentum conservation. Likewise exact, and the Taylor-Green initial +# field has zero net momentum, so it must stay at zero. +# +# 3. Decomposition independence. Collision is per-cell and streaming is a +# permutation, so the field is bitwise independent of how the domain is +# split across ranks; only the order of the final reduction differs. +# Running the same case on 1, 2 and 4 ranks must therefore agree to near +# machine epsilon. This is the most valuable check here: the Cartesian +# decomposition and the 26-neighbour halo exchange are the most intricate +# part of the code, and the part that is hardest to eyeball on a machine +# where you get one job per queue wait. +# +# What is deliberately NOT checked: monotone decay of the kinetic energy. The +# distributions are initialised at equilibrium with an unrelaxed pressure field, +# so the first tens of steps carry an acoustic transient and ke genuinely +# oscillates (9.852 -> 10.219 -> 9.889 at 32^3). Asserting monotonicity would be +# a flaky test that is also wrong about the physics. +# +# Tolerances are ~100x the values measured on Kokkos 5.1.1 OpenMP + Open MPI: +# mass drift within a run 1.2e-14 +# mass across rank counts 1.4e-14 +# kinetic energy across ranks 1.0e-14 +# If a tolerance ever has to be loosened, find out why first -- on these +# invariants, a growing residual is evidence, not noise. + +set -euo pipefail + +TOL_MASS_DRIFT=1e-12 # relative +TOL_MOMENTUM=1e-10 # absolute; the exact initial value is ~1e-14 +TOL_RANK_AGREE=1e-12 # relative +RANKS="1 2 4" + +# Resolve the executable before changing directory. +raw=${1:-./imexlbm} +EXE="$(cd "$(dirname "$raw")" && pwd)/$(basename "$raw")" +if [ ! -x "$EXE" ]; then + echo "FAIL no executable at $EXE" >&2 + exit 1 +fi + +workdir=$(mktemp -d) +trap 'rm -rf "$workdir"' EXIT +cd "$workdir" + +# 32^3 for 20 steps, reporting every step. Small enough to run in seconds on a +# CI runner, large enough that the decomposition is non-trivial at 4 ranks. +printf '1.0 0.1 1600\n0.05 20 1\n32 32 32\n' > input.in + +echo "running $EXE on ${RANKS// /, } rank(s)" +for n in $RANKS; do + if ! OMP_NUM_THREADS=1 OMP_PROC_BIND=false \ + mpirun -n "$n" --oversubscribe "$EXE" > "run.$n.log" 2>&1; then + echo "FAIL solver exited non-zero on $n rank(s)" >&2 + tail -20 "run.$n.log" >&2 + exit 1 + fi + grep '^cons' "run.$n.log" > "cons.$n" || true + if [ ! -s "cons.$n" ]; then + echo "FAIL no 'cons' diagnostic lines on $n rank(s)" >&2 + echo " the binary predates LBM::Conserved, or the run produced no output" >&2 + exit 1 + fi + if grep -qiE 'nan|inf' "cons.$n"; then + echo "FAIL non-finite value in the diagnostics on $n rank(s)" >&2 + grep -iE 'nan|inf' "cons.$n" | head -3 >&2 + exit 1 + fi +done +echo + +status=0 + +# --- 1 and 2: invariants hold within each run ------------------------------ +# Field layout of a 'cons' line under FS=[ |]+ : +# 1 cons 2 step 3 4 mass 5 6 px 7 +# 8 py 9 10 pz 11 12 ke 13 +for n in $RANKS; do + awk -F'[ |]+' -v n="$n" -v tm="$TOL_MASS_DRIFT" -v tp="$TOL_MOMENTUM" ' + NR == 1 { m0 = $5 } + { + d = ($5 - m0) / m0; if (d < 0) d = -d; if (d > dm) dm = d + for (i = 7; i <= 11; i += 2) { p = $i; if (p < 0) p = -p; if (p > dp) dp = p } + } + END { + ok = 1 + if (dm > tm) { printf "FAIL mass not conserved on %s rank(s): max relative drift %.3e > %s\n", n, dm, tm; ok = 0 } + else { printf "ok mass conserved on %s rank(s) (max drift %.3e)\n", n, dm } + if (dp > tp) { printf "FAIL momentum not conserved on %s rank(s): max |p| %.3e > %s\n", n, dp, tp; ok = 0 } + else { printf "ok net momentum stays zero on %s rank(s) (max |p| %.3e)\n", n, dp } + exit ok ? 0 : 1 + }' "cons.$n" || status=1 +done +echo + +# --- 3: the answer does not depend on the decomposition -------------------- +ref=${RANKS%% *} +for n in $RANKS; do + [ "$n" = "$ref" ] && continue + paste "cons.$ref" "cons.$n" | awk -F'[ |\t]+' -v n="$n" -v r="$ref" -v t="$TOL_RANK_AGREE" ' + { + m1 = $5; k1 = $13; m2 = $18; k2 = $26 + d = (m2 - m1) / m1; if (d < 0) d = -d; if (d > dm) dm = d + d = (k2 - k1) / k1; if (d < 0) d = -d; if (d > dk) dk = d + } + END { + ok = 1 + if (dm > t) { printf "FAIL mass differs between %s and %s ranks: %.3e > %s\n", r, n, dm, t; ok = 0 } + else { printf "ok mass agrees between %s and %s ranks (%.3e)\n", r, n, dm } + if (dk > t) { printf "FAIL kinetic energy differs between %s and %s ranks: %.3e > %s\n", r, n, dk, t; ok = 0 } + else { printf "ok kinetic energy agrees between %s and %s ranks (%.3e)\n", r, n, dk } + exit ok ? 0 : 1 + }' || status=1 +done + +echo +if [ "$status" -ne 0 ]; then + echo "regression test FAILED" + exit 1 +fi +echo "regression test passed"