Use KineticEnergy instead of getEnergy

[lucafedeli88]

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?

[roelof-groenewald]

Thanks for the clean-up @lucafedeli88! We do need to force double precision for the cases where getEnergy is used as we have seen incorrect physics results when allowing single precision in that function.

[lucafedeli88]

Thanks for the clean-up @lucafedeli88! We do need to force double precision for the cases where getEnergy is used as we have seen incorrect physics results when allowing single precision in that function.

I see!
Would this solution be OK ?

    template <typename T>
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    T KineticEnergy(
        const T ux, const T uy, const T uz, const T mass)
    {
        using namespace amrex;

        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(one + u2*inv_c2);
        return one/(one + gamma)*mass*u2;
    }

Double precision can be forced like that: Algorithms::KineticEnergy<double>(vx, vy, vz, m)

[roelof-groenewald]

The particle velocity components would be of type ParticleReal. Won't this cause a problem since it now expects all inputs to be of the templated type T (which is specified to be double)?
Could this be:

Suggested change

	const T ux, const T uy, const T uz, const T mass)
	const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const T mass)

[lucafedeli88]

It should be fine if you force double precision, since it would be either a no-op or a promotion. But you are right, it could be an issue if for some reason amrex::ParticleReal is double and you try to force single precision.

[lucafedeli88]

This could be a solution that avoids this issue:

    /**
     * \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 equalt 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
    T KineticEnergy(
        const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::ParticleReal mass)
    {
        using namespace amrex;

        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 = 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;
    }

T is by default equal to ParticleReal. However, you can override this and then calculations are carried out with precision T and the result is of type T. Conversions are carried out using static_cast, so it should work ok.

[roelof-groenewald]

Thanks for the code cleaning!