Skip to content

Conversation

@Guo-astro
Copy link
Contributor

@Guo-astro Guo-astro commented Jan 8, 2026

Summary

Completes the SR-GSPH implementation with Riemann solvers and mode orchestration.

Recovery Methods

  • RecoveryBase: abstract interface for primitive recovery
  • NewtonRaphson: Newton-Raphson iterative recovery algorithm for conservative-to-primitive conversion

Riemann Solvers

  • RiemannBase: abstract SR Riemann solver interface
  • Exact: exact relativistic Riemann solver following Pons et al. (2000)

Mode Implementation

  • SRMode.cpp: complete SR physics mode orchestration
    • Ghost field setup for SR variables (Lorentz factor, momentum, etc.)
    • Primitive recovery from conserved variables
    • Force computation with relativistic corrections
    • Integration of SR equations

Physics Background

  • Implements exact relativistic Riemann solver (Pons 2000)
  • Newton-Raphson primitive recovery for robustness
  • Full SR-GSPH physics loop following Kitajima 2024

Dependencies

Depends on #1525 (SR physics core)

@github-actions
Copy link

github-actions bot commented Jan 8, 2026

Thanks @Guo-astro 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

@gemini-code-assist
Copy link

Summary of Changes

Hello @Guo-astro, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request significantly extends the Special Relativistic Smoothed Particle Hydrodynamics (SR-GSPH) framework by integrating essential components for relativistic fluid simulations. It provides a robust mechanism for interconverting conserved and primitive fluid variables, introduces an exact Riemann solver capable of accurately modeling relativistic wave phenomena, and establishes a comprehensive physics mode to manage the full simulation lifecycle, ensuring both physical accuracy and numerical resilience.

Highlights

  • Primitive Recovery Implementation: Introduced RecoveryBase as an abstract interface for primitive recovery in Special Relativistic Smoothed Particle Hydrodynamics (SR-GSPH). A concrete NewtonRaphson implementation was added to convert conserved variables (momentum, energy, lab-frame density) into primitive variables (Lorentz factor, rest-frame density, pressure, velocities) by solving a quartic equation for the Lorentz factor.
  • Exact Relativistic Riemann Solver: A new RiemannBase interface and an Exact relativistic Riemann solver have been implemented. This solver, based on Pons et al. (2000), uses a hybrid Newton-bisection method to determine the interface states (pressure, normal and tangential velocities) by resolving relativistic shocks and rarefaction waves, including numerical integration via Gauss-Legendre quadrature for rarefactions with tangential velocities.
  • Comprehensive SR Physics Mode: The SRMode class has been added to orchestrate the entire SR-GSPH timestep sequence. This includes managing predictor-corrector steps, handling boundary conditions, performing h-iteration with a volume-based smoothing length, initializing and recovering primitive/conserved variables, computing forces using the new exact Riemann solver, and calculating the CFL timestep.
  • Enhanced Numerical Stability and Diagnostics: Integrated robust checks throughout the SRMode to detect and report NaN/Inf values in primitive recovery and force derivatives, providing early warnings for numerical instabilities. The implementation also includes clamping mechanisms to prevent superluminal velocities, contributing to overall simulation stability.
  • Volume-Based Smoothing Length Iteration: The h-iteration loop for smoothing length calculation has been adapted for SR-GSPH, utilizing a volume-based approach (compute_omega_sr) as described by Kitajima et al. (2025). This method computes the lab-frame baryon density (N) directly from kernel summation and iteratively adjusts the smoothing length based on particle volume.
Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature Command Description
Code Review /gemini review Performs a code review for the current pull request in its current state.
Pull Request Summary /gemini summary Provides a summary of the current pull request in its current state.
Comment @gemini-code-assist Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help /gemini help Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Footnotes

  1. Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution.

Copy link

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

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 a comprehensive Special Relativistic (SR) hydrodynamics mode for GSPH, including new implementations for primitive recovery and an exact Riemann solver. The NewtonRaphson.hpp file provides a Newton-Raphson method to convert conserved variables to primitives, while Exact.hpp implements an iterative SR Riemann solver based on Pons et al. (2000), featuring Gauss-Legendre quadrature for rarefaction waves and robust bracketing for root finding. Corresponding RecoveryBase.hpp and RiemannBase.hpp files define interfaces and configuration structs for these solvers. The SRMode.cpp file orchestrates the SR simulation timestep, managing predictor-corrector steps, volume-based h-iteration, EOS computations, and force calculations using the new SR solvers. Key changes include the introduction of SR-specific conserved fields (momentum 'S', energy 'e') and lab-frame density 'N'.

