Cleanup and more docs

doc/src/physics.rst

@@ -340,7 +340,7 @@ For example, if there is only a single, stationary, particle in the domain or if
 The ``frame`` parameter is typically set to its default value of ``global``.
 This means that the REBOUND simulation is in the global frame of reference.
 Note that this can be a different frame than the |code| frame if the rotating frame is active.
-Setting ``frame=local`` means that |code| artemis will not apply corrections to the REBOUND frame, but caution should be used with this option.
+Setting ``frame=local`` means that |code| will not apply corrections to the REBOUND frame, but caution should be used with this option.
 
 |code| can also make full use of the collision functionality in REBOUND.
 Particles can be given radii defining their collisional cross-section.
@@ -404,13 +404,22 @@ One example is:
    beta  = 0.0     # = torque-free mass removal
 
 There are three softening models available: ``none`` does nothing, ``spline`` is the spline softening used in Gadget (Springel 2001), and ``plummer`` modifies :math:`r^2 \rightarrow r^2 + r_s^2`.
-Spline softening is exactly Kepelerian outside the softening radius, whereas Plummer softening asymptoically approaches Keplerian outside the softening radius.
+Spline softening is exactly Keplerian outside the softening radius, whereas Plummer softening asymptoically approaches Keplerian outside the softening radius.
 Roughly, the Plummer softening radius is :math:`\sim 2.8 \times` the spline softening radius.
 
 Mass accretion is handled by a particle's ``<../sink>`` attribute.
-Two types of removal rates are avaible.
+|code| implements mass accretion by solving:
+
+.. math::
+   \frac{\partial \rho}{\partial t} = - \gamma \rho \\
+   \frac{\partial (\rho v_r)}{\partial t} = - \gamma \rho v_r \\
+   \frac{\partial (\rho v_\theta)}{\partial t} = - \gamma \rho v_\theta \\
+   \frac{\partial (\rho v_\phi)}{\partial t} = - \beta \rho v_\phi \\
+where :math:`(r,\theta,\phi)` refers to a spherical coordinate system coordinate system centered on the particle. 
+
+Two types of removal rates are available.
 The first, ``gamma``, sets the rate of mass removal for cells inside the particle's sink radius.
-The second, ``beta``, controls the amount of momentum removed from the fluid.
+The second, ``beta``, controls the amount of angular momentum removed from the fluid.
 Typically ``beta`` will either be zero -- corresponding to angular-momentum conserving mass removal -- or equal to ``gamma`` -- corresponding to isotropic momentum removal.
 
 
@@ -615,8 +624,45 @@ All planets share the options set in the ``<planet>`` block.
 Note that the central object has been left unspecified.
 An additional particle can be added at the center of the system to represent the star, or a ``<binary>`` block could be added to model a circumbinary planetary sytem, or even a ``<system>`` block for a circum-cluster planetary system.
 
+N-Body Outputs 
+^^^^^^^^^^^^^^
+
+There are two types of outputs that the ``<nbody>`` will produce if ``dt_output`` is defined.
+The first output is the ``.reb`` file. 
+This file lists the state of each particle in the simulation at the output time, as well as the total amount of mass and momentum change since the last output.
+It is laid out as an ascii history file. 
+The header for a 2 particle system looks like
+
+::
+
+  # job.reb
+  # NBody data N = 2
+  # [1]=time  [2]=hash  [3]=active  [4]=mass  [5]=x [6]=y [7]=z [8]=vx  [9]=vy  [10]=vz [11]=dm [12]=dmx_g  [13]=dmy_g  [14]=dmz_g  [15]=dmx_a  [16]=dmy_a  [17]=dmz_a
+
+Note that everything is in Cartesian coordinates -- regardless of the problem geometry. 
+Columns 11-17 hold cumulative quantities (since the last output) of mass accreted (``dm``), momentum change due to the gravitational interaction with the fluids (``dm(x,y,z)_g``), and momentum accreted (``dm(x,y,z)_a``). 
+These quantities are calculated using every timestep of the simulation. 
+Because of this, very accurate time-averages can be obtained from this data. 
+Each particle is written as a new line. 
+Therefore, for this 2 particle example, the first two rows of data correspond to the same time. 
+
+
+The second type of output file is the ``.orb`` file. 
+This file takes all the data that is already in the ``.reb`` file and outputs derived quantities useful for analysis of binaries, e.g., orbital elements, differential forces, and average forces.
+An example header for the ``.orb`` file is:
+
+::
 
