Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
a838583
hyperbolic cleaning
ChrisZYJ Apr 24, 2025
c745398
flux per dim (wrong)
ChrisZYJ Jun 8, 2025
e98fb01
hyper works
ChrisZYJ Jun 9, 2025
6798e4b
format
ChrisZYJ Jun 9, 2025
4d651e8
hardcoded cases & m_start_up hard-coded psi
ChrisZYJ Jul 1, 2025
89a0ac5
Added ellipse marker patch
danieljvickers Dec 12, 2025
ef446f5
Added ellipse IB patch
danieljvickers Dec 12, 2025
2f9b1be
Added an ellipse immersed boundary patch
danieljvickers Dec 12, 2025
f0b7239
validation cases
ChrisZYJ Dec 13, 2025
2c29ab2
Merged master
danieljvickers Jan 23, 2026
7db99b4
Missed one somehow
danieljvickers Jan 23, 2026
f3b6a37
Resolved another commit
danieljvickers Jan 23, 2026
32655d0
Builds!
danieljvickers Jan 23, 2026
edc9f5d
Merge branch 'master' into mhd_hypercleaning
danieljvickers Jan 23, 2026
64f4d96
Formatting
danieljvickers Jan 23, 2026
96ccf10
Remove duplicate values from user_inputs namelist
danieljvickers Jan 24, 2026
f942bb9
Missing comma
danieljvickers Jan 24, 2026
4480d18
Replaced outdated ACC calls with unified macros
danieljvickers Jan 24, 2026
c642669
Adjust comment syntax to make the linter happy
danieljvickers Jan 24, 2026
bd12ee0
Merge branch 'master' into mhd_hypercleaning
danieljvickers Jan 24, 2026
efb124d
Updates to the rotor problem to match original literature
danieljvickers Jan 25, 2026
3011794
Merge branch 'mhd_hypercleaning' of github.com:ChrisZYJ/MFC_PR into m…
danieljvickers Jan 26, 2026
6c61e02
Fixed the issue. Had to change how we damped psi
danieljvickers Jan 26, 2026
780f9bd
fix hyper cleaning
ChrisZYJ Jan 27, 2026
f54854b
remove powell; docs; golden for hyper_cleaning
ChrisZYJ Jan 27, 2026
4e25edb
format; add examples
ChrisZYJ Jan 27, 2026
178a9bb
missed on golden file
danieljvickers Jan 27, 2026
f1862ed
Now deleting analysis cases that won't be added to examples
danieljvickers Jan 27, 2026
3503bb5
Domain boundaries and endings need to be floats or Cray compilers wil…
danieljvickers Jan 27, 2026
059b941
RHS updates need to be performed on the GPU. Also deleted block comme…
danieljvickers Jan 27, 2026
9e04bd7
formatting
danieljvickers Jan 27, 2026
ba275e8
Merge branch 'master' into mhd_hypercleaning
sbryngelson Feb 2, 2026
ba7b4af
Fixing spelling change that got merged in from master
danieljvickers Feb 2, 2026
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .typos.toml
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ TKE = "TKE"
HSA = "HSA"
infp = "infp"
Sur = "Sur"
iz = "iz"

[files]
extend-exclude = ["docs/documentation/references*", "tests/", "toolchain/cce_simulation_workgroup_256.sh"]
23 changes: 13 additions & 10 deletions docs/documentation/case.md
Original file line number Diff line number Diff line change
Expand Up @@ -939,22 +939,25 @@ By convention, positive accelerations in the `x[y,z]` direction are in the posit

### 14. Magnetohydrodynamics (MHD)

