diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu index 7ca48a9fd..46273ef03 100644 --- a/src/dust/dust_cuda.cu +++ b/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; } diff --git a/src/dust/dust_cuda.h b/src/dust/dust_cuda.h index fb72007ac..212901e8a 100644 --- a/src/dust/dust_cuda.h +++ b/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. diff --git a/src/dust/dust_cuda_tests.cpp b/src/dust/dust_cuda_tests.cpp index a1357037a..fb0677edf 100644 --- a/src/dust/dust_cuda_tests.cpp +++ b/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; diff --git a/src/global/global.cpp b/src/global/global.cpp index 552e52c83..1b8dada22 100644 --- a/src/global/global.cpp +++ b/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); diff --git a/src/global/global.h b/src/global/global.h index b7b3d5643..87cc124a4 100644 --- a/src/global/global.h +++ b/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); diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp index d22e1242b..b315d7764 100644 --- a/src/grid/grid3D.cpp +++ b/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 diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index aff94c898..3319822fc 100644 --- a/src/grid/grid3D.h +++ b/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 diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index 89c42aa24..332898f85 100644 --- a/src/grid/initial_conditions.cpp +++ b/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 } } }