Skip to content

Commit

Permalink
Merge branch 'dev' into dev-iss241
Browse files Browse the repository at this point in the history
  • Loading branch information
@evaneschneider
evaneschneider authored Jan 18, 2024
2 parents 5fe262a + 2bd8222 commit e429d3f
Show file tree
Hide file tree
Showing 13 changed files with 217 additions and 33 deletions.
2 changes: 1 addition & 1 deletion cholla-tests-data
Submodule cholla-tests-data updated 1 files
+1 −1 system_tests/tMHDSYSTEMParameterizedMpi_RyuAndJones1aShockTubeCorrectInputExpectCorrectOutput.h5
23 changes: 19 additions & 4 deletions src/global/global.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,25 @@ typedef double Real;
#define TEMP_FLOOR 1e-3
#define DENS_FLOOR 1e-5 // in code units

// Parameter for Enzo dual Energy Condition
#define DE_ETA_1 \
0.001 // Ratio of U to E for which Internal Energy is used to compute the
// Pressure
// Parameters for Enzo dual Energy Condition
// - Prior to GH PR #356, DE_ETA_1 nominally had a value of 0.001 in all
// simulations (in practice, the value of DE_ETA_1 had minimal significance
// in those simulations). In PR #356, we revised the internal-energy
// synchronization to account for the value of DE_ETA_1. This was necessary
// for non-cosmology simulations.
// - In Cosmological simulation, we set DE_ETA_1 to a large number (it doesn't
// really matter what, as long as its >=1) to maintain the older behavior
// - In the future, we run tests and revisit the choice of DE_ETA_1 in
// cosmological simulations
#ifdef COSMOLOGY
#define DE_ETA_1 10.0
#else
#define DE_ETA_1 \
0.001 // Ratio of U to E for which Internal Energy is used to compute the
// Pressure. This also affects when the Internal Energy is used for
// the update.
#endif

#define DE_ETA_2 \
0.035 // Ratio of U to max(E_local) used to select which Internal Energy is
// used for the update.
Expand Down
18 changes: 15 additions & 3 deletions src/hydro/hydro_cuda.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -839,6 +839,7 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho
int imo, ipo;
n_cells = nx;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand All @@ -864,7 +865,10 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho
Emax = fmax(dev_conserved[4 * n_cells + imo], E);
Emax = fmax(Emax, dev_conserved[4 * n_cells + ipo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand All @@ -887,6 +891,7 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i
int imo, ipo, jmo, jpo;
n_cells = nx * ny;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand Down Expand Up @@ -922,7 +927,10 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i
Emax = fmax(Emax, dev_conserved[4 * n_cells + jmo]);
Emax = fmax(Emax, dev_conserved[4 * n_cells + jpo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand All @@ -945,6 +953,7 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i
int imo, ipo, jmo, jpo, kmo, kpo;
n_cells = nx * ny * nz;

Real eta_1 = DE_ETA_1;
Real eta_2 = DE_ETA_2;

// get a global thread ID
Expand Down Expand Up @@ -987,7 +996,10 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i
Emax = fmax(Emax, dev_conserved[4 * n_cells + kmo]);
Emax = fmax(Emax, dev_conserved[4 * n_cells + kpo]);

if (U_total / Emax > eta_2) {
// We only use the "advected" internal energy if both:
// - the thermal energy divided by total energy is a small fraction (smaller than eta_1)
// - AND we aren't masking shock heating (details controlled by Emax & eta_2)
if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
U = U_total;
} else {
U = U_advected;
Expand Down
3 changes: 3 additions & 0 deletions src/reconstruction/plmc_cuda.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,9 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_
reconstruction::Primitive interface_L_iph = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, 1.0);
reconstruction::Primitive interface_R_imh = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, -1.0);

// Limit the interfaces
reconstruction::Plm_Limit_Interfaces(interface_L_iph, interface_R_imh, cell_imo, cell_i, cell_ipo);

#ifndef VL

Real const dtodx = dt / dx;
Expand Down
22 changes: 11 additions & 11 deletions src/reconstruction/plmc_cuda_tests.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,17 +204,17 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
{{21, 0.73640639402573249},
{85, 3.3462413154443715},
{149, 2.1945584994458125},
{213, 0.67418839414138987},
{277, 16.909618487528142},
{341, 2.1533768050263267},
{405, 1.6994195863331925}},
{213, 1.1837630990406585},
{277, 17.570011907061254},
{341, 2.1583975283044725},
{405, 1.7033818819502551}},
{{21, 0.25340904981266843},
{85, 2.0441984720128734},
{149, 1.9959059157695584},
{213, 0.45377591914009824},
{277, 23.677832869261188},
{341, 1.5437923271692418},
{405, 1.8141353672443383}}};
{277, 24.018953780483471},
{341, 1.7033818819502551},
{405, 1.8587936590169301}}};
std::vector<std::unordered_map<int, double>> fiducial_interface_right = {{{20, 0.59023012197434721},
{84, 3.0043379408547275},
{148, 2.6320759184913625},
Expand All @@ -227,18 +227,18 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
{81, 2.5027813113931279},
{145, 2.6371119205792346},
{209, 1.0210845222961809},
{273, 21.360010722689488},
{273, 21.353253570231175},
{337, 2.1634182515826184},
{401, 1.7073441775673177},
{401, 1.7033818819502551},
},
{
{5, 0.92705119413602599},
{69, 1.9592598982258778},
{133, 0.96653490574340428},
{197, 1.3203867992383289},
{261, 8.0057564947791793},
{261, 7.9217487636977353},
{325, 1.8629714367312684},
{389, 1.9034519507895218},
{389, 1.8587936590169301},
}};