+  # job.orb.0_1
+  # NBody Orbit data
+  # [1]=time  [2]=mb  [3]=xc  [4]=yc  [5]=zc  [6]=xb  [7]=yb  [8]=zb  [9]=vxc [10]=vyc  [11]=vzc  [12]=vxb  [13]=vyb  [14]=vzb  [15]=qb [16]=nb [17]=ab [18]=eb [19]=Ib [20]=o  [21]=O  [22]=pomega [23]=f  [24]=h  [25]=ex [26]=ey [27]=ix [28]=iy [29]=dm [30]=Fx_grav_com  [31]=Fy_grav_com  [32]=Fz_grav_com  [33]=Fx_acc_com [34]=Fy_acc_com [35]=Fz_acc_com [36]=Fx_grav_bin  [37]=Fy_grav_bin  [38]=Fz_grav_bin  [39]=Fx_acc_bin [40]=Fy_acc_bin [41]=Fz_acc_bin
+
+At every output time, |code| looks for all pairs of bound particles. Each bound pair is output to their own file, e.g., particles 0 and 1 are put in ``.orb.0_1``. 
+Columns 3-14 provide the center of mass position/velocity and separation position/velocity. 
+Columns 15-28 have the orbital elements (in addition to the mass ratio and mean motion) computed directly from REBOUND. 
+Columns 29-41 combine Columns 11-17 of the ``.reb`` file into forces on the center of mass or the separation. In particular, forces marked with ``bin`` are combined as :math:`\mu_1 f_2 - \mu_2 f_1`, where :math:`\mu_i = m_i/m_b`. 
 
+To disable these outputs, set ``disable_outputs = true`` in the ``<nbody>`` block. 
 
 External Gravity
 ----------------

doc/src/running.rst

@@ -78,6 +78,7 @@ Typically, there are three types of outputs specified here, history, hdf5 snapsh
 
 The next blocks define the simulation mesh dimensions, boundary conditions, and meshblock size.
 This example sets up a 2D cylindrical mesh that spans the full :math:`2 \pi` in azimuth.
+
 ::
 
   <parthenon/mesh>
@@ -161,8 +162,8 @@ For more details see the :ref:`physics` and :ref:`parameters` sections
   alpha = 1e-3
 
   <gravity>
-  type = binary
   gm = 1.0
+  <gravity/binary>
   q = 1e-3
   a = 1.0
   sft2 = .06
@@ -191,7 +192,7 @@ Run |code|
 The exact command to launch it depends on the system it is run on.
 This example will assume a SLURM-like cluster.
 
-To launch a fresh |code| on ``$NPROCS`` CPUs with ``srun``,
+To launch a fresh |code| simulation on ``$NPROCS`` CPUs with ``srun``,
 
 ::
 
@@ -205,14 +206,23 @@ To restart a previous run, use the ``-r`` argument
 
 A modified input file can optionally still be passed with the ``-i`` argument.
 
+When launching |code| on GPUs with ``srun``, the application passed to ``srun`` should not be |code|, but instead a script that 
+determines how to bind the CPU cores to the available GPUs. 
+One such script that comes with ``Kokkos`` is ``hpcbind``. 
+Assuming the |code| source code lives in ``$ARTEMIS_HOME``, the following will launch |code| on the available number of GPUs
+
+::
+
+  srun -n $NPROCS $ARTEMIS_HOME/external/parthenon/external/Kokkos/bin/hpcbind -- artemis -i input.par
+
 Return codes
 ^^^^^^^^^^^^
 
 When using batch submissions, it is possible to set up a self-restarting job.
 The easiest way to do this is to take advantage of SLURM interrupt signals and the |code| return code.
 |code|
 
-An example batch submission script, ``run.sh``, would look like:
+An example CPU batch submission script, ``run.sh``, would look like:
 
 ::
 

inputs/disk/binary_cyl.in

@@ -75,7 +75,6 @@ riemann = hllc
 gamma = 1.00001
 dfloor = 1.0e-10
 siefloor = 1.0e-10
-de_switch = 0.0
 refine_field = pressure
 refine_type = gradient
 refine_thr = 3.0

inputs/disk/disk_axi.in

@@ -72,7 +72,6 @@ reconstruct = plm
 riemann = hlle
 dfloor = 1.0e-10
 pfloor = 1.0e-13
-de_switch = 1.0
 
 <rotating_frame>
 omega = 1.0

inputs/disk/disk_cart.in

@@ -96,7 +96,6 @@ reconstruct = plm
 riemann = hllc
 dfloor = 1.0e-10
 siefloor = 1.0e-15
-de_switch = 1.0
 
 <gas/viscosity>
 type = alpha

inputs/disk/disk_cyl.in

@@ -73,7 +73,6 @@ reconstruct = plm
 riemann = hllc
 dfloor = 1.0e-10
 siefloor = 1.0e-10
-de_switch = 1.0
 
 <rotating_frame>
 omega = 1.0

inputs/disk/disk_sph.in

@@ -72,7 +72,6 @@ reconstruct = plm
 riemann = hlle
 dfloor = 1.0e-10
 pfloor = 1.0e-13
-de_switch = 1.0
 
 <rotating_frame>
 omega = 1.0

src/geometry/axisymmetric.hpp

