[SPH] 1Pn disc-sink interactions - #1923
Conversation
|
Thanks @AugustinDart 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}
There was a problem hiding this comment.
Code Review
This pull request introduces a new 1PN external force model (AddForce1PN) and updates the Lense-Thirring force model to support configurable central positions and velocities. Feedback focuses on preventing potential division-by-zero errors in both force modules when a particle is exactly at the central position. Additionally, it is recommended to update the particle accretion logic to respect the new Lense-Thirring central position, and to fix a minor variable alignment issue.
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.
| Tvec r_a = xyz[gid] - cpos; | ||
| Tvec v_a = vxyz[gid] - cvel; | ||
| Tscal abs_ra = sycl::length(r_a); | ||
| Tscal abs_ra_2 = abs_ra * abs_ra; | ||
| Tscal abs_ra_3 = abs_ra_2 * abs_ra; |
There was a problem hiding this comment.
Similar to the 1PN force, if a particle is exactly at the central position cpos, abs_ra will be 0, leading to a division by zero when calculating omega_a. Adding a guard to skip calculation when abs_ra is extremely small prevents numerical instability. Additionally, since this relative position and distance calculation logic is duplicated across multiple force modules, consider refactoring it into a shared helper function or lambda to improve readability and maintainability.
| Tvec r_a = xyz[gid] - cpos; | |
| Tvec v_a = vxyz[gid] - cvel; | |
| Tscal abs_ra = sycl::length(r_a); | |
| Tscal abs_ra_2 = abs_ra * abs_ra; | |
| Tscal abs_ra_3 = abs_ra_2 * abs_ra; | |
| Tvec r_a = xyz[gid] - cpos; | |
| Tscal abs_ra = sycl::length(r_a); | |
| if (abs_ra < 1e-12) { | |
| return; | |
| } | |
| Tvec v_a = vxyz[gid] - cvel; | |
| Tscal abs_ra_2 = abs_ra * abs_ra; | |
| Tscal abs_ra_3 = abs_ra_2 * abs_ra; |
References
- Refactor duplicated logic into a helper function or lambda to improve readability and maintainability.
| auto set_central_pos | ||
| = register_constant_set<Tvec>(solver_graph, prefix_central_pos, [&]() { | ||
| return Tvec{0, 0, 0}; // no support for offset yet | ||
| return ext_force->central_pos; | ||
| }); |
There was a problem hiding this comment.
Since LenseThirring now supports a non-zero central_pos, please make sure to also update point_mass_accrete_particles() (around line 678 in this file) to use ext_force->central_pos instead of the hardcoded {0, 0, 0}. Otherwise, accretion will still occur at the origin instead of the actual central position.
| Tscal c = edges.constant_c.data; | ||
| Tscal cmass = edges.central_mass.data; | ||
| Tvec cpos = edges.central_pos.data; | ||
| Tvec cvel = edges.central_vel.data; |
There was a problem hiding this comment.
The variable cvel is not aligned with the other variable declarations in this block. Aligning the assignment operator improves readability.
| Tvec cvel = edges.central_vel.data; | |
| Tvec cvel = edges.central_vel.data; |
📝 WalkthroughWalkthroughAdds the ChangesExternal force models
Estimated code review effort: 4 (Complex) | ~45 minutes Sequence Diagram(s)sequenceDiagram
participant PythonConfig
participant SolverConfig
participant ExtForceConfig
participant ExternalForces
participant AddForce1PN
participant AccelerationField
PythonConfig->>SolverConfig: add_ext_force_1pn(central_mass, central_pos, central_vel)
SolverConfig->>ExtForceConfig: add_1pn(...)
ExtForceConfig-->>ExternalForces: PN_1PN variant
ExternalForces->>AddForce1PN: wire constants and field edges
AddForce1PN->>AccelerationField: accumulate 1PN acceleration
Suggested reviewers: 🚥 Pre-merge checks | ✅ 4 | ❌ 1❌ Failed checks (1 warning)
✅ Passed checks (4 passed)
✨ Finishing Touches🧪 Generate unit tests (beta)
Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp`:
- Around line 65-110: Use a single dimensionally consistent radius clamp in the
1PN kernel around the acceleration calculation in AddForce1PN: clamp r once into
r_safe, then derive and consistently use r_safe² and r_safe³ for all
radius-dependent denominators, preserving the existing acceleration behavior. In
AddForceLenseThirring, clamp abs_ra before computing its powers so the
central-position case cannot produce NaN acceleration; apply the change in
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp lines
65-110 and
src/shammodels/common/include/shammodels/common/modules/AddForceLenseThirring.hpp
lines 80-89.
In `@src/shammodels/sph/src/modules/ExternalForces.cpp`:
- Around line 488-500: Update the deferred getters registered in
ExternalForces.cpp at lines 488-500 to capture the external-force data by value
rather than retaining the block-local ext_force reference; apply this to the 1PN
mass, position, and velocity getters. Also update lines 536-541 to value-capture
the Lense–Thirring position and velocity getters, preserving their existing
returned fields.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Repository UI
Review profile: CHILL
Plan: Pro Plus
Run ID: 24f44f40-919e-4158-ace0-8304d1519182
📒 Files selected for processing (6)
src/shammodels/common/include/shammodels/common/ExtForceConfig.hppsrc/shammodels/common/include/shammodels/common/modules/AddForce1PN.hppsrc/shammodels/common/include/shammodels/common/modules/AddForceLenseThirring.hppsrc/shammodels/sph/include/shammodels/sph/SolverConfig.hppsrc/shammodels/sph/src/modules/ExternalForces.cppsrc/shammodels/sph/src/pySPHModel.cpp
| Tscal eps_gr = 1e-10; // small value to avoid division by zero | ||
| Tscal cmass = edges.central_mass.data; | ||
| Tvec cpos = edges.central_pos.data; | ||
| Tvec cvel = edges.central_vel.data; | ||
| Tscal GM = cmass * G; | ||
|
|
||
|
|
||
| sham::distributed_data_kernel_call( | ||
| shamsys::instance::get_compute_scheduler_ptr(), | ||
|
|
||
| sham::DDMultiRef{ | ||
| edges.spans_positions.get_spans(), | ||
| edges.spans_velocities.get_spans() | ||
| }, | ||
|
|
||
| sham::DDMultiRef{ | ||
| edges.spans_accel_ext.get_spans() | ||
| }, | ||
|
|
||
| edges.sizes.indexes, | ||
|
|
||
| [cpos, cvel, GM, c, eps_gr](u32 gid, | ||
| const Tvec *xyz, | ||
| const Tvec *vxyz, | ||
| Tvec *axyz_ext) { | ||
|
|
||
| Tvec r_a = xyz[gid] - cpos; | ||
| Tvec v_a = vxyz[gid] - cvel; | ||
|
|
||
| Tscal r = sycl::length(r_a); | ||
|
|
||
| Tvec r_hat = r_a / (r+eps_gr); | ||
|
|
||
| Tscal v2 = sham::dot(v_a, v_a); | ||
|
|
||
| Tscal vr = sham::dot(v_a, r_hat); | ||
|
|
||
|
|
||
| Tvec acc_1PN = | ||
| -GM / (r * r + eps_gr) | ||
| * | ||
| ( | ||
| ( | ||
| v2 / (c * c) | ||
| - | ||
| 4 * GM / ((r + eps_gr) * c * c) |
There was a problem hiding this comment.
🩺 Stability & Availability | 🟠 Major | ⚡ Quick win
Use one dimensionally consistent radius-clamping strategy. The current kernels either mix length and squared-length epsilon terms or permit division by zero.
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp#L65-L110: clampronce and user_safe,r_safe², andr_safe³consistently.src/shammodels/common/include/shammodels/common/modules/AddForceLenseThirring.hpp#L80-L89: clampabs_rabefore calculating its powers to prevent NaN acceleration at the central position.
📍 Affects 2 files
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp#L65-L110(this comment)src/shammodels/common/include/shammodels/common/modules/AddForceLenseThirring.hpp#L80-L89
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp`
around lines 65 - 110, Use a single dimensionally consistent radius clamp in the
1PN kernel around the acceleration calculation in AddForce1PN: clamp r once into
r_safe, then derive and consistently use r_safe² and r_safe³ for all
radius-dependent denominators, preserving the existing acceleration behavior. In
AddForceLenseThirring, clamp abs_ra before computing its powers so the
central-position case cannot produce NaN acceleration; apply the change in
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp lines
65-110 and
src/shammodels/common/include/shammodels/common/modules/AddForceLenseThirring.hpp
lines 80-89.
| auto set_cmass = register_constant_set<Tscal>(solver_graph, prefix_cmass, [&]() { | ||
| return ext_force->central_mass; | ||
| }); | ||
|
|
||
| auto set_central_pos | ||
| = register_constant_set<Tvec>(solver_graph, prefix_central_pos, [&]() { | ||
| return ext_force->central_pos; | ||
| }); | ||
|
|
||
| auto set_central_vel | ||
| = register_constant_set<Tvec>(solver_graph, prefix_central_vel, [&]() { | ||
| return ext_force->central_vel; | ||
| }); |
There was a problem hiding this comment.
🩺 Stability & Availability | 🔴 Critical | ⚡ Quick win
Capture external-force constants by value. register_constant_set stores each getter for deferred execution, while [&] retains a reference to the block-local ext_force pointer variable after its lifetime ends.
src/shammodels/sph/src/modules/ExternalForces.cpp#L488-L500: value-capture the 1PN mass, position, and velocity.src/shammodels/sph/src/modules/ExternalForces.cpp#L536-L541: value-capture the Lense–Thirring position and velocity.
📍 Affects 1 file
src/shammodels/sph/src/modules/ExternalForces.cpp#L488-L500(this comment)src/shammodels/sph/src/modules/ExternalForces.cpp#L536-L541
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/shammodels/sph/src/modules/ExternalForces.cpp` around lines 488 - 500,
Update the deferred getters registered in ExternalForces.cpp at lines 488-500 to
capture the external-force data by value rather than retaining the block-local
ext_force reference; apply this to the 1PN mass, position, and velocity getters.
Also update lines 536-541 to value-capture the Lense–Thirring position and
velocity getters, preserving their existing returned fields.
There was a problem hiding this comment.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (2)
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp (2)
122-123: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winUse lower_case names for the new TVI C++ identifiers.
The new internal identifiers use
C_1_fluidandC_delta_v, violating the repository naming rule. Rename them toc_1_fluidandc_delta_v; preserve the JSON/Python names if they are external compatibility keys.
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp#L122-L123: rename the newMonofluidTVImembers.src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp#L142-L143: rename the setter parameters and references.src/shammodels/sph/src/pySPHModel.cpp#L273-L274: rename the binding-lambda parameters while retaining publicpy::argnames if required.As per coding guidelines, C++ variables, parameters, and members must use lower_case.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp` around lines 122 - 123, Rename the new MonofluidTVI C++ members in src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp at lines 122-123 to c_1_fluid and c_delta_v, and update the corresponding setter parameters and references at lines 142-143. In src/shammodels/sph/src/pySPHModel.cpp lines 273-274, rename the binding-lambda parameters consistently while preserving any required public py::arg names for JSON/Python compatibility.Source: Coding guidelines
192-197: 🗄️ Data Integrity & Integration | 🟠 Major | ⚡ Quick winKeep TVI defaults optional when loading JSON
mode_from_json()still makes the newmonofluid_tvifields mandatory, so older configs withoutC_1_fluid,C_delta_v,cfl_density_threshold, orensure_s_j_positivitywill fail to load. Read them with the same defaults asset_monofluid_tvi(), or version the format.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp` around lines 192 - 197, The mode_from_json() loader currently requires the new monofluid_tvi fields; make C_1_fluid, C_delta_v, cfl_density_threshold, and ensure_s_j_positivity optional by applying the same defaults defined by set_monofluid_tvi(). Preserve loading of older JSON configurations while retaining explicitly provided values.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Outside diff comments:
In `@src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp`:
- Around line 122-123: Rename the new MonofluidTVI C++ members in
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp at lines 122-123 to
c_1_fluid and c_delta_v, and update the corresponding setter parameters and
references at lines 142-143. In src/shammodels/sph/src/pySPHModel.cpp lines
273-274, rename the binding-lambda parameters consistently while preserving any
required public py::arg names for JSON/Python compatibility.
- Around line 192-197: The mode_from_json() loader currently requires the new
monofluid_tvi fields; make C_1_fluid, C_delta_v, cfl_density_threshold, and
ensure_s_j_positivity optional by applying the same defaults defined by
set_monofluid_tvi(). Preserve loading of older JSON configurations while
retaining explicitly provided values.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Repository UI
Review profile: ASSERTIVE
Plan: Pro Plus
Run ID: 8f53b1d0-0ca9-4250-ae76-4ba99b67ef28
📒 Files selected for processing (2)
src/shammodels/sph/include/shammodels/sph/SolverConfig.hppsrc/shammodels/sph/src/pySPHModel.cpp
There was a problem hiding this comment.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp (1)
102-117: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick winInclude the Newtonian term in this standalone 1PN force.
PN_1PNwiring registers only this node, but this expression adds only the (O(c^{-2})) correction. It therefore omits the leading (-GM\hat r/r^2) acceleration entirely.Proposed fix
( - v2 / (c * c) + Tscal{1} + v2 / (c * c) - 4 * GM * inv_r / (c * c) )🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp` around lines 102 - 117, Update the standalone 1PN acceleration expression in PN_1PN to include the leading Newtonian -GM * inv_r2 * r_hat term in addition to the existing O(c^-2) correction, so the registered node returns the complete acceleration.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Outside diff comments:
In `@src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp`:
- Around line 102-117: Update the standalone 1PN acceleration expression in
PN_1PN to include the leading Newtonian -GM * inv_r2 * r_hat term in addition to
the existing O(c^-2) correction, so the registered node returns the complete
acceleration.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Repository UI
Review profile: ASSERTIVE
Plan: Pro Plus
Run ID: 737fdea0-31ba-470c-a8f6-d237e1760a45
📒 Files selected for processing (1)
src/shammodels/common/include/shammodels/common/modules/AddForce1PN.hpp
3 participants
-1PN for Sink-disk interaction added
-Lense thirring for Sink-Disk interaction modified