ablastr constants : add variable templates to allow for flexibility in using single/double precision

[lucafedeli88]

With this PR the constants in ablastr are first defined using variable templates and then specialized for amrex::Real, e.g., for mathematical constants:

    namespace math
    {
        namespace variable_templates
        {

            template<typename T>
            constexpr auto pi = static_cast<T>(3.14159265358979323846);

            template<typename T>
            constexpr auto tau = static_cast<T>(2.0 * pi<double>);

        } 

        constexpr auto pi = variable_templates::pi<amrex::Real>;

        constexpr auto tau = variable_templates::tau<amrex::Real>;

    } // namespace math

This is done with the intent of allowing the use of double precision constants even if the code is compiled in single precision (or vice versa if you really want it !), e.g. :

using namespace ablastr;
namespace math_vt = math::variable_templates;

const auto a = 2_rt * math::pi; // a is an amrex::Real

const auto b = 2 * math_vt::pi<double>; // b is a double !

I will need this feature to import the QED module from PICSAR into WarpX. I've started migrating the code in #5677 ,but the PR is becoming too big, so I decided to split it into smaller chunks.

The new feature is now used in the template function Algorithms::KineticEnergy . This function is conceived to be able to force double precision in some cases regardless of how WarpX has been compiled. Therefore, using c in double precision in these cases is coherent with its goal.

In case you don't like the long namespace name variable_template, an alternative would be to call it vt.

Note that I've changed the variables from static constexpr to constexpr because, to my understanding, static does not make any difference in this context.

Updates

• Following @dpgrote 's suggestion, I've changed variable_templates to const_templates
• Following @aeriforme 's suggestion, I've added c2, invc and inv_c2 constants, and I have tested them by modifying a couple of source files.

[dpgrote]

Since this is for constants, perhaps call it constant_templates?

[lucafedeli88]

Sorry for not having replied earlier, @dpgrote !
Would const_templates work according to you?

[dpgrote]

Thanks, it's better, but const has a specific meaning in C++ so that could cause confusion.

[lucafedeli88]

I think that you are right, @dpgrote . I've changed the name to constant_templates, as you originally suggested.

[dpgrote]

Would there be much penalty if the constants were all explicitly declared double? This might generate warnings of explicit down casting when single precision is being used, but would there be a performance penalty?

[lucafedeli88]

Would there be much penalty if the constants were all explicitly declared double? This might generate warnings of explicit down casting when single precision is being used, but would there be a performance penalty?

The main issue I see is that, due to our extensive use of the auto keyword, inserting just one constant in double precision may silently transform whole blocks of calculations to double precision. I am afraid that this might have an impact on performances. What I like of the approach followed in this PR is that everything works as before unless you specifically want your constants with a given precision. The verbosity of constant_templates is a feature: you only use it when you really need it!

[aeriforme]

We make massive and inconsistent use of $c^{-2}, c^{-1}, c^2$ and similar.
Should we add them here, so that we don't have to re-calculate them every time?

[lucafedeli88]

We make massive and inconsistent use of c − 2 , c − 1 , c 2 and similar. Should we add them here, so that we don't have to re-calculate them every time?

We often use constexpr when we define these constants in the code, so this should not result into an overhead.
However, we can consider adding these constants.

[lucafedeli88]

We make massive and inconsistent use of c − 2 , c − 1 , c 2 and similar. Should we add them here, so that we don't have to re-calculate them every time?

@aeriforme , I've added c2, invc, and invc2 constants, and I have tested them by modifying a couple of source files. If we like the new constants we can use them everywhere in WarpX, maybe by making the change with another PR.

[lucafedeli88]

ping, @aeriforme ;-)

[EZoni]

Thanks for the PR, @lucafedeli88! I left only a couple of minor comments.

[EZoni]

Maybe we could take this opportunity to rewrite the docstring from

//! xi: nonlinearity parameter of Heisenberg-Euler effective theory = (2.*alpha*alpha*ep0*ep0*hbar*hbar*hbar)/(45.*m_e*m_e*m_e*m_e*c*c*c*c*c)

to

//! xi: nonlinearity parameter of Heisenberg-Euler effective theory = (2.*alpha**2*ep0**2*hbar**3)/(45.*m_e**4*c**5)

[EZoni]

I'm wondering whether we could drop the templates part in the name...

Obviously it does describe what this is, as far as the implementation is concerned, but from the perspective of this being used in the code I wonder if PhysConst::constant_templates::invc2<amrex::Real> is really "needed", as opposed to, say, simply PhysConst::constants::invc2<amrex::Real>?

[lucafedeli88]

@EZoni , what about an abbreviation, such as tmpl (constant is already in PhysConst):
e.g.
PhysConst::tmpl::invc2<amrex::Real>

?

[EZoni]