@@ -59,12 +59,11 @@ class Coords<Coordinates::axisymmetric>
 
   KOKKOS_INLINE_FUNCTION std::array<Real, 3> FaceCenX2(const CellFace f) {
     // <r> = d(r^3/3) / d(r^2/2)
-    std::array<Real, 3> xf{2.0 / 3.0 *
-                               (bnds.x1[0] * bnds.x1[0] + bnds.x1[0] * bnds.x1[1] +
-                                bnds.x1[1] * bnds.x1[1]) /
-                               (bnds.x1[0] + bnds.x1[1]),
-                           bnds.x2[static_cast<int>(f)], 0.5 * (bnds.x3[0] + bnds.x3[1])};
-    return xf;
+    return {2.0 / 3.0 *
+                (bnds.x1[0] * bnds.x1[0] + bnds.x1[0] * bnds.x1[1] +
+                 bnds.x1[1] * bnds.x1[1]) /
+                (bnds.x1[0] + bnds.x1[1]),
+            bnds.x2[static_cast<int>(f)], 0.5 * (bnds.x3[0] + bnds.x3[1])};
   }
 
   KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) {
@@ -100,8 +99,7 @@ class Coords<Coordinates::axisymmetric>
   ConvertCoordsToCart(const std::array<Real, 3> &xi) {
     const Real cp = std::cos(xi[2]);
     const Real sp = std::sin(xi[2]);
-    std::array<Real, 3> xo{xi[0] * cp, xi[0] * sp, xi[1]};
-    return xo;
+    return {xi[0] * cp, xi[0] * sp, xi[1]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCart(const std::array<Real, 3> &xi) {
     const Real cp = std::cos(xi[2]);
@@ -119,8 +117,7 @@ class Coords<Coordinates::axisymmetric>
     const Real R = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]);
     const Real ct = xi[1] / (R + Fuzz<Real>());
     const Real st = xi[0] / (R + Fuzz<Real>());
-    std::array<Real, 3> xo{R, std::acos(ct), xi[2]};
-    return xo;
+    return {R, std::acos(ct), xi[2]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToSph(const std::array<Real, 3> &xi) {
     const Real rsph = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]);
@@ -136,8 +133,7 @@ class Coords<Coordinates::axisymmetric>
 
   KOKKOS_INLINE_FUNCTION std::array<Real, 3>
   ConvertCoordsToCyl(const std::array<Real, 3> &xi) {
-    std::array<Real, 3> xo{xi[0], xi[2], xi[1]};
-    return xo;
+    return {xi[0], xi[2], xi[1]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array<Real, 3> &xi) {
     // clang-format off

src/geometry/cylindrical.hpp

@@ -56,12 +56,11 @@ class Coords<Coordinates::cylindrical>
 
   KOKKOS_INLINE_FUNCTION std::array<Real, 3> FaceCenX3(const CellFace f) {
     //  <r> = d(r^3/3) / d(r^2/2)
-    std::array<Real, 3> xf{2.0 / 3.0 *
-                               (bnds.x1[0] * bnds.x1[0] + bnds.x1[0] * bnds.x1[1] +
-                                bnds.x1[1] * bnds.x1[1]) /
-                               (bnds.x1[0] + bnds.x1[1]),
-                           0.5 * (bnds.x2[0] + bnds.x2[1]), bnds.x3[static_cast<int>(f)]};
-    return xf;
+    return {2.0 / 3.0 *
+                (bnds.x1[0] * bnds.x1[0] + bnds.x1[0] * bnds.x1[1] +
+                 bnds.x1[1] * bnds.x1[1]) /
+                (bnds.x1[0] + bnds.x1[1]),
+            0.5 * (bnds.x2[0] + bnds.x2[1]), bnds.x3[static_cast<int>(f)]};
   }
 
   KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) {
@@ -97,8 +96,7 @@ class Coords<Coordinates::cylindrical>
   ConvertCoordsToCart(const std::array<Real, 3> &xi) {
     const Real cp = std::cos(xi[1]);
     const Real sp = std::sin(xi[1]);
-    std::array<Real, 3> xo{xi[0] * cp, xi[0] * sp, xi[2]};
-    return xo;
+    return {xi[0] * cp, xi[0] * sp, xi[2]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCart(const std::array<Real, 3> &xi) {
     const Real cp = std::cos(xi[1]);
@@ -114,8 +112,7 @@ class Coords<Coordinates::cylindrical>
     const Real R = std::sqrt(xi[0] * xi[0] + xi[2] * xi[2]);
     const Real ct = xi[2] / (R + Fuzz<Real>());
     const Real st = xi[0] / (R + Fuzz<Real>());
-    std::array<Real, 3> xo{R, std::acos(ct), xi[1]};
-    return xo;
+    return {R, std::acos(ct), xi[1]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToSph(const std::array<Real, 3> &xi) {
     const Real rsph = std::sqrt(xi[0] * xi[0] + xi[2] * xi[2]);
@@ -140,8 +137,7 @@ class Coords<Coordinates::cylindrical>
 
   KOKKOS_INLINE_FUNCTION std::array<Real, 3>
   ConvertCoordsToAxi(const std::array<Real, 3> &xi) {
-    std::array<Real, 3> xo{xi[0], xi[2], xi[1]};
-    return xi;
+    return {xi[0], xi[2], xi[1]};
   }
   KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array<Real, 3> &xi) {
     std::array<Real, 3> ex1{1.0, 0.0, 0.0};