Merge branch 'dev' into dev-scalar-floor

builds/make.type.disk

@@ -21,7 +21,6 @@ DFLAGS += -DGRAVITY_5_POINTS_GRADIENT
 
 #DFLAGS += -DSTATIC_GRAV
 
-#DFLAGS    += -DOUTPUT_ALWAYS
 DFLAGS    += -DCUDA
 DFLAGS    += -DMPI_CHOLLA
 DFLAGS    += -DPRECISION=2

builds/make.type.mhd

@@ -51,6 +51,3 @@ DFLAGS    += $(MPI_GPU)
 
 # Limit the number of steps to evolve.
 # DFLAGS += -DN_STEPS_LIMIT=1000
-
-# Output on every time step
-# DFLAGS += -DOUTPUT_ALWAYS

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

@@ -27,9 +27,11 @@ TEST(tDUSTTestSputteringTimescale, CorrectInputExpectCorrectOutput)
   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

@@ -240,6 +240,10 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "out_float32_GasEnergy") == 0) {
     parms->out_float32_GasEnergy = atoi(value);
 #endif  // DE
+  } else if (strcmp(name, "output_always") == 0) {
+    int tmp = atoi(value);
+    CHOLLA_ASSERT((tmp == 0) or (tmp == 1), "output_always must be 1 or 0.");
+    parms->output_always = tmp;
 #ifdef MHD
   } else if (strcmp(name, "out_float32_magnetic_x") == 0) {
     parms->out_float32_magnetic_x = atoi(value);
@@ -453,6 +457,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

@@ -204,6 +204,7 @@ struct Parameters {
 #ifdef DE
   int out_float32_GasEnergy = 0;
 #endif
+  bool output_always = false;
 #ifdef STATIC_GRAV
   int custom_grav = 0;  // flag to set specific static gravity field
 #endif
@@ -341,6 +342,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

@@ -280,6 +280,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;
 }
 
@@ -524,7 +530,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

@@ -268,6 +268,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/main.cpp

@@ -325,9 +325,7 @@ int main(int argc, char *argv[])
         "%9.3f ms   total time = %9.4f s\n\n",
         G.H.n_step, G.H.t, G.H.dt, (stop_step - start_step) * 1000, G.H.t_wall);
 
-#ifdef OUTPUT_ALWAYS
-    G.H.Output_Now = true;
-#endif
+    if (P.output_always) G.H.Output_Now = true;
 
 #ifdef ANALYSIS
     if (G.Analysis.Output_Now) {

src/reconstruction/plmc_cuda.cu

@@ -124,6 +124,9 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_
   reconstruction::Primitive interface_L_iph = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, 1.0);
   reconstruction::Primitive interface_R_imh = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, -1.0);
 
+  // Limit the interfaces
+  reconstruction::Plm_Limit_Interfaces(interface_L_iph, interface_R_imh, cell_imo, cell_i, cell_ipo);
+
 #ifndef VL
 
   Real const dtodx = dt / dx;

src/reconstruction/plmc_cuda_tests.cu

@@ -204,17 +204,17 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
                                                                            {{21, 0.73640639402573249},
                                                                             {85, 3.3462413154443715},
                                                                             {149, 2.1945584994458125},
-                                                                            {213, 0.67418839414138987},
-                                                                            {277, 16.909618487528142},
-                                                                            {341, 2.1533768050263267},
-                                                                            {405, 1.6994195863331925}},
+                                                                            {213, 1.1837630990406585},
+                                                                            {277, 17.570011907061254},
+                                                                            {341, 2.1583975283044725},
+                                                                            {405, 1.7033818819502551}},
                                                                            {{21, 0.25340904981266843},
                                                                             {85, 2.0441984720128734},
                                                                             {149, 1.9959059157695584},
                                                                             {213, 0.45377591914009824},
-                                                                            {277, 23.677832869261188},
-                                                                            {341, 1.5437923271692418},
-                                                                            {405, 1.8141353672443383}}};
+                                                                            {277, 24.018953780483471},
+                                                                            {341, 1.7033818819502551},
+                                                                            {405, 1.8587936590169301}}};
   std::vector<std::unordered_map<int, double>> fiducial_interface_right = {{{20, 0.59023012197434721},
                                                                             {84, 3.0043379408547275},
                                                                             {148, 2.6320759184913625},
@@ -227,18 +227,18 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
                                                                                {81, 2.5027813113931279},
                                                                                {145, 2.6371119205792346},
                                                                                {209, 1.0210845222961809},
-                                                                               {273, 21.360010722689488},
+                                                                               {273, 21.353253570231175},
                                                                                {337, 2.1634182515826184},
-                                                                               {401, 1.7073441775673177},
+                                                                               {401, 1.7033818819502551},
                                                                            },
                                                                            {
                                                                                {5, 0.92705119413602599},
                                                                                {69, 1.9592598982258778},
                                                                                {133, 0.96653490574340428},
                                                                                {197, 1.3203867992383289},
-                                                                               {261, 8.0057564947791793},
+                                                                               {261, 7.9217487636977353},
                                                                                {325, 1.8629714367312684},
-                                                                               {389, 1.9034519507895218},
+                                                                               {389, 1.8587936590169301},
                                                                            }};
 
   // Loop over different directions

