Use KineticEnergy instead of getEnergy (#5786)

We currently have two functions in the code that do almost exactly the
same thing:

```
    /**
     * \brief Return (relativistic) particle energy given velocity and mass.
     * Note the use of `double` since this calculation is prone to error with
     * single precision.
     *
     * @param[in] u2 square of particle speed (i.e. u dot u where u = gamma*v)
     * @param[in] mass particle mass
     * @param[out] energy particle energy in eV
     */
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    void getEnergy ( amrex::ParticleReal const u2, double const mass,
                     double& energy )
    {
        using std::sqrt;
        using namespace amrex::literals;

        constexpr auto c2 = PhysConst::c * PhysConst::c;
        energy = mass * u2 / (sqrt(1.0_rt + u2 / c2) + 1.0_rt) / PhysConst::q_e;
    }
```

and 

```
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    amrex::ParticleReal KineticEnergy(
        const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz,
        const amrex::ParticleReal mass)
    {
        using namespace amrex;

        constexpr auto inv_c2 = 1.0_prt/(PhysConst::c * PhysConst::c);

        // The expression used is derived by reducing the expression
        // (gamma - 1)*(gamma + 1)/(gamma + 1)

        const auto u2 = ux*ux + uy*uy + uz*uz;
        const auto gamma = std::sqrt(1.0_prt + u2*inv_c2);
        return 1.0_prt/(1.0_prt + gamma)*mass*u2;
    }
```
The only differences are that `getEnergy` forces the use of double
precision and returns an energy in `eV`. I think that it would be better
if we could use only one function for clarity. Therefore, in this PR I
replaced `getEnergy` with `KineticEnergy`.

The units of the output are easy to deal with (we can just divide by 1eV
the result). The precision is trickier. @roelof-groenewald , do you
think that we really need to force `double` precision here? If it is the
case, I think that it could be preferable to make `getEnergy` a template
function and force double precision with something like
`KineticEnergy<double>(...)`. What do you think?

---------

Author: lucafedeli88

Source/Particles/Algorithms/KineticEnergy.H

@@ -12,6 +12,7 @@
 
 #include "AMReX_Extension.H"
 #include "AMReX_GpuQualifiers.H"
+#include "AMReX_Math.H"
 #include "AMReX_REAL.H"
 
 #include <cmath>
@@ -22,28 +23,30 @@ namespace Algorithms{
      * \brief Computes the kinetic energy of a particle.
      * This method should not be used with photons.
      *
+     * @tparam T floating-point type to be used in internal calculations (by default equal to amrex::ParticleReal)
      * @param[in] ux x component of the particle momentum (code units)
      * @param[in] uy y component of the particle momentum (code units)
      * @param[in] uz z component of the particle momentum (code units)
      * @param[in] mass mass of the particle (in S.I. units)
      *
      * @return the kinetic energy of the particle (in S.I. units)
      */
+    template <typename T = amrex::ParticleReal>
     AMREX_GPU_HOST_DEVICE AMREX_INLINE
-    amrex::ParticleReal KineticEnergy(
-        const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz,
-        const amrex::ParticleReal mass)
+    T KineticEnergy(
+        const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::ParticleReal mass)
     {
         using namespace amrex;
 
-        constexpr auto inv_c2 = 1.0_prt/(PhysConst::c * PhysConst::c);
+        constexpr auto one = static_cast<T>(1.0);
+        constexpr auto inv_c2 = one/Math::powi<2>(static_cast<T>(PhysConst::c));
 
         // The expression used is derived by reducing the expression
         // (gamma - 1)*(gamma + 1)/(gamma + 1)
 
-        const auto u2 = ux*ux + uy*uy + uz*uz;
-        const auto gamma = std::sqrt(1.0_prt + u2*inv_c2);
-        return 1.0_prt/(1.0_prt + gamma)*mass*u2;
+        const auto u2 = static_cast<T>(ux*ux + uy*uy + uz*uz);
+        const auto gamma = std::sqrt(one + u2*inv_c2);
+        return one/(one + gamma)*static_cast<T>(mass)*u2;
     }
 
     /**

Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp

@@ -7,6 +7,7 @@
 #include "BackgroundMCCCollision.H"
 
 #include "ImpactIonization.H"
+#include "Particles/Algorithms/KineticEnergy.H"
 #include "Particles/ParticleCreation/FilterCopyTransform.H"
 #include "Particles/ParticleCreation/SmartCopy.H"
 #include "Utils/Parser/ParserUtils.H"
@@ -435,8 +436,8 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile
                                   // subtract any energy penalty of the collision from the
                                   // projectile energy
                                   if (scattering_process.m_energy_penalty > 0.0_prt) {
-                                      ParticleUtils::getEnergy(v_coll2, m, E_coll);
-                                      E_coll = (E_coll - scattering_process.m_energy_penalty) * PhysConst::q_e;
+                                      constexpr auto eV = PhysConst::q_e;
+                                      E_coll = (Algorithms::KineticEnergy<double>(vx, vy, vz, m) - scattering_process.m_energy_penalty*eV);
                                       const auto scale_fac = static_cast<amrex::ParticleReal>(
                                         std::sqrt(E_coll * (E_coll + 2.0_prt*mc2) / c2) / m / v_coll);
                                       vx *= scale_fac;

Source/Particles/Collision/BackgroundMCC/ImpactIonization.H

@@ -9,6 +9,7 @@
 
 #include "Particles/Collision/ScatteringProcess.H"
 
+#include "Particles/Algorithms/KineticEnergy.H"
 #include "Utils/ParticleUtils.H"
 #include "Utils/WarpXConst.H"
 
@@ -94,13 +95,14 @@ public:
         const ParticleReal uz = ptd.m_rdata[PIdx::uz][i];
 
         // calculate kinetic energy
-        const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
-        ParticleUtils::getEnergy(u_coll2, m_mass, E_coll);
+        constexpr auto eV = PhysConst::q_e;
+        E_coll = Algorithms::KineticEnergy<double>(ux, uy, uz, m_mass)/eV;
 
         // get collision cross-section
         const ParticleReal sigma_E = m_mcc_process.getCrossSection(static_cast<amrex::ParticleReal>(E_coll));
 
         // calculate normalized collision frequency
+        const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
         const ParticleReal nu_i = n_a * sigma_E * sqrt(u_coll2) / m_nu_max;
 
         // check if this collision should be performed
@@ -198,11 +200,11 @@ public:
         auto& i_uz = dst2.m_rdata[PIdx::uz][i_dst2];
 
         // calculate kinetic energy
-        const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
-        ParticleUtils::getEnergy(u_coll2, m_mass1, E_coll);
+        constexpr auto eV = PhysConst::q_e;
+        E_coll = Algorithms::KineticEnergy<double>(ux, uy, uz, m_mass1)/eV;
 
         // each electron gets half the energy (could change this later)
-        const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) / 2.0_prt * PhysConst::q_e);
+        const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) / 2.0_prt * eV);
 
         // precalculate often used value
         constexpr auto c2 = PhysConst::c * PhysConst::c;

Source/Utils/ParticleUtils.H

@@ -33,26 +33,6 @@ namespace ParticleUtils {
                              amrex::MFIter const & mfi,
                              WarpXParticleContainer::ParticleTileType & ptile);
 
-    /**
-     * \brief Return (relativistic) particle energy given velocity and mass.
-     * Note the use of `double` since this calculation is prone to error with
-     * single precision.
-     *
-     * @param[in] u2 square of particle speed (i.e. u dot u where u = gamma*v)
-     * @param[in] mass particle mass
-     * @param[out] energy particle energy in eV
-     */
-    AMREX_GPU_HOST_DEVICE AMREX_INLINE
-    void getEnergy ( amrex::ParticleReal const u2, double const mass,
-                     double& energy )
-    {
-        using std::sqrt;
-        using namespace amrex::literals;
-
-        constexpr auto c2 = PhysConst::c * PhysConst::c;
-        energy = mass * u2 / (sqrt(1.0_rt + u2 / c2) + 1.0_rt) / PhysConst::q_e;
-    }
-
     /**
      * \brief Return (relativistic) collision energy assuming the target (with
      * mass M) is stationary and the projectile is approaching with the