Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Merge branch '4-write-a-basic-test-simulation-and-user-guide-for-sphe…
Browse files Browse the repository at this point in the history 
…rical-harmonics-use' of github.com:MintonGroup/swiftest into 4-write-a-basic-test-simulation-and-user-guide-for-spherical-harmonics-use
  • Loading branch information
@anand43
anand43 committed Mar 1, 2024
2 parents 1b0804f + d7393ed commit cc3f838
Show file tree
Hide file tree
Showing 4 changed files with 257 additions and 151 deletions.
65 changes: 61 additions & 4 deletions docs/user-guide/basic-simulation/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,12 @@ file. The default value is 10, which means that Swiftest will store 10 outputs i
Setting ``dump_cadence`` to 0 is a a special case that tells Swiftest to store *all* output in memory until the end of the
simulation. This is useful for short simulations, but can be memory intensive for long simulations.

.. note::
Changing the value of ``dump_cadence`` does not change the amount of data that is output by the end of the simulation. It only
changes how often the data is written to disk. Changing the value of ``tstep_out`` (or ``istep_out``) *does* change the amount of
data that is output by the end of the simulation. Intermediate steps between output steps are not saved to disk, and cannot be
recovered after the simulation has finished.

The choice of what values to set for ``tstep_out`` (or ``istep_out``), ``nstep_out``, and ``dump_cadence`` depends on the particular
simulation. Higher values of ``dump_cadence`` are typically useful for simulations with small numbers of bodies and small values
of ```tstep_out`` where frequent writing to disk can severely impact performance. For simulations with large numbers of bodies and
Expand All @@ -87,10 +93,10 @@ so that the memory usage does not become too large. The default value of ``dump_
caes.

We can set these simulation parameters using the :func:`set_parameter <swiftest.Simulation.set_parameter>` method.
Here we have a simulation that runs for 1 My a step size of 0.01 y. We will also save the system every 1000 y and wait until the end
Here we have a simulation that runs for 100,000 y a step size of 0.01 y. We will also save the system every 1000 y and wait until the end
of the simulation to write the simulation data to file using the ``dump_cadence=0`` argument::

sim.set_parameter(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
sim.set_parameter(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)

Once everything is set up, we call the :func:`run <swiftest.Simulation.run>` method to integrate the system forward in time::

Expand All @@ -99,7 +105,7 @@ Once everything is set up, we call the :func:`run <swiftest.Simulation.run>` met
Swiftest is relatively flexible with arguments. You can pass the parameters in when initializing the simulation object, or even later when running.
So the following are all equivalent::

sim = swiftest.Simulation(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
sim = swiftest.Simulation(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)
sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
sim.run()

Expand All @@ -110,7 +116,7 @@ So the following are all equivalent::

sim = swiftest.Simulation()
sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
sim.run(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
sim.run(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)

.. note::
Swiftest uses OpenMP parallelization to help speed up the integration, however the parallelization is most effective when there
Expand All @@ -131,6 +137,57 @@ Once a simulation has been run, its output data is stored in the ``./simdata`` d
default name of ``data.nc``, which is a netCDF file. It is read in and stored as an
`Xarray Dataset <https://docs.xarray.dev/en/stable/>`__ object in the ``sim.data`` attribute.

Here is an example of what the dataset looks like after the above simulation has been run::

In [8]: sim.data
Out[8]:
<xarray.Dataset> Size: 229kB

