MHD Hyperbolic Divergence Cleaning #1086
Conversation
📝 WalkthroughWalkthroughThe PR replaces Powell divergence correction with hyperbolic (Dedner) divergence cleaning for MHD simulations. This involves removing the Powell module, introducing hyper_cleaning feature flags and parameters, extending the state vector with a psi component to track divergence, updating initialization/finalization paths, modifying Riemann solver wave speed constraints, and adjusting RHS computations to propagate cleaning via the psi variable. Changes
Sequence DiagramsequenceDiagram
participant TS as Time Stepper
participant RS as Riemann Solver
participant RHS as RHS Computation
participant State as State Vector
Note over TS,State: Hyper-Cleaning Flow
TS->>RS: Compute fluxes<br/>(constrain c_fast by ±hyper_cleaning_speed)
RS->>RS: Update flux_rs for psi<br/>(propagate B divergence)
RS->>RHS: Return constrained fluxes
RHS->>State: Compute psi RHS<br/>= -psi/hyper_cleaning_tau
RHS->>State: Update rhs_vf(psi_idx)<br/>with divergence & relaxation
State->>State: Advance q_cons_vf(psi_idx)<br/>via standard time integration
TS->>State: q_cons_vf(psi_idx) moderates<br/>∇·B via source term
Estimated code review effort🎯 4 (Complex) | ⏱️ ~45 minutes The changes span multiple critical paths (Riemann solver, RHS, state initialization) across 25+ files, introducing a new physical feature with interdependent logic (hyper_cleaning parameters, psi state variable, wave speed constraints, flux propagation). While the changes are cohesive and scoped to hyperbolic cleaning, they require understanding of MHD discretization, state vector extension, and numerical wave propagation. Possibly related PRs
Suggested labels
Suggested reviewers
Poem
🚥 Pre-merge checks | ✅ 2 | ❌ 1❌ Failed checks (1 warning)
✅ Passed checks (2 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing touches
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
|
beautiful! @ChrisZYJ Maybe something we can look at... |
|
Yeah, this the the hyperbolic solver I told you I was going to add. The one I put in this week is an outdated elliptic solver. I was hoping to update this one as I have time and redo our 2D convergence test that I added previously. |
|
Hey, @ChrisZYJ, I think I almost have this merged with master. I wanted to ask about this variable This appear to work fine if I just remove it, but I wanted to see if there is a purpose to it sticking around |
|
Hi @danieljvickers , thanks for the help with the PR!!
|
|
I don't mind finishing it up. It looks like it's just an issue casting double to complex. You've really done most the heavy lifting already. |
1 similar comment
|
I don't mind finishing it up. It looks like it's just an issue casting double to complex. You've really done most the heavy lifting already. |
|
Thanks! Please let me know if you need me to work on anything. |
|
MacOS failures are due to issues with updated gcc in the latest brew update. I implemented a fix on #1111, which currently is passing all tests. Once that gets merged, we can pull that in and resolve the failure. |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1086 +/- ##
==========================================
- Coverage 44.16% 44.14% -0.03%
==========================================
Files 71 70 -1
Lines 20417 20416 -1
Branches 1991 2002 +11
==========================================
- Hits 9018 9013 -5
+ Misses 10252 10244 -8
- Partials 1147 1159 +12 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
I am going to put this PR down for a couple of days, but I manged to do at least some tracking towards a solution. I tried to also damp in the RHS module by adding lines like this: When I have this term in AND the term in time stepper, I do actually get the rotor problem to converge for m=n=400 and when I drop the hyper_cleaning_tau = dt: |
|
Thanks so much for helping to debug!! The checkerboarding usually means something is completely wrong. Might be something stupid I changed after it first worked. I'll study your latest PR and get the bug fixed. |
|
I've narrowed down the bug:
I'm comparing the versions line by line and should have a fix soon. |
|
Should be fixed. The bug was that In this PR I removed |
|
I've also added docs, tests, and examples. Also deleted Powell's method as it's no longer relevant. If the checks for the new cases fail, increasing the tolerance should fix them (hyperbolic cleaning can be sensitive to machine errors). I don't think all the cases in |
|
Ahh darn, you are right. I even dug through the git diff yesterday and completely missed it. Thanks for taking time to look at it. Another set of eyes was helpful. Since this PR is good to go (I hope), I deleted the |
|
@danieljvickers No problem, glad we caught it! By the way, thanks for moving the damping equations to the RHS so they align better with the literature. I see the Frontier compilation is still failing, but don't want to cause merge conflicts with overlapping changes. Let me know if you want me to take a look! |
|
Yeah, o saw the failure yesterday and restarted it to see if that would fix it. I saw it pass before I added the GPU change, and I didn't touch anything in common or post process. So my bet is that this is spurious failures on frontier. But I'll take a look. |
|
@sbryngelson This PR is now ready for merge, as long as there isn't any other input from @ChrisZYJ |
|
very good, waiting on @anandrdbz for the AMDFlang compiler stuff first. |
PR Reviewer Guide 🔍Here are some key observations to aid the review process:
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
There was a problem hiding this comment.
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.)
| ! hard-coded psi | ||
| if (hyper_cleaning) then | ||
| do j = 0, m | ||
| 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) | ||
| end do | ||
| end do | ||
| end if |
There was a problem hiding this comment.
Suggestion: Relocate the hard-coded psi initialization to prevent it from being overwritten by a subsequent zero-initialization. [possible issue, importance: 9]
| ! hard-coded psi | |
| if (hyper_cleaning) then | |
| do j = 0, m | |
| 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) | |
| end do | |
| end do | |
| end if | |
| ! This block should be moved to 'src/pre_process/m_initial_condition.fpp' | |
| ! after the zero-initialization of psi, or that initialization should be removed/made conditional. | |
| if (hyper_cleaning) then | |
| do j = 0, m | |
| 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) | |
| end do | |
| end do | |
| end if |
| 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 |
There was a problem hiding this comment.
Suggestion: Add an else branch to the MHD Rotor problem to initialize the ambient region's density and momentum. [possible issue, importance: 8]
| 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 |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 5
🤖 Fix all issues with AI agents
In `@src/simulation/m_global_parameters.fpp`:
- Around line 853-855: When hyper_cleaning is true, add validation that
hyper_cleaning_speed and hyper_cleaning_tau are finite and > 0 to avoid
division-by-zero and invalid wave speeds; locate where hyper_cleaning,
hyper_cleaning_speed and hyper_cleaning_tau are initialized in
m_global_parameters.fpp (the routine that sets defaults/reads parameters) and
insert checks that use isfinite() (or the Fortran equivalent) and value > 0, and
if a check fails call the existing error/stop logging path with a clear message
naming hyper_cleaning_speed or hyper_cleaning_tau (so m_rhs.fpp line 993 and
m_riemann_solvers.fpp line 956 won’t receive invalid values), or set safe
fallback values and log a warning if that matches project behavior.
- Around line 294-295: The psi_idx integer (used by GPU kernels in
m_riemann_solvers and m_rhs) is missing from the device mirroring lists; update
the GPU_DECLARE and GPU_UPDATE macros/declarations to include psi_idx so the
variable is copied to/from the device and kept current for kernels. Locate the
GPU_DECLARE and GPU_UPDATE blocks where global parameters are listed (alongside
other indices) and add psi_idx to both the declaration and update lists to
ensure proper device mirroring for GPU kernels that reference psi_idx.
- Around line 1164-1168: The runtime flag hyper_cleaning unconditionally adds
psi_idx to sys_size even when MHD support (mhd) is not compiled in; add a
runtime guard that aborts when hyper_cleaning is true but mhd is false (e.g.,
check hyper_cleaning .and. .not. mhd and call s_mpi_abort with a descriptive
message) before altering psi_idx/sys_size so the extra equation cannot be
enabled without MHD; update the block that sets psi_idx/sys_size to first
validate mhd and call s_mpi_abort if invalid.
In `@src/simulation/m_rhs.fpp`:
- Around line 987-998: The divergence-cleaning damping term is using the stale
passed-in primitive array q_prim_vf instead of the up-to-date primitive computed
during c2p; in the hyper_cleaning block update the subtraction to read from
q_prim_qp%vf(psi_idx)%sf(j,k,l) divided by hyper_cleaning_tau so
rhs_vf(psi_idx)%sf is reduced using the current primitive psi value (keep the
surrounding loop and GPU pragmas and reference psi_idx, rhs_vf, q_prim_vf,
q_prim_qp, and hyper_cleaning_tau when making the change).
In `@src/simulation/m_riemann_solvers.fpp`:
- Around line 946-958: The normal-component magnetic flux
(flux_rs${XYZ}$_vf(..., B_idx%beg + norm_dir - 1)) is being computed before
applying the psi correction which reintroduces HLL dissipation; update the block
inside the hyper_cleaning branch to restore the original gating by applying (1 -
dir_flg(i + 1)) to the B_normal flux (or explicitly set flux_rs${XYZ}$_vf(...,
B_idx%beg + norm_dir - 1)=0._wp when dir_flg(i + 1) is true) before you add the
psi contribution, ensuring the psi update that uses qL_prim_rs${XYZ}$_vf,
qR_prim_rs${XYZ}$_vf and psi_idx only modifies the non-gated value.
🧹 Nitpick comments (1)
src/common/include/2dHardcodedIC.fpp (1)
256-266: Consider removing commented-out code.The blob implementation is correct. However, lines 263-265 contain commented-out alternative formulations that appear to be debug/exploration artifacts rather than documentation.
🧹 Suggested cleanup
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
| integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD | ||
| !> @} |
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
# Find all usages of psi_idx with context
rg -n "psi_idx" -B 2 -A 2 src/simulationRepository: MFlowCode/MFC
Length of output: 4223
🏁 Script executed:
# Check the variable declaration in m_global_parameters.fpp
head -n 300 src/simulation/m_global_parameters.fpp | tail -n 20Repository: MFlowCode/MFC
Length of output: 1770
🏁 Script executed:
# Search for GPU_DECLARE and GPU_UPDATE in the codebase to understand the pattern
rg -n "GPU_DECLARE|GPU_UPDATE" src/simulation | head -20Repository: MFlowCode/MFC
Length of output: 2026
🏁 Script executed:
# Check if psi_idx appears in any GPU_UPDATE declarations
rg -n "GPU_UPDATE.*psi_idx" src/simulationRepository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Get the full GPU_DECLARE and GPU_UPDATE sections in m_global_parameters.fpp
sed -n '290,320p' src/simulation/m_global_parameters.fppRepository: MFlowCode/MFC
Length of output: 1671
🏁 Script executed:
# Check if there are additional GPU_UPDATE declarations related to global parameters elsewhere
rg -n "GPU_UPDATE.*device=.*momxb" src/simulation -A 3Repository: MFlowCode/MFC
Length of output: 500
Add psi_idx to GPU_DECLARE and GPU_UPDATE declarations. The variable is used in GPU device code (e.g., m_riemann_solvers.fpp:954–956 to index array elements in flux calculations, and m_rhs.fpp:992–993 in GPU-enabled loops), but is missing from both GPU_DECLARE and GPU_UPDATE lists. Without device mirroring, GPU kernels access stale or uninitialized values.
GPU declaration fix
- $:GPU_DECLARE(create='[tensor_size,species_idx,c_idx]')
+ $:GPU_DECLARE(create='[tensor_size,species_idx,c_idx,psi_idx]')GPU update fix
$:GPU_UPDATE(device='[momxb,momxe,advxb,advxe,contxb,contxe, &
& bubxb,bubxe,intxb,intxe,sys_size,buff_size,E_idx, &
& alf_idx,n_idx,adv_n,adap_dt,pi_fac,strxb,strxe, &
- & chemxb,chemxe,c_idx,adap_dt_tol,adap_dt_max_iters]')
+ & chemxb,chemxe,c_idx,psi_idx,adap_dt_tol,adap_dt_max_iters]')🤖 Prompt for AI Agents
In `@src/simulation/m_global_parameters.fpp` around lines 294 - 295, The psi_idx
integer (used by GPU kernels in m_riemann_solvers and m_rhs) is missing from the
device mirroring lists; update the GPU_DECLARE and GPU_UPDATE
macros/declarations to include psi_idx so the variable is copied to/from the
device and kept current for kernels. Locate the GPU_DECLARE and GPU_UPDATE
blocks where global parameters are listed (alongside other indices) and add
psi_idx to both the declaration and update lists to ensure proper device
mirroring for GPU kernels that reference psi_idx.
| hyper_cleaning_speed = dflt_real | ||
| hyper_cleaning_tau = dflt_real | ||
|
|
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
rg -n "hyper_cleaning" -C 3 src --type f90 --type fppRepository: MFlowCode/MFC
Length of output: 84
🏁 Script executed:
rg -n "hyper_cleaning" -C 3 srcRepository: MFlowCode/MFC
Length of output: 21189
Add validation for hyper_cleaning_speed and hyper_cleaning_tau when hyper_cleaning = .true. These parameters default to dflt_real and are used in critical calculations without bounds checking. Specifically, hyper_cleaning_tau is used as a divisor in the RHS computation (m_rhs.fpp line 993) and hyper_cleaning_speed is squared in flux calculations (m_riemann_solvers.fpp line 956). Enforce that both are positive and finite to prevent division-by-zero and invalid wave speeds.
🤖 Prompt for AI Agents
In `@src/simulation/m_global_parameters.fpp` around lines 853 - 855, When
hyper_cleaning is true, add validation that hyper_cleaning_speed and
hyper_cleaning_tau are finite and > 0 to avoid division-by-zero and invalid wave
speeds; locate where hyper_cleaning, hyper_cleaning_speed and hyper_cleaning_tau
are initialized in m_global_parameters.fpp (the routine that sets defaults/reads
parameters) and insert checks that use isfinite() (or the Fortran equivalent)
and value > 0, and if a check fails call the existing error/stop logging path
with a clear message naming hyper_cleaning_speed or hyper_cleaning_tau (so
m_rhs.fpp line 993 and m_riemann_solvers.fpp line 956 won’t receive invalid
values), or set safe fallback values and log a warning if that matches project
behavior.
| if (hyper_cleaning) then | ||
| psi_idx = sys_size + 1 | ||
| sys_size = psi_idx | ||
| end if | ||
|
|
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
rg -n "hyper_cleaning" -C3 src/Repository: MFlowCode/MFC
Length of output: 21189
🏁 Script executed:
rg -n "hyper_cleaning.*=" -B5 -A5 src/pre_process/m_start_up.fpp | head -100Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Search for error checks or validations related to hyper_cleaning
rg -n "if.*hyper_cleaning.*mhd|if.*mhd.*hyper_cleaning" -B3 -A3 src/Repository: MFlowCode/MFC
Length of output: 745
🏁 Script executed:
# Search for where hyper_cleaning is read from input file or parsed
rg -n "hyper_cleaning" src/pre_process/m_start_up.fpp -B10 -A10 | grep -A10 -B10 "read\|input\|parse"Repository: MFlowCode/MFC
Length of output: 1420
🏁 Script executed:
# Check if mhd flag exists and how it's validated
rg -n "if.*mhd.*then" src/simulation/m_global_parameters.fpp | head -20Repository: MFlowCode/MFC
Length of output: 205
🏁 Script executed:
# Search for any validation/error checking in the initialization
rg -n "s_mpi_abort" src/simulation/m_global_parameters.fpp -B5 -A2Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Find where input parameters are read (likely a namelist)
rg -n "namelist\|read.*\(" src/pre_process/m_start_up.fpp | head -30Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Look for the actual input parameter reading section around line 155
sed -n '130,200p' src/pre_process/m_start_up.fppRepository: MFlowCode/MFC
Length of output: 3291
🏁 Script executed:
# Check for any validation that checks mhd and hyper_cleaning together
rg -n "mhd\|hyper_cleaning" src/simulation/m_global_parameters.fpp | grep -A2 -B2 "1164\|1165\|1166"Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Look at the section around 1164 in simulation/m_global_parameters.fpp
sed -n '1140,1180p' src/simulation/m_global_parameters.fppRepository: MFlowCode/MFC
Length of output: 1451
🏁 Script executed:
# Check if there's any validation checking mhd vs hyper_cleaning in the entire m_global_parameters.fpp
rg -n "mhd\|hyper_cleaning" src/simulation/m_global_parameters.fpp -C2 | grep -A5 -B5 "if\|validation\|error"Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Look for any error handling around the sys_size computation
sed -n '1100,1200p' src/simulation/m_global_parameters.fppRepository: MFlowCode/MFC
Length of output: 3997
🏁 Script executed:
# Check what model_eqns == 2 or 3 represents and if it's related to MHD
rg -n "model_eqns" src/simulation/m_global_parameters.fpp | head -40Repository: MFlowCode/MFC
Length of output: 512
🏁 Script executed:
# Search for model_eqns definitions/documentation
rg -n "model_eqns.*=" src/simulation/m_global_parameters.fpp -B3 -A3 | head -50Repository: MFlowCode/MFC
Length of output: 1820
🏁 Script executed:
# Look for any input validation that ties hyper_cleaning to mhd
rg -n "hyper_cleaning\|mhd" src/simulation/m_start_up.fpp -B2 -A2 | head -80Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Search for where MHD is enabled (mhd flag setting)
rg -n "mhd.*=" src/simulation/m_global_parameters.fpp -B3 -A3 | head -100Repository: MFlowCode/MFC
Length of output: 916
🏁 Script executed:
# Check if hyper_cleaning has any input validation
rg -n "hyper_cleaning" src/simulation/m_start_up.fpp -B5 -A5Repository: MFlowCode/MFC
Length of output: 741
🏁 Script executed:
# Check if there's any runtime validation that checks hyper_cleaning against mhd
rg -n "if.*hyper_cleaning" src/simulation/m_start_up.fpp -B10 -A10 | head -100Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Look for error checking or validation subroutines in simulation startup
rg -n "subroutine.*validation\|subroutine.*check" src/simulation/m_start_up.fpp -A30Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Check if hyper_cleaning is validated/constrained when reading input
sed -n '150,400p' src/simulation/m_start_up.fppRepository: MFlowCode/MFC
Length of output: 9639
🏁 Script executed:
# Look at the s_check_inputs subroutines to see if there's validation
rg -n "subroutine s_check_inputs" src/simulation/m_start_up.fpp -A50 | head -100Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Search all files for validation of hyper_cleaning with mhd
rg -n "hyper_cleaning.*mhd|mhd.*hyper_cleaning" src/Repository: MFlowCode/MFC
Length of output: 156
🏁 Script executed:
# Look at the entire check_inputs subroutine in simulation startup
sed -n '800,1200p' src/simulation/m_start_up.fpp | grep -A100 "subroutine s_check_inputs"Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Check what mhd really is (parameter vs variable)
sed -n '120,140p' src/simulation/m_global_parameters.fppRepository: MFlowCode/MFC
Length of output: 1922
🏁 Script executed:
# Check if mhd is compile-time or runtime
rg -n "logical.*mhd\|integer.*mhd" src/simulation/m_global_parameters.fpp | head -20Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Look for how hyper_cleaning relates to actual MHD physics
rg -n "hyper_cleaning" src/simulation/m_riemann_solvers.fpp -B5 -A5Repository: MFlowCode/MFC
Length of output: 2561
🏁 Script executed:
# Check if there's any CMake or build-time constraint
fd "CMakeLists.txt" | head -5Repository: MFlowCode/MFC
Length of output: 109
🏁 Script executed:
# Understand what happens when hyper_cleaning is true but mhd is not compiled
# Check if B_idx exists in non-MHD configurations
rg -n "B_idx" src/simulation/m_global_parameters.fpp | head -20Repository: MFlowCode/MFC
Length of output: 515
🏁 Script executed:
# Check what code blocks are guarded by mhd check
rg -n "if.*mhd.*then" src/simulation/m_global_parameters.fpp -A5 | head -50Repository: MFlowCode/MFC
Length of output: 771
Guard hyper_cleaning to require MHD at compile-time. The hyper_cleaning flag is a runtime parameter that unconditionally adds psi_idx to the system (lines 1164–1167) but is only meaningful when mhd is compiled in. While the code at line 649 in m_riemann_solvers.fpp has an inline comment linking hyper_cleaning to MHD, there is no runtime validation preventing hyper_cleaning=true when mhd is not compiled, which would introduce an unused equation. Add a runtime check (e.g., if (hyper_cleaning .and. .not. mhd) call s_mpi_abort(...)) to prevent this configuration.
🤖 Prompt for AI Agents
In `@src/simulation/m_global_parameters.fpp` around lines 1164 - 1168, The runtime
flag hyper_cleaning unconditionally adds psi_idx to sys_size even when MHD
support (mhd) is not compiled in; add a runtime guard that aborts when
hyper_cleaning is true but mhd is false (e.g., check hyper_cleaning .and. .not.
mhd and call s_mpi_abort with a descriptive message) before altering
psi_idx/sys_size so the extra equation cannot be enabled without MHD; update the
block that sets psi_idx/sys_size to first validate mhd and call s_mpi_abort if
invalid.
| if (hyper_cleaning) then | ||
| $:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3) | ||
| do l = 0, p | ||
| do k = 0, n | ||
| do j = 0, m | ||
| rhs_vf(psi_idx)%sf(j, k, l) = rhs_vf(psi_idx)%sf(j, k, l) - & | ||
| q_prim_vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau | ||
| end do | ||
| end do | ||
| end do | ||
| $:END_GPU_PARALLEL_LOOP() | ||
| end if |
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# Check for assignments between q_prim_vf and q_prim_qp in m_rhs.fpp
rg -n -B5 -A5 'q_prim_vf.*=.*q_prim_qp|q_prim_qp.*=.*q_prim_vf' src/simulation/m_rhs.fppRepository: MFlowCode/MFC
Length of output: 646
🏁 Script executed:
#!/bin/bash
# View the context around the hyper_cleaning block and check for earlier q_prim_vf updates
sed -n '950,1090p' src/simulation/m_rhs.fpp | cat -nRepository: MFlowCode/MFC
Length of output: 6817
🏁 Script executed:
#!/bin/bash
# Search for updates to q_prim_vf in the entire file to understand when it's synchronized
# Show line numbers and context
rg -n 'q_prim_vf.*=' src/simulation/m_rhs.fpp | head -30Repository: MFlowCode/MFC
Length of output: 148
🏁 Script executed:
#!/bin/bash
# Find the subroutine definition that contains the hyper_cleaning block
# Look backwards from line 987 to find the subroutine definition
sed -n '1,987p' src/simulation/m_rhs.fpp | grep -n 'subroutine\|function' | tail -5Repository: MFlowCode/MFC
Length of output: 449
🏁 Script executed:
#!/bin/bash
# View the subroutine signature and variable declarations
sed -n '644,750p' src/simulation/m_rhs.fpp | cat -nRepository: MFlowCode/MFC
Length of output: 5819
🏁 Script executed:
#!/bin/bash
# Search for psi_idx definition and usage
rg -n 'psi_idx' src/simulation/m_rhs.fpp | head -20Repository: MFlowCode/MFC
Length of output: 501
Use current primitive psi value in divergence cleaning source term.
The damping term at line 993 reads from q_prim_vf(psi_idx)%sf, which is the passed-in parameter not synchronized with the current state. Since q_prim_qp%vf(psi_idx) is computed at line 71 during the conservative-to-primitive conversion and is already GPU-enabled, use it instead to ensure the damping applies the current divergence cleaning potential.
Change line 993 from:
q_prim_vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau
to:
q_prim_qp%vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau
🤖 Prompt for AI Agents
In `@src/simulation/m_rhs.fpp` around lines 987 - 998, The divergence-cleaning
damping term is using the stale passed-in primitive array q_prim_vf instead of
the up-to-date primitive computed during c2p; in the hyper_cleaning block update
the subtraction to read from q_prim_qp%vf(psi_idx)%sf(j,k,l) divided by
hyper_cleaning_tau so rhs_vf(psi_idx)%sf is reduced using the current primitive
psi value (keep the surrounding loop and GPU pragmas and reference psi_idx,
rhs_vf, q_prim_vf, q_prim_qp, and hyper_cleaning_tau when making the change).
| flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + i) = (s_M*(vel_R(dir_idx(1))*B%R(i + 1) - vel_R(i + 1)*B%R(norm_dir)) - & | ||
| s_P*(vel_L(dir_idx(1))*B%L(i + 1) - vel_L(i + 1)*B%L(norm_dir)) + & | ||
| s_M*s_P*(B%L(i + 1) - B%R(i + 1)))/(s_M - s_P) | ||
| end do | ||
|
|
||
| if (hyper_cleaning) then | ||
| ! propagate magnetic field divergence as a wave | ||
| flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) = flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) + & | ||
| (s_M*qR_prim_rs${XYZ}$_vf(j + 1, k, l, psi_idx) - s_P*qL_prim_rs${XYZ}$_vf(j, k, l, psi_idx))/(s_M - s_P) | ||
|
|
||
| flux_rs${XYZ}$_vf(j, k, l, psi_idx) = (hyper_cleaning_speed**2*(s_M*B%R(norm_dir) - s_P*B%L(norm_dir)) + s_M*s_P*(qL_prim_rs${XYZ}$_vf(j, k, l, psi_idx) - qR_prim_rs${XYZ}$_vf(j + 1, k, l, psi_idx)))/(s_M - s_P) | ||
| else | ||
| flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) = 0._wp ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero |
There was a problem hiding this comment.
Restore normal-B flux gating before adding psi
Line 946-948 now compute the normal-component induction flux before the psi correction, which reintroduces a non-zero HLL dissipation term in B_normal (this matches the rotor failure noted in the PR discussion). Reapply the (1 - dir_flg(i + 1)) gating (or explicitly zero B_normal) before adding the psi term.
🔧 Suggested fix
- flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + i) = (s_M*(vel_R(dir_idx(1))*B%R(i + 1) - vel_R(i + 1)*B%R(norm_dir)) - &
- s_P*(vel_L(dir_idx(1))*B%L(i + 1) - vel_L(i + 1)*B%L(norm_dir)) + &
- s_M*s_P*(B%L(i + 1) - B%R(i + 1)))/(s_M - s_P)
+ flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + i) = (1._wp - dir_flg(i + 1))*( &
+ s_M*(vel_R(dir_idx(1))*B%R(i + 1) - vel_R(i + 1)*B%R(norm_dir)) - &
+ s_P*(vel_L(dir_idx(1))*B%L(i + 1) - vel_L(i + 1)*B%L(norm_dir)) + &
+ s_M*s_P*(B%L(i + 1) - B%R(i + 1)))/(s_M - s_P)🤖 Prompt for AI Agents
In `@src/simulation/m_riemann_solvers.fpp` around lines 946 - 958, The
normal-component magnetic flux (flux_rs${XYZ}$_vf(..., B_idx%beg + norm_dir -
1)) is being computed before applying the psi correction which reintroduces HLL
dissipation; update the block inside the hyper_cleaning branch to restore the
original gating by applying (1 - dir_flg(i + 1)) to the B_normal flux (or
explicitly set flux_rs${XYZ}$_vf(..., B_idx%beg + norm_dir - 1)=0._wp when
dir_flg(i + 1) is true) before you add the psi contribution, ensuring the psi
update that uses qL_prim_rs${XYZ}$_vf, qR_prim_rs${XYZ}$_vf and psi_idx only
modifies the non-gated value.
There was a problem hiding this comment.
Pull request overview
Implements GLM (hyperbolic) divergence cleaning for 2D/3D MHD by introducing a new psi state variable and associated runtime parameters, replacing the prior Powell-based approach. Updates the toolchain schema/validation, simulation/pre/post-processing plumbing, documentation, and example/test cases to support the new capability.
Changes:
- Introduces
hyper_cleaning+hyper_cleaning_speed/hyper_cleaning_tauand threads a newpsi_idxstate variable through pre/sim/post pipelines. - Updates the MHD HLL flux path and adds
psidamping behavior in the RHS; removes the Powell module. - Adds/updates example cases and test harness references; updates docs to cite Dedner et al. (2002).
Reviewed changes
Copilot reviewed 28 out of 31 changed files in this pull request and generated 5 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/test/cases.py | Updates MHD test entry to use the hyper-cleaning Orszag–Tang example. |
| toolchain/mfc/run/case_dicts.py | Adds hyper_cleaning/speed/tau to the toolchain parameter schema; removes Powell key. |
| toolchain/mfc/case_validator.py | Switches validation from powell to hyper_cleaning constraints. |
| tests/F057F8E6/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/EC8F87D9/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/B4DC99F9/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| src/simulation/m_time_steppers.fpp | Allocates/deallocates psi primitive storage when hyper_cleaning is enabled. |
| src/simulation/m_start_up.fpp | Removes Powell module init/finalize wiring; reads hyper-cleaning params. |
| src/simulation/m_riemann_solvers.fpp | Adds hyper-cleaning wave-speed adjustment and GLM flux terms in HLL path. |
| src/simulation/m_rhs.fpp | Adds psi handling in RHS setup and a decay source term. |
| src/simulation/m_mpi_proxy.fpp | Broadcasts hyper_cleaning and its parameters via MPI. |
| src/simulation/m_mhd.fpp | Removes Powell source-term module implementation. |
| src/simulation/m_global_parameters.fpp | Adds hyper_cleaning, psi_idx, and GLM parameter definitions / sys_size accounting. |
| src/pre_process/m_start_up.fpp | Reads hyper_cleaning and introduces a psi initialization hook. |
| src/pre_process/m_mpi_proxy.fpp | Broadcasts hyper_cleaning via MPI in pre_process. |
| src/pre_process/m_initial_condition.fpp | Initializes psi to zero when enabled. |
| src/pre_process/m_global_parameters.fpp | Adds hyper_cleaning and psi_idx sizing in pre_process globals. |
| src/post_process/m_start_up.fpp | Writes psi to formatted database output when enabled. |
| src/post_process/m_mpi_proxy.fpp | Broadcasts hyper_cleaning via MPI in post_process. |
| src/post_process/m_global_parameters.fpp | Adds hyper_cleaning and psi_idx sizing in post_process globals. |
| src/post_process/m_data_output.fpp | Increments DB variable count when psi is included. |
| src/common/m_variables_conversion.fpp | Propagates psi through primitive↔conservative conversion. |
| src/common/include/2dHardcodedIC.fpp | Adds new hardcoded 2D MHD ICs (rotor, Gaussian divergence pulse, tilted shock tube) and reassigns the vortex case id. |
| examples/2D_orszag_tang_hyper_cleaning/case.py | Adds Orszag–Tang example configured for hyperbolic cleaning. |
| examples/2D_mhd_rotor/case.py | Adds a 2D MHD rotor example configured for hyperbolic cleaning. |
| examples/2D_mhd_magnetic_vortex/case.py | Updates hardcoded IC id to match the vortex case renumbering. |
| docs/documentation/references.md | Replaces Powell reference with Dedner et al. (2002). |
| docs/documentation/case.md | Documents hyper_cleaning and its parameters in the case-input reference. |
Comments suppressed due to low confidence (1)
src/simulation/m_global_parameters.fpp:298
psi_idxis used in GPU kernels (e.g., Riemann solver / RHS), but it is not included in the module’sGPU_DECLARElists (and therefore won’t be present underdefault(present)), nor is it included in theGPU_UPDATEindex updates. Addpsi_idxto the appropriateGPU_DECLARE/GPU_UPDATElists, and consider setting it todflt_intwhenhyper_cleaningis false to avoid an uninitialized index.
integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD
!> @}
$:GPU_DECLARE(create='[sys_size,E_idx,n_idx,bub_idx,alf_idx,gamma_idx]')
$:GPU_DECLARE(create='[pi_inf_idx,B_idx,stress_idx,xi_idx,b_size]')
$:GPU_DECLARE(create='[tensor_size,species_idx,c_idx]')
| self.prohibit(powell and fd_order is None, | ||
| "fd_order must be set if Powell's method is enabled") | ||
| self.prohibit(hyper_cleaning and not mhd, | ||
| "Hyperbolic cleaning requires mhd to be enabled") |
There was a problem hiding this comment.
check_mhd_simulation currently allows hyper_cleaning = T with riemann_solver = 4 (HLLD), but the simulation-side hyperbolic cleaning fluxes appear to be implemented only in the HLL path. This can silently run with incorrect (missing) GLM fluxes. Add a constraint to prohibit hyper_cleaning when riemann_solver == 4 (or implement the corresponding HLLD+GLM fluxes).
| "Hyperbolic cleaning requires mhd to be enabled") | |
| "Hyperbolic cleaning requires mhd to be enabled") | |
| self.prohibit(hyper_cleaning and riemann_solver == 4, | |
| "Hyperbolic cleaning is not supported with HLLD (riemann_solver = 4)") |
| self.prohibit(hyper_cleaning and not mhd, | ||
| "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") |
There was a problem hiding this comment.
When hyper_cleaning is enabled, hyper_cleaning_speed and hyper_cleaning_tau are required for stability/correctness (and hyper_cleaning_tau is used as a divisor in the RHS). Add validation to require both parameters to be provided and strictly positive.
| ! hard-coded psi | ||
| if (hyper_cleaning) then | ||
| do j = 0, m | ||
| 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) | ||
| end do | ||
| end do | ||
| end if | ||
|
|
There was a problem hiding this comment.
This hard-coded psi initialization runs for every case with hyper_cleaning = T, overriding whatever the initial condition module produced (including the earlier "psi initialized to zero" logic). This makes psi nonzero by default for unrelated MHD cases. Consider removing this, or gating it behind a specific hard-coded IC/case id/patch parameter so it only applies to the intended validation case(s).
| ! hard-coded psi | |
| if (hyper_cleaning) then | |
| do j = 0, m | |
| 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) | |
| end do | |
| end do | |
| end if |
| 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) |
There was a problem hiding this comment.
The Gaussian psi initialization uses default-kind literals (1d-2, 2.0, 0.05) and only writes l=0, which is inconsistent with the codebase’s real(wp) usage and will not correctly initialize 3D cases (p>0). Use _wp literals and loop over l=0,p (or assign the full 3D array) if this initialization is kept.
| if (hyper_cleaning) then | ||
| $:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3) | ||
| do l = 0, p | ||
| do k = 0, n | ||
| do j = 0, m |
There was a problem hiding this comment.
The psi damping source term is applied inside the dimensional-splitting loop, so in 2D/3D it will be applied num_dims times per RHS evaluation (over-damping by a factor of 2 or 3). Move this source term outside the directional loop, or split it consistently (e.g., divide by num_dims) if that’s the intended operator splitting.
![cubic-dev-ai[bot]](https://avatars.githubusercontent.com/in/1082092?s=60&v=4){kind=small}
There was a problem hiding this comment.
7 issues found across 31 files
Prompt for AI agents (all issues)
Check if these issues are valid — if so, understand the root cause of each and fix them.
<file name="src/pre_process/m_start_up.fpp">
<violation number="1" location="src/pre_process/m_start_up.fpp:851">
P2: Hard-coded Gaussian psi initialization runs for every hyper_cleaning case after IC setup, overwriting any user/loaded psi values and only setting the z=0 plane. This makes non-Gaussian or 3D initial conditions incorrect. Gate this behind a specific case/parameter or move it to a case-specific IC routine.</violation>
</file>
<file name="src/pre_process/m_global_parameters.fpp">
<violation number="1" location="src/pre_process/m_global_parameters.fpp:892">
P2: `hyper_cleaning` can be enabled when `model_eqns` is not 2/3, leaving `psi_idx` uninitialized while other modules still use it under `hyper_cleaning` guards. Add validation or disable hyper_cleaning for non-MHD configurations to avoid invalid indexing.</violation>
</file>
<file name="tests/F057F8E6/golden-metadata.txt">
<violation number="1" location="tests/F057F8E6/golden-metadata.txt:29">
P2: Redact machine-specific absolute paths in committed golden metadata to avoid leaking developer environment details and to keep snapshots portable.</violation>
</file>
<file name="tests/EC8F87D9/golden-metadata.txt">
<violation number="1" location="tests/EC8F87D9/golden-metadata.txt:29">
P3: Avoid committing user-specific absolute paths in golden metadata; redact or replace with a portable placeholder so the test artifact doesn’t leak local environment details.</violation>
</file>
<file name="src/simulation/m_riemann_solvers.fpp">
<violation number="1" location="src/simulation/m_riemann_solvers.fpp:650">
P1: `c_fast%L` is forced negative via `min(..., -hyper_cleaning_speed)` even though `c_fast` is a positive speed, which can invert the left wave speed and destabilize the solver. Use a positive magnitude expansion instead.</violation>
</file>
<file name="src/simulation/m_global_parameters.fpp">
<violation number="1" location="src/simulation/m_global_parameters.fpp:294">
P2: psi_idx is introduced for hyperbolic cleaning but never added to GPU_DECLARE/GPU_UPDATE lists, even though it is used inside GPU_PARALLEL_LOOP kernels. This can leave psi_idx undefined on the device and cause incorrect indexing or GPU runtime failures when hyper_cleaning is enabled.</violation>
</file>
<file name="tests/B4DC99F9/golden-metadata.txt">
<violation number="1" location="tests/B4DC99F9/golden-metadata.txt:29">
P2: Golden metadata embeds a developer-specific absolute path (username and local directory). This leaks local environment details and makes the test artifact machine-specific. Consider redacting or replacing with a portable placeholder before committing.</violation>
</file>
Since this is your first cubic review, here's how it works:
- cubic automatically reviews your code and comments on bugs and improvements
- Teach cubic by replying to its comments. cubic learns from your replies and gets better over time
- Ask questions if you need clarification on any suggestion
Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.
| end if | ||
|
|
||
| if (hyper_cleaning) then ! mhd | ||
| c_fast%L = min(c_fast%L, -hyper_cleaning_speed) |
There was a problem hiding this comment.
P1: c_fast%L is forced negative via min(..., -hyper_cleaning_speed) even though c_fast is a positive speed, which can invert the left wave speed and destabilize the solver. Use a positive magnitude expansion instead.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_riemann_solvers.fpp, line 650:
<comment>`c_fast%L` is forced negative via `min(..., -hyper_cleaning_speed)` even though `c_fast` is a positive speed, which can invert the left wave speed and destabilize the solver. Use a positive magnitude expansion instead.</comment>
<file context>
@@ -646,6 +646,11 @@ contains
end if
+ if (hyper_cleaning) then ! mhd
+ c_fast%L = min(c_fast%L, -hyper_cleaning_speed)
+ c_fast%R = max(c_fast%R, hyper_cleaning_speed)
+ end if
</file context>
| if (hyper_cleaning) then | ||
| do j = 0, m | ||
| 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)) |
There was a problem hiding this comment.
P2: Hard-coded Gaussian psi initialization runs for every hyper_cleaning case after IC setup, overwriting any user/loaded psi values and only setting the z=0 plane. This makes non-Gaussian or 3D initial conditions incorrect. Gate this behind a specific case/parameter or move it to a case-specific IC routine.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/pre_process/m_start_up.fpp, line 851:
<comment>Hard-coded Gaussian psi initialization runs for every hyper_cleaning case after IC setup, overwriting any user/loaded psi values and only setting the z=0 plane. This makes non-Gaussian or 3D initial conditions incorrect. Gate this behind a specific case/parameter or move it to a case-specific IC routine.</comment>
<file context>
@@ -842,6 +844,16 @@ contains
+ if (hyper_cleaning) then
+ do j = 0, m
+ 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)
+ end do
</file context>
| sys_size = damage_idx | ||
| end if | ||
|
|
||
| if (hyper_cleaning) then |
There was a problem hiding this comment.
P2: hyper_cleaning can be enabled when model_eqns is not 2/3, leaving psi_idx uninitialized while other modules still use it under hyper_cleaning guards. Add validation or disable hyper_cleaning for non-MHD configurations to avoid invalid indexing.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/pre_process/m_global_parameters.fpp, line 892:
<comment>`hyper_cleaning` can be enabled when `model_eqns` is not 2/3, leaving `psi_idx` uninitialized while other modules still use it under `hyper_cleaning` guards. Add validation or disable hyper_cleaning for non-MHD configurations to avoid invalid indexing.</comment>
<file context>
@@ -886,6 +889,11 @@ contains
sys_size = damage_idx
end if
+ if (hyper_cleaning) then
+ psi_idx = sys_size + 1
+ sys_size = psi_idx
</file context>
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /home/chris/source/MFC_mhd_hypercleaning/build/venv/bin/fypp |
There was a problem hiding this comment.
P2: Redact machine-specific absolute paths in committed golden metadata to avoid leaking developer environment details and to keep snapshots portable.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At tests/F057F8E6/golden-metadata.txt, line 29:
<comment>Redact machine-specific absolute paths in committed golden metadata to avoid leaking developer environment details and to keep snapshots portable.</comment>
<file context>
@@ -0,0 +1,185 @@
+ OpenACC : OFF
+ OpenMP : OFF
+
+ Fypp : /home/chris/source/MFC_mhd_hypercleaning/build/venv/bin/fypp
+ Doxygen :
+
</file context>
| type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. | ||
| integer :: c_idx !< Index of color function | ||
| integer :: damage_idx !< Index of damage state variable (D) for continuum damage model | ||
| integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD |
There was a problem hiding this comment.
P2: psi_idx is introduced for hyperbolic cleaning but never added to GPU_DECLARE/GPU_UPDATE lists, even though it is used inside GPU_PARALLEL_LOOP kernels. This can leave psi_idx undefined on the device and cause incorrect indexing or GPU runtime failures when hyper_cleaning is enabled.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_global_parameters.fpp, line 294:
<comment>psi_idx is introduced for hyperbolic cleaning but never added to GPU_DECLARE/GPU_UPDATE lists, even though it is used inside GPU_PARALLEL_LOOP kernels. This can leave psi_idx undefined on the device and cause incorrect indexing or GPU runtime failures when hyper_cleaning is enabled.</comment>
<file context>
@@ -290,6 +291,7 @@ module m_global_parameters
type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns.
integer :: c_idx !< Index of color function
integer :: damage_idx !< Index of damage state variable (D) for continuum damage model
+ integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD
!> @}
$:GPU_DECLARE(create='[sys_size,E_idx,n_idx,bub_idx,alf_idx,gamma_idx]')
</file context>
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /home/chris/source/MFC_mhd_hypercleaning/build/venv/bin/fypp |
There was a problem hiding this comment.
P2: Golden metadata embeds a developer-specific absolute path (username and local directory). This leaks local environment details and makes the test artifact machine-specific. Consider redacting or replacing with a portable placeholder before committing.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At tests/B4DC99F9/golden-metadata.txt, line 29:
<comment>Golden metadata embeds a developer-specific absolute path (username and local directory). This leaks local environment details and makes the test artifact machine-specific. Consider redacting or replacing with a portable placeholder before committing.</comment>
<file context>
@@ -0,0 +1,185 @@
+ OpenACC : OFF
+ OpenMP : OFF
+
+ Fypp : /home/chris/source/MFC_mhd_hypercleaning/build/venv/bin/fypp
+ Doxygen :
+
</file context>
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp |
There was a problem hiding this comment.
P3: Avoid committing user-specific absolute paths in golden metadata; redact or replace with a portable placeholder so the test artifact doesn’t leak local environment details.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At tests/EC8F87D9/golden-metadata.txt, line 29:
<comment>Avoid committing user-specific absolute paths in golden metadata; redact or replace with a portable placeholder so the test artifact doesn’t leak local environment details.</comment>
<file context>
@@ -0,0 +1,193 @@
+ OpenACC : OFF
+ OpenMP : OFF
+
+ Fypp : /home/dan/Documents/repos/MFC/build/venv/bin/fypp
+ Doxygen :
+
</file context>
3 participants
User description
Description
Implements the correct hyperbolic divergence cleaning for 2D/3D MHD in place of the current Powell's method. It is fully working and validated (with results below).
It is based on an outdated MFC branch, so it needs to be rebased, cleaned up, and made GPU-ready. @danieljvickers might be interested in collaborating on this; otherwise, I will get this ready as time permits.
Formulation
Type of change
Scope
How Has This Been Tested?
All the test cases are temporarily placed in the
case_newfolder.Orszag-Tang result at final time; comparison across WENO variants. Top row is p. Bottom row is ∇·B. The last column shows the last saved time step before crashing, for a simulation without hyperbolic cleaning. 3Z is a typo; it should be 5Z.
Checklist
docs/)examples/that demonstrate my new feature performing as expected.They run to completion and demonstrate "interesting physics"
./mfc.sh formatbefore committing my codeIf your code changes any code source files (anything in
src/simulation)To make sure the code is performing as expected on GPU devices, I have:
nvtxranges so that they can be identified in profiles./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.PR Type
Enhancement, Bug fix
Description
Replace Powell's method with hyperbolic (GLM) divergence cleaning for MHD
Add new state variable
psi_idxfor tracking magnetic divergenceImplement hyperbolic cleaning in Riemann solver flux calculations
Add three new MHD test cases: rotor, Gaussian pulse, tilted shock tube
Diagram Walkthrough
File Walkthrough
19 files
Add three new MHD hardcoded initial conditionsAdd hyperbolic cleaning parameters and psi indexImplement hyperbolic cleaning in flux calculationsAdd hyperbolic cleaning to MPI broadcast variablesAdd hyperbolic cleaning support to post-processorInitialize psi state variable for hyperbolic cleaningAdd hyperbolic cleaning parameters to pre-processorReplace Powell RHS with hyperbolic cleaning damping termAllocate and deallocate psi state variable arraysAdd hyperbolic cleaning to post-processor initializationAdd psi conversion between primitive and conservativeBroadcast hyperbolic cleaning flag in pre-processorBroadcast hyperbolic cleaning flag in post-processorInitialize psi to zero in initial conditionsAdd psi variable to database output countUpdate Orszag-Tang case to use hyperbolic cleaningUpdate validation to check hyperbolic cleaning constraintsAdd hyperbolic cleaning parameters to case dictionaryUpdate magnetic vortex case to use new case ID1 files
Remove Powell's method, add hyperbolic cleaning imports2 files
New 2D MHD rotor test case with hyperbolic cleaningReplace Powell test with hyperbolic cleaning test2 files
Document hyperbolic cleaning parameters and remove PowellReplace Powell reference with Dedner et al. reference7 files
Summary by CodeRabbit
Release Notes
New Features
Deprecation
Documentation
✏️ Tip: You can customize this high-level summary in your review settings.