| Parameter | Type | Description |
| ---: | :---: | :--- |
| `mhd` | Logical | Enable ideal MHD simulation |
| `relativity` | Logical | Enable relativistic MHD simulation |
| `powell` | Logical | Enable Powell's method for solenoidal constraint |
| `fd_order` | Integer | Finite difference order for Powell's method |
| `Bx[y,z]` | Real | Initial magnetic field in the x[y,z] direction |
| `Bx0` | Real | Constant magnetic field in the x direction (1D only)|
| Parameter | Type | Description |
| ---: | :---: | :--- |
| `mhd` | Logical | Enable ideal MHD simulation |
| `relativity` | Logical | Enable relativistic MHD simulation |
| `hyper_cleaning` | Logical | Enable hyperbolic (GLM) divergence cleaning for `div B` |
| `hyper_cleaning_speed` | Real | Cleaning wave speed `c_h` |
| `hyper_cleaning_tau` | Real | Cleaning damping timescale `tau` |
| `Bx[y,z]` | Real | Initial magnetic field in the x[y,z] direction |
| `Bx0` | Real | Constant magnetic field in the x direction (1D only) |

- `mhd` is currently only available for single-component flows and 5-equation model. Its compatibility with most other features is work in progress.

- `relativity` only works for `mhd` enabled and activates relativistic MHD (RMHD) simulation.

- `powell` activates Powell's eight-wave method to impose the solenoidal constraint in the MHD simulation [Powell (1994)](references.md). It should not be used in conjunction with HLLD (`riemann_solver = 4`).
- `hyper_cleaning` [Dedner et al., 2002](references.md) only works with `mhd` in 2D/3D and reduces numerical `div B` errors by propagation and damping. Currently not compatible with HLLD (`riemann_solver = 4`).

- `hyper_cleaning_speed` sets the propagation speed of divergence-cleaning waves.

- `fd_order` specifies the finite difference order for computing the RHS of the Powell's method. `fd_order = 1`, `2`, and `4` are allowed.
- `hyper_cleaning_tau` sets the decay timescale for divergence-cleaning.

- `Bx0` is only used in 1D simulations to specify the constant magnetic field in the x direction. It must be specified in 1D simulations. `Bx` must not be used in 1D simulations.

Expand Down
2 changes: 1 addition & 1 deletion docs/documentation/references.md
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@

- <a id="Miyoshi05">Miyoishi, T., and Kusano, K. (2005). A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208(1), 315-344.</a>

- <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>
- <a id="Dedner02">Dedner, A., Kemm, F., Kröner, D., Munz, C. D., Schnitzer, T., & Wesenberg, M. (2002). Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics, 175(2), 645-673.</a>

- <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>

Expand Down
2 changes: 1 addition & 1 deletion examples/2D_mhd_magnetic_vortex/case.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@
"parallel_io": "T",
# MHD
"mhd": "T",
"patch_icpp(1)%hcid": 252,
"patch_icpp(1)%hcid": 253,
"patch_icpp(1)%geometry": 3,
"patch_icpp(1)%x_centroid": 0.0,
"patch_icpp(1)%y_centroid": 0.0,
Expand Down
84 changes: 84 additions & 0 deletions examples/2D_mhd_rotor/case.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
#!/usr/bin/env python3
import json
import math

