diff --git a/cholla-tests-data b/cholla-tests-data index dcd73ff52..71eb66d63 160000 --- a/cholla-tests-data +++ b/cholla-tests-data @@ -1 +1 @@ -Subproject commit dcd73ff52b9027627b247c6d888bcdb56840c85e +Subproject commit 71eb66d63622ac15c0844ae96ec9034cd5e4f4d3 diff --git a/src/global/global.h b/src/global/global.h index 7c7b0dc13..f7030b563 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -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. diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu index d84cd7f00..33f739449 100644 --- a/src/hydro/hydro_cuda.cu +++ b/src/hydro/hydro_cuda.cu @@ -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 @@ -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; @@ -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 @@ -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; @@ -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 @@ -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; diff --git a/src/reconstruction/plmc_cuda.cu b/src/reconstruction/plmc_cuda.cu index 41a5ae505..bb31e9904 100644 --- a/src/reconstruction/plmc_cuda.cu +++ b/src/reconstruction/plmc_cuda.cu @@ -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; diff --git a/src/reconstruction/plmc_cuda_tests.cu b/src/reconstruction/plmc_cuda_tests.cu index 0c567b6d0..68f11b396 100644 --- a/src/reconstruction/plmc_cuda_tests.cu +++ b/src/reconstruction/plmc_cuda_tests.cu @@ -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> fiducial_interface_right = {{{20, 0.59023012197434721}, {84, 3.0043379408547275}, {148, 2.6320759184913625}, @@ -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 diff --git a/src/reconstruction/reconstruction.h b/src/reconstruction/reconstruction.h index 07aae21a6..23442a776 100644 --- a/src/reconstruction/reconstruction.h +++ b/src/reconstruction/reconstruction.h @@ -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 diff --git a/src/reconstruction/reconstruction_tests.cu b/src/reconstruction/reconstruction_tests.cu index 6c2e19af7..dc1f10720 100644 --- a/src/reconstruction/reconstruction_tests.cu +++ b/src/reconstruction/reconstruction_tests.cu @@ -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"); @@ -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"); +} diff --git a/src/system_tests/mhd_system_tests.cpp b/src/system_tests/mhd_system_tests.cpp index c7a21aaae..4261797b2 100644 --- a/src/system_tests/mhd_system_tests.cpp +++ b/src/system_tests/mhd_system_tests.cpp @@ -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 diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index c0f783e1c..24caff9f7 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -35,7 +35,7 @@ inline __host__ __device__ Real Calc_Pressure_Primitive(Real const &E, Real cons Real const &vz, Real const &gamma, Real const &magnetic_x = 0.0, Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0) { - Real pressure = (E - 0.5 * d * (vx * vx + ((vy * vy) + (vz * vz)))); + Real pressure = E - 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz); #ifdef MHD pressure -= mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z); @@ -48,7 +48,7 @@ inline __host__ __device__ Real Calc_Pressure_Conserved(Real const &E, Real cons Real const &mz, Real const &gamma, Real const &magnetic_x = 0.0, Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0) { - Real pressure = (E - 0.5 * (mx * mx + my * my + mz * mz) / d); + Real pressure = E - 0.5 * math_utils::SquareMagnitude(mx, my, mz) / d; #ifdef MHD pressure -= mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z); @@ -76,7 +76,7 @@ inline __host__ __device__ Real Calc_Energy_Primitive(Real const &P, Real const Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0) { // Compute and return energy - Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) + 0.5 * d * (vx * vx + vy * vy + vz * vz); + Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) + 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz); #ifdef MHD energy += mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z); @@ -92,7 +92,7 @@ inline __host__ __device__ Real Calc_Energy_Conserved(Real const &P, Real const { // Compute and return energy Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) + - (0.5 / d) * (momentum_x * momentum_x + momentum_y * momentum_y + momentum_z * momentum_z); + (0.5 / d) * math_utils::SquareMagnitude(momentum_x, momentum_y, momentum_z); #ifdef MHD energy += mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z); @@ -130,7 +130,7 @@ inline __host__ __device__ Real Get_Pressure_From_DE(Real const &E, Real const & inline __host__ __device__ Real Calc_Kinetic_Energy_From_Velocity(Real const &d, Real const &vx, Real const &vy, Real const &vz) { - return 0.5 * d * (vx * vx + vy * vy * vz * vz); + return 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz); } /*! @@ -145,7 +145,7 @@ inline __host__ __device__ Real Calc_Kinetic_Energy_From_Velocity(Real const &d, inline __host__ __device__ Real Calc_Kinetic_Energy_From_Momentum(Real const &d, Real const &mx, Real const &my, Real const &mz) { - return (0.5 / d) * (mx * mx + my * my * mz * mz); + return (0.5 / d) * math_utils::SquareMagnitude(mx, my, mz); } /*! diff --git a/src/utils/hydro_utilities_tests.cpp b/src/utils/hydro_utilities_tests.cpp index eda204c76..b200ddd8c 100644 --- a/src/utils/hydro_utilities_tests.cpp +++ b/src/utils/hydro_utilities_tests.cpp @@ -238,7 +238,7 @@ TEST(tHYDROHydroUtilsGetPressureFromDE, CorrectInputExpectCorrectOutput) TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput) { TestParams parameters; - std::vector fiducialEnergies{0.0, 6.307524975350106e-145, 7.3762470327090601e+249}; + std::vector fiducialEnergies{0.0, 6.307524975350106e-145, 1.9018677140549924e+150}; double const coef = 1E-50; for (size_t i = 0; i < parameters.names.size(); i++) { @@ -252,7 +252,7 @@ TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput) TEST(tHYDROtMHDCalcKineticEnergyFromMomentum, CorrectInputExpectCorrectOutput) { TestParams parameters; - std::vector fiducialEnergies{0.0, 0.0, 7.2568536478335773e+147}; + std::vector fiducialEnergies{0.0, 0.0, 3.0042157852278499e+49}; double const coef = 1E-50; for (size_t i = 0; i < parameters.names.size(); i++) { diff --git a/src/utils/math_utilities.h b/src/utils/math_utilities.h index 1480f852c..68d13f19d 100644 --- a/src/utils/math_utilities.h +++ b/src/utils/math_utilities.h @@ -82,4 +82,20 @@ inline __device__ __host__ Real dotProduct(Real const &a1, Real const &a2, Real }; // ========================================================================= +// ========================================================================= +/*! + * \brief Compute the magnitude of a vector + * + * \param[in] v1 The first element of the vector + * \param[in] v2 The second element of the vector + * \param[in] v3 The third element of the vector + * + * \return Real The dot product of a and b + */ +inline __device__ __host__ Real SquareMagnitude(Real const &v1, Real const &v2, Real const &v3) +{ + return dotProduct(v1, v2, v3, v1, v2, v3); +}; +// ========================================================================= + } // namespace math_utils diff --git a/src/utils/math_utilities_tests.cpp b/src/utils/math_utilities_tests.cpp index 83fd1d232..a49cd8a41 100644 --- a/src/utils/math_utilities_tests.cpp +++ b/src/utils/math_utilities_tests.cpp @@ -56,4 +56,22 @@ TEST(tALLDotProduct, CorrectInputExpectCorrectOutput) // Now check results testing_utilities::Check_Results(fiducialDotProduct, testDotProduct, "dot product"); } +// ========================================================================= + +// ========================================================================= +/*! + * \brief Test the math_utils::dotProduct function + * + */ +TEST(tALLSquareMagnitude, CorrectInputExpectCorrectOutput) +{ + std::vector a = {11.503067766457753, 98.316634031589935, 41.12177317622657}; + + double const fiducial_square_magnitude = 11489.481324498336; + + double test_square_magnitude = math_utils::SquareMagnitude(a.at(0), a.at(1), a.at(2)); + + // Now check results + testing_utilities::Check_Results(fiducial_square_magnitude, test_square_magnitude, "dot product"); +} // ========================================================================= \ No newline at end of file diff --git a/src/utils/mhd_utilities.h b/src/utils/mhd_utilities.h index 55ecc6f75..f409fd4b0 100644 --- a/src/utils/mhd_utilities.h +++ b/src/utils/mhd_utilities.h @@ -18,6 +18,7 @@ #include "../grid/grid3D.h" #include "../utils/cuda_utilities.h" #include "../utils/gpu.hpp" +#include "../utils/math_utilities.h" namespace mhd::utils { @@ -74,7 +75,7 @@ inline __host__ __device__ Real _magnetosonicSpeed(Real const &density, Real con inline __host__ __device__ Real computeMagneticEnergy(Real const &magneticX, Real const &magneticY, Real const &magneticZ) { - return 0.5 * (magneticX * magneticX + ((magneticY * magneticY) + (magneticZ * magneticZ))); + return 0.5 * math_utils::SquareMagnitude(magneticX, magneticY, magneticZ); } // ========================================================================= @@ -98,9 +99,7 @@ inline __host__ __device__ Real computeThermalEnergy(Real const &energyTot, Real Real const &magneticX, Real const &magneticY, Real const &magneticZ, Real const &gamma) { - return energyTot - - 0.5 * (momentumX * momentumX + ((momentumY * momentumY) + (momentumZ * momentumZ))) / - fmax(density, TINY_NUMBER) - + return energyTot - 0.5 * math_utils::SquareMagnitude(momentumX, momentumY, momentumZ) / fmax(density, TINY_NUMBER) - computeMagneticEnergy(magneticX, magneticY, magneticZ); } // =========================================================================