Refactor stopping times

src/artemis.cpp

@@ -23,7 +23,9 @@
 #include "nbody/nbody.hpp"
 #include "rotating_frame/rotating_frame.hpp"
 #include "utils/artemis_utils.hpp"
+#include "utils/eos/eos.hpp"
 #include "utils/history.hpp"
+#include "utils/opacity/opacity.hpp"
 #include "utils/units.hpp"
 
 // Jaybenne includes
@@ -104,7 +106,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   if (do_dust) packages.Add(Dust::Initialize(pin.get(), units));
   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_drag) packages.Add(Drag::Initialize(pin.get(), constants, units));
   if (do_radiation) {
     auto eos_h = packages.Get("gas")->Param<EOS>("eos_h");
     auto opacity_h = packages.Get("gas")->Param<MeanOpacity>("opacity_h");

src/artemis.hpp

@@ -21,9 +21,6 @@
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 
-// Singularity-eos includes
-#include <singularity-eos/eos/eos.hpp>
-
 using namespace parthenon;
 using namespace parthenon::driver::prelude;
 using namespace parthenon::package::prelude;

src/drag/drag.cpp

@@ -16,13 +16,16 @@
 #include "artemis.hpp"
 #include "geometry/geometry.hpp"
 #include "utils/eos/eos.hpp"
+#include "utils/units.hpp"
 
 using ArtemisUtils::EOS;
 namespace Drag {
 //----------------------------------------------------------------------------------------
 //! \fn  StateDescriptor Drag::Initialize
 //! \brief Adds intialization function for damping package
-std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
+                                            ArtemisUtils::Constants &constants,
+                                            ArtemisUtils::Units &units) {
   auto drag = std::make_shared<StateDescriptor>("drag");
   Params &params = drag->AllParams();
 
@@ -82,7 +85,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
     // Enforce 1 gas species
 
-    params.Add("stopping_time_params", StoppingTimeParams("dust/stopping_time", pin));
+    params.Add("stopping_time_params",
+               StoppingTimeParams("dust/stopping_time", pin, constants, units));
   }
 
   return drag;
@@ -139,6 +143,10 @@ TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt) {
           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::powerlaw) {
+          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_plaw,
+                                      DragModel::powerlaw, 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_plaw,
                                       DragModel::stokes, GEOM>(
@@ -149,6 +157,10 @@ TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt) {
           return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_alpha,
                                       DragModel::constant, GEOM>(
               md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par);
+        } else if (stop_par.model == DragModel::powerlaw) {
+          return SimpleDragSourceImpl<Diffusion::DiffType::viscosity_alpha,
+                                      DragModel::powerlaw, 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_alpha,
                                       DragModel::stokes, GEOM>(
@@ -162,6 +174,9 @@ TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt) {
       if (stop_par.model == DragModel::constant) {
         return SimpleDragSourceImpl<Diffusion::DiffType::null, DragModel::constant, GEOM>(
             md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par);
+      } else if (stop_par.model == DragModel::powerlaw) {
+        return SimpleDragSourceImpl<Diffusion::DiffType::null, DragModel::powerlaw, 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::null, DragModel::stokes, GEOM>(
             md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par);

src/drag/drag.hpp

@@ -22,6 +22,7 @@
 #include "utils/artemis_utils.hpp"
 #include "utils/diffusion/diffusion_coeff.hpp"
 #include "utils/eos/eos.hpp"
+#include "utils/units.hpp"
 
 using namespace parthenon::package::prelude;
 using ArtemisUtils::EOS;
@@ -54,7 +55,7 @@ namespace Drag {
 */
 
 enum class Coupling { simple_dust, self, null };
-enum class DragModel { constant, stokes, null };
+enum class DragModel { constant, stokes, dp15, powerlaw, null };
 
 inline Coupling ChooseDrag(const std::string choice) {
   if (choice == "self") {
@@ -111,10 +112,11 @@ struct SelfDragParams {
 
 struct StoppingTimeParams {
 
-  Real scale;
+  Real scale, dh, mass_scale, p1, p2, p3, rho_plaw, x_plaw;
   DragModel model;
   ParArray1D<Real> tau;
-  StoppingTimeParams(std::string block_name, ParameterInput *pin) {
+  StoppingTimeParams(std::string block_name, ParameterInput *pin,
+                     ArtemisUtils::Constants &constants, ArtemisUtils::Units &units) {
     const std::string choice = pin->GetString(block_name, "type");
     const int nd = pin->GetOrAddInteger("dust", "nspecies", 1);
     tau = ParArray1D<Real>("tau", nd);
@@ -137,13 +139,147 @@ struct StoppingTimeParams {
         h_tau(n) = scale;
       }
       tau.DeepCopy(h_tau);
+    } else if (choice == "dp15") {
+      model = DragModel::dp15;
+      scale = pin->GetOrAddReal(block_name, "scale", 1.0);
+      auto h_tau = tau.GetHostMirror();
+      for (int n = 0; n < nd; n++) {
+        h_tau(n) = scale;
+      }
+      tau.DeepCopy(h_tau);
+      dh = units.GetLengthPhysicalToCode() *
+           pin->GetOrAddReal(block_name, "gas_diameter", 2.71e-8 /* cm */);
+      p1 = pin->GetOrAddReal(block_name, "p1", 3.07);
+      p2 = pin->GetOrAddReal(block_name, "p2", 0.6688);
+      p3 = pin->GetOrAddReal(block_name, "p3", 0.681);
+      // save this so that we don't have to pass units and constants struct to the Get()
+      // routines
+      mass_scale = units.GetMassPhysicalToCode() * constants.GetAMUCode();
+
+    } else if (choice == "powerlaw") {
+      model = DragModel::constant;
+      scale = pin->GetOrAddReal(block_name, "scale", 1.0);
+      rho_plaw = pin->GetOrAddReal(block_name, "density_plaw", 0.0);
+      x_plaw = pin->GetOrAddReal(block_name, "r_plaw", 0.0);
+      std::vector<Real> taus = pin->GetVector<Real>(block_name, "tau");
+      auto h_tau = tau.GetHostMirror();
+      for (int n = 0; n < nd; n++) {
+        h_tau(n) = scale * taus[n];
+      }
+      tau.DeepCopy(h_tau);
     } else {
       PARTHENON_FAIL("bad type for stopping time model");
     }
   }
 };
 
-std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+template <DragModel DTYP>
+class DragCoeff {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    PARTHENON_FAIL("No default implementation for drag coefficient");
+  }
+};
+
+// null
+template <>
+class DragCoeff<DragModel::null> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    return Big<Real>();
+  }
+};
+
+template <>
+class DragCoeff<DragModel::constant> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    return dp.tau(id);
+  }
+};
+
+template <>
+class DragCoeff<DragModel::powerlaw> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    return dp.tau(id) * std::pow(dg, dp.rho_plaw) * std::pow(std::abs(xc[0]), dp.x_plaw);
+  }
+};
+
+template <>
+class DragCoeff<DragModel::stokes> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    const Real cv = eos.SpecificHeatFromDensityTemperature(dg, Tg);
+    const Real gm1 = eos.GruneisenParamFromDensityTemperature(dg, Tg);
+    //  kb T/ mu = cv * gm1 * T
+    const Real vth = std::sqrt(8.0 / M_PI * gm1 * cv * Tg);
+    return dp.tau(id) * grain_density / dg * size / vth;
+  }
+};
+
+template <>
+class DragCoeff<DragModel::dp15> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  const std::array<Real, 3> xc, const Real dg,
+                                  const Real Tg, const Real u, const Real grain_density,
+                                  const Real size, const EOS &eos) const {
+    const Real cv = eos.SpecificHeatFromDensityTemperature(dg, Tg);
+    const Real gm1 = eos.GruneisenParamFromDensityTemperature(dg, Tg);
+    const Real mu = dp.mass_scale * eos.MeanAtomicMass();
+    const Real gamma = gm1 + 1.;
+    const Real sgam = std::sqrt(gamma);
+    //  kb T/ mu = cv * gm1 * T
+
+    const Real vth = std::sqrt(8.0 / M_PI * gm1 * cv * Tg);
+    // Mach and Reynolds numbers
+    // Mu is non-zero, M can be zero
+    const Real Mu = std::sqrt(8. / (M_PI * gamma)) / vth;
+    const Real M = u * Mu;
+    const Real K = 5. / (32. * std::sqrt(M_PI)) * mu / SQR(dp.dh) / (dg * size * sgam);
+    // Ru is non-zero, R can be zero
+    const Real Ru = Mu / K;
+    const Real R = Ru * u;
+
+    // CD = 2 + (CS - 2) * exp(-p1 sqrt(g) K * G(R)) + CE * exp(-1/(2*K))
+
+    // The two coefficients (multiplied by u)
+    const Real CEu = 1. / (sgam * Mu) * (4.6 / (1. + M) + 1.7 /* srt(Td/Tg) */);
+    const Real CSu =
+        24. / Ru * (1. + 0.15 * std::pow(R, dp.p3)) + 0.407 * R * u / (R + 8710.);
+
+    // weight functions
+    const Real x = std::pow(R / 312., dp.p2);
+    const Real G = std::exp(2.5 * x / (1. + x));
+    const Real ws = std::exp(-dp.p1 * sgam * K * G);
+    const Real we = std::exp(-1. / (2. * K));
+
+    // The final stopping time
+    const Real CDu = 2 * u * (1. - ws) + CSu * ws + CEu * we;
+    const Real alpha = 3. / 8 * CDu * dg / (grain_density * size);
+    return dp.tau(id) / (alpha + Fuzz<Real>());
+  }
+};
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
+                                            ArtemisUtils::Constants &constants,
+                                            ArtemisUtils::Units &units);
 
 template <Coordinates GEOM>
 TaskStatus DragSource(MeshData<Real> *md, const Real time, const Real dt);
