[SPH][NIMHD] Add Non-Ideal MHD terms#1765
Open
y-lapeyre wants to merge 13 commits into
Open
Conversation
Contributor
|
Thanks @y-lapeyre for opening this PR! You can do multiple things directly here: Once the workflow completes a message will appear displaying informations related to the run. Also the PR gets automatically reviewed by gemini, you can: |
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Code Review
This pull request introduces support for Non-Ideal Magnetohydrodynamics (NIMHD) by adding Ohmic resistivity, Hall effect, and Ambipolar Diffusion terms. Key changes include the addition of NIMHD configuration parameters, mathematical kernels for current and heating terms, and a new CFL timestep constraint for Ambipolar Diffusion. Review feedback highlighted a critical architectural issue regarding the pairwise calculation of electric current, which should be a field quantity requiring a two-pass approach. Additionally, several potential division-by-zero bugs were identified in the mathematical kernels and timestep calculations, along with a likely incorrect gradient index in the induction term and misleading argument names in the Python bindings.
Collaborator
Author
|
/gemini review |
Contributor
There was a problem hiding this comment.
Code Review
This pull request introduces support for Non-Ideal Magnetohydrodynamics (NIMHD) in the SPH solver, adding configuration parameters, mathematical formulations for current density and resistivity, and the corresponding derivative updates. However, several critical issues must be addressed: a bug where the current density J_a is passed by value instead of reference in MagCurrentJ_sum (discarding accumulated values), a race condition and event synchronization bug in the CFL timestep calculation due to variable shadowing, and potential out-of-bounds reads in compute_J from iterating over ghost particles that lack neighbor lists. Additionally, minor improvements are recommended, such as removing unused parameters in B_NI_terms, cleaning up a temporary comment, and renaming a Python binding parameter for consistency.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
Contributor
Workflow reportworkflow report corresponding to commit 8b64f74 Pre-commit check reportPre-commit check: ✅ Test pipeline can run. Clang-tidy diff reportDoxygen diff with
|
Collaborator
Author
|
/gemini review |
Codecov Report❌ Patch coverage is 📢 Thoughts on this report? Let us know! |
Contributor
There was a problem hiding this comment.
Code Review
This pull request introduces Non-Ideal Magnetohydrodynamics (NIMHD) support to the SPH solver, including configurations for Ohmic resistivity, Hall effect, and Ambipolar diffusion, as well as the calculation of electric current density and non-ideal heating/induction terms. The review feedback highlights critical issues with SYCL event shadowing that can lead to concurrent kernel execution and data races on the timestep buffer. Other recommendations include evaluating the NIMHD flag at compile-time, skipping self-interactions in SPH neighbor loops, removing commented-out code, and resolving a developer query about the correct kernel gradient.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
| // first, finish the outer kernel and write in cfl_dt | ||
| // then, compute psi_dt and write in cfl_dt | ||
| depends_list.add_event(e); | ||
| auto e = q.submit(depends_list, [&](sycl::handler &cgh) { |
Contributor
There was a problem hiding this comment.
The variable e is shadowed here using auto e = ... inside the if (has_psi_field) block. This prevents the outer e variable (from the first CFL kernel) from being updated. Consequently, when entering the if (do_NIMHD) block, depends_list.add_event(e) will add the first CFL kernel's event instead of the psi cleaning kernel's event. This causes the NIMHD kernel and the psi cleaning kernel to run concurrently, leading to a data race on cfl_dt and potential undefined behavior. Removing auto to reassign the outer e variable resolves this issue.
Suggested change
| auto e = q.submit(depends_list, [&](sycl::handler &cgh) { | |
| e = q.submit(depends_list, [&](sycl::handler &cgh) { |
|
|
||
| if (do_NIMHD) { | ||
| depends_list.add_event(e); | ||
| auto e = q.submit(depends_list, [&](sycl::handler &cgh) { |
Contributor
There was a problem hiding this comment.
Similarly, e is shadowed here using auto e = ... inside the if (do_NIMHD) block. This prevents the outer e from being updated to the NIMHD kernel's event, meaning the subsequent complete_event_state(e) calls on lines 2614-2616 may not wait for the NIMHD kernel to finish. Removing auto to reassign the outer e variable resolves this issue.
Suggested change
| auto e = q.submit(depends_list, [&](sycl::handler &cgh) { | |
| e = q.submit(depends_list, [&](sycl::handler &cgh) { |
| u32 iB_on_rho_interf = ghost_layout.get_field_idx<Tvec>("B/rho"); | ||
| u32 ipsi_on_ch_interf = ghost_layout.get_field_idx<Tscal>("psi/ch"); | ||
|
|
||
| bool do_NIMHD = solver_config.do_NIMHD(); |
Contributor
There was a problem hiding this comment.
do_NIMHD is evaluated at runtime using solver_config.do_NIMHD(). Since MHD_mode is a template parameter of update_derivs_MHD_impl, this can be evaluated at compile-time as constexpr bool do_NIMHD = (MHD_mode == shamrock::sph::mhd::MHDType::NonIdeal);. This allows the compiler to optimize away the Non-Ideal branches when compiling the Ideal MHD version.
Suggested change
| bool do_NIMHD = solver_config.do_NIMHD(); | |
| constexpr bool do_NIMHD = (MHD_mode == shamrock::sph::mhd::MHDType::NonIdeal); |
Comment on lines
+822
to
+825
| particle_looper.for_each_object(id_a, [&](u32 id_b) { | ||
| Tvec dr = xyz_a - xyz[id_b]; | ||
| Tscal rab2 = sycl::dot(dr, dr); | ||
| Tscal h_b = hpart[id_b]; |
Contributor
There was a problem hiding this comment.
In SPH neighbor loops, self-interaction (id_a == id_b) can be skipped early to avoid redundant calculations (since B_a - B_b is zero, the contribution to J_a is mathematically zero).
particle_looper.for_each_object(id_a, [&](u32 id_b) {
if (id_a == id_b) {
return;
}
Tvec dr = xyz_a - xyz[id_b];
Tscal rab2 = sycl::dot(dr, dr);
Tscal h_b = hpart[id_b];
Comment on lines
+1079
to
+1082
| particle_looper.for_each_object(id_a, [&](u32 id_b) { | ||
| Tvec dr = xyz_a - xyz[id_b]; | ||
| Tscal rab2 = sycl::dot(dr, dr); | ||
| Tscal h_b = hpart[id_b]; |
Contributor
There was a problem hiding this comment.
In SPH neighbor loops, self-interaction (id_a == id_b) can be skipped early to avoid redundant calculations.
particle_looper.for_each_object(id_a, [&](u32 id_b) {
if (id_a == id_b) {
return;
}
Tvec dr = xyz_a - xyz[id_b];
Tscal rab2 = sycl::dot(dr, dr);
Tscal h_b = hpart[id_b];
Comment on lines
+461
to
+464
| // Tvec J_a = MagCurrentJ<Tvec, Tscal, MHD_mode>( | ||
| // pmass, B_a, B_b, r_ab_unit * dWab_a, sub_fact_a, mu_0); | ||
| // Tvec J_b = MagCurrentJ<Tvec, Tscal, MHD_mode>( | ||
| // pmass, B_a, B_b, r_ab_unit * dWab_b, sub_fact_b, mu_0); |
Contributor
There was a problem hiding this comment.
Commented-out code for J_a and J_b calculation should be removed.
| += v_ab * dB_on_rho_induction_term(pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_b); | ||
| dB_on_rho_dt += v_ab | ||
| * dB_on_rho_induction_term( | ||
| pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_a); // @@@ dWab_b ? |
Contributor
There was a problem hiding this comment.
The developer left a query // @@@ dWab_b ? regarding whether dWab_a or dWab_b should be used. Since the induction term is evaluated for particle a using particle a's properties (rho_a_sq, B_a, omega_a), the kernel gradient must be evaluated with respect to particle a's smoothing length (h_a), which corresponds to dWab_a. Thus, dWab_a is mathematically correct, and the comment can be removed.
Suggested change
| pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_a); // @@@ dWab_b ? | |
| pmass, rho_a_sq, B_a, omega_a, r_ab_unit * dWab_a); |
tdavidcl
changed the title
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Add and enable the use of no-ideal MHD in the MHD solver, using the structure I previously set out for it.
Tests:
Could be merged as is and fixed with a 2nd PR that adds the wave damping CI test.