# Configuring case dictionary
print(
json.dumps(
{
# Logistics
"run_time_info": "T",
# Computational Domain Parameters
"x_domain%beg": 0.0,
"x_domain%end": 1.0,
"y_domain%beg": 0.0,
"y_domain%end": 1.0,
"m": 399,
"n": 399,
"p": 0,
"dt": 0.0004,
"t_step_start": 0,
"t_step_stop": 375,
"t_step_save": 5,
# Simulation Algorithm Parameters
"num_patches": 1,
"model_eqns": 2,
"alt_soundspeed": "F",
"num_fluids": 1,
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": 3,
"weno_order": 5,
"mapped_weno": "T",
"weno_eps": 1.0e-6,
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": 1,
"wave_speeds": 1,
"avg_state": 2,
"bc_x%beg": -1,
"bc_x%end": -1,
"bc_y%beg": -1,
"bc_y%end": -1,
# Formatted Database Files Structure Parameters
"format": 1,
"precision": 2,
"prim_vars_wrt": "T",
"rho_wrt": "T",
"parallel_io": "T",
# MHD
"mhd": "T",
"hyper_cleaning": "T",
"hyper_cleaning_speed": 2.5,
"hyper_cleaning_tau": 0.004,
# Patch 1 - 2D MHD Rotor Problem
# gamma = 1.4
# Ambient medium (r > 0.1):
# rho = 1
# p = 1
# v = (0, 0, 0)
# B = (1, 0, 0)
# Rotor (r <= 0.1):
# rho = 10
# v has angular velocity of 20
"patch_icpp(1)%hcid": 252,
"patch_icpp(1)%geometry": 3,
"patch_icpp(1)%x_centroid": 0.5,
"patch_icpp(1)%y_centroid": 0.5,
"patch_icpp(1)%length_x": 1.0,
"patch_icpp(1)%length_y": 1.0,
"patch_icpp(1)%vel(1)": 0.0,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%vel(3)": 0.0,
"patch_icpp(1)%pres": 1.0,
"patch_icpp(1)%Bx": 5.0 / math.sqrt(math.pi * 4.0),
"patch_icpp(1)%By": 0.0,
"patch_icpp(1)%Bz": 0.0,
"patch_icpp(1)%alpha_rho(1)": 1.0,
"patch_icpp(1)%alpha(1)": 1.0,
# Fluids Physical Parameters
"fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
"fluid_pp(1)%pi_inf": 0.0,
}
)
)
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@
"m": 512,
"n": 512,
"p": 0,
"dt": 0.0005,
"dt": 0.0004,
"t_step_start": 0,
"t_step_stop": 2000,
"t_step_stop": 2500,
"t_step_save": 100,
# Simulation Algorithm Parameters
"num_patches": 1,
Expand All @@ -28,7 +28,8 @@
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": 3,
"weno_order": 1,
"weno_order": 5,
"mapped_weno": "T",
"weno_eps": 1.0e-6,
"null_weights": "F",
"mp_weno": "F",
Expand All @@ -47,8 +48,9 @@
"parallel_io": "T",
# MHD
"mhd": "T",
"powell": "T",
"fd_order": 2,
"hyper_cleaning": "T",
"hyper_cleaning_speed": 2.5,
"hyper_cleaning_tau": 0.004,
# Patch 1 - Analytical for v and B
# gamma = 5/3
# rho = 25/(36*pi)
Expand Down
101 changes: 99 additions & 2 deletions src/common/include/2dHardcodedIC.fpp
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

High-level Suggestion

Refactor the new hardcoded initial conditions in Fortran files like 2dHardcodedIC.fpp into configurable input files for the test cases. This improves modularity and separates test configuration from core logic. [High-level, importance: 6]

Solution Walkthrough:

Before:

! File: src/common/include/2dHardcodedIC.fpp
...
case (251)
    ...
