Add direction symmetry, MPI consistency, restart roundtrip, and kernel tests

sbryngelson · claude · sbryngelson · commit 8d19a799ca06 · 2026-02-20T23:29:37.000-05:00
Four new test categories to catch edge-case bugs:

1. Direction symmetry: 3D shock tube tests with shock propagating in x
   and y directions (default 3D tests use z), catching direction-specific
   bugs in reconstruction and gradient calculations.

2. MPI consistency: ppn=2 tests for bubbles, viscous flows, and
   hypoelasticity, catching broadcast/reduction bugs in MPI-sensitive
   physics modules.

3. Restart roundtrip: 1D and 3D tests that run a straight simulation,
   then restart from the midpoint and compare the restarted output
   against the straight run output to verify restart I/O fidelity.

4. Kernel golden values: polydisperse bubble test (nb=3 with poly_sigma
   distribution) exercising bubble array indexing code paths.

Co-Authored-By: Claude Sonnet 4.6 &lt;noreply@anthropic.com&gt;
diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py
@@ -108,11 +108,13 @@ class TestCase(case.Case):
     ppn:          int
     trace:        str
     override_tol: Optional[float] = None
+    restart_check: bool = False
 
-    def __init__(self, trace: str, mods: dict, ppn: int = None, override_tol: float = None) -> None:
+    def __init__(self, trace: str, mods: dict, ppn: int = None, override_tol: float = None, restart_check: bool = False) -> None:
         self.trace        = trace
         self.ppn          = ppn or 1
         self.override_tol = override_tol
+        self.restart_check = restart_check
         super().__init__({**BASE_CFG.copy(), **mods})
 
     def run(self, targets: List[Union[str, MFCTarget]], gpus: Set[int]) -> subprocess.CompletedProcess:
