MHD Hyperbolic Divergence Cleaning #1086
Conversation
📝 WalkthroughWalkthroughThis PR introduces hyperbolic cleaning for MHD simulations with a new psi state variable to control magnetic divergence. It adds eight new test case configurations for 2D Orszag-Tang and rotor simulations, renames Powell MHD modules to Clean, and integrates the feature across preprocessing, simulation, and post-processing modules via parameter additions and state variable handling. Changes
Estimated code review effort🎯 4 (Complex) | ⏱️ ~45 minutes 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.08% -0.09%
==========================================
Files 71 71
Lines 20417 20473 +56
Branches 1991 2006 +15
==========================================
+ Hits 9018 9026 +8
- Misses 10252 10287 +35
- Partials 1147 1160 +13 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
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}
| ! 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: Remove the hard-coded, test-specific initialization for the hyperbolic cleaning variable psi. [possible issue, importance: 7]
| ! 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 | |
| ! 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 |
| if (hyper_cleaning) then | ||
| do l = 0, p | ||
| do k = 0, n | ||
| do j = 0, m | ||
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & | ||
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & | ||
| exp(-dt/hyper_cleaning_tau) | ||
| end do | ||
| end do | ||
| end do | ||
| end if |
There was a problem hiding this comment.
Suggestion: Move the hyperbolic cleaning damping term into the RHS flux calculation to ensure it is applied at each stage of the Runge-Kutta time integration, rather than only at the end of the full step. [possible issue, importance: 9]
| if (hyper_cleaning) then | |
| do l = 0, p | |
| do k = 0, n | |
| do j = 0, m | |
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & | |
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & | |
| exp(-dt/hyper_cleaning_tau) | |
| end do | |
| end do | |
| end do | |
| end if | |
| if (hyper_cleaning) then | |
| do l = 0, p | |
| do k = 0, n | |
| do j = 0, m | |
| flux_vf(psi_idx)%sf(j,k,l) = flux_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 if |
| if (hyperelasticity) call s_initialize_hyperelastic_module() | ||
|
|
||
| if (mhd .and. powell) call s_initialize_mhd_powell_module | ||
| if (powell) call s_initialize_mhd_clean_module |
There was a problem hiding this comment.
Suggestion: Correct the condition for initializing the divergence-cleaning module to use the hyper_cleaning flag instead of the powell flag. [general, importance: 8]
| if (powell) call s_initialize_mhd_clean_module | |
| if (hyper_cleaning) call s_initialize_mhd_clean_module |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 7
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (4)
src/simulation/m_mhd.fpp (1)
41-45: Double allocation of arrays will cause runtime error in 3D.Line 41 allocates
dw_dx_mhd,dw_dy_mhd,dw_dz_mhdunconditionally, but line 44 allocates them again whenp > 0. This will cause a double allocation error in 3D simulations.🐛 Proposed fix
- @:ALLOCATE(du_dx_mhd(0:m,0:n,0:p), dv_dx_mhd(0:m,0:n,0:p), dw_dx_mhd(0:m,0:n,0:p)) + @:ALLOCATE(du_dx_mhd(0:m,0:n,0:p), dv_dx_mhd(0:m,0:n,0:p)) @:ALLOCATE(du_dy_mhd(0:m,0:n,0:p), dv_dy_mhd(0:m,0:n,0:p), dw_dy_mhd(0:m,0:n,0:p)) if (p > 0) then - @:ALLOCATE(dw_dx_mhd(0:m,0:n,0:p), dw_dy_mhd(0:m,0:n,0:p), dw_dz_mhd(0:m,0:n,0:p)) + @:ALLOCATE(dw_dx_mhd(0:m,0:n,0:p), dw_dz_mhd(0:m,0:n,0:p)) + @:ALLOCATE(du_dz_mhd(0:m,0:n,0:p), dv_dz_mhd(0:m,0:n,0:p)) end ifNote: Looking at the declared module variables (lines 24-26), you may need to review which derivatives are actually needed for 2D vs 3D and adjust accordingly. The
du_dz_mhd,dv_dz_mhd, anddw_dz_mhdarrays are declared butdu_dz_mhdanddv_dz_mhddon't appear to be allocated anywhere.src/pre_process/m_start_up.fpp (1)
829-855: Hard-coded psi init overrides IC and only initializes l=0.
This clobbers any preexisting/loaded IC and leaves other z-planes untouched for 3D. It also mixes literal kinds. Consider moving this to a case/IC patch or at least initialize all planes with_wpliterals.🐛 Proposed fix (initialize all planes; use wp literals)
- integer :: j, k + integer :: j, k, l ... - ! 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 + ! hard-coded psi (consider moving to case/IC patch) + if (hyper_cleaning) then + 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._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 + end do + end ifsrc/simulation/m_riemann_solvers.fpp (1)
4971-4984: Uninitialized components in Powell flux_gsrc copy.
Whenpowell=.true., only component 1 offlux_gsrc_rs*_vfis assigned (line 962), but lines 4971–4984 copy allsys_sizecomponents intoflux_gsrc_vf. Components 2 throughsys_sizewill contain uninitialized memory. Either zero all components offlux_gsrc_rs*_vfafter allocation or initialize only the required components before copying.src/simulation/m_global_parameters.fpp (1)
294-298: Addpsi_idxanddamage_idxto the GPU global parameter declarations.Both
psi_idxanddamage_idxare state variable indices used in GPU device contexts (GPU_ENTER_DATA in m_rhs.fpp, @:ACC_SETUP_SFs in m_time_steppers.fpp, and Riemann solver kernels). They follow the same pattern as other state variable indices (E_idx, B_idx, stress_idx, etc.) that are explicitly included in theGPU_DECLAREstatements. Add them to the appropriate GPU_DECLARE line for consistency and to ensure availability in GPU device code.
🤖 Fix all issues with AI agents
In `@case_new/gauss_div_pulse_2/case.py`:
- Around line 12-18: The comment claiming a 1D domain ("n=1" and using a
non-ASCII ×) is inconsistent with the actual grid settings where "m": 100 and
"n": 100; either update the comment to accurately describe a 2D computational
domain (e.g., replace the non-ASCII × with 'x' and state m=100, n=100 over
x_domain%beg/end and y_domain%beg/end) or change the grid parameter "n" to 1 if
the intent was a 1D x-only run; locate the grid dict entries ("x_domain%beg",
"x_domain%end", "y_domain%beg", "y_domain%end", "m", "n") and make the
corresponding comment/text consistent with the chosen configuration.
- Around line 56-70: Add the missing hard-coded initial condition ID and
parameters so the Gaussian divergence pulse is actually activated: set hcid to
260 and add the two missing patch_icpp parameters patch_icpp(1)%a(2) and
patch_icpp(1)%a(3) (use the same numeric values used in your validation run),
alongside the existing patch_icpp(1)%a(1) entry so the IC block matches the
hardcoded expectation.
In `@case_new/slanted_shock_tube/case.py`:
- Around line 5-24: The file and meta should be consistent: decide whether this
is a 1D Gaussian divergence pulse or a 2D slanted shock tube and update the
header/comment, hcid, and domain/grid accordingly. If you mean a 1D y-slice
Gaussian pulse, change the header from "Gaussian Divergence Pulse (1D via
uniform y-slice)" to match the filename or vice versa, set "y_domain%end" to a
very small extent (e.g. equal to a single cell), reduce "n" to 1 (or minimal),
keep "m" as the long x resolution, and ensure hcid reflects the Gaussian pulse
case; alternatively, if this is the slanted shock tube test, update the
header/comment and hcid to describe a 2D slanted shock tube and keep the current
2D "y_domain%end" and "n"/"m" values. Locate and edit the JSON keys
"y_domain%end", "n", "m", the header comment block, and the hcid identifier to
make them consistent.
In `@src/common/include/2dHardcodedIC.fpp`:
- Around line 242-254: In the Gaussian divergence pulse case (case 260) validate
patch_icpp(patch_id)%a(3) (sigma) before using it: check for zero or near-zero
and either set a safe minimum (e.g., sigma_min) or raise/log an error and abort
per .fpp guidelines; then compute C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp and use
that safe sigma for the erf denominator and any other divisions so
q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1 + C_mhd*erf((x_cc(i)-0.5_wp)/sigma) cannot
produce NaNs. Ensure the check is colocated with the lines that set eps_mhd,
sigma and C_mhd in the Gaussian Divergence Pulse block.
In `@src/simulation/m_rhs.fpp`:
- Around line 248-253: The pointer association block under "if (hyper_cleaning)
then" must be guarded by the non-IGR path and must avoid leaking an existing
allocation: change the condition to check ".not. igr" (e.g. if (hyper_cleaning
.and. .not. igr) then) so q_cons_qp%vf(psi_idx)%sf is only aliased when it's
allocated, and before doing the association for q_prim_qp%vf(psi_idx)%sf ensure
you either deallocate q_prim_qp%vf(psi_idx)%sf or skip its earlier allocation in
the generic allocation loop so the subsequent pointer assignment does not leak;
keep using the same GPU directives ($:GPU_ENTER_DATA) but only inside this
guarded branch.
In `@src/simulation/m_riemann_solvers.fpp`:
- Around line 649-652: The current hyper_cleaning branch incorrectly sets
c_fast%L = min(c_fast%L, -hyper_cleaning_speed) which forces a negative left
bound and can shrink the HLL envelope; instead ensure the magnitude is applied
as a positive bound on speeds by using max for c_fast%L and min for c_fast%R (or
apply the magnitude directly when computing s_L/s_R). Update the hyper_cleaning
block around c_fast%L, c_fast%R so c_fast%L = max(c_fast%L,
-hyper_cleaning_speed) is not used — rather clamp magnitudes correctly (e.g.,
c_fast%L = -max(abs(c_fast%L), hyper_cleaning_speed) or set c_fast%L =
-hyper_cleaning_speed with c_fast%R = hyper_cleaning_speed) ensuring s_L = vel -
c_fast and s_R = vel + c_fast produce a proper HLL envelope; modify the code
referencing c_fast%L/c_fast%R in this block accordingly.
In `@src/simulation/m_time_steppers.fpp`:
- Around line 637-647: The decay loop runs on the host but q_cons_ts is on the
device in GPU builds; wrap the triple loop in the GPU loop macro (e.g., the
project's DEVICE/GPU_LOOP macro) and precompute the decay factor into a
loop-local variable (e.g., real :: decay = exp(-dt/hyper_cleaning_tau)) before
entering the GPU loop; then replace repeated exp() calls with that decay and
mark any loop-local scalars (like decay and the loop indices if required) as
private so the loops over j,k,l that update q_cons_ts(1)%vf(psi_idx)%sf(j,k,l)
execute correctly and efficiently on the device.
♻️ Duplicate comments (2)
src/simulation/m_riemann_solvers.fpp (2)
5027-5041: Same concern as the y-direction Powell copy.
Line 5027 duplicates the full-array copy forpowell; please ensure the copied components are fully initialized. (See Line 4971 comment.)
5083-5096: Same concern as the y-direction Powell copy.
Line 5083 duplicates the full-array copy forpowell; please ensure the copied components are fully initialized. (See Line 4971 comment.)
🧹 Nitpick comments (7)
case_new/2D_orszag_tang_weno5m_nohyper/case.py (2)
12-15: Minor type inconsistency in domain bounds.
x_domain%begandy_domain%begare integers (0) whilex_domain%endandy_domain%endare floats (1.0). For consistency and to avoid potential type-handling issues downstream, consider using0.0for the begin values.Suggested fix
# Computational Domain Parameters - "x_domain%beg": 0, + "x_domain%beg": 0.0, "x_domain%end": 1.0, - "y_domain%beg": 0, + "y_domain%beg": 0.0, "y_domain%end": 1.0,
32-32: Inconsistent quote style.This line uses single quotes
'T'while all other boolean string values use double quotes"T". While functionally equivalent in Python, consistent formatting improves readability.Suggested fix
- "mapped_weno": 'T', + "mapped_weno": "T",src/simulation/m_mhd.fpp (1)
39-39: Consider updating error message to match renamed module.The error message still references "Powell correction" but the module has been renamed to "clean". Consider updating for consistency.
♻️ Suggested fix
- if (n == 0) call s_mpi_abort('Fatal Error: Powell correction is not applicable for 1D') + if (n == 0) call s_mpi_abort('Fatal Error: MHD divergence cleaning is not applicable for 1D')src/simulation/m_data_output.fpp (1)
1494-1498: Remove commented-out placeholder code or implement the TODO.This commented-out block is a non-functional stub that doesn't reduce any meaningful variable. Either implement the actual hyper_cleaning reduction logic or remove the placeholder to avoid dead code accumulation.
♻️ Suggested fix (remove placeholder)
if (cont_damage) then tmp = damage_state call s_mpi_allreduce_sum(tmp, damage_state) end if - - ! if (hyper_cleaning) then - ! tmp = 0._wp - ! call s_mpi_allreduce_sum(tmp, tmp) ! TODO - ! end if end ifIf hyper_cleaning requires probe output support, would you like me to help implement the actual reduction logic?
src/pre_process/m_global_parameters.fpp (1)
378-379: Initializepsi_idxwith a safe default.
psi_idxis new but never assigned in the defaults block. Setting it todflt_int(or-1) reduces accidental use whenhyper_cleaning = .false..♻️ Proposed fix
- cont_damage = .false. - hyper_cleaning = .false. + cont_damage = .false. + hyper_cleaning = .false. + psi_idx = dflt_intsrc/simulation/m_global_parameters.fpp (1)
1166-1169: Consider adding else clause forpsi_idxinitialization.Unlike the post-processing module (which sets
psi_idx = dflt_intwhenhyper_cleaningis disabled), this simulation module leavespsi_idxuninitialized when the feature is off. While this matches the pattern fordamage_idx, it could lead to undefined behavior if code inadvertently referencespsi_idxwithout first checkinghyper_cleaning.This is consistent with existing code, so it's likely intentional, but worth noting for awareness.
♻️ Optional: Add else clause for consistency
if (hyper_cleaning) then psi_idx = sys_size + 1 sys_size = psi_idx + else + psi_idx = dflt_int end ifcase_new/2D_orszag_tang_weno3m/case.py (1)
5-78: Consider extracting shared case defaults to reduce duplication.
These newcase_new/*/case.pyfiles are nearly identical; a small helper to build a base dict and override per-case knobs would reduce drift and make future updates safer.
| # Computational Domain: [0,1]×[0,1] but n=1 → effectively 1D in x | ||
| "x_domain%beg":-0.5, | ||
| "x_domain%end": 0.5, | ||
| "y_domain%beg":-0.5, | ||
| "y_domain%end": 0.5, | ||
| "m": 100, | ||
| "n": 100, |
There was a problem hiding this comment.
Comment doesn’t match the grid settings (and uses a non-ASCII ×).
Line 12 claims n=1, but n=100; please align the comment with the actual setup (or set n=1 if that was the intent).
✏️ Possible tweak
- # Computational Domain: [0,1]×[0,1] but n=1 → effectively 1D in x
+ # Computational Domain: [-0.5,0.5] x [-0.5,0.5]; IC is x-only so the solution stays 1D🧰 Tools
🪛 Ruff (0.14.13)
12-12: Comment contains ambiguous × (MULTIPLICATION SIGN). Did you mean x (LATIN SMALL LETTER X)?
(RUF003)
🤖 Prompt for AI Agents
In `@case_new/gauss_div_pulse_2/case.py` around lines 12 - 18, The comment
claiming a 1D domain ("n=1" and using a non-ASCII ×) is inconsistent with the
actual grid settings where "m": 100 and "n": 100; either update the comment to
accurately describe a 2D computational domain (e.g., replace the non-ASCII ×
with 'x' and state m=100, n=100 over x_domain%beg/end and y_domain%beg/end) or
change the grid parameter "n" to 1 if the intent was a 1D x-only run; locate the
grid dict entries ("x_domain%beg", "x_domain%end", "y_domain%beg",
"y_domain%end", "m", "n") and make the corresponding comment/text consistent
with the chosen configuration.
| # --- Patch 1: Gaussian Divergence Pulse IC --- | ||
| "patch_icpp(1)%geometry": 3, # Cartesian | ||
| "patch_icpp(1)%x_centroid": 0.0, | ||
| "patch_icpp(1)%y_centroid": 0.0, | ||
| "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)%alpha_rho(1)": 1.0, | ||
| "patch_icpp(1)%alpha(1)": 1.0, | ||
| "patch_icpp(1)%Bx": 0.0, | ||
| "patch_icpp(1)%By": 0.0, | ||
| "patch_icpp(1)%Bz": 0.0, |
There was a problem hiding this comment.
Gaussian divergence pulse IC isn’t activated.
The hardcoded IC expects hcid=260 plus patch_icpp%a(2) and a(3); without them this case stays uniform and won’t validate divergence cleaning. Please wire the hcid and parameters.
🛠️ Suggested wiring (set values to the ones used in your validation)
# --- Patch 1: Gaussian Divergence Pulse IC ---
+ "patch_icpp(1)%hcid": 260,
+ "patch_icpp(1)%a(2)": 0.1, # eps_mhd (∂Bx/∂x amplitude) — adjust to intended value
+ "patch_icpp(1)%a(3)": 0.05, # sigma — adjust to intended value
"patch_icpp(1)%geometry": 3, # Cartesian📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| # --- Patch 1: Gaussian Divergence Pulse IC --- | |
| "patch_icpp(1)%geometry": 3, # Cartesian | |
| "patch_icpp(1)%x_centroid": 0.0, | |
| "patch_icpp(1)%y_centroid": 0.0, | |
| "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)%alpha_rho(1)": 1.0, | |
| "patch_icpp(1)%alpha(1)": 1.0, | |
| "patch_icpp(1)%Bx": 0.0, | |
| "patch_icpp(1)%By": 0.0, | |
| "patch_icpp(1)%Bz": 0.0, | |
| # --- Patch 1: Gaussian Divergence Pulse IC --- | |
| "patch_icpp(1)%hcid": 260, | |
| "patch_icpp(1)%a(2)": 0.1, # eps_mhd (∂Bx/∂x amplitude) — adjust to intended value | |
| "patch_icpp(1)%a(3)": 0.05, # sigma — adjust to intended value | |
| "patch_icpp(1)%geometry": 3, # Cartesian | |
| "patch_icpp(1)%x_centroid": 0.0, | |
| "patch_icpp(1)%y_centroid": 0.0, | |
| "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)%alpha_rho(1)": 1.0, | |
| "patch_icpp(1)%alpha(1)": 1.0, | |
| "patch_icpp(1)%Bx": 0.0, | |
| "patch_icpp(1)%By": 0.0, | |
| "patch_icpp(1)%Bz": 0.0, |
🤖 Prompt for AI Agents
In `@case_new/gauss_div_pulse_2/case.py` around lines 56 - 70, Add the missing
hard-coded initial condition ID and parameters so the Gaussian divergence pulse
is actually activated: set hcid to 260 and add the two missing patch_icpp
parameters patch_icpp(1)%a(2) and patch_icpp(1)%a(3) (use the same numeric
values used in your validation run), alongside the existing patch_icpp(1)%a(1)
entry so the IC block matches the hardcoded expectation.
| # ----------------------------------------------------------------------------- | ||
| # 2.1 Gaussian Divergence Pulse (1D via uniform y-slice) | ||
| # ----------------------------------------------------------------------------- | ||
| print( | ||
| json.dumps( | ||
| { | ||
| # Logistics | ||
| "run_time_info": "T", | ||
| # Computational Domain: | ||
| "x_domain%beg": -1.0, | ||
| "x_domain%end": 2.0, | ||
| "y_domain%beg":-1.0, | ||
| # "y_domain%end": 2.0/256.0, | ||
| "y_domain%end": 2.0, | ||
| "m": 256*3, | ||
| "n": 256*3, | ||
| "dt": 0.08/math.sqrt(5.0)/170.0/2.0, | ||
| "t_step_start": 0, | ||
| "t_step_stop": 170*2, | ||
| "t_step_save": 10*2, |
There was a problem hiding this comment.
Mismatch between filename, description, and domain setup.
This file is named slanted_shock_tube but the header and hcid indicate a Gaussian divergence pulse, and the domain/grid are 2D rather than a 1D y-slice. Please align (a) the filename, (b) the header/comments, and (c) the y-extent + n values with the intended test to avoid running the wrong validation case.
🤖 Prompt for AI Agents
In `@case_new/slanted_shock_tube/case.py` around lines 5 - 24, The file and meta
should be consistent: decide whether this is a 1D Gaussian divergence pulse or a
2D slanted shock tube and update the header/comment, hcid, and domain/grid
accordingly. If you mean a 1D y-slice Gaussian pulse, change the header from
"Gaussian Divergence Pulse (1D via uniform y-slice)" to match the filename or
vice versa, set "y_domain%end" to a very small extent (e.g. equal to a single
cell), reduce "n" to 1 (or minimal), keep "m" as the long x resolution, and
ensure hcid reflects the Gaussian pulse case; alternatively, if this is the
slanted shock tube test, update the header/comment and hcid to describe a 2D
slanted shock tube and keep the current 2D "y_domain%end" and "n"/"m" values.
Locate and edit the JSON keys "y_domain%end", "n", "m", the header comment
block, and the hcid identifier to make them consistent.
| 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 | ||
|
|
||
| ! B-field | ||
| q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp + C_mhd*erf((x_cc(i) - 0.5_wp)/sigma) | ||
|
|
There was a problem hiding this comment.
Guard against σ=0 in the Gaussian divergence pulse.
sigma is used in a denominator; if unset/zero it will produce NaNs. Please validate it before use (per .fpp guidelines).
🛡️ Suggested safety check
eps_mhd = patch_icpp(patch_id)%a(2)
sigma = patch_icpp(patch_id)%a(3)
+ @:ASSERT(sigma > 0._wp, "Gaussian divergence pulse requires patch_icpp%a(3) (sigma) > 0")
C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp🤖 Prompt for AI Agents
In `@src/common/include/2dHardcodedIC.fpp` around lines 242 - 254, In the Gaussian
divergence pulse case (case 260) validate patch_icpp(patch_id)%a(3) (sigma)
before using it: check for zero or near-zero and either set a safe minimum
(e.g., sigma_min) or raise/log an error and abort per .fpp guidelines; then
compute C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp and use that safe sigma for the
erf denominator and any other divisions so q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1
+ C_mhd*erf((x_cc(i)-0.5_wp)/sigma) cannot produce NaNs. Ensure the check is
colocated with the lines that set eps_mhd, sigma and C_mhd in the Gaussian
Divergence Pulse block.
| if (hyper_cleaning) then | ||
| q_prim_qp%vf(psi_idx)%sf => & | ||
| q_cons_qp%vf(psi_idx)%sf | ||
| $:GPU_ENTER_DATA(copyin='[q_prim_qp%vf(psi_idx)%sf]') | ||
| $:GPU_ENTER_DATA(attach='[q_prim_qp%vf(psi_idx)%sf]') | ||
| end if |
There was a problem hiding this comment.
Guard psi aliasing to avoid unallocated access and leaks.
Line 248 currently runs even when igr is true, but q_cons_qp%vf(psi_idx)%sf is only allocated inside the .not. igr path. Also, in the non-IGR path q_prim_qp%vf(psi_idx)%sf is already allocated in the generic allocation loop, so re-aliasing here leaks that allocation. Please guard the block with .not. igr and deallocate (or skip allocating) q_prim_qp%vf(psi_idx)%sf before the pointer association.
🐛 Proposed fix
- if (hyper_cleaning) then
+ if (.not. igr .and. hyper_cleaning) then
+ if (associated(q_prim_qp%vf(psi_idx)%sf)) then
+ @:DEALLOCATE(q_prim_qp%vf(psi_idx)%sf)
+ end if
q_prim_qp%vf(psi_idx)%sf => &
q_cons_qp%vf(psi_idx)%sf
$:GPU_ENTER_DATA(copyin='[q_prim_qp%vf(psi_idx)%sf]')
$:GPU_ENTER_DATA(attach='[q_prim_qp%vf(psi_idx)%sf]')
end if🤖 Prompt for AI Agents
In `@src/simulation/m_rhs.fpp` around lines 248 - 253, The pointer association
block under "if (hyper_cleaning) then" must be guarded by the non-IGR path and
must avoid leaking an existing allocation: change the condition to check ".not.
igr" (e.g. if (hyper_cleaning .and. .not. igr) then) so q_cons_qp%vf(psi_idx)%sf
is only aliased when it's allocated, and before doing the association for
q_prim_qp%vf(psi_idx)%sf ensure you either deallocate q_prim_qp%vf(psi_idx)%sf
or skip its earlier allocation in the generic allocation loop so the subsequent
pointer assignment does not leak; keep using the same GPU directives
($:GPU_ENTER_DATA) but only inside this guarded branch.
| 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 |
There was a problem hiding this comment.
Fix sign of c_fast bounding for hyper_cleaning.
Line 650 drives c_fast%L negative, which flips the left-going wave bound used in s_L = vel - c_fast and can shrink the HLL envelope. If hyper_cleaning_speed is a magnitude, bound by max (or apply to s_L/s_R directly).
🔧 Suggested fix
- 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
+ if (hyper_cleaning) then ! mhd
+ c_fast%L = max(c_fast%L, hyper_cleaning_speed)
+ c_fast%R = max(c_fast%R, hyper_cleaning_speed)
+ end if📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| 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 | |
| if (hyper_cleaning) then ! mhd | |
| c_fast%L = max(c_fast%L, hyper_cleaning_speed) | |
| c_fast%R = max(c_fast%R, hyper_cleaning_speed) | |
| end if |
🤖 Prompt for AI Agents
In `@src/simulation/m_riemann_solvers.fpp` around lines 649 - 652, The current
hyper_cleaning branch incorrectly sets c_fast%L = min(c_fast%L,
-hyper_cleaning_speed) which forces a negative left bound and can shrink the HLL
envelope; instead ensure the magnitude is applied as a positive bound on speeds
by using max for c_fast%L and min for c_fast%R (or apply the magnitude directly
when computing s_L/s_R). Update the hyper_cleaning block around c_fast%L,
c_fast%R so c_fast%L = max(c_fast%L, -hyper_cleaning_speed) is not used — rather
clamp magnitudes correctly (e.g., c_fast%L = -max(abs(c_fast%L),
hyper_cleaning_speed) or set c_fast%L = -hyper_cleaning_speed with c_fast%R =
hyper_cleaning_speed) ensuring s_L = vel - c_fast and s_R = vel + c_fast produce
a proper HLL envelope; modify the code referencing c_fast%L/c_fast%R in this
block accordingly.
| if (hyper_cleaning) then | ||
| do l = 0, p | ||
| do k = 0, n | ||
| do j = 0, m | ||
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & | ||
| q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & | ||
| exp(-dt/hyper_cleaning_tau) | ||
| end do | ||
| end do | ||
| end do | ||
| end if |
There was a problem hiding this comment.
Decay loop must run on device (missing GPU macro).
q_cons_ts lives on the device in GPU builds, so the host-side loop can be incorrect and will definitely be slow. Wrap this in the GPU loop macro and precompute the decay factor. As per coding guidelines, GPU loops should use the macros and declare loop-local variables via private when needed.
✅ Proposed fix (GPU loop + precomputed decay)
- if (hyper_cleaning) then
- do l = 0, p
- do k = 0, n
- do j = 0, m
- q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = &
- q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* &
- exp(-dt/hyper_cleaning_tau)
- end do
- end do
- end do
- end if
+ if (hyper_cleaning) then
+ psi_decay = exp(-dt/hyper_cleaning_tau)
+ $:GPU_PARALLEL_LOOP(collapse=3)
+ do l = 0, p
+ do k = 0, n
+ do j = 0, m
+ q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = &
+ q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) * psi_decay
+ end do
+ end do
+ end do
+ $:END_GPU_PARALLEL_LOOP()
+ end if🤖 Prompt for AI Agents
In `@src/simulation/m_time_steppers.fpp` around lines 637 - 647, The decay loop
runs on the host but q_cons_ts is on the device in GPU builds; wrap the triple
loop in the GPU loop macro (e.g., the project's DEVICE/GPU_LOOP macro) and
precompute the decay factor into a loop-local variable (e.g., real :: decay =
exp(-dt/hyper_cleaning_tau)) before entering the GPU loop; then replace repeated
exp() calls with that decay and mark any loop-local scalars (like decay and the
loop indices if required) as private so the loops over j,k,l that update
q_cons_ts(1)%vf(psi_idx)%sf(j,k,l) execute correctly and efficiently on the
device.
![cubic-dev-ai[bot]](https://avatars.githubusercontent.com/in/1082092?s=60&v=4){kind=small}
There was a problem hiding this comment.
4 issues found across 28 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="case_new/gauss_div_pulse_2/case.py">
<violation number="1" location="case_new/gauss_div_pulse_2/case.py:18">
P2: The case is documented as a 1D y-slice, but n is set to 100, making it 2D. If the intent is a 1D slice, set n to 1 to match the stated setup.</violation>
</file>
<file name="src/pre_process/m_start_up.fpp">
<violation number="1" location="src/pre_process/m_start_up.fpp:847">
P2: The new hyper_cleaning block hard-codes a Gaussian psi and overwrites any psi set by the initial-condition generator or restart files. This makes hyper_cleaning runs ignore user-provided psi values and breaks restart fidelity. Consider removing this block or gating it behind a dedicated input option so default runs preserve the provided initial condition.</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 but never declared/copied to the GPU, yet it is used inside GPU kernels. This can leave psi_idx undefined on device and cause out-of-bounds indexing when hyper_cleaning is enabled. Add psi_idx to the GPU_DECLARE and GPU_UPDATE lists used for scalar indices.</violation>
</file>
<file name="src/simulation/m_rhs.fpp">
<violation number="1" location="src/simulation/m_rhs.fpp:249">
P2: q_prim_qp%vf(psi_idx)%sf is already allocated earlier when hyper_cleaning is enabled, and this pointer reassignment discards that allocation, leaking a full-sized scalar field. This can waste significant memory on large grids. Consider deallocating before aliasing or skipping the initial allocation for psi_idx.</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.
| "y_domain%beg":-0.5, | ||
| "y_domain%end": 0.5, | ||
| "m": 100, | ||
| "n": 100, |
There was a problem hiding this comment.
P2: The case is documented as a 1D y-slice, but n is set to 100, making it 2D. If the intent is a 1D slice, set n to 1 to match the stated setup.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At case_new/gauss_div_pulse_2/case.py, line 18:
<comment>The case is documented as a 1D y-slice, but n is set to 100, making it 2D. If the intent is a 1D slice, set n to 1 to match the stated setup.</comment>
<file context>
@@ -0,0 +1,76 @@
+ "y_domain%beg":-0.5,
+ "y_domain%end": 0.5,
+ "m": 100,
+ "n": 100,
+ "p": 0,
+ "dt": 0.005,
</file context>
|
|
||
| call s_generate_initial_condition() | ||
|
|
||
| ! hard-coded psi |
There was a problem hiding this comment.
P2: The new hyper_cleaning block hard-codes a Gaussian psi and overwrites any psi set by the initial-condition generator or restart files. This makes hyper_cleaning runs ignore user-provided psi values and breaks restart fidelity. Consider removing this block or gating it behind a dedicated input option so default runs preserve the provided initial condition.
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 847:
<comment>The new hyper_cleaning block hard-codes a Gaussian psi and overwrites any psi set by the initial-condition generator or restart files. This makes hyper_cleaning runs ignore user-provided psi values and breaks restart fidelity. Consider removing this block or gating it behind a dedicated input option so default runs preserve the provided initial condition.</comment>
<file context>
@@ -842,6 +844,16 @@ contains
call s_generate_initial_condition()
+ ! hard-coded psi
+ if (hyper_cleaning) then
+ do j = 0, m
</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 but never declared/copied to the GPU, yet it is used inside GPU kernels. This can leave psi_idx undefined on device and cause out-of-bounds indexing when hyper_cleaning is enabled. Add psi_idx to the GPU_DECLARE and GPU_UPDATE lists used for scalar indices.
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 but never declared/copied to the GPU, yet it is used inside GPU kernels. This can leave psi_idx undefined on device and cause out-of-bounds indexing when hyper_cleaning is enabled. Add psi_idx to the GPU_DECLARE and GPU_UPDATE lists used for scalar indices.</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>
| end if | ||
|
|
||
| if (hyper_cleaning) then | ||
| q_prim_qp%vf(psi_idx)%sf => & |
There was a problem hiding this comment.
P2: q_prim_qp%vf(psi_idx)%sf is already allocated earlier when hyper_cleaning is enabled, and this pointer reassignment discards that allocation, leaking a full-sized scalar field. This can waste significant memory on large grids. Consider deallocating before aliasing or skipping the initial allocation for psi_idx.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_rhs.fpp, line 249:
<comment>q_prim_qp%vf(psi_idx)%sf is already allocated earlier when hyper_cleaning is enabled, and this pointer reassignment discards that allocation, leaking a full-sized scalar field. This can waste significant memory on large grids. Consider deallocating before aliasing or skipping the initial allocation for psi_idx.</comment>
<file context>
@@ -245,6 +245,13 @@ contains
end if
+ if (hyper_cleaning) then
+ q_prim_qp%vf(psi_idx)%sf => &
+ q_cons_qp%vf(psi_idx)%sf
+ $:GPU_ENTER_DATA(copyin='[q_prim_qp%vf(psi_idx)%sf]')
</file context>
|
Looks like it passes tests. I stupidly didn't move the new cases to examples. Once I generate more golden files and it passes, we can rerun the test suite and hopefully this will be wrapped up. |
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
Implements hyperbolic divergence cleaning for MHD to replace Powell's method
Adds new state variable
psi_idxfor hyperbolic cleaning throughout codebaseIntroduces three new MHD test cases: Gaussian divergence pulse, tilted shock tube, 2D rotor
Adds Orszag-Tang vortex validation cases with multiple WENO schemes
Includes hard-coded initial condition for
psiwith Gaussian profileModifies Riemann solver to include hyperbolic cleaning flux contributions
Diagram Walkthrough
File Walkthrough
19 files
Add new MHD test cases with hardcoded ICsHandle psi variable in primitive-conservative conversionInclude psi in database output variable countAdd hyper_cleaning flag and psi_idx parameterBroadcast hyper_cleaning flag across MPI processesWrite psi variable to formatted database outputAdd hyper_cleaning flag and psi_idx to pre-processorInitialize psi to zero in initial condition setupBroadcast hyper_cleaning flag in pre-processor MPIHard-code Gaussian psi initialization and read hyper_cleaningAdd psi to simulation data output (commented TODO)Add hyper_cleaning parameters and GPU declarationsRename Powell module to generic cleaning moduleBroadcast hyper_cleaning parameters across MPIAdd hyperbolic cleaning flux and wave speed modificationsInitialize hyperbolic cleaning module and read parametersAllocate psi and apply exponential decay time steppingUpdate magnetic vortex case to use new case IDAdd hyper_cleaning parameters to case dictionary8 files
Orszag-Tang vortex test with WENO3 mapped schemeOrszag-Tang vortex test with WENO5 mapped schemeOrszag-Tang vortex without hyperbolic cleaningOrszag-Tang vortex test with WENO5Z schemeOrszag-Tang vortex test with WENO7 mapped scheme2D MHD rotor problem test case configurationGaussian divergence pulse test case setupTilted MHD shock tube test case configuration1 files
Summary by CodeRabbit
New Features
Updates
✏️ Tip: You can customize this high-level summary in your review settings.