From e5568afebbea0651ad70dc161b5d3ac3348326ea Mon Sep 17 00:00:00 2001 From: Alwin Date: Fri, 21 Jul 2023 23:45:09 -0400 Subject: [PATCH 01/14] reconstruction diff --- src/reconstruction/ppmc_cuda.cu | 36 +++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index d2e9589e5..b2253da36 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -27,11 +27,43 @@ __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); + + int xs, xe, ys, ye, zs, ze; + switch (dir) { + case 0: + xs = 2; + xe = nx - 3; + ys = 0; + ye = ny; + zs = 0; + ze = nz; + break; + case 1: + xs = 0; + xe = nx; + ys = 2; + ye = ny - 3; + zs = 0; + ze = nz; + break; + case 2: + xs = 0; + xe = nx; + ys = 0; + ye = ny; + zs = 2; + ze = nz - 3; + break; + } + if (xid < xs || xid >= xe || yid < ys || yid >= ye || zid < zs || zid >= ze) return; + + /* // 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) { return; } + */ // Compute the total number of cells int const n_cells = nx * ny * nz; @@ -56,6 +88,10 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun break; } + + + + // load the 5-cell stencil into registers // cell i reconstruction::Primitive const cell_i = From 9e97283ee80b1ec1e014941242570e3d014321b2 Mon Sep 17 00:00:00 2001 From: Alwin Date: Thu, 27 Jul 2023 16:57:28 -0400 Subject: [PATCH 02/14] add debug_utilities --- src/utils/debug_utilities.cu | 55 ++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 src/utils/debug_utilities.cu diff --git a/src/utils/debug_utilities.cu b/src/utils/debug_utilities.cu new file mode 100644 index 000000000..98c0672ff --- /dev/null +++ b/src/utils/debug_utilities.cu @@ -0,0 +1,55 @@ +#include + +#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); + } +} + From 7a5874e6f9b00a11e4f4d0f4abcaa2ae864a7397 Mon Sep 17 00:00:00 2001 From: Alwin Date: Thu, 27 Jul 2023 17:06:53 -0400 Subject: [PATCH 03/14] more reversions --- src/grid/boundary_conditions.cpp | 5 ++- src/reconstruction/ppmc_cuda.cu | 53 +++++++++++++++----------------- src/utils/debug_utilities.cu | 1 - 3 files changed, 26 insertions(+), 33 deletions(-) diff --git a/src/grid/boundary_conditions.cpp b/src/grid/boundary_conditions.cpp index 06e7196af..c13327987 100644 --- a/src/grid/boundary_conditions.cpp +++ b/src/grid/boundary_conditions.cpp @@ -153,9 +153,8 @@ 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); + "WARNING: Possibly invalid boundary conditions for direction: %d flag: %d. Must select between 1 (periodic), 2 " + "(reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n", i, flags[i]); } if (flags[i] == 4) { /*custom boundaries*/ diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index b2253da36..aaec9895a 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -27,33 +27,32 @@ __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); - int xs, xe, ys, ye, zs, ze; switch (dir) { - case 0: - xs = 2; - xe = nx - 3; - ys = 0; - ye = ny; - zs = 0; - ze = nz; - break; - case 1: - xs = 0; - xe = nx; - ys = 2; - ye = ny - 3; - zs = 0; - ze = nz; - break; - case 2: - xs = 0; - xe = nx; - ys = 0; - ye = ny; - zs = 2; - ze = nz - 3; - break; + case 0: + xs = 2; + xe = nx - 3; + ys = 0; + ye = ny; + zs = 0; + ze = nz; + break; + case 1: + xs = 0; + xe = nx; + ys = 2; + ye = ny - 3; + zs = 0; + ze = nz; + break; + case 2: + xs = 0; + xe = nx; + ys = 0; + ye = ny; + zs = 2; + ze = nz - 3; + break; } if (xid < xs || xid >= xe || yid < ys || yid >= ye || zid < zs || zid >= ze) return; @@ -88,10 +87,6 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun break; } - - - - // load the 5-cell stencil into registers // cell i reconstruction::Primitive const cell_i = diff --git a/src/utils/debug_utilities.cu b/src/utils/debug_utilities.cu index 98c0672ff..1cb6214a5 100644 --- a/src/utils/debug_utilities.cu +++ b/src/utils/debug_utilities.cu @@ -52,4 +52,3 @@ void Check_For_Nan(Real* device_array, int array_size, int check_num) chexit(-1); } } - From 0caae5780d210d0e0d5905a98b22439717e1eef2 Mon Sep 17 00:00:00 2001 From: Alwin Date: Mon, 31 Jul 2023 20:38:50 -0400 Subject: [PATCH 04/14] format boundaries --- src/grid/boundary_conditions.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/grid/boundary_conditions.cpp b/src/grid/boundary_conditions.cpp index c13327987..d7d332e8a 100644 --- a/src/grid/boundary_conditions.cpp +++ b/src/grid/boundary_conditions.cpp @@ -151,11 +151,15 @@ int Grid3D::Check_Custom_Boundary(int *flags, struct parameters P) } for (int i = 0; i < 6; i++) { + /* Alwin: I am disabling this check because it is needlessly occurring every timestep. if (flags[i] < 1 or flags[i] > 5) { chprintf( - "WARNING: Possibly invalid boundary conditions for direction: %d flag: %d. Must select between 1 (periodic), 2 " - "(reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n", i, flags[i]); + "WARNING: Possibly invalid boundary conditions for direction: %d flag: %d. Must select between 1 (periodic), " + "2 " + "(reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n", + i, flags[i]); } + */ if (flags[i] == 4) { /*custom boundaries*/ return 1; From 20b5f10911b5bffc10b3d6b86c9ea5cd29341ad5 Mon Sep 17 00:00:00 2001 From: Alwin Date: Mon, 31 Jul 2023 20:57:12 -0400 Subject: [PATCH 05/14] initial condition revert --- src/grid/initial_conditions.cpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index eebfbb21a..4452f3dcc 100644 --- a/src/grid/initial_conditions.cpp +++ b/src/grid/initial_conditions.cpp @@ -508,9 +508,9 @@ 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 iend = H.nx - H.n_ghost; - size_t jstart, kstart, jend, kend; + int const istart = H.n_ghost; + int const iend = H.nx - H.n_ghost; + int jstart, kstart, jend, kend; if (H.ny > 1) { jstart = H.n_ghost; jend = H.ny - H.n_ghost; @@ -527,9 +527,9 @@ void Grid3D::Riemann(parameters const &P) } // 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 (int k = kstart - 1; k < kend; k++) { + for (int j = jstart - 1; j < jend; j++) { + for (int i = istart - 1; 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; @@ -582,12 +583,12 @@ void Grid3D::Riemann(parameters const &P) #endif // SCALAR #ifdef DE C.GasEnergy[id] = P.P_r / (gama - 1.0); -#endif // DE - } - } - } - } - } +#endif // DE + } // if diaph + } // if real + } // k + } // j + } // i } /*! \fn void Shu_Osher() @@ -1973,4 +1974,4 @@ void Grid3D::Orszag_Tang_Vortex() } } } -#endif // MHD \ No newline at end of file +#endif // MHD From 9e8602da921cca013f74ed7584157a7738413a88 Mon Sep 17 00:00:00 2001 From: Alwin Date: Tue, 1 Aug 2023 10:01:07 -0400 Subject: [PATCH 06/14] appease clang tidy and disable incorrect test --- src/reconstruction/ppmc_cuda_tests.cu | 2 ++ src/utils/debug_utilities.cu | 12 +++++++++--- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/reconstruction/ppmc_cuda_tests.cu b/src/reconstruction/ppmc_cuda_tests.cu index 7dd9b49e3..79b4aafac 100644 --- a/src/reconstruction/ppmc_cuda_tests.cu +++ b/src/reconstruction/ppmc_cuda_tests.cu @@ -24,6 +24,8 @@ TEST(tHYDROPpmcCTUReconstructor, CorrectInputExpectCorrectOutput) { + // Alwin: skip until this has been fixed + GTEST_SKIP(); // Set up PRNG to use std::mt19937_64 prng(42); std::uniform_real_distribution doubleRand(0.1, 5); diff --git a/src/utils/debug_utilities.cu b/src/utils/debug_utilities.cu index 1cb6214a5..9a1157aca 100644 --- a/src/utils/debug_utilities.cu +++ b/src/utils/debug_utilities.cu @@ -8,7 +8,9 @@ __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; + if (tid >= array_size) { + return; + } kernel_printf("Dump Values: marker %d tid %d value %g \n", marker, tid, device_array[tid]); } @@ -26,8 +28,12 @@ void Dump_Values(Real* device_array, int array_size, int 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; + 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); } From 5a550cc893ca7222d69ec5c63f4e807797876026 Mon Sep 17 00:00:00 2001 From: Alwin Date: Tue, 1 Aug 2023 10:20:43 -0400 Subject: [PATCH 07/14] clang tidy --- src/reconstruction/ppmc_cuda.cu | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index aaec9895a..efe3e51c1 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -54,7 +54,9 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun ze = nz - 3; break; } - if (xid < xs || xid >= xe || yid < ys || yid >= ye || zid < zs || zid >= ze) return; + if (xid < xs || xid >= xe || yid < ys || yid >= ye || zid < zs || zid >= ze) { + return; + } /* // Ensure that we are only operating on cells that will be used From 297225bd5c6d786d278ac652adf14a8c5901a5bd Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 10:21:20 -0400 Subject: [PATCH 08/14] Move check for boundary conditions type This previously was checked on every single time step. Now it is checked once at the beginning in the `Check_Configuration` function. --- src/grid/boundary_conditions.cpp | 9 --------- src/utils/error_handling.cpp | 14 ++++++++++++++ 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/grid/boundary_conditions.cpp b/src/grid/boundary_conditions.cpp index d7d332e8a..eca473fdb 100644 --- a/src/grid/boundary_conditions.cpp +++ b/src/grid/boundary_conditions.cpp @@ -151,15 +151,6 @@ int Grid3D::Check_Custom_Boundary(int *flags, struct parameters P) } for (int i = 0; i < 6; i++) { - /* Alwin: I am disabling this check because it is needlessly occurring every timestep. - if (flags[i] < 1 or flags[i] > 5) { - chprintf( - "WARNING: Possibly invalid boundary conditions for direction: %d flag: %d. Must select between 1 (periodic), " - "2 " - "(reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n", - i, flags[i]); - } - */ if (flags[i] == 4) { /*custom boundaries*/ return 1; diff --git a/src/utils/error_handling.cpp b/src/utils/error_handling.cpp index fd2a59ad7..38374c704 100644 --- a/src/utils/error_handling.cpp +++ b/src/utils/error_handling.cpp @@ -54,6 +54,20 @@ void Check_Configuration(parameters const &P) #error "Only one integrator can be enabled at a time." #endif // Only one integrator check + // Check the boundary conditions + auto Check_Boundary = [](int const &boundary) { + bool is_allowed_bc = boundary >= 0 and boundary <= 4; + assert(is_allowed_bc && + "WARNING: Possibly invalid boundary conditions for direction: %d flag: %d. Must select between (periodic), " + "2 (reflective), 3 (transmissive), 4 (custom), 5 (mpi).\n"); + }; + Check_Boundary(P.xl_bcnd); + Check_Boundary(P.xu_bcnd); + Check_Boundary(P.yl_bcnd); + Check_Boundary(P.yu_bcnd); + Check_Boundary(P.zl_bcnd); + Check_Boundary(P.zu_bcnd); + // warn if error checking is disabled #ifndef CUDA_ERROR_CHECK #warning "CUDA error checking is disabled. Enable it with the CUDA_ERROR_CHECK macro" From df7df09ab36ce0757abec64009eb05af275dbeef Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 10:26:32 -0400 Subject: [PATCH 09/14] Remove commented out code --- src/grid/initial_conditions.cpp | 12 ++++++------ src/reconstruction/ppmc_cuda.cu | 8 -------- 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index 4452f3dcc..26c4cf3ca 100644 --- a/src/grid/initial_conditions.cpp +++ b/src/grid/initial_conditions.cpp @@ -583,12 +583,12 @@ void Grid3D::Riemann(parameters const &P) #endif // SCALAR #ifdef DE C.GasEnergy[id] = P.P_r / (gama - 1.0); -#endif // DE - } // if diaph - } // if real - } // k - } // j - } // i +#endif // DE + } + } + } + } + } } /*! \fn void Shu_Osher() diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index efe3e51c1..ed4af9daa 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -58,14 +58,6 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun return; } - /* - // 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) { - return; - } - */ - // Compute the total number of cells int const n_cells = nx * ny * nz; From a6ace6b2c5605a3978228af4e0946146a3a36af6 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 10:33:46 -0400 Subject: [PATCH 10/14] Replace `int` indices with `size_t` in Riemann IC --- src/grid/initial_conditions.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index 26c4cf3ca..38967e4b7 100644 --- a/src/grid/initial_conditions.cpp +++ b/src/grid/initial_conditions.cpp @@ -508,18 +508,18 @@ void Grid3D::Square_Wave(parameters const &P) * \brief Initialize the grid with a Riemann problem. */ void Grid3D::Riemann(parameters const &P) { - int const istart = H.n_ghost; - int const iend = H.nx - H.n_ghost; - int jstart, kstart, jend, kend; + 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; @@ -527,9 +527,9 @@ void Grid3D::Riemann(parameters const &P) } // set initial values of conserved variables - for (int k = kstart - 1; k < kend; k++) { - for (int j = jstart - 1; j < jend; j++) { - for (int 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; From fe78156530da3dbbda905a60a9df138773e91f3a Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 11:45:44 -0400 Subject: [PATCH 11/14] Add 1D & 2D Sod tests --- cholla-tests-data | 2 +- src/system_tests/hydro_system_tests.cpp | 12 ++++ ...nsionalCorrectInputExpectCorrectOutput.txt | 56 +++++++++++++++++++ ...nsionalCorrectInputExpectCorrectOutput.txt | 56 +++++++++++++++++++ 4 files changed, 125 insertions(+), 1 deletion(-) create mode 100644 src/system_tests/input_files/tHYDROSYSTEMSodShockTube_OneDimensionalCorrectInputExpectCorrectOutput.txt create mode 100644 src/system_tests/input_files/tHYDROSYSTEMSodShockTube_TwoDimensionalCorrectInputExpectCorrectOutput.txt diff --git a/cholla-tests-data b/cholla-tests-data index 321416680..dcd73ff52 160000 --- a/cholla-tests-data +++ b/cholla-tests-data @@ -1 +1 @@ -Subproject commit 321416680f95d97b5d4ccc6f0b83a8b9ecafdaf0 +Subproject commit dcd73ff52b9027627b247c6d888bcdb56840c85e diff --git a/src/system_tests/hydro_system_tests.cpp b/src/system_tests/hydro_system_tests.cpp index 292935813..288690290 100644 --- a/src/system_tests/hydro_system_tests.cpp +++ b/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); diff --git a/src/system_tests/input_files/tHYDROSYSTEMSodShockTube_OneDimensionalCorrectInputExpectCorrectOutput.txt b/src/system_tests/input_files/tHYDROSYSTEMSodShockTube_OneDimensionalCorrectInputExpectCorrectOutput.txt new file mode 100644 index 000000000..dd54ff082 --- /dev/null +++ b/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 diff --git a/src/system_tests/input_files/tHYDROSYSTEMSodShockTube_TwoDimensionalCorrectInputExpectCorrectOutput.txt b/src/system_tests/input_files/tHYDROSYSTEMSodShockTube_TwoDimensionalCorrectInputExpectCorrectOutput.txt new file mode 100644 index 000000000..c89e179be --- /dev/null +++ b/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 From 11804ba2362387f67143ca0f1f135b1a69856a8d Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 14:49:58 -0400 Subject: [PATCH 12/14] Add new reconstruction::Thread_Guard function This function determines if a thread is in the domain to be operated on and returns the appropriate bool for a threadguard; true if it isn't, false if it is. - Enable PLMC test, for some reason it was commented out - Update PPMC tests for new threadguard --- src/reconstruction/plmc_cuda.cu | 4 +- src/reconstruction/plmc_cuda_tests.cu | 506 +++++++++++++------------- src/reconstruction/ppmc_cuda.cu | 32 +- src/reconstruction/ppmc_cuda_tests.cu | 194 +++++----- src/reconstruction/reconstruction.h | 23 ++ 5 files changed, 376 insertions(+), 383 deletions(-) diff --git a/src/reconstruction/plmc_cuda.cu b/src/reconstruction/plmc_cuda.cu index 8db428b82..41a5ae505 100644 --- a/src/reconstruction/plmc_cuda.cu +++ b/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; } diff --git a/src/reconstruction/plmc_cuda_tests.cu b/src/reconstruction/plmc_cuda_tests.cu index 11f859967..3616d2d0a 100644 --- a/src/reconstruction/plmc_cuda_tests.cu +++ b/src/reconstruction/plmc_cuda_tests.cu @@ -1,280 +1,280 @@ -// /*! -// * \file plmc_cuda_tests.cu -// * \brief Tests for the contents of plmc_cuda.h and plmc_cuda.cu -// * -// */ +/*! + * \file plmc_cuda_tests.cu + * \brief Tests for the contents of plmc_cuda.h and plmc_cuda.cu + * + */ -// // STL Includes -// #include -// #include -// #include -// #include +// STL Includes +#include +#include +#include +#include -// // External Includes -// #include // Include GoogleTest and related libraries/headers +// External Includes +#include // Include GoogleTest and related libraries/headers -// // Local Includes -// #include +// Local Includes +#include -// #include "../global/global.h" -// #include "../io/io.h" -// #include "../reconstruction/plmc_cuda.h" -// #include "../utils/DeviceVector.h" -// #include "../utils/hydro_utilities.h" -// #include "../utils/testing_utilities.h" +#include "../global/global.h" +#include "../io/io.h" +#include "../reconstruction/plmc_cuda.h" +#include "../utils/DeviceVector.h" +#include "../utils/hydro_utilities.h" +#include "../utils/testing_utilities.h" -// TEST(tHYDROPlmcReconstructor, CorrectInputExpectCorrectOutput) -// { -// // Set up PRNG to use -// std::mt19937_64 prng(42); -// std::uniform_real_distribution doubleRand(0.1, 5); +TEST(tHYDROPlmcReconstructor, CorrectInputExpectCorrectOutput) +{ + // Set up PRNG to use + std::mt19937_64 prng(42); + std::uniform_real_distribution doubleRand(0.1, 5); -// // Mock up needed information -// size_t const nx = 5; -// size_t const ny = 4; -// size_t const nz = 4; -// size_t const n_fields = 5; -// double const dx = doubleRand(prng); -// double const dt = doubleRand(prng); -// double const gamma = 5.0 / 3.0; + // Mock up needed information + size_t const nx = 5; + size_t const ny = 4; + size_t const nz = 4; + size_t const n_fields = 5; + double const dx = doubleRand(prng); + double const dt = doubleRand(prng); + double const gamma = 5.0 / 3.0; -// // Setup host grid. Fill host grid with random values and randomly assign maximum value -// std::vector host_grid(nx * ny * nz * n_fields); -// for (Real &val : host_grid) { -// val = doubleRand(prng); -// } + // Setup host grid. Fill host grid with random values and randomly assign maximum value + std::vector host_grid(nx * ny * nz * n_fields); + for (Real &val : host_grid) { + val = doubleRand(prng); + } -// // Allocating and copying to device -// cuda_utilities::DeviceVector dev_grid(host_grid.size()); -// dev_grid.cpyHostToDevice(host_grid); + // Allocating and copying to device + cuda_utilities::DeviceVector dev_grid(host_grid.size()); + dev_grid.cpyHostToDevice(host_grid); -// // Fiducial Data -// std::vector> fiducial_interface_left = {{{26, 2.1584359129984056}, -// {27, 0.70033864721549188}, -// {106, 2.2476363309467553}, -// {107, 3.0633780053857027}, -// {186, 2.2245934101106259}, -// {187, 2.1015872413794123}, -// {266, 2.1263341057778309}, -// {267, 3.9675148506537838}, -// {346, 3.3640057502842691}, -// {347, 21.091316282933843}}, -// {{21, 0.72430827309279655}, -// {37, 0.19457128219588618}, -// {101, 5.4739527659741896}, -// {117, 4.4286255636679313}, -// {181, 0.12703829036056602}, -// {197, 2.2851440769830953}, -// {261, 1.5337035731959561}, -// {277, 2.697375839048191}, -// {341, 22.319601655044117}, -// {357, 82.515887983144168}}, -// {{25, 2.2863650183226212}, -// {29, 1.686415421301841}, -// {105, 0.72340346106443465}, -// {109, 5.4713687086831388}, -// {185, 3.929100145230096}, -// {189, 4.9166140516911483}, -// {265, 0.95177493689267167}, -// {269, 0.46056494878491938}, -// {345, 3.6886096301452787}, -// {349, 16.105488797582133}}}; -// std::vector> fiducial_interface_right = {{{25, 3.8877922383184833}, -// {26, 0.70033864721549188}, -// {105, 1.5947787943675635}, -// {106, 3.0633780053857027}, -// {185, 4.0069556576401011}, -// {186, 2.1015872413794123}, -// {265, 1.7883678016935785}, -// {266, 3.9675148506537838}, -// {345, 2.8032969746372527}, -// {346, 21.091316282933843}}, -// {{17, 0.43265217076853835}, -// {33, 0.19457128219588618}, -// {97, 3.2697645945288754}, -// {113, 4.4286255636679313}, -// {177, 0.07588397666718491}, -// {193, 2.2851440769830953}, -// {257, 0.91612950577699748}, -// {273, 2.697375839048191}, -// {337, 13.332201861384396}, -// {353, 82.515887983144168}}, -// {{5, 2.2863650183226212}, -// {9, 1.686415421301841}, -// {85, 0.72340346106443465}, -// {89, 1.7792505446336098}, -// {165, 5.3997753452111859}, -// {169, 1.4379190463124139}, -// {245, 0.95177493689267167}, -// {249, 0.46056494878491938}, -// {325, 6.6889498465051407}, -// {329, 1.6145084086614281}}}; + // Fiducial Data + std::vector> fiducial_interface_left = {{{26, 2.1584359129984056}, + {27, 0.70033864721549188}, + {106, 2.2476363309467553}, + {107, 3.0633780053857027}, + {186, 2.2245934101106259}, + {187, 2.1015872413794123}, + {266, 2.1263341057778309}, + {267, 3.9675148506537838}, + {346, 3.3640057502842691}, + {347, 21.091316282933843}}, + {{21, 0.72430827309279655}, + {37, 0.19457128219588618}, + {101, 5.4739527659741896}, + {117, 4.4286255636679313}, + {181, 0.12703829036056602}, + {197, 2.2851440769830953}, + {261, 1.5337035731959561}, + {277, 2.697375839048191}, + {341, 22.319601655044117}, + {357, 82.515887983144168}}, + {{25, 2.2863650183226212}, + {29, 1.686415421301841}, + {105, 0.72340346106443465}, + {109, 5.4713687086831388}, + {185, 3.929100145230096}, + {189, 4.9166140516911483}, + {265, 0.95177493689267167}, + {269, 0.46056494878491938}, + {345, 3.6886096301452787}, + {349, 16.105488797582133}}}; + std::vector> fiducial_interface_right = {{{25, 3.8877922383184833}, + {26, 0.70033864721549188}, + {105, 1.5947787943675635}, + {106, 3.0633780053857027}, + {185, 4.0069556576401011}, + {186, 2.1015872413794123}, + {265, 1.7883678016935785}, + {266, 3.9675148506537838}, + {345, 2.8032969746372527}, + {346, 21.091316282933843}}, + {{17, 0.43265217076853835}, + {33, 0.19457128219588618}, + {97, 3.2697645945288754}, + {113, 4.4286255636679313}, + {177, 0.07588397666718491}, + {193, 2.2851440769830953}, + {257, 0.91612950577699748}, + {273, 2.697375839048191}, + {337, 13.332201861384396}, + {353, 82.515887983144168}}, + {{5, 2.2863650183226212}, + {9, 1.686415421301841}, + {85, 0.72340346106443465}, + {89, 1.7792505446336098}, + {165, 5.3997753452111859}, + {169, 1.4379190463124139}, + {245, 0.95177493689267167}, + {249, 0.46056494878491938}, + {325, 6.6889498465051407}, + {329, 1.6145084086614281}}}; -// // Loop over different directions -// for (size_t direction = 0; direction < 3; direction++) { -// // Assign the shape -// size_t nx_rot, ny_rot, nz_rot; -// switch (direction) { -// case 0: -// nx_rot = nx; -// ny_rot = ny; -// nz_rot = nz; -// break; -// case 1: -// nx_rot = ny; -// ny_rot = nz; -// nz_rot = nx; -// break; -// case 2: -// nx_rot = nz; -// ny_rot = nx; -// nz_rot = ny; -// break; -// } + // Loop over different directions + for (size_t direction = 0; direction < 3; direction++) { + // Assign the shape + size_t nx_rot, ny_rot, nz_rot; + switch (direction) { + case 0: + nx_rot = nx; + ny_rot = ny; + nz_rot = nz; + break; + case 1: + nx_rot = ny; + ny_rot = nz; + nz_rot = nx; + break; + case 2: + nx_rot = nz; + ny_rot = nx; + nz_rot = ny; + break; + } -// // Allocate device buffers -// cuda_utilities::DeviceVector dev_interface_left(host_grid.size(), true); -// cuda_utilities::DeviceVector dev_interface_right(host_grid.size(), true); + // Allocate device buffers + cuda_utilities::DeviceVector dev_interface_left(host_grid.size(), true); + cuda_utilities::DeviceVector dev_interface_right(host_grid.size(), true); -// // Launch kernel -// hipLaunchKernelGGL(PLMC_cuda, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), -// dev_interface_right.data(), nx_rot, ny_rot, nz_rot, dx, dt, gamma, direction, n_fields); -// CudaCheckError(); -// CHECK(cudaDeviceSynchronize()); + // Launch kernel + hipLaunchKernelGGL(PLMC_cuda, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx_rot, ny_rot, nz_rot, dx, dt, gamma, direction, n_fields); + CudaCheckError(); + CHECK(cudaDeviceSynchronize()); -// // Perform Comparison -// for (size_t i = 0; i < host_grid.size(); i++) { -// // Check the left interface -// double test_val = dev_interface_left.at(i); -// double fiducial_val = -// (fiducial_interface_left.at(direction).find(i) == fiducial_interface_left.at(direction).end()) -// ? 0.0 -// : fiducial_interface_left.at(direction)[i]; + // Perform Comparison + for (size_t i = 0; i < host_grid.size(); i++) { + // Check the left interface + double test_val = dev_interface_left.at(i); + double fiducial_val = + (fiducial_interface_left.at(direction).find(i) == fiducial_interface_left.at(direction).end()) + ? 0.0 + : fiducial_interface_left.at(direction)[i]; -// testingUtilities::checkResults( -// fiducial_val, test_val, -// "left interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); + testingUtilities::checkResults( + fiducial_val, test_val, + "left interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); -// // Check the right interface -// test_val = dev_interface_right.at(i); -// fiducial_val = (fiducial_interface_right.at(direction).find(i) == fiducial_interface_right.at(direction).end()) -// ? 0.0 -// : fiducial_interface_right.at(direction)[i]; + // Check the right interface + test_val = dev_interface_right.at(i); + fiducial_val = (fiducial_interface_right.at(direction).find(i) == fiducial_interface_right.at(direction).end()) + ? 0.0 + : fiducial_interface_right.at(direction)[i]; -// testingUtilities::checkResults( -// fiducial_val, test_val, -// "right interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); -// } -// } -// } + testingUtilities::checkResults( + fiducial_val, test_val, + "right interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); + } + } +} -// TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput) -// { -// // Set up PRNG to use -// std::mt19937_64 prng(42); -// std::uniform_real_distribution doubleRand(0.1, 5); +TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput) +{ + // Set up PRNG to use + std::mt19937_64 prng(42); + std::uniform_real_distribution doubleRand(0.1, 5); -// // Mock up needed information -// size_t const nx = 4, ny = nx, nz = nx; -// size_t const n_fields = 8; -// size_t const n_cells_grid = nx * ny * nz * n_fields; -// size_t const n_cells_interface = nx * ny * nz * (n_fields - 1); -// double const dx = doubleRand(prng); -// double const dt = doubleRand(prng); -// double const gamma = 5.0 / 3.0; + // Mock up needed information + size_t const nx = 4, ny = nx, nz = nx; + size_t const n_fields = 8; + size_t const n_cells_grid = nx * ny * nz * n_fields; + size_t const n_cells_interface = nx * ny * nz * (n_fields - 1); + double const dx = doubleRand(prng); + double const dt = doubleRand(prng); + double const gamma = 5.0 / 3.0; -// // Setup host grid. Fill host grid with random values and randomly assign maximum value -// std::vector host_grid(n_cells_grid); -// for (Real &val : host_grid) { -// val = doubleRand(prng); -// } + // Setup host grid. Fill host grid with random values and randomly assign maximum value + std::vector host_grid(n_cells_grid); + for (Real &val : host_grid) { + val = doubleRand(prng); + } -// // Allocating and copying to device -// cuda_utilities::DeviceVector dev_grid(host_grid.size()); -// dev_grid.cpyHostToDevice(host_grid); + // Allocating and copying to device + cuda_utilities::DeviceVector dev_grid(host_grid.size()); + dev_grid.cpyHostToDevice(host_grid); -// // Fiducial Data -// std::vector> fiducial_interface_left = {{{21, 0.59023012197434721}, -// {85, 3.0043379408547275}, -// {149, 2.6320759184913625}, -// {213, 0.9487867623146744}, -// {277, 18.551193003661723}, -// {341, 1.8587936590169301}, -// {405, 2.1583975283044725}}, -// {{21, 0.73640639402573249}, -// {85, 3.3462413154443715}, -// {149, 2.1945584994458125}, -// {213, 0.67418839414138987}, -// {277, 16.909618487528142}, -// {341, 2.1533768050263267}, -// {405, 1.6994195863331925}}, -// {{21, 0.25340904981266843}, -// {85, 2.0441984720128734}, -// {149, 1.9959059157695584}, -// {213, 0.45377591914009824}, -// {277, 23.677832869261188}, -// {341, 1.5437923271692418}, -// {405, 1.8141353672443383}}}; -// std::vector> fiducial_interface_right = {{{20, 0.59023012197434721}, -// {84, 3.0043379408547275}, -// {148, 2.6320759184913625}, -// {212, 0.9487867623146744}, -// {276, 22.111134849009044}, -// {340, 1.8587936590169301}, -// {404, 2.1583975283044725}}, -// { -// {17, 0.44405384992296193}, -// {81, 2.5027813113931279}, -// {145, 2.6371119205792346}, -// {209, 1.0210845222961809}, -// {273, 21.360010722689488}, -// {337, 2.1634182515826184}, -// {401, 1.7073441775673177}, -// }, -// { -// {5, 0.92705119413602599}, -// {69, 1.9592598982258778}, -// {133, 0.96653490574340428}, -// {197, 1.3203867992383289}, -// {261, 8.0057564947791793}, -// {325, 1.8629714367312684}, -// {389, 1.9034519507895218}, -// }}; + // Fiducial Data + std::vector> fiducial_interface_left = {{{21, 0.59023012197434721}, + {85, 3.0043379408547275}, + {149, 2.6320759184913625}, + {213, 0.9487867623146744}, + {277, 18.551193003661723}, + {341, 1.8587936590169301}, + {405, 2.1583975283044725}}, + {{21, 0.73640639402573249}, + {85, 3.3462413154443715}, + {149, 2.1945584994458125}, + {213, 0.67418839414138987}, + {277, 16.909618487528142}, + {341, 2.1533768050263267}, + {405, 1.6994195863331925}}, + {{21, 0.25340904981266843}, + {85, 2.0441984720128734}, + {149, 1.9959059157695584}, + {213, 0.45377591914009824}, + {277, 23.677832869261188}, + {341, 1.5437923271692418}, + {405, 1.8141353672443383}}}; + std::vector> fiducial_interface_right = {{{20, 0.59023012197434721}, + {84, 3.0043379408547275}, + {148, 2.6320759184913625}, + {212, 0.9487867623146744}, + {276, 22.111134849009044}, + {340, 1.8587936590169301}, + {404, 2.1583975283044725}}, + { + {17, 0.44405384992296193}, + {81, 2.5027813113931279}, + {145, 2.6371119205792346}, + {209, 1.0210845222961809}, + {273, 21.360010722689488}, + {337, 2.1634182515826184}, + {401, 1.7073441775673177}, + }, + { + {5, 0.92705119413602599}, + {69, 1.9592598982258778}, + {133, 0.96653490574340428}, + {197, 1.3203867992383289}, + {261, 8.0057564947791793}, + {325, 1.8629714367312684}, + {389, 1.9034519507895218}, + }}; -// // Loop over different directions -// for (size_t direction = 0; direction < 3; direction++) { -// // Allocate device buffers -// cuda_utilities::DeviceVector dev_interface_left(n_cells_interface, true); -// cuda_utilities::DeviceVector dev_interface_right(n_cells_interface, true); + // Loop over different directions + for (size_t direction = 0; direction < 3; direction++) { + // Allocate device buffers + cuda_utilities::DeviceVector dev_interface_left(n_cells_interface, true); + cuda_utilities::DeviceVector dev_interface_right(n_cells_interface, true); -// // Launch kernel -// hipLaunchKernelGGL(PLMC_cuda, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), -// dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, direction, n_fields); -// CudaCheckError(); -// CHECK(cudaDeviceSynchronize()); + // Launch kernel + hipLaunchKernelGGL(PLMC_cuda, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, direction, n_fields); + CudaCheckError(); + CHECK(cudaDeviceSynchronize()); -// // Perform Comparison -// for (size_t i = 0; i < dev_interface_right.size(); i++) { -// // Check the left interface -// double test_val = dev_interface_left.at(i); -// double fiducial_val = -// (fiducial_interface_left.at(direction).find(i) == fiducial_interface_left.at(direction).end()) -// ? 0.0 -// : fiducial_interface_left.at(direction)[i]; + // Perform Comparison + for (size_t i = 0; i < dev_interface_right.size(); i++) { + // Check the left interface + double test_val = dev_interface_left.at(i); + double fiducial_val = + (fiducial_interface_left.at(direction).find(i) == fiducial_interface_left.at(direction).end()) + ? 0.0 + : fiducial_interface_left.at(direction)[i]; -// testingUtilities::checkResults( -// fiducial_val, test_val, -// "left interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); + testingUtilities::checkResults( + fiducial_val, test_val, + "left interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); -// // Check the right interface -// test_val = dev_interface_right.at(i); -// fiducial_val = (fiducial_interface_right.at(direction).find(i) == fiducial_interface_right.at(direction).end()) -// ? 0.0 -// : fiducial_interface_right.at(direction)[i]; + // Check the right interface + test_val = dev_interface_right.at(i); + fiducial_val = (fiducial_interface_right.at(direction).find(i) == fiducial_interface_right.at(direction).end()) + ? 0.0 + : fiducial_interface_right.at(direction)[i]; -// testingUtilities::checkResults( -// fiducial_val, test_val, -// "right interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); -// } -// } -// } + testingUtilities::checkResults( + fiducial_val, test_val, + "right interface at i=" + std::to_string(i) + ", in direction " + std::to_string(direction)); + } + } +} diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index ed4af9daa..4db993d70 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -27,34 +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); - int xs, xe, ys, ye, zs, ze; - switch (dir) { - case 0: - xs = 2; - xe = nx - 3; - ys = 0; - ye = ny; - zs = 0; - ze = nz; - break; - case 1: - xs = 0; - xe = nx; - ys = 2; - ye = ny - 3; - zs = 0; - ze = nz; - break; - case 2: - xs = 0; - xe = nx; - ys = 0; - ye = ny; - zs = 2; - ze = nz - 3; - break; - } - if (xid < xs || xid >= xe || yid < ys || yid >= ye || zid < zs || zid >= ze) { + if (reconstruction::Thread_Guard<3>(nx, ny, nz, xid, yid, zid)) { return; } @@ -573,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; } diff --git a/src/reconstruction/ppmc_cuda_tests.cu b/src/reconstruction/ppmc_cuda_tests.cu index 79b4aafac..1c7515ec0 100644 --- a/src/reconstruction/ppmc_cuda_tests.cu +++ b/src/reconstruction/ppmc_cuda_tests.cu @@ -24,16 +24,14 @@ TEST(tHYDROPpmcCTUReconstructor, CorrectInputExpectCorrectOutput) { - // Alwin: skip until this has been fixed - GTEST_SKIP(); // Set up PRNG to use std::mt19937_64 prng(42); std::uniform_real_distribution doubleRand(0.1, 5); // Mock up needed information - size_t const nx = 7; - size_t const ny = 7; - size_t const nz = 7; + size_t const nx = 6; + size_t const ny = 6; + size_t const nz = 6; size_t const n_fields = 5; double const dx = doubleRand(prng); double const dt = doubleRand(prng); @@ -50,37 +48,37 @@ TEST(tHYDROPpmcCTUReconstructor, CorrectInputExpectCorrectOutput) dev_grid.cpyHostToDevice(host_grid); // Fiducial Data - std::vector> fiducial_interface_left = {{{171, 1.7598055553475744}, - {514, 3.3921082637175894}, - {857, 3.5866056366266772}, - {1200, 3.4794572581328902}, - {1543, 10.363861270296034}}, - {{171, 1.6206985712721598}, - {514, 3.123972986618837}, - {857, 3.30309596610488}, - {1200, 3.204417323222251}, - {1543, 9.544631281899882}}, - {{171, 1.6206985712721595}, - {514, 5.0316428671215876}, - {857, 2.3915465711497186}, - {1200, 3.2044173232222506}, - {1543, 12.74302824034023}}}; + std::vector> fiducial_interface_left = {{{86, 2.6558981128823214}, + {302, 0.84399195916314151}, + {518, 2.2002498722761787}, + {734, 1.764334292986655}, + {950, 3.3600925565746804}}, + {{86, 2.4950488327292639}, + {302, 0.79287723513518138}, + {518, 1.7614576990062414}, + {734, 1.8238574169157304}, + {950, 3.14294317122161}}, + {{86, 2.6558981128823214}, + {302, 0.84399195916314151}, + {518, 2.0109603398129137}, + {734, 1.764334292986655}, + {950, 3.2100231679403066}}}; - std::vector> fiducial_interface_right = {{{170, 1.7857012385420896}, - {513, 3.4420234152477129}, - {856, 3.6393828329638049}, - {1199, 3.5306577572855762}, - {1542, 10.516366339570284}}, - {{164, 1.6206985712721595}, - {507, 3.1239729866188366}, - {850, 3.3030959661048795}, - {1193, 3.2044173232222506}, - {1536, 9.5446312818998802}}, - {{122, 1.6206985712721595}, - {465, 5.4375307473677061}, - {808, 2.2442413290889327}, - {1151, 3.2044173232222506}, - {1494, 13.843305272338561}}}; + std::vector> fiducial_interface_right = {{{85, 2.6558981128823214}, + {301, 0.84399195916314151}, + {517, 1.8381070277226794}, + {733, 1.764334292986655}, + {949, 3.0847691079841209}}, + {{80, 3.1281603739188069}, + {296, 0.99406757727427164}, + {512, 1.8732124042412865}, + {728, 1.6489758692176784}, + {944, 2.8820015278590443}}, + {{50, 2.6558981128823214}, + {266, 0.84399195916314151}, + {482, 2.0109603398129137}, + {698, 1.764334292986655}, + {914, 3.2100231679403066}}}; // Loop over different directions for (size_t direction = 0; direction < 3; direction++) { @@ -134,9 +132,9 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) std::uniform_real_distribution doubleRand(0.1, 5); // Mock up needed information - size_t const nx = 7; - size_t const ny = 7; - size_t const nz = 7; + size_t const nx = 6; + size_t const ny = 6; + size_t const nz = 6; double const gamma = 5.0 / 3.0; #ifdef MHD size_t const n_fields = 8; @@ -156,70 +154,70 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) // Fiducial Data #ifdef MHD - std::vector> fiducial_interface_left = {{{171, 1.5556846217288991}, - {514, 1.7422005905354798}, - {857, 3.6289199464135558}, - {1200, 2.1487031353407438}, - {1543, 22.988345461909127}, - {1886, 3.1027541330860546}, - {2229, 3.2554981416903335}}, - {{171, 1.7167767631895592}, - {514, 1.8447385381907686}, - {857, 2.9211469103910663}, - {1200, 2.626030390823102}, - {1543, 28.84165870179233}, - {1886, 3.8209152940021962}, - {2229, 2.7248523895714203}}, - {{171, 1.421933695280897}, - {514, 1.2318388818745061}, - {857, 2.8667822907691818}, - {1200, 2.1256773710028964}, - {1543, 15.684026541123352}, - {1886, 2.3642698195433232}, - {2229, 2.9207483994866617}}}; + std::vector> fiducial_interface_left = {{{86, 3.6926886385390683}, + {302, 2.3022467009220993}, + {518, 2.3207781368125389}, + {734, 2.6544338753333747}, + {950, 11.430630157120799}, + {1166, 0.6428577630032507}, + {1382, 4.1406925096276597}}, + {{86, 3.811691682348938}, + {302, 1.4827993897794758}, + {518, 2.3955690789476871}, + {734, 4.06241130448349}, + {950, 10.552876853630949}, + {1166, 3.5147238706385471}, + {1382, 1.2344879085821312}}, + {{86, 3.1608655959160155}, + {302, 1.5377824007725194}, + {518, 0.41798730655927896}, + {734, 2.2721408530383784}, + {950, 5.6329522765789646}, + {1166, 0.84450832590555991}, + {1382, 1.4279317910797107}}}; - std::vector> fiducial_interface_right = {{{170, 1.4838721492695441}, - {513, 1.3797509020377114}, - {856, 3.223172223924883}, - {1199, 2.2593969253004111}, - {1542, 15.634488002075017}, - {1885, 2.7494588681249819}, - {2228, 3.2540533219925698}}, - {{164, 1.4075989434297753}, - {507, 1.34947711631431}, - {850, 3.605198021293794}, - {1193, 1.9244827470895529}, - {1536, 13.52285212927548}, - {1879, 2.9568307038177966}, - {2222, 2.1086380065800636}}, - {{122, 1.9532382085816002}, - {465, 2.6860067041011249}, - {808, 5.1657781029381917}, - {1151, 2.7811084475444732}, - {1494, 24.999993264381686}, - {1837, 2.3090650532529238}, - {2180, 2.8525500781893642}}}; + std::vector> fiducial_interface_right = {{{85, 2.8949509658187838}, + {301, 0.25766140043685887}, + {517, 1.8194165731976308}, + {733, 2.0809921071868756}, + {949, 8.1315538869542046}, + {1165, 0.49708185787322312}, + {1381, 3.2017395511439881}}, + {{80, 2.8600082827930269}, + {296, 0.37343415089084014}, + {512, 1.7974558224423689}, + {728, 0.94369445956099784}, + {944, 7.7011501503138504}, + {1160, 3.5147238706385471}, + {1376, 1.2344879085821312}}, + {{50, 3.1608655959160155}, + {266, 0.32035830490636008}, + {482, 3.1721881746709815}, + {698, 2.2721408530383784}, + {914, 14.017699282483312}, + {1130, 1.5292690020097823}, + {1346, -0.12121484974901264}}}; #else // not MHD std::vector> fiducial_interface_left = { - {{171, 1.5239648818969727}, {514, 1.658831367400063}, {857, 3.3918153400617137}, {1200, 2.4096936604224304}}, - {{171, 1.5239639282226562}, {514, 1.6246850138898132}, {857, 3.391813217514656}, {1200, 2.3220060950058032}}, - {{171, 1.7062816619873047}, {514, 1.3300289077249516}, {857, 3.5599794228554593}, {1200, 2.5175993972231074}}}; + {{86, 4.155160222900312}, {302, 1.1624633361407897}, {518, 1.6379195998743412}, {734, 2.9868746414179093}}, + {{86, 4.1795874335665655}, {302, 2.1094239978455054}, {518, 2.6811988240843849}, {734, 4.2540957888954054}}, + {{86, 2.1772852940944429}, {302, 0.58167501916840214}, {518, 1.3683785996473696}, {734, 0.40276763592716164}}}; - std::vector> fiducial_interface_right = {{{135, 6.5824208227997447}, - {170, 1.5239620208740234}, - {513, 1.5386557138925041}, - {856, 3.3918089724205411}, - {1199, 1.9263881802230425}}, - {{135, 6.4095055796015963}, - {164, 1.5239639282226562}, - {507, 1.5544994569400168}, - {850, 3.391813217514656}, - {1193, 2.1017627061702138}}, - {{122, 1.3893871307373047}, - {135, 6.0894802934332555}, - {465, 2.1518846449159135}, - {808, 3.4792525252435533}, - {1151, 2.0500250813102903}}}; + std::vector> fiducial_interface_right = {{{54, 3.8655260187947502}, + {85, 2.6637168309565289}, + {301, 0.69483650107094164}, + {517, 2.7558388224532218}, + {733, 1.9147729154830744}}, + {{54, 5.7556871317935459}, + {80, 2.6515032256234021}, + {296, 0.39344537106429511}, + {512, 1.6491544916805785}, + {728, 0.85830485311660487}}, + {{50, 2.8254070932730269}, + {54, 2.1884721760267873}, + {266, 0.75482470285166003}, + {482, 1.7757096932649317}, + {698, 3.6101832818706452}}}; #endif // MHD // Loop over different directions diff --git a/src/reconstruction/reconstruction.h b/src/reconstruction/reconstruction.h index 3ed89a8b6..17a2f013e 100644 --- a/src/reconstruction/reconstruction.h +++ b/src/reconstruction/reconstruction.h @@ -78,6 +78,29 @@ struct Characteristic { }; // ===================================================================================================================== +// ===================================================================================================================== +template +bool __device__ __host__ __inline__ Thread_Guard(int const &nx, int const &ny, int const &nz, int const &xid, + int const &yid, int const &zid) +{ + // x check + bool out_of_bounds_thread = xid < order - 1 or xid >= nx - order; + + // y check + if (ny > 1) { + out_of_bounds_thread = yid < order - 1 or yid >= ny - order or out_of_bounds_thread; + } + + // z check + if (nz > 1) { + out_of_bounds_thread = zid < order - 1 or zid >= nz - order or out_of_bounds_thread; + } + out_of_bounds_thread = zid >= nz or out_of_bounds_thread; + + return out_of_bounds_thread; +} +// ===================================================================================================================== + // ===================================================================================================================== /*! * \brief Load the data for reconstruction From 69d43a68228f60e551011f377e719058f4de8840 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 15:08:12 -0400 Subject: [PATCH 13/14] Add test for reconstruction::Thread_Guard --- src/reconstruction/reconstruction_tests.cu | 29 ++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/reconstruction/reconstruction_tests.cu b/src/reconstruction/reconstruction_tests.cu index 62e615b39..5f8000bf8 100644 --- a/src/reconstruction/reconstruction_tests.cu +++ b/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 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(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 From 890310ddfe75e929c2e777c544098bb2c1ffb901 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 8 Aug 2023 15:59:53 -0400 Subject: [PATCH 14/14] Add documentation for reconstruction::Thread_Guard --- src/reconstruction/reconstruction.h | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/reconstruction/reconstruction.h b/src/reconstruction/reconstruction.h index 17a2f013e..07aae21a6 100644 --- a/src/reconstruction/reconstruction.h +++ b/src/reconstruction/reconstruction.h @@ -79,23 +79,41 @@ 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 bool __device__ __host__ __inline__ Thread_Guard(int const &nx, int const &ny, int const &nz, int const &xid, int const &yid, int const &zid) { - // x check + // 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 + // 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 + // 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; } - out_of_bounds_thread = zid >= nz 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; }