Merge branch 'develop' into STS

.github/workflows/ci.yml

@@ -40,11 +40,11 @@ jobs:
           python-version: '3.x'  # Specify the Python version you need
       - name: Install dependencies
         run: |
+          pip install clang-format==12.0.1.2
           pip install black
       - name: Run format check
         run: |
-          source env/bash
-          VERBOSE=1 ./style/format.sh
+          VERBOSE=1 CFM=clang-format ./style/format.sh
           git diff --exit-code --ignore-submodules
 
   cpu:

doc/src/physics.rst

@@ -673,10 +673,7 @@ Specifically, |code|, adds the body force,
 .. math::
    \rho \frac{D \mathbf{v}}{D t} = \rho \mathbf{g}
 
-To activate external gravity, there must be a ``<gravity>`` node and ``gravity = true`` under the ``<physics>`` node.
-The ``<gravity>`` requires setting the parameter ``gm``, typically this is set to ``gm = 1.0``.
-In addition to ``gm``, the ``tstart`` and ``tstop`` parameters control when gravity is active.
-
+To activate external gravity, there must be a subnode of ``<gravity>`` (e.g., ``<gravity/point>``) and ``gravity = true`` under the ``<physics>`` node.
 The specific model for the gravitational acceleration, :math:`\mathbf{g}`, is controlled by adding the appropriate subnode.
 Available options are:
 
@@ -687,8 +684,6 @@ Available options are:
 
   ::
 
-   <gravity>
-   gm = 1.0
    <gravity/constant>
    gx1 = 0.0
    gx2 = 0.0
@@ -703,9 +698,8 @@ Available options are:
 
   ::
 
-   <gravity>
-   gm = 1.0
    <gravity/point>
+   mass = 1.0
    x = 0.0
    y = 0.0
    z = 0.0
@@ -721,9 +715,8 @@ Available options are:
 
   ::
 
-   <gravity>
-   gm = 1.0            # Binary total GM
    <gravity/binary>
+   mass = 1.0          # Total binary mass
    x = 0.0             # Binary x center of mass
    y = 0.0             # Binary y center of mass
    z = 0.0             # Binary z center of mass
@@ -745,14 +738,11 @@ Available options are:
 * ``<gravity/nbody>``
 
   This indicates that the gravitational force will be calculated by the N-body system defined in the ``<nbody>`` input block.
-  The only parameter required is the ``gm`` parameter.
-  Note that the total mass of the system defined in the ``<nbody>`` block will be rescaled to ``gm``.
-  An example input block would thus read:
+  Unlike the other gravity nodes, ``<gravity/nbody>`` has no parameters. Instead, all parameters live under the ``<nbody>`` node.
+  To activate nbody gravity, simply add:
 
   ::
 
-   <gravity>
-   gm = 1.0
    <gravity/nbody>
 
 See `N-Body Dynamics`_ for a description of how to set up the N-body system.

inputs/diffusion/alpha_disk.in

@@ -83,9 +83,8 @@ beta0 = 0.0
 tcyl = .0025
 cyl_plaw = -1.0
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <problem>
 r0 = 1.0

inputs/diffusion/conduction.in

@@ -71,7 +71,7 @@ dfloor = 1.0e-10
 siefloor = 1.0e-15
 
 <gas/conductivity>
-type = constant 
+type = conductivity 
 cond = 1.0e-1
 
 <gravity/uniform>

inputs/diffusion/gaussian_bump.in

@@ -80,7 +80,7 @@ type = constant
 nu = 5.0e-3
 
 <gas/conductivity>
-type = constant
+type = conductivity 
 cond = 5.0e-3
 
 <problem>

inputs/disk/binary_cyl.in

@@ -86,9 +86,8 @@ omega = 1.0
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/binary>
+mass = 1.0
 q = 1.0e-5
 a = 1.0
 e = 0.0

inputs/disk/binary_nbody_cyl.in

