Skip to content

[SPH] Add PN (1, 1.5, 2, 2.5) terms for sink binaries - #1904

Open
AugustinDart wants to merge 23 commits into
Shamrock-code:mainfrom
AugustinDart:SPIN/RR
Open

[SPH] Add PN (1, 1.5, 2, 2.5) terms for sink binaries#1904
AugustinDart wants to merge 23 commits into
Shamrock-code:mainfrom
AugustinDart:SPIN/RR

Conversation

@AugustinDart

@AugustinDart AugustinDart commented Jul 6, 2026

Copy link
Copy Markdown
Contributor

1PN added for sink-sink interaction (orbital precession)
1.5 PN added for sink-sink interaction (Lense thirring frame dragging effect)
2PN added for Sink-sink interaction (we only added the Spin-spin interaction and not the term 2PN of the orbital precession)
2.5 PN added (Radiation Reaction or emission of GW)

We added the equation of Spins evolution for the binary in SinkParticleUpdate.cpp

Run_sph_binary_orbit.py is also modified in order to see the different effects (with the visualisation of orbital elements --> a,e, i and w). Possible to observe the evolution of spins too.

@github-actions

github-actions Bot commented Jul 6, 2026

Copy link
Copy Markdown
Contributor

Thanks @AugustinDart for opening this PR!

You can do multiple things directly here:
1 - Comment pre-commit.ci run to run pre-commit checks.
2 - Comment pre-commit.ci autofix to apply fixes.
3 - Add label autofix.ci to fix authorship & pre-commit for every commit made.
4 - Add label light-ci to only trigger a reduced & faster version of the CI (need the full one before merge).
5 - Add label trigger-ci to create an empty commit to trigger the CI.

Once the workflow completes a message will appear displaying informations related to the run.

Also the PR gets automatically reviewed by gemini, you can:
1 - Comment /gemini review to trigger a review
2 - Comment /gemini summary for a summary
3 - Tag it using @gemini-code-assist either in the PR or in review comments on files

@AugustinDart
AugustinDart marked this pull request as draft July 6, 2026 15:30

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Code Review

This pull request introduces Post-Newtonian (PN) corrections (including orbit precession, spin-orbit, spin-spin, and radiation reaction terms) and spin precession evolution for sink particles, along with new test and visualization scripts. The review feedback highlights a critical bug in compute_ext_forces where the force accumulation is incorrectly scaled inside the loop, as well as a performance issue caused by printing configuration logs to std::cout on every timestep.

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.

Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
@tdavidcl tdavidcl changed the title Spin/rr [SPH] Add PN (1, 1.5, 2, 2.5) terms for sink binaries Jul 6, 2026
@y-lapeyre
y-lapeyre self-requested a review July 6, 2026 15:34
@y-lapeyre y-lapeyre added the draft label Jul 6, 2026
@AugustinDart
AugustinDart marked this pull request as ready for review July 7, 2026 11:38
@github-actions

github-actions Bot commented Jul 9, 2026

Copy link
Copy Markdown
Contributor

Workflow report

workflow report corresponding to commit 0e0697a
Commiter email is 66853113+pre-commit-ci[bot]@users.noreply.github.com
You are using github private e-mail. This prevent proper tracing of who contributed what, please disable it (see Keep my email addresses private).

Pre-commit check report

Some failures were detected in base source checks checks.
Check the On PR / Linting / Base source checks (pull_request) job in the tests for more detailed output

Suggested changes

Detailed changes :

Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
Comment thread orbit.png Outdated
Comment thread orbit.png Outdated
Comment thread spins.png Outdated
Comment thread vitesse.png Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
@y-lapeyre

Copy link
Copy Markdown
Collaborator

