From 541de4a1bb67390450a5836e51c0a12042e09eba Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:17:17 -0500 Subject: [PATCH 01/10] 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 1db57b99570174ad134393b95db1926ecdbfe56a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:18:06 -0500 Subject: [PATCH 02/10] 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 e012f54273142d30d5d929d98b0e733d454ca86a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:20:05 -0500 Subject: [PATCH 03/10] 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 a8e0cc8932babff80e3b43371e6062409a934b74 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:02:35 -0500 Subject: [PATCH 04/10] 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 08e6ce7794b543b127eb157cfac660ae18ff6756 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:02:40 -0500 Subject: [PATCH 05/10] 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 0daa425bc7362106ac25740b0eb2ed1748cc833f Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 00:11:44 -0500 Subject: [PATCH 06/10] 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 b63f9e23e4a6c78a602be02bd0d5771bc5bfab61 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 21 Feb 2026 01:13:27 -0500 Subject: [PATCH 07/10] 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 cfed440fa985321cb28641cff10cbaf79d2a3a33 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 20 Feb 2026 22:20:41 -0500 Subject: [PATCH 08/10] 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 ea6715dbf600f66dfeaf0ded584cef780b04f698 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 09:19:06 -0500 Subject: [PATCH 09/10] 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 9276e0bbc9cda971a67d7c2d8dd85d7cbba56aa4 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Mon, 23 Feb 2026 09:42:00 -0500 Subject: [PATCH 10/10] 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