@@ -140,6 +142,36 @@ def run(self, targets: List[Union[str, MFCTarget]], gpus: Set[int]) -> subproces
 
         return common.system(command, print_cmd=False, text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
 
+    def run_restart(self, targets, gpus):
+        """Run a restart roundtrip: simulate to midpoint, then restart to end."""
+        mid_step = self.params['t_step_stop'] // 2
+        orig = dict(self.params)
+
+        try:
+            self.delete_output()
+
+            # Phase 1: Run to midpoint (generates restart data)
+            self.params = {**orig, 't_step_stop': mid_step, 't_step_save': mid_step}
+            self.create_directory()
+            result1 = self.run(targets, gpus)
+            if result1.returncode != 0:
+                return result1
+
+            # Delete output data but keep restart_data/
+            dirpath = self.get_dirpath()
+            common.delete_directory(os.path.join(dirpath, "D"))
+            common.delete_directory(os.path.join(dirpath, "p_all"))
+            common.delete_directory(os.path.join(dirpath, "silo_hdf5"))
+
+            # Phase 2: Restart from midpoint
+            self.params = {**orig, 'old_ic': 'T', 'old_grid': 'T',
+                           't_step_start': mid_step,
+                           't_step_save': orig['t_step_stop'] - mid_step}
+            self.create_directory()
+            return self.run(targets, gpus)
+        finally:
+            self.params = orig
+
     def get_trace(self) -> str:
         return self.trace
 
@@ -278,6 +310,7 @@ class TestCaseBuilder:
     ppn:          int
     functor:      Optional[Callable]
     override_tol: Optional[float] = None
+    restart_check: bool = False
 
     def get_uuid(self) -> str:
         return trace_to_uuid(self.trace)
@@ -302,7 +335,7 @@ def to_case(self) -> TestCase:
         if self.functor:
             self.functor(dictionary)
 
-        return TestCase(self.trace, dictionary, self.ppn, self.override_tol)
+        return TestCase(self.trace, dictionary, self.ppn, self.override_tol, self.restart_check)
 
 
 @dataclasses.dataclass
@@ -330,7 +363,7 @@ def define_case_f(trace: str, path: str, args: List[str] = None, ppn: int = None
 
 
 # pylint: disable=too-many-arguments, too-many-positional-arguments
-def define_case_d(stack: CaseGeneratorStack, newTrace: str, newMods: dict, ppn: int = None, functor: Callable = None, override_tol: float = None) -> TestCaseBuilder:
+def define_case_d(stack: CaseGeneratorStack, newTrace: str, newMods: dict, ppn: int = None, functor: Callable = None, override_tol: float = None, restart_check: bool = False) -> TestCaseBuilder:
     mods: dict = {}
 
     for mod in stack.mods:
@@ -346,7 +379,7 @@ def define_case_d(stack: CaseGeneratorStack, newTrace: str, newMods: dict, ppn:
         if not common.isspace(trace):
             traces.append(trace)
 
-    return TestCaseBuilder(' -> '.join(traces), mods, None, None, ppn or 1, functor, override_tol)
+    return TestCaseBuilder(' -> '.join(traces), mods, None, None, ppn or 1, functor, override_tol, restart_check)
 
 def input_bubbles_lagrange(self):
     if "lagrange_bubblescreen" in self.trace:
diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py
@@ -1109,6 +1109,227 @@ def chemistry_cases():
 
     chemistry_cases()
 
+    def direction_symmetry_tests():
+        """3D tests with shock propagating in x and y directions.
+
+        Default 3D tests have the shock along z. These test x and y
+        code paths to catch direction-specific bugs in reconstruction,
+        Riemann solvers, and gradient calculations.
+        """
+        for direction in ['x', 'y']:
+            others = [d for d in ['x', 'y', 'z'] if d != direction]
+            mods = {
+                'm': 24, 'n': 24, 'p': 24,
+                'x_domain%beg': 0.E+00, 'x_domain%end': 1.E+00,
+                'y_domain%beg': 0.E+00, 'y_domain%end': 1.E+00,
+                'z_domain%beg': 0.E+00, 'z_domain%end': 1.E+00,
+                'bc_x%beg': -3, 'bc_x%end': -3,
+                'bc_y%beg': -3, 'bc_y%end': -3,
+                'bc_z%beg': -3, 'bc_z%end': -3,
+            }
+
+            centroids = [0.05, 0.45, 0.9]
+            lengths   = [0.1,  0.7,  0.2]
+
+            for patchID in range(1, 4):
+                mods[f'patch_icpp({patchID})%geometry'] = 9
+                mods[f'patch_icpp({patchID})%vel(1)'] = 0.0
+                mods[f'patch_icpp({patchID})%vel(2)'] = 0.0
+                mods[f'patch_icpp({patchID})%vel(3)'] = 0.0
+                mods[f'patch_icpp({patchID})%{direction}_centroid'] = centroids[patchID - 1]
+                mods[f'patch_icpp({patchID})%length_{direction}']   = lengths[patchID - 1]
+                for od in others:
+                    mods[f'patch_icpp({patchID})%{od}_centroid'] = 0.5
+                    mods[f'patch_icpp({patchID})%length_{od}']   = 1
+
+            stack.push(f'3D Direction Symmetry -> Shock in {direction.upper()}', mods)
+            cases.append(define_case_d(stack, '', {}))
+            stack.pop()
+
+    direction_symmetry_tests()
+
+    def mpi_consistency_tests():
+        """ppn=2 tests for physics sensitive to MPI decomposition.
+
+        Exercises bubble dynamics, viscous flows, and hypoelasticity
+        with 2 MPI ranks to catch broadcast/reduction bugs.
+        """
+        base_3d = {
+            'm': 29, 'n': 29, 'p': 49,
+            'x_domain%beg': 0.E+00, 'x_domain%end': 1.E+00,
+            'y_domain%beg': 0.E+00, 'y_domain%end': 1.E+00,
+            'z_domain%beg': 0.E+00, 'z_domain%end': 1.E+00,
+            'bc_x%beg': -3, 'bc_x%end': -3,
+            'bc_y%beg': -3, 'bc_y%end': -3,
+            'bc_z%beg': -3, 'bc_z%end': -3,
+        }
+
+        for patchID in range(1, 4):
+            base_3d[f'patch_icpp({patchID})%geometry'] = 9
+            base_3d[f'patch_icpp({patchID})%vel(1)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%vel(2)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%vel(3)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%x_centroid'] = 0.5
+            base_3d[f'patch_icpp({patchID})%length_x'] = 1
+            base_3d[f'patch_icpp({patchID})%y_centroid'] = 0.5
+            base_3d[f'patch_icpp({patchID})%length_y'] = 1
+        base_3d.update({
+            'patch_icpp(1)%z_centroid': 0.05, 'patch_icpp(1)%length_z': 0.1,
+            'patch_icpp(2)%z_centroid': 0.45, 'patch_icpp(2)%length_z': 0.7,
+            'patch_icpp(3)%z_centroid': 0.9,  'patch_icpp(3)%length_z': 0.2,
+        })
+
+        # Bubbles with 2 MPI ranks
+        stack.push('MPI Consistency -> 3D -> Bubbles', {**base_3d,
+            'dt': 1e-06,
+            'bubbles_euler': 'T', 'nb': 1, 'polytropic': 'T', 'bubble_model': 2,
+            'fluid_pp(1)%gamma': 0.16, 'fluid_pp(1)%pi_inf': 3515.0,
+            'bub_pp%R0ref': 1.0, 'bub_pp%p0ref': 1.0, 'bub_pp%rho0ref': 1.0, 'bub_pp%T0ref': 1.0,
+            'bub_pp%ss': 0.07179866765358993, 'bub_pp%pv': 0.02308216136195411,
+            'bub_pp%vd': 0.2404125083932959,
+            'bub_pp%mu_l': 0.009954269975623244, 'bub_pp%mu_v': 8.758168074360729e-05,
+            'bub_pp%mu_g': 0.00017881922111898042, 'bub_pp%gam_v': 1.33, 'bub_pp%gam_g': 1.4,
+            'bub_pp%M_v': 18.02, 'bub_pp%M_g': 28.97, 'bub_pp%k_v': 0.5583395141263873,
+            'bub_pp%k_g': 0.7346421281308791, 'bub_pp%R_v': 1334.8378710170155,
+            'bub_pp%R_g': 830.2995663005393,
+            'patch_icpp(1)%alpha_rho(1)': 0.96, 'patch_icpp(1)%alpha(1)': 4e-02,
+            'patch_icpp(2)%alpha_rho(1)': 0.96, 'patch_icpp(2)%alpha(1)': 4e-02,
+            'patch_icpp(3)%alpha_rho(1)': 0.96, 'patch_icpp(3)%alpha(1)': 4e-02,
+            'patch_icpp(1)%pres': 1.0, 'patch_icpp(2)%pres': 1.0, 'patch_icpp(3)%pres': 1.0,
+        })
+        cases.append(define_case_d(stack, '', {}, ppn=2))
+        stack.pop()
+
+        # Viscous with 2 MPI ranks
+        stack.push('MPI Consistency -> 3D -> Viscous', {**base_3d,
+            'dt': 1e-11,
+            'fluid_pp(1)%Re(1)': 0.0001, 'viscous': 'T',
+            'patch_icpp(1)%vel(1)': 1.0,
+            'patch_icpp(2)%vel(1)': 1.0,
+            'patch_icpp(3)%vel(1)': 1.0,
+        })
+        cases.append(define_case_d(stack, '', {}, ppn=2))
+        stack.pop()
+
+        # Hypoelasticity with 2 MPI ranks
+        stack.push('MPI Consistency -> 3D -> Hypoelasticity', {**base_3d,
+            'hypoelasticity': 'T', 'riemann_solver': 1, 'fd_order': 4,
+            'fluid_pp(1)%gamma': 0.3, 'fluid_pp(1)%pi_inf': 7.8E+05,
+            'fluid_pp(1)%G': 1.E+05,
+            'patch_icpp(1)%pres': 1.E+06, 'patch_icpp(1)%alpha_rho(1)': 1000.E+00,
+            'patch_icpp(2)%pres': 1.E+05, 'patch_icpp(2)%alpha_rho(1)': 1000.E+00,
+            'patch_icpp(3)%pres': 5.E+05, 'patch_icpp(3)%alpha_rho(1)': 1000.E+00,
+            'patch_icpp(1)%tau_e(1)': 0.E+00, 'patch_icpp(2)%tau_e(1)': 0.E+00,
+            'patch_icpp(3)%tau_e(1)': 0.E+00,
+            'patch_icpp(1)%tau_e(2)': 0.E+00, 'patch_icpp(1)%tau_e(3)': 0.E+00,
+            'patch_icpp(2)%tau_e(2)': 0.E+00, 'patch_icpp(2)%tau_e(3)': 0.E+00,
+            'patch_icpp(3)%tau_e(2)': 0.E+00, 'patch_icpp(3)%tau_e(3)': 0.E+00,
+            'patch_icpp(1)%tau_e(4)': 0.E+00, 'patch_icpp(1)%tau_e(5)': 0.E+00,
+            'patch_icpp(1)%tau_e(6)': 0.E+00,
+            'patch_icpp(2)%tau_e(4)': 0.E+00, 'patch_icpp(2)%tau_e(5)': 0.E+00,
+            'patch_icpp(2)%tau_e(6)': 0.E+00,
+            'patch_icpp(3)%tau_e(4)': 0.E+00, 'patch_icpp(3)%tau_e(5)': 0.E+00,
+            'patch_icpp(3)%tau_e(6)': 0.E+00,
+        })
+        cases.append(define_case_d(stack, '', {}, ppn=2))
+        stack.pop()
+
+    mpi_consistency_tests()
+
+    def restart_roundtrip_tests():
+        """Tests that verify save-restart roundtrip fidelity.
+
+        Each test runs a straight simulation, then a restart from the
+        midpoint. The restarted output is compared against the straight
+        run output to verify restart I/O doesn't introduce drift.
+        """
+        # 1D restart
+        stack.push('Restart Roundtrip -> 1D', {
+            'm': 299, 'n': 0, 'p': 0,
+            'x_domain%beg': 0.E+00, 'x_domain%end': 1.E+00,
+            'bc_x%beg': -3, 'bc_x%end': -3,
+            'patch_icpp(1)%geometry': 1, 'patch_icpp(2)%geometry': 1,
+            'patch_icpp(3)%geometry': 1,
+            'patch_icpp(1)%x_centroid': 0.05, 'patch_icpp(1)%length_x': 0.1,
+            'patch_icpp(2)%x_centroid': 0.45, 'patch_icpp(2)%length_x': 0.7,
+            'patch_icpp(3)%x_centroid': 0.9,  'patch_icpp(3)%length_x': 0.2,
+            'patch_icpp(1)%vel(1)': 0.0, 'patch_icpp(2)%vel(1)': 0.0,
+            'patch_icpp(3)%vel(1)': 0.0,
+        })
+        cases.append(define_case_d(stack, '', {}, restart_check=True))
+        stack.pop()
+
+        # 3D restart
+        base_3d = {
+            'm': 24, 'n': 24, 'p': 24,
+            'x_domain%beg': 0.E+00, 'x_domain%end': 1.E+00,
+            'y_domain%beg': 0.E+00, 'y_domain%end': 1.E+00,
+            'z_domain%beg': 0.E+00, 'z_domain%end': 1.E+00,
+            'bc_x%beg': -3, 'bc_x%end': -3,
+            'bc_y%beg': -3, 'bc_y%end': -3,
+            'bc_z%beg': -3, 'bc_z%end': -3,
+        }
+        for patchID in range(1, 4):
+            base_3d[f'patch_icpp({patchID})%geometry'] = 9
+            base_3d[f'patch_icpp({patchID})%vel(1)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%vel(2)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%vel(3)'] = 0.0
+            base_3d[f'patch_icpp({patchID})%x_centroid'] = 0.5
+            base_3d[f'patch_icpp({patchID})%length_x'] = 1
+            base_3d[f'patch_icpp({patchID})%y_centroid'] = 0.5
+            base_3d[f'patch_icpp({patchID})%length_y'] = 1
+        base_3d.update({
+            'patch_icpp(1)%z_centroid': 0.05, 'patch_icpp(1)%length_z': 0.1,
+            'patch_icpp(2)%z_centroid': 0.45, 'patch_icpp(2)%length_z': 0.7,
+            'patch_icpp(3)%z_centroid': 0.9,  'patch_icpp(3)%length_z': 0.2,
+        })
+        stack.push('Restart Roundtrip -> 3D', base_3d)
+        cases.append(define_case_d(stack, '', {}, restart_check=True))
+        stack.pop()
+
+    restart_roundtrip_tests()
+
+    def kernel_golden_tests():
+        """Focused golden-value tests for specific physics kernels.
+
+        Polydisperse bubbles: exercises bubble array indexing with
+        multiple bubble sizes and weight distributions.
+        """
+        # Polydisperse bubbles (nb=3 with poly_sigma distribution)
+        stack.push('Kernel -> 1D -> Polydisperse Bubbles', {
+            'm': 299, 'n': 0, 'p': 0,
+            'x_domain%beg': 0.E+00, 'x_domain%end': 1.E+00,
+            'bc_x%beg': -3, 'bc_x%end': -3,
+            'patch_icpp(1)%geometry': 1, 'patch_icpp(2)%geometry': 1,
+            'patch_icpp(3)%geometry': 1,
+            'patch_icpp(1)%x_centroid': 0.05, 'patch_icpp(1)%length_x': 0.1,
+            'patch_icpp(2)%x_centroid': 0.45, 'patch_icpp(2)%length_x': 0.7,
+            'patch_icpp(3)%x_centroid': 0.9,  'patch_icpp(3)%length_x': 0.2,
+            'patch_icpp(1)%vel(1)': 0.0, 'patch_icpp(2)%vel(1)': 0.0,
+            'patch_icpp(3)%vel(1)': 0.0,
+            'dt': 1e-07,
+            'bubbles_euler': 'T', 'nb': 3, 'polydisperse': 'T', 'poly_sigma': 0.3,
+            'polytropic': 'T', 'bubble_model': 2,
+            'fluid_pp(1)%gamma': 0.16, 'fluid_pp(1)%pi_inf': 3515.0,
+            'bub_pp%R0ref': 1.0, 'bub_pp%p0ref': 1.0,
+            'bub_pp%rho0ref': 1.0, 'bub_pp%T0ref': 1.0,
+            'bub_pp%ss': 0.07179866765358993, 'bub_pp%pv': 0.02308216136195411,
+            'bub_pp%vd': 0.2404125083932959,
+            'bub_pp%mu_l': 0.009954269975623244, 'bub_pp%mu_v': 8.758168074360729e-05,
+            'bub_pp%mu_g': 0.00017881922111898042, 'bub_pp%gam_v': 1.33, 'bub_pp%gam_g': 1.4,
+            'bub_pp%M_v': 18.02, 'bub_pp%M_g': 28.97, 'bub_pp%k_v': 0.5583395141263873,
+            'bub_pp%k_g': 0.7346421281308791, 'bub_pp%R_v': 1334.8378710170155,
+            'bub_pp%R_g': 830.2995663005393,
+            'patch_icpp(1)%alpha_rho(1)': 0.96, 'patch_icpp(1)%alpha(1)': 4e-02,
+            'patch_icpp(2)%alpha_rho(1)': 0.96, 'patch_icpp(2)%alpha(1)': 4e-02,
+            'patch_icpp(3)%alpha_rho(1)': 0.96, 'patch_icpp(3)%alpha(1)': 4e-02,
+            'patch_icpp(1)%pres': 1.0, 'patch_icpp(2)%pres': 1.0, 'patch_icpp(3)%pres': 1.0,
+        })
+        cases.append(define_case_d(stack, '', {}))
+        stack.pop()
+
+    kernel_golden_tests()
+
     # Sanity Check 1
     if stack.size() != 0:
         raise common.MFCException("list_cases: stack isn't fully pop'ed")
diff --git a/toolchain/mfc/test/test.py b/toolchain/mfc/test/test.py
@@ -384,6 +384,34 @@ def _handle_case(case: TestCase, devices: typing.Set[int]):
                 if msg is not None:
                     raise MFCException(f"Test {case}: {msg}")
 
+        # Restart roundtrip verification: run to midpoint, restart,
+        # and compare restarted output against the straight run.
+        if getattr(case, 'restart_check', False) and not ARG("add_new_variables"):
+            straight_pack = pack
+
+            if timeout_flag.is_set():
+                raise TestTimeoutError("Test case exceeded 1 hour timeout")
+
+            restart_result = case.run_restart([PRE_PROCESS, SIMULATION], devices)
+            out_filepath_restart = os.path.join(case.get_dirpath(), "out_restart.txt")
+            common.file_write(out_filepath_restart, restart_result.stdout)
+
+            if restart_result.returncode != 0:
+                cons.print(restart_result.stdout)
+                raise MFCException(f"Test {case}: Restart roundtrip run failed.")
+
+            restart_pack, restart_err = packer.pack(case.get_dirpath())
+            if restart_err is not None:
+                raise MFCException(f"Test {case}: Restart pack error: {restart_err}")
+
+            if restart_pack.has_NaNs():
+                raise MFCException(f"Test {case}: NaNs detected in restarted output.")
+
+            _, restart_msg = packtol.compare(
+                restart_pack, straight_pack, packtol.Tolerance(tol, tol))
+            if restart_msg is not None:
+                raise MFCException(f"Test {case}: Restart roundtrip mismatch: {restart_msg}")
+
         if ARG("test_all"):
             case.delete_output()
             # Check timeout before launching the (potentially long) post-process run
