refactor AU units

- add conversion au units -> si - add conversion au unis -> pic ci: picongpu
ComputationalRadiationPhysics · Sep 23, 2024 · 6f9b7a8 · 6f9b7a8
1 parent a839023
commit 6f9b7a8
Show file tree

Hide file tree

Showing 14 changed files with 282 additions and 222 deletions.
diff --git a/docs/source/usage/workflows/ionization.rst b/docs/source/usage/workflows/ionization.rst
@@ -69,13 +69,13 @@ To include charge-state-only simulations in your simulations you must:
                 float_X,
                 7, // number charge states
                 Nitrogen, // name you may reference this by later, remember to prepend the namespace and append _t!
-                14.53413 * UNITCONV_eV_to_AU,
-                29.60125 * UNITCONV_eV_to_AU,
-                47.4453 * UNITCONV_eV_to_AU,
-                77.4735 * UNITCONV_eV_to_AU,
-                97.89013 * UNITCONV_eV_to_AU,
-                552.06731 * UNITCONV_eV_to_AU,
-                667.04609 * UNITCONV_eV_to_AU);
+                sim.si.conv.ev2auEnergy(14.53413),
+                sim.si.conv.ev2auEnergy(29.60125),
+                sim.si.conv.ev2auEnergy(47.4453),
+                sim.si.conv.ev2auEnergy(77.4735),
+                sim.si.conv.ev2auEnergy(97.89013),
+                sim.si.conv.ev2auEnergy(552.06731),
+                sim.si.conv.ev2auEnergy(667.04609));
         }; // namespace picongpu::ionization::energies::AU
 
     .. note::

diff --git a/include/picongpu/defines.hpp b/include/picongpu/defines.hpp
@@ -40,5 +40,4 @@ namespace picongpu
 #include "picongpu/param/speciesConstants.param"
 #include "picongpu/param/simulation.param"
 #include "picongpu/unitless/simulation.unitless"
-#include "picongpu/unitless/physicalConstants.unitless"
 // clang-format on
diff --git a/include/picongpu/param/ionizationEnergies.param b/include/picongpu/param/ionizationEnergies.param
diff --git a/include/picongpu/param/physicalConstants.param b/include/picongpu/param/physicalConstants.param
@@ -46,26 +46,5 @@ namespace picongpu
         constexpr float_64 ELECTRON_MASS_SI = 9.1093837139e-31;
         //! unit: C,    2022 CODATA value, https://physics.nist.gov/cgi-bin/cuu/Value?e
         constexpr float_64 ELECTRON_CHARGE_SI = -1.602176634e-19;
-
-        /* atomic unit for energy:
-         * 2 Rydberg = 27.21 eV --> converted to Joule
-         */
-        constexpr float_64 ATOMIC_UNIT_ENERGY = 4.36e-18;
-
-        /* atomic unit for electric field in V/m:
-         * field strength between electron and core in ground state hydrogen
-         */
-        constexpr float_64 ATOMIC_UNIT_EFIELD = 5.14e11;
-
-        /* atomic unit for time in s:
-         * 150 attoseconds (classical electron orbit time in hydrogen) / 2 PI
-         */
-        constexpr float_64 ATOMIC_UNIT_TIME = 2.4189e-17;
-
     } // namespace SI
-
-    /* 1 atomic unit of energy is equal to 1 Hartree or 2 Rydberg
-     * which is twice the ground state binding energy of atomic hydrogen */
-    constexpr float_64 UNITCONV_AU_to_eV = 27.21139;
-    constexpr float_64 UNITCONV_eV_to_AU = (1.0 / UNITCONV_AU_to_eV);
 } // namespace picongpu
diff --git a/include/picongpu/particles/atomicPhysics/kernel/DecelerateElectrons.kernel b/include/picongpu/particles/atomicPhysics/kernel/DecelerateElectrons.kernel
@@ -142,7 +142,7 @@ namespace picongpu::particles::atomicPhysics::kernel
                     // sqrt(E' * (E' + 2* m*c^2)) / c
                     float_X const newElectronMomentum
                         = math::sqrt(newEnergyElectron * (newEnergyElectron + 2 * mcSquaredElectron)) * scalingFactor;
-                    // AU = ATOMIC_UNIT_ENERGY
+                    // AU = sim.si.conv.auEnergy2Joule()
                     // sqrt(eV * (eV + eV))/(sim.unit.length()/sim.unit.time()) * (sim.unit.energy()/eV)
                     //  = sqrt((eV)^2)/(eV) * sim.unit.time()/sim.unit.length() * sim.unit.energy()
                     //  = sim.unit.mass() * sim.unit.length()^2/sim.unit.time()^2 * sim.unit.time()/sim.unit.length()

