From 6c4171a6f4608cca4fa81cdd07f91b8ec3ae78a3 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 21 May 2024 15:42:15 -0400 Subject: [PATCH] Promote DSMC calculations to double precision to allow use with single-precision particles (#4941) * promote DSMC calculations to double precision * remove commented-out check * double -> auto for clang-tidy --- .../BinaryCollision/BinaryCollisionUtils.H | 62 ++++++++++--------- .../BinaryCollision/DSMC/DSMCFunc.cpp | 3 - 2 files changed, 34 insertions(+), 31 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H index 08ddcd27174..2e24a93a35e 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H @@ -62,60 +62,66 @@ namespace BinaryCollisionUtils{ using namespace amrex::literals; using namespace amrex::Math; - constexpr auto one_pr = amrex::ParticleReal(1.); - constexpr auto inv_four_pr = amrex::ParticleReal(1./4.); constexpr double c_sq = PhysConst::c * PhysConst::c; constexpr double inv_csq = 1.0 / c_sq; - const amrex::ParticleReal m1_sq = m1*m1; - const amrex::ParticleReal m2_sq = m2*m2; + // Cast input parameters to double before computing collision properties + // This is needed to avoid errors when using single-precision particles + const auto m1_dbl = static_cast(m1); + const auto m2_dbl = static_cast(m2); + const auto u1x_dbl = static_cast(u1x); + const auto u1y_dbl = static_cast(u1y); + const auto u1z_dbl = static_cast(u1z); + const auto u2x_dbl = static_cast(u2x); + const auto u2y_dbl = static_cast(u2y); + const auto u2z_dbl = static_cast(u2z); + + const double m1_sq = m1_dbl*m1_dbl; + const double m2_sq = m2_dbl*m2_dbl; // Compute Lorentz factor gamma in the lab frame - const double g1 = std::sqrt( 1.0 + static_cast(u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq ); - const double g2 = std::sqrt( 1.0 + static_cast(u2x*u2x+u2y*u2y+u2z*u2z)*inv_csq ); + const double g1 = std::sqrt( 1.0 + (u1x_dbl*u1x_dbl + u1y_dbl*u1y_dbl + u1z_dbl*u1z_dbl)*inv_csq ); + const double g2 = std::sqrt( 1.0 + (u2x_dbl*u2x_dbl + u2y_dbl*u2y_dbl + u2z_dbl*u2z_dbl)*inv_csq ); // Compute momenta - const amrex::ParticleReal p1x = u1x * m1; - const amrex::ParticleReal p1y = u1y * m1; - const amrex::ParticleReal p1z = u1z * m1; - const amrex::ParticleReal p2x = u2x * m2; - const amrex::ParticleReal p2y = u2y * m2; - const amrex::ParticleReal p2z = u2z * m2; + const double p1x = u1x_dbl * m1_dbl; + const double p1y = u1y_dbl * m1_dbl; + const double p1z = u1z_dbl * m1_dbl; + const double p2x = u2x_dbl * m2_dbl; + const double p2y = u2y_dbl * m2_dbl; + const double p2z = u2z_dbl * m2_dbl; + // Square norm of the total (sum between the two particles) momenta in the lab frame - const auto p_total_sq = static_cast( - powi<2>(p1x + p2x) + powi<2>(p1y + p2y) + powi<2>(p1z + p2z) - ); + const double p_total_sq = powi<2>(p1x + p2x) + powi<2>(p1y + p2y) + powi<2>(p1z + p2z); // Total energy in the lab frame // Note the use of `double` for energy since this calculation is // prone to error with single precision. - const auto m1_dbl = static_cast(m1); - const auto m2_dbl = static_cast(m2); - const double E_lab = (m1_dbl * g1 + m2_dbl * g2) * c_sq; + const double E_lab = (m1_dbl*g1 + m2_dbl*g2) * c_sq; // Total energy squared in the center of mass frame, calculated using the Lorentz invariance // of the four-momentum norm const double E_star_sq = E_lab*E_lab - c_sq*p_total_sq; // Kinetic energy in the center of mass frame const double E_star = std::sqrt(E_star_sq); + + // Cast back to chosen precision for output E_kin_COM = static_cast(E_star - (m1_dbl + m2_dbl)*c_sq); // Square of the norm of the momentum of one of the particles in the center of mass frame // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle // The expression below is specifically written in a form that avoids returning // small negative numbers due to machine precision errors, for low-energy particles - const auto E_ratio = static_cast(E_star/((m1 + m2)*c_sq)); - const auto p_star_sq = static_cast( - m1*m2*c_sq * ( powi<2>(E_ratio) - one_pr ) - + powi<2>(m1 - m2)*c_sq*inv_four_pr * powi<2>( E_ratio - 1._prt/E_ratio) - ); + const double E_ratio = E_star/((m1_dbl + m2_dbl)*c_sq); + const double p_star_sq = m1_dbl*m2_dbl*c_sq * ( powi<2>(E_ratio) - 1.0 ) + + powi<2>(m1_dbl - m2_dbl)*c_sq/4.0 * powi<2>( E_ratio - 1.0/E_ratio); // Lorentz factors in the center of mass frame - const auto g1_star = std::sqrt(one_pr + p_star_sq / static_cast(m1_sq*c_sq)); - const auto g2_star = std::sqrt(one_pr + p_star_sq / static_cast(m2_sq*c_sq)); + const double g1_star = std::sqrt(1.0 + p_star_sq / (m1_sq*c_sq)); + const double g2_star = std::sqrt(1.0 + p_star_sq / (m2_sq*c_sq)); - // relative velocity in the center of mass frame - v_rel_COM = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + one_pr/(m2*g2_star)); + // relative velocity in the center of mass frame, cast back to chosen precision + v_rel_COM = static_cast(std::sqrt(p_star_sq) * (1.0/(m1_dbl*g1_star) + 1.0/(m2_dbl*g2_star))); // Cross sections and relative velocity are computed in the center of mass frame. // On the other hand, the particle densities (weight over volume) in the lab frame are used. @@ -124,7 +130,7 @@ namespace BinaryCollisionUtils{ // COM frame and the Lorentz factors in the lab frame (see // Perez et al., Phys.Plasmas.19.083104 (2012)). The correction factor // is calculated here. - lab_to_COM_lorentz_factor = g1_star*g2_star/static_cast(g1*g2); + lab_to_COM_lorentz_factor = static_cast(g1_star*g2_star/(g1*g2)); } /** diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp index ca0818a6041..9ee676bf002 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -23,9 +23,6 @@ DSMCFunc::DSMCFunc ( { using namespace amrex::literals; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (std::is_same::value), - "Particle precision must be double for DSMC collisions."); - const amrex::ParmParse pp_collision_name(collision_name); // query for a list of collision processes