From 497df2eeedc50803bbbc0197ad150824e160d735 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 18 Jun 2024 10:35:15 -0400 Subject: [PATCH] Remove `_specific` from struct member variables --- src/reconstruction/pcm_cuda.h | 4 +- src/reconstruction/plmc_cuda.cu | 309 ++++++++++++++++++ src/reconstruction/ppmc_cuda.cu | 64 ++-- src/reconstruction/reconstruction_internals.h | 35 +- src/riemann_solvers/exact_cuda.cu | 16 +- src/riemann_solvers/hll_cuda.cu | 38 +-- src/riemann_solvers/hllc_cuda.cu | 26 +- src/riemann_solvers/hlld_cuda.cu | 10 +- src/riemann_solvers/roe_cuda.cu | 28 +- src/utils/basic_structs.h | 14 +- src/utils/hydro_utilities.h | 8 +- src/utils/hydro_utilities_tests.cpp | 3 +- 12 files changed, 424 insertions(+), 131 deletions(-) create mode 100644 src/reconstruction/plmc_cuda.cu diff --git a/src/reconstruction/pcm_cuda.h b/src/reconstruction/pcm_cuda.h index d1c593ff2..be21cf950 100644 --- a/src/reconstruction/pcm_cuda.h +++ b/src/reconstruction/pcm_cuda.h @@ -63,11 +63,11 @@ reconstruction::InterfaceState __device__ __host__ inline PCM_Reconstruction(Rea interface_state.magnetic = conserved_data.magnetic; #endif // MHD #ifdef DE - interface_state.gas_energy_specific = primitive_data.gas_energy_specific; + interface_state.gas_energy = primitive_data.gas_energy; #endif // DE #ifdef SCALAR for (size_t i = 0; i < grid_enum::nscalars; i++) { - interface_state.scalar_specific[i] = primitive_data.scalar_specific[i]; + interface_state.scalar[i] = primitive_data.scalar[i]; } #endif // SCALAR diff --git a/src/reconstruction/plmc_cuda.cu b/src/reconstruction/plmc_cuda.cu new file mode 100644 index 000000000..0fb7a88fe --- /dev/null +++ b/src/reconstruction/plmc_cuda.cu @@ -0,0 +1,309 @@ +/*! \file plmc_cuda.cu + * \brief Definitions of the piecewise linear reconstruction functions with + limiting applied in the characteristic variables, as described + in Stone et al., 2008. */ + +#include + +#include "../global/global.h" +#include "../global/global_cuda.h" +#include "../reconstruction/plmc_cuda.h" +#include "../reconstruction/reconstruction_internals.h" +#include "../utils/cuda_utilities.h" +#include "../utils/gpu.hpp" + +#ifdef DE // PRESSURE_DE + #include "../utils/hydro_utilities.h" +#endif // DE + +/*! \fn __global__ 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) + * \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 n_fields) +{ + // get a thread ID + int const thread_id = threadIdx.x + blockIdx.x * blockDim.x; + int xid, yid, zid; + cuda_utilities::compute3DIndices(thread_id, nx, ny, xid, yid, zid); + + // Ensure that we are only operating on cells that will be used + if (reconstruction::Thread_Guard<2>(nx, ny, nz, xid, yid, zid)) { + return; + } + + // Compute the total number of cells + int const n_cells = nx * ny * nz; + + // Set the field indices for the various directions + int o1, o2, o3; + switch (dir) { + case 0: + o1 = grid_enum::momentum_x; + o2 = grid_enum::momentum_y; + o3 = grid_enum::momentum_z; + break; + case 1: + o1 = grid_enum::momentum_y; + o2 = grid_enum::momentum_z; + o3 = grid_enum::momentum_x; + break; + case 2: + o1 = grid_enum::momentum_z; + o2 = grid_enum::momentum_x; + o3 = grid_enum::momentum_y; + break; + } + + // load the 3-cell stencil into registers + // cell i + hydro_utilities::Primitive const cell_i = + 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 = 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 = 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); + Real const sound_speed_squared = sound_speed * sound_speed; + +// Compute the eigenvectors +#ifdef MHD + reconstruction::EigenVecs const eigenvectors = + reconstruction::Compute_Eigenvectors(cell_i, sound_speed, sound_speed_squared, gamma); +#else + reconstruction::EigenVecs eigenvectors; +#endif // MHD + + // Compute the left, right, centered, and van Leer differences of the + // primitive variables Note that here L and R refer to locations relative to + // the cell center + + // left + hydro_utilities::Primitive const del_L = reconstruction::Compute_Slope(cell_imo, cell_i); + + // right + hydro_utilities::Primitive const del_R = reconstruction::Compute_Slope(cell_i, cell_ipo); + + // centered + hydro_utilities::Primitive const del_C = reconstruction::Compute_Slope(cell_imo, cell_ipo, 0.5); + + // Van Leer + hydro_utilities::Primitive const del_G = reconstruction::Compute_Van_Leer_Slope(del_L, del_R); + + // Project the left, right, centered and van Leer differences onto the + // characteristic variables Stone Eqn 37 (del_a are differences in + // characteristic variables, see Stone for notation) Use the eigenvectors + // given in Stone 2008, Appendix A + reconstruction::Characteristic const del_a_L = + reconstruction::Primitive_To_Characteristic(cell_i, del_L, eigenvectors, sound_speed, sound_speed_squared, gamma); + + reconstruction::Characteristic const del_a_R = + reconstruction::Primitive_To_Characteristic(cell_i, del_R, eigenvectors, sound_speed, sound_speed_squared, gamma); + + reconstruction::Characteristic const del_a_C = + reconstruction::Primitive_To_Characteristic(cell_i, del_C, eigenvectors, sound_speed, sound_speed_squared, gamma); + + reconstruction::Characteristic const del_a_G = + reconstruction::Primitive_To_Characteristic(cell_i, del_G, eigenvectors, sound_speed, sound_speed_squared, gamma); + + // Apply monotonicity constraints to the differences in the characteristic variables and project the monotonized + // difference in the characteristic variables back onto the primitive variables Stone Eqn 39 + hydro_utilities::Primitive del_m_i = reconstruction::Monotonize_Characteristic_Return_Primitive( + cell_i, del_L, del_R, del_C, del_G, del_a_L, del_a_R, del_a_C, del_a_G, eigenvectors, sound_speed, + sound_speed_squared, gamma); + + // Compute the left and right interface values using the monotonized difference in the primitive variables + hydro_utilities::Primitive interface_L_iph = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, 1.0); + hydro_utilities::Primitive interface_R_imh = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, -1.0); + +#ifndef VL + + Real const dtodx = dt / dx; + + // Compute the eigenvalues of the linearized equations in the + // primitive variables using the cell-centered primitive variables + Real const lambda_m = cell_i.velocity.x() - sound_speed; + Real const lambda_0 = cell_i.velocity.x(); + Real const lambda_p = cell_i.velocity.x() + sound_speed; + + // Integrate linear interpolation function over domain of dependence + // defined by max(min) eigenvalue + Real qx = -0.5 * fmin(lambda_m, 0.0) * dtodx; + interface_R_imh.density = interface_R_imh.density + qx * del_m_i.density; + interface_R_imh.velocity.x() = interface_R_imh.velocity.x() + qx * del_m_i.velocity.x(); + interface_R_imh.velocity.y() = interface_R_imh.velocity.y() + qx * del_m_i.velocity.y(); + interface_R_imh.velocity.z() = interface_R_imh.velocity.z() + qx * del_m_i.velocity.z(); + interface_R_imh.pressure = interface_R_imh.pressure + qx * del_m_i.pressure; + + qx = 0.5 * fmax(lambda_p, 0.0) * dtodx; + interface_L_iph.density = interface_L_iph.density - qx * del_m_i.density; + interface_L_iph.velocity.x() = interface_L_iph.velocity.x() - qx * del_m_i.velocity.x(); + interface_L_iph.velocity.y() = interface_L_iph.velocity.y() - qx * del_m_i.velocity.y(); + interface_L_iph.velocity.z() = interface_L_iph.velocity.z() - qx * del_m_i.velocity.z(); + interface_L_iph.pressure = interface_L_iph.pressure - qx * del_m_i.pressure; + + #ifdef DE + interface_R_imh.gas_energy = interface_R_imh.gas_energy + qx * del_m_i.gas_energy; + interface_L_iph.gas_energy = interface_L_iph.gas_energy - qx * del_m_i.gas_energy; + #endif // DE + + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + interface_R_imh.scalar[i] = interface_R_imh.scalar[i] + qx * del_m_i.scalar[i]; + interface_L_iph.scalar[i] = interface_L_iph.scalar[i] - qx * del_m_i.scalar[i]; + } + #endif // SCALAR + + // Perform the characteristic tracing + // Stone Eqns 42 & 43 + + // left-hand interface value, i+1/2 + Real sum_0 = 0.0, sum_1 = 0.0, sum_2 = 0.0, sum_3 = 0.0, sum_4 = 0.0; + #ifdef DE + Real sum_ge = 0; + #endif // DE + #ifdef SCALAR + Real sum_scalar[NSCALARS]; + for (int i = 0; i < NSCALARS; i++) { + sum_scalar[i] = 0.0; + } + #endif // SCALAR + if (lambda_m >= 0) { + Real lamdiff = lambda_p - lambda_m; + + sum_0 += lamdiff * (-cell_i.density * del_m_i.velocity.x() / (2 * sound_speed) + + del_m_i.pressure / (2 * sound_speed_squared)); + sum_1 += lamdiff * (del_m_i.velocity.x() / 2.0 - del_m_i.pressure / (2 * sound_speed * cell_i.density)); + sum_4 += lamdiff * (-cell_i.density * del_m_i.velocity.x() * sound_speed / 2.0 + del_m_i.pressure / 2.0); + } + if (lambda_0 >= 0) { + Real lamdiff = lambda_p - lambda_0; + + sum_0 += lamdiff * (del_m_i.density - del_m_i.pressure / (sound_speed_squared)); + sum_2 += lamdiff * del_m_i.velocity.y(); + sum_3 += lamdiff * del_m_i.velocity.z(); + #ifdef DE + sum_ge += lamdiff * del_m_i.gas_energy; + #endif // DE + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + sum_scalar[i] += lamdiff * del_m_i.scalar[i]; + } + #endif // SCALAR + } + if (lambda_p >= 0) { + Real lamdiff = lambda_p - lambda_p; + + sum_0 += lamdiff * + (cell_i.density * del_m_i.velocity.x() / (2 * sound_speed) + del_m_i.pressure / (2 * sound_speed_squared)); + sum_1 += lamdiff * (del_m_i.velocity.x() / 2.0 + del_m_i.pressure / (2 * sound_speed * cell_i.density)); + sum_4 += lamdiff * (cell_i.density * del_m_i.velocity.x() * sound_speed / 2.0 + del_m_i.pressure / 2.0); + } + + // add the corrections to the initial guesses for the interface values + interface_L_iph.density += 0.5 * dtodx * sum_0; + interface_L_iph.velocity.x() += 0.5 * dtodx * sum_1; + interface_L_iph.velocity.y() += 0.5 * dtodx * sum_2; + interface_L_iph.velocity.z() += 0.5 * dtodx * sum_3; + interface_L_iph.pressure += 0.5 * dtodx * sum_4; + #ifdef DE + interface_L_iph.gas_energy += 0.5 * dtodx * sum_ge; + #endif // DE + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + interface_L_iph.scalar[i] += 0.5 * dtodx * sum_scalar[i]; + } + #endif // SCALAR + + // right-hand interface value, i-1/2 + sum_0 = sum_1 = sum_2 = sum_3 = sum_4 = 0; + #ifdef DE + sum_ge = 0; + #endif // DE + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + sum_scalar[i] = 0; + } + #endif // SCALAR + if (lambda_m <= 0) { + Real lamdiff = lambda_m - lambda_m; + + sum_0 += lamdiff * (-cell_i.density * del_m_i.velocity.x() / (2 * sound_speed) + + del_m_i.pressure / (2 * sound_speed_squared)); + sum_1 += lamdiff * (del_m_i.velocity.x() / 2.0 - del_m_i.pressure / (2 * sound_speed * cell_i.density)); + sum_4 += lamdiff * (-cell_i.density * del_m_i.velocity.x() * sound_speed / 2.0 + del_m_i.pressure / 2.0); + } + if (lambda_0 <= 0) { + Real lamdiff = lambda_m - lambda_0; + + sum_0 += lamdiff * (del_m_i.density - del_m_i.pressure / (sound_speed_squared)); + sum_2 += lamdiff * del_m_i.velocity.y(); + sum_3 += lamdiff * del_m_i.velocity.z(); + #ifdef DE + sum_ge += lamdiff * del_m_i.gas_energy; + #endif // DE + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + sum_scalar[i] += lamdiff * del_m_i.scalar[i]; + } + #endif // SCALAR + } + if (lambda_p <= 0) { + Real lamdiff = lambda_m - lambda_p; + + sum_0 += lamdiff * + (cell_i.density * del_m_i.velocity.x() / (2 * sound_speed) + del_m_i.pressure / (2 * sound_speed_squared)); + sum_1 += lamdiff * (del_m_i.velocity.x() / 2.0 + del_m_i.pressure / (2 * sound_speed * cell_i.density)); + sum_4 += lamdiff * (cell_i.density * del_m_i.velocity.x() * sound_speed / 2.0 + del_m_i.pressure / 2.0); + } + + // add the corrections + interface_R_imh.density += 0.5 * dtodx * sum_0; + interface_R_imh.velocity.x() += 0.5 * dtodx * sum_1; + interface_R_imh.velocity.y() += 0.5 * dtodx * sum_2; + interface_R_imh.velocity.z() += 0.5 * dtodx * sum_3; + interface_R_imh.pressure += 0.5 * dtodx * sum_4; + #ifdef DE + interface_R_imh.gas_energy += 0.5 * dtodx * sum_ge; + #endif // DE + #ifdef SCALAR + for (int i = 0; i < NSCALARS; i++) { + interface_R_imh.scalar[i] += 0.5 * dtodx * sum_scalar[i]; + } + #endif // SCALAR +#endif // CTU + + // apply minimum constraints + interface_R_imh.density = fmax(interface_R_imh.density, (Real)TINY_NUMBER); + interface_L_iph.density = fmax(interface_L_iph.density, (Real)TINY_NUMBER); + interface_R_imh.pressure = fmax(interface_R_imh.pressure, (Real)TINY_NUMBER); + interface_L_iph.pressure = fmax(interface_L_iph.pressure, (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 + size_t id = cuda_utilities::compute1DIndex(xid, yid, zid, nx, ny); + reconstruction::Write_Data(interface_L_iph, dev_bounds_L, dev_conserved, id, n_cells, o1, o2, o3, gamma); + + 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/ppmc_cuda.cu b/src/reconstruction/ppmc_cuda.cu index c13b78b91..4945ac69e 100644 --- a/src/reconstruction/ppmc_cuda.cu +++ b/src/reconstruction/ppmc_cuda.cu @@ -273,17 +273,15 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun Real const p_6 = 6.0 * (cell_i.pressure - 0.5 * (interface_R_imh.pressure + interface_L_iph.pressure)); #ifdef DE - del_m_i.gas_energy_specific = interface_L_iph.gas_energy_specific - interface_R_imh.gas_energy_specific; - Real const ge_6 = 6.0 * (cell_i.gas_energy_specific - - 0.5 * (interface_R_imh.gas_energy_specific + interface_L_iph.gas_energy_specific)); + del_m_i.gas_energy = interface_L_iph.gas_energy - interface_R_imh.gas_energy; + Real const ge_6 = 6.0 * (cell_i.gas_energy - 0.5 * (interface_R_imh.gas_energy + interface_L_iph.gas_energy)); #endif // DE #ifdef SCALAR Real scalar_6[NSCALARS]; for (int i = 0; i < NSCALARS; i++) { - del_m_i.scalar_specific[i] = interface_L_iph.scalar_specific[i] - interface_R_imh.scalar_specific[i]; - scalar_6[i] = 6.0 * (cell_i.scalar_specific[i] - - 0.5 * (interface_R_imh.scalar_specific[i] + interface_L_iph.scalar_specific[i])); + del_m_i.scalar[i] = interface_L_iph.scalar[i] - interface_R_imh.scalar[i]; + scalar_6[i] = 6.0 * (cell_i.scalar[i] - 0.5 * (interface_R_imh.scalar[i] + interface_L_iph.scalar[i])); } #endif // SCALAR @@ -342,24 +340,22 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun lambda_min * (0.5 * dtodx) * (del_m_i.pressure + (1.0 + (2.0 / 3.0) * lambda_min * dtodx) * p_6); #ifdef DE - interface_L_iph.gas_energy_specific = - interface_L_iph.gas_energy_specific - - lambda_max * (0.5 * dtodx) * (del_m_i.gas_energy_specific - (1.0 - (2.0 / 3.0) * lambda_max * dtodx) * ge_6); - interface_R_imh.gas_energy_specific = - interface_R_imh.gas_energy_specific - - lambda_min * (0.5 * dtodx) * (del_m_i.gas_energy_specific + (1.0 + (2.0 / 3.0) * lambda_min * dtodx) * ge_6); + interface_L_iph.gas_energy = + interface_L_iph.gas_energy - + lambda_max * (0.5 * dtodx) * (del_m_i.gas_energy - (1.0 - (2.0 / 3.0) * lambda_max * dtodx) * ge_6); + interface_R_imh.gas_energy = + interface_R_imh.gas_energy - + lambda_min * (0.5 * dtodx) * (del_m_i.gas_energy + (1.0 + (2.0 / 3.0) * lambda_min * dtodx) * ge_6); #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - interface_L_iph.scalar_specific[i] = - interface_L_iph.scalar_specific[i] - - lambda_max * (0.5 * dtodx) * - (del_m_i.scalar_specific[i] - (1.0 - (2.0 / 3.0) * lambda_max * dtodx) * scalar_6[i]); - interface_R_imh.scalar_specific[i] = - interface_R_imh.scalar_specific[i] - - lambda_min * (0.5 * dtodx) * - (del_m_i.scalar_specific[i] + (1.0 + (2.0 / 3.0) * lambda_min * dtodx) * scalar_6[i]); + interface_L_iph.scalar[i] = + interface_L_iph.scalar[i] - + lambda_max * (0.5 * dtodx) * (del_m_i.scalar[i] - (1.0 - (2.0 / 3.0) * lambda_max * dtodx) * scalar_6[i]); + interface_R_imh.scalar[i] = + interface_R_imh.scalar[i] - + lambda_min * (0.5 * dtodx) * (del_m_i.scalar[i] + (1.0 + (2.0 / 3.0) * lambda_min * dtodx) * scalar_6[i]); } #endif // SCALAR @@ -404,11 +400,11 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun Real const chi_4 = A * (del_m_i.velocity.z() - vz_6) + B * vz_6; Real const chi_5 = A * (del_m_i.pressure - p_6) + B * p_6; #ifdef DE - chi_ge = A * (del_m_i.gas_energy_specific - ge_6) + B * ge_6; + chi_ge = A * (del_m_i.gas_energy - ge_6) + B * ge_6; #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - chi_scalar[i] = A * (del_m_i.scalar_specific[i] - scalar_6[i]) + B * scalar_6[i]; + chi_scalar[i] = A * (del_m_i.scalar[i] - scalar_6[i]) + B * scalar_6[i]; } #endif // SCALAR @@ -446,11 +442,11 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun interface_L_iph.velocity.z() += sum_4; interface_L_iph.pressure += sum_5; #ifdef DE - interface_L_iph.gas_energy_specific += sum_ge; + interface_L_iph.gas_energy += sum_ge; #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - interface_L_iph.scalar_specific[i] += sum_scalar[i]; + interface_L_iph.scalar[i] += sum_scalar[i]; } #endif // SCALAR @@ -492,11 +488,11 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun Real const chi_4 = C * (del_m_i.velocity.z() + vz_6) + D * vz_6; Real const chi_5 = C * (del_m_i.pressure + p_6) + D * p_6; #ifdef DE - chi_ge = C * (del_m_i.gas_energy_specific + ge_6) + D * ge_6; + chi_ge = C * (del_m_i.gas_energy + ge_6) + D * ge_6; #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - chi_scalar[i] = C * (del_m_i.scalar_specific[i] + scalar_6[i]) + D * scalar_6[i]; + chi_scalar[i] = C * (del_m_i.scalar[i] + scalar_6[i]) + D * scalar_6[i]; } #endif // SCALAR @@ -534,11 +530,11 @@ __global__ void PPMC_CTU(Real *dev_conserved, Real *dev_bounds_L, Real *dev_boun interface_R_imh.velocity.z() += sum_4; interface_R_imh.pressure += sum_5; #ifdef DE - interface_R_imh.gas_energy_specific += sum_ge; + interface_R_imh.gas_energy += sum_ge; #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - interface_R_imh.scalar_specific[i] += sum_scalar[i]; + interface_R_imh.scalar[i] += sum_scalar[i]; } #endif // SCALAR @@ -680,17 +676,13 @@ __global__ __launch_bounds__(TPB) void PPMC_VL(Real *dev_conserved, Real *dev_bo // Compute the interfaces for the variables that don't have characteristics #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); + reconstruction::PPM_Single_Variable(cell_im2.gas_energy, cell_im1.gas_energy, cell_i.gas_energy, cell_ip1.gas_energy, + cell_ip2.gas_energy, interface_L_iph.gas_energy, interface_R_imh.gas_energy); #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]); + reconstruction::PPM_Single_Variable(cell_im2.scalar[i], cell_im1.scalar[i], cell_i.scalar[i], cell_ip1.scalar[i], + cell_ip2.scalar[i], interface_L_iph.scalar[i], interface_R_imh.scalar[i]); } #endif // SCALAR diff --git a/src/reconstruction/reconstruction_internals.h b/src/reconstruction/reconstruction_internals.h index dbb624e04..f8e041081 100644 --- a/src/reconstruction/reconstruction_internals.h +++ b/src/reconstruction/reconstruction_internals.h @@ -201,7 +201,7 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Load_Data( #endif // MHD loaded_data.pressure = hydro_utilities::Get_Pressure_From_DE(energy, energy - energy_non_thermal, gas_energy, gamma); - loaded_data.gas_energy_specific = gas_energy / loaded_data.density; + loaded_data.gas_energy = gas_energy / loaded_data.density; #else // not DE #ifdef MHD loaded_data.pressure = hydro_utilities::Calc_Pressure_Primitive( @@ -217,7 +217,7 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Load_Data( #ifdef SCALAR for (size_t i = 0; i < grid_enum::nscalars; i++) { - loaded_data.scalar_specific[i] = dev_conserved[(grid_enum::scalar + i) * n_cells + id] / loaded_data.density; + loaded_data.scalar[i] = dev_conserved[(grid_enum::scalar + i) * n_cells + id] / loaded_data.density; } #endif // SCALAR @@ -302,12 +302,12 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Compute_Slope(hydro_ut #endif // MHD #ifdef DE - slopes.gas_energy_specific = coef * (right.gas_energy_specific - left.gas_energy_specific); + slopes.gas_energy = coef * (right.gas_energy - left.gas_energy); #endif // DE #ifdef SCALAR for (size_t i = 0; i < grid_enum::nscalars; i++) { - slopes.scalar_specific[i] = coef * (right.scalar_specific[i] - left.scalar_specific[i]); + slopes.scalar[i] = coef * (right.scalar[i] - left.scalar[i]); } #endif // SCALAR @@ -348,12 +348,12 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Compute_Van_Leer_Slope #endif // MHD #ifdef DE - vl_slopes.gas_energy_specific = Calc_Vl_Slope(left_slope.gas_energy_specific, right_slope.gas_energy_specific); + vl_slopes.gas_energy = Calc_Vl_Slope(left_slope.gas_energy, right_slope.gas_energy); #endif // DE #ifdef SCALAR for (size_t i = 0; i < grid_enum::nscalars; i++) { - vl_slopes.scalar_specific[i] = Calc_Vl_Slope(left_slope.scalar_specific[i], right_slope.scalar_specific[i]); + vl_slopes.scalar[i] = Calc_Vl_Slope(left_slope.scalar[i], right_slope.scalar[i]); } #endif // SCALAR @@ -727,13 +727,13 @@ void __device__ __host__ __inline__ Monotonize_Parabolic_Interface(hydro_utiliti #endif // MHD #ifdef DE - Monotonize(cell_i.gas_energy_specific, cell_im1.gas_energy_specific, cell_ip1.gas_energy_specific, - interface_L_iph.gas_energy_specific, interface_R_imh.gas_energy_specific); + Monotonize(cell_i.gas_energy, cell_im1.gas_energy, cell_ip1.gas_energy, interface_L_iph.gas_energy, + interface_R_imh.gas_energy); #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - Monotonize(cell_i.scalar_specific[i], cell_im1.scalar_specific[i], cell_ip1.scalar_specific[i], - interface_L_iph.scalar_specific[i], interface_R_imh.scalar_specific[i]); + Monotonize(cell_i.scalar[i], cell_im1.scalar[i], cell_ip1.scalar[i], interface_L_iph.scalar[i], + interface_R_imh.scalar[i]); } #endif // SCALAR } @@ -767,11 +767,11 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Calc_Interface_Linear( #endif // MHD #ifdef DE - output.gas_energy_specific = interface(primitive.gas_energy_specific, slopes.gas_energy_specific); + output.gas_energy = interface(primitive.gas_energy, slopes.gas_energy); #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - output.scalar_specific[i] = interface(primitive.scalar_specific[i], slopes.scalar_specific[i]); + output.scalar[i] = interface(primitive.scalar[i], slopes.scalar[i]); } #endif // SCALAR @@ -817,13 +817,11 @@ hydro_utilities::Primitive __device__ __host__ __inline__ Calc_Interface_Parabol #endif // MHD #ifdef DE - output.gas_energy_specific = interface(cell_i.gas_energy_specific, cell_im1.gas_energy_specific, - slopes_i.gas_energy_specific, slopes_im1.gas_energy_specific); + output.gas_energy = interface(cell_i.gas_energy, cell_im1.gas_energy, slopes_i.gas_energy, slopes_im1.gas_energy); #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - output.scalar_specific[i] = interface(cell_i.scalar_specific[i], cell_im1.scalar_specific[i], - slopes_i.scalar_specific[i], slopes_im1.scalar_specific[i]); + output.scalar[i] = interface(cell_i.scalar[i], cell_im1.scalar[i], slopes_i.scalar[i], slopes_im1.scalar[i]); } #endif // SCALAR @@ -1021,12 +1019,11 @@ void __device__ __host__ __inline__ Write_Data(hydro_utilities::Primitive const #endif // MHD #ifdef DE - dev_interface[grid_enum::GasEnergy * n_cells + id] = interface_state.density * interface_state.gas_energy_specific; + dev_interface[grid_enum::GasEnergy * n_cells + id] = interface_state.density * interface_state.gas_energy; #endif // DE #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - dev_interface[(grid_enum::scalar + i) * n_cells + id] = - interface_state.density * interface_state.scalar_specific[i]; + dev_interface[(grid_enum::scalar + i) * n_cells + id] = interface_state.density * interface_state.scalar[i]; } #endif // SCALAR } diff --git a/src/riemann_solvers/exact_cuda.cu b/src/riemann_solvers/exact_cuda.cu index 736580ec5..0626e676a 100644 --- a/src/riemann_solvers/exact_cuda.cu +++ b/src/riemann_solvers/exact_cuda.cu @@ -80,11 +80,11 @@ __global__ void Calculate_Exact_Fluxes_CUDA(Real const *dev_conserved, Real cons left_state.pressure = fmax(left_state.pressure, (Real)TINY_NUMBER); #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - left_state.scalar_specific[i] = dev_bounds_L[(5 + i) * n_cells + tid] / left_state.density; + left_state.scalar[i] = dev_bounds_L[(5 + i) * n_cells + tid] / left_state.density; } #endif #ifdef DE - left_state.gas_energy_specific = dge / left_state.density; + left_state.gas_energy = dge / left_state.density; #endif right_state.density = dev_bounds_R[tid]; right_state.velocity.x() = dev_bounds_R[o1 * n_cells + tid] / right_state.density; @@ -109,11 +109,11 @@ __global__ void Calculate_Exact_Fluxes_CUDA(Real const *dev_conserved, Real cons right_state.pressure = fmax(right_state.pressure, (Real)TINY_NUMBER); #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - right_state.scalar_specific[i] = dev_bounds_R[(5 + i) * n_cells + tid] / right_state.density; + right_state.scalar[i] = dev_bounds_R[(5 + i) * n_cells + tid] / right_state.density; } #endif #ifdef DE - right_state.gas_energy_specific = dge / right_state.density; + right_state.gas_energy = dge / right_state.density; #endif } // compute sounds speeds in left and right regions @@ -144,11 +144,11 @@ __global__ void Calculate_Exact_Fluxes_CUDA(Real const *dev_conserved, Real cons dev_flux[o3 * n_cells + tid] = ds * vs * left_state.velocity.z(); #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - dev_flux[(5 + i) * n_cells + tid] = ds * vs * left_state.scalar_specific[i]; + dev_flux[(5 + i) * n_cells + tid] = ds * vs * left_state.scalar[i]; } #endif #ifdef DE - dev_flux[(n_fields - 1) * n_cells + tid] = ds * vs * left_state.gas_energy_specific; + dev_flux[(n_fields - 1) * n_cells + tid] = ds * vs * left_state.gas_energy; #endif Es = (ps / (gamma - 1.0)) + 0.5 * ds * (vs * vs + left_state.velocity.y() * left_state.velocity.y() + @@ -158,11 +158,11 @@ __global__ void Calculate_Exact_Fluxes_CUDA(Real const *dev_conserved, Real cons dev_flux[o3 * n_cells + tid] = ds * vs * right_state.velocity.z(); #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - dev_flux[(5 + i) * n_cells + tid] = ds * vs * right_state.scalar_specific[i]; + dev_flux[(5 + i) * n_cells + tid] = ds * vs * right_state.scalar[i]; } #endif #ifdef DE - dev_flux[(n_fields - 1) * n_cells + tid] = ds * vs * right_state.gas_energy_specific; + dev_flux[(n_fields - 1) * n_cells + tid] = ds * vs * right_state.gas_energy; #endif Es = (ps / (gamma - 1.0)) + 0.5 * ds * (vs * vs + right_state.velocity.y() * right_state.velocity.y() + diff --git a/src/riemann_solvers/hll_cuda.cu b/src/riemann_solvers/hll_cuda.cu index ad8adcb98..5073d8595 100644 --- a/src/riemann_solvers/hll_cuda.cu +++ b/src/riemann_solvers/hll_cuda.cu @@ -21,8 +21,6 @@ __global__ void Calculate_HLL_Fluxes_CUDA(Real const *dev_conserved, Real const int xid, yid, zid; cuda_utilities::compute3DIndices(tid, nx, ny, xid, yid, zid); - // Note that in only this riemann solver neither the scalar nor the dual energy that are carried through are NOT the - // specific versions despite the name in the struct reconstruction::InterfaceState left_state, right_state; Real f_d_l, f_mx_l, f_my_l, f_mz_l, f_E_l; @@ -71,11 +69,11 @@ __global__ void Calculate_HLL_Fluxes_CUDA(Real const *dev_conserved, Real const left_state.energy = dev_bounds_L[4 * n_cells + tid]; #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - left_state.scalar_specific[i] = dev_bounds_L[(5 + i) * n_cells + tid]; + left_state.scalar[i] = dev_bounds_L[(5 + i) * n_cells + tid]; } #endif #ifdef DE - left_state.gas_energy_specific = dev_bounds_L[(n_fields - 1) * n_cells + tid]; + left_state.gas_energy = dev_bounds_L[(n_fields - 1) * n_cells + tid]; #endif right_state.density = dev_bounds_R[tid]; @@ -85,11 +83,11 @@ __global__ void Calculate_HLL_Fluxes_CUDA(Real const *dev_conserved, Real const right_state.energy = dev_bounds_R[4 * n_cells + tid]; #ifdef SCALAR for (int i = 0; i < NSCALARS; i++) { - right_state.scalar_specific[i] = dev_bounds_R[(5 + i) * n_cells + tid]; + right_state.scalar[i] = dev_bounds_R[(5 + i) * n_cells + tid]; } #endif #ifdef DE - right_state.gas_energy_specific = dev_bounds_R[(n_fields - 1) * n_cells + tid]; + right_state.gas_energy = dev_bounds_R[(n_fields - 1) * n_cells + tid]; #endif // calculate primitive variables @@ -102,7 +100,7 @@ __global__ void Calculate_HLL_Fluxes_CUDA(Real const *dev_conserved, Real const (left_state.velocity.x() * left_state.velocity.x() + left_state.velocity.y() * left_state.velocity.y() + left_state.velocity.z() * left_state.velocity.z()); left_state.pressure = hydro_utilities::Get_Pressure_From_DE(left_state.energy, left_state.energy - E_kin, - left_state.gas_energy_specific, gamma); + left_state.gas_energy, gamma); #else left_state.pressure = (left_state.energy - 0.5 * left_state.density * (left_state.velocity.x() * left_state.velocity.x() + @@ -113,11 +111,11 @@ __global__ void Calculate_HLL_Fluxes_CUDA(Real const *dev_conserved, Real const left_state.pressure = fmax(left_state.pressure, (Real)TINY_NUMBER); // #ifdef SCALAR // for (int i=0; i= 0.0) { - dev_flux[(5 + i) * n_cells + tid] = dev_flux[tid] * left_state.scalar_specific[i]; + dev_flux[(5 + i) * n_cells + tid] = dev_flux[tid] * left_state.scalar[i]; } else { - dev_flux[(5 + i) * n_cells + tid] = dev_flux[tid] * right_state.scalar_specific[i]; + dev_flux[(5 + i) * n_cells + tid] = dev_flux[tid] * right_state.scalar[i]; } } #endif #ifdef DE if (dev_flux[tid] >= 0.0) { - dev_flux[(n_fields - 1) * n_cells + tid] = dev_flux[tid] * left_state.gas_energy_specific; + dev_flux[(n_fields - 1) * n_cells + tid] = dev_flux[tid] * left_state.gas_energy; } else { - dev_flux[(n_fields - 1) * n_cells + tid] = dev_flux[tid] * right_state.gas_energy_specific; + dev_flux[(n_fields - 1) * n_cells + tid] = dev_flux[tid] * right_state.gas_energy; } #endif } diff --git a/src/utils/basic_structs.h b/src/utils/basic_structs.h index bdef682f9..b3467b84f 100644 --- a/src/utils/basic_structs.h +++ b/src/utils/basic_structs.h @@ -123,11 +123,11 @@ struct Primitive { #ifdef DE /// The specific thermal energy in the gas - Real gas_energy_specific; + Real gas_energy; #endif // DE #ifdef SCALAR - Real scalar_specific[grid_enum::nscalars]; + Real scalar[grid_enum::nscalars]; #endif // SCALAR /// Default constructor, should init everything to zero @@ -135,7 +135,7 @@ struct Primitive { /// Manual constructor, mostly used for testing and doesn't init all members. The `in_` prefix stands for input, /// mostly to avoid name collision with the member variables Primitive(Real const in_density, VectorXYZ const& in_velocity, Real const in_pressure, - VectorXYZ const& in_magnetic = {0, 0, 0}, Real const in_gas_energy_specific = 0.0) + VectorXYZ const& in_magnetic = {0, 0, 0}, Real const in_gas_energy = 0.0) : density(in_density), velocity(in_velocity), pressure(in_pressure) { #ifdef MHD @@ -143,7 +143,7 @@ struct Primitive { #endif // mhd #ifdef DE - gas_energy_specific = in_gas_energy_specific; + gas_energy = in_gas_energy; #endif // DE }; }; @@ -166,11 +166,11 @@ struct InterfaceState { #endif // MHD #ifdef DE - Real gas_energy_specific; + Real gas_energy; #endif // DE #ifdef SCALAR - Real scalar_specific[grid_enum::nscalars]; + Real scalar[grid_enum::nscalars]; #endif // SCALAR // Define the constructors @@ -191,7 +191,7 @@ struct InterfaceState { total_pressure = in_total_pressure; #endif // MHD #ifdef DE - gas_energy_specific = 0.0; + gas_energy = 0.0; #endif // DE }; }; diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index ec9759d72..ddeecd608 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -308,12 +308,12 @@ __inline__ __host__ __device__ Primitive Conserved_2_Primitive(Conserved const & #endif // MHD #ifdef DE - output.gas_energy_specific = conserved_in.gas_energy / conserved_in.density; + output.gas_energy = 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; + output.scalar[i] = conserved_in.scalar[i] / conserved_in.density; } #endif // SCALAR @@ -368,12 +368,12 @@ __inline__ __host__ __device__ Conserved Primitive_2_Conserved(Primitive const & #endif // MHD #ifdef DE - output.gas_energy = primitive_in.gas_energy_specific * primitive_in.density; + output.gas_energy = primitive_in.gas_energy * 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; + output.scalar[i] = primitive_in.scalar[i] * primitive_in.density; } #endif // SCALAR diff --git a/src/utils/hydro_utilities_tests.cpp b/src/utils/hydro_utilities_tests.cpp index ea04d7549..d369255a9 100644 --- a/src/utils/hydro_utilities_tests.cpp +++ b/src/utils/hydro_utilities_tests.cpp @@ -425,8 +425,7 @@ TEST(tALLConserved2Primitive, CorrectInputExpectCorrectOutput) 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"); + testing_utilities::Check_Results(fiducial_data.gas_energy, test_data.gas_energy, "gas_energy"); #endif // DE }