Merge pull request #43 from lanl/dempsey/dust_update

Convert dust params from cgs
lanl · Jan 14, 2025 · b9e1093 · b9e1093
2 parents c98a205 + 72984f6
commit b9e1093
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 14 deletions.
diff --git a/src/artemis.cpp b/src/artemis.cpp
@@ -101,7 +101,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   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_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()));

diff --git a/src/dust/dust.cpp b/src/dust/dust.cpp
@@ -25,14 +25,16 @@
 #include "utils/artemis_utils.hpp"
 #include "utils/fluxes/fluid_fluxes.hpp"
 #include "utils/history.hpp"
+#include "utils/units.hpp"
 
 using ArtemisUtils::VI;
 
 namespace Dust {
 //----------------------------------------------------------------------------------------
 //! \fn  StateDescriptor Dust::Initialize
 //! \brief Adds intialization function for dust hydrodynamics package
-std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
+                                            ArtemisUtils::Units &units) {
   auto dust = std::make_shared<StateDescriptor>("dust");
   Params &params = dust->AllParams();
 
@@ -99,10 +101,13 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   // Dust sizes
   const auto size_dist = pin->GetOrAddString("dust", "size_input", "direct");
+  const Real length_conv = units.GetLengthPhysicalToCode();
+  const Real rho_conv = units.GetMassDensityPhysicalToCode();
+
   if (size_dist == "linspace") {
     // uniform
-    auto min_size = pin->GetReal("dust", "min_size");
-    auto max_size = pin->GetReal("dust", "max_size");
+    auto min_size = length_conv * pin->GetReal("dust", "min_size");
+    auto max_size = length_conv * pin->GetReal("dust", "max_size");
 
     ParArray1D<Real> sizes("sizes", nspecies);
     auto h_sizes = sizes.GetHostMirror();
@@ -117,8 +122,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   } else if (size_dist == "logspace") {
     // uniform in log-space
-    const auto lmin = std::log10(pin->GetReal("dust", "min_size"));
-    const auto lmax = std::log10(pin->GetReal("dust", "max_size"));
+    const auto lmin = std::log10(length_conv * pin->GetReal("dust", "min_size"));
+    const auto lmax = std::log10(length_conv * pin->GetReal("dust", "max_size"));
 
     ParArray1D<Real> sizes("sizes", nspecies);
     auto h_sizes = sizes.GetHostMirror();
@@ -138,7 +143,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     ParArray1D<Real> sizes("sizes", nspecies);
     auto h_sizes = sizes.GetHostMirror();
     for (int n = 0; n < nspecies; n++) {
-      h_sizes(n) = sizes_v[n];
+      h_sizes(n) = length_conv * sizes_v[n];
     }
     sizes.DeepCopy(h_sizes);
     params.Add("sizes", sizes);
@@ -153,11 +158,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     auto h_sizes = sizes.GetHostMirror();
     if (data.size() == 1) {
       for (int n = 0; n < nspecies; n++) {
-        h_sizes(n) = data[0][n];
+        h_sizes(n) = length_conv * data[0][n];
       }
     } else {
       for (int n = 0; n < nspecies; n++) {
-        h_sizes(n) = data[n][0];
+        h_sizes(n) = length_conv * data[n][0];
       }
     }
     sizes.DeepCopy(h_sizes);
@@ -168,7 +173,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   }
 
   // Dust density
-  params.Add("grain_density", pin->GetOrAddReal("dust", "grain_density", 1.0));
+  params.Add("grain_density", rho_conv * pin->GetOrAddReal("dust", "grain_density", 1.0));
 
   // Scratch for dust flux
   const int scr_level = pin->GetOrAddInteger("dust", "scr_level", 0);

diff --git a/src/dust/dust.hpp b/src/dust/dust.hpp
@@ -13,17 +13,16 @@
 #ifndef DUST_DUST_HPP_
 #define DUST_DUST_HPP_
 
-// Parthenon includes
-#include <parthenon/package.hpp>
-
 // Artemis includes
 #include "artemis.hpp"
+#include "utils/units.hpp"
 
 using namespace parthenon::package::prelude;
 
 namespace Dust {
 
-std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
+                                            ArtemisUtils::Units &units);
 
 template <Coordinates GEOM>
 Real EstimateTimestepMesh(MeshData<Real> *md);