Implement Reviewer suggestions

Author: Marco Garten

src/picongpu/include/particles/ionization/byCollision/ThomasFermi/AlgorithmThomasFermi.hpp

@@ -20,12 +20,11 @@
 
 #pragma once
 
-#include "pmacc_types.hpp"
 #include "simulation_defines.hpp"
 #include "particles/traits/GetAtomicNumbers.hpp"
 #include "traits/attribute/GetChargeState.hpp"
+
 #include "algorithms/math/floatMath/floatingPoint.tpp"
-#include "particles/ionization/byCollision/ThomasFermi/TFFittingParameters.def"
 
 /** \file AlgorithmThomasFermi.hpp
  *
@@ -48,11 +47,13 @@ namespace ionization
      * \brief calculation for the Thomas-Fermi pressure ionization model
      *
      * This model uses local density and "temperature" values as input
-     * parameters. To be able to speak of a "temperature" an equilibrium state
-     * would be required. Typical high power laser-plasma interaction is highly
+     * parameters. A physical temperature requires a defined equilibrium state.
+     * Typical high power laser-plasma interaction is highly
      * non-equilibrated, though. The name "temperature" is kept to illustrate
      * the origination from the Thomas-Fermi model. It is nevertheless
-     * more accurate to think of it as an averaged kinetic energy.
+     * more accurate to think of it as an averaged kinetic energy
+     * which is not backed by the model and should therefore only be used with
+     * a certain suspicion in such Non-LTE scenarios.
      */
     struct AlgorithmThomasFermi
     {
@@ -68,86 +69,87 @@ namespace ionization
          * \param randNr random number
          */
         template<typename KinEnergyDensityType, typename DensityType, typename ParticleType >
-        HDINLINE void
+        HDINLINE float_X
         operator()( const KinEnergyDensityType kinEnergyDensity, const DensityType density, ParticleType& parentIon, float_X randNr )
         {
 
             /* @TODO replace the float_64 with float_X and make sure the values are scaled to PIConGPU units */
-            const float_64 protonNumber = GetAtomicNumbers<ParticleType>::type::numberOfProtons;
-            const float_64 neutronNumber = GetAtomicNumbers<ParticleType>::type::numberOfNeutrons;
+            float_64 constexpr protonNumber = GetAtomicNumbers<ParticleType>::type::numberOfProtons;
+            float_64 constexpr neutronNumber = GetAtomicNumbers<ParticleType>::type::numberOfNeutrons;
             float_64 chargeState = attribute::getChargeState(parentIon);
 
             /* ionization condition */
             if (chargeState < protonNumber)
             {
                 /* atomic mass number (usually A) A = N + Z */
-                float_64 massNumber = neutronNumber + protonNumber;
+                float_64 constexpr massNumber = neutronNumber + protonNumber;
 
                 /** @TODO replace the static_cast<float_64> by casts to float_X
                  * or leave out entirely and compute everything in PIConGPU scaled units
                  */
-                float_64 const densityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::Density().getUnit()[0]);
-                float_64 const kinEnergyDensityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::EnergyDensity().getUnit()[0]);
+                float_64 constexpr densityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::Density().getUnit()[0]);
+                float_64 constexpr kinEnergyDensityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::EnergyDensity().getUnit()[0]);
                 /* convert from kinetic energy density to average kinetic energy per particle */
-                float_64 kinEnergy = (kinEnergyDensity / density) * (kinEnergyDensityUnit / densityUnit);
+                float_64 constexpr kinEnergyUnit = kinEnergyDensityUnit / densityUnit;
+                float_64 const kinEnergy = (kinEnergyDensity / density) * kinEnergyUnit;
                 /** convert kinetic energy in J to "temperature" in eV by assuming an ideal electron gas
                  * E_kin = 3/2 k*T
                  */