src/reconstruction/reconstruction.h

@@ -673,6 +673,56 @@ Primitive __device__ __host__ __inline__ Calc_Interface_Linear(Primitive const &
 }
 // =====================================================================================================================
 
+// =====================================================================================================================
+/*!
+ * \brief Apply limiting the the primitive interfaces in PLM reconstructions
+ *
+ * \param[in,out] interface_L_iph The unlimited left plus 1/2 interface
+ * \param[in,out] interface_R_imh The unlimited right minus 1/2 interface
+ * \param[in] cell_imo The cell centered values at i-1
+ * \param[in] cell_i The cell centered values at i
+ * \param[in] cell_ipo The cell centered values at i+1
+ */
+void __device__ __host__ __inline__ Plm_Limit_Interfaces(Primitive &interface_L_iph, Primitive &interface_R_imh,
+                                                         Primitive const &cell_imo, Primitive const &cell_i,
+                                                         Primitive const &cell_ipo)
+{
+  auto limiter = [](Real &l_iph, Real &r_imh, Real const &val_imo, Real const &val_i, Real const &val_ipo) {
+    r_imh = fmax(fmin(val_i, val_imo), r_imh);
+    r_imh = fmin(fmax(val_i, val_imo), r_imh);
+    l_iph = fmax(fmin(val_i, val_ipo), l_iph);
+    l_iph = fmin(fmax(val_i, val_ipo), l_iph);
+  };
+
+  limiter(interface_L_iph.density, interface_R_imh.density, cell_imo.density, cell_i.density, cell_ipo.density);
+  limiter(interface_L_iph.velocity_x, interface_R_imh.velocity_x, cell_imo.velocity_x, cell_i.velocity_x,
+          cell_ipo.velocity_x);
+  limiter(interface_L_iph.velocity_y, interface_R_imh.velocity_y, cell_imo.velocity_y, cell_i.velocity_y,
+          cell_ipo.velocity_y);
+  limiter(interface_L_iph.velocity_z, interface_R_imh.velocity_z, cell_imo.velocity_z, cell_i.velocity_z,
+          cell_ipo.velocity_z);
+  limiter(interface_L_iph.pressure, interface_R_imh.pressure, cell_imo.pressure, cell_i.pressure, cell_ipo.pressure);
+
+#ifdef MHD
+  limiter(interface_L_iph.magnetic_y, interface_R_imh.magnetic_y, cell_imo.magnetic_y, cell_i.magnetic_y,
+          cell_ipo.magnetic_y);
+  limiter(interface_L_iph.magnetic_z, interface_R_imh.magnetic_z, cell_imo.magnetic_z, cell_i.magnetic_z,
+          cell_ipo.magnetic_z);
+#endif  // MHD
+
+#ifdef DE
+  limiter(interface_L_iph.gas_energy, interface_R_imh.gas_energy, cell_imo.gas_energy, cell_i.gas_energy,
+          cell_ipo.gas_energy);
+#endif  // DE
+#ifdef SCALAR
+  for (int i = 0; i < NSCALARS; i++) {
+    limiter(interface_L_iph.scalar[i], interface_R_imh.scalar[i], cell_imo.scalar[i], cell_i.scalar[i],
+            cell_ipo.scalar[i]);
+  }
+#endif  // SCALAR
+}
+// =====================================================================================================================
+
 // =====================================================================================================================
 /*!
  * \brief Compute the interface state for the CTU version fo the reconstructor from the slope and cell centered state