Dimensions: (time: 101, space: 3, name: 9)
Coordinates:
* time (time) float64 808B 0.0 1e+03 2e+03 ... 9.9e+04 1e+05
* space (space) <U1 12B 'x' 'y' 'z'
* name (name) <U32 1kB 'Sun' 'Mercury' ... 'Uranus' 'Neptune'
Data variables:
id (name) int64 72B 0 1 2 3 4 5 6 7 8
status (time, name) int64 7kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
npl (time) int64 808B 8 8 8 8 8 8 8 8 8 8 ... 8 8 8 8 8 8 8 8 8
ntp (time) int64 808B 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
nplm (time) int64 808B 8 8 8 8 8 8 8 8 8 8 ... 8 8 8 8 8 8 8 8 8
particle_type (name) <U32 1kB 'Central Body' ... 'Massive Body'
rh (time, name, space) float64 22kB nan nan ... -15.76 -0.4033
vh (time, name, space) float64 22kB nan nan ... -0.03418
gr_pseudo_vh (time, name, space) float64 22kB nan nan ... -0.03418
a (time, name) float64 7kB nan 0.3871 0.7233 ... 19.16 30.17
e (time, name) float64 7kB nan 0.2056 ... 0.03379 0.008614
inc (time, name) float64 7kB nan 7.003 3.394 ... 1.306 1.884
capom (time, name) float64 7kB nan 48.3 76.6 ... 150.3 124.7
omega (time, name) float64 7kB nan 29.2 54.96 ... 135.7 310.1
capm (time, name) float64 7kB nan 338.3 200.5 ... 212.3 254.7
varpi (time, name) float64 7kB nan 77.5 131.6 ... 286.0 74.82
lam (time, name) float64 7kB nan 55.84 332.0 ... 138.3 329.5
f (time, name) float64 7kB nan 327.0 200.2 ... 210.3 253.8
cape (time, name) float64 7kB nan 333.0 200.3 ... 211.3 254.3
Gmass (time, name) float64 7kB 39.48 6.554e-06 ... 0.002033
mass (time, name) float64 7kB 1.0 1.66e-07 ... 5.15e-05
rhill (time, name) float64 7kB nan nan nan nan ... nan nan nan
radius (time, name) float64 7kB 0.00465 1.631e-05 ... 0.0001646
origin_time (name) float64 72B 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
origin_type (name) <U32 1kB 'Initial conditions' ... 'Initial condit...
origin_rh (name, space) float64 216B 0.0 0.0 0.0 ... 2.045 -0.7287
origin_vh (name, space) float64 216B 0.0 0.0 0.0 ... 1.149 -0.02168
collision_id (name) int64 72B 0 0 0 0 0 0 0 0 0
discard_time (name) float64 72B 0.0 1.798e+308 ... 1.798e+308 1.798e+308
discard_rh (name, space) float64 216B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
discard_vh (name, space) float64 216B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
discard_body_id (name) int64 72B -2147483647 -2147483647 ... -2147483647
Ip (time, name, space) float64 22kB 0.0 0.0 0.07 ... 0.0 0.23
rot (time, name, space) float64 22kB 642.2 ... 1.721e+05
rotphase (time) float64 808B 0.0 281.4 197.1 ... 176.4 126.8 344.4
j2rp2 (time) float64 808B 4.754e-12 4.754e-12 ... 4.754e-12
j4rp4 (time) float64 808B -2.247e-18 -2.247e-18 ... -2.247e-18



.. toctree::
.. :maxdepth: 2
Expand Down
46 changes: 27 additions & 19 deletions src/netcdf_io/netcdf_io_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,25 +64,6 @@ module netcdf_io
!! ID for the space variable
character(len=1), dimension(3) :: space_coords = ["x","y","z"]
!! The space dimension coordinate labels
character(NAMELEN) :: sign_dimname = "sign"
!! name of the sign dimension for c_lm
integer(I4B) :: sign_dimid
!! ID for sign dimension
integer(I4B) :: sign_varid
!! ID for sign variable
character(NAMELEN) :: l_dimname = "l"
!! name of l dimension for c_lm
integer(I4B) :: l_dimid
!! ID for the l dimension for c_lm
integer(I4B) :: l_varid
!! ID for the l variable
character(NAMELEN) :: m_dimname = "m"
!! name of m dimension for c_lm
integer(I4B) :: m_dimid
!! ID for the m dimension for c_lm
integer(I4B) :: m_varid
!! ID for the m variable


! Non-dimension ids and variable names
character(NAMELEN) :: id_varname = "id"
Expand Down Expand Up @@ -291,8 +272,35 @@ module netcdf_io
!! ID for the id of the other body involved in the discard
logical :: lpseudo_vel_exists = .false.
!! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file.