Please translate all comments to English. Non ascii characters ("é"'s) will make the pre-commit test fail :)

@coderabbitai

coderabbitai Bot commented Jul 16, 2026

Copy link
Copy Markdown

Review Change Stack

Note

Reviews paused

It looks like this branch is under active development. To avoid overwhelming you with review comments due to an influx of new commits, CodeRabbit has automatically paused this review. You can configure this behavior by changing the reviews.auto_review.auto_pause_after_reviewed_commits setting.

Use the following commands to manage reviews:

  • @coderabbitai resume to resume automatic reviews.
  • @coderabbitai review to trigger a single review.

Use the checkboxes below for quick actions:

  • ▶️ Resume reviews
  • 🔍 Trigger review
📝 Walkthrough

Walkthrough

Adds configurable post-Newtonian sink forces and spin updates, persists PN solver flags, preserves sink angular momentum, and introduces three runnable binary-orbit examples with SPH evolution, orbital diagnostics, plotting, and disk rendering.

Changes

Post-Newtonian sink dynamics

Layer / File(s) Summary
Sink state and PN configuration contracts
src/shammodels/sph/include/shammodels/sph/Model.hpp, src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp, src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp, src/shammodels/sph/src/pySPHModel.cpp
Sink insertion accepts angular momentum, solver configuration stores and serializes four PN flags, Python exposes their setters, and sink force APIs accept a timestep.
Timestep-aware PN force pipeline
src/shammodels/sph/src/Solver.cpp, src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
The solver forwards dt; sink updates accumulate Newtonian and enabled PN pairwise terms and update sink spins.
Binary sink setup and orbital analysis
examples/physics/run_binary_sink_only.py
A runnable binary-sink example initializes rotated orbits and spins, records evolution, computes orbital elements, plots diagnostics, and exports CSV data.
Binary sink disk workflow
examples/physics/run_binary_sink_disk.py
A runnable disk example builds and recenters an SPH disk with binary sinks, evolves the system, records snapshots, and renders orbit and density animations.
Sink integration test workflow
examples/physics/run_sink_integ_test_problems.py
A standalone example builds a circular sink binary, advances it with configurable stepping, records snapshots, and plots trajectories.

Estimated code review effort: 4 (Complex) | ~60 minutes

Sequence Diagram(s)

sequenceDiagram
  participant BinaryExample
  participant Solver
  participant SinkParticlesUpdate
  participant SinkPairInteraction
  participant OrbitalAnalysis
  BinaryExample->>Solver: evolve_once_override_time(dt)
  Solver->>SinkParticlesUpdate: compute_ext_forces(dt)
  SinkParticlesUpdate->>SinkPairInteraction: accumulate Newtonian and enabled PN terms
  SinkPairInteraction-->>SinkParticlesUpdate: sink accelerations
  SinkParticlesUpdate->>SinkParticlesUpdate: update_sink_spins(dt)
  Solver-->>BinaryExample: updated sink state
  BinaryExample->>OrbitalAnalysis: compute orbital and spin diagnostics
Loading
🚥 Pre-merge checks | ✅ 4 | ❌ 1

❌ Failed checks (1 warning)

Check name Status Explanation Resolution
Docstring Coverage ⚠️ Warning Docstring coverage is 7.25% which is insufficient. The required threshold is 80.00%. Write docstrings for the functions missing them to satisfy the coverage threshold.
✅ Passed checks (4 passed)
Check name Status Explanation
Title check ✅ Passed The title accurately summarizes the main change: adding PN terms for sink binaries.
Description check ✅ Passed The description matches the changes by describing the added PN terms, spin evolution, and orbit visualization.
Linked Issues check ✅ Passed Check skipped because no linked issues were found for this pull request.
Out of Scope Changes check ✅ Passed Check skipped because no linked issues were found for this pull request.
✨ Finishing Touches
🧪 Generate unit tests (beta)
  • Create PR with unit tests

Comment @coderabbitai help to get the list of available commands.

@y-lapeyre

Copy link
Copy Markdown
Collaborator

/gemini review

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Code Review

This pull request introduces Post-Newtonian (PN) developments for binary orbit functions and sink particles, adding configuration flags for orbital precession, spin-orbit, spin-spin, and radiation reaction terms, along with several example scripts. The review feedback highlights critical issues in the C++ implementation of these PN terms, including integer division bugs in the radiation reaction term, potential division-by-zero errors when computing unit vectors, a missing factor of three in the spin-spin term, and an update-in-place bug in spin precession that violates angular momentum conservation. Additionally, the feedback recommends replacing hardcoded absolute local file paths in the Python examples with relative paths.

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.

Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment thread src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp Outdated
Comment on lines +515 to 561
void shammodels::sph::modules::SinkParticlesUpdate<Tvec, SPHKernel>::update_sink_spins(Tscal dt) {

//Definition of the constants G and c for the calculations of spin precession (the same as in the compute_ext_forces function)
Tscal G = solver_config.get_constant_G(); //G=4*pi*2
Tscal c = solver_config.get_constant_c(); //c= 63 241.077 AU/year

if (storage.sinks.is_empty()) {
return;
}

std::vector<Sink> &sink_parts = storage.sinks.get();


Tscal epsilon_spin = 1e-9;

for (Sink &s1 : sink_parts) {
Tvec dS = {};

for (Sink &s2 : sink_parts) {
if (&s1 == &s2) {
continue;
}
Tvec rij = s1.pos - s2.pos;
Tvec vij = s1.velocity - s2.velocity;
Tscal m1=s1.mass;
Tscal m2=s2.mass;
Tscal M = m1 + m2;
Tscal nu = m1 * m2 / M;
Tvec L = nu * sycl::cross(rij, vij);


Tscal rij_scal = sycl::length(rij) + epsilon_spin;
Tvec nij = rij / rij_scal;
Tvec S1 = s1.angular_momentum;
Tvec S2 = s2.angular_momentum;
Tscal prefactor = G / (c * c * rij_scal * rij_scal * rij_scal);
// Simple spin precession structure.


Tvec Omega_prec = prefactor * ((2 + 3 * m2 / (2 * m1)) * L - S2 + 3*sycl::dot(nij, S2) * nij) ;

dS += sycl::cross(Omega_prec, S1);
}
s1.ext_acceleration = sum;

s1.angular_momentum += dS * dt;
}
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

high

There are two major issues in update_sink_spins:

  1. Update-in-place coupled integration bug: s1.angular_momentum is updated immediately inside the loop. In the next iteration, when s1 becomes s2, the calculation uses the already updated spin at $t + dt$ instead of the state at $t$. This breaks symmetry and violates conservation of total angular momentum. Spin derivatives should be accumulated first and applied in a separate loop.
  2. Numerical instability: rij_scal is defined as sycl::length(rij) + epsilon_spin. If length(rij) is 0, rij_scal is 1e-9, and rij_scal^3 is 1e-27. Dividing by 1e-27 in prefactor causes an extreme numerical spike and overflow. It is much safer to add epsilon_spin directly to the denominator of prefactor.
void shammodels::sph::modules::SinkParticlesUpdate<Tvec, SPHKernel>::update_sink_spins(Tscal dt) {

//Definition of the constants G and c for the calculations of spin precession (the same as in the compute_ext_forces function)
Tscal G = solver_config.get_constant_G();       //G=4*pi*2
Tscal c = solver_config.get_constant_c();     //c= 63 241.077 AU/year

    if (storage.sinks.is_empty()) {
        return;
    }

    std::vector<Sink> &sink_parts = storage.sinks.get();

    Tscal epsilon_spin = 1e-9;
    std::vector<Tvec> dS_list(sink_parts.size(), Tvec{});

    for (size_t i = 0; i < sink_parts.size(); ++i) {
        Sink &s1 = sink_parts[i];
        Tvec dS = {};

        for (size_t j = 0; j < sink_parts.size(); ++j) {
            Sink &s2 = sink_parts[j];
            if (i == j) {
                continue;
            }
            Tvec rij = s1.pos - s2.pos;
            Tvec vij = s1.velocity - s2.velocity;
            Tscal m1=s1.mass;
            Tscal m2=s2.mass;
            Tscal M = m1 + m2;
            Tscal nu = m1 * m2 / M;
            Tvec L = nu * sycl::cross(rij, vij);

            Tscal rij_scal = sycl::length(rij);
            Tvec nij = rij / (rij_scal + epsilon_spin);
            Tvec S1 = s1.angular_momentum;
            Tvec S2 = s2.angular_momentum;
            Tscal prefactor = G / (c * c * (rij_scal * rij_scal * rij_scal + epsilon_spin));
            // Simple spin precession structure.
            

            Tvec Omega_prec = prefactor * ((2 + 3 * m2 / (2 * m1)) * L - S2 + 3*sycl::dot(nij, S2) * nij) ;
            
            dS += sycl::cross(Omega_prec, S1);
        }
        dS_list[i] = dS;
    }

    for (size_t i = 0; i < sink_parts.size(); ++i) {
        sink_parts[i].angular_momentum += dS_list[i] * dt;
    }
}
References
  1. If a test for a function's property (like symmetry) is failing, investigate the root cause in the test or function itself instead of altering the test to make it pass trivially.

Comment thread examples/physics/run_binary_sink_only.py Outdated
Comment thread examples/physics/run_binary_sink_only.py Outdated

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Actionable comments posted: 19

🧹 Nitpick comments (4)
examples/physics/run_sink_integ_test_problems.py (1)

141-141: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Resolve the reported unused unpacking.

nstep and ndim are unused, triggering Ruff RUF059. Replace them with _nstep/_ndim or unpack only the needed dimension.

🤖 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 `@examples/physics/run_sink_integ_test_problems.py` at line 141, Update the
sinks_positions shape unpacking to avoid unused variables flagged by Ruff
RUF059: retain nsink for downstream use and replace unused nstep and ndim with
underscore-prefixed names or unpack only the required dimension.

Source: Linters/SAST tools

src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp (1)

401-404: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low value

New physics-notation variables use uppercase naming.

OP, SO, SS, RR, G, c, M, S1, S2, S, Delta are all new local variables that don't follow lower_case. This matches existing convention elsewhere in the codebase (e.g. get_constant_G()), so it's a minor/pre-existing-style issue rather than a new regression, but it's still a guideline violation in the changed lines.

As per coding guidelines, "use lower_case for functions, variables, parameters, and members" for **/*.{cpp,cc,cxx,h,hpp,hh}.

Also applies to: 419-420, 455-459

🤖 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/SinkParticlesUpdate.cpp` around lines 401 -
404, Rename the newly introduced physics-notation locals in the affected logic
to lower_case identifiers, including OP, SO, SS, RR, G, c, M, S1, S2, S, and
Delta, and update all references consistently. Preserve their existing meanings
and behavior.

Source: Coding guidelines

src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp (1)

50-51: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low value

compute_ext_forces(dt) signature update looks correct; consider update_sink_spins visibility.

update_sink_spins is only called internally from compute_ext_forces (see SinkParticlesUpdate.cpp line 511) — no external caller is visible in Solver.cpp. Consider making it private unless it's intended to be unit-testable/callable independently.

🤖 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/modules/SinkParticlesUpdate.hpp`
around lines 50 - 51, Change update_sink_spins(Tscal dt) in the
SinkParticlesUpdate class from public to private, since compute_ext_forces is
its only caller. Keep compute_ext_forces publicly accessible and preserve the
existing internal call.
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp (1)

421-430: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Document the PN flag acronyms (OP/SO/SS/RR) at both the C++ definition and Python binding. Both sites expose the same undocumented acronyms with no explanation of what they control.

  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp#L421-L430: add @brief doc comments on compute_OP/compute_SO/compute_SS/compute_RR and their setters explaining Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and Radiation Reaction.
  • src/shammodels/sph/src/pySPHModel.cpp#L74-L77: add a short docstring/kwarg description to set_compute_OP/set_compute_SO/set_compute_SS/set_compute_RR, similar to the R"pbdoc(...)" block used for make_generator_disc_mc in the same file.
🤖 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 421
- 430, Document the PN flag acronyms at both affected sites: in
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp lines 421-430, add
`@brief` comments to compute_OP, compute_SO, compute_SS, compute_RR and their
setters identifying Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and
Radiation Reaction; in src/shammodels/sph/src/pySPHModel.cpp lines 74-77, add
matching short kwarg/docstring descriptions for set_compute_OP, set_compute_SO,
set_compute_SS, and set_compute_RR using the existing R"pbdoc(...)"
documentation style.
🤖 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 `@examples/physics/run_binary_sink_disk.py`:
- Around line 526-551: Move the ani.save call in the animation setup flow to
after the ax.set_title, axis-label, limit, aspect, fig.colorbar, and
plt.tight_layout configuration calls, so the generated GIF includes all figure
decorations.
- Around line 366-385: Update the render block in the simulation loop so it does
not append every full-resolution rho array to render_frames; stream each
generated frame directly to the animation writer, or otherwise reduce retained
frame data by using a larger render_stride or lower nx/ny resolution while
preserving rendering behavior.
- Around line 253-296: Move the momentum and barycenter correction,
recalculation, logging, and assertions into the existing generate_disk
conditional. When generate_disk is false, skip the entire correction block so
sink-only mode performs no disk analysis.
- Around line 579-591: Update the render_disk_and_orbit call in the
run_binary_orbit_PN flow to pass the same 15 AU spatial extent used by
render_ext (rout * 1.5). Preserve the existing frame rendering while ensuring
density images and sink coordinates share the render extent.

In `@examples/physics/run_sink_integ_test_problems.py`:
- Around line 196-202: Update the get_binary_rotated call in the circular-orbit
setup to pass zero eccentricity via e=0.0, while leaving the other orbital
parameters unchanged.
- Around line 66-69: Before the sink-creation loop in the relevant setup,
validate that positions, velocities, masses, and accretion_radii have equal
lengths and reject mismatches instead of allowing zip() to truncate. Also update
the “Circular orbit” setup’s eccentricity argument from e=0.2 to e=0 so its
initial condition is circular.

In `@examples/physics/run_sph_binary_sink.py`:
- Line 660: Update the default filename in save_orbital_elements to use a
portable relative ASCII-only path, and apply the same correction to the
corresponding filename default at the other reported occurrence. Remove the
author-specific absolute directory and accented characters while preserving the
existing output behavior.
- Around line 31-36: Translate all French comments, labels, and other
user-facing text in run_sph_binary_sink.py—including the referenced ranges—into
clear English, and replace accented or other non-ASCII characters with ASCII
equivalents. Preserve the existing code behavior and numeric values while
ensuring the entire affected text is ASCII-only.
- Around line 39-43: Update the physical setup comments near the spin parameters
and periapsis configuration: describe the spins as tilted by theta = pi/3 rather
than aligned with orbital angular momentum, and correct the periapsis value to A
* (1 - E) = 90 AU. Change comments only; preserve the existing parameter values
and calculations.
- Around line 77-92: Both binary initializer functions expose nu but always
initialize periapsis conditions; update binary_initial_conditions in
examples/physics/run_sph_binary_sink.py (lines 77-92) and the corresponding
initializer in examples/physics/run_binary_sink_disk.py (lines 135-150) to use
anomaly-dependent position and velocity calculations for nonzero nu, or remove
nu from both APIs and all callers consistently.
- Around line 39-47: Keep only dimensionless spin parameters and the shared spin
axis in the module-level setup, avoiding mass-dependent spin vectors. In
examples/physics/run_sph_binary_sink.py lines 167-180, update the parameterized
builder to derive spin vectors from its supplied m1 and m2 values or accept
explicit vectors; in lines 623-628, pass or accept m1 and m2 rather than relying
on main-block globals. Apply the same changes in
examples/physics/run_binary_sink_disk.py lines 97-105 and 229-242 so each
builder derives or receives vectors scaled for its supplied masses.
- Around line 609-610: Update the omega processing near the periapsis-angle
plotting logic to unwrap the angular values rather than converting between
degrees and radians. Use a continuous-angle unwrapping operation with the
correct unit so jumps at 360 degrees are removed while preserving the plotted
angle values.
- Around line 570-574: Update the eccentricity-vector calculation in the
relevant function to return the raw e_vec from the cross-product expression
without dividing by its norm or replacing zero-magnitude vectors. Preserve the
existing component values so e_x, e_y, and e_z retain the actual
eccentricity-vector magnitude.
- Around line 738-750: Move the orbital plotting, save_orbital_elements,
model.get_sinks, and plot_spins calls into the existing main execution guard
that initializes snapshots, m1, m2, and model. Keep their current order and
arguments so importing the module performs no post-processing or model access.

In `@src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp`:
- Around line 407-411: Move the five Post-Newtonian status logs from the
per-timestep force-evaluation path in compute_ext_forces to a one-time
initialization or SolverConfig::check_config() location. Preserve the reported
configuration values while ensuring they are emitted only once per solver run.
- Around line 497-505: Update the RR term calculation in the visible term4
expression so every fractional coefficient uses floating-point arithmetic,
especially the 8/5 prefactor and both 2/3 gravitational sub-terms. Preserve the
existing formula while ensuring these coefficients evaluate as 1.6 and
approximately 0.667 rather than integer divisions.
- Around line 447-450: Guard the unit-vector calculation in the sink-particle
update around rij and nij by adding the established epsilon_spin to rij_scal
before dividing. Preserve the existing rij distance calculation and ensure exact
position coincidence cannot produce 0/0 or propagate NaNs into PN-term
ext_acceleration.
- Around line 399-404: Guard the PN-term setup in the acceleration calculation
around OP, SO, SS, and RR so configurations with more than two sinks cannot
apply these binary-only corrections. When sink_parts.size() exceeds two and any
PN term is enabled, emit an appropriate warning and disable or otherwise prevent
those PN contributions while preserving Newtonian pairwise acceleration.
- Line 258: Remove the update_sink_spins(dt) call from compute_ext_forces so
that method remains force-only. Invoke update_sink_spins(dt) exactly once per
timestep from the solver’s timestep flow, outside both predictor_step and later
force re-evaluation calls, while preserving the existing dt argument and
accumulation behavior.

---

Nitpick comments:
In `@examples/physics/run_sink_integ_test_problems.py`:
- Line 141: Update the sinks_positions shape unpacking to avoid unused variables
flagged by Ruff RUF059: retain nsink for downstream use and replace unused nstep
and ndim with underscore-prefixed names or unpack only the required dimension.

In `@src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp`:
- Around line 50-51: Change update_sink_spins(Tscal dt) in the
SinkParticlesUpdate class from public to private, since compute_ext_forces is
its only caller. Keep compute_ext_forces publicly accessible and preserve the
existing internal call.

In `@src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp`:
- Around line 421-430: Document the PN flag acronyms at both affected sites: in
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp lines 421-430, add
`@brief` comments to compute_OP, compute_SO, compute_SS, compute_RR and their
setters identifying Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and
Radiation Reaction; in src/shammodels/sph/src/pySPHModel.cpp lines 74-77, add
matching short kwarg/docstring descriptions for set_compute_OP, set_compute_SO,
set_compute_SS, and set_compute_RR using the existing R"pbdoc(...)"
documentation style.

In `@src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp`:
- Around line 401-404: Rename the newly introduced physics-notation locals in
the affected logic to lower_case identifiers, including OP, SO, SS, RR, G, c, M,
S1, S2, S, and Delta, and update all references consistently. Preserve their
existing meanings and behavior.
🪄 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: cfdac8b8-d3a8-4892-ad47-2177cf9d0296

📥 Commits

Reviewing files that changed from the base of the PR and between d12b8c7 and b6bfe2f.

📒 Files selected for processing (9)
  • examples/physics/run_binary_sink_disk.py
  • examples/physics/run_sink_integ_test_problems.py
  • examples/physics/run_sph_binary_sink.py
  • src/shammodels/sph/include/shammodels/sph/Model.hpp
  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
  • src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp
  • src/shammodels/sph/src/Solver.cpp
  • src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
  • src/shammodels/sph/src/pySPHModel.cpp

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Caution

Inline review comments failed to post. This is likely due to GitHub's internal server error or limits when posting large numbers of comments. If you are seeing this consistently it is likely a permissions issue. Please check "Moderation" -> "Code review limits" under your organization settings.

Actionable comments posted: 19

🧹 Nitpick comments (4)
examples/physics/run_sink_integ_test_problems.py (1)

141-141: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Resolve the reported unused unpacking.

nstep and ndim are unused, triggering Ruff RUF059. Replace them with _nstep/_ndim or unpack only the needed dimension.

🤖 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 `@examples/physics/run_sink_integ_test_problems.py` at line 141, Update the
sinks_positions shape unpacking to avoid unused variables flagged by Ruff
RUF059: retain nsink for downstream use and replace unused nstep and ndim with
underscore-prefixed names or unpack only the required dimension.

Source: Linters/SAST tools

src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp (1)

401-404: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low value

New physics-notation variables use uppercase naming.

OP, SO, SS, RR, G, c, M, S1, S2, S, Delta are all new local variables that don't follow lower_case. This matches existing convention elsewhere in the codebase (e.g. get_constant_G()), so it's a minor/pre-existing-style issue rather than a new regression, but it's still a guideline violation in the changed lines.

As per coding guidelines, "use lower_case for functions, variables, parameters, and members" for **/*.{cpp,cc,cxx,h,hpp,hh}.

Also applies to: 419-420, 455-459

🤖 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/SinkParticlesUpdate.cpp` around lines 401 -
404, Rename the newly introduced physics-notation locals in the affected logic
to lower_case identifiers, including OP, SO, SS, RR, G, c, M, S1, S2, S, and
Delta, and update all references consistently. Preserve their existing meanings
and behavior.

Source: Coding guidelines

src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp (1)

50-51: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low value

compute_ext_forces(dt) signature update looks correct; consider update_sink_spins visibility.

update_sink_spins is only called internally from compute_ext_forces (see SinkParticlesUpdate.cpp line 511) — no external caller is visible in Solver.cpp. Consider making it private unless it's intended to be unit-testable/callable independently.

🤖 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/modules/SinkParticlesUpdate.hpp`
around lines 50 - 51, Change update_sink_spins(Tscal dt) in the
SinkParticlesUpdate class from public to private, since compute_ext_forces is
its only caller. Keep compute_ext_forces publicly accessible and preserve the
existing internal call.
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp (1)

421-430: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Document the PN flag acronyms (OP/SO/SS/RR) at both the C++ definition and Python binding. Both sites expose the same undocumented acronyms with no explanation of what they control.

  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp#L421-L430: add @brief doc comments on compute_OP/compute_SO/compute_SS/compute_RR and their setters explaining Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and Radiation Reaction.
  • src/shammodels/sph/src/pySPHModel.cpp#L74-L77: add a short docstring/kwarg description to set_compute_OP/set_compute_SO/set_compute_SS/set_compute_RR, similar to the R"pbdoc(...)" block used for make_generator_disc_mc in the same file.
🤖 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 421
- 430, Document the PN flag acronyms at both affected sites: in
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp lines 421-430, add
`@brief` comments to compute_OP, compute_SO, compute_SS, compute_RR and their
setters identifying Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and
Radiation Reaction; in src/shammodels/sph/src/pySPHModel.cpp lines 74-77, add
matching short kwarg/docstring descriptions for set_compute_OP, set_compute_SO,
set_compute_SS, and set_compute_RR using the existing R"pbdoc(...)"
documentation style.
🤖 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 `@examples/physics/run_binary_sink_disk.py`:
- Around line 526-551: Move the ani.save call in the animation setup flow to
after the ax.set_title, axis-label, limit, aspect, fig.colorbar, and
plt.tight_layout configuration calls, so the generated GIF includes all figure
decorations.
- Around line 366-385: Update the render block in the simulation loop so it does
not append every full-resolution rho array to render_frames; stream each
generated frame directly to the animation writer, or otherwise reduce retained
frame data by using a larger render_stride or lower nx/ny resolution while
preserving rendering behavior.
- Around line 253-296: Move the momentum and barycenter correction,
recalculation, logging, and assertions into the existing generate_disk
conditional. When generate_disk is false, skip the entire correction block so
sink-only mode performs no disk analysis.
- Around line 579-591: Update the render_disk_and_orbit call in the
run_binary_orbit_PN flow to pass the same 15 AU spatial extent used by
render_ext (rout * 1.5). Preserve the existing frame rendering while ensuring
density images and sink coordinates share the render extent.

In `@examples/physics/run_sink_integ_test_problems.py`:
- Around line 196-202: Update the get_binary_rotated call in the circular-orbit
setup to pass zero eccentricity via e=0.0, while leaving the other orbital
parameters unchanged.
- Around line 66-69: Before the sink-creation loop in the relevant setup,
validate that positions, velocities, masses, and accretion_radii have equal
lengths and reject mismatches instead of allowing zip() to truncate. Also update
the “Circular orbit” setup’s eccentricity argument from e=0.2 to e=0 so its
initial condition is circular.

In `@examples/physics/run_sph_binary_sink.py`:
- Line 660: Update the default filename in save_orbital_elements to use a
portable relative ASCII-only path, and apply the same correction to the
corresponding filename default at the other reported occurrence. Remove the
author-specific absolute directory and accented characters while preserving the
existing output behavior.
- Around line 31-36: Translate all French comments, labels, and other
user-facing text in run_sph_binary_sink.py—including the referenced ranges—into
clear English, and replace accented or other non-ASCII characters with ASCII
equivalents. Preserve the existing code behavior and numeric values while
ensuring the entire affected text is ASCII-only.
- Around line 39-43: Update the physical setup comments near the spin parameters
and periapsis configuration: describe the spins as tilted by theta = pi/3 rather
than aligned with orbital angular momentum, and correct the periapsis value to A
* (1 - E) = 90 AU. Change comments only; preserve the existing parameter values
and calculations.
- Around line 77-92: Both binary initializer functions expose nu but always
initialize periapsis conditions; update binary_initial_conditions in
examples/physics/run_sph_binary_sink.py (lines 77-92) and the corresponding
initializer in examples/physics/run_binary_sink_disk.py (lines 135-150) to use
anomaly-dependent position and velocity calculations for nonzero nu, or remove
nu from both APIs and all callers consistently.
- Around line 39-47: Keep only dimensionless spin parameters and the shared spin
axis in the module-level setup, avoiding mass-dependent spin vectors. In
examples/physics/run_sph_binary_sink.py lines 167-180, update the parameterized
builder to derive spin vectors from its supplied m1 and m2 values or accept
explicit vectors; in lines 623-628, pass or accept m1 and m2 rather than relying
on main-block globals. Apply the same changes in
examples/physics/run_binary_sink_disk.py lines 97-105 and 229-242 so each
builder derives or receives vectors scaled for its supplied masses.
- Around line 609-610: Update the omega processing near the periapsis-angle
plotting logic to unwrap the angular values rather than converting between
degrees and radians. Use a continuous-angle unwrapping operation with the
correct unit so jumps at 360 degrees are removed while preserving the plotted
angle values.
- Around line 570-574: Update the eccentricity-vector calculation in the
relevant function to return the raw e_vec from the cross-product expression
without dividing by its norm or replacing zero-magnitude vectors. Preserve the
existing component values so e_x, e_y, and e_z retain the actual
eccentricity-vector magnitude.
- Around line 738-750: Move the orbital plotting, save_orbital_elements,
model.get_sinks, and plot_spins calls into the existing main execution guard
that initializes snapshots, m1, m2, and model. Keep their current order and
arguments so importing the module performs no post-processing or model access.

In `@src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp`:
- Around line 407-411: Move the five Post-Newtonian status logs from the
per-timestep force-evaluation path in compute_ext_forces to a one-time
initialization or SolverConfig::check_config() location. Preserve the reported
configuration values while ensuring they are emitted only once per solver run.
- Around line 497-505: Update the RR term calculation in the visible term4
expression so every fractional coefficient uses floating-point arithmetic,
especially the 8/5 prefactor and both 2/3 gravitational sub-terms. Preserve the
existing formula while ensuring these coefficients evaluate as 1.6 and
approximately 0.667 rather than integer divisions.
- Around line 447-450: Guard the unit-vector calculation in the sink-particle
update around rij and nij by adding the established epsilon_spin to rij_scal
before dividing. Preserve the existing rij distance calculation and ensure exact
position coincidence cannot produce 0/0 or propagate NaNs into PN-term
ext_acceleration.
- Around line 399-404: Guard the PN-term setup in the acceleration calculation
around OP, SO, SS, and RR so configurations with more than two sinks cannot
apply these binary-only corrections. When sink_parts.size() exceeds two and any
PN term is enabled, emit an appropriate warning and disable or otherwise prevent
those PN contributions while preserving Newtonian pairwise acceleration.
- Line 258: Remove the update_sink_spins(dt) call from compute_ext_forces so
that method remains force-only. Invoke update_sink_spins(dt) exactly once per
timestep from the solver’s timestep flow, outside both predictor_step and later
force re-evaluation calls, while preserving the existing dt argument and
accumulation behavior.

---

Nitpick comments:
In `@examples/physics/run_sink_integ_test_problems.py`:
- Line 141: Update the sinks_positions shape unpacking to avoid unused variables
flagged by Ruff RUF059: retain nsink for downstream use and replace unused nstep
and ndim with underscore-prefixed names or unpack only the required dimension.

In `@src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp`:
- Around line 50-51: Change update_sink_spins(Tscal dt) in the
SinkParticlesUpdate class from public to private, since compute_ext_forces is
its only caller. Keep compute_ext_forces publicly accessible and preserve the
existing internal call.

In `@src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp`:
- Around line 421-430: Document the PN flag acronyms at both affected sites: in
src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp lines 421-430, add
`@brief` comments to compute_OP, compute_SO, compute_SS, compute_RR and their
setters identifying Orbital Precession (1PN), Spin-Orbit, Spin-Spin, and
Radiation Reaction; in src/shammodels/sph/src/pySPHModel.cpp lines 74-77, add
matching short kwarg/docstring descriptions for set_compute_OP, set_compute_SO,
set_compute_SS, and set_compute_RR using the existing R"pbdoc(...)"
documentation style.

In `@src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp`:
- Around line 401-404: Rename the newly introduced physics-notation locals in
the affected logic to lower_case identifiers, including OP, SO, SS, RR, G, c, M,
S1, S2, S, and Delta, and update all references consistently. Preserve their
existing meanings and behavior.
🪄 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: cfdac8b8-d3a8-4892-ad47-2177cf9d0296

📥 Commits

Reviewing files that changed from the base of the PR and between d12b8c7 and b6bfe2f.

📒 Files selected for processing (9)
  • examples/physics/run_binary_sink_disk.py
  • examples/physics/run_sink_integ_test_problems.py
  • examples/physics/run_sph_binary_sink.py
  • src/shammodels/sph/include/shammodels/sph/Model.hpp
  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
  • src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp
  • src/shammodels/sph/src/Solver.cpp
  • src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
  • src/shammodels/sph/src/pySPHModel.cpp
🛑 Comments failed to post (19)
examples/physics/run_binary_sink_disk.py (4)

253-296: 🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

Skip disk corrections when generate_disk=False.

The function advertises a sink-only mode but still computes and asserts momentum and barycenter values for an empty disk. Keep the correction block inside the generate_disk branch.

🧰 Tools
🪛 Ruff (0.15.21)

[warning] 287-287: Unpacked variable disc_mass_value is never used

Prefix it with an underscore or any other dummy variable pattern

(RUF059)

🤖 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 `@examples/physics/run_binary_sink_disk.py` around lines 253 - 296, Move the
momentum and barycenter correction, recalculation, logging, and assertions into
the existing generate_disk conditional. When generate_disk is false, skip the
entire correction block so sink-only mode performs no disk analysis.

366-385: 🚀 Performance & Scalability | 🟠 Major | 🏗️ Heavy lift

Avoid retaining every full-resolution density frame in memory.

Each render appends a dense rho array. With the current caller's 1024-by-1024 resolution and roughly 200 steps, this alone retains about 1.6 GiB. Stream frames to the animation writer or enforce a larger stride/lower resolution.

🤖 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 `@examples/physics/run_binary_sink_disk.py` around lines 366 - 385, Update the
render block in the simulation loop so it does not append every full-resolution
rho array to render_frames; stream each generated frame directly to the
animation writer, or otherwise reduce retained frame data by using a larger
render_stride or lower nx/ny resolution while preserving rendering behavior.

526-551: 🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win

Configure the figure before saving the animation.

The GIF is saved before its title, labels, limits, aspect ratio, and colorbar are added, so those decorations are absent from the generated file. Move ani.save(...) after the figure configuration.

🤖 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 `@examples/physics/run_binary_sink_disk.py` around lines 526 - 551, Move the
ani.save call in the animation setup flow to after the ax.set_title, axis-label,
limit, aspect, fig.colorbar, and plt.tight_layout configuration calls, so the
generated GIF includes all figure decorations.

579-591: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Use the same spatial extent when rendering and displaying frames.

Frames are generated with render_ext = 15 AU, but render_disk_and_orbit defaults to 2.5 AU. The density image and sink coordinates are therefore plotted on inconsistent scales.

Proposed fix
-    render_disk_and_orbit(render_frames)
+    render_disk_and_orbit(render_frames, ext=rout * 1.5)
📝 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.

    snapshots, render_frames = run_binary_orbit_PN(
        model,
        render=True,
        render_ext=rout * 1.5,
        nx=1024,
        ny=1024,
    )

    for snapshot in snapshots[:3]:
        print("time", snapshot["time"], "positions", snapshot["positions"])

    plot_orbit_trajectory(snapshots)
    render_disk_and_orbit(render_frames, ext=rout * 1.5)
🤖 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 `@examples/physics/run_binary_sink_disk.py` around lines 579 - 591, Update the
render_disk_and_orbit call in the run_binary_orbit_PN flow to pass the same 15
AU spatial extent used by render_ext (rout * 1.5). Preserve the existing frame
rendering while ensuring density images and sink coordinates share the render
extent.
examples/physics/run_sink_integ_test_problems.py (2)

66-69: 🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win

🧩 Analysis chain

🏁 Script executed:

#!/bin/bash
fd -HI 'pyproject.toml|setup.cfg|tox.ini|.python-version|runtime.txt' . -x sh -c 'echo "--- $1"; rg -n "requires-python|python_requires|target-version|python_version" "$1"' sh {}

Repository: Shamrock-code/Shamrock

Length of output: 160


🏁 Script executed:

#!/bin/bash
set -euo pipefail

echo "== relevant files =="
git ls-files 'examples/physics/run_sink_integ_test_problems.py' 'pyproject.toml' 'setup.cfg' 'tox.ini' '.python-version' 'runtime.txt' '.github/workflows/*' | sed 's#^`#-` #'

echo
echo "== pyproject/setup config snippets =="
for f in pyproject.toml setup.cfg tox.ini .python-version runtime.txt; do
  if [ -f "$f" ]; then
    echo "--- $f"
    cat -n "$f" | sed -n '1,220p'
  fi
done

echo
echo "== target file excerpt =="
cat -n examples/physics/run_sink_integ_test_problems.py | sed -n '1,220p'

Repository: Shamrock-code/Shamrock

Length of output: 8481


🏁 Script executed:

#!/bin/bash
set -euo pipefail

echo "== repo root files =="
git ls-files | sed -n '1,220p'

echo
echo "== Python version hints =="
rg -n --hidden --glob '!*.git/*' 'requires-python|python_requires|target-version|python_version|Python 3\.[0-9]+' .

echo
echo "== target file around the flagged lines =="
sed -n '1,220p' examples/physics/run_sink_integ_test_problems.py | nl -ba | sed -n '1,220p'

Repository: Shamrock-code/Shamrock

Length of output: 8687


🏁 Script executed:

#!/bin/bash
set -euo pipefail

echo "== examples/physics/run_sink_integ_test_problems.py =="
cat -n examples/physics/run_sink_integ_test_problems.py | sed -n '1,220p'

Repository: Shamrock-code/Shamrock

Length of output: 7278


Reject mismatched sink-input lengths.

  • zip() truncates silently, so add an explicit length check before adding sinks.
  • The “Circular orbit” setup still passes e=0.2, so the initial condition isn’t circular.
🧰 Tools
🪛 Ruff (0.15.21)

[warning] 66-68: zip() without an explicit strict= parameter

Add explicit value for parameter strict=

(B905)

🤖 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 `@examples/physics/run_sink_integ_test_problems.py` around lines 66 - 69,
Before the sink-creation loop in the relevant setup, validate that positions,
velocities, masses, and accretion_radii have equal lengths and reject mismatches
instead of allowing zip() to truncate. Also update the “Circular orbit” setup’s
eccentricity argument from e=0.2 to e=0 so its initial condition is circular.

Source: Linters/SAST tools


196-202: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Use zero eccentricity for the circular-orbit integration case.

This scenario is labeled and plotted as circular, but e=0.2 initializes an elliptic orbit. Set e=0.0 so the example exercises the stated circular initial conditions.

Proposed fix
-    m1=1.0, m2=m2, a=a, e=0.2, nu=0.0, G=G, roll=0.0, pitch=0.0, yaw=0.0
+    m1=1.0, m2=m2, a=a, e=0.0, nu=0.0, G=G, roll=0.0, pitch=0.0, yaw=0.0
📝 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.

# Circular orbit
m1 = 1.0
m2 = 0.2
a = 1.0
_x1, _x2, _v1, _v2 = shamrock.phys.get_binary_rotated(
    m1=1.0, m2=m2, a=a, e=0.0, nu=0.0, G=G, roll=0.0, pitch=0.0, yaw=0.0
)
🤖 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 `@examples/physics/run_sink_integ_test_problems.py` around lines 196 - 202,
Update the get_binary_rotated call in the circular-orbit setup to pass zero
eccentricity via e=0.0, while leaving the other orbital parameters unchanged.
examples/physics/run_sph_binary_sink.py (8)

31-36: 📐 Maintainability & Code Quality | 🟠 Major | ⚡ Quick win

Translate the added text and remove all non-ASCII characters.

These ranges still contain French comments, labels, and accented characters, which the PR summary says will fail pre-commit. Use English ASCII text throughout.

Also applies to: 53-55, 90-90, 384-390, 405-407, 428-434, 451-453, 490-497, 544-546, 591-593, 615-617, 642-658, 668-699, 736-737

🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 31 - 36, Translate all
French comments, labels, and other user-facing text in
run_sph_binary_sink.py—including the referenced ranges—into clear English, and
replace accented or other non-ASCII characters with ASCII equivalents. Preserve
the existing code behavior and numeric values while ensuring the entire affected
text is ASCII-only.

39-43: 📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win

Correct the physical setup comments.

The spin is tilted by theta = pi/3, not aligned, and the configured periapsis is A * (1 - E) = 90 AU, not approximately 0.7 AU.

Also applies to: 346-346

🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 39 - 43, Update the
physical setup comments near the spin parameters and periapsis configuration:
describe the spins as tilted by theta = pi/3 rather than aligned with orbital
angular momentum, and correct the periapsis value to A * (1 - E) = 90 AU. Change
comments only; preserve the existing parameter values and calculations.

39-47: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Derive spin vectors from each builder's mass arguments.

The parameterized builders attach spin vectors scaled using fixed module-level M1 and M2. Calling either builder with different masses therefore produces the wrong dimensionless spins.

  • examples/physics/run_sph_binary_sink.py#L39-L47: retain dimensionless spin parameters and the axis rather than precomputing mass-dependent vectors globally.
  • examples/physics/run_sph_binary_sink.py#L167-L180: compute vectors from m1 and m2, or accept explicit vectors.
  • examples/physics/run_sph_binary_sink.py#L623-L628: accept m1 and m2 instead of relying on main-block globals.
  • examples/physics/run_binary_sink_disk.py#L97-L105: avoid precomputing mass-dependent vectors globally.
  • examples/physics/run_binary_sink_disk.py#L229-L242: derive or receive spin vectors for the supplied masses.
📍 Affects 2 files
  • examples/physics/run_sph_binary_sink.py#L39-L47 (this comment)
  • examples/physics/run_sph_binary_sink.py#L167-L180
  • examples/physics/run_sph_binary_sink.py#L623-L628
  • examples/physics/run_binary_sink_disk.py#L97-L105
  • examples/physics/run_binary_sink_disk.py#L229-L242
🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 39 - 47, Keep only
dimensionless spin parameters and the shared spin axis in the module-level
setup, avoiding mass-dependent spin vectors. In
examples/physics/run_sph_binary_sink.py lines 167-180, update the parameterized
builder to derive spin vectors from its supplied m1 and m2 values or accept
explicit vectors; in lines 623-628, pass or accept m1 and m2 rather than relying
on main-block globals. Apply the same changes in
examples/physics/run_binary_sink_disk.py lines 97-105 and 229-242 so each
builder derives or receives vectors scaled for its supplied masses.

77-92: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Honor or remove the nu orbital-anomaly parameter.

Both initializers expose nu, but always place the binary at periapsis. Nonzero values silently produce incorrect initial conditions.

  • examples/physics/run_sph_binary_sink.py#L77-L92: restore anomaly-dependent position and velocity calculations, or remove nu.
  • examples/physics/run_binary_sink_disk.py#L135-L150: apply the same contract consistently.
📍 Affects 2 files
  • examples/physics/run_sph_binary_sink.py#L77-L92 (this comment)
  • examples/physics/run_binary_sink_disk.py#L135-L150
🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 77 - 92, Both binary
initializer functions expose nu but always initialize periapsis conditions;
update binary_initial_conditions in examples/physics/run_sph_binary_sink.py
(lines 77-92) and the corresponding initializer in
examples/physics/run_binary_sink_disk.py (lines 135-150) to use
anomaly-dependent position and velocity calculations for nonzero nu, or remove
nu from both APIs and all callers consistently.

570-574: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Return the eccentricity vector without normalizing it.

Normalization discards its magnitude, so plots and CSV columns named e_x, e_y, and e_z do not contain the actual eccentricity-vector components.

Proposed fix
     e_vec = np.cross(v_vec, h)/mu - r_vec/r
-    e=np.linalg.norm(e_vec)
-    e_vec = e_vec / e if e != 0 else np.zeros_like(e_vec)  `#vecteur` excentricité normé ou de Laplace-Runge-Lenz
 
     return e_vec
📝 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.

    e_vec = np.cross(v_vec, h)/mu - r_vec/r

    return e_vec
🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 570 - 574, Update the
eccentricity-vector calculation in the relevant function to return the raw e_vec
from the cross-product expression without dividing by its norm or replacing
zero-magnitude vectors. Preserve the existing component values so e_x, e_y, and
e_z retain the actual eccentricity-vector magnitude.

609-610: 🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win

Actually unwrap the periapsis angle.

Converting degrees to radians and immediately back is an identity operation, so the plotted angle still jumps at 360 degrees.

Proposed fix
-    # retire les sauts de 360°
-    omega = np.degrees((np.radians(omega)))
+    # Remove discontinuities at 360 degrees.
+    omega = np.degrees(np.unwrap(np.radians(omega)))
📝 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.

    # Remove discontinuities at 360 degrees.
    omega = np.degrees(np.unwrap(np.radians(omega)))
🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 609 - 610, Update the
omega processing near the periapsis-angle plotting logic to unwrap the angular
values rather than converting between degrees and radians. Use a
continuous-angle unwrapping operation with the correct unit so jumps at 360
degrees are removed while preserving the plotted angle values.

660-660: 🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

Use a portable ASCII output path.

The author-specific absolute directory will usually not exist, causing np.savetxt to fail. The accented path also violates the ASCII-only requirement.

Proposed fix
-def save_orbital_elements(snapshots, m1, m2, filename="/home/dartencet/Chamrock/CodePythons/tableaux éléments orbitaux/ elements_orbitaux_sham.csv"):
+def save_orbital_elements(snapshots, m1, m2, filename="orbital_elements_sham.csv"):
...
-save_orbital_elements(snapshots, m1, m2, filename="/home/dartencet/Chamrock/CodePythons/tableaux éléments orbitaux/ elements_orbitaux_sham.csv")
+save_orbital_elements(snapshots, m1, m2)

Also applies to: 744-744

🤖 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 `@examples/physics/run_sph_binary_sink.py` at line 660, Update the default
filename in save_orbital_elements to use a portable relative ASCII-only path,
and apply the same correction to the corresponding filename default at the other
reported occurrence. Remove the author-specific absolute directory and accented
characters while preserving the existing output behavior.

738-750: 🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

Guard post-processing against execution during import.

These calls are at module scope, but snapshots, m1, m2, and model are only initialized by the earlier main block. Importing this example therefore raises NameError.

Proposed fix
-plot_eccentricity(snapshots, m1, m2)
-...
-plot_spins(snapshots)
+if __name__ == "__main__":
+    plot_eccentricity(snapshots, m1, m2)
+    # Run the remaining analysis and plotting calls here.
+    plot_spins(snapshots, m1, m2)
🤖 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 `@examples/physics/run_sph_binary_sink.py` around lines 738 - 750, Move the
orbital plotting, save_orbital_elements, model.get_sinks, and plot_spins calls
into the existing main execution guard that initializes snapshots, m1, m2, and
model. Keep their current order and arguments so importing the module performs
no post-processing or model access.
src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp (5)

258-258: 🗄️ Data Integrity & Integration | 🔴 Critical | ⚡ Quick win

🧩 Analysis chain

🏁 Script executed:

rg -n 'compute_ext_forces\(|update_sink_spins\(' src/shammodels/sph/src

Repository: Shamrock-code/Shamrock

Length of output: 729


🏁 Script executed:

#!/bin/bash
set -euo pipefail

# Inspect the relevant code regions with line numbers.
sed -n '230,275p' src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
printf '\n====\n'
sed -n '380,565p' src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp
printf '\n====\n'
sed -n '1688,1735p' src/shammodels/sph/src/Solver.cpp

Repository: Shamrock-code/Shamrock

Length of output: 9078


Move update_sink_spins(dt) out of compute_ext_forces(dt)
compute_ext_forces(dt) runs from both predictor_step and the later force re-evaluation in Solver.cpp, but update_sink_spins(dt) is an accumulating update (+=). That advances spin precession twice per timestep and skews the new spin-evolution path. Keep compute_ext_forces(dt) force-only and call the spin update once per timestep from the solver.

🤖 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/SinkParticlesUpdate.cpp` at line 258, Remove
the update_sink_spins(dt) call from compute_ext_forces so that method remains
force-only. Invoke update_sink_spins(dt) exactly once per timestep from the
solver’s timestep flow, outside both predictor_step and later force
re-evaluation calls, while preserving the existing dt argument and accumulation
behavior.

399-404: 🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win

No guard against >2 sinks with PN terms enabled.

The comment explicitly notes these PN terms "are only true for binary (two sinks)", but pairwise summation is applied regardless of sink_parts.size(). For N>2 sinks with any of OP/SO/SS/RR enabled, this silently produces physically invalid results with no warning.

🛡️ Proposed fix
+    if (sink_parts.size() > 2 && (OP || SO || SS || RR)) {
+        logger::warn_ln(
+            "SinkParticleUpdate",
+            "Post-Newtonian terms are only valid for exactly two sinks; results with more "
+            "sinks will not reproduce correct N-body PN dynamics.");
+    }
📝 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.

    //In the following part of the code, we calculate the acceleration depending of the solver config( Orbital precession, Spin-Orbit, Spin-Spin, Radiation Reaction)
    //Note that all these terms (except for the Newton) are only true for binary (two sinks)
    bool OP = solver_config.compute_OP;
    bool SO = solver_config.compute_SO;
    bool SS = solver_config.compute_SS;
    bool RR = solver_config.compute_RR;
    if (sink_parts.size() > 2 && (OP || SO || SS || RR)) {
        logger::warn_ln(
            "SinkParticleUpdate",
            "Post-Newtonian terms are only valid for exactly two sinks; results with more "
            "sinks will not reproduce correct N-body PN dynamics.");
    }
🤖 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/SinkParticlesUpdate.cpp` around lines 399 -
404, Guard the PN-term setup in the acceleration calculation around OP, SO, SS,
and RR so configurations with more than two sinks cannot apply these binary-only
corrections. When sink_parts.size() exceeds two and any PN term is enabled, emit
an appropriate warning and disable or otherwise prevent those PN contributions
while preserving Newtonian pairwise acceleration.

407-411: 🚀 Performance & Scalability | 🟠 Major | ⚡ Quick win

Per-timestep status logging is called on every force evaluation.

These 5 logger::info_ln calls run unconditionally each time compute_ext_forces is invoked (at least twice per timestep, per the double-invocation issue above), flooding logs over long runs with static configuration values that never change during a run. This is the same block flagged in a prior review ("Avoid if/else spaghetti") and is functionally unchanged (only the field names were renamed).

♻️ Proposed fix

Move this diagnostic to a one-time location (e.g., solver init / SolverConfig::check_config()) instead of the per-timestep force computation:

-    logger::info_ln("-------- SinkParticleUpdate: Post-Newtonian terms --------");
-    logger::info_ln("1PN", solver_config.compute_OP);
-    logger::info_ln("SO", solver_config.compute_SO);
-    logger::info_ln("SS", solver_config.compute_SS);
-    logger::info_ln("RR", solver_config.compute_RR);
🤖 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/SinkParticlesUpdate.cpp` around lines 407 -
411, Move the five Post-Newtonian status logs from the per-timestep
force-evaluation path in compute_ext_forces to a one-time initialization or
SolverConfig::check_config() location. Preserve the reported configuration
values while ensuring they are emitted only once per solver run.

447-450: 🩺 Stability & Availability | 🟡 Minor | ⚡ Quick win

nij unit vector unguarded against exact coincidence (rij_scal == 0).

Unlike update_sink_spins, which adds epsilon_spin before dividing (rij_scal = sycl::length(rij) + epsilon_spin;), nij = rij / rij_scal; here uses the raw distance, so exact position coincidence produces 0/0 (NaN) that would then poison ext_acceleration for all PN terms.

🛡️ Proposed fix
-            Tvec nij = rij / rij_scal;
+            Tvec nij = rij / (rij_scal + epsilon_grav_sink);
📝 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.

            Tvec rij = s1.pos - s2.pos;
            Tscal rij_scal = sycl::length(rij);

            Tvec nij = rij / (rij_scal + epsilon_grav_sink);
🤖 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/SinkParticlesUpdate.cpp` around lines 447 -
450, Guard the unit-vector calculation in the sink-particle update around rij
and nij by adding the established epsilon_spin to rij_scal before dividing.
Preserve the existing rij distance calculation and ensure exact position
coincidence cannot produce 0/0 or propagate NaNs into PN-term ext_acceleration.

497-505: 🎯 Functional Correctness | 🔴 Critical | ⚡ Quick win

Integer-division bug zeroes/truncates the RR (2.5PN) coefficients.

8/5 and 2/3 are both int-literal divisions in C++, truncating to 1 and 0 respectively instead of 1.6 and 0.667. This silently corrupts the radiation-reaction coefficient and drops the 2/3*G*M/r sub-term entirely.

🐛 Proposed fix
-                term4 = 8/5*G*G*eta* M *M/(c*c*c*c*c*(rij_scal*rij_scal*rij_scal+epsilon_grav_sink))
+                term4 = (8.0/5.0)*G*G*eta* M *M/(c*c*c*c*c*(rij_scal*rij_scal*rij_scal+epsilon_grav_sink))
                 * (
-                    vij_nij*nij*(18*v2 + 2/3*G*M/(rij_scal + epsilon_grav_sink)-25*vij_nij*vij_nij) 
+                    vij_nij*nij*(18*v2 + (2.0/3.0)*G*M/(rij_scal + epsilon_grav_sink)-25*vij_nij*vij_nij) 
                     - (6*v2 - 2*G*M/(rij_scal + epsilon_grav_sink)-15*vij_nij*vij_nij)*vij
                 );
📝 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(RR){
                term4 = (8.0/5.0)*G*G*eta* M *M/(c*c*c*c*c*(rij_scal*rij_scal*rij_scal+epsilon_grav_sink))
                * (
                    vij_nij*nij*(18*v2 + (2.0/3.0)*G*M/(rij_scal + epsilon_grav_sink)-25*vij_nij*vij_nij) 
                    - (6*v2 - 2*G*M/(rij_scal + epsilon_grav_sink)-15*vij_nij*vij_nij)*vij
                );
                
                sum += s2.mass/M*term4;
            } 
🤖 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/SinkParticlesUpdate.cpp` around lines 497 -
505, Update the RR term calculation in the visible term4 expression so every
fractional coefficient uses floating-point arithmetic, especially the 8/5
prefactor and both 2/3 gravitational sub-terms. Preserve the existing formula
while ensuring these coefficients evaluate as 1.6 and approximately 0.667 rather
than integer divisions.

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Actionable comments posted: 8

🧹 Nitpick comments (2)
src/shammodels/sph/src/pySPHModel.cpp (1)

274-287: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Use lower_case for the C++ lambda parameters.

Rename C_1_fluid and C_delta_v to c_1_fluid and c_delta_v locally, while preserving the existing py::arg(...) names if those Python keywords are public API.

As per coding guidelines, C++ functions, 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/src/pySPHModel.cpp` around lines 274 - 287, Rename the
lambda parameters C_1_fluid and C_delta_v in the set_monofluid_tvi binding to
c_1_fluid and c_delta_v, and update their uses in the call to
self.dust_config.set_monofluid_tvi. Preserve any existing py::arg names for the
public Python API.

Source: Coding guidelines

examples/physics/run_binary_sink_only.py (1)

205-238: 🚀 Performance & Scalability | 🔵 Trivial | ⚡ Quick win

Make per-step diagnostics optional.

The defaults emit roughly 45,000 lines. Add a verbose or log_every option so console I/O does not dominate the example.

Proposed refactor
-def run_binary_orbit_PN(model, n_steps=n_steps, dt=dt):
+def run_binary_orbit_PN(model, n_steps=n_steps, dt=dt, log_every=None):
...
-    for _ in range(n_steps):
+    for step in range(n_steps):
...
-        print(f"t = {current_time:.4f}, dt = {next_dt:.6f}, distance = {distance:.6f}")
+        if log_every and (step + 1) % log_every == 0:
+            print(f"t = {current_time:.4f}, dt = {next_dt:.6f}, distance = {distance:.6f}")
🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 205 - 238, Make the
per-step diagnostic output in run_binary_orbit_PN optional by adding a verbose
or log_every parameter and guarding the timestep/distance print accordingly.
Preserve the existing default simulation behavior while preventing console
output on every iteration unless explicitly requested.
🤖 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 `@examples/physics/run_binary_sink_only.py`:
- Line 660: Update save_orbital_elements to remove the contributor-specific
absolute default path and use a portable local CSV filename instead; ensure the
corresponding call site also supports the resulting output destination,
optionally allowing a CLI-provided path.
- Line 346: Correct the periapsis distance in the comment near the racc setting
to approximately 90 AU, consistent with A = 100 and E = 0.1, and remove the
incorrect 0.7 AU value.
- Around line 44-47: Update the sink-spin construction in build_binary_sph_model
to derive spin magnitudes from its m1 and m2 parameters rather than the global
M1 and M2 values. Ensure the vectors passed to add_sink for both sinks use these
builder-local masses, including the corresponding logic in the other affected
sink setup.
- Around line 623-628: Update plot_spins to accept m1 and m2 as explicit
parameters, and use those parameters in the spin scaling calculations instead of
relying on script-level globals. Update every call to plot_spins, including the
usage around the additionally referenced location, to pass the corresponding
masses while preserving existing plotting behavior.
- Around line 31-40: Translate all remaining French comments, plot labels,
messages, and path text in run_binary_sink_only.py into English, and replace
every non-ASCII character, including degree symbols, with ASCII equivalents.
Preserve behavior while using terms such as “years,” “deg,” and “orbital
elements” consistently across the referenced sections.
- Around line 340-366: Move the executable setup currently under the `__main__`
guard, including assignments to `m1`, `m2`, `a`, and `e`, model construction via
`build_binary_sph_model`, orbit execution via `run_binary_orbit_PN`, snapshot
printing, and `plot_orbit_trajectory`, into one final `if __name__ ==
"__main__":` block after all function definitions. Ensure no post-processing
code executes on import and that `snapshots`, `m1`, `m2`, and `model` are only
referenced after being initialized.
- Around line 675-677: Update the output column headers associated with the
dimensionless spin values computed as spin1 and spin2 so they are named a1_x
through a1_z and a2_x through a2_z, rather than angular-momentum S labels. Apply
the same header correction in the additional output block around the
corresponding later rows, while leaving the spin calculations unchanged.
- Around line 600-610: Update plot_omega to unwrap the computed periapsis angle
across snapshots before plotting, replacing the identity degrees/radians
conversion with an actual phase-unwrapping operation while keeping the resulting
values in degrees.

---

Nitpick comments:
In `@examples/physics/run_binary_sink_only.py`:
- Around line 205-238: Make the per-step diagnostic output in
run_binary_orbit_PN optional by adding a verbose or log_every parameter and
guarding the timestep/distance print accordingly. Preserve the existing default
simulation behavior while preventing console output on every iteration unless
explicitly requested.

In `@src/shammodels/sph/src/pySPHModel.cpp`:
- Around line 274-287: Rename the lambda parameters C_1_fluid and C_delta_v in
the set_monofluid_tvi binding to c_1_fluid and c_delta_v, and update their uses
in the call to self.dust_config.set_monofluid_tvi. Preserve any existing py::arg
names for the public Python API.
🪄 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: f111491b-b2a9-4892-bb44-9aa4fa9162a8

📥 Commits

Reviewing files that changed from the base of the PR and between b6bfe2f and b23608a.

📒 Files selected for processing (4)
  • examples/physics/run_binary_sink_only.py
  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
  • src/shammodels/sph/src/Solver.cpp
  • src/shammodels/sph/src/pySPHModel.cpp
🚧 Files skipped from review as they are similar to previous changes (2)
  • src/shammodels/sph/src/Solver.cpp
  • src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp

Comment thread examples/physics/run_binary_sink_only.py Outdated
Comment on lines +44 to +47
spin_mag_1 = a1 * G * M1 * M1 / c
spin_mag_2 = a2 * G * M2 * M2 / c
spin_vec_1 = spin_mag_1 * spin_axis
spin_vec_2 = spin_mag_2 * spin_axis

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Derive sink spins from the builder's masses.

add_sink receives vectors precomputed from global M1/M2, so build_binary_sph_model(m1, m2, ...) produces incorrect dimensionless spins whenever its arguments differ from those globals.

Proposed fix
-        tuple(spin_vec_1.tolist()),
+        tuple((a1 * G * m1**2 / c * spin_axis).tolist()),
...
-        tuple(spin_vec_2.tolist()),
+        tuple((a2 * G * m2**2 / c * spin_axis).tolist()),

Also applies to: 167-180

🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 44 - 47, Update the
sink-spin construction in build_binary_sph_model to derive spin magnitudes from
its m1 and m2 parameters rather than the global M1 and M2 values. Ensure the
vectors passed to add_sink for both sinks use these builder-local masses,
including the corresponding logic in the other affected sink setup.

Comment on lines +340 to +366
if __name__ == "__main__":
m1 = M1
m2 = M2
a = A
e = E

# racc=0.001 AU is much smaller than binary separation (~0.7 AU at periapsis)
ctx, model = build_binary_sph_model(
m1,
m2,
a,
e,
roll=0.0,
pitch=0.0,
yaw=0.0,
racc=0.001,
compute_op=True,
compute_so=True,
compute_ss=True,
compute_rr=True,
)
snapshots = run_binary_orbit_PN(model)

for snapshot in snapshots[:3]:
print("time", snapshot["time"], "positions", snapshot["positions"])

plot_orbit_trajectory(snapshots)

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

Move all executable calls into one final main guard.

On import, snapshots, m1, m2, and model are undefined when Lines 738-750 execute. Place the existing setup and all post-processing under a single if __name__ == "__main__": block after the function definitions.

Also applies to: 738-750

🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 340 - 366, Move the
executable setup currently under the `__main__` guard, including assignments to
`m1`, `m2`, `a`, and `e`, model construction via `build_binary_sph_model`, orbit
execution via `run_binary_orbit_PN`, snapshot printing, and
`plot_orbit_trajectory`, into one final `if __name__ == "__main__":` block after
all function definitions. Ensure no post-processing code executes on import and
that `snapshots`, `m1`, `m2`, and `model` are only referenced after being
initialized.

a = A
e = E

# racc=0.001 AU is much smaller than binary separation (~0.7 AU at periapsis)

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win

Correct the stated periapsis distance.

With A = 100 and E = 0.1, periapsis is 90 AU, not approximately 0.7 AU.

🤖 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 `@examples/physics/run_binary_sink_only.py` at line 346, Correct the periapsis
distance in the comment near the racc setting to approximately 90 AU, consistent
with A = 100 and E = 0.1, and remove the incorrect 0.7 AU value.

Comment on lines +600 to +610
def plot_omega(snapshots, m1, m2):

times = np.array([snap["time"] for snap in snapshots])

omega = np.array([
compute_omega(snap, m1, m2)
for snap in snapshots
])

# retire les sauts de 360°
omega = np.degrees((np.radians(omega)))

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Actually unwrap the periapsis angle.

np.degrees(np.radians(omega)) leaves the values unchanged, so the plot retains 360-degree jumps.

Proposed fix
-    omega = np.degrees((np.radians(omega)))
+    omega = np.degrees(np.unwrap(np.radians(omega)))
📝 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.

Suggested change
def plot_omega(snapshots, m1, m2):
times = np.array([snap["time"] for snap in snapshots])
omega = np.array([
compute_omega(snap, m1, m2)
for snap in snapshots
])
# retire les sauts de 360°
omega = np.degrees((np.radians(omega)))
def plot_omega(snapshots, m1, m2):
times = np.array([snap["time"] for snap in snapshots])
omega = np.array([
compute_omega(snap, m1, m2)
for snap in snapshots
])
# retire les sauts de 360°
omega = np.degrees(np.unwrap(np.radians(omega)))
🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 600 - 610, Update
plot_omega to unwrap the computed periapsis angle across snapshots before
plotting, replacing the identity degrees/radians conversion with an actual
phase-unwrapping operation while keeping the resulting values in degrees.

Comment on lines +623 to +628
def plot_spins(snapshots):

times = np.array([s["time"] for s in snapshots])

a1 = c/(G*m1*m1)*np.array([s["spins"][0] for s in snapshots])
a2 = c/(G*m2*m2)*np.array([s["spins"][1] for s in snapshots])

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Pass masses explicitly to plot_spins.

The function depends on globals created only during script execution, making it unusable after import.

Proposed fix
-def plot_spins(snapshots):
+def plot_spins(snapshots, m1, m2):
...
-plot_spins(snapshots)
+plot_spins(snapshots, m1, m2)

Also applies to: 750-750

🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 623 - 628, Update
plot_spins to accept m1 and m2 as explicit parameters, and use those parameters
in the spin scaling calculations instead of relying on script-level globals.
Update every call to plot_spins, including the usage around the additionally
referenced location, to pass the corresponding masses while preserving existing
plotting behavior.

Comment thread examples/physics/run_binary_sink_only.py Outdated
Comment on lines +675 to +677
# spins adimensionnés
spin1 = c/(G*m1*m1) * np.array(snap["spins"][0])
spin2 = c/(G*m2*m2) * np.array(snap["spins"][1])

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

🗄️ Data Integrity & Integration | 🟡 Minor | ⚡ Quick win

Label dimensionless spin columns as a, not S.

The rows contain cS/(Gm^2), but the header labels them as angular momentum components. Rename them to a1_x through a2_z.

Also applies to: 718-723

🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 675 - 677, Update the
output column headers associated with the dimensionless spin values computed as
spin1 and spin2 so they are named a1_x through a1_z and a2_x through a2_z,
rather than angular-momentum S labels. Apply the same header correction in the
additional output block around the corresponding later rows, while leaving the
spin calculations unchanged.

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Actionable comments posted: 1

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (1)
examples/physics/run_binary_sink_only.py (1)

152-176: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick win

Rotate spins with the orbital frame.

binary_initial_conditions rotates positions and velocities for nonzero roll, pitch, or yaw, but the spin vectors passed to add_sink remain unrotated. This breaks the stated spin alignment and gives PN spin terms the wrong initial orientation. Apply the same rotation to both spin vectors before adding the sinks.

🤖 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 `@examples/physics/run_binary_sink_only.py` around lines 152 - 176, Apply the
same orbital-frame rotation used by binary_initial_conditions to spin_vec_1 and
spin_vec_2 before the model.add_sink calls. Update the spin vectors in the
run_binary_sink_only flow so both tuple(spin_vec_*.tolist()) arguments contain
the rotated orientations while preserving their existing magnitudes and sink
setup.
🤖 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 `@examples/physics/run_binary_sink_only.py`:
- Around line 656-657: Add the missing colon to the `save_orbital_elements`
function definition so the `def` statement is syntactically valid; leave the
commented-out filename line unchanged.

---

Outside diff comments:
In `@examples/physics/run_binary_sink_only.py`:
- Around line 152-176: Apply the same orbital-frame rotation used by
binary_initial_conditions to spin_vec_1 and spin_vec_2 before the model.add_sink
calls. Update the spin vectors in the run_binary_sink_only flow so both
tuple(spin_vec_*.tolist()) arguments contain the rotated orientations while
preserving their existing magnitudes and sink setup.
🪄 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: ASSERTIVE

Plan: Pro Plus

Run ID: 01832c53-e811-4f54-887f-5773ca8ad591

📥 Commits

Reviewing files that changed from the base of the PR and between b23608a and 7a3e582.

📒 Files selected for processing (1)
  • examples/physics/run_binary_sink_only.py

Comment thread examples/physics/run_binary_sink_only.py Outdated

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Actionable comments posted: 4

🤖 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 `@examples/physics/run_binary_sink_only.py`:
- Line 612: Update the plot title in the plotting section to use “omega” instead
of “w,” matching the y-axis terminology while preserving the existing title
wording and formatting.
- Line 403: Update the user-facing y-axis label in the diagnostic plot from
“Excentricity e” to the standard spelling “Eccentricity e” at the plt.ylabel
call.
- Line 640: Update the title passed to plt.title in the spin plotting code to
use the grammatically correct text “Evolution of the spin components”.
- Line 384: Update the comment near the h variable in the binary sink setup from
“Specificcinetic momentum” to “Specific angular momentum,” correcting the
spelling and terminology without changing the implementation.
🪄 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: ASSERTIVE

Plan: Pro Plus

Run ID: 9be3febb-b6fd-4889-b889-13f72b0aaf68

📥 Commits

Reviewing files that changed from the base of the PR and between 7a3e582 and d5a1409.

📒 Files selected for processing (1)
  • examples/physics/run_binary_sink_only.py

Comment thread examples/physics/run_binary_sink_only.py
Comment thread examples/physics/run_binary_sink_only.py
Comment thread examples/physics/run_binary_sink_only.py
plt.xlabel("Time")
plt.ylabel("Spin")

plt.title("Evolution of the spins components")

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win

Fix the spin plot title grammar.

Use Evolution of the spin components.

Proposed fix
-    plt.title("Evolution of the spins components")
+    plt.title("Evolution of the spin components")
📝 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.

Suggested change
plt.title("Evolution of the spins components")
plt.title("Evolution of the spin components")
🤖 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 `@examples/physics/run_binary_sink_only.py` at line 640, Update the title
passed to plt.title in the spin plotting code to use the grammatically correct
text “Evolution of the spin components”.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants