From f64512790d805a4d300f947a5ad0adabf004c2f0 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:17:17 -0500 Subject: [PATCH 01/13] Fix moncon_cutoff declared as integer, truncating 1e-8 to 0 moncon_cutoff is assigned 1e-8_wp but declared as integer, so Fortran silently truncates it to 0. This makes all THINC monotonicity constraint comparisons always true, completely disabling the constraint. Co-Authored-By: Claude Opus 4.6 --- src/common/m_constants.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index f926824cb8..1a87f0a912 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -44,7 +44,7 @@ module m_constants ! Interface Compression real(wp), parameter :: dflt_ic_eps = 1e-4_wp !< Ensure compression is only applied to surface cells in THINC real(wp), parameter :: dflt_ic_beta = 1.6_wp !< Sharpness parameter's default value used in THINC - integer, parameter :: moncon_cutoff = 1e-8_wp !< Monotonicity constraint's limiter to prevent extremas in THINC + real(wp), parameter :: moncon_cutoff = 1e-8_wp !< Monotonicity constraint's limiter to prevent extremas in THINC ! Chemistry real(wp), parameter :: dflt_T_guess = 1200._wp ! Default guess for temperature (when a previous value is not available) From ad4667d95cac30a7619877dcb2d771ccfc12de37 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:18:06 -0500 Subject: [PATCH 02/13] Fix grid stretching using wrong array bounds for y_cc and z_cc y_cc(0:m) and z_cc(0:m) use the x-dimension size m instead of the y-dimension n and z-dimension p respectively. Corrupts cell-center coordinates when grid stretching is enabled and m != n or m != p. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_grid.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pre_process/m_grid.f90 b/src/pre_process/m_grid.f90 index 92bf4aa052..6ea71bd592 100644 --- a/src/pre_process/m_grid.f90 +++ b/src/pre_process/m_grid.f90 @@ -131,7 +131,7 @@ impure subroutine s_generate_serial_grid end do y_cb = y_cb*length - y_cc(0:m) = (y_cb(0:n) + y_cb(-1:n - 1))/2._wp + y_cc(0:n) = (y_cb(0:n) + y_cb(-1:n - 1))/2._wp dy = minval(y_cb(0:n) - y_cb(-1:n - 1)) @@ -168,7 +168,7 @@ impure subroutine s_generate_serial_grid end do z_cb = z_cb*length - z_cc(0:m) = (z_cb(0:p) + z_cb(-1:p - 1))/2._wp + z_cc(0:p) = (z_cb(0:p) + z_cb(-1:p - 1))/2._wp dz = minval(z_cb(0:p) - z_cb(-1:p - 1)) From 6bd9bf7964cae215dc3bd318ccd3f7758ae91d37 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:20:05 -0500 Subject: [PATCH 03/13] Remove duplicate R3bar accumulation line that doubles bubble perturbation Character-for-character identical duplicate line causes R3bar to always be 2x the intended value for QBMM cases. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_assign_variables.fpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index e41428bf90..65809485ba 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -233,7 +233,6 @@ contains if (qbmm) then do i = 1, nb R3bar = R3bar + weight(i)*0.5_wp*(q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l))**3._wp - R3bar = R3bar + weight(i)*0.5_wp*(q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l))**3._wp end do else do i = 1, nb From 0a7b4a391a0c4dbb5fddfffa0f785c5e0fb28c25 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:02:35 -0500 Subject: [PATCH 04/13] Fix y-velocity perturbation using already-modified x-velocity Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_perturbation.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pre_process/m_perturbation.fpp b/src/pre_process/m_perturbation.fpp index 18fb89cd34..91b1f339ed 100644 --- a/src/pre_process/m_perturbation.fpp +++ b/src/pre_process/m_perturbation.fpp @@ -83,8 +83,8 @@ contains perturb_alpha = q_prim_vf(E_idx + perturb_flow_fluid)%sf(i, j, k) call random_number(rand_real) rand_real = rand_real*perturb_flow_mag - q_prim_vf(mom_idx%beg)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(mom_idx%beg)%sf(i, j, k) q_prim_vf(mom_idx%end)%sf(i, j, k) = rand_real*q_prim_vf(mom_idx%beg)%sf(i, j, k) + q_prim_vf(mom_idx%beg)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(mom_idx%beg)%sf(i, j, k) if (bubbles_euler) then q_prim_vf(alf_idx)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(alf_idx)%sf(i, j, k) end if From 680716e6415519740caafd7cf43a61a1e5fc92ce Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:02:40 -0500 Subject: [PATCH 05/13] Fix hyper_cleaning psi initialization only covering k=0 plane Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_start_up.fpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 3b45547c29..b631211ef5 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -767,7 +767,7 @@ contains real(wp), intent(inout) :: start, finish - integer :: j, k + integer :: j, k, l ! Setting up the grid and the initial condition. If the grid is read in from ! preexisting grid data files, it is checked for consistency. If the grid is @@ -787,10 +787,12 @@ contains ! hard-coded psi if (hyper_cleaning) then - do j = 0, m + do l = 0, p do k = 0, n - q_cons_vf(psi_idx)%sf(j, k, 0) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2)/(2.0*0.05**2)) - q_prim_vf(psi_idx)%sf(j, k, 0) = q_cons_vf(psi_idx)%sf(j, k, 0) + do j = 0, m + q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0*0.05**2)) + q_prim_vf(psi_idx)%sf(j, k, l) = q_cons_vf(psi_idx)%sf(j, k, l) + end do end do end do end if From dbbe35b50f0fd637989477d82a2462a6dfdf6bcb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:11:44 -0500 Subject: [PATCH 06/13] Guard z_cc access in hyper_cleaning init for 2D case (p=0) z_cc is only allocated when p > 0. The previous commit accessed z_cc(l) unconditionally, which would crash or read garbage in 2D runs. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_start_up.fpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index b631211ef5..f02dc27f65 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -790,7 +790,11 @@ contains do l = 0, p do k = 0, n do j = 0, m - q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0*0.05**2)) + if (p > 0) then + q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0*0.05**2)) + else + q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2)/(2.0*0.05**2)) + end if q_prim_vf(psi_idx)%sf(j, k, l) = q_cons_vf(psi_idx)%sf(j, k, l) end do end do From 9ab9c37f272765b9a7f1d491993f09b1febac7d0 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 01:13:27 -0500 Subject: [PATCH 07/13] Fix mixed-precision literals in hyper_cleaning psi initialization Replace `1d-2` (hard double precision) and bare `2.0`/`0.05` (default real) with `_wp`-suffixed literals so the Gaussian initialization is consistent with the working precision across single and double builds. Co-Authored-By: Claude Sonnet 4.6 --- src/pre_process/m_start_up.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index f02dc27f65..1f349ab139 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -791,9 +791,9 @@ contains do k = 0, n do j = 0, m if (p > 0) then - q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0*0.05**2)) + q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0_wp*0.05_wp**2)) else - q_cons_vf(psi_idx)%sf(j, k, l) = 1d-2*exp(-(x_cc(j)**2 + y_cc(k)**2)/(2.0*0.05**2)) + q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2)/(2.0_wp*0.05_wp**2)) end if q_prim_vf(psi_idx)%sf(j, k, l) = q_cons_vf(psi_idx)%sf(j, k, l) end do From 097f6c7757be1b46170f57d9a1dde15525287bc6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:20:41 -0500 Subject: [PATCH 08/13] Fix IB patch vel/angular_vel/angles broadcast inside wrong loop patch_ib vel, angular_vel, and angles were broadcast inside the num_bc_patches_max loop instead of the num_patches_max loop where the other patch_ib fields are broadcast. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_mpi_proxy.fpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 30ef061689..0a074ba616 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -76,10 +76,6 @@ contains call MPI_BCAST(patch_bc(i)%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor - #:for VAR in ['vel', 'angular_vel', 'angles'] - call MPI_BCAST(patch_ib(i)%${VAR}$, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - #:endfor - call MPI_BCAST(patch_bc(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:for VAR in ['centroid', 'length'] @@ -121,6 +117,9 @@ contains call MPI_BCAST(patch_icpp(i)%Y, size(patch_icpp(i)%Y), mpi_p, 0, MPI_COMM_WORLD, ierr) end if ! Broadcast IB variables + #:for VAR in ['vel', 'angular_vel', 'angles'] + call MPI_BCAST(patch_ib(i)%${VAR}$, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) + #:endfor call MPI_BCAST(patch_ib(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(patch_ib(i)%model_filepath, len(patch_ib(i)%model_filepath), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(patch_ib(i)%model_threshold, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) From 9d7d4abfe066e3891096e505e48423957a094d0b Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 09:19:06 -0500 Subject: [PATCH 09/13] Simplify hyper_cleaning init and use size() for IB broadcast - Remove redundant if(p>0) branch in psi initialization: z_cc(0)=0 in 2D, so the 3D formula produces identical results in both cases. - Replace hardcoded 3 with size(patch_ib(i)%VAR) in MPI_BCAST for vel/angular_vel/angles, consistent with other IB broadcasts. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_mpi_proxy.fpp | 2 +- src/pre_process/m_start_up.fpp | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 0a074ba616..d88981b019 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -118,7 +118,7 @@ contains end if ! Broadcast IB variables #:for VAR in ['vel', 'angular_vel', 'angles'] - call MPI_BCAST(patch_ib(i)%${VAR}$, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(patch_ib(i)%${VAR}$, size(patch_ib(i)%${VAR}$), mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(patch_ib(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(patch_ib(i)%model_filepath, len(patch_ib(i)%model_filepath), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 1f349ab139..115b893237 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -790,11 +790,7 @@ contains do l = 0, p do k = 0, n do j = 0, m - if (p > 0) then - q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0_wp*0.05_wp**2)) - else - q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2)/(2.0_wp*0.05_wp**2)) - end if + q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0_wp*0.05_wp**2)) q_prim_vf(psi_idx)%sf(j, k, l) = q_cons_vf(psi_idx)%sf(j, k, l) end do end do From 22304a19ab17c514365e2504849b1b764819aecc Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 09:42:00 -0500 Subject: [PATCH 10/13] Guard y_cc and z_cc access in hyper_cleaning init for 1D/2D safety Build the Gaussian r^2 incrementally, guarding y_cc(k) behind n > 0 and z_cc(l) behind p > 0, since these arrays are only allocated in 2D+ and 3D respectively. Prevents undefined behavior if hyper_cleaning is enabled in a 1D MHD configuration. Co-Authored-By: Claude Opus 4.6 --- src/pre_process/m_start_up.fpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 115b893237..259a56a6d7 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -768,6 +768,7 @@ contains real(wp), intent(inout) :: start, finish integer :: j, k, l + real(wp) :: r2 ! Setting up the grid and the initial condition. If the grid is read in from ! preexisting grid data files, it is checked for consistency. If the grid is @@ -790,7 +791,10 @@ contains do l = 0, p do k = 0, n do j = 0, m - q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-(x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)/(2.0_wp*0.05_wp**2)) + r2 = x_cc(j)**2 + if (n > 0) r2 = r2 + y_cc(k)**2 + if (p > 0) r2 = r2 + z_cc(l)**2 + q_cons_vf(psi_idx)%sf(j, k, l) = 1.0e-2_wp*exp(-r2/(2.0_wp*0.05_wp**2)) q_prim_vf(psi_idx)%sf(j, k, l) = q_cons_vf(psi_idx)%sf(j, k, l) end do end do From a8f0d84771a3086431a75ed9d48071a5120b2ddb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 17:31:43 -0500 Subject: [PATCH 11/13] Prohibit hyper_cleaning in 2D simulations (p == 0) The psi field used by hyperbolic cleaning is initialized with a 3D loop that accesses z_cc, which is out-of-bounds when p == 0. The 1D prohibition already existed; extend it to 2D. Co-Authored-By: Claude Sonnet 4.6 --- toolchain/mfc/case_validator.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 2c29c75e7f..39b7209831 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1029,6 +1029,7 @@ def check_mhd_simulation(self): relativity = self.get('relativity', 'F') == 'T' hyper_cleaning = self.get('hyper_cleaning', 'F') == 'T' n = self.get('n', 0) + p = self.get('p', 0) self.prohibit(mhd and riemann_solver is not None and riemann_solver not in [1, 4], "MHD simulations require riemann_solver = 1 (HLL) or riemann_solver = 4 (HLLD)") @@ -1040,6 +1041,8 @@ def check_mhd_simulation(self): "Hyperbolic cleaning requires mhd to be enabled") self.prohibit(hyper_cleaning and n is not None and n == 0, "Hyperbolic cleaning is not supported for 1D simulations") + self.prohibit(hyper_cleaning and p is not None and p == 0, + "Hyperbolic cleaning is not supported for 2D simulations") def check_igr_simulation(self): # pylint: disable=too-many-locals From 4997be7769ae91a639e389e7182d04013418fd5a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 18:16:44 -0500 Subject: [PATCH 12/13] Fix hyper_cleaning 2D prohibition double-triggering for 1D cases In 1D simulations n==0 and p==0, so both the 1D and 2D prohibition checks were firing simultaneously. Gate the 2D check on n > 0 so it only triggers for configurations that are actually 2D. Co-Authored-By: Claude Sonnet 4.6 --- toolchain/mfc/case_validator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 39b7209831..a34d4d242d 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1041,7 +1041,7 @@ def check_mhd_simulation(self): "Hyperbolic cleaning requires mhd to be enabled") self.prohibit(hyper_cleaning and n is not None and n == 0, "Hyperbolic cleaning is not supported for 1D simulations") - self.prohibit(hyper_cleaning and p is not None and p == 0, + self.prohibit(hyper_cleaning and n is not None and n > 0 and p is not None and p == 0, "Hyperbolic cleaning is not supported for 2D simulations") From d27fe4b27b2812f61cf0b05e39e5e8d547262f19 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 25 Feb 2026 09:36:16 -0500 Subject: [PATCH 13/13] Remove incorrect 2D hyper_cleaning prohibition GLM divergence cleaning is valid in 2D. The z_cc access was already guarded with `if (p > 0)` in 22304a19. The prohibition broke the existing 2D hyper_cleaning, mhd_rotor, and orszag_tang tests. Co-Authored-By: Claude Opus 4.6 --- toolchain/mfc/case_validator.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index a34d4d242d..2c29c75e7f 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1029,7 +1029,6 @@ def check_mhd_simulation(self): relativity = self.get('relativity', 'F') == 'T' hyper_cleaning = self.get('hyper_cleaning', 'F') == 'T' n = self.get('n', 0) - p = self.get('p', 0) self.prohibit(mhd and riemann_solver is not None and riemann_solver not in [1, 4], "MHD simulations require riemann_solver = 1 (HLL) or riemann_solver = 4 (HLLD)") @@ -1041,8 +1040,6 @@ def check_mhd_simulation(self): "Hyperbolic cleaning requires mhd to be enabled") self.prohibit(hyper_cleaning and n is not None and n == 0, "Hyperbolic cleaning is not supported for 1D simulations") - self.prohibit(hyper_cleaning and n is not None and n > 0 and p is not None and p == 0, - "Hyperbolic cleaning is not supported for 2D simulations") def check_igr_simulation(self): # pylint: disable=too-many-locals