-                float_64 temperature = kinEnergy * UNITCONV_Joule_to_keV * float_64(1.e3) * float_64(2./3.);
+                float_64 constexpr convKinEnergyToTemperature = UNITCONV_Joule_to_keV * float_64(1.e3) * float_64(2./3.);
+                float_64 const temperature = kinEnergy * convKinEnergyToTemperature;
 
-                float_64 T_0 = temperature/math::pow(protonNumber,float_64(4./3.));
+                float_64 const T_0 = temperature/math::pow(protonNumber,float_64(4./3.));
 
-                float_64 T_F = T_0 / (float_64(1.) + T_0);
+                float_64 const T_F = T_0 / (float_64(1.) + T_0);
 
-                /* for all the fitting parameters @see TFFittingParameters.def */
+                /* for all the fitting parameters @see ionizerConfig.param */
 
                 /** this is weird - I have to define temporary variables because
                  * otherwise the math::pow function won't recognize those at the
                  * exponent position */
-                float_64 TFA2_temp = TFA2;
-                float_64 TFA4_temp = TFA4;
-                float_64 TFBeta_temp = TFBeta;
+                float_64 constexpr TFA2_temp = TFA2;
+                float_64 constexpr TFA4_temp = TFA4;
+                float_64 constexpr TFBeta_temp = TFBeta;
 
-                float_64 A = TFA1 * math::pow(T_0,TFA2_temp) + TFA3 * math::pow(T_0,TFA4_temp);
+                float_64 const A = TFA1 * math::pow(T_0,TFA2_temp) + TFA3 * math::pow(T_0,TFA4_temp);
 
-                float_64 B = -math::exp(TFB0 + TFB1*T_F + TFB2*math::pow(T_F,float_64(7.)));
+                float_64 const B = -math::exp(TFB0 + TFB1*T_F + TFB2*math::pow(T_F,float_64(7.)));
 
-                float_64 C = TFC1 * T_F + TFC2;
+                float_64 const C = TFC1 * T_F + TFC2;
 
-                /** requires mass density in g/cm^3
-                 * @TODO relocate constants to param file or leave out and calculate unitless
-                 */
-                float_64 const nAvogadro = 6.022e23;
-                float_64 const convM3ToCM3 = 1.e6;
+                /* requires mass density in g/cm^3 */
+                float_64 constexpr nAvogadro = SI::N_AVOGADRO;
+                float_64 constexpr convM3ToCM3 = 1.e6;
 
-                float_64 massDensity = density * densityUnit * massNumber / nAvogadro / convM3ToCM3;
-                float_64 R = massDensity/(protonNumber * massNumber);
+                float_64 constexpr convToMassDensity = densityUnit * massNumber / nAvogadro / convM3ToCM3;
+                float_64 const massDensity = density * convToMassDensity;
 
-                float_64 Q_1 = A * math::pow(R,B);
+                float_64 constexpr atomicTimesMassNumber = protonNumber * massNumber;
+                float_64 const R = massDensity/atomicTimesMassNumber;
 
-                float_64 Q = math::pow(math::pow(R,C) + math::pow(Q_1,C), float_64(1.)/C);
+                float_64 const Q_1 = A * math::pow(R,B);
 
-                float_64 x = TFAlpha * math::pow(Q,TFBeta_temp);
+                float_64 const Q = math::pow(math::pow(R,C) + math::pow(Q_1,C), float_64(1.)/C);
+
+                float_64 const x = TFAlpha * math::pow(Q,TFBeta_temp);
 
                 /* Thomas-Fermi average ionization state */
-                float_64 ZStar = protonNumber * x / (float_64(1.) + x + math::sqrt(float_64(1.) + float_64(2.)*x));
+                float_64 const ZStar = protonNumber * x / (float_64(1.) + x + math::sqrt(float_64(1.) + float_64(2.)*x));
 
                 /* integral part of the charge state */
                 float_64 intZStar;
                 /* fractional part of the charge state */
                 float_X fracZStar = static_cast<float_X>(math::modf(ZStar,&intZStar));
 
                 /* determine charge state */
