diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 000000000..656359398 --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,23 @@ +on: [push] + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v3 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: paper/paper.md + - name: Upload + uses: actions/upload-artifact@v1 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: paper/paper.pdf \ No newline at end of file diff --git a/.gitignore b/.gitignore index 19b8038d1..3b52904fb 100644 --- a/.gitignore +++ b/.gitignore @@ -19,11 +19,13 @@ swiftest_driver.sh !README.swifter dump* !**/.gitignore -!.github/workflows/build_wheels.yml +!.github/workflows/*.yml !setup.py !examples/** !swiftest/** !tests/** +!.zenodo.json +!CITATION.cff *ipynb_checkpoints **/.DS_Store !version.txt diff --git a/.zenodo.json b/.zenodo.json new file mode 100644 index 000000000..db7b69cbf --- /dev/null +++ b/.zenodo.json @@ -0,0 +1,58 @@ +{ + "description": "Official release for the Journal of Open Source Software.", + "license" : "GPLv3+", + "title": "profminton/swiftest: v2023.10.2", + "version": "v2023.10.2", + "upload_type": "software", + "publication_date": "2023-10-07", + "creators": [ + { + "name": "Carlisle Wishard", + "affiliation": [ + "Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA", + "Evansville Museum of Arts, History & Science, USA" + ], + "orcid": "0009-0001-0733-3268" + }, + { + "name": "Jennifer Pouplin", + "affiliation": [ + "Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA", + "School of Aeronautics and Astronautics, Purdue University, USA" + ] + }, + { + "name": "Jacob Elliott", + "affiliation": "Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA" + }, + { + "name": "Dana Singh", + "affiliation": "SAIC, Princeton, NJ, USA" + }, + { + "name": "Kaustub Anand", + "affiliation": [ + "Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA", + "Department of Physics and Astronomy, Purdue University, USA" + ] + }, + { + "name": "David Minton", + "affiliation": "Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA", + "orcid": "0000-0003-1656-9704" + } + ], + "access_right": "open", + "related_identifiers": [ + { + "scheme": "url", + "identifier": "https://github.com/profminton/swiftest/tree/v2023.10.2", + "relation": "isSupplementTo" + }, + { + "scheme": "doi", + "identifier": "10.5281/zenodo.8417147", + "relation": "isVersionOf" + } + ] +} diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 000000000..adac9fea2 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,22 @@ +cff-version: 1.1.0 +message: "If you use this software, please cite it as below." +authors: + - name: Wishard + given-names: Carlislse + orcid: https://orcid.org/0009-0001-0733-3268 + - family-names: Pouplin + given-names: Jennifer + - family-names: Elliott + given-names: Jacob + - family-names: Dana + given-names: Singh + - family-names: Anand + given-names: Kaustub + - family-names: Minton + given-names: David + affiliation: "1" + orcid: https://orcid.org/0000-0003-1656-9704 +title: profminton/swiftest: v2023.10.2 +version: v2023.10.2 +date-released: 2023-10-07 + \ No newline at end of file diff --git a/cmake/Modules/SetSwiftestFlags.cmake b/cmake/Modules/SetSwiftestFlags.cmake index 8cbc0de29..8a5390908 100644 --- a/cmake/Modules/SetSwiftestFlags.cmake +++ b/cmake/Modules/SetSwiftestFlags.cmake @@ -563,8 +563,10 @@ IF (CMAKE_BUILD_TYPE STREQUAL "RELEASE" OR CMAKE_BUILD_TYPE STREQUAL "PROFILE") ) # Tells the compiler to link to certain libraries in the Intel oneAPI Math Kernel Library (oneMKL). SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-qmkl=cluster" # Intel - "-qmkl" # Intel + Fortran "-mkl=cluster" + "-mkl" + "-qmkl=cluster" + "-qmkl" ) # Enables additional interprocedural optimizations for a single file compilation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 4905e05f6..b7154f4bb 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -14,22 +14,24 @@ Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation outputs are stored in the subdirectory named after their collisional regime. +**NOTE: You must have ffmpeg installed on your system before running this script. For instance, on MacOS: + +```brew install ffmpeg``` + +on Ubuntu: + +```sudo apt-get install ffmpeg``` + + Inputs _______ -None. +User selects which collisional movie to generate Output ------ -collisions.log : An ASCII file containing the information of any collisional events that occured. -collisions.nc : A NetCDF file containing the collision output. -data.nc : A NetCDF file containing the simulation output. -encounters.nc : A NetCDF file containing the encounter output. -init_cond.nc : A NetCDF file containing the initial conditions for the simulation. -param.00...0.in : A series of parameter input files containing the parameters for the simulation at every output stage. -param.in : An ASCII file containing the inital parameters for the simulation. -param.restart.in : An ASCII file containing the parameters for the simulation at the last output. -swiftest.log : An ASCII file containing the information on the status of the simulation as it runs. -collision.mp4 : A movie file named after the collisional regime depicting the collision. + +Files are stored in directories for each collision type name requested +[COLLISION_TYPE].mp4 : A movie file named after the collisional regime depicting the collision. """ @@ -48,6 +50,22 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) num_movie_frames = 1000 + +# ---------------------------------------------------------------------------------------------------------------------- +# To increase the number of bodies generated in each collision type, decrease the value of the corresponding nfrag_reduction number +# ---------------------------------------------------------------------------------------------------------------------- +nfrag_reduction = {"disruption_headon" : 10.0, + "disruption_off_axis" : 10.0, + "supercatastrophic_headon" : 10.0, + "supercatastrophic_off_axis" : 10.0, + "hitandrun_disrupt" : 10.0, + "hitandrun_pure" : 1.0, + "merge" : 1.0, + "merge_spinner" : 1.0, + } + + + # These initial conditions were generated by trial and error names = ["Target","Projectile"] pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), @@ -124,16 +142,6 @@ "merge_spinner" : 5.0e-3, } -nfrag_reduction = {"disruption_headon" : 1.0, - "disruption_off_axis" : 1.0, - "supercatastrophic_headon" : 1.0, - "supercatastrophic_off_axis" : 10.0, - "hitandrun_disrupt" : 1.0, - "hitandrun_pure" : 1.0, - "merge" : 1.0, - "merge_spinner" : 1.0, - } - density = 3000 * swiftest.AU2M**3 / swiftest.MSun GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3 body_radius = body_Gmass.copy() @@ -358,7 +366,11 @@ def vec_props(self, c): # Set fragmentation parameters minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.50 * body_Gmass[style][1] - sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style]) + sim.set_parameter(collision_model="fraggle", + encounter_save="both", + gmtiny=gmtiny, + minimum_fragment_gmass=minimum_fragment_gmass, + nfrag_reduction=nfrag_reduction[style]) sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/examples/Multibody_Fragmentation/Multibody_Movie.py b/examples/Multibody_Fragmentation/Multibody_Movie.py index 7c3d28359..dde262355 100755 --- a/examples/Multibody_Fragmentation/Multibody_Movie.py +++ b/examples/Multibody_Fragmentation/Multibody_Movie.py @@ -14,6 +14,15 @@ Generates, runs, and processes a set of initial conditions for a multi-body super-catastrophic distruption collisional event. All Swiftest output files are stored in the /supercatastrophic_multi subdirectory. +**NOTE: You must have ffmpeg installed on your system before running this script. For instance, on MacOS: + +```brew install ffmpeg``` + +on Ubuntu: + +```sudo apt-get install ffmpeg``` + + Inputs _______ None. @@ -47,6 +56,16 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) num_movie_frames = 1000 + +# ---------------------------------------------------------------------------------------------------------------------- +# To increase the number of bodies generated in each collision type, decrease the value of the corresponding nfrag_reduction number +# ---------------------------------------------------------------------------------------------------------------------- +nfrag_reduction = { + "supercatastrophic_multi" : 5.0, + } + + + # These initial conditions were generated by trial and error names = ["Body1","Body2","Body3","Body4"] @@ -287,7 +306,12 @@ def vec_props(self, c): # Set fragmentation parameters minimum_fragment_gmass = 0.05 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] - sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="fraggle", + encounter_save="both", + gmtiny=gmtiny, + minimum_fragment_gmass=minimum_fragment_gmass, + nfrag_reduction=nfrag_reduction[style], + verbose=False) sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/paper/paper.bib b/paper/paper.bib index fddf9bdaf..6e6260224 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -1,32 +1,174 @@ -@article{Levison:1994, - Author = {{Levison}, H. and {Duncan}, M.}, - Journal = {Icarus}, - Title = {{The Long-Term Dynamical Behavior of Short-Period Comets}}, - Year = 1994, - Month = mar, - Volume = 108, - DOI = {10.1006/icar.1994.1039}, - url = {https://www.sciencedirect.com/science/article/pii/S0019103584710396?via%3Dihub} +@article{Chambers:1999, + author = {{Chambers}, J.~E.}, + title = "{A hybrid symplectic integrator that permits close encounters between massive bodies}", + journal = {Monthly Notices of the Royal Astronomical Society}, + keywords = {ACCRETION, ACCRETION DISCS, METHODS: NUMERICAL, CELESTIAL MECHANICS, STELLAR DYNAMICS, SOLAR SYSTEM: GENERAL}, + year = 1999, + month = apr, + volume = {304}, + number = {4}, + pages = {793-799}, + doi = {10.1046/j.1365-8711.1999.02379.x}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1999MNRAS.304..793C}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @article{Duncan:1998, - Author = {{Duncan}, M., {Levison}, H., and {Lee}, M. H.}, - Journal = {The Astronomical Journal}, - Title = {{A Multiple Time Step Symplectic Algorithm for Integrating Close Encounters}}, - Year = 1998, - Month = oct, - Volume = 116, - DOI = {10.1086/300541}, - url = {https://iopscience.iop.org/article/10.1086/300541} + title = {A {Multiple} {Time} {Step} {Symplectic} {Algorithm} for {Integrating} {Close} {Encounters}}, + volume = {116}, + url = {http://stacks.iop.org/1538-3881/116/i=4/a=2067}, + doi = {10.1086/300541}, + number = {4}, + journal = {The Astronomical Journal}, + author = {Duncan, Martin J. and Levison, Harold~F. and Lee, Man Hoi}, + month = oct, + year = {1998}, + pages = {2067--2077}, +} + + +@article{Grimm:2022, + title = {{GENGA}. {II}. {GPU} {Planetary} {N}-body {Simulations} with {Non}-{Newtonian} {Forces} and {High} {Number} of {Particles}}, + volume = {932}, + issn = {0004-637X}, + url = {https://dx.doi.org/10.3847/1538-4357/ac6dd2}, + doi = {10.3847/1538-4357/ac6dd2}, + language = {en}, + number = {2}, + urldate = {2023-10-06}, + journal = {The Astrophysical Journal}, + author = {Grimm, Simon L. and Stadel, Joachim G. and Brasser, Ramon and Meier, Matthias M. M. and Mordasini, Christoph}, + month = jun, + year = {2022}, + pages = {124}, +} + + +@misc{kaufmann_swifter_nodate, + title = {Swifter}, + url = {https://www.boulder.swri.edu/swifter/}, + abstract = {The Swifter software package, written by David E. Kaufmann, is a completely redesigned and improved version of the Swift package of Hal Levison and Martin Duncan. Like Swift, Swifter is designed to integrate a set of mutually gravitationally interacting bodies together with a group of massless test particles that feel the gravitational influence of the massive bodies but do not affect each other or the massive bodies. In addition, the SyMBA integrator supports a second class of massive bodies whose masses are less than some user-specified value. These bodies gravitationally interact with the more massive bodies, but do not interact with each other. Seven integration techniques are included in the current beta release of Swifter:}, + urldate = {2023-10-06}, + journal = {Swifter — an improved solar system integration software package}, + author = {Kaufmann, David E.}, } @article{Leinhardt:2012, - Author = {{Leinhardt}, Z. and {Stewart}, S.}, - Journal = {The Astronomical Journal}, - Title = {{Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling Laws.}}, - Year = 2012, - Month = dec, - Volume = 745, - DOI = {10.1088/0004-637X/745/1/79}, - url = {https://iopscience.iop.org/article/10.1088/0004-637X/745/1/79} + title = {Collisions between {Gravity}-dominated {Bodies}. {I}. {Outcome} {Regimes} and {Scaling} {Laws}}, + volume = {745}, + url = {http://stacks.iop.org/0004-637X/745/i=1/a=79}, + doi = {10.1088/0004-637X/745/1/79}, + number = {1}, + journal = {The Astrophysical Journal}, + author = {Leinhardt, Zoë M and Stewart, Sarah T}, + month = jan, + year = {2012}, + pages = {79}, +} + + +@article{Levison:1994, + author = {{Levison}, Harold F. and {Duncan}, M.}, + title = {{The Long-Term Dynamical Behavior of Short-Period Comets}}, + journal = {Icarus}, + year = 1994, + month = mar, + volume = 108, + doi = {10.1006/icar.1994.1039}, + url = {https://www.sciencedirect.com/science/article/pii/S0019103584710396?via%3Dihub} +} + + +@article{Levison:2000, + title = {Symplectically {Integrating} {Close} {Encounters} with the {Sun}}, + volume = {120}, + issn = {0004-6256}, + url = {http://dx.doi.org/10.1086/301553}, + doi = {10.1086/301553}, + number = {4}, + journal = {The Astronomical Journal}, + author = {{Levison}, Harold F. and {Duncan}, Martin J.}, + month = oct, + year = {2000}, + pages = {2117--2123}, +} + +@article{Moore:2011, + title = {{QYMSYM}: {A} {GPU}-accelerated hybrid symplectic integrator that permits close encounters}, + volume = {16}, + url = {http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2011NewA...16..445M&link_type=ABSTRACT}, + doi = {10.1016/j.newast.2011.03.009}, + number = {7}, + journal = {New Astronomy}, + author = {Moore, Alexander and Quillen, Alice C.}, + month = nov, + year = {2011}, + pages = {445--455}, +} + +@article{Rein:2012, + title = {{REBOUND}: an open-source multi-purpose {N}-body code for collisional dynamics}, + volume = {537}, + url = {http://www.aanda.org/10.1051/0004-6361/201118085}, + doi = {10.1051/0004-6361/201118085}, + journal = {Astronomy \& Astrophysics}, + author = {{Rein}, Hanno and {Liu}, S. F.}, + month = jan, + year = {2012}, + pages = {A128}, +} + + +@article{Rein:2015, + title = {WHFAST: a fast and unbiased implementation of a symplectic {Wisdom}–{Holman} integrator for long-term gravitational simulations}, + volume = {452}, + issn = {0035-8711}, + shorttitle = {whfast}, + url = {https://doi.org/10.1093/mnras/stv1257}, + doi = {10.1093/mnras/stv1257}, + number = {1}, + urldate = {2023-10-06}, + journal = {Monthly Notices of the Royal Astronomical Society}, + author = {{Rein}, Hanno and {Tamayo}, Daniel}, + month = sep, + year = {2015}, + pages = {376--388}, +} + +@article{Saha:1994, + title = {Long-term planetary integration with individual time steps}, + volume = {108}, + url = {http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1994AJ....108.1962S&link_type=ABSTRACT}, + doi = {10.1086/117210}, + journal = {The Astronomical Journal}, + author = {Saha, Prasenjit and Tremaine, Scott}, + month = nov, + year = {1994}, + pages = {1962--1969}, +} + +@article{Stewart:2012, + title = {Collisions between {Gravity}-dominated {Bodies}. {II}. {The} {Diversity} of {Impact} {Outcomes} during the {End} {Stage} of {Planet} {Formation}}, + volume = {751}, + url = {http://stacks.iop.org/0004-637X/751/i=1/a=32}, + doi = {10.1088/0004-637X/751/1/32}, + number = {1}, + journal = {The Astrophysical Journal}, + author = {Stewart, Sarah T and Leinhardt, Zoë M}, + month = jan, + year = {2012}, + pages = {32}, +} + + +@article{Wisdom:1991, + title = {Symplectic maps for the n-body problem}, + volume = {102}, + url = {http://adsabs.harvard.edu/cgi-bin/bib_query?1991AJ....102.1528W}, + doi = {10.1086/115978}, + journal = {The Astronomical Journal}, + author = {Wisdom, Jack and Holman, Matthew}, + month = oct, + year = {1991}, + pages = {1528--1538}, } diff --git a/paper/paper.md b/paper/paper.md index f5c8c6310..402c75632 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -9,48 +9,77 @@ tags: - Planetary Systems authors: - name: Carlisle Wishard + affiliation: "1,4" orcid: 0009-0001-0733-3268 equal-contrib: true - corresponding: true - affiliation: 1 - - name: David Minton - equal-contrib: true - affiliation: 1 - name: Jennifer Pouplin - affiliation: "1" - - name: Jake Elliott + affiliation: "1,3" + - name: Jacob Elliott affiliation: "1" - name: Dana Singh + affiliation: "5" + - name: Kaustub Anand + affiliation: "1, 2" + - name: David Minton affiliation: "1" + orcid: 0000-0003-1656-9704 + equal-contrib: true + corresponding: true affiliations: - - name: Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA - index: 1 -date: 05 April 2023 + - name: Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, USA + index: 1 + - name: School of Aeronautics and Astronautics, Purdue University, USA + index: 2 + - name: Department of Physics and Astronomy, Purdue University, USA + index: 3 + - name: Evansville Museum of Arts, History & Science, USA + index: 4 + - name : SAIC, Princeton, NJ, USA + index: 5 + +date: 06 October 2023 bibliography: paper.bib --- # Summary - -The dynamical evolution of planetary systems is dominated by gravitational interactions between massive bodies. Determining the orbits of massive bodies over long time scales is the first step towards understanding the formation and evolution of planets, moons, asteroids, comets, and more. To model these systems, which often include hundreds or thousands of gravitationally interacting bodies, a numerical tool called an \textit{n}-body integrator is often employed. +`Swiftest` is a software package designed to model the long-term dynamics of system of bodies in orbit around a dominant central body, such a planetary system around a star, or a satellite system around a planet. The main body of the program is written in Modern Fortran, taking advantage of the object-oriented capabilities included with Fortran 2003 and the parallel capabilities included with Fortran 2008 and Fortran 2018. `Swiftest` also includes a Python package that allows the user to quickly generate input, run simulations, and process output from the simulations. `Swiftest` uses a NetCDF output file format which makes data analysis with the `Swiftest` Python package a streamlined and flexible process for the user. # Statement of Need +Building off a strong legacy, including its predecessors `Swifter` [@kaufmann_swifter_nodate] and `Swift` [@Levison:1994], `Swiftest` takes the next step in modeling the dynamics of planetary systems by improving the performance and ease of use of software, and by introducing a new collisional fragmentation model. Currently, `Swiftest` includes the four main symplectic integrators included in its predecessors: -`Swiftest` is a software package designed to model gravitationally dominated systems. The main body of the program is written in Modern Fortran, taking advantage of the object-oriented capabilities included with Fortran 2003 and the parallel capabilities included with Fortran 2008 and Fortran 2018. `Swiftest` also includes a Python package that allows the user to quickly generate input and process output from the main integrator. `Swiftest` uses a NetCDF output file format which makes data analysis with the `Swiftest` Python package a streamlined and flexible process for the user. +WHM + : Wisdom-Holman method. `WHM` is a symplectic integrator that is best suited for cases in which the orbiting bodies do not have close encounters with each other. For details see @Wisdom:1991. -Building off a strong legacy, including its predecessors `Swifter` [@Duncan:1998] and `Swift` [@Levison:1994], `Swiftest` takes the next step in modeling gravitationally dominated systems by including collisional fragmentation. Our collisional fragmentation algorithm, `Fraggle` (based on the work of @Leinhardt:2012), is designed to resolve collisions between massive bodies and generate collisional debris. `Swiftest` fully incorporates this debris into the gravitational system, evolving these new bodies along with pre-existing bodies. This allows for a more complete model of the orbital evolution of the system and the growth of massive bodies. +RMVS + : Regularized Mixed Variable Symplectic. `RMVS` is an extension of `WHM` that handles close approaches between test particles and massive bodies. For details, see @Levison:1994. -The combination of modern programming practices, flexible data processing tools, and the latest advances in the field of collisional dynamics make `Swiftest` the ideal tool for studying the formation of planetary systems, the growth of planetary moons, the evolution of asteroid families, and beyond. +HELIO + : Democratic Heliocentric method. This is a basic symplectic integrator that uses democratic heliocentric coordinates instead of the Jacobi coordinates used by `WHM`. Like `WHM` it is not suited for simulating close encounters between orbiting bodies. For details, see @Duncan:1998. -# Performance +SyMBA + : Symplectic Massive Body Algorithm. This is an extension of `HELIO` that handles close approaches between massive bodies and any of the other objects in the simulation. It also includes semi-interacting massive bodies that can gravitationally influence fully massive bodies but not each other. This algorithm is described in the @Duncan:1998. See also @Levison:2000. + +All of the integrators that are currently implemented are based symplectic integrators, which are designed to model massive bodies (e.g. planets) or test particles (e.g. asteroids) in orbit of a single dominant central body (e.g. the Sun) for very long periods of times (e.g. Gy). The core components of `Swiftest` were inherited from the `Swift`/`Swifter` codebase, but have been somewhat modified for performance. The `Swift`-family of integrators are most comparable to other symplectic integrators, which are among the most popular numerical methods used to study the long-term dynamics of planetary systems. + +The `SyMBA` integrator included in `Swiftest` is most similar to the hybrid symplectic integrator `MERCURY6` [@Chambers:1999], the `MERCURIUS` integrator of `REBOUND` [@Rein:2012;@Rein:2015], and the GPU-enabled hybrid symplectic integrators such as `QYMSYM` [@Moore:2011] and `GENGA II` [@Grimm:2022], with some important distinctions. The hybrid symplectic integrators typically employ a symplectic method, such as the original WHM method in Jacobi coordinates or the modified method that uses the Democratic Heliocentric coordinates, only when bodies are far from each other relative to their gravitational spheres of influence (some multiple of a Hill's sphere). When bodies approach each other, the integration is smoothly switched to a non-symplectic method, such as Bulirsch-Stoer or IAS15. In contrast, `SyMBA` is a multi-step method that recursively subdivides the step size of bodies undergoing close approach with each other. +In addition `Swiftest` contains a number of significant enhancements relative to its predecessors: + +- It includes a fully symplectic implementation of the post-Newtonian correction (General Relativity) as described in @Saha:1994. +- Output data is stored in NetCDF4 (HDF5) format, making data analysis, cross-platform compatibility, and portability far more robust than the flat binary file formats used by its predecessors. +- It makes use of CPU-based parallelization via OpenMP. It also makes use of SIMD instructions, but these are available only when compiled for target host machine rather than when installed from PyPI via pip. +- Simulations can be run by means of a standalone binary driver program, similar to its predecessors, or from a compiled Python module. Both the standalone driver and Python module link to the same compiled library, so simulations run by either method are identical. +- It comes with an extensive set of Python scripts to help generate simulation initial conditions and post-process simulation results. +- It includes the `Fraggle` collisional fragmentation model that is used to generate collisional fragments on trajectories that conserve linear and angular momentum and lose the appropriate amount of collisional energy for inelastic collisions between massive bodies in `SyMBA` simulations. The collisional outcome is determined using standard methods based on the work of @Leinhardt:2012 and @Stewart:2012. + +# Performance Modeling the behavior of thousands of fully interacting bodies over long timescales is computationally expensive, with typical runs taking weeks or months to complete. The addition of collisional fragmentation can quickly generate hundreds or thousands of new bodies in a short time period, creating further computational challenges for traditional \textit{n}-body integrators. As a result, enhancing computational performance was a key aspect of the development of `Swiftest`. Here we show a comparison between the performance of `Swift`, `Swifter-OMP` (a parallel version of `Swifter`), and `Swiftest` on simulations with 1k, 2k, 8k, and 16k fully interacting bodies. The number of cores dedicated to each run is varied from 1 to 24 to test the parallel performance of each program. \autoref{fig:performance} shows the results of this performance test. We can see that `Swiftest` outperforms `Swifter-OMP` and `Swift` in each simulation set, even when run in serial. When run in parallel, `Swiftest` shows a significant performance boost when the number of bodies is increased. The improved performance of `Swiftest` compared to `Swifter-OMP` and `Swift` is a critical step forward in \textit{n}-body modeling, providing a powerful tool for modeling the dynamical evolution of planetary systems. -![Performance testing of `Swiftest` on systems of (a) 1k, (b) 2k, (c) 8k, and (d) 16k fully interacting massive bodies. All simulations were run using the \textit{SyMBA} integrator included in `Swift`, `Swifter-OMP`, and `Swiftest`. Speedup is measured relative to `Swift` (dashed), with an ideal 1:1 speedup relative to `Swiftest` in serial shown as an upper limit (dotted). The performance of `Swifter-OMP` is shown in green while the performance of `Swiftest` is shown in blue. All simulations were run on the Purdue University Rosen Center for Advanced Computing Brown Community Cluster. Brown contains 550 Dell compute nodes, with each node containing 2 12-core Intel Xeon Gold Sky Lake processors, resulting in 24 cores per node. Each node has 96 GB of memory. \label{fig:performance}](performance.png) +![Performance testing of `Swiftest` on systems of (a) 1k, (b) 2k, (c) 8k, and (d) 16k fully interacting massive bodies. All simulations were run using the \textit{SyMBA} integrator included in `Swift`, `Swifter-OMP` (and earlier attempt to parallelize `Swifter`), and `Swiftest`. Speedup is measured relative to `Swiftest` for 1 CPU (dashed), with an ideal 1:1 speedup shown as an upper limit (dotted). The performance of `Swifter-OMP` is shown in green while the performance of `Swiftest` is shown in blue. All simulations were run on the Purdue University Rosen Center for Advanced Computing Brown Community Cluster. Brown contains 550 Dell compute nodes, with each node containing 2 12-core Intel Xeon Gold Sky Lake processors, resulting in 24 cores per node. Each node has 96 GB of memory. \label{fig:performance}](performance.png) # Acknowledgements - `Swiftest` was developed at Purdue University and was funded under the NASA Emerging Worlds and Solar System Workings programs. Active development by the Purdue Swiftest Team is ongoing and contributions from the community are highly encouraged. # References diff --git a/paper/performance.png b/paper/performance.png index 7a425e840..7f95d781e 100644 Binary files a/paper/performance.png and b/paper/performance.png differ diff --git a/pyproject.toml b/pyproject.toml index 5e76778a6..0e7fe3767 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "swiftest" -version = "2023.10.0" +version = "2023.10.2" authors=[ {name = 'David A. Minton', email='daminton@purdue.edu'}, {name = 'Carlisle Wishard'}, @@ -22,7 +22,7 @@ classifiers=[ 'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)', 'Programming Language :: Python :: 3', ] -keywords=['astronomy','astrophysics', 'planetary', 'nbody integrator', 'symplectic', 'wisdom-holman'] +keywords=['astronomy','astrophysics', 'planetary', 'n-body', 'integrator', 'symplectic', 'wisdom-holman', 'symba'] dependencies = [ 'numpy>=1.24.3', 'scipy>=1.10.1', @@ -32,13 +32,12 @@ dependencies = [ 'bottleneck>=1.3', 'h5netcdf>=1.1', 'h5py>=3.9', - 'netcdf4>=1.6.4', + 'netcdf4>=1.6', 'matplotlib>=3.7', 'astropy>=5.2', 'astroquery>=0.4.6', 'tqdm>=4.66', 'cython>=3.0.0', - 'pyshtools>=4.10' ] [project.urls] @@ -99,14 +98,14 @@ netCDF-Fortran_DIR="${PREFIX}/lib/cmake/netCDF" [tool.cibuildwheel.macos] before-all = [ - "brew install coreutils shtools", - "LIBS=\"\" buildscripts/build_dependencies.sh -p ${PREFIX} -d ${PREFIX}/build -m ${MACOSX_DEPLOYMENT_TARGET}" + "brew install coreutils", + "LIBS=\"\" buildscripts/build_dependencies.sh -p ${PREFIX} -d ${HOME}/Downloads -m ${MACOSX_DEPLOYMENT_TARGET}" ] [tool.cibuildwheel.linux] skip = "cp312-* pp* -manylinux_i686* *-musllinux*" before-all = [ - "yum install doxygen libxml2-devel libcurl-devel fftw-devel blas-devel lapack-devel -y", + "yum install doxygen libxml2-devel libcurl-devel -y", "buildscripts/build_dependencies.sh -p /usr/local" ] [tool.cibuildwheel.linux.environment] diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index aeeae27d9..746e768d2 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -635,7 +635,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, vdisp, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new real(DP), dimension(NDIM) :: dL_metric real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_avail, ke_remove, dE_best, fscale - real(DP) :: dE_metric, mfrag, rn, dL1_mag, dE_conv, vumag + real(DP) :: dE_metric, mfrag, rn, dL1_mag, dE_conv, vumag, L_mag_factor integer(I4B), dimension(:), allocatable :: vsign real(DP), dimension(:), allocatable :: vscale real(DP), dimension(:), allocatable :: dLi_mag @@ -644,8 +644,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess - integer(I4B), parameter :: MAXLOOP = 25 - integer(I4B), parameter :: MAXTRY = 10 + integer(I4B), parameter :: MAXLOOP = 50 + integer(I4B), parameter :: MAXTRY = 50 integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -754,9 +754,10 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Try to put residual angular momentum into the spin, but if this would go past the spin barrier, then put it into ! velocity shear instead call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) + L_mag_factor = .mag.(collider_local%L_total(:,1) + collider_local%L_total(:,2)) + L_residual(:) = (collider_local%L_total(:,2) / L_mag_factor - collider_local%L_total(:,1) / L_mag_factor) L_residual_unit(:) = .unit. L_residual(:) - if (nsteps == 1) L_residual_best(:) = L_residual(:) + if (nsteps == 1) L_residual_best(:) = L_residual(:) * L_mag_factor ! Use equipartition of spin kinetic energy to distribution spin angular momentum #ifdef DOCONLOC @@ -768,7 +769,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu (fragments%radius(i) / fragments%radius(istart))**2 * & (fragments%Ip(3,i) / fragments%Ip(3,istart)))**(1.5_DP) end do - dL1_mag = .mag.L_residual(:) / sum(dLi_mag(istart:fragments%nbody)) + dL1_mag = .mag.L_residual(:) * L_mag_factor / sum(dLi_mag(istart:fragments%nbody)) do i = istart,fragments%nbody dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:) @@ -789,20 +790,22 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i) end if end if - L_residual(:) = L_residual(:) + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + L_residual(:) = L_residual(:) + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 & + / L_mag_factor end do ! Put any remaining residual into velocity shear angmtm: do j = 1, MAXANGMTM if (j == MAXANGMTM) exit inner call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric(:) = abs(L_residual(:)) / .mag.(collider_local%L_total(:,1)) / MOMENTUM_SUCCESS_METRIC + L_mag_factor = .mag.(collider_local%L_total(:,1) + collider_local%L_total(:,2)) + L_residual(:) = (collider_local%L_total(:,2) / L_mag_factor - collider_local%L_total(:,1)/L_mag_factor) + dL_metric(:) = abs(L_residual(:)) / MOMENTUM_SUCCESS_METRIC if (all(dL_metric(:) <= 1.0_DP)) exit angmtm do i = istart, fragments%nbody - dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody)) + dL(:) = -L_residual(:) * L_mag_factor * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody)) call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) end do call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) @@ -817,8 +820,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual_last = E_residual E_residual = dE + impactors%Qloss - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - dL_metric(:) = abs(L_residual(:)) / .mag.collider_local%L_total(:,1) / MOMENTUM_SUCCESS_METRIC + L_mag_factor = .mag.(collider_local%L_total(:,1) + collider_local%L_total(:,2)) + L_residual(:) = (collider_local%L_total(:,2) / L_mag_factor - collider_local%L_total(:,1) / L_mag_factor) + dL_metric(:) = abs(L_residual(:)) / MOMENTUM_SUCCESS_METRIC dE_conv = abs(E_residual - E_residual_last) / abs(E_residual_last) ! Check if we've converged on our constraints @@ -826,7 +830,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity E_residual_best = E_residual - L_residual_best(:) = L_residual(:) + L_residual_best(:) = L_residual(:) * L_mag_factor dE_best = dE nsteps_best = nsteps @@ -877,8 +881,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu // trim(adjustl(message)) // " steps.") call collider%get_energy_and_momentum(nbody_system, param, phase="after") - L_residual(:) = (collider%L_total(:,2) - collider%L_total(:,1)) - call collision_util_velocity_torque(-L_residual(:), collider%fragments%mtot, impactors%rbcom, impactors%vbcom) + L_mag_factor = .mag.(collider%L_total(:,1) + collider%L_total(:,2)) + L_residual(:) = (collider%L_total(:,2) / L_mag_factor - collider%L_total(:,1)) / L_mag_factor + call collision_util_velocity_torque(-L_residual(:) * L_mag_factor, collider%fragments%mtot, impactors%rbcom, impactors%vbcom) #ifdef DOCONLOC do concurrent(i = 1:nfrag) shared(collider, impactors) diff --git a/src/globals/globals_module.f90 b/src/globals/globals_module.f90 index d1a078329..d7d00e085 100644 --- a/src/globals/globals_module.f90 +++ b/src/globals/globals_module.f90 @@ -48,7 +48,7 @@ module globals integer(I4B), parameter :: UPPERCASE_OFFSET = iachar('A') - iachar('a') !! ASCII character set parameter for lower to upper !! conversion - offset between upper and lower - character(*), parameter :: VERSION = "2023.10.0" !! Swiftest version + character(*), parameter :: VERSION = "2023.10.2" !! Swiftest version !> Symbolic name for integrator types character(*), parameter :: UNKNOWN_INTEGRATOR = "UKNOWN INTEGRATOR" diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index e885708df..bce34f3aa 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -128,8 +128,8 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) real(DP), dimension(NDIM) :: L_total_now, L_orbit_now, L_spin_now real(DP) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now, be_cb_now, be_cb_orig, te_now real(DP) :: GMtot_now - character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE_orbit/|E0| = ", ES12.5, ";"' & - //'" DE_total/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)' + character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE_orbit/|E0| = ", ES12.5,' & + //'"; DE_total/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)' associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit) diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index d0700d321..cfd9ef54c 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -191,7 +191,7 @@ module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, ahj(:,:) = 0.0_DP !$omp parallel do default(private) schedule(static)& !$omp shared(npl, nplm, r, Gmass, radius) & - !$omp reduction(+:ahi,j) + !$omp reduction(+:ahi,ahj) do i = 1, nplm #ifdef DOCONLOC do concurrent(j = i+1:npl) shared(i,r,radius,ahi,ahj,Gmass) local(rx,ry,rz,rji2,rlim2) diff --git a/version.txt b/version.txt index 3bae60812..d0d0c6122 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2023.10.0 \ No newline at end of file +2023.10.2 \ No newline at end of file