More thread crash improvements

src/grid/grid3D.cpp

@@ -8,15 +8,17 @@
 #endif
 #include "../global/global.h"
 #include "../grid/grid3D.h"
-#include "../grid/grid_enum.h"    // provides grid_enum
-#include "../hydro/hydro_cuda.h"  // provides Calc_dt_GPU
+#include "../grid/grid_enum.h"       // provides grid_enum
+#include "../hydro/average_cells.h"  // provides Average_Slow_Cells and SlowCellConditionChecker
+#include "../hydro/hydro_cuda.h"     // provides Calc_dt_GPU
 #include "../integrators/VL_1D_cuda.h"
 #include "../integrators/VL_2D_cuda.h"
 #include "../integrators/VL_3D_cuda.h"
 #include "../integrators/simple_1D_cuda.h"
 #include "../integrators/simple_2D_cuda.h"
 #include "../integrators/simple_3D_cuda.h"
 #include "../io/io.h"
+#include "../utils/DeviceVector.h"
 #include "../utils/error_handling.h"
 #ifdef MPI_CHOLLA
   #include <mpi.h>
@@ -152,7 +154,9 @@ void Grid3D::Initialize(struct Parameters *P)
 
 #ifdef AVERAGE_SLOW_CELLS
   H.min_dt_slow = 1e-100;  // Initialize the minumum dt to a tiny number
-#endif                     // AVERAGE_SLOW_CELLS
+#else
+  H.min_dt_slow = -1.0;
+#endif  // AVERAGE_SLOW_CELLS
 
 #ifndef MPI_CHOLLA
 
@@ -422,6 +426,9 @@ void Grid3D::Execute_Hydro_Integrator(void)
   Timer.Hydro_Integrator.Start();
 #endif  // CPU_TIME
 
+  // this buffer holds 1 element that is initialized to 0
+  cuda_utilities::DeviceVector<int> error_code_buffer(1, true);
+
   // Run the hydro integrator on the grid
   if (H.nx > 1 && H.ny == 1 && H.nz == 1)  // 1D
   {
@@ -446,18 +453,24 @@ void Grid3D::Execute_Hydro_Integrator(void)
 #ifdef VL
     VL_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
                          H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, H.density_floor,
-                         C.Grav_potential);
+                         C.Grav_potential, SlowCellConditionChecker(1.0 / H.min_dt_slow, H.dx, H.dy, H.dz),
+                         error_code_buffer.data());
 #endif  // VL
 #ifdef SIMPLE
     Simple_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
                              H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, H.density_floor,
-                             C.Grav_potential);
+                             C.Grav_potential, SlowCellConditionChecker(1.0 / H.min_dt_slow, H.dx, H.dy, H.dz),
+                             error_code_buffer.data());
 #endif  // SIMPLE
   } else {
     chprintf("Error: Grid dimensions nx: %d  ny: %d  nz: %d  not supported.\n", H.nx, H.ny, H.nz);
     chexit(-1);
   }
 
+  if (error_code_buffer[0] != 0) {
+    CHOLLA_ERROR("An error occurred during the hydro calculation.");
+  }
+
 #ifdef CPU_TIME
   Timer.Hydro_Integrator.End(true);
 #endif  // CPU_TIME
@@ -578,8 +591,9 @@ Real Grid3D::Update_Hydro_Grid()
   ny_off = ny_local_start;
   nz_off = nz_local_start;
   #endif
-  Average_Slow_Cells(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dx, H.dy, H.dz, gama, max_dti_slow, H.xbound,
-                     H.ybound, H.zbound, nx_off, ny_off, nz_off);
+  Average_Slow_Cells(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, gama,
+                     SlowCellConditionChecker(max_dti_slow, H.dx, H.dy, H.dz), H.xbound, H.ybound, H.zbound, nx_off,
+                     ny_off, nz_off);
 #endif  // AVERAGE_SLOW_CELLS
 
   // ==Calculate the next time step using Calc_dt_GPU from hydro/hydro_cuda.h==

src/grid/grid3D.h

@@ -210,9 +210,10 @@ struct Header {
    *  \brief Length of the current timestep */
   Real dt;
 
-#ifdef AVERAGE_SLOW_CELLS
+  /*! \brief Cells that introduce timesteps shorter than will be averaged with
+   *         neighboring cells. Should be a negative value when the
+   *         AVERAGE_SLOW_CELLS macro isn't defined. */
   Real min_dt_slow;
-#endif
 
   /*! \var t_wall
    *  \brief Wall time */

src/hydro/average_cells.h

@@ -0,0 +1,56 @@
+/*! \file average_cells.h
+ *  \brief Definitions of functions and classes that implement logic related to averaging cells with
+ *         neighbors. */
+
+#ifndef AVERAGE_CELLS_H
+#define AVERAGE_CELLS_H
+
+#include <math.h>
+
+#include "../global/global.h"
+
+/*! \brief Object that checks whether a given cell meets the conditions for slow-cell averaging.
+ *          The main motivation for creating this class is reducing ifdef statements (and allow to modify the
+ *          actual slow-cell-condition. */
+struct SlowCellConditionChecker {
+// omit data-members if they aren't used for anything
+#ifdef AVERAGE_SLOW_CELLS
+  Real max_dti_slow, dx, dy, dz;
+#endif
+
+  /*! \brief Construct a new object. */
+  __host__ __device__ SlowCellConditionChecker(Real max_dti_slow, Real dx, Real dy, Real dz)
+#ifdef AVERAGE_SLOW_CELLS
+      : max_dti_slow{max_dti_slow}, dx{dx}, dy{dy}, dz{dz}
+#endif
+  {
+  }
+
+  /*! \brief Returns whether the cell meets the condition for being considered a slow cell that must
+   *  be averaged. */
+  template <bool verbose = false>
+  __device__ bool is_slow(Real total_energy, Real density, Real density_inv, Real velocity_x, Real velocity_y,
+                          Real velocity_z, Real gamma) const
+  {
+    return this->max_dti_if_slow(total_energy, density, density_inv, velocity_x, velocity_y, velocity_z, gamma) > 0.0;
+  }
+
+  /*! \brief Returns the max inverse timestep of the specified cell, if it meets the criteria for being
+   *  a slow cell. If it doesn't, return a negative value instead.
+   */
+  __device__ Real max_dti_if_slow(Real total_energy, Real density, Real density_inv, Real velocity_x, Real velocity_y,
+                                  Real velocity_z, Real gamma) const;
+};
+
+#ifdef AVERAGE_SLOW_CELLS
+
+void Average_Slow_Cells(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma,
+                        SlowCellConditionChecker slow_check, Real xbound, Real ybound, Real zbound, int nx_offset,
+                        int ny_offset, int nz_offset);
+
+__global__ void Average_Slow_Cells_3D(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
+                                      Real gamma, SlowCellConditionChecker slow_check, Real xbound, Real ybound,
+                                      Real zbound, int nx_offset, int ny_offset, int nz_offset);
+#endif
+
+#endif /* AVERAGE_CELLS_H */