diff --git a/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp b/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp
@@ -92,8 +92,8 @@ namespace picongpu::particles::atomicPhysics::rateCalculation
                 chargeStateDataBox);
             // unitless
             float_X const effectivePrincipalQuantumNumber
-                = screenedCharge / math::sqrt(2._X * ionizationEnergy * picongpu::UNITCONV_eV_to_AU);
-            float_X const eFieldNorm_AU = eFieldNorm / ATOMIC_UNIT_EFIELD;
+                = screenedCharge / math::sqrt(2._X * ionizationEnergy * sim.si.conv.ev2auEnergy(1.0));
+            float_X const eFieldNorm_AU = sim.pic.conv.eField2auEField(eFieldNorm);
             float_X const screenedChargeCubed = pmacc::math::cPow(screenedCharge, 3u);
             float_X const dBase = 4.0_X * math::exp(1._X) * screenedChargeCubed
                 / (eFieldNorm_AU * pmacc::math::cPow(effectivePrincipalQuantumNumber, 4u));
@@ -103,7 +103,7 @@ namespace picongpu::particles::atomicPhysics::rateCalculation
             constexpr float_X pi = pmacc::math::Pi<float_X>::value;
             float_X const nEffCubed = pmacc::math::cPow(effectivePrincipalQuantumNumber, 3u);
 
-            // 1/ATOMIC_UNIT_TIME
+            // 1/atomicTime
             float_X rateADK_AU = eFieldNorm_AU * pmacc::math::cPow(dFromADK, 2u) / (8._X * pi * screenedCharge)
                 * math::exp(-2._X * screenedChargeCubed / (3._X * nEffCubed * eFieldNorm_AU));
 
@@ -112,7 +112,7 @@ namespace picongpu::particles::atomicPhysics::rateCalculation
                 u32(T_ADKLaserPolarization) == u32(atomicPhysics::enums::ADKLaserPolarization::linearPolarization))
                 rateADK_AU *= math::sqrt(3._X * nEffCubed * eFieldNorm_AU / (pi * screenedChargeCubed));
 
-            return rateADK_AU / ATOMIC_UNIT_TIME;
+            return rateADK_AU / sim.pic.conv.auTime2Time(1.0);
         }
     };
 } // namespace picongpu::particles::atomicPhysics::rateCalculation
diff --git a/include/picongpu/particles/ionization/byField/ADK/AlgorithmADK.hpp b/include/picongpu/particles/ionization/byField/ADK/AlgorithmADK.hpp
@@ -84,7 +84,7 @@ namespace picongpu
 
                         constexpr float_X pi = pmacc::math::Pi<float_X>::value;
                         /* electric field in atomic units - only absolute value */
-                        float_X const eInAU = pmacc::math::l2norm(eField) / ATOMIC_UNIT_EFIELD;
+                        float_X const eInAU = sim.pic.conv.eField2auEField(pmacc::math::l2norm(eField));
 
                         /* the charge that attracts the electron that is to be ionized:
                          * equals `protonNumber - #allInnerElectrons`
@@ -113,7 +113,7 @@ namespace picongpu
                         }
 
                         /* simulation time step in atomic units */
-                        auto const timeStepAU = float_X(sim.pic.getDt() / ATOMIC_UNIT_TIME);
+                        auto const timeStepAU = sim.pic.conv.time2auTime(sim.pic.getDt());
                         /* ionization probability
                          *
                          * probability = rate * time step

diff --git a/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSI.hpp b/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSI.hpp
@@ -77,7 +77,7 @@ namespace picongpu
                         /* critical field strength in atomic units */
                         float_X const critField = iEnergy * iEnergy / (float_X(4.0) * effectiveCharge);
                         /* ionization condition */
