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,11 +106,12 @@ 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<Opacity>("opacity_h");
-    auto scattering_h = packages.Get("gas")->Param<Scattering>("scattering_h");
+    auto eos_h = packages.Get("gas")->Param<ArtemisUtils::EOS>("eos_h");
+    auto opacity_h = packages.Get("gas")->Param<ArtemisUtils::Opacity>("opacity_h");
+    auto scattering_h =
+        packages.Get("gas")->Param<ArtemisUtils::Scattering>("scattering_h");
     packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h));
     PARTHENON_REQUIRE(coords == Coordinates::cartesian,
                       "Jaybenne currently supports only Cartesian coordinates!");

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;

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, 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;
   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,125 @@ 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 {
       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 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 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 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::stokes> {
+ public:
+  KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id,
+                                  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 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 +512,8 @@ TaskStatus SimpleDragSourceImpl(MeshData<Real> *md, const Real time, const Real
         Real fvd[3] = {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
         const Real vdt[3] = {0.0};
@@ -415,11 +524,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, 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 +559,13 @@ 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, 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

@@ -103,17 +103,28 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
   if (eos_name == "ideal") {
     const Real gamma = pin->GetOrAddReal("gas", "gamma", 1.66666666667);
     auto cv = Null<Real>();
+    auto mu = Null<Real>();
+
     if (pin->DoesParameterExist("gas", "cv")) {
-      PARTHENON_REQUIRE(!pin->DoesParameterExist("gas", "mmw"),
-                        "Cannot specify both cv and mmw");
+      PARTHENON_REQUIRE(!pin->DoesParameterExist("gas", "mu"),
+                        "Cannot specify both cv and mu");
       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 {
-      const Real mu = pin->GetOrAddReal("gas", "mu", 1.);
+      mu = pin->GetOrAddReal("gas", "mu", 1.);
       PARTHENON_REQUIRE(mu > 0, "Only positive mean molecular weight allowed!");
       cv = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * mu);
     }
-    EOS eos_host = singularity::IdealGas(gamma - 1., cv);
+    // 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::IdealGas(gamma - 1., cv, {mu, zbar});
     EOS eos_device = eos_host.GetOnDevice();
     params.Add("eos_h", eos_host);
     params.Add("eos_d", eos_device);

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 {