diff --git a/doc/src/physics.rst b/doc/src/physics.rst index 7e2b187..3c6b204 100644 --- a/doc/src/physics.rst +++ b/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 ```` 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 ```` block could be added to model a circumbinary planetary sytem, or even a ```` block for a circum-cluster planetary system. +N-Body Outputs +^^^^^^^^^^^^^^ + +There are two types of outputs that the ```` 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 ```` block. External Gravity ---------------- diff --git a/doc/src/running.rst b/doc/src/running.rst index cca63e7..edddc9e 100644 --- a/doc/src/running.rst +++ b/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. + :: @@ -161,8 +162,8 @@ For more details see the :ref:`physics` and :ref:`parameters` sections alpha = 1e-3 - type = binary gm = 1.0 + 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,6 +206,15 @@ 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 ^^^^^^^^^^^^ @@ -212,7 +222,7 @@ 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: :: diff --git a/inputs/disk/binary_cyl.in b/inputs/disk/binary_cyl.in index b530961..7f4bec8 100644 --- a/inputs/disk/binary_cyl.in +++ b/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 diff --git a/inputs/disk/disk_axi.in b/inputs/disk/disk_axi.in index e035ebf..27b01b2 100644 --- a/inputs/disk/disk_axi.in +++ b/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 omega = 1.0 diff --git a/inputs/disk/disk_cart.in b/inputs/disk/disk_cart.in index c587fb4..fcb8c0a 100644 --- a/inputs/disk/disk_cart.in +++ b/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 type = alpha diff --git a/inputs/disk/disk_cyl.in b/inputs/disk/disk_cyl.in index 5aec57d..9cd634a 100644 --- a/inputs/disk/disk_cyl.in +++ b/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 omega = 1.0 diff --git a/inputs/disk/disk_sph.in b/inputs/disk/disk_sph.in index 56dc2a9..dfa17c4 100644 --- a/inputs/disk/disk_sph.in +++ b/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 omega = 1.0 diff --git a/src/geometry/axisymmetric.hpp b/src/geometry/axisymmetric.hpp index b31afe2..c205c44 100644 --- a/src/geometry/axisymmetric.hpp +++ b/src/geometry/axisymmetric.hpp @@ -59,12 +59,11 @@ class Coords KOKKOS_INLINE_FUNCTION std::array FaceCenX2(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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(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(f)], 0.5 * (bnds.x3[0] + bnds.x3[1])}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -100,8 +99,7 @@ class Coords ConvertCoordsToCart(const std::array &xi) { const Real cp = std::cos(xi[2]); const Real sp = std::sin(xi[2]); - std::array 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 &xi) { const Real cp = std::cos(xi[2]); @@ -119,8 +117,7 @@ class Coords const Real R = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]); const Real ct = xi[1] / (R + Fuzz()); const Real st = xi[0] / (R + Fuzz()); - std::array xo{R, std::acos(ct), xi[2]}; - return xo; + return {R, std::acos(ct), xi[2]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToSph(const std::array &xi) { const Real rsph = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]); @@ -136,8 +133,7 @@ class Coords KOKKOS_INLINE_FUNCTION std::array ConvertCoordsToCyl(const std::array &xi) { - std::array xo{xi[0], xi[2], xi[1]}; - return xo; + return {xi[0], xi[2], xi[1]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array &xi) { // clang-format off diff --git a/src/geometry/cylindrical.hpp b/src/geometry/cylindrical.hpp index 4bef543..b624d34 100644 --- a/src/geometry/cylindrical.hpp +++ b/src/geometry/cylindrical.hpp @@ -56,12 +56,11 @@ class Coords KOKKOS_INLINE_FUNCTION std::array FaceCenX3(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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(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(f)]}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -97,8 +96,7 @@ class Coords ConvertCoordsToCart(const std::array &xi) { const Real cp = std::cos(xi[1]); const Real sp = std::sin(xi[1]); - std::array 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 &xi) { const Real cp = std::cos(xi[1]); @@ -114,8 +112,7 @@ class Coords const Real R = std::sqrt(xi[0] * xi[0] + xi[2] * xi[2]); const Real ct = xi[2] / (R + Fuzz()); const Real st = xi[0] / (R + Fuzz()); - std::array xo{R, std::acos(ct), xi[1]}; - return xo; + return {R, std::acos(ct), xi[1]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToSph(const std::array &xi) { const Real rsph = std::sqrt(xi[0] * xi[0] + xi[2] * xi[2]); @@ -140,8 +137,7 @@ class Coords KOKKOS_INLINE_FUNCTION std::array ConvertCoordsToAxi(const std::array &xi) { - std::array xo{xi[0], xi[2], xi[1]}; - return xi; + return {xi[0], xi[2], xi[1]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array &xi) { std::array ex1{1.0, 0.0, 0.0}; diff --git a/src/geometry/geometry.hpp b/src/geometry/geometry.hpp index fc50437..430164f 100644 --- a/src/geometry/geometry.hpp +++ b/src/geometry/geometry.hpp @@ -184,21 +184,18 @@ class CoordsBase { KOKKOS_INLINE_FUNCTION std::array FaceCenX1(const CellFace f) { // The centroid value of the X1 face - std::array xf{bnds.x1[static_cast(f)], static_cast(this)->x2v(), - static_cast(this)->x3v()}; - return xf; + return {bnds.x1[static_cast(f)], static_cast(this)->x2v(), + static_cast(this)->x3v()}; } KOKKOS_INLINE_FUNCTION std::array FaceCenX2(const CellFace f) { // The centroid value of the X2 face - std::array xf{static_cast(this)->x1v(), bnds.x2[static_cast(f)], - static_cast(this)->x3v()}; - return xf; + return {static_cast(this)->x1v(), bnds.x2[static_cast(f)], + static_cast(this)->x3v()}; } KOKKOS_INLINE_FUNCTION std::array FaceCenX3(const CellFace f) { // The centroid value of the X3 face - std::array xf{static_cast(this)->x1v(), static_cast(this)->x2v(), - bnds.x3[static_cast(f)]}; - return xf; + return {static_cast(this)->x1v(), static_cast(this)->x2v(), + bnds.x3[static_cast(f)]}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -269,8 +266,7 @@ class CoordsBase { const Real ct = xi[2] / (r + Fuzz()); const Real cp = xi[0] / (R + Fuzz()); const Real sp = xi[1] / (R + Fuzz()); - std::array xo{r, std::acos(ct), std::atan2(sp, cp)}; - return xo; + return {r, std::acos(ct), std::atan2(sp, cp)}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToSph(const std::array &xi) { // Convert the input vector from the problem coordinate system to spherical @@ -292,8 +288,7 @@ class CoordsBase { Real R = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]); const Real cp = xi[0] / (R + Fuzz()); const Real sp = xi[1] / (R + Fuzz()); - std::array xo{R, std::atan2(sp, cp), xi[2]}; - return xo; + return {R, std::atan2(sp, cp), xi[2]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array &xi) { // Convert the input vector from the problem coordinate system to cylindrical @@ -312,8 +307,7 @@ class CoordsBase { Real R = std::sqrt(xi[0] * xi[0] + xi[1] * xi[1]); const Real cp = xi[0] / (R + Fuzz()); const Real sp = xi[1] / (R + Fuzz()); - std::array xo{R, xi[2], std::atan2(sp, cp)}; - return xo; + return {R, xi[2], std::atan2(sp, cp)}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array &xi) { // Convert the input vector from the problem coordinate system to axisymmetric @@ -354,44 +348,37 @@ class CoordsBase { // Return all cell widths const Real xv[3] = {static_cast(this)->x1v(), static_cast(this)->x2v(), static_cast(this)->x3v()}; - const std::array dx{ - static_cast(this)->hx1(xv[0], xv[1], xv[2]) * (bnds.x1[1] - bnds.x1[0]), - static_cast(this)->hx2(xv[0], xv[1], xv[2]) * (bnds.x2[1] - bnds.x2[0]), - static_cast(this)->hx3(xv[0], xv[1], xv[2]) * (bnds.x3[1] - bnds.x3[0])}; - return dx; + return {static_cast(this)->hx1(xv[0], xv[1], xv[2]) * (bnds.x1[1] - bnds.x1[0]), + static_cast(this)->hx2(xv[0], xv[1], xv[2]) * (bnds.x2[1] - bnds.x2[0]), + static_cast(this)->hx3(xv[0], xv[1], xv[2]) * (bnds.x3[1] - bnds.x3[0])}; } KOKKOS_INLINE_FUNCTION std::array GetCellCenter() { // Get the cell centroid - std::array xv{static_cast(this)->x1v(), static_cast(this)->x2v(), - static_cast(this)->x3v()}; - return xv; + return {static_cast(this)->x1v(), static_cast(this)->x2v(), + static_cast(this)->x3v()}; } KOKKOS_INLINE_FUNCTION std::array GetScaleFactors() { // Get the volume averaged scale factors - std::array hx{static_cast(this)->hx1v(), static_cast(this)->hx2v(), - static_cast(this)->hx3v()}; - return hx; + return {static_cast(this)->hx1v(), static_cast(this)->hx2v(), + static_cast(this)->hx3v()}; } KOKKOS_INLINE_FUNCTION std::array GetFaceAreaX1() { // Get the lower and upper face areas in the X1 direction - const std::array areas{static_cast(this)->AreaX1(bnds.x1[0]), - static_cast(this)->AreaX1(bnds.x1[1])}; - return areas; + return {static_cast(this)->AreaX1(bnds.x1[0]), + static_cast(this)->AreaX1(bnds.x1[1])}; } KOKKOS_INLINE_FUNCTION std::array GetFaceAreaX2() { // Get the lower and upper face areas in the X2 direction - const std::array areas{static_cast(this)->AreaX2(bnds.x2[0]), - static_cast(this)->AreaX2(bnds.x2[1])}; - return areas; + return {static_cast(this)->AreaX2(bnds.x2[0]), + static_cast(this)->AreaX2(bnds.x2[1])}; } KOKKOS_INLINE_FUNCTION std::array GetFaceAreaX3() { // Get the lower and upper face areas in the X3 direction - const std::array areas{static_cast(this)->AreaX3(bnds.x3[0]), - static_cast(this)->AreaX3(bnds.x3[1])}; - return areas; + return {static_cast(this)->AreaX3(bnds.x3[0]), + static_cast(this)->AreaX3(bnds.x3[1])}; } template @@ -422,24 +409,18 @@ class CoordsBase { KOKKOS_INLINE_FUNCTION std::array GetConnX1() { // { dh1/dx1, dh2/dx1, dh3/dx1 } - const std::array Gamma{static_cast(this)->dh1dx1(), - static_cast(this)->dh2dx1(), - static_cast(this)->dh3dx1()}; - return Gamma; + return {static_cast(this)->dh1dx1(), static_cast(this)->dh2dx1(), + static_cast(this)->dh3dx1()}; } KOKKOS_INLINE_FUNCTION std::array GetConnX2() { // { dh1/dx2, dh2/dx2, dh3/dx2 } - const std::array Gamma{static_cast(this)->dh1dx2(), - static_cast(this)->dh2dx2(), - static_cast(this)->dh3dx2()}; - return Gamma; + return {static_cast(this)->dh1dx2(), static_cast(this)->dh2dx2(), + static_cast(this)->dh3dx2()}; } KOKKOS_INLINE_FUNCTION std::array GetConnX3() { // { dh1/dx3, dh2/dx3, dh3/dx3 } - const std::array Gamma{static_cast(this)->dh1dx3(), - static_cast(this)->dh2dx3(), - static_cast(this)->dh3dx3()}; - return Gamma; + return {static_cast(this)->dh1dx3(), static_cast(this)->dh2dx3(), + static_cast(this)->dh3dx3()}; } KOKKOS_INLINE_FUNCTION Mat3x3 GetConns() { diff --git a/src/geometry/spherical.hpp b/src/geometry/spherical.hpp index b6fac32..7d8e802 100644 --- a/src/geometry/spherical.hpp +++ b/src/geometry/spherical.hpp @@ -87,22 +87,20 @@ class Coords KOKKOS_INLINE_FUNCTION std::array FaceCenX2(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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(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(f)], 0.5 * (bnds.x3[0] + bnds.x3[1])}; } KOKKOS_INLINE_FUNCTION std::array FaceCenX3(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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(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(f)]}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -176,8 +174,7 @@ class Coords const Real sp = std::sin(xi[2]); const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; - return xo; + return {xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCart(const std::array &xi) { const Real cp = std::cos(xi[2]); @@ -209,8 +206,7 @@ class Coords const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st, xi[2], xi[0] * ct}; - return xo; + return {xi[0] * st, xi[2], xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array &xi) { // rhat = st * Rhat + ct * Zhat = (st, 0, ct) @@ -228,8 +224,7 @@ class Coords ConvertCoordsToAxi(const std::array &xi) { const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st, xi[0] * ct, xi[2]}; - return xo; + return {xi[0] * st, xi[0] * ct, xi[2]}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array &xi) { const Real ct = std::cos(xi[1]); @@ -296,22 +291,20 @@ class Coords KOKKOS_INLINE_FUNCTION std::array FaceCenX2(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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(f)], 0.0}; - 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(f)], 0.0}; } KOKKOS_INLINE_FUNCTION std::array FaceCenX3(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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]), 0.0}; - 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]), 0.0}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -382,8 +375,7 @@ class Coords const Real sp = 0.0; const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; - return xo; + return {xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCart(const std::array &xi) { const Real cp = 1.0; @@ -415,8 +407,7 @@ class Coords const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st, 0.0, xi[0] * ct}; - return xo; + return {xi[0] * st, 0.0, xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array &xi) { // rhat = st * Rhat + ct * Zhat = (st, 0, ct) @@ -434,8 +425,7 @@ class Coords ConvertCoordsToAxi(const std::array &xi) { const Real ct = std::cos(xi[1]); const Real st = std::sin(xi[1]); - std::array xo{xi[0] * st, xi[0] * ct, 0.0}; - return xo; + return {xi[0] * st, xi[0] * ct, 0.0}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array &xi) { const Real ct = std::cos(xi[1]); @@ -472,22 +462,20 @@ class Coords KOKKOS_INLINE_FUNCTION std::array FaceCenX2(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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]), - M_PI * 0.5, 0.0}; - 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]), + M_PI * 0.5, 0.0}; } KOKKOS_INLINE_FUNCTION std::array FaceCenX3(const CellFace f) { // = d(r^3/3) / d(r^2/2) - std::array 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]), - M_PI * 0.5, 0.0}; - 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]), + M_PI * 0.5, 0.0}; } KOKKOS_INLINE_FUNCTION Real AreaX1(const Real x1f) { @@ -543,8 +531,7 @@ class Coords const Real sp = 0.0; const Real ct = 0.0; const Real st = 1.0; - std::array xo{xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; - return xo; + return {xi[0] * st * cp, xi[0] * st * sp, xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCart(const std::array &xi) { const Real cp = 1.0; @@ -576,8 +563,7 @@ class Coords const Real ct = 0.0; const Real st = 1.0; - std::array xo{xi[0] * st, 0.0, xi[0] * ct}; - return xo; + return {xi[0] * st, 0.0, xi[0] * ct}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToCyl(const std::array &xi) { // rhat = st * Rhat + ct * Zhat = (st, 0, ct) @@ -595,8 +581,7 @@ class Coords ConvertCoordsToAxi(const std::array &xi) { const Real ct = 0.0; const Real st = 1.0; - std::array xo{xi[0] * st, xi[0] * ct, 0.0}; - return xo; + return {xi[0] * st, xi[0] * ct, 0.0}; } KOKKOS_INLINE_FUNCTION Mat3x3 ConvertVecToAxi(const std::array &xi) { const Real ct = 0.0; diff --git a/src/rotating_frame/rotating_frame.hpp b/src/rotating_frame/rotating_frame.hpp index 44ae672..9e7daa9 100644 --- a/src/rotating_frame/rotating_frame.hpp +++ b/src/rotating_frame/rotating_frame.hpp @@ -40,8 +40,7 @@ KOKKOS_INLINE_FUNCTION std::array RotationVelocity(const std::array vrot{ex1[1] * vp, ex2[1] * vp, ex3[1] * vp}; - return vrot; + return {ex1[1] * vp, ex2[1] * vp, ex3[1] * vp}; } } // namespace RotatingFrame diff --git a/tst/scripts/disk/disk.py b/tst/scripts/disk/disk.py index 1f97404..3db7622 100644 --- a/tst/scripts/disk/disk.py +++ b/tst/scripts/disk/disk.py @@ -68,6 +68,7 @@ def run(**kwargs): g, int(10 * gam), b ), "problem/polytropic_index={:.2f}".format(gam), + "gas/de_switch=" + str(0.0 if g != "sph" else 1e-2), ] + geom_args, )