Add a physical regression test and the diagnostics it needs

.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

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<double>(lmass), Kokkos::Sum<double>(lmx),
+        Kokkos::Sum<double>(lmy), Kokkos::Sum<double>(lmz),
+        Kokkos::Sum<double>(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

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);
 

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);
             }
         }

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 <it>  4 mass  5 <mass>  6 px  7 <px>
+#   8 py  9 <py>  10 pz  11 <pz>  12 ke  13 <ke>
+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"