Set dust grain size from input file

src/dust/dust_cuda.cu

@@ -25,17 +25,20 @@
   #include "../utils/gpu.hpp"
   #include "../utils/hydro_utilities.h"
 
-void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma)
+void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                 Real grain_radius)
 {
   int n_cells = nx * ny * nz;
   int ngrid   = (n_cells + TPB - 1) / TPB;
   dim3 dim1dGrid(ngrid, 1, 1);
   dim3 dim1dBlock(TPB, 1, 1);
-  hipLaunchKernelGGL(Dust_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, dt, gamma);
+  hipLaunchKernelGGL(Dust_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, dt, gamma,
+                     grain_radius);
   GPU_Error_Check();
 }
 
-__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma)
+__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                            Real grain_radius)
 {
   // get grid indices
   int n_cells = nx * ny * nz;
@@ -99,8 +102,8 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
     // if dual energy is turned on use temp from total internal energy
     temperature = temperature_init;
 
-    Real tau_sp =
-        Calc_Sputtering_Timescale(number_density, temperature) / TIME_UNIT;  // sputtering timescale, kyr (sim units)
+    Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) /
+                  TIME_UNIT;  // sputtering timescale, kyr (sim units)
 
     dd_dt = Calc_dd_dt(density_dust, tau_sp);  // rate of change in dust density at current timestep
     dd    = dd_dt * dt;                        // change in dust density at current timestep
@@ -126,17 +129,18 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
 }
 
 // McKinnon et al. (2017) sputtering timescale
-__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature)
+__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature, Real grain_radius)
 {
-  Real grain_radius  = 1;          // dust grain size in units of 0.1 micrometers
-  Real temperature_0 = 2e6;        // temp above which the sputtering rate is ~constant in K
-  Real omega         = 2.5;        // controls the low-temperature scaling of the sputtering rate
-  Real A             = 5.3618e15;  // 0.17 Gyr in s
+  Real a             = grain_radius;  // dust grain size in units of 0.1 micrometers
+  Real temperature_0 = 2e6;           // temp above which the sputtering rate is ~constant in K
+  Real omega         = 2.5;           // controls the low-temperature scaling of the sputtering rate
+  Real A             = 5.3618e15;     // 0.17 Gyr in s
 
   number_density /= (6e-4);  // gas number density in units of 10^-27 g/cm^3
 
   // sputtering timescale, s
-  Real tau_sp = A * (grain_radius / number_density) * (pow(temperature_0 / temperature, omega) + 1);
+  printf("%e\n", grain_radius);
+  Real tau_sp = A * (a / number_density) * (pow(temperature_0 / temperature, omega) + 1);
 
   return tau_sp;
 }

src/dust/dust_cuda.h

@@ -27,7 +27,8 @@
  * \param[in] dt Simulation timestep
  * \param[in] gamma Specific heat ratio
  */
-void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma);
+void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                 Real grain_radius);
 
 /*!
  * \brief Compute the change in dust density for a cell and update its value in dev_conserved.
@@ -42,8 +43,8 @@ void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n
  * \param[in] dt Simulation timestep
  * \param[in] gamma Specific heat ratio
  */
-__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt,
-                            Real gamma);
+__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                            Real grain_radius);
 
 /*!
  * \brief Compute the sputtering timescale based on a cell's density and temperature.
@@ -53,7 +54,7 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
  *
  * \return Real Sputtering timescale in seconds (McKinnon et al. 2017)
  */
-__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature);
+__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature, Real grain_radius);
 
 /*!
  * \brief Compute the rate of change in dust density based on the current dust density and sputtering timescale.

src/dust/dust_cuda_tests.cpp

@@ -28,9 +28,11 @@ TEST(tDUSTTestSputteringTimescale,
   Real YR_IN_S                     = 3.154e7;
   Real const k_test_number_density = 1;
   Real const k_test_temperature    = pow(10, 5.0);
+  Real const k_test_grain_radius   = 1;
   Real const k_fiducial_num        = 182565146.96398282;
 
-  Real test_num = Calc_Sputtering_Timescale(k_test_number_density, k_test_temperature) / YR_IN_S;  // yr
+  Real test_num =
+      Calc_Sputtering_Timescale(k_test_number_density, k_test_temperature, k_test_grain_radius) / YR_IN_S;  // yr
 
   double abs_diff;
   int64_t ulps_diff;

src/global/global.cpp

@@ -445,6 +445,12 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "skewersdir") == 0) {
     strncpy(parms->skewersdir, value, MAXLEN);
   #endif
+#endif
+#ifdef SCALAR
+  #ifdef DUST
+  } else if (strcmp(name, "grain_radius") == 0) {
+    parms->grain_radius = atoi(value);
+  #endif
 #endif
   } else if (!Is_Param_Valid(name)) {
     chprintf("WARNING: %s/%s: Unknown parameter/value pair!\n", name, value);

src/global/global.h

@@ -337,6 +337,11 @@ struct Parameters {
   char skewersdir[MAXLEN];
   #endif
 #endif
+#ifdef SCALAR
+  #ifdef DUST
+  Real grain_radius;
+  #endif
+#endif
 };
 
 /*! \fn void parse_params(char *param_file, struct Parameters * parms);

src/grid/grid3D.cpp

@@ -274,6 +274,12 @@ void Grid3D::Initialize(struct Parameters *P)
   H.OUTPUT_SCALE_FACOR = not(P->scale_outputs_file[0] == '\0');
 #endif
 
+#ifdef SCALAR
+  #ifdef DUST
+  H.grain_radius = P->grain_radius;
+  #endif
+#endif
+
   H.Output_Initial = true;
 }
 
@@ -518,7 +524,7 @@ Real Grid3D::Update_Hydro_Grid()
 
   #ifdef DUST
   // ==Apply dust from dust/dust_cuda.h==
-  Dust_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama);
+  Dust_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama, H.grain_radius);
   #endif  // DUST
 
 #endif  // CUDA

src/grid/grid3D.h

@@ -267,6 +267,12 @@ struct Header {
    *  \brief Flag set to true when all the data will  be written to file
    * (Restart File ) */
   bool Output_Complete_Data;
+
+#ifdef SCALAR
+  #ifdef DUST
+  Real grain_radius;
+  #endif
+#endif
 };
 
 /*! \class Grid3D

src/grid/initial_conditions.cpp

@@ -1415,9 +1415,11 @@ void Grid3D::Clouds()
             C.GasEnergy[id] = p_cl / (gama - 1.0);
 #endif  // DE
 
-#ifdef DUST
+#ifdef SCALAR
+  #ifdef DUST
             C.host[id + H.n_cells * grid_enum::dust_density] = rho_cl * 1e-2;
-#endif  // DUST
+  #endif  // DUST
+#endif
           }
         }
       }