diff --git a/src/integrators/VL_1D_cuda.cu b/src/integrators/VL_1D_cuda.cu index f2ad520b8..48ead68f2 100644 --- a/src/integrators/VL_1D_cuda.cu +++ b/src/integrators/VL_1D_cuda.cu @@ -92,7 +92,7 @@ void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Rea n_fields); #endif #ifdef PLMC - hipLaunchKernelGGL(PLMC_cuda, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0, + hipLaunchKernelGGL(PLMC_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, n_fields); #endif #ifdef PLMP @@ -104,7 +104,7 @@ void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Rea gama, 0, n_fields); #endif #ifdef PPMC - hipLaunchKernelGGL(PPMC_VL, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama, 0); + hipLaunchKernelGGL(PPMC_VL<0>, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama); #endif GPU_Error_Check(); diff --git a/src/integrators/VL_2D_cuda.cu b/src/integrators/VL_2D_cuda.cu index 3c8992d71..f5f305b28 100644 --- a/src/integrators/VL_2D_cuda.cu +++ b/src/integrators/VL_2D_cuda.cu @@ -102,10 +102,10 @@ void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_of dt, gama, 1, n_fields); #endif #ifdef PLMC - hipLaunchKernelGGL(PLMC_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, - 0, n_fields); - hipLaunchKernelGGL(PLMC_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, - 1, n_fields); + hipLaunchKernelGGL(PLMC_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, + gama, n_fields); + hipLaunchKernelGGL(PLMC_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, + gama, n_fields); #endif #ifdef PPMP hipLaunchKernelGGL(PPMP_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, @@ -114,8 +114,8 @@ void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_of dt, gama, 1, n_fields); #endif // PPMP #ifdef PPMC - hipLaunchKernelGGL(PPMC_VL, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama, 0); - hipLaunchKernelGGL(PPMC_VL, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, gama, 1); + hipLaunchKernelGGL(PPMC_VL<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama); + hipLaunchKernelGGL(PPMC_VL<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, gama); #endif // PPMC GPU_Error_Check(); diff --git a/src/integrators/VL_3D_cuda.cu b/src/integrators/VL_3D_cuda.cu index 1d30cc0e2..5dde15443 100644 --- a/src/integrators/VL_3D_cuda.cu +++ b/src/integrators/VL_3D_cuda.cu @@ -145,8 +145,7 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int // Step 2: Calculate first-order upwind fluxes #ifdef EXACT - cuda_utilities::AutomaticLaunchParams static const exact_launch_params(Calculate_Exact_Fluxes_CUDA, - n_cellsCalculate_Exact_Fluxes_CUDA); + cuda_utilities::AutomaticLaunchParams static const exact_launch_params(Calculate_Exact_Fluxes_CUDA, n_cells); hipLaunchKernelGGL(Calculate_Exact_Fluxes_CUDA, exact_launch_params.get_numBlocks(), exact_launch_params.get_threadsPerBlock(), 0, 0, Q_Lx, Q_Rx, F_x, nx, ny, nz, n_ghost, gama, 0, n_fields); @@ -250,13 +249,13 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, n_ghost, dz, dt, gama, 2, n_fields); #endif // PLMP #ifdef PLMC - cuda_utilities::AutomaticLaunchParams static const plmc_vl_launch_params(PLMC_cuda, n_cells); - hipLaunchKernelGGL(PLMC_cuda, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), 0, - 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0, n_fields); - hipLaunchKernelGGL(PLMC_cuda, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), 0, - 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, 1, n_fields); - hipLaunchKernelGGL(PLMC_cuda, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), 0, - 0, dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama, 2, n_fields); + cuda_utilities::AutomaticLaunchParams static const plmc_vl_launch_params(PLMC_cuda<0>, n_cells); + hipLaunchKernelGGL(PLMC_cuda<0>, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), + 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, n_fields); + hipLaunchKernelGGL(PLMC_cuda<1>, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), + 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, n_fields); + hipLaunchKernelGGL(PLMC_cuda<2>, plmc_vl_launch_params.get_numBlocks(), plmc_vl_launch_params.get_threadsPerBlock(), + 0, 0, dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama, n_fields); #endif // PLMC #ifdef PPMP cuda_utilities::AutomaticLaunchParams static const ppmp_launch_params(PPMP_cuda, n_cells); @@ -268,13 +267,13 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, n_ghost, dz, dt, gama, 2, n_fields); #endif // PPMP #ifdef PPMC - cuda_utilities::AutomaticLaunchParams static const ppmc_vl_launch_params(PPMC_VL, n_cells); - hipLaunchKernelGGL(PPMC_VL, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, 0, - dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama, 0); - hipLaunchKernelGGL(PPMC_VL, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, 0, - dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, gama, 1); - hipLaunchKernelGGL(PPMC_VL, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, 0, - dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, gama, 2); + cuda_utilities::AutomaticLaunchParams static const ppmc_vl_launch_params(PPMC_VL<0>, n_cells); + hipLaunchKernelGGL(PPMC_VL<0>, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, + 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, gama); + hipLaunchKernelGGL(PPMC_VL<1>, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, + 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, gama); + hipLaunchKernelGGL(PPMC_VL<2>, ppmc_vl_launch_params.get_numBlocks(), ppmc_vl_launch_params.get_threadsPerBlock(), 0, + 0, dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, gama); #endif // PPMC GPU_Error_Check(); diff --git a/src/integrators/simple_1D_cuda.cu b/src/integrators/simple_1D_cuda.cu index 80f26021a..e8a8af067 100644 --- a/src/integrators/simple_1D_cuda.cu +++ b/src/integrators/simple_1D_cuda.cu @@ -64,7 +64,7 @@ void Simple_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, GPU_Error_Check(); #endif #ifdef PLMC - hipLaunchKernelGGL(PLMC_cuda, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0, + hipLaunchKernelGGL(PLMC_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, n_fields); GPU_Error_Check(); #endif @@ -74,7 +74,7 @@ void Simple_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, GPU_Error_Check(); #endif #ifdef PPMC - hipLaunchKernelGGL(PPMC_CTU, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0); + hipLaunchKernelGGL(PPMC_CTU<0>, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); GPU_Error_Check(); #endif diff --git a/src/integrators/simple_2D_cuda.cu b/src/integrators/simple_2D_cuda.cu index 97d435c51..33b54e9b7 100644 --- a/src/integrators/simple_2D_cuda.cu +++ b/src/integrators/simple_2D_cuda.cu @@ -65,9 +65,9 @@ void Simple_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int gama, 1, n_fields); #endif #ifdef PLMC - hipLaunchKernelGGL(PLMC_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0, + hipLaunchKernelGGL(PLMC_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, n_fields); - hipLaunchKernelGGL(PLMC_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, 1, + hipLaunchKernelGGL(PLMC_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, n_fields); #endif #ifdef PPMP @@ -77,8 +77,8 @@ void Simple_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int gama, 1, n_fields); #endif #ifdef PPMC - hipLaunchKernelGGL(PPMC_CTU, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0); - hipLaunchKernelGGL(PPMC_CTU, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, 1); + hipLaunchKernelGGL(PPMC_CTU<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + hipLaunchKernelGGL(PPMC_CTU<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); #endif GPU_Error_Check(); diff --git a/src/integrators/simple_3D_cuda.cu b/src/integrators/simple_3D_cuda.cu index 528eab04f..c2ce8e1e1 100644 --- a/src/integrators/simple_3D_cuda.cu +++ b/src/integrators/simple_3D_cuda.cu @@ -99,11 +99,11 @@ void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, gama, 2, n_fields); #endif // PLMP #ifdef PLMC - hipLaunchKernelGGL(PLMC_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0, + hipLaunchKernelGGL(PLMC_cuda<0>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, n_fields); - hipLaunchKernelGGL(PLMC_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, 1, + hipLaunchKernelGGL(PLMC_cuda<1>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, n_fields); - hipLaunchKernelGGL(PLMC_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama, 2, + hipLaunchKernelGGL(PLMC_cuda<2>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama, n_fields); #endif #ifdef PPMP @@ -115,9 +115,9 @@ void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, gama, 2, n_fields); #endif // PPMP #ifdef PPMC - hipLaunchKernelGGL(PPMC_CTU, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama, 0); - hipLaunchKernelGGL(PPMC_CTU, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama, 1); - hipLaunchKernelGGL(PPMC_CTU, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama, 2); + hipLaunchKernelGGL(PPMC_CTU<0>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + hipLaunchKernelGGL(PPMC_CTU<1>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); + hipLaunchKernelGGL(PPMC_CTU<2>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama); GPU_Error_Check(); #endif // PPMC diff --git a/src/reconstruction/plmc_cuda.cu b/src/reconstruction/plmc_cuda.cu index c670723dd..83004d0c7 100644 --- a/src/reconstruction/plmc_cuda.cu +++ b/src/reconstruction/plmc_cuda.cu @@ -21,8 +21,9 @@ gamma, int dir) * \brief When passed a stencil of conserved variables, returns the left and right boundary values for the interface calculated using plm. */ +template __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real dx, Real dt, Real gamma, int dir, int n_fields) + int ny, int nz, Real dx, Real dt, Real gamma, int n_fields) { // get a thread ID int const thread_id = threadIdx.x + blockIdx.x * blockDim.x; @@ -60,15 +61,15 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_ // load the 3-cell stencil into registers // cell i hydro_utilities::Primitive const cell_i = - reconstruction::Load_Data(dev_conserved, xid, yid, zid, nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Load_Cell_Primitive(dev_conserved, xid, yid, zid, nx, ny, n_cells, gamma); // cell i-1. The equality checks the direction and will subtract one from the correct direction - hydro_utilities::Primitive const cell_imo = reconstruction::Load_Data( - dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_imo = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, gamma); // cell i+1. The equality checks the direction and add one to the correct direction - hydro_utilities::Primitive const cell_ipo = reconstruction::Load_Data( - dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_ipo = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, gamma); // calculate the adiabatic sound speed in cell i Real const sound_speed = hydro_utilities::Calc_Sound_Speed(cell_i.pressure, cell_i.density, gamma); @@ -295,3 +296,14 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_ id = cuda_utilities::compute1DIndex(xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny); reconstruction::Write_Data(interface_R_imh, dev_bounds_R, dev_conserved, id, n_cells, o1, o2, o3, gamma); } + +// Instantiate the relevant template specifications +template __global__ __launch_bounds__(TPB) void PLMC_cuda<0>(Real *dev_conserved, Real *dev_bounds_L, + Real *dev_bounds_R, int nx, int ny, int nz, Real dx, + Real dt, Real gamma, int n_fields); +template __global__ __launch_bounds__(TPB) void PLMC_cuda<1>(Real *dev_conserved, Real *dev_bounds_L, + Real *dev_bounds_R, int nx, int ny, int nz, Real dx, + Real dt, Real gamma, int n_fields); +template __global__ __launch_bounds__(TPB) void PLMC_cuda<2>(Real *dev_conserved, Real *dev_bounds_L, + Real *dev_bounds_R, int nx, int ny, int nz, Real dx, + Real dt, Real gamma, int n_fields); diff --git a/src/reconstruction/plmc_cuda.h b/src/reconstruction/plmc_cuda.h index c2d25df84..d50f493a1 100644 --- a/src/reconstruction/plmc_cuda.h +++ b/src/reconstruction/plmc_cuda.h @@ -15,7 +15,8 @@ gamma, int dir) * \brief When passed a stencil of conserved variables, returns the left and right boundary values for the interface calculated using plm. */ +template __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real dx, Real dt, Real gamma, int dir, int n_fields); + int ny, int nz, Real dx, Real dt, Real gamma, int n_fields); #endif // PLMC_CUDA_H diff --git a/src/reconstruction/plmc_cuda_tests.cu b/src/reconstruction/plmc_cuda_tests.cu index dd0f35891..6b8e95021 100644 --- a/src/reconstruction/plmc_cuda_tests.cu +++ b/src/reconstruction/plmc_cuda_tests.cu @@ -148,8 +148,21 @@ TEST(tHYDROPlmcReconstructor, CorrectInputExpectCorrectOutput) 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); + std::cout << "direction = " << direction << std::endl; + switch (direction) { + case 0: + hipLaunchKernelGGL(PLMC_cuda<0>, 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, n_fields); + break; + case 1: + hipLaunchKernelGGL(PLMC_cuda<1>, 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, n_fields); + break; + case 2: + hipLaunchKernelGGL(PLMC_cuda<2>, 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, n_fields); + break; + } GPU_Error_Check(); GPU_Error_Check(cudaDeviceSynchronize()); @@ -261,8 +274,20 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput) 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); + switch (direction) { + case 0: + hipLaunchKernelGGL(PLMC_cuda<0>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, n_fields); + break; + case 1: + hipLaunchKernelGGL(PLMC_cuda<1>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, n_fields); + break; + case 2: + hipLaunchKernelGGL(PLMC_cuda<2>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, n_fields); + break; + } GPU_Error_Check(); GPU_Error_Check(cudaDeviceSynchronize()); diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index c98435075..b610b9902 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -19,8 +19,9 @@ /*! * \brief When passed a stencil of conserved variables, returns the left and right boundary values for the interface calculated using ppm. */ +template __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, int nz, Real dx, - Real dt, Real gamma, int dir) + Real dt, Real gamma) { // get a thread ID int const thread_id = threadIdx.x + blockIdx.x * blockDim.x; @@ -57,29 +58,27 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun // load the 5-cell stencil into registers // cell i hydro_utilities::Primitive const cell_i = - reconstruction::Load_Data(dev_conserved, xid, yid, zid, nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Load_Cell_Primitive(dev_conserved, xid, yid, zid, nx, ny, n_cells, gamma); // cell i-1. The equality checks check the direction and subtracts one from the direction // im1 stands for "i minus 1" - hydro_utilities::Primitive const cell_im1 = reconstruction::Load_Data( - dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_im1 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, gamma); // cell i+1. The equality checks check the direction and adds one to the direction // ip1 stands for "i plus 1" - hydro_utilities::Primitive const cell_ip1 = reconstruction::Load_Data( - dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_ip1 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, gamma); // cell i-2. The equality checks check the direction and subtracts one from the direction // im2 stands for "i minus 2" - hydro_utilities::Primitive const cell_im2 = - reconstruction::Load_Data(dev_conserved, xid - 2 * int(dir == 0), yid - 2 * int(dir == 1), - zid - 2 * int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_im2 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid - 2 * int(dir == 0), yid - 2 * int(dir == 1), zid - 2 * int(dir == 2), nx, ny, n_cells, gamma); // cell i+2. The equality checks check the direction and adds one to the direction // ip2 stands for "i plus 2" - hydro_utilities::Primitive const cell_ip2 = - reconstruction::Load_Data(dev_conserved, xid + 2 * int(dir == 0), yid + 2 * int(dir == 1), - zid + 2 * int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_ip2 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid + 2 * int(dir == 0), yid + 2 * int(dir == 1), zid + 2 * int(dir == 2), nx, ny, n_cells, gamma); // Steps 2 - 5 are repeated for cell i-1, i, and i+1 @@ -541,8 +540,9 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun // ===================================================================================================================== // ===================================================================================================================== +template __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real gamma, int dir) + int ny, int nz, Real gamma) { // get a thread ID int const thread_id = threadIdx.x + blockIdx.x * blockDim.x; @@ -580,29 +580,27 @@ __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bo // load the 5-cell stencil into registers // cell i hydro_utilities::Primitive const cell_i = - reconstruction::Load_Data(dev_conserved, xid, yid, zid, nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Load_Cell_Primitive(dev_conserved, xid, yid, zid, nx, ny, n_cells, gamma); // cell i-1. The equality checks the direction and will subtract one from the correct direction // im1 stands for "i minus 1" - hydro_utilities::Primitive const cell_im1 = reconstruction::Load_Data( - dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_im1 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny, n_cells, gamma); // cell i+1. The equality checks the direction and add one to the correct direction // ip1 stands for "i plus 1" - hydro_utilities::Primitive const cell_ip1 = reconstruction::Load_Data( - dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_ip1 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid + int(dir == 0), yid + int(dir == 1), zid + int(dir == 2), nx, ny, n_cells, gamma); // cell i-2. The equality checks the direction and will subtract two from the correct direction // im2 stands for "i minus 2" - hydro_utilities::Primitive const cell_im2 = - reconstruction::Load_Data(dev_conserved, xid - 2 * int(dir == 0), yid - 2 * int(dir == 1), - zid - 2 * int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_im2 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid - 2 * int(dir == 0), yid - 2 * int(dir == 1), zid - 2 * int(dir == 2), nx, ny, n_cells, gamma); // cell i+2. The equality checks the direction and add two to the correct direction // ip2 stands for "i plus 2" - hydro_utilities::Primitive const cell_ip2 = - reconstruction::Load_Data(dev_conserved, xid + 2 * int(dir == 0), yid + 2 * int(dir == 1), - zid + 2 * int(dir == 2), nx, ny, n_cells, o1, o2, o3, gamma); + hydro_utilities::Primitive const cell_ip2 = hydro_utilities::Load_Cell_Primitive( + dev_conserved, xid + 2 * int(dir == 0), yid + 2 * int(dir == 1), zid + 2 * int(dir == 2), nx, ny, n_cells, gamma); // Convert to the characteristic variables Real const sound_speed = hydro_utilities::Calc_Sound_Speed(cell_i.pressure, cell_i.density, gamma); @@ -701,4 +699,18 @@ __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bo id = cuda_utilities::compute1DIndex(xid - int(dir == 0), yid - int(dir == 1), zid - int(dir == 2), nx, ny); reconstruction::Write_Data(interface_R_imh, dev_bounds_R, dev_conserved, id, n_cells, o1, o2, o3, gamma); } +// Instantiate the relevant template specifications +template __global__ void PPMC_CTU<0>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, + int nz, Real dx, Real dt, Real gamma); +template __global__ void PPMC_CTU<1>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, + int nz, Real dx, Real dt, Real gamma); +template __global__ void PPMC_CTU<2>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, + int nz, Real dx, Real dt, Real gamma); + +template __global__ __launch_bounds__(TPB) void PPMC_VL<0>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, + int nx, int ny, int nz, Real gamma); +template __global__ __launch_bounds__(TPB) void PPMC_VL<1>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, + int nx, int ny, int nz, Real gamma); +template __global__ __launch_bounds__(TPB) void PPMC_VL<2>(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, + int nx, int ny, int nz, Real gamma); // ===================================================================================================================== diff --git a/src/reconstruction/ppmc_cuda.h b/src/reconstruction/ppmc_cuda.h index 916853874..a820b7e0b 100644 --- a/src/reconstruction/ppmc_cuda.h +++ b/src/reconstruction/ppmc_cuda.h @@ -14,6 +14,7 @@ * in the characteristic variables to monotonize the slopes followed by limiting the interface states using the limiter * from Colella & Woodward 1984. * + * \tparam dir The direction to reconstruct. 0=X, 1=Y, 2=Z * \param[in] dev_conserved The conserved variable array * \param[out] dev_bounds_L The array of left interfaces * \param[out] dev_bounds_R The array of right interfaces @@ -23,10 +24,10 @@ * \param[in] dx The length of the cells in the `dir` direction * \param[in] dt The time step * \param[in] gamma The adiabatic index - * \param[in] dir The direction to reconstruct. 0=X, 1=Y, 2=Z */ +template __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, int nz, Real dx, - Real dt, Real gamma, int dir); + Real dt, Real gamma); /*! * \brief Computes the left and right interface states using PPM with limiting in the characteristic variables. Used for @@ -38,6 +39,7 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun * oscillatory, and faster than the method described in Stone et al. 2008 which is used in PPMC_CTU. The difference is * most pronounced in the Brio & Wu shock tube where the PPM oscillations are much smaller using this method. * + * \tparam dir The direction to reconstruct. 0=X, 1=Y, 2=Z * \param[in] dev_conserved The conserved variable array * \param[out] dev_bounds_L The array of left interfaces * \param[out] dev_bounds_R The array of right interfaces @@ -45,9 +47,9 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun * \param[in] ny The number of cells in the Y-direction * \param[in] nz The number of cells in the Z-direction * \param[in] gamma The adiabatic index - * \param[in] dir The direction to reconstruct. 0=X, 1=Y, 2=Z */ +template __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real gamma, int dir); + int ny, int nz, Real gamma); #endif // PPMC_CUDA_H diff --git a/src/reconstruction/ppmc_cuda_tests.cu b/src/reconstruction/ppmc_cuda_tests.cu index 9e9b11140..bcacb9869 100644 --- a/src/reconstruction/ppmc_cuda_tests.cu +++ b/src/reconstruction/ppmc_cuda_tests.cu @@ -87,8 +87,20 @@ TEST(tHYDROPpmcCTUReconstructor, CorrectInputExpectCorrectOutput) cuda_utilities::DeviceVector dev_interface_right(host_grid.size(), true); // Launch kernel - hipLaunchKernelGGL(PPMC_CTU, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), - dev_interface_right.data(), nx, ny, nz, dx, dt, gamma, direction); + switch (direction) { + case 0: + hipLaunchKernelGGL(PPMC_CTU<0>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma); + break; + case 1: + hipLaunchKernelGGL(PPMC_CTU<1>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma); + break; + case 2: + hipLaunchKernelGGL(PPMC_CTU<2>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, dx, dt, gamma); + break; + } GPU_Error_Check(); GPU_Error_Check(cudaDeviceSynchronize()); @@ -227,8 +239,20 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) cuda_utilities::DeviceVector dev_interface_right(nx * ny * nz * (n_fields - 1), true); // Launch kernel - hipLaunchKernelGGL(PPMC_VL, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), - dev_interface_right.data(), nx, ny, nz, gamma, direction); + switch (direction) { + case 0: + hipLaunchKernelGGL(PPMC_VL<0>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, gamma); + break; + case 1: + hipLaunchKernelGGL(PPMC_VL<1>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, gamma); + break; + case 2: + hipLaunchKernelGGL(PPMC_VL<2>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + dev_interface_right.data(), nx, ny, nz, gamma); + break; + } GPU_Error_Check(); GPU_Error_Check(cudaDeviceSynchronize()); diff --git a/src/reconstruction/reconstruction_internals.h b/src/reconstruction/reconstruction_internals.h index c60fff7fa..1d2eb7f63 100644 --- a/src/reconstruction/reconstruction_internals.h +++ b/src/reconstruction/reconstruction_internals.h @@ -48,7 +48,7 @@ enum Kind { chosen = plmc #elif defined(PPMP) chosen = ppmp -#elif defined(PLMC) +#elif defined(PPMC) chosen = ppmc #else #error "no reconstruction selected" diff --git a/src/utils/basic_structs.h b/src/utils/basic_structs.h index f51b733ae..ae3288e19 100644 --- a/src/utils/basic_structs.h +++ b/src/utils/basic_structs.h @@ -87,6 +87,22 @@ struct Conserved { #ifdef SCALAR Real scalar[grid_enum::nscalars]; #endif // SCALAR + + /// Default constructor, should init everything to zero + Conserved() = default; + /// Manual constructor, mostly used for testing and doesn't init all members + Conserved(Real const in_density, VectorXYZ const& in_momentum, Real const in_energy, + VectorXYZ const& in_magnetic = {0, 0, 0}, Real const in_gas_energy = 0.0) + : density(in_density), momentum(in_momentum), energy(in_energy) + { +#ifdef MHD + magnetic = in_magnetic; +#endif // mhd + +#ifdef DE + gas_energy = in_gas_energy; +#endif // DE + }; }; // ===================================================================================================================== diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index 5e61f45b8..57e91b3ec 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -211,4 +211,215 @@ inline __host__ __device__ Real Calc_Sound_Speed(Real const &P, Real const &d, R return sqrt(gamma * P / d); } +// ===================================================================================================================== +/*! + * \brief Load the conserved variables from a single cell. Note that with MHD this returns cell-centered magnetic + * fields, not face centered fields. + * + * \tparam dir The direction to load the data in. i.e. which direction is the cell x, y, and z. If dir = 0 then local + * xyz is the same as the true xyz. If dir = 1 then local xyz is true yzx, for dir = 2 then local xyz is true zxy. This + * option is so that the reconstructors, and any other split kernels, can load data that is in the appropriated + * direction for that kernel. + * \param[in] dev_conserved The pointer to the conserved variable array + * \param[in] xid The index for the x location + * \param[in] yid The index for the y location + * \param[in] zid The index for the z location + * \param[in] nx The total number of cells in the x direction + * \param[in] ny The total number of cells in the y direction + * \param[in] n_cells The total number of cells + * \return Conserved The cell centered conserved variables in the cell at location xid, yid, zid. + */ +template +inline __host__ __device__ Conserved Load_Cell_Conserved(Real const *dev_conserved, size_t const xid, size_t const yid, + size_t const zid, size_t const nx, size_t const ny, + size_t const n_cells) +{ + // First, check that our direction is correct + static_assert((0 <= dir) and (dir <= 2), "dir is not in the proper range"); + + // Compute index + size_t const cell_id = cuda_utilities::compute1DIndex(xid, yid, zid, nx, ny); + + // Load all the data + Conserved loaded_data; + + // Hydro variables + loaded_data.density = dev_conserved[cell_id + n_cells * grid_enum::density]; + loaded_data.momentum.x() = dev_conserved[cell_id + n_cells * grid_enum::momentum_x]; + loaded_data.momentum.y() = dev_conserved[cell_id + n_cells * grid_enum::momentum_y]; + loaded_data.momentum.z() = dev_conserved[cell_id + n_cells * grid_enum::momentum_z]; + loaded_data.energy = dev_conserved[cell_id + n_cells * grid_enum::Energy]; + +#ifdef MHD + // These are all cell centered values + loaded_data.magnetic = mhd::utils::cellCenteredMagneticFields(dev_conserved, cell_id, xid, yid, zid, n_cells, nx, ny); +#endif // MHD + +#ifdef DE + loaded_data.gas_energy = dev_conserved[cell_id + n_cells * grid_enum::GasEnergy]; +#endif // DE + +#ifdef SCALAR + for (size_t i = 0; i < grid_enum::nscalars; i++) { + loaded_data.scalar[i] = dev_conserved[cell_id + n_cells * (grid_enum::scalar + i)]; + } +#endif // SCALAR + + // Now that all the data is loaded, let's sort out the direction + // if constexpr(dir == 0) in this case everything is already set so we'll skip this case + if constexpr (dir == 1) { + math_utils::Cyclic_Permute_Once(loaded_data.momentum); +#ifdef MHD + math_utils::Cyclic_Permute_Once(loaded_data.magnetic); +#endif // MHD + } else if constexpr (dir == 2) { + math_utils::Cyclic_Permute_Twice(loaded_data.momentum); +#ifdef MHD + math_utils::Cyclic_Permute_Twice(loaded_data.magnetic); +#endif // MHD + } + + return loaded_data; +} +// ===================================================================================================================== + +// ===================================================================================================================== +/*! + * \brief Convert Conserved cell centered variables to primitive variables + * + * \param conserved_in The conserved variables to convert + * \param gamma The adiabatic index + * \return Primitive The cell centered primitive variables + */ +__inline__ __host__ __device__ Primitive Conserved_2_Primitive(Conserved const &conserved_in, Real const gamma) +{ + Primitive output; + + // First the easy ones + output.density = conserved_in.density; + output.velocity.x() = conserved_in.momentum.x() / conserved_in.density; + output.velocity.y() = conserved_in.momentum.y() / conserved_in.density; + output.velocity.z() = conserved_in.momentum.z() / conserved_in.density; + +#ifdef MHD + output.magnetic.x() = conserved_in.magnetic.x(); + output.magnetic.y() = conserved_in.magnetic.y(); + output.magnetic.z() = conserved_in.magnetic.z(); +#endif // MHD + +#ifdef DE + output.gas_energy_specific = conserved_in.gas_energy / conserved_in.density; +#endif // DE + +#ifdef SCALAR + for (size_t i = 0; i < grid_enum::nscalars; i++) { + output.scalar_specific[i] = conserved_in.scalar[i] / conserved_in.density; + } +#endif // SCALAR + +// Now that the easy ones are done let's figure out the pressure +#ifdef DE // DE + Real E_non_thermal = hydro_utilities::Calc_Kinetic_Energy_From_Velocity(output.density, output.velocity.x(), + output.velocity.y(), output.velocity.z()); + + #ifdef MHD + E_non_thermal += mhd::utils::computeMagneticEnergy(output.magnetic.x(), output.magnetic.y(), output.magnetic.z()); + #endif // MHD + + output.pressure = hydro_utilities::Get_Pressure_From_DE(conserved_in.energy, conserved_in.energy - E_non_thermal, + conserved_in.gas_energy, gamma); +#else // not DE + #ifdef MHD + output.pressure = hydro_utilities::Calc_Pressure_Primitive( + conserved_in.energy, conserved_in.density, output.velocity.x(), output.velocity.y(), output.velocity.z(), gamma, + output.magnetic.x(), output.magnetic.y(), output.magnetic.z()); + #else // not MHD + output.pressure = hydro_utilities::Calc_Pressure_Primitive( + conserved_in.energy, conserved_in.density, output.velocity.x(), output.velocity.y(), output.velocity.z(), gamma); + #endif // MHD +#endif // DE + + return output; +} +// ===================================================================================================================== + +// ===================================================================================================================== +/*! + * \brief Convert primitive cell centered variables to conserved variables + * + * \param primitive_in The primitive variables to convert + * \param gamma The adiabatic index + * \return Conserved The cell centered conserved variables + */ +__inline__ __host__ __device__ Conserved Primitive_2_Conserved(Primitive const &primitive_in, Real const gamma) +{ + Conserved output; + + // First the easy ones + output.density = primitive_in.density; + output.momentum.x() = primitive_in.velocity.x() * primitive_in.density; + output.momentum.y() = primitive_in.velocity.y() * primitive_in.density; + output.momentum.z() = primitive_in.velocity.z() * primitive_in.density; + +#ifdef MHD + output.magnetic.x() = primitive_in.magnetic.x(); + output.magnetic.y() = primitive_in.magnetic.y(); + output.magnetic.z() = primitive_in.magnetic.z(); +#endif // MHD + +#ifdef DE + output.gas_energy = primitive_in.gas_energy_specific * primitive_in.density; +#endif // DE + +#ifdef SCALAR + for (size_t i = 0; i < grid_enum::nscalars; i++) { + output.scalar[i] = primitive_in.scalar_specific[i] * primitive_in.density; + } +#endif // SCALAR + +// Now that the easy ones are done let's figure out the energy +#ifdef MHD + output.energy = hydro_utilities::Calc_Energy_Primitive(primitive_in.pressure, primitive_in.density, + primitive_in.velocity.x(), primitive_in.velocity.y(), + primitive_in.velocity.z(), gamma, primitive_in.magnetic.x(), + primitive_in.magnetic.y(), primitive_in.magnetic.z()); +#else // not MHD + output.energy = + hydro_utilities::Calc_Energy_Primitive(primitive_in.pressure, primitive_in.density, primitive_in.velocity.x(), + primitive_in.velocity.y(), primitive_in.velocity.z(), gamma); +#endif // MHD + + return output; +} +// ===================================================================================================================== + +// ===================================================================================================================== +/*! + * \brief Load the primitive variables from a single cell. Note that with MHD this returns cell-centered magnetic + * fields, not face centered fields. + * + * \tparam dir The direction to load the data in. i.e. which direction is the cell x, y, and z. If dir = 0 then local + * xyz is the same as the true xyz. If dir = 1 then local xyz is true yzx, for dir = 2 then local xyz is true zxy. This + * option is so that the reconstructors, and any other split kernels, can load data that is in the appropriated + * direction for that kernel. + * \param[in] dev_conserved The pointer to the conserved variable array + * \param[in] xid The index for the x location + * \param[in] yid The index for the y location + * \param[in] zid The index for the z location + * \param[in] nx The total number of cells in the x direction + * \param[in] ny The total number of cells in the y direction + * \param[in] n_cells The total number of cells + * \param[in] gamma The adiabatic index + * \return Primitive The cell centered conserved variables in the cell at location xid, yid, zid. + */ +template +inline __host__ __device__ Primitive Load_Cell_Primitive(Real const *dev_conserved, size_t const xid, size_t const yid, + size_t const zid, size_t const nx, size_t const ny, + size_t const n_cells, Real const gamma) +{ + Conserved const conserved_cell = Load_Cell_Conserved(dev_conserved, xid, yid, zid, nx, ny, n_cells); + return Conserved_2_Primitive(conserved_cell, gamma); +} +// ===================================================================================================================== + } // namespace hydro_utilities diff --git a/src/utils/hydro_utilities_tests.cpp b/src/utils/hydro_utilities_tests.cpp index fe0cbe9e6..b6467fbd3 100644 --- a/src/utils/hydro_utilities_tests.cpp +++ b/src/utils/hydro_utilities_tests.cpp @@ -261,4 +261,189 @@ TEST(tHYDROtMHDCalcKineticEnergyFromMomentum, CorrectInputExpectCorrectOutput) testing_utilities::Check_Results(fiducialEnergies.at(i), testEnergy, parameters.names.at(i)); } +} + +TEST(tALLLoadCellPrimitive, CorrectInputExpectCorrectOutput) +{ + // Set up test and mock up grid + size_t const nx = 3, ny = 3, nz = 3; + size_t const n_cells = nx * ny * nz; + size_t const xid = 1, yid = 1, zid = 1; + size_t const o1 = grid_enum::momentum_x, o2 = grid_enum::momentum_y, o3 = grid_enum::momentum_z; + Real const gamma = 5. / 3.; + + std::vector conserved(n_cells * grid_enum::num_fields); + std::iota(conserved.begin(), conserved.end(), 0.0); + + // Up the energy part of the grid to avoid negative pressure + for (size_t i = grid_enum::Energy * n_cells; i < (grid_enum::Energy + 1) * n_cells; i++) { + conserved.at(i) *= 5.0E2; + } + + for (int direction = 0; direction < 3; direction++) { + // Get test data + hydro_utilities::Primitive test_data; + + // Get the test data and permute the vector quantities back to the original order + switch (direction) { + case 0: + test_data = hydro_utilities::Load_Cell_Primitive<0>(conserved.data(), xid, yid, zid, nx, ny, n_cells, gamma); + break; + case 1: + test_data = hydro_utilities::Load_Cell_Primitive<1>(conserved.data(), xid, yid, zid, nx, ny, n_cells, gamma); + math_utils::Cyclic_Permute_Twice(test_data.velocity); +#ifdef MHD + math_utils::Cyclic_Permute_Twice(test_data.magnetic); +#endif // MHD + break; + case 2: + test_data = hydro_utilities::Load_Cell_Primitive<2>(conserved.data(), xid, yid, zid, nx, ny, n_cells, gamma); + math_utils::Cyclic_Permute_Once(test_data.velocity); +#ifdef MHD + math_utils::Cyclic_Permute_Once(test_data.magnetic); +#endif // MHD + break; + } + +// Check results +#ifdef MHD + hydro_utilities::Primitive const fiducial_data{ + 13, {3.0769230769230771, 5.1538461538461542, 7.2307692307692308}, 9662.3910256410272, {147.5, 173.5, 197.5}}; + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.velocity.x(), test_data.velocity.x(), "velocity.x()"); + testing_utilities::Check_Results(fiducial_data.velocity.y(), test_data.velocity.y(), "velocity.y()"); + testing_utilities::Check_Results(fiducial_data.velocity.z(), test_data.velocity.z(), "velocity.z()"); + testing_utilities::Check_Results(fiducial_data.pressure, test_data.pressure, "pressure"); + testing_utilities::Check_Results(fiducial_data.magnetic.x(), test_data.magnetic.x(), "magnetic.x()"); + testing_utilities::Check_Results(fiducial_data.magnetic.y(), test_data.magnetic.y(), "magnetic.y()"); + testing_utilities::Check_Results(fiducial_data.magnetic.z(), test_data.magnetic.z(), "magnetic.z()"); +#else // MHD + hydro_utilities::Primitive fiducial_data{ + 13, {3.0769230769230771, 5.1538461538461542, 7.2307692307692308}, 39950.641025641031}; + #ifdef DE + fiducial_data.pressure = 39950.641025641031; + #endif // DE + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.velocity.x(), test_data.velocity.x(), "velocity.x"); + testing_utilities::Check_Results(fiducial_data.velocity.y(), test_data.velocity.y(), "velocity.y"); + testing_utilities::Check_Results(fiducial_data.velocity.z(), test_data.velocity.z(), "velocity.z"); + testing_utilities::Check_Results(fiducial_data.pressure, test_data.pressure, "pressure"); +#endif // MHD + } +} + +TEST(tALLLoadCellConserved, CorrectInputExpectCorrectOutput) +{ + // Set up test and mock up grid + size_t const nx = 3, ny = 3, nz = 3; + size_t const n_cells = nx * ny * nz; + size_t const xid = 1, yid = 1, zid = 1; + size_t const o1 = grid_enum::momentum_x, o2 = grid_enum::momentum_y, o3 = grid_enum::momentum_z; + Real const gamma = 5. / 3.; + + std::vector conserved(n_cells * grid_enum::num_fields); + std::iota(conserved.begin(), conserved.end(), 0.0); + + // Up the energy part of the grid to avoid negative pressure + for (size_t i = grid_enum::Energy * n_cells; i < (grid_enum::Energy + 1) * n_cells; i++) { + conserved.at(i) *= 5.0E2; + } + + for (int direction = 0; direction < 3; direction++) { + // Get test data + hydro_utilities::Conserved test_data; + + // Get the test data and permute the vector quantities back to the original order + switch (direction) { + case 0: + test_data = hydro_utilities::Load_Cell_Conserved<0>(conserved.data(), xid, yid, zid, nx, ny, n_cells); + break; + case 1: + test_data = hydro_utilities::Load_Cell_Conserved<1>(conserved.data(), xid, yid, zid, nx, ny, n_cells); + math_utils::Cyclic_Permute_Twice(test_data.momentum); +#ifdef MHD + math_utils::Cyclic_Permute_Twice(test_data.magnetic); +#endif // MHD + break; + case 2: + test_data = hydro_utilities::Load_Cell_Conserved<2>(conserved.data(), xid, yid, zid, nx, ny, n_cells); + math_utils::Cyclic_Permute_Once(test_data.momentum); +#ifdef MHD + math_utils::Cyclic_Permute_Once(test_data.magnetic); +#endif // MHD + break; + } + + // Check results + hydro_utilities::Conserved fiducial_data{13, {40, 67, 94}, 60500, {147.5, 173.5, 197.5}}; +#ifdef MHD + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.momentum.x(), test_data.momentum.x(), "momentum.x"); + testing_utilities::Check_Results(fiducial_data.momentum.y(), test_data.momentum.y(), "momentum.y"); + testing_utilities::Check_Results(fiducial_data.momentum.z(), test_data.momentum.z(), "momentum.z"); + testing_utilities::Check_Results(fiducial_data.energy, test_data.energy, "energy"); + testing_utilities::Check_Results(fiducial_data.magnetic.x(), test_data.magnetic.x(), "magnetic.x"); + testing_utilities::Check_Results(fiducial_data.magnetic.y(), test_data.magnetic.y(), "magnetic.y"); + testing_utilities::Check_Results(fiducial_data.magnetic.z(), test_data.magnetic.z(), "magnetic.z"); +#else // MHD + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.momentum.x(), test_data.momentum.x(), "momentum.x"); + testing_utilities::Check_Results(fiducial_data.momentum.y(), test_data.momentum.y(), "momentum.y"); + testing_utilities::Check_Results(fiducial_data.momentum.z(), test_data.momentum.z(), "momentum.z"); + testing_utilities::Check_Results(fiducial_data.energy, test_data.energy, "energy"); +#endif // MHD + } +} + +TEST(tALLConserved2Primitive, CorrectInputExpectCorrectOutput) +{ + Real const gamma = 5. / 3.; + hydro_utilities::Conserved input_data{2, {2, 3, 4}, 90, {6, 7, 8}, 9}; + hydro_utilities::Primitive test_data = hydro_utilities::Conserved_2_Primitive(input_data, gamma); + + hydro_utilities::Primitive fiducial_data{2, {1, 1.5, 2}, 55.166666666666671, {6, 7, 8}, 4.5}; +#ifdef MHD + fiducial_data.pressure = 5.5000000000000009; +#endif // MHD + + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.velocity.x(), test_data.velocity.x(), "velocity.x"); + testing_utilities::Check_Results(fiducial_data.velocity.y(), test_data.velocity.y(), "velocity.y"); + testing_utilities::Check_Results(fiducial_data.velocity.z(), test_data.velocity.z(), "velocity.z"); + testing_utilities::Check_Results(fiducial_data.pressure, test_data.pressure, "pressure"); +#ifdef MHD + testing_utilities::Check_Results(fiducial_data.magnetic.x(), test_data.magnetic.x(), "magnetic.x"); + testing_utilities::Check_Results(fiducial_data.magnetic.y(), test_data.magnetic.y(), "magnetic.y"); + testing_utilities::Check_Results(fiducial_data.magnetic.z(), test_data.magnetic.z(), "magnetic.z"); +#endif // MHD +#ifdef DE + testing_utilities::Check_Results(fiducial_data.gas_energy_specific, test_data.gas_energy_specific, + "gas_energy_specific"); +#endif // DE +} + +TEST(tALLPrimitive2Conserved, CorrectInputExpectCorrectOutput) +{ + Real const gamma = 5. / 3.; + hydro_utilities::Primitive input_data{2, {2, 3, 4}, 90, {6, 7, 8}, 9}; + hydro_utilities::Conserved test_data = hydro_utilities::Primitive_2_Conserved(input_data, gamma); + + hydro_utilities::Conserved fiducial_data{2, {4, 6, 8}, 163.99999999999997, {6, 7, 8}, 18}; +#ifdef MHD + fiducial_data.energy = 238.49999999999997; +#endif // MHD + + testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density"); + testing_utilities::Check_Results(fiducial_data.momentum.x(), test_data.momentum.x(), "momentum.x"); + testing_utilities::Check_Results(fiducial_data.momentum.y(), test_data.momentum.y(), "momentum.y"); + testing_utilities::Check_Results(fiducial_data.momentum.z(), test_data.momentum.z(), "momentum.z"); + testing_utilities::Check_Results(fiducial_data.energy, test_data.energy, "energy"); +#ifdef MHD + testing_utilities::Check_Results(fiducial_data.magnetic.x(), test_data.magnetic.x(), "magnetic.x"); + testing_utilities::Check_Results(fiducial_data.magnetic.y(), test_data.magnetic.y(), "magnetic.y"); + testing_utilities::Check_Results(fiducial_data.magnetic.z(), test_data.magnetic.z(), "magnetic.z"); +#endif // MHD +#ifdef DE + testing_utilities::Check_Results(fiducial_data.gas_energy, test_data.gas_energy, "gas_energy"); +#endif // DE } \ No newline at end of file diff --git a/src/utils/math_utilities.h b/src/utils/math_utilities.h index 68d13f19d..b23fe5c9a 100644 --- a/src/utils/math_utilities.h +++ b/src/utils/math_utilities.h @@ -16,6 +16,7 @@ // Local Includes #include "../global/global.h" #include "../global/global_cuda.h" +#include "../utils/basic_structs.h" #include "../utils/gpu.hpp" namespace math_utils @@ -98,4 +99,34 @@ inline __device__ __host__ Real SquareMagnitude(Real const &v1, Real const &v2, }; // ========================================================================= +// ===================================================================================================================== +/*! + * \brief Cyclically permute a Vector once. i.e. (x,y,z) becomes (y,z,x) + * + * \param[in,out] vec The vector to permute + */ +inline __device__ __host__ void Cyclic_Permute_Once(hydro_utilities::VectorXYZ &vec) +{ + Real temp = vec.x(); + vec.x() = vec.y(); + vec.y() = vec.z(); + vec.z() = temp; +} +// ===================================================================================================================== + +// ===================================================================================================================== +/*! + * \brief Cyclically permute a Vector twice. i.e. (x,y,z) becomes (z,x,y) + * + * \param[in,out] vec The vector to permute + */ +inline __device__ __host__ void Cyclic_Permute_Twice(hydro_utilities::VectorXYZ &vec) +{ + Real temp = vec.y(); + vec.y() = vec.x(); + vec.x() = vec.z(); + vec.z() = temp; +} +// ===================================================================================================================== + } // namespace math_utils diff --git a/src/utils/math_utilities_tests.cpp b/src/utils/math_utilities_tests.cpp index a49cd8a41..ee0da7ba8 100644 --- a/src/utils/math_utilities_tests.cpp +++ b/src/utils/math_utilities_tests.cpp @@ -13,6 +13,7 @@ // Local Includes #include "../global/global.h" +#include "../utils/basic_structs.h" #include "../utils/math_utilities.h" #include "../utils/testing_utilities.h" @@ -74,4 +75,40 @@ TEST(tALLSquareMagnitude, CorrectInputExpectCorrectOutput) // Now check results testing_utilities::Check_Results(fiducial_square_magnitude, test_square_magnitude, "dot product"); } -// ========================================================================= \ No newline at end of file +// ========================================================================= + +// ========================================================================= +/*! + * \brief Test the math_utils::Cyclic_Permute_Once function + * + */ +TEST(tALLCyclicPermuteOnce, CorrectInputExpectCorrectOutput) +{ + hydro_utilities::VectorXYZ test_vec{1, 2, 3}; + + math_utils::Cyclic_Permute_Once(test_vec); + + // Now check results + testing_utilities::Check_Results(2, test_vec.x(), "Failure in x term"); + testing_utilities::Check_Results(3, test_vec.y(), "Failure in y term"); + testing_utilities::Check_Results(1, test_vec.z(), "Failure in z term"); +} +// ========================================================================= + +// ========================================================================= +/*! + * \brief Test the math_utils::Cyclic_Permute_Twice function + * + */ +TEST(tALLCyclicPermuteTwice, CorrectInputExpectCorrectOutput) +{ + hydro_utilities::VectorXYZ test_vec{1, 2, 3}; + + math_utils::Cyclic_Permute_Twice(test_vec); + + // Now check results + testing_utilities::Check_Results(3, test_vec.x(), "Failure in x term"); + testing_utilities::Check_Results(1, test_vec.y(), "Failure in y term"); + testing_utilities::Check_Results(2, test_vec.z(), "Failure in z term"); +} +// =========================================================================