-                        if(pmacc::math::l2norm(eField) / ATOMIC_UNIT_EFIELD >= critField)
+                        if(sim.pic.conv.eField2auEField(pmacc::math::l2norm(eField)) >= critField)
                         {
                             /* return ionization energy and number of macro electrons to produce */
                             return IonizerReturn{iEnergy, 1u};

diff --git a/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSIEffectiveZ.hpp b/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSIEffectiveZ.hpp
@@ -76,7 +76,7 @@ namespace picongpu
                         float_X critField = iEnergy * iEnergy / (float_X(4.0) * ZEff);
 
                         /* ionization condition */
-                        if(pmacc::math::l2norm(eField) / ATOMIC_UNIT_EFIELD >= critField)
+                        if(sim.pic.conv.eField2auEField(pmacc::math::l2norm(eField)) >= critField)
                         {
                             /* return ionization energy and number of macro electrons to produce */
                             return IonizerReturn{iEnergy, 1u};

diff --git a/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSIStarkShifted.hpp b/include/picongpu/particles/ionization/byField/BSI/AlgorithmBSIStarkShifted.hpp
@@ -75,7 +75,7 @@ namespace picongpu
                             = (math::sqrt(float_X(2.)) - float_X(1.)) * math::pow(iEnergy, float_X(3. / 2.));
 
                         /* ionization condition */
-                        if(pmacc::math::l2norm(eField) / ATOMIC_UNIT_EFIELD >= critField)
+                        if(sim.pic.conv.eField2auEField(pmacc::math::l2norm(eField)) >= critField)
                         {
                             /* return ionization energy number of electrons to produce */
                             return IonizerReturn{iEnergy, 1u};

diff --git a/include/picongpu/particles/ionization/byField/IonizationCurrent/IonizationCurrent.hpp b/include/picongpu/particles/ionization/byField/IonizationCurrent/IonizationCurrent.hpp
@@ -61,8 +61,7 @@ namespace picongpu
                      * ionization current. */
                     if(retValue.ionizationEnergy != 0.0_X)
                     {
-                        auto ionizationEnergy = weighting * retValue.ionizationEnergy * SI::ATOMIC_UNIT_ENERGY
-                            / sim.unit.energy(); // convert to PIConGPU units
+                        auto ionizationEnergy = sim.pic.conv.auEnergy2Joule(weighting * retValue.ionizationEnergy);
                         /* calculate ionization current at particle position */
                         float3_X jIonizationPar = JIonizationCalc{}(ionizationEnergy, eField);
                         /* assign ionization current to grid points */

diff --git a/include/picongpu/particles/ionization/byField/Keldysh/AlgorithmKeldysh.hpp b/include/picongpu/particles/ionization/byField/Keldysh/AlgorithmKeldysh.hpp
@@ -76,7 +76,7 @@ namespace picongpu
 
                         constexpr float_X pi = pmacc::math::Pi<float_X>::value;
                         /* electric field in atomic units - only absolute value */
-                        float_X eInAU = pmacc::math::l2norm(eField) / ATOMIC_UNIT_EFIELD;
+                        float_X eInAU = sim.pic.conv.eField2auEField(pmacc::math::l2norm(eField));
 
                         /* factor two avoid calculation math::pow(2,5./4.); */
                         const float_X twoToFiveQuarters = 2.3784142300054;
@@ -89,7 +89,7 @@ namespace picongpu
                             * math::sqrt(float_X(1.) / charExpArg) * math::exp(-float_X(2. / 3.) * charExpArg);
 
                         /* simulation time step in atomic units */
-                        const auto timeStepAU = float_X(sim.pic.getDt() / ATOMIC_UNIT_TIME);
+                        const auto timeStepAU = sim.pic.conv.time2auTime(sim.pic.getDt());
                         /* ionization probability
                          *
                          * probability = rate * time step

diff --git a/include/picongpu/unitless/simulation.unitless b/include/picongpu/unitless/simulation.unitless
@@ -127,6 +127,35 @@ namespace picongpu
             {
                 return static_cast<T_Type>(si.getClassicalElectronRadius() / unit.length());
             }
+
+            struct Convert
+            {
+                template<typename T_Type = float_X>
+                constexpr T_Type eField2auEField(T_Type eField = 1.0) const
+                {
+                    return static_cast<T_Type>(si.conv.eField2auEField(eField * unit.eField()));
+                }
+
+                template<typename T_Type = float_X>
+                constexpr T_Type auEnergy2Joule(T_Type auEnergy = 1.0) const
+                {
+                    return static_cast<T_Type>(si.conv.auEnergy2Joule(auEnergy) / unit.energy());
+                }
+
+                template<typename T_Type = float_X>
+                constexpr T_Type auTime2Time(T_Type auTime = 1.0) const
+                {
+                    return static_cast<T_Type>(si.conv.auTime2Time(auTime) / unit.time());
+                }
+
+                template<typename T_Type = float_X>
+                constexpr T_Type time2auTime(T_Type time = 1.0) const
+                {
+                    return static_cast<T_Type>(si.conv.time2auTime(time * unit.time()));
+                }
+            };
+
+            static constexpr Convert conv{};
         };
         struct SiUnits
         {
@@ -289,6 +318,55 @@ namespace picongpu
                 {
                     return joule / ev2Joule(1.0);
                 }
+
+                /* atomic unit for energy:
+                 * 2 Rydberg = 27.21 eV --> converted to Joule
+                 */
+                constexpr float_64 auEnergy2Joule(float_64 auEnergy = 1.0) const
+                {
+                    return auEnergy * 2.0 * si.getRydbergEnergy();
+                }
+
+                constexpr float_64 joule2auEnergy(float_64 joule = 1.0) const
+                {
+                    return joule / auEnergy2Joule(1.0);
+                }
+
+                constexpr float_64 auEnergy2ev(float_64 auEnergy = 1.0) const
+                {
+                    return joule2ev(auEnergy2Joule(auEnergy));
+                }
+
+                constexpr float_64 ev2auEnergy(float_64 ev = 1.0) const
+                {
+                    return ev / auEnergy2ev(1.0);
+                }
+
+                /* atomic unit for electric field in V/m:
+                 * field strength between electron and core in ground state hydrogen
+                 */
+                constexpr float_64 auEField2EField(float_64 auEField = 1.0) const
+                {
+                    return auEField * 5.14e11;
+                }
+
+                constexpr float_64 eField2auEField(float_64 eField = 1.0) const
+                {
+                    return eField / auEField2EField();
+                }
+
+                /* atomic unit for time in s:
+                 * 150 attoseconds (classical electron orbit time in hydrogen) / 2 PI
+                 */
+                constexpr float_64 auTime2Time(float_64 auTime = 1.0) const
+                {
+                    return auTime * 2.4189e-17;
+                }
+
+                constexpr float_64 time2auTime(float_64 time = 1.0) const
+                {
+                    return time / auTime2Time(1.0);
+                }
             };
 
             static constexpr Convert conv{};

diff --git a/share/picongpu/tests/compileAtomicPhysics/include/picongpu/param/ionizationEnergies.param b/share/picongpu/tests/compileAtomicPhysics/include/picongpu/param/ionizationEnergies.param
@@ -39,7 +39,7 @@ namespace picongpu
          * Please follow these rules for defining ionization energies of atomic species,
          * unless your chosen ionization model requires a different unit system than `AU::`
          *     - input of values in either atomic units or converting eV or Joule to them
-         *     -> use either UNITCONV_eV_to_AU or SI::ATOMIC_UNIT_ENERGY for that purpose
+         *     -> use either sim.si.conv.ev2auEnergy() or sim.si.conv.joule2auEnergy() for that purpose
          *     - use `float_X` as the preferred data type
          *
          * example:
@@ -76,24 +76,24 @@ namespace picongpu
                     float_X,
                     18,
                     Argon,
-                    15.7596119 * UNITCONV_eV_to_AU,
-                    27.62967 * UNITCONV_eV_to_AU,
-                    40.735 * UNITCONV_eV_to_AU,
-                    59.58 * UNITCONV_eV_to_AU,
-                    74.84 * UNITCONV_eV_to_AU,
-                    91.290 * UNITCONV_eV_to_AU,
-                    124.41 * UNITCONV_eV_to_AU,
-                    143.4567 * UNITCONV_eV_to_AU,
-                    422.60 * UNITCONV_eV_to_AU,
-                    479.76 * UNITCONV_eV_to_AU,
-                    540.4 * UNITCONV_eV_to_AU,
-                    619.0 * UNITCONV_eV_to_AU,
-                    685.5 * UNITCONV_eV_to_AU,
-                    755.13 * UNITCONV_eV_to_AU,
-                    855.5 * UNITCONV_eV_to_AU,
-                    918.375 * UNITCONV_eV_to_AU,
-                    4120.6657 * UNITCONV_eV_to_AU,
-                    4426.2229 * UNITCONV_eV_to_AU);
+                    sim.si.conv.ev2auEnergy(15.7596119),
+                    sim.si.conv.ev2auEnergy(27.62967),
+                    sim.si.conv.ev2auEnergy(40.735),
+                    sim.si.conv.ev2auEnergy(59.58),
+                    sim.si.conv.ev2auEnergy(74.84),
+                    sim.si.conv.ev2auEnergy(91.290),
+                    sim.si.conv.ev2auEnergy(124.41),
+                    sim.si.conv.ev2auEnergy(143.4567),
+                    sim.si.conv.ev2auEnergy(422.60),
+                    sim.si.conv.ev2auEnergy(479.76),
+                    sim.si.conv.ev2auEnergy(540.4),
+                    sim.si.conv.ev2auEnergy(619.0),
+                    sim.si.conv.ev2auEnergy(685.5),
+                    sim.si.conv.ev2auEnergy(755.13),
+                    sim.si.conv.ev2auEnergy(855.5),
+                    sim.si.conv.ev2auEnergy(918.375),
+                    sim.si.conv.ev2auEnergy(4120.6657),
+                    sim.si.conv.ev2auEnergy(4426.2229));
 
             } // namespace AU
         } // namespace energies