@@ -92,9 +92,8 @@ omega = 1.0
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/nbody>
+mtot = 1.0
 
 <nbody>
 dt_reb = 0.01

inputs/disk/cb_disk.in

@@ -81,9 +81,8 @@ siefloor = 1.0e-6
 type = constant
 nu = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/nbody>
+mtot = 1.0
 
 <nbody>
 integrator = ias15

inputs/disk/disk_alpha.in

@@ -74,9 +74,8 @@ siefloor = 1.0e-10
 type = alpha
 alpha = 1e-2
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <cooling>
 type = beta

inputs/disk/disk_axi.in

@@ -80,9 +80,8 @@ omega = 1.0
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <problem>
 r0 = 1.0

inputs/disk/disk_cart.in

@@ -110,9 +110,8 @@ riemann = hlle
 dfloor = 1.0e-10
 pfloor = 1.0e-15
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <problem>
 r0 = 1.0

inputs/disk/disk_collision.in

@@ -79,9 +79,8 @@ riemann = hllc
 dfloor = 1.0e-10
 siefloor = 1.0e-10
 
-<gravity>
-mass_tot = 1.0
 <gravity/nbody>
+mtot = 1.0
 
 <gas/damping>
 inner_x1 = 0.45

inputs/disk/disk_cyl.in

@@ -81,9 +81,8 @@ omega = 1.0
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <problem>
 r0 = 1.0

inputs/disk/disk_nbody_cyl.in

@@ -78,8 +78,6 @@ siefloor = 1.0e-10
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/nbody>
 
 <nbody>

inputs/disk/disk_sph.in

@@ -80,9 +80,8 @@ omega = 1.0
 type = alpha
 alpha = 1e-3
 
-<gravity>
-mass_tot = 1.0
 <gravity/point>
+mass = 1.0
 
 <problem>
 r0 = 1.0

inputs/ssheet/ssheet.in

@@ -76,9 +76,8 @@ refine_type = magnitude
 refine_thr = 3.0
 deref_thr = 0.8
 
-<gravity>
-mass_tot = 1.0e-5
 <gravity/point>
+mass = 1.0e-5
 soft = 0.03
 x = 0.0
 y = 0.0

src/artemis.cpp

@@ -108,14 +108,15 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   artemis->AddParam("coord_sys", sys);
 
   // Call package initializers here
-  if (do_gas) packages.Add(Gas::Initialize(pin.get(), units, constants));
+  if (do_nbody) packages.Add(NBody::Initialize(pin.get(), constants));
+  if (do_gravity) packages.Add(Gravity::Initialize(pin.get(), constants, packages));
+  if (do_gas) packages.Add(Gas::Initialize(pin.get(), units, constants, packages));
   if (do_dust) packages.Add(Dust::Initialize(pin.get()));
-  if (do_gravity) packages.Add(Gravity::Initialize(pin.get(), constants));
   if (do_rotating_frame) packages.Add(RotatingFrame::Initialize(pin.get()));
   if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get()));
   if (do_drag) packages.Add(Drag::Initialize(pin.get()));
-  if (do_nbody) packages.Add(NBody::Initialize(pin.get(), constants));
   if (do_sts) packages.Add(STS::Initialize(pin.get()));
