diff --git a/src/integrators/VL_1D_cuda.cu b/src/integrators/VL_1D_cuda.cu index 421351fba..fe91819dc 100644 --- a/src/integrators/VL_1D_cuda.cu +++ b/src/integrators/VL_1D_cuda.cu @@ -13,8 +13,7 @@ #include "../integrators/VL_1D_cuda.h" #include "../io/io.h" #include "../reconstruction/plm_cuda.h" - #include "../reconstruction/ppmc_cuda.h" - #include "../reconstruction/ppmp_cuda.h" + #include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hllc_cuda.h" @@ -85,12 +84,8 @@ void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Rea #if defined(PLMP) or defined(PLMC) hipLaunchKernelGGL(PLM_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); #endif // PLMP or PLMC - #ifdef PPMP - hipLaunchKernelGGL(PPMP_cuda, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, dt, - gama, 0, n_fields); - #endif - #ifdef PPMC - hipLaunchKernelGGL(PPMC_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + #if defined(PPMP) or defined(PPMC) + hipLaunchKernelGGL(PPM_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); #endif GPU_Error_Check(); diff --git a/src/integrators/VL_2D_cuda.cu b/src/integrators/VL_2D_cuda.cu index c993d695c..e510d3e55 100644 --- a/src/integrators/VL_2D_cuda.cu +++ b/src/integrators/VL_2D_cuda.cu @@ -11,8 +11,7 @@ #include "../hydro/hydro_cuda.h" #include "../integrators/VL_2D_cuda.h" #include "../reconstruction/plm_cuda.h" - #include "../reconstruction/ppmc_cuda.h" - #include "../reconstruction/ppmp_cuda.h" + #include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hllc_cuda.h" @@ -97,16 +96,10 @@ void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_of hipLaunchKernelGGL(PLM_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); #endif // PLMP or PLMC - #ifdef PPMP - hipLaunchKernelGGL(PPMP_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, - dt, gama, 0, n_fields); - hipLaunchKernelGGL(PPMP_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, n_ghost, dy, - dt, gama, 1, n_fields); - #endif // PPMP - #ifdef PPMC - hipLaunchKernelGGL(PPMC_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, + #if defined(PPMP) or defined(PPMC) + hipLaunchKernelGGL(PPM_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, + hipLaunchKernelGGL(PPM_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); #endif // PPMC GPU_Error_Check(); diff --git a/src/integrators/VL_3D_cuda.cu b/src/integrators/VL_3D_cuda.cu index ea31d09c0..09cf2144a 100644 --- a/src/integrators/VL_3D_cuda.cu +++ b/src/integrators/VL_3D_cuda.cu @@ -18,8 +18,7 @@ #include "../mhd/ct_electric_fields.h" #include "../mhd/magnetic_update.h" #include "../reconstruction/plm_cuda.h" - #include "../reconstruction/ppmc_cuda.h" - #include "../reconstruction/ppmp_cuda.h" + #include "../reconstruction/ppm_cuda.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hll_cuda.h" #include "../riemann_solvers/hllc_cuda.h" @@ -245,23 +244,14 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int hipLaunchKernelGGL(PLM_cuda<2>, plm_vl_launch_params.get_numBlocks(), plm_vl_launch_params.get_threadsPerBlock(), 0, 0, dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama); #endif // PLMP or PLMC - #ifdef PPMP - cuda_utilities::AutomaticLaunchParams static const ppmp_launch_params(PPMP_cuda, n_cells); - hipLaunchKernelGGL(PPMP_cuda, ppmp_launch_params.get_numBlocks(), ppmp_launch_params.get_threadsPerBlock(), 0, 0, - dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, dt, gama, 0, n_fields); - hipLaunchKernelGGL(PPMP_cuda, ppmp_launch_params.get_numBlocks(), ppmp_launch_params.get_threadsPerBlock(), 0, 0, - dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, n_ghost, dy, dt, gama, 1, n_fields); - hipLaunchKernelGGL(PPMP_cuda, ppmp_launch_params.get_numBlocks(), ppmp_launch_params.get_threadsPerBlock(), 0, 0, - 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_cuda<0>, n_cells); - hipLaunchKernelGGL(PPMC_cuda<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, dx, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<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, dy, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<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, dz, dt, gama); + #if defined(PPMP) or defined(PPMC) + cuda_utilities::AutomaticLaunchParams static const ppm_launch_params(PPM_cuda<0>, n_cells); + hipLaunchKernelGGL(PPM_cuda<0>, ppm_launch_params.get_numBlocks(), ppm_launch_params.get_threadsPerBlock(), 0, 0, + dev_conserved_half, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + hipLaunchKernelGGL(PPM_cuda<1>, ppm_launch_params.get_numBlocks(), ppm_launch_params.get_threadsPerBlock(), 0, 0, + dev_conserved_half, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); + hipLaunchKernelGGL(PPM_cuda<2>, ppm_launch_params.get_numBlocks(), ppm_launch_params.get_threadsPerBlock(), 0, 0, + dev_conserved_half, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama); #endif // PPMC GPU_Error_Check(); diff --git a/src/integrators/simple_1D_cuda.cu b/src/integrators/simple_1D_cuda.cu index 3137a2f27..f83aff576 100644 --- a/src/integrators/simple_1D_cuda.cu +++ b/src/integrators/simple_1D_cuda.cu @@ -11,8 +11,7 @@ #include "../integrators/simple_1D_cuda.h" #include "../io/io.h" #include "../reconstruction/plm_cuda.h" -#include "../reconstruction/ppmc_cuda.h" -#include "../reconstruction/ppmp_cuda.h" +#include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hllc_cuda.h" @@ -56,13 +55,8 @@ void Simple_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, hipLaunchKernelGGL(PLM_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); GPU_Error_Check(); #endif // PLMP or PLMC -#ifdef PPMP - hipLaunchKernelGGL(PPMP_cuda, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, dt, gama, - 0, n_fields); - GPU_Error_Check(); -#endif -#ifdef PPMC - hipLaunchKernelGGL(PPMC_cuda<0>, dimGrid, dimBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); +#if defined(PPMP) or defined(PPMC) + hipLaunchKernelGGL(PPM_cuda<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 bdfe1e25e..9873379c1 100644 --- a/src/integrators/simple_2D_cuda.cu +++ b/src/integrators/simple_2D_cuda.cu @@ -9,8 +9,7 @@ #include "../hydro/hydro_cuda.h" #include "../integrators/simple_2D_cuda.h" #include "../reconstruction/plm_cuda.h" -#include "../reconstruction/ppmc_cuda.h" -#include "../reconstruction/ppmp_cuda.h" +#include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hllc_cuda.h" @@ -57,15 +56,9 @@ void Simple_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int hipLaunchKernelGGL(PLM_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); hipLaunchKernelGGL(PLM_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); #endif // PLMP or PLMC -#ifdef PPMP - hipLaunchKernelGGL(PPMP_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, dt, - gama, 0, n_fields); - hipLaunchKernelGGL(PPMP_cuda, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, n_ghost, dy, dt, - gama, 1, n_fields); -#endif -#ifdef PPMC - hipLaunchKernelGGL(PPMC_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<1>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); +#if defined(PPMP) or defined(PPMC) + hipLaunchKernelGGL(PPM_cuda<0>, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + hipLaunchKernelGGL(PPM_cuda<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 671e61768..56c575b81 100644 --- a/src/integrators/simple_3D_cuda.cu +++ b/src/integrators/simple_3D_cuda.cu @@ -13,8 +13,7 @@ #include "../integrators/simple_3D_cuda.h" #include "../io/io.h" #include "../reconstruction/plm_cuda.h" - #include "../reconstruction/ppmc_cuda.h" - #include "../reconstruction/ppmp_cuda.h" + #include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction.h" #include "../riemann_solvers/exact_cuda.h" #include "../riemann_solvers/hll_cuda.h" @@ -89,18 +88,10 @@ void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, hipLaunchKernelGGL(PLM_cuda<1>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); hipLaunchKernelGGL(PLM_cuda<2>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama); #endif // PLMP or PLMC - #ifdef PPMP - hipLaunchKernelGGL(PPMP_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, n_ghost, dx, dt, - gama, 0, n_fields); - hipLaunchKernelGGL(PPMP_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, n_ghost, dy, dt, - gama, 1, n_fields); - hipLaunchKernelGGL(PPMP_cuda, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, n_ghost, dz, dt, - gama, 2, n_fields); - #endif // PPMP - #ifdef PPMC - hipLaunchKernelGGL(PPMC_cuda<0>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<1>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); - hipLaunchKernelGGL(PPMC_cuda<2>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lz, Q_Rz, nx, ny, nz, dz, dt, gama); + #if defined(PPMP) or defined(PPMC) + hipLaunchKernelGGL(PPM_cuda<0>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, nx, ny, nz, dx, dt, gama); + hipLaunchKernelGGL(PPM_cuda<1>, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Ly, Q_Ry, nx, ny, nz, dy, dt, gama); + hipLaunchKernelGGL(PPM_cuda<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/plm_cuda.h b/src/reconstruction/plm_cuda.h index c214661a9..e0f43ec89 100644 --- a/src/reconstruction/plm_cuda.h +++ b/src/reconstruction/plm_cuda.h @@ -2,8 +2,8 @@ * \brief Declarations of the cuda plm kernels */ -#ifndef PLMC_CUDA_H -#define PLMC_CUDA_H +#ifndef PLM_CUDA_H +#define PLM_CUDA_H #include "../global/global.h" #include "../grid/grid_enum.h" @@ -314,4 +314,4 @@ auto __device__ __inline__ PLM_Reconstruction(Real *dev_conserved, int const xid } } // namespace reconstruction -#endif // PLMC_CUDA_H +#endif // PLM_CUDA_H diff --git a/src/reconstruction/ppmc_cuda.cu b/src/reconstruction/ppm_cuda.cu similarity index 84% rename from src/reconstruction/ppmc_cuda.cu rename to src/reconstruction/ppm_cuda.cu index a57ea923f..2a99ee83a 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppm_cuda.cu @@ -1,4 +1,4 @@ -/*! \file ppmc_cuda.cu +/*! \file ppm_cuda.cu * \brief Functions definitions for the ppm kernels, using characteristic tracing. Written following Stone et al. 2008. */ @@ -6,7 +6,7 @@ #include "../global/global.h" #include "../global/global_cuda.h" -#include "../reconstruction/ppmc_cuda.h" +#include "../reconstruction/ppm_cuda.h" #include "../reconstruction/reconstruction_internals.h" #include "../utils/gpu.hpp" #include "../utils/hydro_utilities.h" @@ -325,8 +325,8 @@ void __device__ __host__ __inline__ PPM_Characteristic_Evolution(hydro_utilities // ===================================================================================================================== template -__global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real dx, Real dt, Real gamma) +__global__ __launch_bounds__(TPB) void PPM_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, + int ny, int nz, Real dx, Real dt, Real gamma) { // get a thread ID int const thread_id = threadIdx.x + blockIdx.x * blockDim.x; @@ -382,6 +382,7 @@ __global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_ 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); +#ifdef PPMC // Compute the eigenvectors reconstruction::EigenVecs const eigenvectors = reconstruction::Compute_Eigenvectors(cell_i, gamma); @@ -406,32 +407,9 @@ __global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_ reconstruction::Primitive_To_Characteristic(cell_i, cell_ip2, eigenvectors, gamma); // Compute the interface states for each field - reconstruction::Characteristic interface_R_imh_characteristic, interface_L_iph_characteristic; - - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a0, cell_im1_characteristic.a0, cell_i_characteristic.a0, - cell_ip1_characteristic.a0, cell_ip2_characteristic.a0, - interface_L_iph_characteristic.a0, interface_R_imh_characteristic.a0); - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a1, cell_im1_characteristic.a1, cell_i_characteristic.a1, - cell_ip1_characteristic.a1, cell_ip2_characteristic.a1, - interface_L_iph_characteristic.a1, interface_R_imh_characteristic.a1); - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a2, cell_im1_characteristic.a2, cell_i_characteristic.a2, - cell_ip1_characteristic.a2, cell_ip2_characteristic.a2, - interface_L_iph_characteristic.a2, interface_R_imh_characteristic.a2); - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a3, cell_im1_characteristic.a3, cell_i_characteristic.a3, - cell_ip1_characteristic.a3, cell_ip2_characteristic.a3, - interface_L_iph_characteristic.a3, interface_R_imh_characteristic.a3); - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a4, cell_im1_characteristic.a4, cell_i_characteristic.a4, - cell_ip1_characteristic.a4, cell_ip2_characteristic.a4, - interface_L_iph_characteristic.a4, interface_R_imh_characteristic.a4); - -#ifdef MHD - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a5, cell_im1_characteristic.a5, cell_i_characteristic.a5, - cell_ip1_characteristic.a5, cell_ip2_characteristic.a5, - interface_L_iph_characteristic.a5, interface_R_imh_characteristic.a5); - reconstruction::PPM_Single_Variable(cell_im2_characteristic.a6, cell_im1_characteristic.a6, cell_i_characteristic.a6, - cell_ip1_characteristic.a6, cell_ip2_characteristic.a6, - interface_L_iph_characteristic.a6, interface_R_imh_characteristic.a6); -#endif // MHD + auto const [interface_L_iph_characteristic, interface_R_imh_characteristic] = + reconstruction::PPM_Interfaces(cell_im2_characteristic, cell_im1_characteristic, cell_i_characteristic, + cell_ip1_characteristic, cell_ip2_characteristic); // Convert back to primitive variables hydro_utilities::Primitive interface_L_iph = @@ -440,20 +418,24 @@ __global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_ reconstruction::Characteristic_To_Primitive(cell_i, interface_R_imh_characteristic, eigenvectors, gamma); // Compute the interfaces for the variables that don't have characteristics -#ifdef DE + #ifdef DE reconstruction::PPM_Single_Variable(cell_im2.gas_energy_specific, cell_im1.gas_energy_specific, cell_i.gas_energy_specific, cell_ip1.gas_energy_specific, cell_ip2.gas_energy_specific, interface_L_iph.gas_energy_specific, interface_R_imh.gas_energy_specific); -#endif // DE -#ifdef SCALAR + #endif // DE + #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { reconstruction::PPM_Single_Variable(cell_im2.scalar_specific[i], cell_im1.scalar_specific[i], cell_i.scalar_specific[i], cell_ip1.scalar_specific[i], cell_ip2.scalar_specific[i], interface_L_iph.scalar_specific[i], interface_R_imh.scalar_specific[i]); } -#endif // SCALAR + #endif // SCALAR +#else // PPMC + auto [interface_L_iph, interface_R_imh] = + reconstruction::PPM_Interfaces(cell_im2, cell_im1, cell_i, cell_ip1, cell_ip2); +#endif // PPMC // Do the characteristic tracing #ifndef VL @@ -477,13 +459,10 @@ __global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_ 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 PPMC_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); -template __global__ __launch_bounds__(TPB) void PPMC_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); -template __global__ __launch_bounds__(TPB) void PPMC_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); +template __global__ __launch_bounds__(TPB) void PPM_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); +template __global__ __launch_bounds__(TPB) void PPM_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); +template __global__ __launch_bounds__(TPB) void PPM_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); // ===================================================================================================================== diff --git a/src/reconstruction/ppm_cuda.h b/src/reconstruction/ppm_cuda.h new file mode 100644 index 000000000..c84936c5c --- /dev/null +++ b/src/reconstruction/ppm_cuda.h @@ -0,0 +1,32 @@ +/*! \file ppm_cuda.h + * \brief Declarations of the cuda ppm kernels, characteristic reconstruction + * version. */ + +#ifndef PPM_CUDA_H +#define PPM_CUDA_H + +#include "../global/global.h" + +/*! + * \brief Computes the left and right interface states using PPM with limiting in the primitive or characteristic + * variables. This uses the PPM method described in Felker & Stone 2018 "A fourth-order accurate finite volume method + * for ideal MHD via upwind constrained transport". This method computes the 3rd order interface then applies a mixture + * of monoticity constraints from from Colella & Sekora 2008, McCorquodale & Colella 2011, and Colella et al. 2011. We + * found that this newer method and limiters was more stable, less oscillatory, and faster than the method described in + * Stone et al. 2008 which was previously used. 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 + * \param[in] nx The number of cells in the X-direction + * \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 + */ +template +__global__ __launch_bounds__(TPB) void PPM_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, + int ny, int nz, Real dx, Real dt, Real gamma); + +#endif // PPM_CUDA_H diff --git a/src/reconstruction/ppmc_cuda_tests.cu b/src/reconstruction/ppm_cuda_tests.cu similarity index 57% rename from src/reconstruction/ppmc_cuda_tests.cu rename to src/reconstruction/ppm_cuda_tests.cu index 03c4c9223..c43fc2d0f 100644 --- a/src/reconstruction/ppmc_cuda_tests.cu +++ b/src/reconstruction/ppm_cuda_tests.cu @@ -1,6 +1,6 @@ /*! - * \file ppmc_cuda_tests.cu - * \brief Tests for the contents of ppmc_cuda.h and ppmc_cuda.cu + * \file ppm_cuda_tests.cu + * \brief Tests for the contents of ppm_cuda.h and ppm_cuda.cu * */ @@ -17,12 +17,12 @@ // Local Includes #include "../global/global.h" #include "../io/io.h" -#include "../reconstruction/ppmc_cuda.h" +#include "../reconstruction/ppm_cuda.h" #include "../utils/DeviceVector.h" #include "../utils/hydro_utilities.h" #include "../utils/testing_utilities.h" -TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) +TEST(tALLPpmReconstructor, CorrectInputExpectCorrectOutput) { #ifdef DE /// This test doesn't support Dual Energy. It wouldn't be that hard to add support for DE but the DE parts of the @@ -58,6 +58,7 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) // Fiducial Data #ifdef MHD + #ifdef PPMC std::vector> fiducial_interface_left = {{{86, 3.6926886385390683}, {302, 2.3022467009220993}, {518, 2.3207781368125389}, @@ -101,7 +102,53 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) {914, 14.017699282483312}, {1130, 1.5292690020097823}, {1346, -0.12121484974901264}}}; -#else // not MHD + #else // PPMC + std::vector> fiducial_interface_left = {{{86, 3.1608646282711232}, + {302, 0.84444422521258167}, + {518, 1.2459789393105685}, + {734, 2.2721401574613527}, + {950, 7.7508629541568022}, + {1166, 0.54567382624989913}, + {1382, 3.5147238706385462}}, + {{86, 3.6292858956631076}, + {302, 1.8316886259802778}, + {518, 2.2809308293670103}, + {734, 3.6939841768696002}, + {950, 10.405768833830281}, + {1166, 3.5147238706385462}, + {1382, 1.234487908582131}}, + {{86, 3.1608646282711232}, + {302, 0.84444422521258167}, + {518, 1.9865377887960551}, + {734, 1.1540870822905045}, + {950, 4.8971025794015812}, + {1166, 1.234487908582131}, + {1382, 0.54567382624989913}}}; + + std::vector> fiducial_interface_right = {{{301, 0.84444422521258167}, + {85, 3.1608646282711232}, + {733, 2.2721401574613527}, + {517, 3.2701799807980008}, + {949, 10.497902459040514}, + {1165, 0.54567382624989913}, + {1381, 3.5147238706385462}}, + {{80, 2.245959460360242}, + {296, 0.33326844362749702}, + {512, 1.4115388872411132}, + {728, 0.72702830835784316}, + {944, 7.5422056995631559}, + {1160, 3.5147238706385462}, + {1376, 1.234487908582131}}, + {{50, 3.1608646282711232}, + {266, 0.84444422521258167}, + {482, 1.9865377887960551}, + {698, 4.1768690252280765}, + {914, 14.823997016980297}, + {1130, 1.234487908582131}, + {1346, 0.54567382624989913}}}; + #endif // PPMC +#else // not MHD + #ifdef PPMC std::vector> fiducial_interface_left = { {{86, 4.155160222900312}, {302, 1.1624633361407897}, {518, 1.6379195998743412}, {734, 2.9868746414179093}}, {{86, 4.1795874335665655}, {302, 2.1094239978455054}, {518, 2.6811988240843849}, {734, 4.2540957888954054}}, @@ -122,7 +169,29 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) {266, 0.75482470285166003}, {482, 1.7757096932649317}, {698, 3.6101832818706452}}}; -#endif // MHD + #else // PPMC + std::vector> fiducial_interface_left = { + {{86, 3.1608646282711232}, {302, 0.84444422521258167}, {518, 1.2459789393105685}, {734, 2.2721401574613527}}, + {{86, 3.6292858956631076}, {302, 1.8316886259802778}, {518, 2.2809308293670103}, {734, 3.6939841768696002}}, + {{86, 3.1608646282711232}, {302, 0.84444422521258167}, {518, 1.9865377887960551}, {734, 1.1540870822905045}}}; + + std::vector> fiducial_interface_right = {{{54, 3.4283787020401455}, + {85, 3.1608646282711232}, + {301, 0.84444422521258167}, + {517, 3.2701799807980008}, + {733, 2.2721401574613527}}, + {{54, 5.3122571267813665}, + {80, 2.245959460360242}, + {296, 0.33326844362749702}, + {512, 1.4115388872411132}, + {728, 0.72702830835784316}}, + {{50, 3.1608646282711232}, + {54, 3.2010935757366896}, + {266, 0.84444422521258167}, + {482, 1.9865377887960551}, + {698, 4.1768690252280765}}}; + #endif // PPMC +#endif // MHD // Loop over different directions for (size_t direction = 0; direction < 3; direction++) { @@ -133,15 +202,15 @@ TEST(tALLPpmcVLReconstructor, CorrectInputExpectCorrectOutput) // Launch kernel switch (direction) { case 0: - hipLaunchKernelGGL(PPMC_cuda<0>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + hipLaunchKernelGGL(PPM_cuda<0>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), dev_interface_right.data(), nx, ny, nz, 0, 0, gamma); break; case 1: - hipLaunchKernelGGL(PPMC_cuda<1>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + hipLaunchKernelGGL(PPM_cuda<1>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), dev_interface_right.data(), nx, ny, nz, 0, 0, gamma); break; case 2: - hipLaunchKernelGGL(PPMC_cuda<2>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), + hipLaunchKernelGGL(PPM_cuda<2>, dev_grid.size(), 1, 0, 0, dev_grid.data(), dev_interface_left.data(), dev_interface_right.data(), nx, ny, nz, 0, 0, gamma); break; } diff --git a/src/reconstruction/ppmc_cuda.h b/src/reconstruction/ppmc_cuda.h deleted file mode 100644 index 23d1f28db..000000000 --- a/src/reconstruction/ppmc_cuda.h +++ /dev/null @@ -1,32 +0,0 @@ -/*! \file ppmc_cuda.h - * \brief Declarations of the cuda ppm kernels, characteristic reconstruction - * version. */ - -#ifndef PPMC_CUDA_H -#define PPMC_CUDA_H - -#include "../global/global.h" - -/*! - * \brief Computes the left and right interface states using PPM with limiting in the characteristic variables. This - * uses the PPM method described in Felker & Stone 2018 "A fourth-order accurate finite volume method for ideal MHD via - * upwind constrained transport". This method computes the 3rd order interface then applies a mixture of monoticity - * constraints from from Colella & Sekora 2008, McCorquodale & Colella 2011, and Colella et al. 2011. We found that this - * newer method and limiters was more stable, less oscillatory, and faster than the method described in Stone et al. - * 2008 which was previously used. 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 - * \param[in] nx The number of cells in the X-direction - * \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 - */ -template -__global__ __launch_bounds__(TPB) void PPMC_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, - int ny, int nz, Real dx, Real dt, Real gamma); - -#endif // PPMC_CUDA_H diff --git a/src/reconstruction/ppmp_cuda.cu b/src/reconstruction/ppmp_cuda.cu deleted file mode 100644 index 2038f215a..000000000 --- a/src/reconstruction/ppmp_cuda.cu +++ /dev/null @@ -1,775 +0,0 @@ -/*! \file ppmp_cuda.cu - * \brief Definitions of the piecewise parabolic reconstruction (Fryxell 2000) - functions with limiting in the primitive variables. */ - -#ifdef PPMP - - #include - - #include "../global/global.h" - #include "../global/global_cuda.h" - #include "../reconstruction/ppmp_cuda.h" - #include "../utils/gpu.hpp" - - #ifdef DE // PRESSURE_DE - #include "../utils/hydro_utilities.h" - #endif - -// #define STEEPENING -// #define FLATTENING -// Note: Errors when using FLATTENING, need to check the ghost cells - -/*! \fn __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real - *dev_bounds_R, int nx, int ny, int nz, int n_ghost, Real gamma, int dir, int - n_fields) - * \brief When passed a stencil of conserved variables, returns the left and - right boundary values for the interface calculated using ppm with limiting in - the primitive variables. */ -__global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, int nz, - int n_ghost, Real dx, Real dt, Real gamma, int dir, int n_fields) -{ - int n_cells = nx * ny * nz; - int o1, o2, o3; - if (dir == 0) { - o1 = 1; - o2 = 2; - o3 = 3; - } - if (dir == 1) { - o1 = 2; - o2 = 3; - o3 = 1; - } - if (dir == 2) { - o1 = 3; - o2 = 1; - o3 = 2; - } - - // declare primitive variables in the stencil - Real d_i, vx_i, vy_i, vz_i, p_i; - Real d_imo, vx_imo, vy_imo, vz_imo, p_imo; - Real d_ipo, vx_ipo, vy_ipo, vz_ipo, p_ipo; - Real d_imt, vx_imt, vy_imt, vz_imt, p_imt; - Real d_ipt, vx_ipt, vy_ipt, vz_ipt, p_ipt; - #ifdef FLATTENING - Real p_imth, p_ipth; - #endif // FLATTENING - - // declare left and right interface values - Real d_L, vx_L, vy_L, vz_L, p_L; - Real d_R, vx_R, vy_R, vz_R, p_R; - - // declare other variables - Real del_q_imo, del_q_i, del_q_ipo; - - #ifndef VL - // #ifdef CTU - Real cs, cl, cr; // sound speed in cell i, and at left and right boundaries - Real del_d, del_vx, del_vy, del_vz, del_p; // "slope" accross cell i - Real d_6, vx_6, vy_6, vz_6, p_6; - Real beta_m, beta_0, beta_p; - Real alpha_m, alpha_0, alpha_p; - Real lambda_m, lambda_0, lambda_p; // speed of characteristics - Real dL_m, vxL_m, pL_m; - Real dL_0, vyL_0, vzL_0, pL_0; - Real vxL_p, pL_p; - Real vxR_m, pR_m; - Real dR_0, vyR_0, vzR_0, pR_0; - Real dR_p, vxR_p, pR_p; - Real chi_L_m, chi_L_0, chi_L_p; - Real chi_R_m, chi_R_0, chi_R_p; - #endif // CTU - - #ifdef DE - Real ge_i, ge_imo, ge_ipo, ge_imt, ge_ipt, ge_L, ge_R, E_kin, E, dge; - #ifndef VL - // #ifdef CTU - Real del_ge, ge_6, geL_0, geR_0; - #endif // CTU - #endif // DE - - #ifdef SCALAR - Real scalar_i[NSCALARS], scalar_imo[NSCALARS], scalar_ipo[NSCALARS], scalar_imt[NSCALARS], scalar_ipt[NSCALARS]; - Real scalar_L[NSCALARS], scalar_R[NSCALARS]; - #ifndef VL - // #ifdef CTU - Real del_scalar[NSCALARS], scalar_6[NSCALARS], scalarL_0[NSCALARS], scalarR_0[NSCALARS]; - #endif // CTU - #endif // SCALAR - - // get a thread ID - int blockId = blockIdx.x + blockIdx.y * gridDim.x; - int tid = threadIdx.x + blockId * blockDim.x; - int id; - int zid = tid / (nx * ny); - int yid = (tid - zid * nx * ny) / nx; - int xid = tid - zid * nx * ny - yid * nx; - - int xs, xe, ys, ye, zs, ze; - - // - // if (dir == 0) { - // xs = 3; xe = nx-4; - // ys = 0; ye = ny; - // zs = 0; ze = nz; - // } - // if (dir == 1) { - // xs = 0; xe = nx; - // ys = 3; ye = ny-4; - // zs = 0; ze = nz; - // } - // if (dir == 2) { - // xs = 0; xe = nx; - // ys = 0; ye = ny; - // zs = 3; ze = nz-4; - // } - - // Ignore only the 2 ghost cells on each side ( instead of ignoring 3 ghost - // cells on each side ) - if (dir == 0) { - xs = 2; - xe = nx - 3; - ys = 0; - ye = ny; - zs = 0; - ze = nz; - } - if (dir == 1) { - xs = 0; - xe = nx; - ys = 2; - ye = ny - 3; - zs = 0; - ze = nz; - } - if (dir == 2) { - xs = 0; - xe = nx; - ys = 0; - ye = ny; - zs = 2; - ze = nz - 3; - } - - if (xid >= xs && xid < xe && yid >= ys && yid < ye && zid >= zs && zid < ze) { - // load the 5-cell stencil into registers - // cell i - id = xid + yid * nx + zid * nx * ny; - d_i = dev_conserved[id]; - vx_i = dev_conserved[o1 * n_cells + id] / d_i; - vy_i = dev_conserved[o2 * n_cells + id] / d_i; - vz_i = dev_conserved[o3 * n_cells + id] / d_i; - #ifdef DE // PRESSURE_DE - E = dev_conserved[4 * n_cells + id]; - E_kin = 0.5 * d_i * (vx_i * vx_i + vy_i * vy_i + vz_i * vz_i); - dge = dev_conserved[(n_fields - 1) * n_cells + id]; - p_i = hydro_utilities::Get_Pressure_From_DE(E, E - E_kin, dge, gamma); - #else - p_i = (dev_conserved[4 * n_cells + id] - 0.5 * d_i * (vx_i * vx_i + vy_i * vy_i + vz_i * vz_i)) * (gamma - 1.0); - #endif // PRESSURE_DE - p_i = fmax(p_i, (Real)TINY_NUMBER); - #ifdef DE - ge_i = dge / d_i; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_i[i] = dev_conserved[(5 + i) * n_cells + id] / d_i; - } - #endif // SCALAR - // cell i-1 - if (dir == 0) id = xid - 1 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid - 1) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid - 1) * nx * ny; - d_imo = dev_conserved[id]; - vx_imo = dev_conserved[o1 * n_cells + id] / d_imo; - vy_imo = dev_conserved[o2 * n_cells + id] / d_imo; - vz_imo = dev_conserved[o3 * n_cells + id] / d_imo; - #ifdef DE // PRESSURE_DE - E = dev_conserved[4 * n_cells + id]; - E_kin = 0.5 * d_imo * (vx_imo * vx_imo + vy_imo * vy_imo + vz_imo * vz_imo); - dge = dev_conserved[(n_fields - 1) * n_cells + id]; - p_imo = hydro_utilities::Get_Pressure_From_DE(E, E - E_kin, dge, gamma); - #else - p_imo = (dev_conserved[4 * n_cells + id] - 0.5 * d_imo * (vx_imo * vx_imo + vy_imo * vy_imo + vz_imo * vz_imo)) * - (gamma - 1.0); - #endif // PRESSURE_DE - p_imo = fmax(p_imo, (Real)TINY_NUMBER); - #ifdef DE - ge_imo = dge / d_imo; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_imo[i] = dev_conserved[(5 + i) * n_cells + id] / d_imo; - } - #endif // SCALAR - // cell i+1 - if (dir == 0) id = xid + 1 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid + 1) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid + 1) * nx * ny; - d_ipo = dev_conserved[id]; - vx_ipo = dev_conserved[o1 * n_cells + id] / d_ipo; - vy_ipo = dev_conserved[o2 * n_cells + id] / d_ipo; - vz_ipo = dev_conserved[o3 * n_cells + id] / d_ipo; - #ifdef DE // PRESSURE_DE - E = dev_conserved[4 * n_cells + id]; - E_kin = 0.5 * d_ipo * (vx_ipo * vx_ipo + vy_ipo * vy_ipo + vz_ipo * vz_ipo); - dge = dev_conserved[(n_fields - 1) * n_cells + id]; - p_ipo = hydro_utilities::Get_Pressure_From_DE(E, E - E_kin, dge, gamma); - #else - p_ipo = (dev_conserved[4 * n_cells + id] - 0.5 * d_ipo * (vx_ipo * vx_ipo + vy_ipo * vy_ipo + vz_ipo * vz_ipo)) * - (gamma - 1.0); - #endif // PRESSURE_DE - p_ipo = fmax(p_ipo, (Real)TINY_NUMBER); - #ifdef DE - ge_ipo = dge / d_ipo; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_ipo[i] = dev_conserved[(5 + i) * n_cells + id] / d_ipo; - } - #endif // SCALAR - // cell i-2 - if (dir == 0) id = xid - 2 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid - 2) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid - 2) * nx * ny; - d_imt = dev_conserved[id]; - vx_imt = dev_conserved[o1 * n_cells + id] / d_imt; - vy_imt = dev_conserved[o2 * n_cells + id] / d_imt; - vz_imt = dev_conserved[o3 * n_cells + id] / d_imt; - #ifdef DE // PRESSURE_DE - E = dev_conserved[4 * n_cells + id]; - E_kin = 0.5 * d_imt * (vx_imt * vx_imt + vy_imt * vy_imt + vz_imt * vz_imt); - dge = dev_conserved[(n_fields - 1) * n_cells + id]; - p_imt = hydro_utilities::Get_Pressure_From_DE(E, E - E_kin, dge, gamma); - #else - p_imt = (dev_conserved[4 * n_cells + id] - 0.5 * d_imt * (vx_imt * vx_imt + vy_imt * vy_imt + vz_imt * vz_imt)) * - (gamma - 1.0); - #endif // PRESSURE_DE - p_imt = fmax(p_imt, (Real)TINY_NUMBER); - #ifdef DE - ge_imt = dge / d_imt; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_imt[i] = dev_conserved[(5 + i) * n_cells + id] / d_imt; - } - #endif // SCALAR - // cell i+2 - if (dir == 0) id = xid + 2 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid + 2) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid + 2) * nx * ny; - d_ipt = dev_conserved[id]; - vx_ipt = dev_conserved[o1 * n_cells + id] / d_ipt; - vy_ipt = dev_conserved[o2 * n_cells + id] / d_ipt; - vz_ipt = dev_conserved[o3 * n_cells + id] / d_ipt; - #ifdef DE // PRESSURE_DE - E = dev_conserved[4 * n_cells + id]; - E_kin = 0.5 * d_ipt * (vx_ipt * vx_ipt + vy_ipt * vy_ipt + vz_ipt * vz_ipt); - dge = dev_conserved[(n_fields - 1) * n_cells + id]; - p_ipt = hydro_utilities::Get_Pressure_From_DE(E, E - E_kin, dge, gamma); - #else - p_ipt = (dev_conserved[4 * n_cells + id] - 0.5 * d_ipt * (vx_ipt * vx_ipt + vy_ipt * vy_ipt + vz_ipt * vz_ipt)) * - (gamma - 1.0); - #endif // PRESSURE_DE - p_ipt = fmax(p_ipt, (Real)TINY_NUMBER); - #ifdef DE - ge_ipt = dge / d_ipt; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_ipt[i] = dev_conserved[(5 + i) * n_cells + id] / d_ipt; - } - #endif // SCALAR - #ifdef FLATTENING - // cell i-3 - if (dir == 0) id = xid - 3 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid - 3) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid - 3) * nx * ny; - p_imth = - (dev_conserved[4 * n_cells + id] - 0.5 * - (dev_conserved[o1 * n_cells + id] * dev_conserved[o1 * n_cells + id] + - dev_conserved[o2 * n_cells + id] * dev_conserved[o2 * n_cells + id] + - dev_conserved[o3 * n_cells + id] * dev_conserved[o3 * n_cells + id]) / - dev_conserved[id]) * - (gamma - 1.0); - p_imth = fmax(p_imth, (Real)TINY_NUMBER); - // cell i+3 - if (dir == 0) id = xid + 3 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid + 3) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid + 3) * nx * ny; - p_ipth = - (dev_conserved[4 * n_cells + id] - 0.5 * - (dev_conserved[o1 * n_cells + id] * dev_conserved[o1 * n_cells + id] + - dev_conserved[o2 * n_cells + id] * dev_conserved[o2 * n_cells + id] + - dev_conserved[o3 * n_cells + id] * dev_conserved[o3 * n_cells + id]) / - dev_conserved[id]) * - (gamma - 1.0); - p_ipth = fmax(p_imth, (Real)TINY_NUMBER); - #endif // FLATTENING - - // use ppm routines to set cell boundary values (see Fryxell Sec. 3.1.1) - - // Calculate the monotonized slopes for cells imo, i, ipo (density) - del_q_imo = Calculate_Slope(d_imt, d_imo, d_i); - del_q_i = Calculate_Slope(d_imo, d_i, d_ipo); - del_q_ipo = Calculate_Slope(d_i, d_ipo, d_ipt); - - // Calculate the interface values for density - Interface_Values_PPM(d_imo, d_i, d_ipo, del_q_imo, del_q_i, del_q_ipo, &d_L, &d_R); - - // Calculate the monotonized slopes for cells imo, i, ipo (x-velocity) - del_q_imo = Calculate_Slope(vx_imt, vx_imo, vx_i); - del_q_i = Calculate_Slope(vx_imo, vx_i, vx_ipo); - del_q_ipo = Calculate_Slope(vx_i, vx_ipo, vx_ipt); - - // Calculate the interface values for x-velocity - Interface_Values_PPM(vx_imo, vx_i, vx_ipo, del_q_imo, del_q_i, del_q_ipo, &vx_L, &vx_R); - - // Calculate the monotonized slopes for cells imo, i, ipo (y-velocity) - del_q_imo = Calculate_Slope(vy_imt, vy_imo, vy_i); - del_q_i = Calculate_Slope(vy_imo, vy_i, vy_ipo); - del_q_ipo = Calculate_Slope(vy_i, vy_ipo, vy_ipt); - - // Calculate the interface values for y-velocity - Interface_Values_PPM(vy_imo, vy_i, vy_ipo, del_q_imo, del_q_i, del_q_ipo, &vy_L, &vy_R); - - // Calculate the monotonized slopes for cells imo, i, ipo (z-velocity) - del_q_imo = Calculate_Slope(vz_imt, vz_imo, vz_i); - del_q_i = Calculate_Slope(vz_imo, vz_i, vz_ipo); - del_q_ipo = Calculate_Slope(vz_i, vz_ipo, vz_ipt); - - // Calculate the interface values for z-velocity - Interface_Values_PPM(vz_imo, vz_i, vz_ipo, del_q_imo, del_q_i, del_q_ipo, &vz_L, &vz_R); - - // Calculate the monotonized slopes for cells imo, i, ipo (pressure) - del_q_imo = Calculate_Slope(p_imt, p_imo, p_i); - del_q_i = Calculate_Slope(p_imo, p_i, p_ipo); - del_q_ipo = Calculate_Slope(p_i, p_ipo, p_ipt); - - // Calculate the interface values for pressure - Interface_Values_PPM(p_imo, p_i, p_ipo, del_q_imo, del_q_i, del_q_ipo, &p_L, &p_R); - - #ifdef DE - // Calculate the monotonized slopes for cells imo, i, ipo (internal energy) - del_q_imo = Calculate_Slope(ge_imt, ge_imo, ge_i); - del_q_i = Calculate_Slope(ge_imo, ge_i, ge_ipo); - del_q_ipo = Calculate_Slope(ge_i, ge_ipo, ge_ipt); - - // Calculate the interface values for internal energy - Interface_Values_PPM(ge_imo, ge_i, ge_ipo, del_q_imo, del_q_i, del_q_ipo, &ge_L, &ge_R); - #endif // DE - - #ifdef SCALAR - // Calculate the monotonized slopes for cells imo, i, ipo (passive scalars) - for (int i = 0; i < NSCALARS; i++) { - del_q_imo = Calculate_Slope(scalar_imt[i], scalar_imo[i], scalar_i[i]); - del_q_i = Calculate_Slope(scalar_imo[i], scalar_i[i], scalar_ipo[i]); - del_q_ipo = Calculate_Slope(scalar_i[i], scalar_ipo[i], scalar_ipt[i]); - - // Calculate the interface values for the passive scalars - Interface_Values_PPM(scalar_imo[i], scalar_i[i], scalar_ipo[i], del_q_imo, del_q_i, del_q_ipo, &scalar_L[i], - &scalar_R[i]); - } - #endif // SCALAR - - #ifdef STEEPENING - Real d2_rho_imo, d2_rho_ipo, eta_i; - // check for contact discontinuities & steepen if necessary (see Fryxell - // Sec 3.1.2) if condition 4 (Fryxell Eqn 37) (Colella Eqn 1.16.5) is true, - // check further conditions, otherwise do nothing - if ((fabs(d_ipo - d_imo) / fmin(d_ipo, d_imo)) > 0.01) { - // calculate the second derivative of the density in the imo and ipo cells - d2_rho_imo = calc_d2_rho(d_imt, d_imo, d_i, dx); - d2_rho_ipo = calc_d2_rho(d_i, d_ipo, d_ipt, dx); - // if condition 1 (Fryxell Eqn 38) (Colella Eqn 1.16.5) is true, check - // further conditions, otherwise do nothing - if ((d2_rho_imo * d2_rho_ipo) < 0) { - // calculate condition 5, pressure vs density jumps (Fryxell Eqn 39) - // (Colella Eqn 3.2) if c5 is true, set value of eta for discontinuity - // steepening - if ((fabs(p_ipo - p_imo) / fmin(p_ipo, p_imo)) < 0.1 * gamma * (fabs(d_ipo - d_imo) / fmin(d_ipo, d_imo))) { - // calculate first eta value (Fryxell Eqn 36) (Colella Eqn 1.16.5) - eta_i = calc_eta(d2_rho_imo, d2_rho_ipo, dx, d_imo, d_ipo); - // calculate steepening coefficient (Fryxell Eqn 40) (Colella - // Eqn 1.16) - eta_i = fmax(0, fmin(20 * (eta_i - 0.05), 1)); - - // calculate new left and right interface variables using monotonized - // slopes - del_q_imo = Calculate_Slope(d_imt, d_imo, d_i); - del_q_ipo = Calculate_Slope(d_i, d_ipo, d_ipt); - - // replace left and right interface values of density (Colella - // Eqn 1.14, 1.15) - d_L = d_L * (1 - eta_i) + (d_imo + 0.5 * del_q_imo) * eta_i; - d_R = d_R * (1 - eta_i) + (d_ipo - 0.5 * del_q_ipo) * eta_i; - } - } - } - #endif // STEEPENING - - #ifdef FLATTENING - Real F_imo, F_i, F_ipo; - // flatten shock fronts that are too narrow (see Fryxell Sec 3.1.3) - // calculate the shock steepness parameters (Fryxell Eqn 43) - // calculate the dimensionless flattening coefficients (Fryxell Eqn 45) - F_imo = fmax(0, fmin(1, 10 * (((p_i - p_imt) / (p_ipo - p_imth)) - 0.75))); - F_i = fmax(0, fmin(1, 10 * (((p_ipo - p_imo) / (p_ipt - p_imt)) - 0.75))); - F_ipo = fmax(0, fmin(1, 10 * (((p_ipt - p_i) / (p_ipth - p_imo)) - 0.75))); - // ensure that we are encountering a shock (Fryxell Eqns 46 & 47) - if (fabs(p_i - p_imt) / fmin(p_i, p_imt) < 1. / 3.) { - F_imo = 0; - } - if (fabs(p_ipo - p_imo) / fmin(p_ipo, p_imo) < 1. / 3.) { - F_i = 0; - } - if (fabs(p_ipt - p_i) / fmin(p_ipt, p_i) < 1. / 3.) { - F_ipo = 0; - } - if (vx_i - vx_imt > 0) { - F_imo = 0; - } - if (vx_ipo - vx_imo > 0) { - F_i = 0; - } - if (vx_ipt - vx_i > 0) { - F_ipo = 0; - } - // set the flattening coefficient (Fryxell Eqn 48) - if (p_ipo - p_imo < 0) { - F_i = fmax(F_i, F_ipo); - } else { - F_i = fmax(F_i, F_imo); - } - // modify the interface values - d_L = F_i * d_i + (1 - F_i) * d_L; - vx_L = F_i * vx_i + (1 - F_i) * vx_L; - vy_L = F_i * vy_i + (1 - F_i) * vy_L; - vz_L = F_i * vz_i + (1 - F_i) * vz_L; - p_L = F_i * p_i + (1 - F_i) * p_L; - #ifdef DE - ge_L = F_i * ge_i + (1 - F_i) * ge_L; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_L[i] = F_i * scalar_i[i] + (1 - F_i) * scalar_L[i]; - } - #endif // SCALAR - d_R = F_i * d_i + (1 - F_i) * d_R; - vx_R = F_i * vx_i + (1 - F_i) * vx_R; - vy_R = F_i * vy_i + (1 - F_i) * vy_R; - vz_R = F_i * vz_i + (1 - F_i) * vz_R; - p_R = F_i * p_i + (1 - F_i) * p_R; - #ifdef DE - ge_R = F_i * ge_i + (1 - F_i) * ge_R; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_R[i] = F_i * scalar_i[i] + (1 - F_i) * scalar_R[i]; - } - #endif // SCALAR - #endif // FLATTENING - - #ifndef VL - // #ifdef CTU - // compute sound speed in cell i - cs = sqrt(gamma * p_i / d_i); - - // compute a first guess at the left and right states by taking the average - // under the characteristic on each side that has the largest speed - - // recompute slope across cell for each variable Fryxell Eqn 29 - del_d = d_R - d_L; - del_vx = vx_R - vx_L; - del_vy = vy_R - vy_L; - del_vz = vz_R - vz_L; - del_p = p_R - p_L; - #ifdef DE - del_ge = ge_R - ge_L; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - del_scalar[i] = scalar_R[i] - scalar_L[i]; - } - #endif // SCALAR - - d_6 = 6.0 * (d_i - 0.5 * (d_L + d_R)); // Fryxell Eqn 30 - vx_6 = 6.0 * (vx_i - 0.5 * (vx_L + vx_R)); // Fryxell Eqn 30 - vy_6 = 6.0 * (vy_i - 0.5 * (vy_L + vy_R)); // Fryxell Eqn 30 - vz_6 = 6.0 * (vz_i - 0.5 * (vz_L + vz_R)); // Fryxell Eqn 30 - p_6 = 6.0 * (p_i - 0.5 * (p_L + p_R)); // Fryxell Eqn 30 - #ifdef DE - ge_6 = 6.0 * (ge_i - 0.5 * (ge_L + ge_R)); // Fryxell Eqn 30 - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_6[i] = 6.0 * (scalar_i[i] - 0.5 * (scalar_L[i] + scalar_R[i])); // Fryxell Eqn 30 - } - #endif // SCALAR - - // set speed of characteristics (v-c, v, v+c) using average values of v and - // c - lambda_m = vx_i - cs; - lambda_0 = vx_i; - lambda_p = vx_i + cs; - - // calculate betas (for left state guesses) - beta_m = fmax((lambda_m * dt / dx), 0.0); // Fryxell Eqn 59 - beta_0 = fmax((lambda_0 * dt / dx), 0.0); // Fryxell Eqn 59 - beta_p = fmax((lambda_p * dt / dx), 0.0); // Fryxell Eqn 59 - - // calculate alphas (for right state guesses) - alpha_m = fmax((-lambda_m * dt / dx), 0.0); // Fryxell Eqn 61 - alpha_0 = fmax((-lambda_0 * dt / dx), 0.0); // Fryxell Eqn 61 - alpha_p = fmax((-lambda_p * dt / dx), 0.0); // Fryxell Eqn 61 - - // average values under characteristics for left interface (Fryxell Eqn 60) - dL_m = d_L + 0.5 * alpha_m * (del_d + d_6 * (1 - (2. / 3.) * alpha_m)); - vxL_m = vx_L + 0.5 * alpha_m * (del_vx + vx_6 * (1 - (2. / 3.) * alpha_m)); - pL_m = p_L + 0.5 * alpha_m * (del_p + p_6 * (1 - (2. / 3.) * alpha_m)); - dL_0 = d_L + 0.5 * alpha_0 * (del_d + d_6 * (1 - (2. / 3.) * alpha_0)); - vyL_0 = vy_L + 0.5 * alpha_0 * (del_vy + vy_6 * (1 - (2. / 3.) * alpha_0)); - vzL_0 = vz_L + 0.5 * alpha_0 * (del_vz + vz_6 * (1 - (2. / 3.) * alpha_0)); - #ifdef DE - geL_0 = ge_L + 0.5 * alpha_0 * (del_ge + ge_6 * (1 - (2. / 3.) * alpha_0)); - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalarL_0[i] = scalar_L[i] + 0.5 * alpha_0 * (del_scalar[i] + scalar_6[i] * (1 - (2. / 3.) * alpha_0)); - } - #endif // SCALAR - pL_0 = p_L + 0.5 * alpha_0 * (del_p + p_6 * (1 - (2. / 3.) * alpha_0)); - vxL_p = vx_L + 0.5 * alpha_p * (del_vx + vx_6 * (1 - (2. / 3.) * alpha_p)); - pL_p = p_L + 0.5 * alpha_p * (del_p + p_6 * (1 - (2. / 3.) * alpha_p)); - - // average values under characteristics for right interface (Fryxell Eqn 58) - vxR_m = vx_R - 0.5 * beta_m * (del_vx - vx_6 * (1 - (2. / 3.) * beta_m)); - pR_m = p_R - 0.5 * beta_m * (del_p - p_6 * (1 - (2. / 3.) * beta_m)); - dR_0 = d_R - 0.5 * beta_0 * (del_d - d_6 * (1 - (2. / 3.) * beta_0)); - vyR_0 = vy_R - 0.5 * beta_0 * (del_vy - vy_6 * (1 - (2. / 3.) * beta_0)); - vzR_0 = vz_R - 0.5 * beta_0 * (del_vz - vz_6 * (1 - (2. / 3.) * beta_0)); - #ifdef DE - geR_0 = ge_R - 0.5 * beta_0 * (del_ge - ge_6 * (1 - (2. / 3.) * beta_0)); - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalarR_0[i] = scalar_R[i] - 0.5 * beta_0 * (del_scalar[i] - scalar_6[i] * (1 - (2. / 3.) * beta_0)); - } - #endif // SCALAR - pR_0 = p_R - 0.5 * beta_0 * (del_p - p_6 * (1 - (2. / 3.) * beta_0)); - dR_p = d_R - 0.5 * beta_p * (del_d - d_6 * (1 - (2. / 3.) * beta_p)); - vxR_p = vx_R - 0.5 * beta_p * (del_vx - vx_6 * (1 - (2. / 3.) * beta_p)); - pR_p = p_R - 0.5 * beta_p * (del_p - p_6 * (1 - (2. / 3.) * beta_p)); - - // as a first guess, use characteristics with the largest speeds - // for transverse velocities, use the 0 characteristic - // left - d_L = dL_m; - vx_L = vxL_m; - vy_L = vyL_0; - vz_L = vzL_0; - p_L = pL_m; - #ifdef DE - ge_L = geL_0; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_L[i] = scalarL_0[i]; - } - #endif // SCALAR - // right - d_R = dR_p; - vx_R = vxR_p; - vy_R = vyR_0; - vz_R = vzR_0; - p_R = pR_p; - #ifdef DE - ge_R = geR_0; - #endif // DE - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - scalar_R[i] = scalarR_0[i]; - } - #endif // SCALAR - - // correct these initial guesses by taking into account the number of - // characteristics on each side of the interface - - // calculate the 'guess' sound speeds - cl = sqrt(gamma * p_L / d_L); - cr = sqrt(gamma * p_R / d_L); - - // calculate the chi values (Fryxell Eqns 62 & 63) - chi_L_m = 1. / (2 * d_L * cl) * (vx_L - vxL_m - (p_L - pL_m) / (d_L * cl)); - chi_L_p = -1. / (2 * d_L * cl) * (vx_L - vxL_p + (p_L - pL_p) / (d_L * cl)); - chi_L_0 = (p_L - pL_0) / (d_L * d_L * cl * cl) + 1. / d_L - 1. / dL_0; - chi_R_m = 1. / (2 * d_R * cr) * (vx_R - vxR_m - (p_R - pR_m) / (d_R * cr)); - chi_R_p = -1. / (2 * d_R * cr) * (vx_R - vxR_p + (p_R - pR_p) / (d_R * cr)); - chi_R_0 = (p_R - pR_0) / (d_R * d_R * cr * cr) + 1. / d_R - 1. / dR_0; - - // set chi to 0 if characteristic velocity has the wrong sign (Fryxell Eqn - // 64) - if (lambda_m >= 0) { - chi_L_m = 0; - } - if (lambda_0 >= 0) { - chi_L_0 = 0; - } - if (lambda_p >= 0) { - chi_L_p = 0; - } - if (lambda_m <= 0) { - chi_R_m = 0; - } - if (lambda_0 <= 0) { - chi_R_0 = 0; - } - if (lambda_p <= 0) { - chi_R_p = 0; - } - - // use the chi values to correct the initial guesses and calculate final - // input states - p_L = p_L + (d_L * d_L * cl * cl) * (chi_L_p + chi_L_m); - vx_L = vx_L + d_L * cl * (chi_L_p - chi_L_m); - d_L = pow(((1.0 / d_L) - (chi_L_m + chi_L_0 + chi_L_p)), -1); - p_R = p_L + (d_R * d_R * cr * cr) * (chi_R_p + chi_R_m); - vx_R = vx_R + d_R * cr * (chi_R_p - chi_R_m); - d_R = pow(((1.0 / d_R) - (chi_R_m + chi_R_0 + chi_R_p)), -1); - #endif // CTU - - // Apply mimimum constraints - d_L = fmax(d_L, (Real)TINY_NUMBER); - d_R = fmax(d_R, (Real)TINY_NUMBER); - p_L = fmax(p_L, (Real)TINY_NUMBER); - p_R = fmax(p_R, (Real)TINY_NUMBER); - - // Convert the left and right states in the primitive to the conserved - // variables send final values back from kernel bounds_R refers to the right - // side of the i-1/2 interface - if (dir == 0) id = xid - 1 + yid * nx + zid * nx * ny; - if (dir == 1) id = xid + (yid - 1) * nx + zid * nx * ny; - if (dir == 2) id = xid + yid * nx + (zid - 1) * nx * ny; - dev_bounds_R[id] = d_L; - dev_bounds_R[o1 * n_cells + id] = d_L * vx_L; - dev_bounds_R[o2 * n_cells + id] = d_L * vy_L; - dev_bounds_R[o3 * n_cells + id] = d_L * vz_L; - dev_bounds_R[4 * n_cells + id] = p_L / (gamma - 1.0) + 0.5 * d_L * (vx_L * vx_L + vy_L * vy_L + vz_L * vz_L); - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - dev_bounds_R[(5 + i) * n_cells + id] = d_L * scalar_L[i]; - } - #endif // SCALAR - #ifdef DE - dev_bounds_R[(n_fields - 1) * n_cells + id] = d_L * ge_L; - #endif // DE - // bounds_L refers to the left side of the i+1/2 interface - id = xid + yid * nx + zid * nx * ny; - dev_bounds_L[id] = d_R; - dev_bounds_L[o1 * n_cells + id] = d_R * vx_R; - dev_bounds_L[o2 * n_cells + id] = d_R * vy_R; - dev_bounds_L[o3 * n_cells + id] = d_R * vz_R; - dev_bounds_L[4 * n_cells + id] = p_R / (gamma - 1.0) + 0.5 * d_R * (vx_R * vx_R + vy_R * vy_R + vz_R * vz_R); - #ifdef SCALAR - for (int i = 0; i < NSCALARS; i++) { - dev_bounds_L[(5 + i) * n_cells + id] = d_R * scalar_R[i]; - } - #endif // SCALAR - #ifdef DE - dev_bounds_L[(n_fields - 1) * n_cells + id] = d_R * ge_R; - #endif // DE - } -} - -/*! \fn __device__ Real Calculate_Slope(Real q_imo, Real q_i, Real q_ipo) - * \brief Calculates the limited slope across a cell.*/ -__device__ Real Calculate_Slope(Real q_imo, Real q_i, Real q_ipo) -{ - Real del_q_L, del_q_R, del_q_C, del_q_G; - Real lim_slope_a, lim_slope_b, del_q_m; - - // Compute the left, right, and centered differences of the primitive - // variables Note that here L and R refer to locations relative to the cell - // center - - // left - del_q_L = q_i - q_imo; - // right - del_q_R = q_ipo - q_i; - // centered - del_q_C = 0.5 * (q_ipo - q_imo); - // Van Leer - if (del_q_L * del_q_R > 0.0) { - del_q_G = 2.0 * del_q_L * del_q_R / (del_q_L + del_q_R); - } else { - del_q_G = 0.0; - } - - // Monotonize the differences - lim_slope_a = fmin(fabs(del_q_L), fabs(del_q_R)); - lim_slope_b = fmin(fabs(del_q_C), fabs(del_q_G)); - - // Minmod limiter - // del_q_m = sgn_CUDA(del_q_C)*fmin(2.0*lim_slope_a, fabs(del_q_C)); - - // Van Leer limiter - del_q_m = sgn_CUDA(del_q_C) * fmin((Real)2.0 * lim_slope_a, lim_slope_b); - - return del_q_m; -} - -/*! \fn __device__ void Interface_Values_PPM(Real q_imo, Real q_i, Real q_ipo, - Real del_q_imo, Real del_q_i, Real del_q_ipo, Real *q_L, Real *q_R) - * \brief Calculates the left and right interface values for a cell using - parabolic reconstruction in the primitive variables with limited slopes - provided. Applies further monotonicity constraints.*/ -__device__ void Interface_Values_PPM(Real q_imo, Real q_i, Real q_ipo, Real del_q_imo, Real del_q_i, Real del_q_ipo, - Real *q_L, Real *q_R) -{ - // Calculate the left and right interface values using the limited slopes - *q_L = 0.5 * (q_i + q_imo) - (1.0 / 6.0) * (del_q_i - del_q_imo); - *q_R = 0.5 * (q_ipo + q_i) - (1.0 / 6.0) * (del_q_ipo - del_q_i); - - // Apply further monotonicity constraints to ensure interface values lie - // between neighboring cell-centered values - - // local maximum or minimum criterion (Fryxell Eqn 52, Fig 11) - if ((*q_R - q_i) * (q_i - *q_L) <= 0) *q_L = *q_R = q_i; - - // steep gradient criterion (Fryxell Eqn 53, Fig 12) - if (6.0 * (*q_R - *q_L) * (q_i - 0.5 * (*q_L + *q_R)) > (*q_R - *q_L) * (*q_R - *q_L)) { - *q_L = 3.0 * q_i - 2.0 * (*q_R); - } - if (6.0 * (*q_R - *q_L) * (q_i - 0.5 * (*q_L + *q_R)) < -(*q_R - *q_L) * (*q_R - *q_L)) { - *q_R = 3.0 * q_i - 2.0 * (*q_L); - } - - *q_L = fmax(fmin(q_i, q_imo), *q_L); - *q_L = fmin(fmax(q_i, q_imo), *q_L); - *q_R = fmax(fmin(q_i, q_ipo), *q_R); - *q_R = fmin(fmax(q_i, q_ipo), *q_R); -} - -/*! \fn calc_d2_rho - * \brief Returns the second derivative of rho across zone i. (Fryxell Eqn 35) - */ -__device__ Real calc_d2_rho(Real rho_imo, Real rho_i, Real rho_ipo, Real dx) -{ - return (1. / (6 * dx * dx)) * (rho_ipo - 2 * rho_i + rho_imo); -} - -/*! \fn calc_eta - * \brief Returns a dimensionless quantity relating the 1st and 3rd derivatives - See Fryxell Eqn 36. */ -__device__ Real calc_eta(Real d2rho_imo, Real d2rho_ipo, Real dx, Real rho_imo, Real rho_ipo) -{ - Real A, B; - - A = (d2rho_ipo - d2rho_imo) * dx * dx; - B = 1.0 / (rho_ipo - rho_imo); - - return -A * B; -} - -#endif // PPMP diff --git a/src/reconstruction/ppmp_cuda.h b/src/reconstruction/ppmp_cuda.h deleted file mode 100644 index 064d328fa..000000000 --- a/src/reconstruction/ppmp_cuda.h +++ /dev/null @@ -1,40 +0,0 @@ -/*! \file ppmp_cuda.h - * \brief Declarations of the cuda ppmp kernels. */ - -#ifndef PPMP_CUDA_H -#define PPMP_CUDA_H - -#include "../global/global.h" - -/*! \fn __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real - *dev_bounds_R, int nx, int ny, int nz, int n_ghost, Real dx, Real dt, Real - gamma, int dir, int n_fields) - * \brief When passed a stencil of conserved variables, returns the left and - right boundary values for the interface calculated using ppm with limiting in - the primitive variables. */ -__global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bounds_R, int nx, int ny, int nz, - int n_ghost, Real dx, Real dt, Real gamma, int dir, int n_fields); - -/*! \fn __device__ Real Calculate_Slope(Real q_imo, Real q_i, Real q_ipo) - * \brief Calculates the limited slope across a cell.*/ -__device__ Real Calculate_Slope(Real q_imo, Real q_i, Real q_ipo); - -/*! \fn __device__ void Interface_Values_PPM(Real q_imo, Real q_i, Real q_ipo, - Real *q_L, Real *q_R) - * \brief Calculates the left and right interface values for a cell using - parabolic reconstruction in the primitive variables with limited slopes - provided. Applies further monotonicity constraints.*/ -__device__ void Interface_Values_PPM(Real q_imo, Real q_i, Real q_ipo, Real del_q_imo, Real del_q_i, Real del_q_ipo, - Real *q_L, Real *q_R); - -/*! \fn calc_d2_rho - * \brief Returns the second derivative of rho across zone i. (Fryxell Eqn 35) - */ -__device__ Real calc_d2_rho(Real rho_imo, Real rho_i, Real rho_ipo, Real dx); - -/*! \fn calc_eta - * \brief Returns a dimensionless quantity relating the 1st and 3rd derivatives - See Fryxell Eqn 36. */ -__device__ Real calc_eta(Real d2rho_imo, Real d2rho_ipo, Real dx, Real rho_imo, Real rho_ipo); - -#endif // PPMP_CUDA_H diff --git a/src/reconstruction/reconstruction_internals.h b/src/reconstruction/reconstruction_internals.h index ffeb03103..43120e73e 100644 --- a/src/reconstruction/reconstruction_internals.h +++ b/src/reconstruction/reconstruction_internals.h @@ -735,6 +735,107 @@ void __device__ __host__ __inline__ PPM_Single_Variable(Real const &cell_im2, Re } // ===================================================================================================================== +// ===================================================================================================================== +/*! + * \brief Compute the primitive PPM interfaces. Calls PPM_Single_Variable on each field + * + * \param[in] cell_im2 The state of the cell at i-2 + * \param[in] cell_im1 The state of the cell at i-1 + * \param[in] cell_i The state of the cell at i + * \param[in] cell_ip1 The state of the cell at i+1 + * \param[in] cell_ip2 The state of the cell at i+2 + * \return auto The left interface at i+1/2 and the right interface at i-1/2 in that order + */ +auto __device__ __host__ __inline__ PPM_Interfaces(hydro_utilities::Primitive const &cell_im2, + hydro_utilities::Primitive const &cell_im1, + hydro_utilities::Primitive const &cell_i, + hydro_utilities::Primitive const &cell_ip1, + hydro_utilities::Primitive const &cell_ip2) +{ + hydro_utilities::Primitive interface_R_imh, interface_L_iph; + + reconstruction::PPM_Single_Variable(cell_im2.density, cell_im1.density, cell_i.density, cell_ip1.density, + cell_ip2.density, interface_L_iph.density, interface_R_imh.density); + reconstruction::PPM_Single_Variable(cell_im2.velocity.x, cell_im1.velocity.x, cell_i.velocity.x, cell_ip1.velocity.x, + cell_ip2.velocity.x, interface_L_iph.velocity.x, interface_R_imh.velocity.x); + reconstruction::PPM_Single_Variable(cell_im2.velocity.y, cell_im1.velocity.y, cell_i.velocity.y, cell_ip1.velocity.y, + cell_ip2.velocity.y, interface_L_iph.velocity.y, interface_R_imh.velocity.y); + reconstruction::PPM_Single_Variable(cell_im2.velocity.z, cell_im1.velocity.z, cell_i.velocity.z, cell_ip1.velocity.z, + cell_ip2.velocity.z, interface_L_iph.velocity.z, interface_R_imh.velocity.z); + reconstruction::PPM_Single_Variable(cell_im2.pressure, cell_im1.pressure, cell_i.pressure, cell_ip1.pressure, + cell_ip2.pressure, interface_L_iph.pressure, interface_R_imh.pressure); + +#ifdef MHD + reconstruction::PPM_Single_Variable(cell_im2.magnetic.y, cell_im1.magnetic.y, cell_i.magnetic.y, cell_ip1.magnetic.y, + cell_ip2.magnetic.y, interface_L_iph.magnetic.y, interface_R_imh.magnetic.y); + reconstruction::PPM_Single_Variable(cell_im2.magnetic.z, cell_im1.magnetic.z, cell_i.magnetic.z, cell_ip1.magnetic.z, + cell_ip2.magnetic.z, interface_L_iph.magnetic.z, interface_R_imh.magnetic.z); +#endif // MHD + +#ifdef DE + reconstruction::PPM_Single_Variable(cell_im2.gas_energy_specific, cell_im1.gas_energy_specific, + cell_i.gas_energy_specific, cell_ip1.gas_energy_specific, + cell_ip2.gas_energy_specific, interface_L_iph.gas_energy_specific, + interface_R_imh.gas_energy_specific); +#endif // DE +#ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + reconstruction::PPM_Single_Variable(cell_im2.scalar_specific[i], cell_im1.scalar_specific[i], + cell_i.scalar_specific[i], cell_ip1.scalar_specific[i], + cell_ip2.scalar_specific[i], interface_L_iph.scalar_specific[i], + interface_R_imh.scalar_specific[i]); + } +#endif // DE + + struct LocalReturnStruct { + hydro_utilities::Primitive left, right; + }; + return LocalReturnStruct{interface_L_iph, interface_R_imh}; +} +// ===================================================================================================================== + +// ===================================================================================================================== +/*! + * \brief Compute the characteristic PPM interfaces. Calls PPM_Single_Variable on each field + * + * \param[in] cell_im2 The state of the cell at i-2 + * \param[in] cell_im1 The state of the cell at i-1 + * \param[in] cell_i The state of the cell at i + * \param[in] cell_ip1 The state of the cell at i+1 + * \param[in] cell_ip2 The state of the cell at i+2 + * \return auto The left interface at i+1/2 and the right interface at i-1/2 in that order + */ +auto __device__ __host__ __inline__ PPM_Interfaces(Characteristic const &cell_im2, Characteristic const &cell_im1, + Characteristic const &cell_i, Characteristic const &cell_ip1, + Characteristic const &cell_ip2) +{ + Characteristic interface_R_imh, interface_L_iph; + + reconstruction::PPM_Single_Variable(cell_im2.a0, cell_im1.a0, cell_i.a0, cell_ip1.a0, cell_ip2.a0, interface_L_iph.a0, + interface_R_imh.a0); + reconstruction::PPM_Single_Variable(cell_im2.a1, cell_im1.a1, cell_i.a1, cell_ip1.a1, cell_ip2.a1, interface_L_iph.a1, + interface_R_imh.a1); + reconstruction::PPM_Single_Variable(cell_im2.a2, cell_im1.a2, cell_i.a2, cell_ip1.a2, cell_ip2.a2, interface_L_iph.a2, + interface_R_imh.a2); + reconstruction::PPM_Single_Variable(cell_im2.a3, cell_im1.a3, cell_i.a3, cell_ip1.a3, cell_ip2.a3, interface_L_iph.a3, + interface_R_imh.a3); + reconstruction::PPM_Single_Variable(cell_im2.a4, cell_im1.a4, cell_i.a4, cell_ip1.a4, cell_ip2.a4, interface_L_iph.a4, + interface_R_imh.a4); + +#ifdef MHD + reconstruction::PPM_Single_Variable(cell_im2.a5, cell_im1.a5, cell_i.a5, cell_ip1.a5, cell_ip2.a5, interface_L_iph.a5, + interface_R_imh.a5); + reconstruction::PPM_Single_Variable(cell_im2.a6, cell_im1.a6, cell_i.a6, cell_ip1.a6, cell_ip2.a6, interface_L_iph.a6, + interface_R_imh.a6); +#endif // MHD + + struct LocalReturnStruct { + Characteristic left, right; + }; + return LocalReturnStruct{interface_L_iph, interface_R_imh}; +} +// ===================================================================================================================== + // ===================================================================================================================== /*! * \brief Write the interface data to the appropriate arrays diff --git a/src/utils/error_handling.cpp b/src/utils/error_handling.cpp index e55cbc951..db0bf676a 100644 --- a/src/utils/error_handling.cpp +++ b/src/utils/error_handling.cpp @@ -95,11 +95,6 @@ void Check_Configuration(Parameters const& P) #error "MHD only supports the HLLD Riemann Solver" #endif //! HLLD or EXACT or ROE or HLL or HLLC - // May only use certain reconstructions - #if ((defined(PCM) + defined(PLMP) + defined(PLMC) + defined(PPMC)) != 1) || defined(PPMP) - #error "MHD does not support PPMP reconstruction." - #endif // Reconstruction check - // must have HDF5 #if defined(OUTPUT) and (not defined(HDF5)) #error "MHD only supports HDF5 output"