! Gravitational harmonics ids and variable names
logical :: lc_lm_exists = .false.
!! Logical flag to indicate whether or not the c_lm array was present in an old file.
character(NAMELEN) :: sign_dimname = "sign"
!! name of the sign dimension for c_lm
integer(I4B) :: sign_dimid
!! ID for sign dimension
integer(I4B) :: sign_varid
!! ID for sign variable
integer(I4B), dimension(2) :: sign_coords = [-1,1]
!! The sign dimension coordinate labels
character(NAMELEN) :: l_dimname = "l"
!! name of l dimension for c_lm
integer(I4B) :: l_dimid
!! ID for the l dimension for c_lm
integer(I4B) :: l_varid
!! ID for the l variable
character(NAMELEN) :: m_dimname = "m"
!! name of m dimension for c_lm
integer(I4B) :: m_dimid
!! ID for the m dimension for c_lm
integer(I4B) :: m_varid
!! ID for the m variable
integer(I4B) :: m_dim_max
!! Maximum value of the m dimension
integer(I4B) :: l_dim_max
!! Maximum value of the l dimension

contains
procedure :: close => netcdf_io_close
!! Closes an open NetCDF file
Expand Down
78 changes: 52 additions & 26 deletions src/swiftest/swiftest_gr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,12 @@ pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, para
!! Adapted from David A. Minton's Swifter routine routine gr_kick_getaccb_ns.f90
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system
!! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
real(DP) :: rmag, rdotv, vmag2
integer(I4B) :: i
Expand Down Expand Up @@ -62,12 +65,18 @@ pure module subroutine swiftest_gr_kick_getacch(mu, x, lmask, n, inv_c2, agr)
!! Adapted from David A. Minton's Swifter routine routine gr_whm_kick_getacch.f90
implicit none
! Arguments
real(DP), dimension(:), intent(in) :: mu !! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x !! Position vectors
logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n !! Total number of bodies
real(DP), intent(in) :: inv_c2 !! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr !! Accelerations
real(DP), dimension(:), intent(in) :: mu
!! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x
!! Position vectors
logical, dimension(:), intent(in) :: lmask
!! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n
!! Total number of bodies
real(DP), intent(in) :: inv_c2
!! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr
!! Accelerations
! Internals
integer(I4B) :: i
real(DP) :: beta, rjmag4
Expand Down Expand Up @@ -100,10 +109,14 @@ pure elemental module subroutine swiftest_gr_p4_pos_kick(inv_c2, rx, ry, rz, vx,
!! Adapted from David A. Minton's Swifter routine gr_whm_p4.f90
implicit none
! Arguments
real(DP), intent(in) :: inv_c2 !! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz !! Position vector
real(DP), intent(in) :: vx, vy, vz !! Velocity vector
real(DP), intent(in) :: dt !! Step size
real(DP), intent(in) :: inv_c2
!! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz
!! Position vector
real(DP), intent(in) :: vx, vy, vz
!! Velocity vector
real(DP), intent(in) :: dt
!! Step size
! Internals
real(DP) :: drx, dry, drz
real(DP) :: vmag2
Expand Down Expand Up @@ -133,11 +146,16 @@ pure module subroutine swiftest_gr_pseudovel2vel(param, mu, rh, pv, vh)
!! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh !! Heliocentric velocity vector
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv
!! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh
!! Heliocentric velocity vector
! Internals
real(DP) :: vmag2, rmag, grterm

Expand Down Expand Up @@ -191,11 +209,16 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
!! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh !! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh
!! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv
!! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
! Internals
real(DP) :: v2, G, pv2, rterm, det
real(DP), dimension(NDIM,NDIM) :: J,Jinv
Expand Down Expand Up @@ -264,11 +287,14 @@ pure module subroutine swiftest_gr_vh2pv_body(self, param)
!! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest particle object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest particle object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
integer(I4B) :: i
real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion
real(DP), dimension(:,:), allocatable :: pv
!! Temporary holder of pseudovelocity for in-place conversion

associate(n => self%nbody)
if (n == 0) return
Expand Down
Loading

0 comments on commit cc3f838

Please sign in to comment.