@@ -398,13 +534,8 @@ TaskStatus SimpleDragSourceImpl(MeshData<Real> *md, const Real time, const Real
         std::array<Real, 3> fvd{0.0, 0.0, 0.0};
         const auto nspecies = vmesh.GetSize(b, dust::cons::density());
 
-        [[maybe_unused]] auto &grain_density_ = grain_density;
-        [[maybe_unused]] Real vth = Null<Real>();
-
-        if constexpr (DRAG == DragModel::stokes) {
-          const Real gm1 = eos_d.GruneisenParamFromDensityInternalEnergy(dg, sieg);
-          vth = std::sqrt(8.0 / M_PI * gm1 * sieg);
-        }
+        DragCoeff<DRAG> drag_coeff;
+        Real Tg = eos_d.TemperatureFromDensityInternalEnergy(dg, sieg);
 
         // First pass to collect \sum rho' and \sum rho' v and compute new vg
         std::array<Real, 3> vdt{0.0, 0.0, 0.0};
@@ -415,11 +546,15 @@ TaskStatus SimpleDragSourceImpl(MeshData<Real> *md, const Real time, const Real
               vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * dens),
               vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * dens),
               vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * dens)};
-          Real tc = tp.tau(id);
-          [[maybe_unused]] auto &sizes_ = sizes;
-          if constexpr (DRAG == DragModel::stokes) {
-            tc = tp.scale * grain_density_ / dg * sizes_(id) / vth;
-          }
+
+          // relative speed
+          const Real u =
+              std::sqrt(SQR(vd[0] - vg[0]) + SQR(vd[1] - vg[1]) + SQR(vd[2] - vg[2]));
+
+          // Get the stopping time
+          Real tc =
+              drag_coeff.Get(tp, id, xv, dg, Tg, u, grain_density, sizes(id), eos_d);
+
           const Real alpha = dt * ((tc <= 0.0) ? Big<Real>() : 1.0 / tc);
           for (int d = 0; d < 3; d++) {
             const Real rhop = dens * alpha / (1.0 + alpha + bd[d]);
@@ -447,11 +582,14 @@ TaskStatus SimpleDragSourceImpl(MeshData<Real> *md, const Real time, const Real
               vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * dens),
               vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * dens),
               vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * dens)};
