Skip to content

Commit

Permalink
Fix bug in calculation of kinetic energy
Browse files Browse the repository at this point in the history 
Two of the terms were multiplied instead of added. Added a new function
for computing the square of the magnitude
  • Loading branch information
@bcaddy
bcaddy committed Jan 12, 2024
1 parent 4c548a7 commit c38c115
Show file tree
Hide file tree
Showing 8 changed files with 35 additions and 20 deletions.
2 changes: 1 addition & 1 deletion cholla-tests-data
Submodule cholla-tests-data updated 1 files
+1 −1 system_tests/tMHDSYSTEMParameterizedMpi_RyuAndJones1aShockTubeCorrectInputExpectCorrectOutput.h5
10 changes: 5 additions & 5 deletions src/reconstruction/plmc_cuda_tests.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,18 +226,18 @@ 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},
},
{
{5, 0.92705119413602599},
{69, 1.9592598982258778},
{133, 0.96653490574340428},
{197, 1.3203867992383289},
{261, 7.6371723945376493},
{325, 1.7033818819502551},
{261, 7.9217487636977353},
{325, 1.8629714367312684},
{389, 1.8587936590169301},
}};

Expand Down
2 changes: 1 addition & 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
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
12 changes: 6 additions & 6 deletions src/utils/hydro_utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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);
}

/*!
Expand All @@ -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);
}

/*!
Expand Down
4 changes: 2 additions & 2 deletions src/utils/hydro_utilities_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@ TEST(tHYDROHydroUtilsGetPressureFromDE, CorrectInputExpectCorrectOutput)
TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput)
{
TestParams parameters;
std::vector<double> fiducialEnergies{0.0, 6.307524975350106e-145, 7.3762470327090601e+249};
std::vector<double> fiducialEnergies{0.0, 6.307524975350106e-145, 1.9018677140549924e+150};
double const coef = 1E-50;

for (size_t i = 0; i < parameters.names.size(); i++) {
Expand All @@ -252,7 +252,7 @@ TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput)
TEST(tHYDROtMHDCalcKineticEnergyFromMomentum, CorrectInputExpectCorrectOutput)
{
TestParams parameters;
std::vector<double> fiducialEnergies{0.0, 0.0, 7.2568536478335773e+147};
std::vector<double> fiducialEnergies{0.0, 0.0, 3.0042157852278499e+49};
double const coef = 1E-50;

for (size_t i = 0; i < parameters.names.size(); i++) {
Expand Down
16 changes: 16 additions & 0 deletions src/utils/math_utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
7 changes: 3 additions & 4 deletions src/utils/mhd_utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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);
}
// =========================================================================

Expand All @@ -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);
}
// =========================================================================
Expand Down

0 comments on commit c38c115

Please sign in to comment.