Promote DSMC calculations to double precision to allow use with singl…

…e-precision particles (ECP-WarpX#4941)

* promote DSMC calculations to double precision

* remove commented-out check

* double -> auto for clang-tidy

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<double>(m1);
+        const auto m2_dbl = static_cast<double>(m2);
+        const auto u1x_dbl = static_cast<double>(u1x);
+        const auto u1y_dbl = static_cast<double>(u1y);
+        const auto u1z_dbl = static_cast<double>(u1z);
+        const auto u2x_dbl = static_cast<double>(u2x);
+        const auto u2y_dbl = static_cast<double>(u2y);
+        const auto u2z_dbl = static_cast<double>(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<double>(u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq );
-        const double g2 = std::sqrt( 1.0 + static_cast<double>(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<double>(
-            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<double>(m1);
-        const auto m2_dbl = static_cast<double>(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<amrex::ParticleReal>(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<amrex::ParticleReal>(E_star/((m1 + m2)*c_sq));
-        const auto p_star_sq = static_cast<amrex::ParticleReal>(
-            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<amrex::ParticleReal>(m1_sq*c_sq));
-        const auto g2_star = std::sqrt(one_pr + p_star_sq / static_cast<amrex::ParticleReal>(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<amrex::ParticleReal>(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<amrex::ParticleReal>(g1*g2);
+        lab_to_COM_lorentz_factor = static_cast<amrex::ParticleReal>(g1_star*g2_star/(g1*g2));
     }
 
     /**

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<amrex::ParticleReal, double>::value),
-    "Particle precision must be double for DSMC collisions.");
-
     const amrex::ParmParse pp_collision_name(collision_name);
 
     // query for a list of collision processes