-          Real tc = tp.tau(id);
-          [[maybe_unused]] auto &sizes_ = sizes;
-          if constexpr (DRAG == DragModel::stokes) {
-            tc = tp.scale * grain_density_ / dg * sizes_(id) / vth;
-          }
+          // relative speed
+          const Real u =
+              std::sqrt(SQR(vd[0] - vg[0]) + SQR(vd[1] - vg[1]) + SQR(vd[2] - vg[2]));
+
+          // Get the stopping time
+          Real tc =
+              drag_coeff.Get(tp, id, xv, dg, Tg, u, grain_density, sizes(id), eos_d);
+
           const Real alpha = dt * ((tc <= 0.0) ? Big<Real>() : 1.0 / tc);
           for (int d = 0; d < 3; d++) {
             Real delta_d = 0.;

src/dust/params.yaml

@@ -120,5 +120,8 @@ dust:
       stokes:
         _type: opt
         _description: "Stopping times between the gas and dust fluid are in the Stokes regime"
+      dp15:
+        _type: opt
+        _description: "Stopping times between the gas and dust fluid use the model described in Appendix A of D'Angelo & Podolok (2015)."
 
 

src/gas/gas.cpp

@@ -110,15 +110,22 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
       cv = pin->GetReal("gas", "cv");
       PARTHENON_REQUIRE(cv > 0, "Only positive cv allowed!");
       mu = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * cv);
+      if (pin->DoesParameterExist("gas", "mu")) {
+        PARTHENON_FAIL("Cannot set both mu and cv in the input file. One is calculated "
+                       "from the other!");
+      }
     } else {
       mu = pin->GetOrAddReal("gas", "mu", 1.);
       PARTHENON_REQUIRE(mu > 0, "Only positive mean molecular weight allowed!");
       cv = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * mu);
     }
-    params.Add("mu", mu);
+    // Store these so they are available in outputs
     params.Add("cv", cv);
+    params.Add("mu", mu);
+    auto zbar = pin->GetOrAddReal("gas", "zbar", 1.0);
     EOS eos_host = singularity::UnitSystem<singularity::IdealGas>(
-        singularity::IdealGas(gamma - 1., cv * units.GetSpecificHeatCodeToPhysical()),
+        singularity::IdealGas(gamma - 1., cv * units.GetSpecificHeatCodeToPhysical(),
+                              {mu, zbar}),
         singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimeCodeToPhysical(),
         units.GetMassCodeToPhysical(), units.GetLengthCodeToPhysical(),
         units.GetTemperatureCodeToPhysical());

src/utils/eos/eos.hpp

@@ -13,6 +13,8 @@
 #ifndef UTILS_EOS_HPP_
 #define UTILS_EOS_HPP_
 
+#include <singularity-eos/eos/eos.hpp>
+
 #include "artemis.hpp"
 
 namespace ArtemisUtils {