APR with GR#679
Conversation
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Summary of Changes
Hello @Fitz-Hu, 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 focuses on enabling and improving the use of Adaptive Particle Refinement (APR) in simulations employing General Relativity (GR). The primary changes involve implementing particle splitting along geodesics in GR, ensuring correct handling of APR particle masses within the GR framework's primitive/conservative conversions and force calculations, and introducing a dedicated module for geodesic integration. Additionally, it includes updates to analysis tools and setup options.
Highlights
- APR-GR Integration: Integrated the Adaptive Particle Refinement (APR) feature with the General Relativity (GR) capabilities of the code.
- Geodesic Particle Splitting: Implemented a new method for splitting particles along geodesics when GR is enabled, replacing the standard tangential plane splitting.
- GR-compatible APR Mass Handling: Fixed bugs related to using APR with GR by ensuring particle mass, particularly the APR-level-dependent mass, is correctly used in primitive-to-conservative conversions and GR force calculations.
- New Geodesic Integration Module: Added a new module (
geodesic.f90) containing a subroutine (integrate_geodesic) to integrate particle motion along geodesics under GR. - Analysis Tool Updates: Modified the radio TDE analysis tool to include internal energy and gravitational potential energy in the shock analysis output.
- Setup File Enhancements: Updated the setup module with new runtime parameters and logic for adding spheres and turbulence around discs.
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 in your pull request via creating an issue comment (i.e. comment on the pull request page) using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands.
| 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 issue 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 is currently in preview and 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 to provide feedback.
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
-
Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configureGemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩
There was a problem hiding this comment.
Code Review
This pull request focuses on integrating Adaptive Particle Refinement (APR) with General Relativity (GR) functionalities. Key changes include:
- A new module
integrate_geodesic.f90for splitting particles along geodesics when GR is active. - Updates to
apr.f90to handle particle splitting differently under GR, utilizing the new geodesic integration. - Modifications across several modules (
cons2prim.f90,substepping.f90,utils_gr.f90) to correctly use particle mass (pmassi) when APR is enabled, ensuring calculations like density and forces are accurate for varying particle resolutions. - Adjustments to analysis scripts (
analysis_radiotde.f90,phantomanalysis.f90) to account for APR and to provide more detailed shock analysis including internal energy. - Updates to build files and dependency lists to incorporate the new GR and APR functionalities.
Overall, the changes appear to be well-targeted for the stated objectives. The review identified a potential build issue in the Makefile, a need for careful initialization of pmassi in one GR path, and a logic change in inject_sim.f90 that needs verification. The new geodesic integration routine is complex and warrants thorough testing. Minor formatting and documentation improvements were also suggested.
| ! | ||
| next_time = -1. | ||
| next_dump = getnextfilename(start_dump) | ||
| next_dump = start_dump |
There was a problem hiding this comment.
The logic for setting next_dump has changed from getnextfilename(start_dump) to start_dump. This implies that the simulation will now start by reading the start_dump itself, rather than the one immediately following it. Ensure this change in behavior is intended and correctly reflects how the initial dump for injection should be handled.
| write(iunit,'(21(A,1x))') '#', 'time', 'rad_min[cm]', 'rad_max[cm]', 'velocity[c]', 'mass[Msun]', 'ki_ene[erg]', 'u[erg]', & !'density[g/cm-3]' | ||
| 'grav[erg]', 'rad_min_tde[cm]', 'rad_max_tde[cm]', 'vel_tde[c]', 'mass_tde[Msun]', & | ||
| 'ki_ene_tde[erg]', 'u_tde[erg]', & | ||
| 'rad_min_cnm[cm]', 'rad_max_cnm[cm]', 'vel_cnm[c]', 'mass_cnm[Msun]', & | ||
| 'ki_ene_cnm[erg]', 'u_cnm[erg]' | ||
| endif | ||
| if (rad_max > 0.) then | ||
| write(iunit,'(16(es18.10,1x))') & | ||
| write(iunit,'(20(es18.10,1x))') & | ||
| time*todays, & | ||
| rad_min*udist, rad_max*udist, shock_v, shock_m*umass/solarm, shock_e*unit_energ, & | ||
| rad_min_tde*udist, rad_max_tde*udist, shock_v_tde, shock_m_tde*umass/solarm, shock_e_tde*unit_energ, & | ||
| rad_min_cnm*udist, rad_max_cnm*udist, shock_v_cnm, shock_m_cnm*umass/solarm, shock_e_cnm*unit_energ !shock_rho*unit_density | ||
| rad_min*udist, rad_max*udist, shock_v, shock_m*umass/solarm, shock_e*unit_energ, shock_u*unit_energ, & | ||
| shock_g*unit_energ, rad_min_tde*udist, rad_max_tde*udist, shock_v_tde, & | ||
| shock_m_tde*umass/solarm, shock_e_tde*unit_energ, shock_u_tde*unit_energ, & | ||
| rad_min_cnm*udist, rad_max_cnm*udist, shock_v_cnm, & | ||
| shock_m_cnm*umass/solarm, shock_e_cnm*unit_energ, shock_u_cnm*unit_energ !shock_rho*unit_density |
There was a problem hiding this comment.
The output format for the 'shock' file has been significantly changed to include 'u[erg]' and 'grav[erg]' for the total shock, and 'u_tde[erg]' and 'u_cnm[erg]' for the TDE and CNM components respectively. Ensure the write statement (line 220-226) correctly maps the new variables (shock_u, shock_g, shock_u_tde, shock_u_cnm) to these new columns and that the format specifier (20(es18.10,1x)) matches the increased number of output variables (now 20, previously 16).
| v = sqrt(dot_product(vxyz,vxyz)) | ||
| tend = dist/v |
There was a problem hiding this comment.
If v (the magnitude of vxyz) is zero or very small, tend could become very large or result in a division by zero if dist is non-zero. Consider adding a check for v > tiny(v) before the division, or handling the case where v is zero (e.g., if dist is also zero, tend is 0; if dist is non-zero and v is zero, this might be an invalid state or require special handling).
| !--------------------------------------------------------------------------! | ||
| ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! | ||
| ! Copyright (c) 2007-2025 The Authors (see AUTHORS) ! | ||
| ! See LICENCE file for usage and distribution conditions ! | ||
| ! http://phantomsph.github.io/ ! | ||
| !--------------------------------------------------------------------------! | ||
| module geodesic | ||
| ! | ||
| ! geodesic | ||
| ! | ||
| ! :References: None | ||
| ! | ||
| ! :Owner: Fitz) Hu | ||
| ! | ||
| ! :Runtime parameters: None | ||
| ! | ||
| ! :Dependencies: cons2primsolver, eos, extern_gr, io, metric_tools, part, | ||
| ! timestep | ||
| ! | ||
|
|
||
| implicit none | ||
|
|
||
| public :: integrate_geodesic | ||
|
|
||
| private | ||
|
|
||
| contains | ||
|
|
||
| subroutine integrate_geodesic(pmass,xyzh,vxyzu,dens,pr,gamma,temp,pxyzu,dist,time) | ||
| use extern_gr, only:get_grforce | ||
| use metric_tools, only:pack_metric,pack_metricderivs | ||
| use eos, only:ieos,equationofstate | ||
| use cons2primsolver,only:conservative2primitive | ||
| use io, only:warning,fatal | ||
| use timestep, only:bignumber,xtol,ptol | ||
| use part, only:rhoh,ien_type | ||
| real, intent(inout) :: xyzh(:),vxyzu(:),pxyzu(:) | ||
| real, intent(inout) :: dens,pr,gamma,temp,pmass | ||
| real, intent(in), optional :: dist,time | ||
| real :: metrics(0:3,0:3,2),metricderivs(0:3,0:3,3),fext(3) | ||
| real :: t,tend,v,dt,hdt | ||
| integer, parameter :: itsmax = 50 | ||
| integer :: its,ierr | ||
| real :: xyz(3),pxyz(3),eni,vxyz(1:3),uui,rho,spsoundi,pondensi | ||
| real :: vxyz_star(3),xyz_prev(3),pprev(3),fstar(3) | ||
| real :: pmom_err,x_err | ||
| logical :: converged | ||
|
|
||
| xyz = xyzh(1:3) | ||
| pxyz = pxyzu(1:3) | ||
| eni = pxyzu(4) | ||
| vxyz = vxyzu(1:3) | ||
| uui = vxyzu(4) | ||
| rho = rhoh(xyzh(4),pmass) | ||
| if (present(dist)) then | ||
| v = sqrt(dot_product(vxyz,vxyz)) | ||
| tend = dist/v | ||
| elseif (present(time)) then | ||
| tend = time | ||
| else | ||
| tend = 0. | ||
| endif | ||
|
|
||
| call pack_metric(xyz,metrics(:,:,:)) | ||
| call pack_metricderivs(xyz,metricderivs(:,:,:)) | ||
| call get_grforce(xyzh(:),metrics(:,:,:),metricderivs(:,:,:),vxyz,dens,uui,pr,fext(1:3),dt) | ||
|
|
||
| t = 0. | ||
| dt = min(0.1,tend*0.1,dt) | ||
| hdt = 0.5*dt | ||
|
|
||
| do while (t <= tend) | ||
| t = t + dt | ||
| pxyz = pxyz + hdt*fext | ||
|
|
||
|
|
||
| ! Note: grforce needs derivatives of the metric, | ||
| ! which do not change between pmom iterations | ||
| its = 0 | ||
| converged = .false. | ||
| pmom_iterations: do while (its <= itsmax .and. .not. converged) | ||
| its = its + 1 | ||
|
|
||
| pprev = pxyz | ||
| call conservative2primitive(xyz,metrics(:,:,:),vxyz,dens,uui,pr,& | ||
| temp,gamma,rho,pxyz,eni,ierr,ien_type) | ||
| if (ierr > 0) call warning('cons2primsolver [in integrate_geodesic (a)]','did not converge') | ||
| call get_grforce(xyzh(:),metrics(:,:,:),metricderivs(:,:,:),vxyz,dens,uui,pr,fstar) | ||
|
|
||
| pxyz = pprev + hdt*(fstar - fext) | ||
| pmom_err = maxval(abs(pxyz - pprev)) | ||
| if (pmom_err < ptol) converged = .true. | ||
| fext = fstar | ||
| enddo pmom_iterations | ||
| if (its > itsmax ) call warning('integrate_geodesic',& | ||
| 'max # of pmom iterations',var='pmom_err',val=pmom_err) | ||
|
|
||
| call conservative2primitive(xyz,metrics(:,:,:),vxyz,dens,uui,pr,temp,& | ||
| gamma,rho,pxyz,eni,ierr,ien_type) | ||
| if (ierr > 0) call warning('cons2primsolver [in integrate_geodesic (b)]','did not converge') | ||
| xyz = xyz + dt*vxyz | ||
| call pack_metric(xyz,metrics(:,:,:)) | ||
|
|
||
| its = 0 | ||
| converged = .false. | ||
| vxyz_star = vxyz | ||
|
|
||
| xyz_iterations: do while (its <= itsmax .and. .not. converged) | ||
| its = its+1 | ||
| xyz_prev = xyz | ||
|
|
||
| call conservative2primitive(xyz,metrics(:,:,:),vxyz_star,dens,uui,& | ||
| pr,temp,gamma,rho,pxyz,eni,ierr,ien_type) | ||
| if (ierr > 0) call warning('cons2primsolver [in integrate_geodesic (c)]','did not converge') | ||
| xyz = xyz_prev + hdt*(vxyz_star - vxyz) | ||
| x_err = maxval(abs(xyz-xyz_prev)) | ||
| if (x_err < xtol) converged = .true. | ||
| vxyz = vxyz_star | ||
| ! UPDATE METRIC HERE | ||
| call pack_metric(xyz,metrics(:,:,:)) | ||
| enddo xyz_iterations | ||
| call pack_metricderivs(xyz,metricderivs(:,:,:)) | ||
| if (its > itsmax ) call warning('integrate_geodesic','Reached max number of x iterations. x_err ',val=x_err) | ||
|
|
||
| call equationofstate(ieos,pondensi,spsoundi,dens,xyzh(1),xyzh(2),xyzh(3),temp,uui) | ||
| pr = pondensi*dens | ||
|
|
||
| call get_grforce(xyzh(:),metrics(:,:,:),metricderivs(:,:,:),vxyzu(1:3),dens,vxyzu(4),pr,fext(1:3),dt) | ||
| pxyzu(1:3) = pxyzu(1:3) + hdt*fext(1:3) | ||
| call conservative2primitive(xyz,metrics(:,:,:),vxyz,dens,uui,& | ||
| pr,temp,gamma,rho,pxyz,eni,ierr,ien_type) | ||
| xyzh(1:3) = xyz(1:3) | ||
| pxyzu(1:3) = pxyz(1:3) | ||
| vxyzu(1:3) = vxyz(1:3) | ||
| vxyzu(4) = uui | ||
| enddo | ||
| end subroutine integrate_geodesic | ||
|
|
||
| end module geodesic |
There was a problem hiding this comment.
This new module geodesic with the integrate_geodesic subroutine is a significant addition for GR simulations.
Key areas to scrutinize:
- Algorithm Correctness: The leapfrog-like integration scheme (predictor-corrector for momentum and position) needs careful validation. Ensure the equations of motion for geodesics in the chosen metric are correctly implemented. The iterative loops for
pmom_iterationsandxyz_iterationssuggest a self-consistent field approach or handling of complex dependencies. Verify convergence criteria (ptol,xtol) and maximum iteration counts (itsmax). - Metric Handling: Calls to
pack_metric,pack_metricderivs, andget_grforceare central. Ensure these are used correctly with updated particle states within the loops. - Variable Updates: Confirm that all relevant particle properties (
xyzh,vxyzu,pxyzu,dens,pr,gamma,temp) are consistently updated and reflect the state after geodesic integration. - Optional Arguments: The handling of
distandtimefor determiningtendis clear. Ensure the defaulttend = 0.is safe if neither is provided (though the loopdo while (t <= tend)would not execute). - Error Handling:
warningcalls are used for convergence failures. Consider iffatalmight be more appropriate in some cases, or if there are recovery mechanisms. - Clarity and Comments: The code is reasonably structured. Adding more comments explaining the physical basis of certain steps or the roles of key variables (e.g.,
fext,fstar) could improve maintainability, especially within the iteration loops.
| call integrate_geodesic(pmass,xyzh(:,npartnew),vxyzu(:,npartnew),dens(npartnew),eos_vars(igasP,npartnew), & | ||
| eos_vars(igamma,npartnew),eos_vars(itemp,npartnew),pxyzu(:,npartnew),dist=sep) | ||
| vxyzu(4,npartnew) = uold ! reset u | ||
| call pack_metric(xyzh(1:3,npartnew),metrics(:,:,:,npartnew)) | ||
| call pack_metricderivs(xyzh(1:3,npartnew),metricderivs(:,:,:,npartnew)) | ||
| call get_grforce(xyzh(:,npartnew),metrics(:,:,:,npartnew),metricderivs(:,:,:,npartnew), & | ||
| vxyzu(1:3,npartnew),dens(npartnew),vxyzu(4,npartnew),eos_vars(igasP,npartnew),fext(1:3,npartnew)) | ||
| if (ind_timesteps) call put_in_smallest_bin(npartnew) |
There was a problem hiding this comment.
The calls to pack_metric, pack_metricderivs, and get_grforce are crucial for GR calculations. Ensure that the arguments passed to these subroutines are correct and in the expected units/frame. Specifically, verify that xyzh(:,npartnew), vxyzu(:,npartnew), dens(npartnew), eos_vars(igasP,npartnew), etc., are all consistent with the requirements of these GR utility functions after the integrate_geodesic call.
| use io, only:fatal,master,id | ||
| use units, only:umass,udist,utime,set_units,unit_density | ||
| use timestep, only:dtmax,tmax | ||
| use timestep, only:dtmax,tmax,idtmax_frac,dtmax_ifactor,idtmax_n |
There was a problem hiding this comment.
The timestep module import now includes idtmax_frac, dtmax_ifactor, and idtmax_n. This indicates changes to how maximum timestep (dtmax) is handled or logged, possibly allowing for fractional or factor-based adjustments.
| dtmax_ifactor = 0 | ||
| idtmax_frac = 0 ! so don't write to .restart | ||
| idtmax_n = 1 |
There was a problem hiding this comment.
New timestep control parameters dtmax_ifactor, idtmax_frac, and idtmax_n are initialized here. The comment for idtmax_frac suggests it's related to restart file writing. Ensure these new parameters are used correctly elsewhere if they influence timestepping behavior or output frequency.
| ! | ||
| ! :Dependencies: boundary, centreofmass, dim, eos, part, physcon, | ||
| ! prompting, readwrite_dumps_fortran, timestep, units | ||
| ! prompting, readwrite_dumps, timestep, units |
There was a problem hiding this comment.
The dependency readwrite_dumps_fortran has been changed to readwrite_dumps. This likely reflects a renaming or consolidation of modules related to dump file I/O. Ensure the new module readwrite_dumps provides all necessary functionality previously offered by readwrite_dumps_fortran for this utility.
| ! :Dependencies: analysis, apr, dim, eos, externalforces, fileutils, | ||
| ! infile_utils, io, kernel, part, readwrite_dumps | ||
| ! | ||
| use dim, only:tagline,do_nucleation,inucleation | ||
| use part, only:xyzh,hfact,massoftype,vxyzu,npart !,npartoftype | ||
| use dim, only:tagline,do_nucleation,inucleation,use_apr | ||
| use part, only:xyzh,hfact,massoftype,vxyzu,npart,apr_level !,npartoftype |
There was a problem hiding this comment.
The phantomanalysis program now imports use_apr from dim and apr_level from part. It also imports init_apr from the apr module. These changes are necessary to correctly initialize and handle APR-related data when analyzing dump files.
| ! initialise apr if it is being used | ||
| if (use_apr) then | ||
| call init_apr(apr_level,ierr) | ||
| else | ||
| apr_level(:) = 1 | ||
| endif |
There was a problem hiding this comment.
This block correctly initializes APR-related data structures if use_apr is true by calling init_apr. If APR is not used, apr_level is defaulted to 1. This is good practice to ensure apr_level is always initialized before potential use in analysis routines.
|
The remaining test failure is genuine and comes from a seg fault in phantomanalysis (with SETUP=star) |
|
|
||
| end subroutine put_in_smallest_bin | ||
|
|
||
| !----------------------------------------------------------------------- |
There was a problem hiding this comment.
as previously discussed, it seems very specific to push a particle forward on a geodesic trajectory by integrating the ODEs rather than just using the local force measurement to push the particle forwards in time. This routine looks ok in that it doesn't copy-and-paste large blocks of code from elsewhere, but are you sure it is necessary compared to just using the local acceleration estimate to extrapolate the velocity?
That is, how far forwards in "time" or "dist" does the integrate_geodesic_gr routine actually push the particle? Is this much more than would occur in one timestep or just a short extrapolation?
|
|
||
| xyzh(1:3) = xyz(1:3) | ||
| vxyzu(1:3) = vxyz(1:3) | ||
| end subroutine integrate_geodesic |
There was a problem hiding this comment.
one small issue with both of these routines is that they are not tested, if they are really needed [see above] then there should be unit tests for them... e.g. they could be compared to the output from calling substep_gr
3 participants
Type of PR:
Bug fix / modification to existing code
Description:
Split particles along geodesics when using gr
Fix bugs with use apr with gr
Testing:
APR of TDE
Did you run the bots? yes
Did you update relevant documentation in the docs directory? no