Merge pull request #315 from alwinm/2023-july-21

Revert minor changes so 1-D examples can run and add debug functions

src/grid/boundary_conditions.cpp

@@ -151,12 +151,6 @@ int Grid3D::Check_Custom_Boundary(int *flags, struct parameters P)
   }
 
   for (int i = 0; i < 6; i++) {
-    if (flags[i] < 1 or flags[i] > 5) {
-      chprintf(
-          "Invalid boundary conditions. Must select between 1 (periodic), 2 "
-          "(reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n");
-      chexit(-1);
-    }
     if (flags[i] == 4) {
       /*custom boundaries*/
       return 1;

src/grid/initial_conditions.cpp

@@ -508,28 +508,28 @@ void Grid3D::Square_Wave(parameters const &P)
  *  \brief Initialize the grid with a Riemann problem. */
 void Grid3D::Riemann(parameters const &P)
 {
-  size_t const istart = H.n_ghost;
+  size_t const istart = H.n_ghost - 1;
   size_t const iend   = H.nx - H.n_ghost;
   size_t jstart, kstart, jend, kend;
   if (H.ny > 1) {
-    jstart = H.n_ghost;
+    jstart = H.n_ghost - 1;
     jend   = H.ny - H.n_ghost;
   } else {
     jstart = 0;
     jend   = H.ny;
   }
   if (H.nz > 1) {
-    kstart = H.n_ghost;
+    kstart = H.n_ghost - 1;
     kend   = H.nz - H.n_ghost;
   } else {
     kstart = 0;
     kend   = H.nz;
   }
 
   // set initial values of conserved variables
-  for (size_t k = kstart - 1; k < kend; k++) {
-    for (size_t j = jstart - 1; j < jend; j++) {
-      for (size_t i = istart - 1; i < iend; i++) {
+  for (size_t k = kstart; k < kend; k++) {
+    for (size_t j = jstart; j < jend; j++) {
+      for (size_t i = istart; i < iend; i++) {
         // get cell index
         size_t const id = i + j * H.nx + k * H.nx * H.ny;
 
@@ -540,6 +540,7 @@ void Grid3D::Riemann(parameters const &P)
 #ifdef MHD
         // Set the magnetic field including the rightmost ghost cell on the
         // left side which is really the left face of the first grid cell
+        // WARNING: Only correct in 3-D
         if (x_pos < P.diaph) {
           C.magnetic_x[id] = P.Bx_l;
           C.magnetic_y[id] = P.By_l;
@@ -1973,4 +1974,4 @@ void Grid3D::Orszag_Tang_Vortex()
     }
   }
 }
-#endif  // MHD
+#endif  // MHD

src/reconstruction/plmc_cuda.cu

@@ -29,8 +29,8 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_
   int xid, yid, zid;
   cuda_utilities::compute3DIndices(thread_id, nx, ny, xid, yid, zid);
 
-  // Thread guard to prevent overrun
-  if (xid < 1 or xid >= nx - 2 or yid < 1 or yid >= ny - 2 or zid < 1 or zid >= nz - 2) {
+  // Ensure that we are only operating on cells that will be used
+  if (reconstruction::Thread_Guard<2>(nx, ny, nz, xid, yid, zid)) {
     return;
   }
 

src/reconstruction/ppmc_cuda.cu

@@ -27,9 +27,7 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun
   int xid, yid, zid;
   cuda_utilities::compute3DIndices(thread_id, nx, ny, xid, yid, zid);
 
-  // Ensure that we are only operating on cells that will be used
-  if (size_t const min = 3, max = 3;
-      xid < min or xid >= nx - max or yid < min or yid >= ny - max or zid < min or zid >= nz - max) {
+  if (reconstruction::Thread_Guard<3>(nx, ny, nz, xid, yid, zid)) {
     return;
   }
 
@@ -548,8 +546,7 @@ __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bo
   cuda_utilities::compute3DIndices(thread_id, nx, ny, xid, yid, zid);
 
   // Ensure that we are only operating on cells that will be used
-  if (size_t const min = 3, max = 3;
-      xid < min or xid >= nx - max or yid < min or yid >= ny - max or zid < min or zid >= nz - max) {
+  if (reconstruction::Thread_Guard<3>(nx, ny, nz, xid, yid, zid)) {
     return;
   }
 

src/reconstruction/reconstruction.h

@@ -78,6 +78,47 @@ struct Characteristic {
 };
 // =====================================================================================================================
 
+// =====================================================================================================================
+/*!
+ * \brief Determine if a thread is within the allowed range
+ *
+ * \tparam order The order of the reconstruction. 2 for PLM, 3 for PPM
+ * \param nx The number of cells in the X-direction
+ * \param ny The number of cells in the Y-direction
+ * \param nz The number of cells in the Z-direction
+ * \param xid The X thread index
+ * \param yid The Y thread index
+ * \param zid The Z thread index
+ * \return true The thread is NOT in the allowed range
+ * \return false The thread is in the allowed range
+ */
+template <int order>
+bool __device__ __host__ __inline__ Thread_Guard(int const &nx, int const &ny, int const &nz, int const &xid,
+                                                 int const &yid, int const &zid)
+{
+  // These checks all make sure that the xid is such that the thread won't try to load any memory that is out of bounds
+
+  // X check
+  bool out_of_bounds_thread = xid < order - 1 or xid >= nx - order;
+
+  // Y check, only used for 2D and 3D
+  if (ny > 1) {
+    out_of_bounds_thread = yid < order - 1 or yid >= ny - order or out_of_bounds_thread;
+  }
+
+  // z check, only used for 3D
+  if (nz > 1) {
+    out_of_bounds_thread = zid < order - 1 or zid >= nz - order or out_of_bounds_thread;
+  }
+  // This is needed in the case that nz == 1 to avoid overrun
+  else {
+    out_of_bounds_thread = zid >= nz or out_of_bounds_thread;
+  }
+
+  return out_of_bounds_thread;
+}
+// =====================================================================================================================
+
 // =====================================================================================================================
 /*!
  * \brief Load the data for reconstruction

src/reconstruction/reconstruction_tests.cu

@@ -18,6 +18,7 @@
 #include "../io/io.h"
 #include "../reconstruction/reconstruction.h"
 #include "../utils/DeviceVector.h"
+#include "../utils/cuda_utilities.h"
 #include "../utils/gpu.hpp"
 #include "../utils/testing_utilities.h"
 
@@ -174,6 +175,34 @@ TEST(tMHDReconstructionComputeEigenvectors, CorrectInputExpectCorrectOutput)
 }
 #endif  // MHD
 
+TEST(tALLReconstructionThreadGuard, CorrectInputExpectCorrectOutput)
+{
+  // Test parameters
+  int const order = 3;
+  int const nx    = 6;
+  int const ny    = 6;
+  int const nz    = 6;
+
+  // fiducial data
+  std::vector<int> fiducial_vals(nx * ny * nz, 1);
+  fiducial_vals.at(86) = 0;
+
+  // loop through all values of the indices and check them
+  for (int xid = 0; xid < nx; xid++) {
+    for (int yid = 0; yid < ny; yid++) {
+      for (int zid = 0; zid < nz; zid++) {
+        // Get the test value
+        bool test_val = reconstruction::Thread_Guard<order>(nx, ny, nz, xid, yid, zid);
+
+        // Compare
+        int id = cuda_utilities::compute1DIndex(xid, yid, zid, nx, ny);
+        ASSERT_EQ(test_val, fiducial_vals.at(id))
+            << "Test value not equal to fiducial value at id = " << id << std::endl;
+      }
+    }
+  }
+}
+
 TEST(tALLReconstructionLoadData, CorrectInputExpectCorrectOutput)
 {
   // Set up test and mock up grid

src/system_tests/hydro_system_tests.cpp

@@ -68,6 +68,18 @@ INSTANTIATE_TEST_SUITE_P(CorrectInputExpectCorrectOutput, tHYDROtMHDSYSTEMSodSho
 /// @}
 // =============================================================================
 
+TEST(tHYDROSYSTEMSodShockTube, OneDimensionalCorrectInputExpectCorrectOutput)
+{
+  systemTest::SystemTestRunner sodTest;
+  sodTest.runTest();
+}
+
+TEST(tHYDROSYSTEMSodShockTube, TwoDimensionalCorrectInputExpectCorrectOutput)
+{
+  systemTest::SystemTestRunner sodTest;
+  sodTest.runTest();
+}
+
 TEST(tHYDROtMHDSYSTEMConstant, CorrectInputExpectCorrectOutput)
 {
   systemTest::SystemTestRunner testObject(false, false, false);

src/system_tests/input_files/tHYDROSYSTEMSodShockTube_OneDimensionalCorrectInputExpectCorrectOutput.txt

@@ -0,0 +1,56 @@
+#
+# Parameter File for 1D Sod Shock tube
+#
+
+################################################
+# number of grid cells in the x dimension
+nx=64
+# number of grid cells in the y dimension
+ny=1
+# number of grid cells in the z dimension
+nz=1
+# final output time
+tout=0.2
+# time interval for output
+outstep=0.2
+# name of initial conditions
+init=Riemann
+# domain properties
+xmin=0.0
+ymin=0.0
+zmin=0.0
+xlen=1.0
+ylen=1.0
+zlen=1.0
+# type of boundary conditions
+xl_bcnd=3
+xu_bcnd=3
+yl_bcnd=3
+yu_bcnd=3
+zl_bcnd=3
+zu_bcnd=3
+# path to output directory
+outdir=./
+
+#################################################
+# Parameters for 1D Riemann problems
+# density of left state
+rho_l=1.0
+# velocity of left state
+vx_l=0.0
+vy_l=0.0
+vz_l=0.0
+# pressure of left state
+P_l=1.0
+# density of right state
+rho_r=0.1
+# velocity of right state
+vx_r=0.0
+vy_r=0.0
+vz_r=0.0
+# pressure of right state
+P_r=0.1
+# location of initial discontinuity
+diaph=0.5
+# value of gamma
+gamma=1.4

src/system_tests/input_files/tHYDROSYSTEMSodShockTube_TwoDimensionalCorrectInputExpectCorrectOutput.txt

@@ -0,0 +1,56 @@
+#
+# Parameter File for 1D Sod Shock tube
+#
+
+################################################
+# number of grid cells in the x dimension
+nx=64
+# number of grid cells in the y dimension
+ny=64
+# number of grid cells in the z dimension
+nz=1
+# final output time
+tout=0.2
+# time interval for output
+outstep=0.2
+# name of initial conditions
+init=Riemann
+# domain properties
+xmin=0.0
+ymin=0.0
+zmin=0.0
+xlen=1.0
+ylen=1.0
+zlen=1.0
+# type of boundary conditions
+xl_bcnd=3
+xu_bcnd=3
+yl_bcnd=3
+yu_bcnd=3
+zl_bcnd=3
+zu_bcnd=3
+# path to output directory
+outdir=./
+
+#################################################
+# Parameters for 1D Riemann problems
+# density of left state
+rho_l=1.0
+# velocity of left state
+vx_l=0.0
+vy_l=0.0
+vz_l=0.0
+# pressure of left state
+P_l=1.0
+# density of right state
+rho_r=0.1
+# velocity of right state
+vx_r=0.0
+vy_r=0.0
+vz_r=0.0
+# pressure of right state
+P_r=0.1
+# location of initial discontinuity
+diaph=0.5
+# value of gamma
+gamma=1.4

src/utils/debug_utilities.cu

@@ -0,0 +1,60 @@
+#include <math.h>
+
+#include "../global/global.h"
+#include "../global/global_cuda.h"
+#include "../io/io.h"                 // provides chprintf
+#include "../utils/error_handling.h"  // provides chexit
+
+__global__ void Dump_Values_Kernel(Real* device_array, int array_size, int marker)
+{
+  int tid = threadIdx.x + blockIdx.x * blockDim.x;
+  if (tid >= array_size) {
+    return;
+  }
+  kernel_printf("Dump Values: marker %d tid %d value %g \n", marker, tid, device_array[tid]);
+}
+
+/*
+  Prints out all values of a device_array
+ */
+void Dump_Values(Real* device_array, int array_size, int marker)
+{
+  int ngrid = (array_size + TPB - 1) / TPB;
+  dim3 dim1dGrid(ngrid, 1, 1);
+  dim3 dim1dBlock(TPB, 1, 1);
+  hipLaunchKernelGGL(Dump_Values_Kernel, dim1dGrid, dim1dBlock, 0, 0, device_array, array_size, marker);
+}
+
+__global__ void Check_For_Nan_Kernel(Real* device_array, int array_size, int check_num, bool* out_bool)
+{
+  int tid = threadIdx.x + blockIdx.x * blockDim.x;
+  if (tid >= array_size) {
+    return;
+  }
+  if (device_array[tid] == device_array[tid]) {
+    return;
+  }
+  out_bool[0] = true;
+  kernel_printf("Check_For_Nan_Kernel found Nan Checknum: %d Thread: %d\n", check_num, tid);
+}
+
+/*
+  Checks a device_array for NaN and prints/exits if found
+ */
+void Check_For_Nan(Real* device_array, int array_size, int check_num)
+{
+  bool host_out_bool[1] = {false};
+  bool* out_bool;
+  cudaMalloc((void**)&out_bool, sizeof(bool));
+  cudaMemcpy(out_bool, host_out_bool, sizeof(bool), cudaMemcpyHostToDevice);
+  int ngrid = (array_size + TPB - 1) / TPB;
+  dim3 dim1dGrid(ngrid, 1, 1);
+  dim3 dim1dBlock(TPB, 1, 1);
+  hipLaunchKernelGGL(Check_For_Nan_Kernel, dim1dGrid, dim1dBlock, 0, 0, device_array, array_size, check_num, out_bool);
+  cudaMemcpy(host_out_bool, out_bool, sizeof(bool), cudaMemcpyDeviceToHost);
+  cudaFree(out_bool);
+
+  if (host_out_bool[0]) {
+    chexit(-1);
+  }
+}