// Loop over different directions
Expand Down
50 changes: 50 additions & 0 deletions src/reconstruction/reconstruction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -673,6 +673,56 @@ Primitive __device__ __host__ __inline__ Calc_Interface_Linear(Primitive const &
}
// =====================================================================================================================

// =====================================================================================================================
/*!
* \brief Apply limiting the the primitive interfaces in PLM reconstructions
*
* \param[in,out] interface_L_iph The unlimited left plus 1/2 interface
* \param[in,out] interface_R_imh The unlimited right minus 1/2 interface
* \param[in] cell_imo The cell centered values at i-1
* \param[in] cell_i The cell centered values at i
* \param[in] cell_ipo The cell centered values at i+1
*/
void __device__ __host__ __inline__ Plm_Limit_Interfaces(Primitive &interface_L_iph, Primitive &interface_R_imh,
Primitive const &cell_imo, Primitive const &cell_i,
Primitive const &cell_ipo)
{
auto limiter = [](Real &l_iph, Real &r_imh, Real const &val_imo, Real const &val_i, Real const &val_ipo) {
r_imh = fmax(fmin(val_i, val_imo), r_imh);
r_imh = fmin(fmax(val_i, val_imo), r_imh);
l_iph = fmax(fmin(val_i, val_ipo), l_iph);
l_iph = fmin(fmax(val_i, val_ipo), l_iph);
};

limiter(interface_L_iph.density, interface_R_imh.density, cell_imo.density, cell_i.density, cell_ipo.density);
limiter(interface_L_iph.velocity_x, interface_R_imh.velocity_x, cell_imo.velocity_x, cell_i.velocity_x,
cell_ipo.velocity_x);
limiter(interface_L_iph.velocity_y, interface_R_imh.velocity_y, cell_imo.velocity_y, cell_i.velocity_y,
cell_ipo.velocity_y);
limiter(interface_L_iph.velocity_z, interface_R_imh.velocity_z, cell_imo.velocity_z, cell_i.velocity_z,
cell_ipo.velocity_z);
limiter(interface_L_iph.pressure, interface_R_imh.pressure, cell_imo.pressure, cell_i.pressure, cell_ipo.pressure);

#ifdef MHD
limiter(interface_L_iph.magnetic_y, interface_R_imh.magnetic_y, cell_imo.magnetic_y, cell_i.magnetic_y,
cell_ipo.magnetic_y);
limiter(interface_L_iph.magnetic_z, interface_R_imh.magnetic_z, cell_imo.magnetic_z, cell_i.magnetic_z,
cell_ipo.magnetic_z);
#endif // MHD

#ifdef DE
limiter(interface_L_iph.gas_energy, interface_R_imh.gas_energy, cell_imo.gas_energy, cell_i.gas_energy,
cell_ipo.gas_energy);
#endif // DE
#ifdef SCALAR
for (int i = 0; i < NSCALARS; i++) {
limiter(interface_L_iph.scalar[i], interface_R_imh.scalar[i], cell_imo.scalar[i], cell_i.scalar[i],
cell_ipo.scalar[i]);
}
#endif // SCALAR
}
// =====================================================================================================================