! case 252 is for the 2D MHD Rotor problem
case (252) ! 2D MHD Rotor Problem
    r_sq = (x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2
    if (r_sq <= 0.1**2) then
        q_prim_vf(contxb)%sf(i, j, 0) = 10._wp
        ...
    else if (r_sq <= 0.115**2) then
        ...
    end if
case (253) ! MHD Smooth Magnetic Vortex
    ...
case (260) ! Gaussian Divergence Pulse
    ...

After:

! File: src/common/include/2dHardcodedIC.fpp
...
case (251)
    ...
! case (252) ! MHD Smooth Magnetic Vortex (ID shifted)
!    ...
! Other hardcoded cases removed and logic moved to python case files.
...

! File: examples/2D_mhd_rotor/case.py
...
# Instead of hcid, use python functions to define ICs
# and pass them to the simulation, for example via a file
# that the Fortran code can read.
# (This framework does not seem to support python-defined ICs directly,
# so the best practice would be to generalize the existing IC system
# if possible, rather than adding more hardcoded cases.)

Original file line number Diff line number Diff line change
@@ -1,9 +1,14 @@
#:def Hardcoded2DVariables()
! Place any declaration of intermediate variables here
real(wp) :: eps
real(wp) :: eps, eps_mhd, C_mhd
real(wp) :: r, rmax, gam, umax, p0
real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph
real(wp) :: factor
real(wp) :: r0, alpha, r2
real(wp) :: sinA, cosA

real(wp) :: r_sq

! # 207
real(wp) :: sigma, gauss1, gauss2
! # 208
Expand Down Expand Up @@ -183,7 +188,43 @@
q_prim_vf(E_idx)%sf(i, j, 0) = 3.e-5_wp
end if

case (252) ! MHD Smooth Magnetic Vortex
! case 252 is for the 2D MHD Rotor problem
case (252) ! 2D MHD Rotor Problem
! Ambient conditions are set in the JSON file.
! This case imposes the dense, rotating cylinder.
!
! gamma = 1.4
! Ambient medium (r > 0.1):
! rho = 1, p = 1, v = 0, B = (1,0,0)
! Rotor (r <= 0.1):
! rho = 10, p = 1
! v has angular velocity w=20, giving v_tan=2 at r=0.1

! Calculate distance squared from the center
r_sq = (x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2

! inner radius of 0.1
if (r_sq <= 0.1**2) then
! -- Inside the rotor --
! Set density uniformly to 10
q_prim_vf(contxb)%sf(i, j, 0) = 10._wp

! Set vup constant rotation of rate v=2
! v_x = -omega * (y - y_c)
! v_y = omega * (x - x_c)
q_prim_vf(momxb)%sf(i, j, 0) = -20._wp*(y_cc(j) - 0.5_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = 20._wp*(x_cc(i) - 0.5_wp)

! taper width of 0.015
else if (r_sq <= 0.115**2) then
Comment on lines +207 to +219
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: In the rotor initial condition, the comparison radii use default-kind real literals (0.1**2, 0.115**2) while r_sq is real(wp), causing mixed‑kind arithmetic and slight precision inconsistencies that can misclassify cells very close to the intended radius thresholds. Use _wp suffixed literals to keep everything in the same kind. [type error]

Severity Level: Major ⚠️
- ⚠️ Rotor initial condition (hcid=252) borderline cell classification.
- ⚠️ Small differences in initial rotor velocity/density assignments.
Suggested change
if (r_sq <= 0.1**2) then
! -- Inside the rotor --
! Set density uniformly to 10
q_prim_vf(contxb)%sf(i, j, 0) = 10._wp
! Set vup constant rotation of rate v=2
! v_x = -omega * (y - y_c)
! v_y = omega * (x - x_c)
q_prim_vf(momxb)%sf(i, j, 0) = -20._wp*(y_cc(j) - 0.5_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = 20._wp*(x_cc(i) - 0.5_wp)
! taper width of 0.015
else if (r_sq <= 0.115**2) then
if (r_sq <= 0.1_wp**2) then
! -- Inside the rotor --
! Set density uniformly to 10
q_prim_vf(contxb)%sf(i, j, 0) = 10._wp
! Set vup constant rotation of rate v=2
! v_x = -omega * (y - y_c)
! v_y = omega * (x - x_c)
q_prim_vf(momxb)%sf(i, j, 0) = -20._wp*(y_cc(j) - 0.5_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = 20._wp*(x_cc(i) - 0.5_wp)
! taper width of 0.015
else if (r_sq <= 0.115_wp**2) then
Steps of Reproduction ✅
1. Configure a patch with hardcoded initial condition id hcid=252 (2D MHD Rotor) in the
input (Hardcoded2D is defined in src/common/include/2dHardcodedIC.fpp).

2. Run the solver to reach initialization where the Hardcoded2D macro executes; the rotor
branch starts at src/common/include/2dHardcodedIC.fpp lines around 192..227 and computes
r_sq at line 204: "r_sq = (x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2".

3. Initialization evaluates the radius tests at lines 207 ("if (r_sq <= 0.1**2) then") and
219 ("else if (r_sq <= 0.115**2) then") to classify inside/transition/outside cells.

4. Because the literals 0.1**2 and 0.115**2 are default‑kind reals while r_sq and other
variables use kind _wp, cells with sqrt(r_sq) numerically extremely close to 0.1 or 0.115
can be classified differently than when using _wp‑suffixed literals; this reproduces the
misclassification at initialization for boundary grid cells. (This is a marginal numerical
precision issue localized to hcid=252; the code path and lines above show the exact
location.)
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/include/2dHardcodedIC.fpp
**Line:** 207:219
**Comment:**
	*Type Error: In the rotor initial condition, the comparison radii use default-kind real literals (`0.1**2`, `0.115**2`) while `r_sq` is `real(wp)`, causing mixed‑kind arithmetic and slight precision inconsistencies that can misclassify cells very close to the intended radius thresholds. Use `_wp` suffixed literals to keep everything in the same kind.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

! linearly smooth the function between r = 0.1 and 0.115
q_prim_vf(contxb)%sf(i, j, 0) = 1._wp + 9._wp*(0.115_wp - sqrt(r_sq))/(0.015_wp)

q_prim_vf(momxb)%sf(i, j, 0) = -(2._wp/sqrt(r_sq))*(y_cc(j) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = (2._wp/sqrt(r_sq))*(x_cc(i) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
end if
Comment on lines +192 to +225
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Add an else branch to the MHD Rotor problem to initialize the ambient region's density and momentum. [possible issue, importance: 8]

Suggested change
case (252) ! 2D MHD Rotor Problem
! Ambient conditions are set in the JSON file.
! This case imposes the dense, rotating cylinder.
!
! gamma = 1.4
! Ambient medium (r > 0.1):
! rho = 1, p = 1, v = 0, B = (1,0,0)
! Rotor (r <= 0.1):
! rho = 10, p = 1
! v has angular velocity w=20, giving v_tan=2 at r=0.1
! Calculate distance squared from the center
r_sq = (x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2
! inner radius of 0.1
if (r_sq <= 0.1**2) then
! -- Inside the rotor --
! Set density uniformly to 10
q_prim_vf(contxb)%sf(i, j, 0) = 10._wp
! Set vup constant rotation of rate v=2
! v_x = -omega * (y - y_c)
! v_y = omega * (x - x_c)
q_prim_vf(momxb)%sf(i, j, 0) = -20._wp*(y_cc(j) - 0.5_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = 20._wp*(x_cc(i) - 0.5_wp)
! taper width of 0.015
else if (r_sq <= 0.115**2) then
! linearly smooth the function between r = 0.1 and 0.115
q_prim_vf(contxb)%sf(i, j, 0) = 1._wp + 9._wp*(0.115_wp - sqrt(r_sq))/(0.015_wp)
q_prim_vf(momxb)%sf(i, j, 0) = -(2._wp/sqrt(r_sq))*(y_cc(j) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
q_prim_vf(momxb + 1)%sf(i, j, 0) = (2._wp/sqrt(r_sq))*(x_cc(i) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
end if
case (252) ! 2D MHD Rotor Problem
...
if (r_sq <= 0.1_wp**2) then
...
else if (r_sq <= 0.115_wp**2) then
...
else
! -- Ambient medium --
q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
q_prim_vf(momxb)%sf(i, j, 0) = 0._wp
q_prim_vf(momxb + 1)%sf(i, j, 0) = 0._wp
end if


case (253) ! MHD Smooth Magnetic Vortex
! Section 5.2 of
! Implicit hybridized discontinuous Galerkin methods for compressible magnetohydrodynamics
! C. Ciuca, P. Fernandez, A. Christophe, N.C. Nguyen, J. Peraire
Expand All @@ -199,6 +240,62 @@
! pressure
q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp + (1 - 2._wp*(x_cc(i)**2 + y_cc(j)**2))*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/((2._wp*pi)**3)

case (260) ! Gaussian Divergence Pulse
! Bx(x) = 1 + C * erf((x-0.5)/σ)
! ⇒ ∂Bx/∂x = C * (2/√π) * exp[-((x-0.5)/σ)**2] * (1/σ)
! Choose C = ε * σ * √π / 2 ⇒ ∂Bx/∂x = ε * exp[-((x-0.5)/σ)**2]
! ψ is initialized to zero everywhere.

eps_mhd = patch_icpp(patch_id)%a(2)
sigma = patch_icpp(patch_id)%a(3)
C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: In the Gaussian divergence pulse case, sigma is read directly from input and then used as a divisor without validation; if the user forgets to set it or sets it to a non-positive value, this will cause a division by zero or mathematically invalid configuration at runtime. [possible bug]

Severity Level: Major ⚠️
- ❌ Initialization for Gaussian divergence pulse (hcid=260) fails.
- ⚠️ Produces NaN/INF B‑field if sigma is zero or negative.
Suggested change
C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp
if (sigma <= 0._wp) then
call s_mpi_abort("Gaussian Divergence Pulse (hcid=260): sigma must be positive.")
end if
Steps of Reproduction ✅
1. Configure a patch to use hardcoded initial condition hcid=260 (Gaussian Divergence
Pulse) in the input so the case (260) branch in src/common/include/2dHardcodedIC.fpp is
executed.

2. Provide patch parameters with a(3)=0 or omit setting a(3) so patch_icpp(patch_id)%a(3)
resolves to zero; the code assigns sigma at src/common/include/2dHardcodedIC.fpp line 249:
"sigma = patch_icpp(patch_id)%a(3)".

3. Initialization computes the B‑field using erf((x_cc(i) - 0.5_wp)/sigma) at line 254.
With sigma == 0 the division produces an invalid argument (INF/NaN) or a runtime fault
depending on compiler/FP settings, observable as NaNs in B arrays or an initialization
failure.

4. Observe the failure by running the solver to initialization and inspecting the
initialized q_prim_vf(B_idx) array or checking for abnormal program termination. Adding
the explicit check (sigma <= 0) at line 251 and aborting with a clear message prevents the
invalid initialization and makes the failure mode explicit.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/include/2dHardcodedIC.fpp
**Line:** 251:251
**Comment:**
	*Possible Bug: In the Gaussian divergence pulse case, `sigma` is read directly from input and then used as a divisor without validation; if the user forgets to set it or sets it to a non-positive value, this will cause a division by zero or mathematically invalid configuration at runtime.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.


! B-field
q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp + C_mhd*erf((x_cc(i) - 0.5_wp)/sigma)

case (261) ! Blob
r0 = 1._wp/sqrt(8._wp)
r2 = x_cc(i)**2 + y_cc(j)**2
r = sqrt(r2)
alpha = r/r0
if (alpha < 1) then
q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp/sqrt(4._wp*pi)*(alpha**8 - 2._wp*alpha**4 + 1._wp)
! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/sqrt(4000._wp*pi) * (4096._wp*r2**4 - 128._wp*r2**2 + 1._wp)
! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp)
! q_prim_vf(E_idx)%sf(i,j,0) = 6._wp - q_prim_vf(B_idx%beg)%sf(i,j,0)**2/2._wp
end if

case (262) ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°)
! rotate by α = atan(2)
alpha = atan(2._wp)
cosA = cos(alpha)
sinA = sin(alpha)
! projection along shock normal
r = x_cc(i)*cosA + y_cc(j)*sinA

if (r <= 0.5_wp) then
! LEFT state: ρ=1, v∥=+10, v⊥=0, p=20, B∥=B⊥=5/√(4π)
q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
q_prim_vf(momxb)%sf(i, j, 0) = 10._wp*cosA
q_prim_vf(momxb + 1)%sf(i, j, 0) = 10._wp*sinA
q_prim_vf(E_idx)%sf(i, j, 0) = 20._wp
q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
- (5._wp/sqrt(4._wp*pi))*sinA
q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
+ (5._wp/sqrt(4._wp*pi))*cosA
else
! RIGHT state: ρ=1, v∥=−10, v⊥=0, p=1, B∥=B⊥=5/√(4π)
q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
q_prim_vf(momxb)%sf(i, j, 0) = -10._wp*cosA
q_prim_vf(momxb + 1)%sf(i, j, 0) = -10._wp*sinA
q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp
q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
- (5._wp/sqrt(4._wp*pi))*sinA
q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
+ (5._wp/sqrt(4._wp*pi))*cosA
end if
! v^z and B^z remain zero by default

case (270)
! This hardcoded case extrudes a 1D profile to initialize a 2D simulation domain
@: HardcodedReadValues()
Expand Down
3 changes: 3 additions & 0 deletions src/common/m_variables_conversion.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -882,6 +882,7 @@ contains

if (cont_damage) qK_prim_vf(damage_idx)%sf(j, k, l) = qK_cons_vf(damage_idx)%sf(j, k, l)

if (hyper_cleaning) qK_prim_vf(psi_idx)%sf(j, k, l) = qK_cons_vf(psi_idx)%sf(j, k, l)
#ifdef MFC_POST_PROCESS
if (bubbles_lagrange) qK_prim_vf(beta_idx)%sf(j, k, l) = qK_cons_vf(beta_idx)%sf(j, k, l)
#endif
Expand Down Expand Up @@ -1152,6 +1153,8 @@ contains

if (cont_damage) q_cons_vf(damage_idx)%sf(j, k, l) = q_prim_vf(damage_idx)%sf(j, k, l)

if (hyper_cleaning) q_cons_vf(psi_idx)%sf(j, k, l) = q_prim_vf(psi_idx)%sf(j, k, l)

end do
end do
end do
Expand Down
3 changes: 3 additions & 0 deletions src/post_process/m_data_output.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -363,6 +363,9 @@ contains
! Damage state variable
if (cont_damage) dbvars = dbvars + 1

! Hyperbolic cleaning for MHD
if (hyper_cleaning) dbvars = dbvars + 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: The hyperbolic cleaning scalar is counted as an output variable whenever hyper_cleaning is true, even if MHD is disabled; this can make the declared number of binary output variables (dbvars) larger than the number actually written, corrupting the binary file layout for non‑MHD runs where hyper_cleaning might be (mis)enabled. [logic error]

Severity Level: Major ⚠️
- ❌ Binary formatted database header mismatches written data.
- ❌ Binary post-processing output files become unreadable.
- ⚠️ Visualization tools reading Binary may mis-interpret data.
- ⚠️ Regression tests comparing Binary dumps may fail.
Suggested change
if (hyper_cleaning) dbvars = dbvars + 1
if (mhd .and. hyper_cleaning) dbvars = dbvars + 1
Steps of Reproduction ✅
1. Configure a run that requests Binary formatted post-processing output (format == 2) and
enable hyperbolic cleaning while MHD is disabled. The dbvars counter is computed inside
s_initialize_data_output_module in src/post_process/m_data_output.fpp; the hyper_cleaning
increment is at the lines shown above (lines 366-367).

2. Start the application so that s_initialize_data_output_module runs and sets dbvars (the
code path that computes dbvars includes the two lines at 366-367). Because hyper_cleaning
is true, dbvars is incremented unconditionally at 367 even though mhd is false.

3. The code path that writes the Binary file header in s_open_formatted_database_file
(src/post_process/m_data_output.fpp) writes the declared number of variables into the file
header: write (dbfile) m, n, p, dbvars. This write uses the dbvars value computed above.

4. Later the routines that actually write per-variable binary blocks
(s_write_variable_to_formatted_database_file and related code in the same file) skip
MHD-related variables when mhd is false, so fewer variable blocks are emitted than the
header advertised. The result is a mismatch between header (contains extra hyper_cleaning
variable) and actual data blocks; binary readers or downstream post-processing will
mis-interpret file layout, causing corrupted reads or failures.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/post_process/m_data_output.fpp
**Line:** 367:367
**Comment:**
	*Logic Error: The hyperbolic cleaning scalar is counted as an output variable whenever `hyper_cleaning` is true, even if MHD is disabled; this can make the declared number of binary output variables (`dbvars`) larger than the number actually written, corrupting the binary file layout for non‑MHD runs where `hyper_cleaning` might be (mis)enabled.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.


! Magnetic field
if (mhd) then
if (n == 0) then
Expand Down
10 changes: 10 additions & 0 deletions src/post_process/m_global_parameters.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ module m_global_parameters
integer :: b_size !< Number of components in the b tensor
integer :: tensor_size !< Number of components in the nonsymmetric tensor
logical :: cont_damage !< Continuum damage modeling
logical :: hyper_cleaning !< Hyperbolic cleaning for MHD
logical :: igr !< enable IGR
integer :: igr_order !< IGR reconstruction order
logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
Expand All @@ -143,6 +144,7 @@ module m_global_parameters
integer :: c_idx !< Index of color function
type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns.
integer :: damage_idx !< Index of damage state variable (D) for continuum damage model
integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD
!> @}

! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
Expand Down Expand Up @@ -407,6 +409,7 @@ contains
b_size = dflt_int
tensor_size = dflt_int
cont_damage = .false.
hyper_cleaning = .false.
igr = .false.

bc_x%beg = dflt_int; bc_x%end = dflt_int
Expand Down Expand Up @@ -812,6 +815,13 @@ contains
damage_idx = dflt_int
end if

if (hyper_cleaning) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: The hyperbolic cleaning variable is appended to the system state vector whenever hyper_cleaning is true, regardless of whether MHD is enabled or the problem is multidimensional, which can create a state layout that does not match the assumptions of the MHD hyperbolic cleaning solver (implemented only for 2D/3D MHD), leading to inconsistent indexing and potential downstream runtime errors. [logic error]

Severity Level: Major ⚠️
- ❌ State-vector layout corruption during accidental hyper_cleaning usage.
- ⚠️ MPI/IO allocation and indexing (MPI_IO_DATA) affected.
Suggested change
if (hyper_cleaning) then
if (hyper_cleaning .and. mhd .and. n > 0) then
Steps of Reproduction ✅
1. Configure a run where model_eqns selects the branch that executes the state-vector
annotation code in impure subroutine s_initialize_global_parameters_module (the large
routine that sets up cont_idx, mom_idx, B_idx, etc.). The psi_idx insertion block is at
lines 818..823 in that routine and will be executed when control reaches the sys_size
annotation logic inside the model_eqns handling.

2. Enable hyper_cleaning in the case input (so the module variable hyper_cleaning becomes
true by the time s_initialize_global_parameters_module runs). The module does not enforce
that hyper_cleaning implies mhd or multidimensional geometry anywhere in this file
(hyper_cleaning is an independent logical declared at line ~120).

3. With hyper_cleaning = .true. and mhd = .false. (or with a 1D run where n == 0), the
existing code at lines 818..823 will append psi_idx to sys_size even though the MHD state
(B_idx) or multidimensionality may not be present — changing the global sys_size and
therefore the expected state-vector layout.

4. The mismatch becomes observable when downstream code that assembles or solves the
system (or any IO/serialization expecting a specific ordering) assumes no psi state unless
MHD/multidimensionality is active. In practice this misconfiguration is realistic (user
sets hyper_cleaning without toggling mhd), so the current guardless insertion can produce
incorrect indexing. Guarding the insertion with ".and. mhd .and. n > 0" prevents adding
the psi field when MHD or multidimensional geometry is not active.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/post_process/m_global_parameters.fpp
**Line:** 818:818
**Comment:**
	*Logic Error: The hyperbolic cleaning variable is appended to the system state vector whenever `hyper_cleaning` is true, regardless of whether MHD is enabled or the problem is multidimensional, which can create a state layout that does not match the assumptions of the MHD hyperbolic cleaning solver (implemented only for 2D/3D MHD), leading to inconsistent indexing and potential downstream runtime errors.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

psi_idx = sys_size + 1
sys_size = psi_idx
else
psi_idx = dflt_int
end if

end if

if (chemistry) then
Expand Down
2 changes: 1 addition & 1 deletion src/post_process/m_mpi_proxy.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ contains
& 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', &
& 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', &
& 'output_partial_domain', 'relativity', 'cont_damage', 'bc_io', &
& 'down_sample','fft_wrt' ]
& 'down_sample','fft_wrt', 'hyper_cleaning' ]
call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
#:endfor

Expand Down
Loading
Loading