From 35e1b8dec99820f0c47a9bfaa5e6219295bb857e Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Fri, 12 Jan 2024 14:15:24 -0500 Subject: [PATCH] Fix bug in calculation of kinetic energy Two of the terms were multiplied instead of added. Added a new function for computing the square of the magnitude --- cholla-tests-data | 2 +- src/reconstruction/plmc_cuda_tests.cu | 10 +++++----- src/reconstruction/reconstruction_tests.cu | 2 +- src/system_tests/mhd_system_tests.cpp | 2 +- src/utils/hydro_utilities.h | 12 ++++++------ src/utils/hydro_utilities_tests.cpp | 4 ++-- src/utils/math_utilities.h | 16 ++++++++++++++++ src/utils/mhd_utilities.h | 7 +++---- 8 files changed, 35 insertions(+), 20 deletions(-) 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/reconstruction/plmc_cuda_tests.cu b/src/reconstruction/plmc_cuda_tests.cu index 429bb7a89..68f11b396 100644 --- a/src/reconstruction/plmc_cuda_tests.cu +++ b/src/reconstruction/plmc_cuda_tests.cu @@ -226,9 +226,9 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput) {17, 0.44405384992296193}, {81, 2.5027813113931279}, {145, 2.6371119205792346}, - {209, 0.71381042558869023}, - {273, 20.742152413015724}, - {337, 2.1583975283044725}, + {209, 1.0210845222961809}, + {273, 21.353253570231175}, + {337, 2.1634182515826184}, {401, 1.7033818819502551}, }, { @@ -236,8 +236,8 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput) {69, 1.9592598982258778}, {133, 0.96653490574340428}, {197, 1.3203867992383289}, - {261, 7.6371723945376493}, - {325, 1.7033818819502551}, + {261, 7.9217487636977353}, + {325, 1.8629714367312684}, {389, 1.8587936590169301}, }}; diff --git a/src/reconstruction/reconstruction_tests.cu b/src/reconstruction/reconstruction_tests.cu index 008ce6752..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"); diff --git a/src/system_tests/mhd_system_tests.cpp b/src/system_tests/mhd_system_tests.cpp index c7a21aaae..89be3059d 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, 1.8E-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/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); } // =========================================================================