Hm, I think I'm a bit confused. Can we try to summarize all namespaces we have currently in development?

I'm seeing

warpx/Source/Utils/WarpXConst.H

Lines 14 to 22 in b4009ed

	namespace MathConst
	{
	using namespace ablastr::constant::math;
	}
	namespace PhysConst
	{
	using namespace ablastr::constant::SI;
	}

and also

warpx/Source/ablastr/constant.H

Lines 15 to 89 in b4009ed

	namespace ablastr::constant
	{
	/** Mathematical constants */
	namespace math
	{
	using namespace amrex::literals;
	//! ratio of a circle's circumference to its diameter
	static constexpr amrex::Real pi = 3.14159265358979323846_rt;
	//! https://tauday.com/tau-manifesto
	static constexpr amrex::Real tau = 2.0_rt * pi;
	} // namespace math
	/** Physical constants
	*
	* Values are the 2022 CODATA recommended values
	* https://physics.nist.gov/cuu/Constants/index.html
	*
	* New additions here should also be considered for addition to
	* `warpx_constants` in WarpXUtil.cpp's `makeParser`, so that they're
	* available in parsing and evaluation of PICMI expressions, as well
	* as the corresponding Python definitions
	*/
	namespace SI
	{
	using namespace amrex::literals;
	//! vacuum speed of light [m/s]
	static constexpr auto c = 299'792'458._rt;
	//! vacuum permittivity: dielectric permittivity of vacuum [F/m]
	static constexpr auto ep0 = 8.8541878188e-12_rt;
	//! vacuum permeability: magnetic permeability of vacuum = 4.0e-7 * pi [H/m]
	//! Note that this is adjusted from the CODATA2022 value, 1.25663706127e-06
	//! So that the relation between mu0, ep0, and c is exact.
	static constexpr auto mu0 = 1.2566370612685e-06_rt;
	//! elementary charge [C]
	static constexpr auto q_e = 1.602176634e-19_rt;
	//! electron mass [kg]
	static constexpr auto m_e = 9.1093837139e-31_rt;
	//! proton mass [kg]
	static constexpr auto m_p = 1.67262192595e-27_rt;
	//! dalton: unified atomic mass unit [kg]
	static constexpr auto m_u = 1.66053906892e-27_rt;
	//! reduced Planck Constant = h / tau [J*s]
	static constexpr auto hbar = 1.0545718176461565e-34_rt;
	//! fine-structure constant = mu0/(4*pi)*q_e*q_e*c/hbar [dimensionless]
	//! The value is calculated from the expression. This differs
	//! slightly from the CODATA2022 value 0.0072973525643.
	static constexpr auto alpha = 0.0072973525643330135_rt;
	//! classical electron radius = 1./(4*pi*ep0) * q_e*q_e/(m_e*c*c) [m]
	static constexpr auto r_e = 2.8179403205e-15_rt;
	//! xi: nonlinearity parameter of Heisenberg-Euler effective theory = (2.*alpha*alpha*ep0*ep0*hbar*hbar*hbar)/(45.*m_e*m_e*m_e*m_e*c*c*c*c*c)
	static constexpr double xi = 1.3050122383827227e-52;
	//! xi times c2 = xi*c*c. This should be usable for single precision instead of xi; very close to smallest float32 number possible (1.2e-38)
	static constexpr auto xi_c2 = 1.1728865075613984e-35_rt;
	//! Boltzmann constant (exact) [J/K]
	static constexpr auto kb = 1.380649e-23_rt;
	//! 1 eV in [J]
	static constexpr auto eV = q_e;
	//! 1 MeV in [J]
	static constexpr auto MeV = q_e * 1e6_rt;
	//! 1 eV/c in [kg*m/s]
	static constexpr auto eV_invc = eV / c;
	//! 1 MeV/c in [kg*m/s]
	static constexpr auto MeV_invc = MeV / c;
	//! 1 eV/c^2 in [kg]
	static constexpr auto eV_invc2 = eV / (c * c);
	//! 1 MeV/c^2 in [kg]
	static constexpr auto MeV_invc2 = MeV / (c * c);
	} // namespace SI
	} // namespace ablastr::constant

So, we can do, e.g.,

PhysConst::c

but we cannot do

PhysConst::constant::c

right?

Where does the conflict with, say, PhysConst::constant::c<amrex::Real> (although I would use constants) come from? When I say "conflict" I refer to what you mentioned about constant being already in PhysConst.

Is the argument that PhysConst::constant_templates::c<amrex::Real> is more intuitive than PhysConst::constant::c<amrex::Real> because with the latter one doesn't immediately foresee that they should pass the template parameter?

[EZoni]

@lucafedeli88

I think we are close to finalizing this PR, as long as we resolve the latest questions raised in #5828 (comment).