-                float_X chargeState = static_cast<float_X>(intZStar) + float_X(1.0)*(randNr < fracZStar);
+                float_X const chargeState = static_cast<float_X>(intZStar) + float_X(1.0)*(randNr < fracZStar);
 
                 /** determine the new number of bound electrons from the TF ionization state
                  * @TODO introduce partial macroparticle ionization / ionization distribution at some point
                  */
-                float_X newBoundElectrons = protonNumber - chargeState;
-                /* safety check to avoid double counting since recombination is not yet implemented */
-                if (newBoundElectrons < parentIon[boundElectrons_])
-                    /* update the particle attribute only if more free electrons are to be created */
-                    parentIon[boundElectrons_] = newBoundElectrons;
+                float_X const newBoundElectrons = protonNumber - chargeState;
+
+                return newBoundElectrons;
             }
 
         }

src/picongpu/include/particles/ionization/byCollision/ThomasFermi/TFFittingParameters.def

@@ -40,7 +40,6 @@ namespace ionization
     struct ThomasFermi_Impl;
 
     /** Thomas-Fermi impact ionization model
-     *
      *
      * This ionization model describes the atom inside the Thomas-Fermi framework
      * in a self-consistent way. There the electrons are modeled as a density

src/picongpu/include/particles/ionization/byCollision/ThomasFermi/ThomasFermi.def

@@ -20,27 +20,26 @@
 
 #pragma once
 
-#include <boost/type_traits/integral_constant.hpp>
-
 #include "simulation_defines.hpp"
 #include "traits/Resolve.hpp"
 #include "traits/UsesRNG.hpp"
-#include "mappings/kernel/AreaMapping.hpp"
 
 #include "fields/FieldTmp.hpp"
 
 #include "particles/ionization/byCollision/ThomasFermi/ThomasFermi.def"
 #include "particles/ionization/byCollision/ThomasFermi/AlgorithmThomasFermi.hpp"
 #include "particles/ionization/ionization.hpp"
-
-#include "compileTime/conversion/TypeToPointerPair.hpp"
-#include "memory/boxes/DataBox.hpp"
-
 #include "particles/ionization/ionizationMethods.hpp"
 
 #include "random/methods/XorMin.hpp"
 #include "random/distributions/Uniform.hpp"
 #include "random/RNGProvider.hpp"
+#include "dataManagement/DataConnector.hpp"
+#include "compileTime/conversion/TypeToPointerPair.hpp"
+#include "memory/boxes/DataBox.hpp"
+#include "mappings/kernel/AreaMapping.hpp"
+
+#include <boost/type_traits/integral_constant.hpp>
 
 namespace picongpu
 {
@@ -72,8 +71,8 @@ namespace ionization
     struct ThomasFermi_Impl
     {
 
-        typedef T_DestSpecies DestSpecies;
-        typedef T_SrcSpecies  SrcSpecies;
+        using DestSpecies = T_DestSpecies;
+        using SrcSpecies = T_SrcSpecies;
 
         typedef typename SrcSpecies::FrameType FrameType;
 
@@ -106,7 +105,7 @@ namespace ionization
             typedef typename RNGFactory::GetRandomType<Distribution>::type RandomGen;
             RandomGen randomGen;
 
-            typedef MappingDesc::SuperCellSize TVec;
+            using SuperCellSize = MappingDesc::SuperCellSize;
 
             typedef FieldTmp::ValueType ValueType_Rho;
             typedef FieldTmp::ValueType ValueType_Ene;
@@ -234,10 +233,10 @@ namespace ionization
                 /* alias for the single macro-particle */
                 auto particle = ionFrame[localIdx];
                 /* particle position, used for field-to-particle interpolation */