// =====================================================================================================================
/*!
* \brief Compute the interface state for the CTU version fo the reconstructor from the slope and cell centered state
Expand Down
73 changes: 72 additions & 1 deletion src/reconstruction/reconstruction_tests.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@ TEST(tALLReconstructionLoadData, CorrectInputExpectCorrectOutput)
reconstruction::Primitive fiducial_data{13, 3.0769230769230771, 5.1538461538461542, 7.2307692307692308,
39950.641025641031};
#ifdef DE
fiducial_data.pressure = 34274.282506448195;
fiducial_data.pressure = 39950.641025641031;
#endif // DE
testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density");
testing_utilities::Check_Results(fiducial_data.velocity_x, test_data.velocity_x, "velocity_x");
Expand Down Expand Up @@ -609,3 +609,74 @@ TEST(tALLReconstructionWriteData, CorrectInputExpectCorrectOutput)
testing_utilities::Check_Results(fiducial_val, test_val, "Interface at i=" + std::to_string(i));
}
}

TEST(tHYDROReconstructionPlmLimitInterfaces, CorrectInputExpectCorrectOutput)
{
// Set up values to test
reconstruction::Primitive interface_l_iph, interface_r_imh;
reconstruction::Primitive cell_im1, cell_i, cell_ip1;
interface_r_imh.density = -1.94432878387898625e+14;
interface_r_imh.velocity_x = 1.42049955114756404e-04;
interface_r_imh.velocity_y = -2.61311412306644180e-06;
interface_r_imh.velocity_z = -1.99429361865204601e-07;
interface_r_imh.pressure = -2.01130121665840250e-14;
interface_l_iph.density = 1.94433200621991188e+14;
interface_l_iph.velocity_x = 1.42025407335853601e-04;
interface_l_iph.velocity_y = -2.61311412306644180e-06;
interface_l_iph.velocity_z = -6.01154878659959398e-06;
interface_l_iph.pressure = 2.01130321665840277e-14;

cell_im1.density = 1.61101072114153951e+08;
cell_i.density = 1.61117046279133737e+08;
cell_ip1.density = 1.61011252191243321e+08;
cell_im1.velocity_x = 1.42067642369120116e-04;
cell_i.velocity_x = 1.42037681225305003e-04;
cell_ip1.velocity_x = 1.41901817571928041e-04;
cell_im1.velocity_y = -2.61228250783092252e-06;
cell_i.velocity_y = -2.61311412306644180e-06;
cell_ip1.velocity_y = -2.61155204131260820e-06;
cell_im1.velocity_z = 2.71420653365757378e-06;
cell_i.velocity_z = -3.10548907423239929e-06;
cell_ip1.velocity_z = -8.91005201578514336e-06;
cell_im1.pressure = 9.99999999999999945e-21;
cell_i.pressure = 9.99999999999999945e-21;
cell_ip1.pressure = 4.70262856027679407e-03;

// Set fiducial values
reconstruction::Primitive interface_r_imh_fiducial, interface_l_iph_fiducial;
interface_r_imh_fiducial.density = 161101072.11415395;
interface_r_imh_fiducial.velocity_x = 1.42049955114756404e-04;
interface_r_imh_fiducial.velocity_y = -2.61311412306644180e-06;
interface_r_imh_fiducial.velocity_z = -1.99429361865204601e-07;
interface_r_imh_fiducial.pressure = 9.99999999999999945e-21;
interface_l_iph_fiducial.density = 1.61117046279133737e+08;
interface_l_iph_fiducial.velocity_x = 1.42025407335853601e-04;
interface_l_iph_fiducial.velocity_y = -2.61311412306644180e-06;
interface_l_iph_fiducial.velocity_z = -6.01154878659959398e-06;
interface_l_iph_fiducial.pressure = 2.0113032166584028e-14;

// Run function
reconstruction::Plm_Limit_Interfaces(interface_l_iph, interface_r_imh, cell_im1, cell_i, cell_ip1);

// Check values
testing_utilities::Check_Results(interface_l_iph_fiducial.density, interface_l_iph.density,
"Mismatch in l_iph density");
testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_x, interface_l_iph.velocity_x,
"Mismatch in l_iph velocity_x");
testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_y, interface_l_iph.velocity_y,
"Mismatch in l_iph velocity_y");
testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_z, interface_l_iph.velocity_z,
"Mismatch in l_iph velocity_z");
testing_utilities::Check_Results(interface_l_iph_fiducial.pressure, interface_l_iph.pressure,
"Mismatch in l_iph pressure");
testing_utilities::Check_Results(interface_r_imh_fiducial.density, interface_r_imh.density,
"Mismatch in r_imh density");
testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_x, interface_r_imh.velocity_x,
"Mismatch in r_imh velocity_x");
testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_y, interface_r_imh.velocity_y,
"Mismatch in r_imh velocity_y");
testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_z, interface_r_imh.velocity_z,
"Mismatch in r_imh velocity_z");
testing_utilities::Check_Results(interface_r_imh_fiducial.pressure, interface_r_imh.pressure,
"Mismatch in r_imh pressure");
}
2 changes: 1 addition & 1 deletion src/system_tests/mhd_system_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -765,7 +765,7 @@ TEST_P(tMHDSYSTEMParameterizedMpi, AdvectingFieldLoopCorrectInputExpectCorrectOu
test_runner.numMpiRanks = GetParam();

// Only do the L2 Norm test. The regular cell-to-cell comparison is brittle for this test across systems
test_runner.runTest(true, 3.9E-8, 1.6E-6);
test_runner.runTest(true, 3.9E-8, 2.25E-6);
}

/// Test the MHD Blast Wave
Expand Down
Loading

0 comments on commit e429d3f

Please sign in to comment.