diff --git a/docs/_static/basic_simulation_a_vs_t_plot.png b/docs/_static/basic_simulation_a_vs_t_plot.png new file mode 100644 index 000000000..6fec4c451 Binary files /dev/null and b/docs/_static/basic_simulation_a_vs_t_plot.png differ diff --git a/docs/user-guide/basic-simulation/index.rst b/docs/user-guide/basic-simulation/index.rst index 5958406e7..32bd6ffaa 100644 --- a/docs/user-guide/basic-simulation/index.rst +++ b/docs/user-guide/basic-simulation/index.rst @@ -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 @@ -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 ` 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 ` method to integrate the system forward in time:: @@ -99,7 +105,7 @@ Once everything is set up, we call the :func:`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() @@ -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 @@ -131,7 +137,69 @@ 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 `__ object in the ``sim.data`` attribute. +Here is an example of what the dataset looks like after the above simulation has been run:: + + In [5]: sim.data + Out[5]: + 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) `__. section of the Xarray documentation. Xarray +Datasets are very powerful and flexible, and can be used to analyze and visualize the simulation data in a variety of ways. +Here is an example where we can generate a simple plot of the semimajor axis vs. time history of all the planets in the system:: + + sim.data['a'].where(sim.data.particle_type != 'Central Body', drop=True).plot(x='time',hue='name') + +.. image:: ../../_static/basic_simulation_a_vs_t_plot.png + +This is just a simple example of what you can do with the simulation data. Xarray has a large number of built-in plotting and +data processing functions. For more information, see the `Xarray documentation `__. -.. toctree:: .. :maxdepth: 2 .. :hidden: diff --git a/docs/user-guide/gravitational-harmonics/index.rst b/docs/user-guide/gravitational-harmonics/index.rst new file mode 100644 index 000000000..f4452666c --- /dev/null +++ b/docs/user-guide/gravitational-harmonics/index.rst @@ -0,0 +1,116 @@ +########################## +Gravitational Harmonics +########################## + +Here, we show how to use Swiftest's Gravitational Harmonics capabilities. This is based on ``/spherical_harmonics_cb`` +in ``swiftest/examples``. Swiftest uses `SHTOOLS `__ to compute the gravitational +harmonics coefficients for a given body and calculate it's associated acceleration kick. The conventions and formulae used +to calculate the additional kick are described `here `__. The gravitational +potential is given by the following equation: + +.. math:: + + U(r) = \frac{GM}{r} \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \left( \frac{R_0}{r} \right)^l C_{lm} Y_{lm} (\theta, \phi) \label{grav_pot} + +Gravitational potential :math:`U` at a point :math:`\vec{r}`; :math:`\theta` is the polar angle; :math:`\phi` is the azimuthal angle; +:math:`R_0` is the central body radius; :math:`G` is the gravitational constant; :math:`Y_{lm}` is the spherical harmonic function at +degree :math:`l` and order :math:`m`; :math:`C_{lm}` is the corresponding coefficient. + +Gravitational Harmonics coefficients +===================================== + +Swiftest adopts the :math:`4\pi` or geodesy normalization for the gravitational harmonics coefficients described +in `Weiczorek et al. (2015) `__. + +The coefficients can be computed in a number of ways: + +- Using the axes measurements of the body. (:func:`clm_from_ellipsoid `) + +- Using a surface relief grid (:func:`clm_from_relief `). *Note: This function is still in development and may not work as expected.* + +- Manually entering the coefficients when adding the central body. (:func:`add_body `) + +Set up a Simulation +==================== + +Let's start with setting up the simulation object with units of `km`, `days`, and `kg`. :: + + import swiftest + + sim = swiftest.Simulation(DU2M = 1e3, TU = 'd', MU = 'kg', integrator = 'symba') + sim.clean() + +Computing coefficients from axes measurements +=============================================== + +Given the axes measurements of a body, the gravitational harmonics coefficients can be computed in a straightforward +manner. Here we use Chariklo as an example body and refer to Jacobi Ellipsoid model from +`Leiva et al. (2017) `__ for the axes measurements. :: + + # Define the central body parameters. + cb_mass = 6.1e18 + cb_radius = 123 + cb_a = 157 + cb_b = 139 + cb_c = 86 + cb_volume = 4.0 / 3 * np.pi * cb_radius**3 + cb_density = cb_mass / cb_volume + cb_T_rotation = 7.004 # hours + cb_T_rotation/= 24.0 # converting to julian days (TU) + cb_rot = [[0, 0, 360.0 / cb_T_rotation]] # degrees/TU + +Once the central body parameters are defined, we can compute the gravitational harmonics coefficients (:math:`C_{lm}`). +The output coefficients are already correctly normalized. :: + + c_lm, cb_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lmax = 6, lref_radius = True) + +Add the central body to the simulation along with the coefficients. :: + + sim.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, c_lm = c_lm) + +Now the user can set up the rest of the simulation as they prefer. + +Final Steps for Running the Simulation +======================================= + +Add other bodies to the simulation. :: + + # Add user-defined massive bodies + npl = 5 + density_pl = cb_density + + name_pl = ["SemiBody_01", "SemiBody_02", "SemiBody_03", "SemiBody_04", "SemiBody_05"] + a_pl = rng.uniform(250, 400, npl) + e_pl = rng.uniform(0.0, 0.05, npl) + inc_pl = rng.uniform(0.0, 10, npl) + capom_pl = rng.uniform(0.0, 360.0, npl) + omega_pl = rng.uniform(0.0, 360.0, npl) + capm_pl = rng.uniform(0.0, 360.0, npl) + R_pl = np.array([0.5, 1.0, 1.2, 0.75, 0.8]) + M_pl = 4.0 / 3 * np.pi * R_pl**3 * density_pl + Ip_pl = np.full((npl,3),0.4,) + rot_pl = np.zeros((npl,3)) + mtiny = 1.1 * np.max(M_pl) + + sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, mass=M_pl, radius=R_pl, Ip=Ip_pl, rot=rot_pl) + +Set the parameters for the simulation and run. :: + + sim.set_parameter(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny) + + # Run the simulation + sim.run() + +Setting a reference radius for the coefficients +================================================== + +The coefficients can be computed with respect to a reference radius. This is useful when the user wants to explicitly set the reference radius. :: + + c_lm, cb_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lmax = 6, lref_radius = True, ref_radius = cb_radius) + + + + +.. .. toctree:: +.. :maxdepth: 2 +.. :hidden: \ No newline at end of file diff --git a/docs/user-guide/index.rst b/docs/user-guide/index.rst index f3a941582..faa81fdc3 100644 --- a/docs/user-guide/index.rst +++ b/docs/user-guide/index.rst @@ -8,7 +8,7 @@ In this user guide, you will find detailed descriptions and examples that descri - A more in-depth :doc:`Detailed Simulation Setup ` -- Understanding the :doc:`Spherical Harmonics capabilities ` +- Understanding the :doc:`Gravitational Harmonics capabilities ` .. toctree:: :maxdepth: 2 @@ -16,4 +16,4 @@ In this user guide, you will find detailed descriptions and examples that descri Basic Simulation Detailed Simulation Setup - Spherical Harmonics \ No newline at end of file + Gravitational Harmonics \ No newline at end of file diff --git a/docs/user-guide/spherical-harmonics/index.rst b/docs/user-guide/spherical-harmonics/index.rst deleted file mode 100644 index 06380bd4a..000000000 --- a/docs/user-guide/spherical-harmonics/index.rst +++ /dev/null @@ -1,12 +0,0 @@ -################### -Spherical Harmonics -################### - -Here, we show how to use Swiftest's Spherical Harmonics capabilities. -This is based on ``/spherical_harmonics_cb`` in ``swiftest/examples``. - - - -.. .. toctree:: -.. :maxdepth: 2 -.. :hidden: \ No newline at end of file