Review comments primarily focus on improving robustness, consistency, and error handling. Several comments highlight inconsistencies between the Result structs defined in the base interface headers (RecoveryBase.hpp, RiemannBase.hpp) and their concrete implementations (NewtonRaphson.hpp, Exact.hpp), particularly regarding converged and iterations flags, and the lack of proper Config object usage in the solver functions. A major point of feedback concerns the extensive use of printf for error reporting, recommending replacement with structured logging or exceptions for better system integration and debugging. Numerical stability is also addressed, with suggestions to use sycl::fmax to prevent sqrt(negative) and to handle edge cases (e.g., very small denominators, non-positive mass flux) more explicitly rather than silently returning default values. Additionally, the complexity of the first-timestep initialization logic in SRMode.cpp and the direct modification of global configuration are noted as areas for potential streamlining and improved design.

@Guo-astro Guo-astro changed the title [GSPH] Add SR Riemann solver and mode implementation [GSPH 4/7] SR Riemann solver and mode implementation Jan 8, 2026
@Guo-astro
Copy link
Contributor Author

/gemini review

Copy link

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

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 a significant feature: a complete Special Relativistic (SR) Godunov SPH (GSPH) implementation. This includes a Newton-Raphson primitive recovery scheme, an exact SR Riemann solver based on Pons et al. (2000), and the main SRMode to orchestrate the simulation loop. The code is well-structured and includes many robustness checks, which is excellent for numerical hydrodynamics. My review focuses on improving consistency, robustness, and maintainability. Key suggestions include aligning implementations with their abstract interfaces, ensuring consistent use of physical units, replacing printf in device code with more appropriate error handling, and some minor code simplifications.

@Guo-astro Guo-astro force-pushed the pr/gsph-sr-riemann-and-mode branch from 45e892a to e51d1f4 Compare January 9, 2026 07:07
Core Infrastructure:
- PhysicsMode base class for strategy pattern implementation
- ForceKernelBase for common force computation interface
- PhysicsModeFactory for creating physics mode instances
- FieldNames.hpp as SSOT for field naming

Modular Components:
- BoundaryHandler: boundary condition processing
- BuildTrees: tree construction
- ComputeCFL: CFL timestep calculation
- ComputeGradients: gradient computation
- ComputeOmega: omega factor computation
- GhostCommunicator: MPI ghost communication
- IterateSmoothingLengthVolume: h-iteration
- NeighbourCache: caching neighbor interactions
- FunctorNode: generic functor nodes

Refactors the monolithic GSPH Solver into modular components,
preparing for physics mode decoupling (Newtonian vs SR).
The ForceKernelBase template class was designed as a Template Method pattern
base but was never used - NewtonianForceKernel and SRForceKernel are standalone
implementations with their own buffer management appropriate for their physics.

The hardcoded CommonBuffers design (buf_density, buf_pressure, etc.) does not
accommodate SR physics which requires distinct lab-frame vs rest-frame field
naming (N_LABFRAME, LORENTZ_FACTOR, ENTHALPY, etc.).
These fields have physics-specific meanings:
- Newtonian: single frame quantities
- SR: lab-frame vs rest-frame distinction

Each physics mode now defines its own field constants with
appropriate semantic names in their respective FieldNames headers.
Headers:
- NewtonianConfig: solver configuration for Newtonian hydro
- NewtonianEOS: ideal gas equation of state interface
- NewtonianFieldNames: field naming constants
- NewtonianForceKernel: force computation interface
- NewtonianMode: PhysicsMode implementation for Newtonian hydro
- NewtonianTimestepper: CFL-based timestepping
- forces.hpp: force math functions
- ReconstructConfig: interface reconstruction settings
- RiemannConfig: Riemann solver configuration

Riemann Solvers:
- RiemannBase: abstract Riemann solver interface
- HLL: Harten-Lax-van Leer approximate Riemann solver
- Iterative: exact iterative Riemann solver

Sources:
- Complete implementations for all Newtonian components

This implements the Newtonian physics mode for GSPH using
the strategy pattern for physics decoupling.
Energy fields have physics-specific meanings and are now defined
directly in NewtonianFieldNames.hpp rather than imported from
the common FieldNames.hpp.
Headers:
- SRConfig: solver configuration for SR hydro
- SREOS: relativistic equation of state interface
- SRFieldNames: field naming constants for SR fields
- SRForceKernel: force computation interface for SR
- SRMode: PhysicsMode implementation for SR hydro
- SRPrimitiveRecovery: conservative to primitive variable recovery
- SRTimestepper: relativistic CFL-based timestepping
- forces.hpp: SR force math functions

Sources:
- SREOS.cpp: relativistic EOS implementation
- SRForceKernel.cpp: SR force kernel implementation
- SRPrimitiveRecovery.cpp: Newton-Raphson primitive recovery
- SRTimestepper.cpp: SR timestep computation

This implements the core SR physics components for GSPH.
SR physics defines its own XYZ, VXYZ, AXYZ, UINT, DUINT constants
directly rather than importing from the common FieldNames.hpp.
This clarifies the physical meaning (lab-frame quantities for SR).
Recovery Methods:
- RecoveryBase: abstract interface for primitive recovery
- NewtonRaphson: Newton-Raphson iterative recovery algorithm

