diff --git a/src/artemis.cpp b/src/artemis.cpp index 8c674300..411ded8f 100644 --- a/src/artemis.cpp +++ b/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 &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_h"); auto opacity_h = packages.Get("gas")->Param("opacity_h"); diff --git a/src/artemis.hpp b/src/artemis.hpp index 8e845aa3..e4a68c7a 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -21,9 +21,6 @@ #include #include -// Singularity-eos includes -#include - using namespace parthenon; using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; diff --git a/src/drag/drag.cpp b/src/drag/drag.cpp index 0d33968a..d51c435a 100644 --- a/src/drag/drag.cpp +++ b/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 Initialize(ParameterInput *pin) { +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Constants &constants, + ArtemisUtils::Units &units) { auto drag = std::make_shared("drag"); Params ¶ms = drag->AllParams(); @@ -82,7 +85,8 @@ std::shared_ptr 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 *md, const Real time, const Real dt) { return SimpleDragSourceImpl( md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + } else if (stop_par.model == DragModel::powerlaw) { + return SimpleDragSourceImpl( + md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( @@ -149,6 +157,10 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { return SimpleDragSourceImpl( md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + } else if (stop_par.model == DragModel::powerlaw) { + return SimpleDragSourceImpl( + md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( @@ -162,6 +174,9 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { if (stop_par.model == DragModel::constant) { return SimpleDragSourceImpl( md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + } else if (stop_par.model == DragModel::powerlaw) { + return SimpleDragSourceImpl( + md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index 5b6fc724..213c0534 100644 --- a/src/drag/drag.hpp +++ b/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 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("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 taus = pin->GetVector(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 Initialize(ParameterInput *pin); +template +class DragCoeff { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array 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 { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array xc, const Real dg, + const Real Tg, const Real u, const Real grain_density, + const Real size, const EOS &eos) const { + return Big(); + } +}; + +template <> +class DragCoeff { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array 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 { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array 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 { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array 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 { + public: + KOKKOS_INLINE_FUNCTION Real Get(const StoppingTimeParams &dp, const int id, + const std::array 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()); + } +}; + +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Constants &constants, + ArtemisUtils::Units &units); template TaskStatus DragSource(MeshData *md, const Real time, const Real dt); @@ -398,13 +534,8 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real std::array 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(); - - if constexpr (DRAG == DragModel::stokes) { - const Real gm1 = eos_d.GruneisenParamFromDensityInternalEnergy(dg, sieg); - vth = std::sqrt(8.0 / M_PI * gm1 * sieg); - } + DragCoeff 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 vdt{0.0, 0.0, 0.0}; @@ -415,11 +546,15 @@ TaskStatus SimpleDragSourceImpl(MeshData *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() : 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 *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() : 1.0 / tc); for (int d = 0; d < 3; d++) { Real delta_d = 0.; diff --git a/src/dust/params.yaml b/src/dust/params.yaml index e0dedbca..1b5afebe 100644 --- a/src/dust/params.yaml +++ b/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)." diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 290f169c..3264317a 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -110,15 +110,22 @@ std::shared_ptr 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(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()); diff --git a/src/utils/eos/eos.hpp b/src/utils/eos/eos.hpp index d1703fab..8417cb81 100644 --- a/src/utils/eos/eos.hpp +++ b/src/utils/eos/eos.hpp @@ -13,6 +13,8 @@ #ifndef UTILS_EOS_HPP_ #define UTILS_EOS_HPP_ +#include + #include "artemis.hpp" namespace ArtemisUtils {