Infall moddump from Calcino et al. 2024/subm, and Pineda streamer model #785
Conversation
Summary of ChangesHello @joshcalcino, 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 introduces significant new capabilities for simulating astrophysical phenomena, particularly focusing on the injection and manipulation of material. It provides a robust framework for modeling infall scenarios with precise control over initial conditions and orbital dynamics, and integrates a specialized streamer model for specific astrophysical contexts. The changes also include refinements to existing setup routines, a suite of new utility moddumps for various simulation adjustments, and the crucial addition of HDF5 support for dump files, which will improve data management and interoperability within the simulation workflow. 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. 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 significant new functionality, including new infall models (Calcino et al. 2024/subm, Pineda streamer model) and various moddump utilities. The changes are extensive and add valuable capabilities. However, my review has identified several critical and high-severity issues, primarily stemming from copy-paste errors, incorrect calculations, and unit mismatches in the new modules. I have also pointed out some medium-severity issues related to code duplication and clarity. Addressing these points will improve the robustness and maintainability of the new code.
| n0 = (sqrt(xp(1)**2 + xp(2)**2))/dma | ||
| pf = 2*dma | ||
| n1 = (xp(2)-y1)/dma | ||
| m0 = 2*sqrt(n0-n1-1.0) |
There was a problem hiding this comment.
The calculation m0 = 2*sqrt(n0-n1-1.0) could fail with a domain error if n0-n1-1.0 is negative. You should add a check to ensure the argument of sqrt is non-negative before calling it.
| x0 = -m0*dma | ||
| y0 = dma*(1.0-(x0/pf)**2) | ||
| dx = xyzh_add(1, i) - xp(1) | ||
| y_prime = 4*dma/x0 *dx |
There was a problem hiding this comment.
There is a potential for division by zero here if x0 is zero. x0 is calculated from m0, which can be zero. You should add a check to handle the case where x0 is zero to prevent a runtime error.
| dustfrac_tmp = sum(dustfrac)/npartoftype(igas) | ||
| dustfrac_tmp = sum(dustfrac(1,:))/npartoftype(igas) |
There was a problem hiding this comment.
The calculation for dustfrac_tmp is performed twice. The first line, dustfrac_tmp = sum(dustfrac)/npartoftype(igas), is incorrect for a 2D array and seems to be a leftover. The second line is also suspicious as sum(dustfrac(1,:)) sums over the entire allocated size of the array, not just the active particles, which could lead to incorrect average values. You should sum only over active particles, for example, sum(dustfrac(1, 1:npart)) / npart if you want the average over all active particles, or sum(dustfrac(1, 1:npartoftype(igas))) / npartoftype(igas) for just gas particles.
dustfrac_tmp = sum(dustfrac(1,1:npartoftype(igas)))/npartoftype(igas)
src/main/inject_streamer_pineda.f90
Outdated
| dustfrac_tmp = sum(dustfrac)/npartoftype(igas) | ||
| dustfrac_tmp = sum(dustfrac(1,:))/npartoftype(igas) |
There was a problem hiding this comment.
The calculation for dustfrac_tmp is performed twice. The first line, dustfrac_tmp = sum(dustfrac)/npartoftype(igas), is incorrect for a 2D array and seems to be a leftover. The second line is also suspicious as sum(dustfrac(1,:)) sums over the entire allocated size of the array, not just the active particles, which could lead to incorrect average values. You should sum only over active particles, for example, sum(dustfrac(1, 1:npart)) / npart if you want the average over all active particles, or sum(dustfrac(1, 1:npartoftype(igas))) / npartoftype(igas) for just gas particles.
dustfrac_tmp = sum(dustfrac(1,1:npartoftype(igas)))/npartoftype(igas)
src/main/inject_streamer_pineda.f90
Outdated
| case('sym_stream') | ||
| read(valstring,*,iostat=ierr) sym_stream | ||
| ngot = ngot + 1 | ||
| if (sym_stream < 0 .or. sym_stream > 2) call fatal(label,'sym_stream < 0 or > 2 in input options') |
There was a problem hiding this comment.
The validation for sym_stream only allows values up to 2, but the code in inject_particles handles a value of 3. This will cause a fatal error for a valid option. The check should be updated to allow sym_stream to be 3.
if (sym_stream < 0 .or. sym_stream > 3) call fatal(label,'sym_stream < 0 or > 3 in input options')
| ! - incy : *inclination about the x-axis (degrees)* | ||
| ! - incz : *inclination about the x-axis (degrees)* |
There was a problem hiding this comment.
The descriptions for incy and incz appear to be copy-pasted from incx and incorrectly state they are for inclination about the x-axis. They should be updated to refer to the y-axis and z-axis, respectively.
! - incy : *inclination about the y-axis (degrees)*
! - incz : *inclination about the z-axis (degrees)*
| ! need a better root finder, so right now r_inj must be r0 | ||
| !call theta_at_r(r_inj/rc, theta0m, ecc, xi0, theta) | ||
| theta = theta0m | ||
| phi = phi0m + acos( tan(theta0m) / tan(theta) ) ! Ulrich 1976, eq. (15) |
There was a problem hiding this comment.
The call to acos could result in a domain error if abs(tan(theta0m) / tan(theta)) is greater than 1. While theta is currently set to theta0m, the commented-out call to theta_at_r suggests this might change. It would be safer to add a check to prevent runtime errors.
| print "(7(a5,1x,a,1x,f13.7,/))",& | ||
| 'x','=',xyzmh_ptmass(1,i),& | ||
| 'y','=',xyzmh_ptmass(2,i),& | ||
| 'z','=',xyzmh_ptmass(3,i),& | ||
| 'mass','=',xyzmh_ptmass(4,i),& | ||
| 'h','=',xyzmh_ptmass(ihsoft,i),& | ||
| 'hacc','=',xyzmh_ptmass(ihacc,i),& | ||
| 'Lnuc','=',xyzmh_ptmass(ilum,i),& | ||
| '???','=',xyzmh_ptmass(12,i) |
There was a problem hiding this comment.
The print statement for sink properties appears to print the luminosity twice with different labels. xyzmh_ptmass(ilum,i) (where ilum is 12) is printed with the label Lnuc, and xyzmh_ptmass(12,i) is printed with the label ???. This is redundant and confusing. The line for ??? should be removed.
print "(6(a5,1x,a,1x,f13.7,/))",&
'x','=',xyzmh_ptmass(1,i),&
'y','=',xyzmh_ptmass(2,i),&
'z','=',xyzmh_ptmass(3,i),&
'mass','=',xyzmh_ptmass(4,i),&
'h','=',xyzmh_ptmass(ihsoft,i),&
'hacc','=',xyzmh_ptmass(ihacc,i),&
'Lnuc','=',xyzmh_ptmass(ilum,i)
| xyzit = xyzi+x0 | ||
| vxyzt = vxyz+v0 | ||
|
|
||
| call rotatevec(xyzit,(/1.,0.,0./),incx/180.*pi) | ||
| call rotatevec(vxyzt,(/1.,0.,0./),incx/180.*pi) | ||
|
|
||
| call rotatevec(xyzit,(/0.,1.,0./),incy/180.*pi) | ||
| call rotatevec(vxyzt,(/0.,1.,0./),incy/180.*pi) | ||
|
|
||
| call rotatevec(xyzit,(/0.,0.,1./),incz/180.*pi) | ||
| call rotatevec(vxyzt,(/0.,0.,1./),incz/180.*pi) | ||
|
|
||
| i_part = npart + 1! all particles are new | ||
| call add_or_update_particle(igas, xyzi+x0, vxyz+v0, hguess, u, i_part, npart, npartoftype, xyzh, vxyzu) | ||
| call add_or_update_particle(igas, xyzit, vxyzt, hguess, u, i_part, npart, npartoftype, xyzh, vxyzu) | ||
|
|
||
| xyzit = -xyzi+x0 | ||
| vxyzt = -vxyz+v0 | ||
|
|
||
| call rotatevec(xyzit,(/1.,0.,0./),incx/180.*pi) | ||
| call rotatevec(vxyzt,(/1.,0.,0./),incx/180.*pi) | ||
|
|
||
| call rotatevec(xyzit,(/0.,1.,0./),incy/180.*pi) | ||
| call rotatevec(vxyzt,(/0.,1.,0./),incy/180.*pi) | ||
|
|
||
| call rotatevec(xyzit,(/0.,0.,1./),incz/180.*pi) | ||
| call rotatevec(vxyzt,(/0.,0.,1./),incz/180.*pi) | ||
| i_part = npart + 1! all particles are new | ||
| call add_or_update_particle(igas, -xyzi+x0, -vxyz+v0, hguess, u, i_part, npart, npartoftype, xyzh, vxyzu) | ||
| call add_or_update_particle(igas, xyzit, vxyzt, hguess, u, i_part, npart, npartoftype, xyzh, vxyzu) |
There was a problem hiding this comment.
The block of code for rotating the particle's position and velocity is duplicated for the positive and negative particle pairs. This duplicated code can be refactored into a helper subroutine to improve maintainability and reduce redundancy. The subroutine would take the position and velocity vectors as arguments and apply the series of rotations.
| enddo | ||
| xyzmh_ptmass(1:3, 1) = xyzmh_ptmass(1:3, 1) - xyzmh_ptmass(1:3, 1) | ||
| ! | ||
| print*, 'Fuck.' |
There was a problem hiding this comment.
This print statement contains unprofessional language and should be removed.
1 participant
Description:
Changes include adding a moddump to add in spheres/ellipses on parabolic/hyperbolic orbits following the flyby prescription (for parabolic) and Dullemond+2019 (hyperbolic). Also the Pineda+2009 streamer model used to fit the streamer in Gupta+, Speedie+ papers.
Components modified:
Type of change:
New injection routines, moddump.
Testing:
Did you run the bots? no
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? yes/no
Checked the code compiles and output looks sensible.