feat: FP stability CI (Verrou) + reduced-energy scheme + HLLC xi-factor fix

.github/workflows/fp-stability.yml

@@ -0,0 +1,92 @@
+name: FP Stability
+
+# Runs the Verrou-based floating-point stability suite.
+#
+# What it tests (./mfc.sh fp-stability):
+#   sod_standard     1-D standard Sod, p_L/p_R=10, ideal gas (well-conditioned baseline)
+#                    25 cells, 5 steps, WENO5 + HLLC
+#                    Threshold 1e-13 — should always pass
+#
+#   sod_strong       1-D Sod shock, p_L/p_R=100,000, ideal gas
+#                    50 cells, 10 steps, WENO5 + HLLC
+#                    Threshold 1e-10
+#                    Probes: HLLC xi-factor cancellation near sonic contact
+#                            (s_L - vel_L)/(s_L - s_S) when s_L ≈ s_S
+#
+#   water_stiffened  1-D water shock, stiffened EOS (pi_inf=4046)
+#                    50 cells, 10 steps, WENO5 + HLLC
+#                    Threshold 1e-8 (loosened; tightens to 1e-10 once Etilde scheme merges)
+#                    Probes: pressure-recovery cancellation p=(E-pi_inf)/gamma
+#                            loses ~4 decimal digits when pi_inf/p_right ~ 40,000
+#
+# For each case: 1 nearest-rounding reference run + N random-rounding runs.
+# PASS if max L∞ deviation across all N samples stays below threshold.
+# On FAIL: verrou_dd_sym runs to identify the responsible function symbols.
+# Logs are uploaded as CI artifacts.
+#
+# Verrou (Valgrind 3.26.0 + edf-hpc/verrou@a58d434) is built once and cached.
+# Build takes ~20 min uncached; cached runs restore in ~30 s.
+
+on:
+  push:
+  workflow_dispatch:
+
+jobs:
+  fp-stability:
+    name: Floating-Point Stability (Verrou)
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Clone
+        uses: actions/checkout@v4
+
+      - name: Cache Verrou
+        id: cache-verrou
+        uses: actions/cache@v4
+        with:
+          path: ~/.local/verrou
+          key: verrou-a58d434-valgrind-3.26.0-${{ runner.os }}
+
+      - name: Install system dependencies
+        run: |
+          sudo apt-get update -y
+          sudo apt-get install -y \
+            build-essential automake python3 python3-numpy libc6-dbg \
+            cmake gfortran
+
+      - name: Build Verrou
+        if: steps.cache-verrou.outputs.cache-hit != 'true'
+        run: |
+          cd /tmp
+          wget -q https://sourceware.org/pub/valgrind/valgrind-3.26.0.tar.bz2
+          tar xf valgrind-3.26.0.tar.bz2
+
+          git clone https://github.com/edf-hpc/verrou.git
+          git -C verrou checkout a58d434
+
+          # Merge Verrou into Valgrind source tree and patch
+          cp -r verrou valgrind-3.26.0/verrou
+          cd valgrind-3.26.0
+          cat verrou/valgrind.*diff | patch -p1
+
+          ./autogen.sh
+          ./configure --enable-only64bit --prefix="$HOME/.local/verrou"
+          make -j"$(nproc)"
+          make install
+
+      - name: Verify Verrou
+        run: ~/.local/verrou/bin/valgrind --version
+
+      - name: Build MFC (debug, serial)
+        run: ./mfc.sh build -t pre_process simulation --no-mpi --debug -j"$(nproc)"
+
+      - name: Run FP Stability Suite
+        run: ./mfc.sh fp-stability -N 5
+
+      - name: Upload FP stability logs
+        if: always()
+        uses: actions/upload-artifact@v4
+        with:
+          name: fp-stability-logs
+          path: fp-stability-logs/
+          if-no-files-found: ignore

.gitignore

@@ -37,6 +37,7 @@ docs/documentation/parameters.md
 /tests/*/**
 !/tests/*/golden.txt
 !/tests/*/golden-metadata.txt
+!/tests/fp_stability/**
 
 # NVIDIA Nsight Compute
 *.nsys-rep

src/simulation/m_riemann_solvers.fpp

@@ -2040,8 +2040,8 @@ contains
                                 s_M = min(0._wp, s_L); s_P = max(0._wp, s_R)
 
                                 ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
-                                xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S)
-                                xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S)
+                                xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
+                                xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
 
                                 ! goes with numerical star velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
                                 xi_M = (5.e-1_wp + sign(0.5_wp, s_S))
@@ -2331,8 +2331,8 @@ contains
                                 s_M = min(0._wp, s_L); s_P = max(0._wp, s_R)
 
                                 ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
-                                xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S)
-                                xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S)
+                                xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
+                                xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
 
                                 ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
                                 xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))
@@ -2688,8 +2688,8 @@ contains
                                 s_M = min(0._wp, s_L); s_P = max(0._wp, s_R)
 
                                 ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
-                                xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S)
-                                xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S)
+                                xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
+                                xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
 
                                 ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
                                 xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))
@@ -3127,8 +3127,8 @@ contains
                                 s_M = min(0._wp, s_L); s_P = max(0._wp, s_R)
 
                                 ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
-                                xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S)
-                                xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S)
+                                xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
+                                xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
 
                                 ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
                                 xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))