From 861f72ea269cc8ee272b41a133349e261fc0fcaf Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 9 Oct 2026 23:06:07 -0400 Subject: [PATCH 1/4] Reject characteristic BCs with MUSCL reconstruction s_cbc computes the characteristic boundary fluxes only when recon_type == WENO. With MUSCL the routine copies the fluxes in and out unchanged, so a characteristic BC (-5 to -12) silently behaves like extrapolation (-3): a MUSCL run with bc = -5 or -6 is bit-identical to the same run with bc = -3. Prohibit the combination in check_muscl_simulation, sharing the characteristic-BC scan with the existing IGR check, and add unit tests. 2D_shockdroplet_muscl and 3D_shockdroplet_muscl used bc_x = -6 with MUSCL; switch them to -3, which they were already getting. Their outputs are unchanged (bit-identical p_all and D on the 2D test case; tests BE796E10 and 23573861 pass against the existing goldens). (cherry picked from commit 8367254f08bd029e8534ee3e4d23c6776659c3a5) --- examples/2D_shockdroplet_muscl/case.py | 4 ++-- examples/3D_shockdroplet_muscl/case.py | 4 ++-- toolchain/mfc/case_validator.py | 18 +++++++++++------- toolchain/mfc/test_case_validator.py | 15 +++++++++++++++ 4 files changed, 30 insertions(+), 11 deletions(-) diff --git a/examples/2D_shockdroplet_muscl/case.py b/examples/2D_shockdroplet_muscl/case.py index 61b904ec85..4367cd079a 100755 --- a/examples/2D_shockdroplet_muscl/case.py +++ b/examples/2D_shockdroplet_muscl/case.py @@ -64,8 +64,8 @@ "riemann_solver": "hllc", "wave_speeds": "direct", "avg_state": "arithmetic", - "bc_x%beg": -6, # 11, - "bc_x%end": -6, # 12 + "bc_x%beg": -3, + "bc_x%end": -3, "bc_y%beg": -2, "bc_y%end": -3, # Formatted Database Files Structure Parameters diff --git a/examples/3D_shockdroplet_muscl/case.py b/examples/3D_shockdroplet_muscl/case.py index 8213b85105..67bd4954e5 100644 --- a/examples/3D_shockdroplet_muscl/case.py +++ b/examples/3D_shockdroplet_muscl/case.py @@ -206,8 +206,8 @@ "riemann_solver": "hllc", "wave_speeds": "direct", "avg_state": "arithmetic", - "bc_x%beg": -6, - "bc_x%end": -6, + "bc_x%beg": -3, + "bc_x%end": -3, "bc_y%beg": -2, "bc_y%end": -3, "bc_z%beg": -2, diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 36fce10a00..befe499f30 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -374,6 +374,11 @@ def _check_order_fits_grid(self, order: int, param_name: str) -> None: def _get_recon_type(self) -> int: return self.get("recon_type", 1) + def _characteristic_bcs(self) -> list: + """Names of the boundaries set to a characteristic BC, BC_CHAR_SLIP_WALL (-5) to BC_CHAR_SUP_OUTFLOW (-12)""" + names = [f"bc_{d}%{b}" for d in ["x", "y", "z"] for b in ["beg", "end"]] + return [name for name in names if self.get(name) is not None and -12 <= self.get(name) <= -5] + def check_parameter_types(self): """Validate parameter types before other checks. @@ -1390,6 +1395,10 @@ def check_muscl_simulation(self): muscl_lim = self.get("muscl_lim") muscl_eps = self.get("muscl_eps") + # s_cbc computes characteristic fluxes for WENO only, so MUSCL would silently ignore a CBC + for name in self._characteristic_bcs(): + self.prohibit(True, f"Characteristic boundary condition {name} is not compatible with recon_type = 2 (MUSCL)") + if muscl_order is None: return @@ -1681,13 +1690,8 @@ def check_igr_simulation(self): self.prohibit(probe_wrt, "IGR does not support probe writes") self.prohibit(int_comp > 0, "IGR does not support int_comp > 0") - # Check BCs - IGR does not support characteristic BCs - # Characteristic BCs are BC_CHAR_SLIP_WALL (-5) through BC_CHAR_SUP_OUTFLOW (-12) - for dir in ["x", "y", "z"]: - for bound in ["beg", "end"]: - bc = self.get(f"bc_{dir}%{bound}") - if bc is not None: - self.prohibit(-12 <= bc <= -5, f"Characteristic boundary condition bc_{dir}%{bound} is not compatible with IGR") + for name in self._characteristic_bcs(): + self.prohibit(True, f"Characteristic boundary condition {name} is not compatible with IGR") def check_acoustic_source(self): """Checks acoustic source parameters (simulation)""" diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index c0170452d9..1000ffbbfa 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -704,3 +704,18 @@ def test_rejects_isothermal_wall_without_any_heat_path(self): {**BASE, "bc_x%beg": -16, "bc_x%end": -16, "bc_x%isothermal_in": "T", "bc_x%Twall_in": 300.0}, "requires a heat-conduction path", ) + + +class TestMusclCharacteristicBcs(ConstraintTestCase): + """s_cbc computes characteristic fluxes for WENO only, so MUSCL silently ignored a CBC.""" + + MUSCL = {**{k: v for k, v in BASE.items() if not k.startswith("weno")}, "recon_type": 2, "muscl_order": 2, "muscl_lim": 1} + + def test_rejects_cbc_with_muscl(self): + self.assertRejects({**self.MUSCL, "bc_x%end": -6}, "bc_x%end is not compatible with recon_type = 2 (MUSCL)") + + def test_accepts_non_characteristic_bc_with_muscl(self): + self.assertAccepts({**self.MUSCL, "bc_x%beg": -3, "bc_x%end": -3}) + + def test_accepts_cbc_with_weno(self): + self.assertAccepts({**BASE, "bc_x%beg": -6, "bc_x%end": -6}) From 15e78cdcdf7555a75ca545d55de2e9a5e13ec35d Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 10 Oct 2026 00:38:50 -0400 Subject: [PATCH 2/4] Restore input checks for patch geometries 16-19 #942 (febf1315) replaced the hard-coded geometries 7, 13 and 15 with hcid and, in the same edit, dropped the print-only branches for geometries 16-19 from s_check_patches. Those geometries (1D bubble pulse, spiral, 2D and 3D varcircle) are still dispatched in s_apply_icpp_patches and documented in case.md, but since then they fall through to "geometry must be between 1 and 21" and cannot be used. Give each a check routine that enforces its dimensionality (the dispatch silently skips a mismatched one) and the parameters case.md lists. case.md called the varcircle thickness `thickness`; the code reads `epsilon`, so name that instead. Also fix the swapped circle messages: n == 0 now reports "n must be greater than zero" and p > 0 "p must be zero". Verified: a 1D geometry-16 and a 2D geometry-18 case now pass pre_process (master rejects them), and geometry 16 in 2D is rejected with "1D bubble pulse patch 2: n must be zero". (cherry picked from commit 877d0b00650b064b70e94c0ff9b3c2ba2466d4f5) --- docs/documentation/case.md | 4 +- src/pre_process/m_check_patches.fpp | 77 ++++++++++++++++++++++++++++- 2 files changed, 77 insertions(+), 4 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 4efcc06e03..7a8f2b4e74 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -1428,8 +1428,8 @@ This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` | 15 | N/A | N/A | N/A | No longer exists. Empty. | | 16 | 1D bubble pulse | 1 | N | Requires `x_centroid`, `length_x` | | 17 | Spiral | 2 | N | Requires `[x,y]_centroid` | -| 18 | 2D Varcircle | 2 | Y | Requires `[x,y]_centroid`, `radius`, and `thickness` | -| 19 | 3D Varcircle | 3 | Y | Requires `[x,y,z]_centroid`, `length_z`, `radius`, and `thickness` | +| 18 | 2D Varcircle | 2 | Y | Requires `[x,y]_centroid`, `radius`, and `epsilon` (thickness) | +| 19 | 3D Varcircle | 3 | Y | Requires `[x,y,z]_centroid`, `length_z`, `radius`, and `epsilon` (thickness) | | 20 | 2D Taylor-Green Vortex | 2 | N | Requires `[x,y]_centroid`, `length_x`, `length_y`, `vel(1)`, and `vel(2)` | | 21 | Model | 2 & 3 | Y | Imports a Model (STL/OBJ). Requires `model_id`. | diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp index 5d31c8dc09..6bba99bd65 100644 --- a/src/pre_process/m_check_patches.fpp +++ b/src/pre_process/m_check_patches.fpp @@ -77,6 +77,14 @@ contains call s_check_2d_modal_patch_geometry(i) else if (patch_icpp(i)%geometry == 14) then call s_check_3d_spherical_harmonic_patch_geometry(i) + else if (patch_icpp(i)%geometry == 16) then + call s_check_1d_bubble_pulse_patch_geometry(i) + else if (patch_icpp(i)%geometry == 17) then + call s_check_spiral_patch_geometry(i) + else if (patch_icpp(i)%geometry == 18) then + call s_check_varcircle_patch_geometry(i) + else if (patch_icpp(i)%geometry == 19) then + call s_check_3d_varcircle_patch_geometry(i) else if (patch_icpp(i)%geometry == 20) then call s_check_2D_TaylorGreen_vortex_patch_geometry(i) else if (patch_icpp(i)%geometry == 21) then @@ -148,8 +156,8 @@ contains call s_int_to_str(patch_id, iStr) - @:PROHIBIT(n == 0, "Circle patch "//trim(iStr)//": n must be zero") - @:PROHIBIT(p > 0, "Circle patch "//trim(iStr)//": p must be greater than zero") + @:PROHIBIT(n == 0, "Circle patch "//trim(iStr)//": n must be greater than zero") + @:PROHIBIT(p > 0, "Circle patch "//trim(iStr)//": p must be zero") @:PROHIBIT(patch_icpp(patch_id)%radius <= 0._wp, "Circle patch "//trim(iStr)//": radius must be greater than zero") @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "Circle patch "//trim(iStr)//": x_centroid must be set") @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "Circle patch "//trim(iStr)//": y_centroid must be set") @@ -229,6 +237,71 @@ contains end subroutine s_check_2D_TaylorGreen_vortex_patch_geometry + !> Check the 1D bubble pulse patch input + impure subroutine s_check_1d_bubble_pulse_patch_geometry(patch_id) + + integer, intent(in) :: patch_id + + call s_int_to_str(patch_id, iStr) + + @:PROHIBIT(n > 0, "1D bubble pulse patch "//trim(iStr)//": n must be zero") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), & + & "1D bubble pulse patch " // trim(iStr) // ": x_centroid must be set") + @:PROHIBIT(patch_icpp(patch_id)%length_x <= 0._wp, & + & "1D bubble pulse patch " // trim(iStr) // ": length_x must be greater than zero") + + end subroutine s_check_1d_bubble_pulse_patch_geometry + + !> Check the spiral patch input + impure subroutine s_check_spiral_patch_geometry(patch_id) + + integer, intent(in) :: patch_id + + call s_int_to_str(patch_id, iStr) + + @:PROHIBIT(n == 0, "Spiral patch "//trim(iStr)//": n must be greater than zero") + @:PROHIBIT(p > 0, "Spiral patch "//trim(iStr)//": p must be zero") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "Spiral patch "//trim(iStr)//": x_centroid must be set") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "Spiral patch "//trim(iStr)//": y_centroid must be set") + + end subroutine s_check_spiral_patch_geometry + + !> Check the 2D varcircle patch input; epsilon is the ring thickness + impure subroutine s_check_varcircle_patch_geometry(patch_id) + + integer, intent(in) :: patch_id + + call s_int_to_str(patch_id, iStr) + + @:PROHIBIT(n == 0, "Varcircle patch "//trim(iStr)//": n must be greater than zero") + @:PROHIBIT(p > 0, "Varcircle patch "//trim(iStr)//": p must be zero") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "Varcircle patch "//trim(iStr)//": x_centroid must be set") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "Varcircle patch "//trim(iStr)//": y_centroid must be set") + @:PROHIBIT(patch_icpp(patch_id)%radius <= 0._wp, "Varcircle patch "//trim(iStr)//": radius must be greater than zero") + @:PROHIBIT(patch_icpp(patch_id)%epsilon <= 0._wp, "Varcircle patch "//trim(iStr)//": epsilon must be greater than zero") + + end subroutine s_check_varcircle_patch_geometry + + !> Check the 3D varcircle patch input; epsilon is the ring thickness + impure subroutine s_check_3d_varcircle_patch_geometry(patch_id) + + integer, intent(in) :: patch_id + + call s_int_to_str(patch_id, iStr) + + @:PROHIBIT(p == 0, "3D varcircle patch "//trim(iStr)//": p must be greater than zero") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "3D varcircle patch "//trim(iStr)//": x_centroid must be set") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "3D varcircle patch "//trim(iStr)//": y_centroid must be set") + @:PROHIBIT(f_is_default(patch_icpp(patch_id)%z_centroid), "3D varcircle patch "//trim(iStr)//": z_centroid must be set") + @:PROHIBIT(patch_icpp(patch_id)%length_z <= 0._wp, & + & "3D varcircle patch " // trim(iStr) // ": length_z must be greater than zero") + @:PROHIBIT(patch_icpp(patch_id)%radius <= 0._wp, & + & "3D varcircle patch " // trim(iStr) // ": radius must be greater than zero") + @:PROHIBIT(patch_icpp(patch_id)%epsilon <= 0._wp, & + & "3D varcircle patch " // trim(iStr) // ": epsilon must be greater than zero") + + end subroutine s_check_3d_varcircle_patch_geometry + !> Check the model patch input impure subroutine s_check_sphere_patch_geometry(patch_id) From b81c71a93e6b99f82e4cffc13b32672942c71ffc Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 10 Oct 2026 01:01:30 -0400 Subject: [PATCH 3/4] Check every IB patch slot in s_check_ib_patches The loop ran to num_patches_max (10, the ICPP patch limit) although patch_ib holds num_ib_patches_max_namelist entries, so IB patches past the tenth were never validated and inactive slots past it were never checked for stray settings. Loop over all patch_ib slots; unused slots keep their defaults, so the inactive-patch checks do not fire. Also set iStr before branching so the inactive-patch message names the right patch instead of a stale index. Verified with an 11-circle case whose 11th IB has no radius: pre_process now aborts with "in circle IB patch 11" (master accepted it), the same case with a radius passes, and the 61 IBM tests pass. (cherry picked from commit 1c0423e175a5ddf4fed5b301cc12b64792d9df96) --- src/pre_process/m_check_ib_patches.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pre_process/m_check_ib_patches.fpp b/src/pre_process/m_check_ib_patches.fpp index 9d037c42d7..5035aba6ef 100644 --- a/src/pre_process/m_check_ib_patches.fpp +++ b/src/pre_process/m_check_ib_patches.fpp @@ -34,9 +34,9 @@ contains integer :: i - do i = 1, num_patches_max + do i = 1, num_ib_patches_max_namelist + call s_int_to_str(i, iStr) if (i <= num_ibs) then - call s_int_to_str(i, iStr) @:PROHIBIT(patch_ib(i)%geometry == dflt_int, "IB patch undefined. patch_ib("//trim(iStr)//")%geometry must be set.") ! Constraints on the geometric initial condition patch parameters From d2822b0711509829c308e3b3a19e0c467940676a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 10 Oct 2026 00:23:38 -0400 Subject: [PATCH 4/4] Restrict sim_data to 3D two-fluid 5/6-equation cases post_process's sim_data writers are hard-wired to one configuration: s_write_energy_data_file multiplies by dz (unallocated in 1D/2D), reads adv(2) and alpha_rho(2) and the volume fractions at eqn_idx%E + 1 and eqn_idx%E + 2, and s_write_intf_data_file thresholds q_prim_vf(eqn_idx%E + 2). Any other case read unallocated memory or the wrong variables. Reject sim_data unless p > 0, num_fluids = 2 and model_eqns is 2 or 3, and cover the check in test_case_validator.py. (cherry picked from commit fa17c3a3f10a5fff30dd9d47b96a0ae6b3132c5e) --- toolchain/mfc/case_validator.py | 11 +++++++ toolchain/mfc/lint_docs.py | 1 + toolchain/mfc/test_case_validator.py | 47 ++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index befe499f30..971ac0b7a7 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -2673,6 +2673,16 @@ def check_surface_tension_post(self): self.prohibit(cf_wrt and not surface_tension, "cf_wrt can only be enabled if surface_tension is enabled") + def check_sim_data(self): + """Checks interface/energy data output constraints (post-process)""" + if self.get("sim_data", "F") != "T": + return + + # s_write_energy_data_file uses dz and s_write_intf_data_file reads alpha_2 at eqn_idx%E + 2 + self.prohibit(self.get("p", 0) == 0, "sim_data requires p > 0 (3D)") + self.prohibit(self.get("num_fluids") != 2, "sim_data requires num_fluids = 2") + self.prohibit(self.get("model_eqns") not in [2, 3], "sim_data requires model_eqns = 2 or 3") + def check_no_flow_variables(self): """Checks that at least one flow variable is selected (post-process)""" rho_wrt = self.get("rho_wrt", "F") == "T" @@ -3110,6 +3120,7 @@ def validate_post_process(self): self.check_liutex_post() self.check_schlieren() self.check_surface_tension_post() + self.check_sim_data() self.check_no_flow_variables() def validate(self, stage: str = "simulation"): diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py index 6eff9185da..65f43c9866 100644 --- a/toolchain/mfc/lint_docs.py +++ b/toolchain/mfc/lint_docs.py @@ -446,6 +446,7 @@ def check_physics_docs_coverage(repo_root: Path) -> list[str]: "check_partial_domain", # output format settings "check_perturb_density", # parameter pairing validation "check_qm", # output dimension requirements + "check_sim_data", # output dimension requirements "check_chemistry", # runtime Cantera validation only # Awaiting proper physics documentation (math, references, explanation) "check_adaptive_time_stepping", diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index 1000ffbbfa..28df8e67c6 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -719,3 +719,50 @@ def test_accepts_non_characteristic_bc_with_muscl(self): def test_accepts_cbc_with_weno(self): self.assertAccepts({**BASE, "bc_x%beg": -6, "bc_x%end": -6}) + + +class TestSimData(ConstraintTestCase): + """sim_data's interface and energy writers assume a 3D two-fluid 5- or 6-equation case.""" + + GOOD = { + **BASE_2D, + **{k: v for k, v in TWO_FLUID.items() if k not in BASE}, + "num_fluids": 2, + "p": 50, + "bc_z%beg": -1, + "bc_z%end": -1, + "z_domain%beg": 0.0, + "z_domain%end": 1.0, + "patch_icpp(1)%z_centroid": 0.5, + "patch_icpp(1)%length_z": 1.0, + "patch_icpp(1)%vel(3)": 0.0, + "rho_wrt": "T", + "sim_data": "T", + } + + def errors_for(self, params) -> str: + """Validate at the post_process stage, where sim_data is read.""" + validator = CaseValidator(dict(params)) + try: + validator.validate("post_process") + except CaseConstraintError as exc: + return str(exc) + return "" + + def test_accepts_3d_two_fluid(self): + self.assertAccepts(self.GOOD) + + def test_accepts_six_equation(self): + self.assertAccepts({**self.GOOD, "model_eqns": 3}) + + def test_rejects_2d(self): + self.assertRejects({**self.GOOD, "p": 0}, "sim_data requires p > 0 (3D)") + + def test_rejects_one_fluid(self): + self.assertRejects({**self.GOOD, "num_fluids": 1}, "sim_data requires num_fluids = 2") + + def test_rejects_gamma_law(self): + self.assertRejects({**self.GOOD, "model_eqns": 1}, "sim_data requires model_eqns = 2 or 3") + + def test_not_tripped_when_disabled(self): + self.assertAccepts({**BASE_2D, "rho_wrt": "T", "sim_data": "F"})