Riemann Solvers:
- RiemannBase: abstract SR Riemann solver interface
- Exact: exact relativistic Riemann solver (Pons 2000)

Mode Implementation:
- SRMode.cpp: complete SR physics mode orchestration
  - Ghost field setup for SR variables
  - Primitive recovery from conserved variables
  - Force computation with relativistic corrections
  - Integration of SR equations

This completes the Special Relativistic GSPH implementation
following Kitajima 2024 formulation.
Solver Changes:
- Solver.hpp/cpp: Refactored to use PhysicsMode strategy pattern
- SolverConfig.hpp/cpp: Updated config for physics mode selection
- Model.hpp/cpp: Updated model registration

Storage & IO:
- SolverStorage.hpp: Updated field storage for physics modes
- VTKDump.hpp/cpp: Physics-aware VTK output

Python Bindings:
- pyGSPHModel.cpp: Extended bindings for SR configuration
  - physics_mode selection (newtonian/sr)
  - SR-specific parameters (gamma, initial conditions)

Build System:
- CMakeLists.txt: Updated for new physics module structure

This integrates the modular physics modes into the GSPH solver,
enabling runtime selection between Newtonian and SR physics.
Newtonian Tests:
- sod_tube_gsph.py: Sod shock tube validation
- blast_wave_gsph.py: Extreme blast wave test

SR Tests (Kitajima 2024 benchmark suite):
- problem1_sod.py: Relativistic Sod shock tube
- problem2_blast.py: Relativistic blast wave
- problem3_strong_blast.py: Strong relativistic blast
- problem4_ultra_relativistic.py: Ultra-relativistic regime
- problem5_tangent_velocity.py: Tangential velocity test
- problem6_2d_sod.py: 2D relativistic Sod tube
- problem7_kh_instability.py: Kelvin-Helmholtz instability

Common:
- sr/__init__.py: SR test utilities
- kitajima_plotting.py: Plotting helpers for Kitajima benchmarks

All tests use:
- ctx.collect_data() for direct memory access (no pyvista)
- Strict tolerances (~1e-8) for regression testing
- Analytic solutions for validation
Unit Tests:
- GSPHForceTests.cpp: Update for new physics structure
- GSPHRiemannTests.cpp: Update Riemann solver tests

SPH Module Fixes:
- IterateSmoothingLengthDensity: Improve h-iteration logging
- BasicSPHGhosts.cpp: Fix ghost handling

Math:
- sphkernels.hpp: Minor kernel fixes

MHD Placeholder:
- MHDConfig.hpp: Placeholder for future MHD physics mode
- NewtonianMode: add compute_omega_newtonian() using standard SPH (no c_smooth)
- SRMode: use SRIterateSmoothingLength with Kitajima volume-based approach
- Remove shared ComputeOmega module (each mode now owns this)
- Remove legacy UpdateDerivs (replaced by NewtonianForceKernel)
- Move IterateSmoothingLengthVolume to physics/sr/SRIterateSmoothingLength
- Update CMakeLists.txt sources

This fixes the density/pressure error regression caused by c_smooth=1.2
being incorrectly applied to Newtonian mode.
…ult (1.0)

SR's volume-based h-iteration (Kitajima Eq. 232-233) requires c_smooth > 1
to smooth h variation across discontinuities. The SR-specific value was
defined in SRConfig but not transferred to the shared config.
The Riemann solver and force computation code has been moved to
the physics/newtonian/ and physics/sr/ directories. The math/ folder
contained duplicate/orphaned code that is no longer used.
Remove deprecated include of math/riemann/iterative.hpp (now deleted)
and update hllc_solver call to use solve_hll from the new location
in physics/newtonian/riemann/HLL.hpp.
@Guo-astro Guo-astro force-pushed the pr/gsph-sr-riemann-and-mode branch from f4fb4fa to d613a83 Compare January 9, 2026 07:48
Guo-astro and others added 3 commits January 9, 2026 19:56
Run buildbot/update_authors.py to add --no git blame-- annotations
to author headers as required by CI checks.
The kernel summation ρ = ν × ΣW gives the density based on particle
positions in the lab frame. For moving particles, the true lab-frame
baryon density is N = γ × ρ due to Lorentz contraction.

Without this fix, pressure was wrong by factor ~1/γ for particles with
non-zero tangent velocity (e.g., P=436 instead of P=1000 for v_t=0.9).

Co-Authored-By: Claude <noreply@anthropic.com>
Remove SolverCallbacks struct and have PhysicsMode evaluate nodes
directly via storage.solver_graph. This aligns with solvergraph design:
- Branching happens at init time (node registration)
- Flow is visible as graph structure
- No runtime callback creation

Changes:
- Register Solver method nodes in init_solver_graph()
- Delete SolverCallbacks struct from PhysicsMode.hpp
- Update evolve_timestep() signature (remove callbacks param)
- NewtonianMode/SRMode evaluate nodes directly via storage.solver_graph

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
@Guo-astro Guo-astro closed this Jan 9, 2026
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.

2 participants