+
   if (do_radiation) {
     auto eos_h = packages.Get("gas")->Param<EOS>("eos_h");
     auto opacity_h = packages.Get("gas")->Param<Opacity>("opacity_h");

src/artemis_params.yaml

@@ -23,13 +23,11 @@ artemis:
     _type: bool
     _default: "false"
     _description: "Add user-defined AMR criterion."
-  units:
-    _description: "Unit conversions between code and cgs."
   physical_units:
     _type: "string"
     _default: "scalefree"
     _description: "What unit system to use for physical units"
-    code:
+    scalefree:
       _type: opt
       _description: "Native code units; no conversions"
     cgs:
@@ -47,12 +45,15 @@ artemis:
       _description: "AU, Year/(2 pi), solar mass units for protoplanetary disks"
   length:
     _type: "Real"
+    _default: "1.0"
     _description: "Physical units value of length equal to 1 code unit"
   time:
     _type: "Real"
+    _default: "1.0"
     _description: "Physical units value of time equal to 1 code unit"
   mass:
     _type: "Real"
+    _default: "1.0"
     _description: "Physical units value of mass equal to 1 code unit"
 
 

src/drag/drag.cpp

@@ -110,8 +110,8 @@ TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt) {
       auto &gas_pkg = pm->packages.Get("gas");
       const auto &dp = gas_pkg->template Param<Diffusion::DiffCoeffParams>("visc_params");
       const auto &eos_d = gas_pkg->template Param<EOS>("eos_d");
-      if (dp.type == Diffusion::DiffType::viscosity_const) {
-        return SelfDragSourceImpl<Diffusion::DiffType::viscosity_const, GEOM>(
+      if (dp.type == Diffusion::DiffType::viscosity_plaw) {
+        return SelfDragSourceImpl<Diffusion::DiffType::viscosity_plaw, GEOM>(
             md, time, dt, dp, eos_d, gas_self_par, dust_self_par);
       } else if (dp.type == Diffusion::DiffType::viscosity_alpha) {
         return SelfDragSourceImpl<Diffusion::DiffType::viscosity_alpha, GEOM>(
@@ -134,13 +134,13 @@ TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt) {
         drag_pkg->template Param<StoppingTimeParams>("stopping_time_params");
     if (gas_self_par.damp_to_visc) {
       auto &dp = gas_pkg->template Param<Diffusion::DiffCoeffParams>("visc_params");
-      if (dp.type == Diffusion::DiffType::viscosity_const) {
+      if (dp.type == Diffusion::DiffType::viscosity_plaw) {
         if (stop_par.model == DragModel::constant) {
-          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_const,
+          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_plaw,
                                       DragModel::constant, GEOM>(
               md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par);
         } else if (stop_par.model == DragModel::stokes) {
-          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_const,
+          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_plaw,
                                       DragModel::stokes, GEOM>(
               md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par);
         }

src/gas/gas.cpp

@@ -40,7 +40,8 @@ namespace Gas {
 //! \brief Adds intialization function for gas hydrodynamics package
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
                                             ArtemisUtils::Units &units,
-                                            ArtemisUtils::Constants &constants) {
+                                            ArtemisUtils::Constants &constants,
+                                            Packages_t &packages) {
   using namespace singularity::photons;
 
   auto gas = std::make_shared<StateDescriptor>("gas");
@@ -195,11 +196,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
   params.Add("do_diffusion", do_diffusion);
 
   if (do_viscosity) {
-    Diffusion::DiffCoeffParams dp("gas/viscosity", "viscosity", pin, constants);
+    Diffusion::DiffCoeffParams dp("gas/viscosity", "viscosity", pin, constants, packages);
     params.Add("visc_params", dp);
   }
   if (do_conduction) {
-    Diffusion::DiffCoeffParams dp("gas/conductivity", "conductivity", pin, constants);
+    Diffusion::DiffCoeffParams dp("gas/conductivity", "conductivity", pin, constants,
+                                  packages);
     params.Add("cond_params", dp);
   }
 
@@ -444,9 +446,9 @@ Real EstimateTimestepMesh(MeshData<Real> *md) {
   const auto do_viscosity = params.template Get<bool>("do_viscosity");
   if (do_viscosity) {
     auto dp = params.template Get<Diffusion::DiffCoeffParams>("visc_params");
-    if (dp.type == Diffusion::DiffType::viscosity_const) {
+    if (dp.type == Diffusion::DiffType::viscosity_plaw) {
       visc_dt = Diffusion::EstimateTimestep<GEOM, Fluid::gas,
-                                            Diffusion::DiffType::viscosity_const>(
+                                            Diffusion::DiffType::viscosity_plaw>(
           md, dp, gas_pkg, eos_d, vmesh);
     } else if (dp.type == Diffusion::DiffType::viscosity_alpha) {
       visc_dt = Diffusion::EstimateTimestep<GEOM, Fluid::gas,
@@ -459,13 +461,13 @@ Real EstimateTimestepMesh(MeshData<Real> *md) {
   const auto do_conduction = params.template Get<bool>("do_conduction");
   if (do_conduction) {
     auto dp = params.template Get<Diffusion::DiffCoeffParams>("cond_params");
-    if (dp.type == Diffusion::DiffType::conductivity_const) {
+    if (dp.type == Diffusion::DiffType::conductivity_plaw) {
       cond_dt = Diffusion::EstimateTimestep<GEOM, Fluid::gas,
-                                            Diffusion::DiffType::conductivity_const>(
+                                            Diffusion::DiffType::conductivity_plaw>(
           md, dp, gas_pkg, eos_d, vmesh);
-    } else if (dp.type == Diffusion::DiffType::thermaldiff_const) {
+    } else if (dp.type == Diffusion::DiffType::thermaldiff_plaw) {
       cond_dt = Diffusion::EstimateTimestep<GEOM, Fluid::gas,
-                                            Diffusion::DiffType::thermaldiff_const>(
+                                            Diffusion::DiffType::thermaldiff_plaw>(
           md, dp, gas_pkg, eos_d, vmesh);
     }
   }
@@ -566,10 +568,10 @@ TaskStatus ViscousFlux(MeshData<Real> *md) {
 
   if (dp.type == Diffusion::DiffType::null) {
     return TaskStatus::complete;
-  } else if (dp.type == Diffusion::DiffType::viscosity_const) {
+  } else if (dp.type == Diffusion::DiffType::viscosity_plaw) {
     return Diffusion::MomentumFluxImpl<GEOM, Fluid::gas,
-                                       Diffusion::DiffType::viscosity_const>(md, dp, pkg,
-                                                                             vprim, vf);
+                                       Diffusion::DiffType::viscosity_plaw>(md, dp, pkg,
+                                                                            vprim, vf);
   } else if (dp.type == Diffusion::DiffType::viscosity_alpha) {
     return Diffusion::MomentumFluxImpl<GEOM, Fluid::gas,
                                        Diffusion::DiffType::viscosity_alpha>(md, dp, pkg,
@@ -604,14 +606,14 @@ TaskStatus ThermalFlux(MeshData<Real> *md) {
 
   if (dp.type == Diffusion::DiffType::null) {
     return TaskStatus::complete;
-  } else if (dp.type == Diffusion::DiffType::conductivity_const) {
-    return Diffusion::ThermalFluxImpl<GEOM, Fluid::gas,
-                                      Diffusion::DiffType::conductivity_const>(
-        md, dp, pkg, vprim, vf);
-  } else if (dp.type == Diffusion::DiffType::thermaldiff_const) {
+  } else if (dp.type == Diffusion::DiffType::conductivity_plaw) {
     return Diffusion::ThermalFluxImpl<GEOM, Fluid::gas,
-                                      Diffusion::DiffType::thermaldiff_const>(md, dp, pkg,
+                                      Diffusion::DiffType::conductivity_plaw>(md, dp, pkg,
                                                                               vprim, vf);
+  } else if (dp.type == Diffusion::DiffType::thermaldiff_plaw) {
+    return Diffusion::ThermalFluxImpl<GEOM, Fluid::gas,
+                                      Diffusion::DiffType::thermaldiff_plaw>(md, dp, pkg,
+                                                                             vprim, vf);
   } else {
     PARTHENON_FAIL("Invalid conductivity type");
   }

src/gas/gas.hpp

@@ -20,7 +20,8 @@ namespace Gas {
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
                                             ArtemisUtils::Units &units,
-                                            ArtemisUtils::Constants &constants);
+                                            ArtemisUtils::Constants &constants,
+                                            Packages_t &packages);
 
 template <Coordinates GEOM>
 Real EstimateTimestepMesh(MeshData<Real> *md);