-                floatD_X pos = particle[position_];
-                const int particleCellIdx = particle[localCellIdx_];
+                floatD_X const pos = particle[position_];
+                int const particleCellIdx = particle[localCellIdx_];
                 /* multi-dim coordinate of the local cell inside the super cell */
-                DataSpace<TVec::dim> localCell(DataSpaceOperations<TVec::dim>::template map<TVec > (particleCellIdx));
+                DataSpace<SuperCellSize::dim> localCell(DataSpaceOperations<SuperCellSize::dim>::template map<SuperCellSize > (particleCellIdx));
                 /* interpolation of density */
                 const fieldSolver::numericalCellType::traits::FieldPosition<FieldTmp> fieldPosRho;
                 ValueType_Rho densityV = Field2ParticleInterpolation()
@@ -257,14 +256,19 @@ namespace ionization
 
                 /* this is the point where actual ionization takes place */
                 IonizationAlgorithm ionizeAlgo;
-                ionizeAlgo(
+                float_X newBoundElectrons = ionizeAlgo(
                     kinEnergyDensity, density,
                     particle, this->randomGen()
                     );
 
-                /* determine number of new macro electrons to be created */
-                newMacroElectrons = prevBoundElectrons - particle[boundElectrons_];
+                particle[boundElectrons_] = newBoundElectrons;
 
+                /* safety check to avoid double counting since recombination is not yet implemented */
+                if (prevBoundElectrons > newBoundElectrons)
+                {
+                    /* determine number of new macro electrons to be created */
+                    newMacroElectrons = static_cast<unsigned int>(prevBoundElectrons - newBoundElectrons);
+                }
             }
 
     };

src/picongpu/include/particles/ionization/byCollision/ThomasFermi/ThomasFermi_Impl.hpp

@@ -35,8 +35,6 @@ namespace particles
 namespace ionization
 {
 
-    struct AlgorithmNone;
-
     struct AlgorithmThomasFermi;
 
 } // namespace ionization

src/picongpu/include/particles/ionization/byCollision/collisionalIonizationCalc.def

@@ -59,4 +59,34 @@ namespace picongpu
         6.665,
         6.665
         );
-}
+
+namespace particles
+{
+namespace ionization
+{
+
+    /* Fitting parameters to average ionization degree Z* = 4/3*pi*R_0^3 * n(R_0)
+     * as an extension towards arbitrary atoms and temperatures
+     *
+     * See table IV of
+     * \url http://www.sciencedirect.com/science/article/pii/S0065219908601451
+     * doi:10.1016/S0065-2199(08)60145-1
+     */
+    float_X constexpr TFAlpha = 14.3139;
+    float_X constexpr TFBeta  =  0.6624;
+
+    float_X constexpr TFA1    =  3.323e-3;
+    float_X constexpr TFA2    =  9.718e-1;
+    float_X constexpr TFA3    =  9.26148e-5;
+    float_X constexpr TFA4    =  3.10165;
+
+    float_X constexpr TFB0    = -1.7630;
+    float_X constexpr TFB1    =  1.43175;
+    float_X constexpr TFB2    =  0.31546;
+
+    float_X constexpr TFC1    = -0.366667;
+    float_X constexpr TFC2    =  0.983333;
+
+} // namespace ionization
+} // namespace particles
+} // namespace picongpu

src/picongpu/include/simulation_defines/param/ionizerConfig.param

@@ -63,6 +63,15 @@ namespace picongpu
          * 150 attoseconds (classical electron orbit time in hydrogen)  / 2 PI
          */
         constexpr float_64 ATOMIC_UNIT_TIME = 2.4189e-17;
+
+        /** Avogadro number
+         * unit: mol^-1
+         *
+         * Y. Azuma et al. Improved measurement results for the Avogadro
+         * constant using a 28-Si-enriched crystal, Metrologie 52, 2015, 360-375
+         * doi:10.1088/0026-1394/52/2/360
+         */
+        constexpr float_64 N_AVOGADRO = 6.02214076e23;
     }
 
     // converts