Collision between a star (sink particle) and a local section of an accretion disc#801
Conversation
Summary of ChangesHello @tajjankovic, 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 expands the simulation capabilities by introducing new setups and supporting utilities for modeling star-disc collisions, particularly relevant for Quasi-Periodic Eruptions (QPEs). It enables the generation of local accretion disk sections and the simulation of their interaction with a central star (represented as a sink particle), offering both hydrodynamic and radiation-hydrodynamic variants. Highlights
🧠 New Feature in Public Preview: You can now enable Memory to help Gemini Code Assist learn from your team's feedback. This makes future code reviews more consistent and personalized to your project's style. Click here to enable Memory in your admin console. Changelog
Activity
Using Gemini Code AssistThe 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
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 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
|
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request introduces new setups for simulating star-disc collisions, a valuable addition for modeling QPEs, with new modules for setup generation, particle injection, and dump file modification. However, a critical security concern has been identified: two potential buffer overflow vulnerabilities exist where user-supplied or state-dependent counts (number of particles and number of sink particles) are used to index arrays without verifying them against maxp and maxptmass. This could lead to memory corruption or application crashes. Additionally, the review highlights areas for improving maintainability and robustness, such as addressing code duplication in the Makefile, removing unused code, and replacing hardcoded values (especially for sink particle properties) with constants or configurable parameters.
src/utils/moddump_addsink.f90
Outdated
| msink = solarm/umass ! mass of the perturber - 1Msun | ||
| !msink = 0.925 ! for Sun-like polytrope with Rsink = 0.8Rstar | ||
| z0 = 5. ! initial distance of the perturber | ||
| incl = 90./180*3.14192 ! inclination | ||
| vz0 = 1 ! 0 for core | ||
| rsink = 1 !0.8 for core |
There was a problem hiding this comment.
The properties of the sink particle (msink, z0, incl, vz0, rsink) are hardcoded. This makes the utility inflexible. The PR description mentions "Allows user control over sink position and velocity", which suggests these should be configurable. Consider reading these values from a setup file, similar to how set_localdisk.f90 handles its parameters. This would greatly improve the usability of this utility.
| print *,'-----Sink particle distance alog z-direction is:', z0 | ||
| print *,'-----Sink particle velocity is:', vz0 | ||
|
|
||
| nptmass = nptmass + 1 |
There was a problem hiding this comment.
The modify_dump subroutine increments nptmass and uses it to index the xyzmh_ptmass and vxyz_ptmass arrays without verifying if it exceeds the maximum allowed number of sink particles (maxptmass). If the limit is already reached, this will lead to an out-of-bounds memory access and potential memory corruption. It is recommended to import maxptmass from the part module and add a check before incrementing nptmass.
| hfact = 1.2 | ||
|
|
||
| !get the number of particles | ||
| npart = Nparts |
There was a problem hiding this comment.
The setpart subroutine reads the number of particles (Nparts) from the setup file and assigns it to npart without checking if it exceeds the maximum allowed particles (maxp). The subsequent loop from 1 to npart indexes the xyzh and vxyzu arrays, which can lead to a buffer overflow if Nparts is greater than maxp. It is recommended to import maxp from the part module and verify that Nparts <= maxp before proceeding.
build/Makefile_setups
Outdated
| ifeq ($(SETUP), qpe_hyd) | ||
| # hydrodynamic QPE simulation | ||
| SETUPFILE= setup_empty.f90 | ||
| KNOWN_SETUP=yes | ||
| ISOTHERMAL=no | ||
| GRAVITY=no | ||
| SRCINJECT=inject_disk.f90 | ||
| FPPFLAGS= -DLIGHTCURVE | ||
| MODFILE= moddump_addsink.f90 | ||
| endif | ||
|
|
||
| ifeq ($(SETUP), qpe_rh) | ||
| # rad.-hydrodnymic QPE simulation | ||
| RADIATION=yes | ||
| SETUPFILE= setup_empty.f90 | ||
| KNOWN_SETUP=yes | ||
| ISOTHERMAL=no | ||
| GRAVITY=no | ||
| SRCINJECT=inject_disk.f90 | ||
| FPPFLAGS= -DLIGHTCURVE | ||
| MODFILE= moddump_addsink.f90 | ||
| endif |
There was a problem hiding this comment.
The qpe_hyd and qpe_rh targets are nearly identical, with the only difference being the RADIATION=yes flag in qpe_rh. This duplication can be reduced to improve maintainability. Consider using a common variable to share settings between them. For example:
ifeq ($(SETUP), qpe_hyd)
QPE_BASE_SETTINGS = yes
endif
ifeq ($(SETUP), qpe_rh)
QPE_BASE_SETTINGS = yes
RADIATION = yes
endif
ifeq ($(QPE_BASE_SETTINGS),yes)
# (rad.-)hydrodynamic QPE simulation
SETUPFILE= setup_empty.f90
KNOWN_SETUP=yes
ISOTHERMAL=no
GRAVITY=no
SRCINJECT=inject_disk.f90
FPPFLAGS= -DLIGHTCURVE
MODFILE= moddump_addsink.f90
endif
src/main/inject_disk.f90
Outdated
| use part, only:igas,kill_particle,shuffle_part,ihsoft | ||
| use partinject,only:add_or_update_particle | ||
| use io, only:master | ||
| use sortutils, only:indexx | ||
| use vectorutils, only:rotatevec | ||
| use physcon, only:pi,solarm,years,c | ||
| use random, only:ran2 |
There was a problem hiding this comment.
There are several unused modules and variables imported in the inject_particles subroutine. For example, igas, kill_particle, shuffle_part from part, and modules like partinject, io, sortutils, vectorutils, physcon, and random are not used within the subroutine. Removing unused imports improves code clarity and reduces dependencies. The only import that seems to be used is ihsoft from the part module.
use part, only:ihsoft
src/main/inject_disk.f90
Outdated
| ! drift remaining time - new method | ||
| rel_r = rel_r + rel_v*(dt - thit) ! r_new | ||
| rmag = norm2(rel_r) | ||
| if (rmag < star_r) rel_r = rel_r * (star_r/rmag) * 1.00001 ! if still inside (or numerically on surface), push slightly outside |
There was a problem hiding this comment.
The magic number 1.00001 is used to push particles slightly outside the star's surface. It's better to define this as a named constant to improve readability and make it easier to modify if needed.
if (rmag < star_r) rel_r = rel_r * (star_r/rmag) * 1.00001_8 ! if still inside (or numerically on surface), push slightly outside
src/main/inject_disk.f90
Outdated
| endif | ||
|
|
||
| enddo | ||
| print *,'Did bounce (properly):',nbounce,'min_d=',min_d, 'dt=',dtlast |
There was a problem hiding this comment.
This print statement appears to be for debugging purposes. It should be removed from the final code to avoid cluttering the output during production runs.
src/setup/set_localdisk.f90
Outdated
|
|
||
| ! units | ||
| r_unit = solarr | ||
| v_unit = 0.1*3e10 ! 0.1*c - typical velocity in QPEs |
There was a problem hiding this comment.
The speed of light c is hardcoded as 3e10. It's better to use the constant c from the physcon module for better maintainability and accuracy. You'll need to add c to the use physcon statement on line 42, like so: use physcon, only:au,solarm,pi,solarr,c.
v_unit = 0.1*c ! 0.1*c - typical velocity in QPEs
src/utils/moddump_addsink.f90
Outdated
| msink = solarm/umass ! mass of the perturber - 1Msun | ||
| !msink = 0.925 ! for Sun-like polytrope with Rsink = 0.8Rstar | ||
| z0 = 5. ! initial distance of the perturber | ||
| incl = 90./180*3.14192 ! inclination |
There was a problem hiding this comment.
The value of pi is hardcoded as 3.14192. The physcon module, which is already in use, provides a more precise pi constant. Using the constant from physcon will improve accuracy and maintainability.
incl = 90./180*pi ! inclination
There was a problem hiding this comment.
can use the deg_to_rad constant from physcon here
|
thanks for the p-r @tajjankovic, most of the failures here are compiler warnings that need fixing before merge. Can check/eliminate these by compiling with NOWARN=yes see also gemini comments, most are straightforward I think. Let me know if you need help! |
|
@danieljprice Thanks! I think I’ve addressed all the comments (bounds checks, configurable Separate design question (before I go too far): I’m trying to make the star–disc collision more accurate by modelling the stellar envelope as gas while replacing the interior (e.g. r < 0.8 R*) with a single softened core/sink. Current workflow: relax a Sun-like polytrope (~1e7 particles), delete r < 0.78 R*, insert a single particle with M = M_deleted and hsoft ≈ 0.78 R*, and keep a ~1h shell near R_env with frozen properties to prevent inward migration. However, I’m not confident that the approach I’m taking in particle freezing is the best within Phantom. Would you prefer I open a separate draft PR with the work-in-progress for discussion, or discuss the design here first and then implement? |
|
@tajjankovic remaining test failure is a warning when compiling phantomsetup with SETUP=qpe: For the putting-a-core-in-a-star question, we normally handle this automatically inside the set_star procedure -- there's an option to add a sink particle core, currently it only works with setup of MESA or Kepler stars but with a couple of lines of code could easily be made to work for a polytrope as well. The procedure there is to remap the stellar profile to a flattened profile that accounts for the softened potential from the sink particle. This profile is then written to a file and read back to construct the "new" star which will then be close to being in hydrostatic equilibrium. Trying to freeze particles is something we've also tried but it's a tricky to maintain conservation properties. Good project for next phantom workshop in Sydney! |
|
also worth a look is that (for similar reasons to you, to simulate a star crashing through a disc) I added the ability to add a star made of gas to the "grdisc" setup. This uses the set_star procedure so is quite general, suggest to have a look in setup_grdisc.f90 for the relevant code... |
|
@danieljprice I’ve pushed a small fix for the remaining CI failure: dtinject is now assigned a default value at entry to inject_particles. I can’t reproduce the warning locally with gfortran, but hopefully this should address the ifort/ifx check. Thanks for the suggestions. Also, thanks for the Sydney workshop note! I’d completely overlooked that. |
|
remaining CI failure is in the build of phantomtest can be reproduced with (I think you just need a blank set_default_options_inject in your inject module). Thanks again for the contributions! |
2 participants
Description:
Adds new setup targets for constructing a local section of an accretion disc and simulating star-disc collisions (as a model for QPEs). Introduces a local-disk setup generator, a star-disc collision setup wrapper (hydro and radiation-hydro variants), a first-pass injection/collision module for disc–sink interactions, and a moddump utility for adding/initialising a sink particle.
Components modified:
Type of change:
Bug fix
Physics improvements
Better initial conditions
Performance improvements
Documentation update
Better testing
Code cleanup / refactor
Other (please describe)
Summary of changes:
build/Makefile_setups:
Adds three new SETUP targets:
qpe_* uses SRCINJECT=inject_disk.f90, MODFILE=moddump_addsink.f90 and defines -DLIGHTCURVE.
src/setup/set_localdisk.f90:
src/main/inject_disk.f90:
src/utils/moddump_addsink.f90:
Testing:
Did you run the bots? no (CI will run on PR)
Did you update relevant documentation in the docs directory? no
Did you add comments such that the purpose of the code is understandable? yes
Is there a unit test that could be added for this feature/bug? no
Related issues: (none)