Add new entries in Physical Species properties (+ refactoring of the SpeciesPhysicalProperties file)

Docs/source/usage/parameters.rst

@@ -570,8 +570,19 @@ Particle initialization
 
 * ``<species_name>.species_type`` (`string`) optional (default `unspecified`)
     Type of physical species.
-    Currently, the accepted species are ``"electron"``, ``"positron"``, ``"photon"``, ``"hydrogen"`` (or equivalently ``"proton"``), ``"helium"`` (or equivalently ``"alpha"``), ``"boron"``, ``"carbon"``, ``"oxygen"``, ``"nitrogen"``, ``"argon"``, ``"copper"`` and ``"xenon"``.
-    Either this or both ``mass`` and ``charge`` have to be specified.
+    Currently, the accepted species are
+    ``"electron"``, ``"positron"``, ``"muon"``, ``"antimuon"``, ``"photon"``, ``"proton"`` , ``"alpha"``,
+    ``"hydrogen"``, ``"protium"``, ``"deuterium"``, ``"tritium"``, ``"helium"``, ``"helium3"``, ``"helium4"``,
+    ``"lithium"``, ``"lithium6"``, ``"lithium7"``, ``"beryllium"``, ``"boron"``, ``"boron10"``, ``"boron11"``,
+    ``"carbon"``, ``"carbon12"``, ``"carbon13"``, ``"nitrogen"``, ``"nitrogen14"``, ``"nitrogen15"``,
+    ``"oxygen"``, ``"oxygen16"``, ``"oxygen17"``, ``"oxygen18"``, ``"fluorine"``, ``"neon"``, ``"neon20"``,
+    ``"neon21"``, ``"neon22"``, ``"aluminium"``, ``"argon"``, ``"copper"``, ``"xenon"`` and ``"gold"``.
+    The difference between ``"proton"`` and ``"protium"`` is that the mass of the latter includes also the mass
+    of the bound electron (same for ``"alpha"`` and ``"helium4"``). When just the name of an element is specified, the mass
+    is a weighted average of the masses of the stable isotopes. For all the elements with ``Z < 11`` we provide
+    also the stable isotopes as an option for ``species_type`` (e.g., ``"helium3"`` and ``"helium4"``), as well as all the isotopes
+    of ``"hydrogen"`` (``"protium"``, ``"deuterium"``, ``"tritium"``).
+    Either ``species_type`` or both ``mass`` and ``charge`` have to be specified.
 
 * ``<species_name>.charge`` (`float`) optional (default `NaN`)
     The charge of one `physical` particle of this species.

Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.cpp

@@ -5,9 +5,12 @@
  * License: BSD-3-Clause-LBNL
  */
 #include "BackTransformParticleFunctor.H"
+
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/WarpXParticleContainer.H"
+#include "Utils/WarpXConst.H"
 #include "WarpX.H"
+
 #include <AMReX.H>
 #include <AMReX_Print.H>
 #include <AMReX_BaseFwd.H>

Source/Initialization/PlasmaInjector.cpp

@@ -125,7 +125,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
     std::string physical_species_s;
     bool species_is_specified = pp_species_name.query("species_type", physical_species_s);
     if (species_is_specified){
-        physical_species = species::from_string( physical_species_s );
+        const auto physical_species_from_string = species::from_string( physical_species_s );
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(physical_species_from_string,
+            physical_species_s + " does not exist!");
+        physical_species = physical_species_from_string.value();
         charge = species::get_charge( physical_species );
         mass = species::get_mass( physical_species );
     }

Source/Particles/CMakeLists.txt

@@ -8,6 +8,7 @@ target_sources(WarpX
     WarpXParticleContainer.cpp
     LaserParticleContainer.cpp
     ParticleBoundaryBuffer.cpp
+    SpeciesPhysicalProperties.cpp
 )
 
 #add_subdirectory(Algorithms)

Source/Particles/Deposition/CurrentDeposition.H

@@ -12,6 +12,7 @@
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/ShapeFactors.H"
 #include "Utils/WarpXAlgorithmSelection.H"
+#include "Utils/WarpXConst.H"
 #ifdef WARPX_DIM_RZ
 #   include "Utils/WarpX_Complex.H"
 #endif

Source/Particles/Gather/GetExternalFields.H

@@ -4,6 +4,7 @@
 #include "Particles/Pusher/GetAndSetPosition.H"
 
 #include "Particles/WarpXParticleContainer_fwd.H"
+#include "Utils/WarpXConst.H"
 
 #include <AMReX.H>
 #include <AMReX_Array.H>

Source/Particles/Make.package

@@ -6,6 +6,7 @@ CEXE_sources += PhotonParticleContainer.cpp
 CEXE_sources += LaserParticleContainer.cpp
 CEXE_sources += ParticleBoundaryBuffer.cpp
 CEXE_sources += ParticleBoundaries.cpp
+CEXE_sources += SpeciesPhysicalProperties.cpp
 
 include $(WARPX_HOME)/Source/Particles/Algorithms/Make.package
 include $(WARPX_HOME)/Source/Particles/Pusher/Make.package

Source/Particles/SpeciesPhysicalProperties.H

@@ -8,180 +8,50 @@
 #ifndef WARPX_SPECIESPHYSICALPROPERTIES_H_
 #define WARPX_SPECIESPHYSICALPROPERTIES_H_
 
-#include "Utils/WarpXConst.H"
-
-#include <AMReX_AmrCore.H>
 #include <AMReX_REAL.H>
 
-#include <limits>
-#include <map>
+#include <optional>
 #include <string>
 
-enum struct PhysicalSpecies{unspecified=0, electron, positron, photon, hydrogen, helium, boron,
-                            boron10, boron11, carbon, nitrogen, oxygen, argon, copper, xenon};
+enum struct PhysicalSpecies{
+    unspecified=0, electron, positron, muon, antimuon, photon, hydrogen, protium, deuterium, tritium, proton,
+    helium, helium3, helium4, alpha, lithium, lithium6, lithium7, beryllium, boron, boron10, boron11, carbon,
+    carbon12, carbon13, nitrogen, nitrogen14, nitrogen15, oxygen, oxygen16, oxygen17, oxygen18, fluorine,
+    neon, neon20, neon21, neon22, aluminium, argon, copper, xenon, gold};
 
 namespace species
 {
-    AMREX_FORCE_INLINE
-    PhysicalSpecies from_string(std::string species)
-    {
-        if( species=="unspecified" )
-            return PhysicalSpecies::unspecified;
-        if( species=="electron" )
-            return PhysicalSpecies::electron;
-        if( species=="positron" )
-            return PhysicalSpecies::positron;
-        if( species=="photon" )
-            return PhysicalSpecies::photon;
-        if( species=="hydrogen" )
-            return PhysicalSpecies::hydrogen;
-        if( species=="proton" )
-            return PhysicalSpecies::hydrogen;
-        if( species=="helium" )
-            return PhysicalSpecies::helium;
-        if( species=="alpha" )
-            return PhysicalSpecies::helium;
-        if( species=="boron" )
-            return PhysicalSpecies::boron;
-        if( species=="boron10" )
-            return PhysicalSpecies::boron10;
-        if( species=="boron11" )
-            return PhysicalSpecies::boron11;
-        if( species=="carbon" )
-            return PhysicalSpecies::carbon;
-        if( species=="nitrogen" )
-            return PhysicalSpecies::nitrogen;
-        if( species=="oxygen" )
-            return PhysicalSpecies::oxygen;
-        if( species=="argon" )
-            return PhysicalSpecies::argon;
-        if( species=="copper" )
-            return PhysicalSpecies::copper;
-        if( species=="xenon" )
-            return PhysicalSpecies::xenon;
-        amrex::Abort("unknown PhysicalSpecies");
-        return PhysicalSpecies::unspecified;
-    }
-
-    AMREX_FORCE_INLINE
-    amrex::Real get_charge (PhysicalSpecies ps)
-    {
-        switch(ps) {
-        case PhysicalSpecies::unspecified:
-            return std::numeric_limits<amrex::Real>::quiet_NaN();
-        case PhysicalSpecies::electron:
-            return -PhysConst::q_e;
-        case PhysicalSpecies::positron:
-            return PhysConst::q_e;
-        case PhysicalSpecies::photon:
-            return 0.;
-        case PhysicalSpecies::hydrogen:
-            return PhysConst::q_e;
-        case PhysicalSpecies::helium:
-            return PhysConst::q_e * amrex::Real(2.0);
-        case PhysicalSpecies::boron:
-            return PhysConst::q_e * amrex::Real(5.0);
-        case PhysicalSpecies::boron10:
-            return PhysConst::q_e * amrex::Real(5.0);
-        case PhysicalSpecies::boron11:
-            return PhysConst::q_e * amrex::Real(5.0);
-        case PhysicalSpecies::carbon:
-            return PhysConst::q_e * amrex::Real(6.0);
-        case PhysicalSpecies::nitrogen:
-            return PhysConst::q_e * amrex::Real(7.0);
-        case PhysicalSpecies::oxygen:
-            return PhysConst::q_e * amrex::Real(8.0);
-        case PhysicalSpecies::argon:
-            return PhysConst::q_e * amrex::Real(18.0);
-        case PhysicalSpecies::copper:
-            return PhysConst::q_e * amrex::Real(29.0);
-        case PhysicalSpecies::xenon:
-            return PhysConst::q_e * amrex::Real(54.0);
-        default:
-            amrex::Abort("unknown PhysicalSpecies");
-            return 0.;
-        }
-    }
-
-    AMREX_FORCE_INLINE
-    amrex::Real get_mass (PhysicalSpecies ps)
-    {
-        switch(ps) {
-        case PhysicalSpecies::unspecified:
-            return std::numeric_limits<amrex::Real>::quiet_NaN();
-        case PhysicalSpecies::electron:
-            return PhysConst::m_e;
-        case PhysicalSpecies::positron:
-            return PhysConst::m_e;
-        case PhysicalSpecies::photon:
-            return 0.;
-        case PhysicalSpecies::hydrogen:
-            return PhysConst::m_p;
-        case PhysicalSpecies::helium:
-            return PhysConst::m_p * amrex::Real(3.97369);
-        case PhysicalSpecies::boron:
-            return PhysConst::m_p * amrex::Real(10.7319);
-        case PhysicalSpecies::boron10:
-            return PhysConst::m_p * amrex::Real(9.94060);
-        case PhysicalSpecies::boron11:
-            return PhysConst::m_p * amrex::Real(10.9298);
-        case PhysicalSpecies::carbon:
-            return PhysConst::m_e * amrex::Real(22032.0);
-        case PhysicalSpecies::nitrogen:
-            return PhysConst::m_e * amrex::Real(25716.9);
-        case PhysicalSpecies::oxygen:
-            return PhysConst::m_p * amrex::Real(15.8834);
-        case PhysicalSpecies::argon:
-            return PhysConst::m_p * amrex::Real(39.9480);
-        case PhysicalSpecies::copper:
-            return PhysConst::m_p * amrex::Real(63.0864);
-        case PhysicalSpecies::xenon:
-            return PhysConst::m_p * amrex::Real(131.293);
-        default:
-            amrex::Abort("unknown PhysicalSpecies");
-            return 0.;
-        }
-    }
-
-    AMREX_FORCE_INLINE
-    std::string get_name (PhysicalSpecies ps)
-    {
-        switch(ps) {
-        case PhysicalSpecies::unspecified:
-            return "unspecified";
-        case PhysicalSpecies::electron:
-            return "electron";
-        case PhysicalSpecies::positron:
-            return "positron";
-        case PhysicalSpecies::photon:
-            return "photon";
-        case PhysicalSpecies::hydrogen:
-            return "hydrogen";
-        case PhysicalSpecies::helium:
-            return "helium";
-        case PhysicalSpecies::boron:
-            return "boron";
-        case PhysicalSpecies::boron10:
-            return "boron10";
-        case PhysicalSpecies::boron11:
-            return "boron11";
-        case PhysicalSpecies::carbon:
-            return "carbon";
-        case PhysicalSpecies::nitrogen:
-            return "nitrogen";
-        case PhysicalSpecies::oxygen:
-            return "oxygen";
-        case PhysicalSpecies::argon:
-            return "argon";
-        case PhysicalSpecies::copper:
-            return "copper";
-        case PhysicalSpecies::xenon:
-            return "xenon";
-        default:
-            amrex::Abort("unknown PhysicalSpecies");
-            return "";
-        }
-    }
+    /**
+     * \brief Returns the PhysicalSpecies associated to a given name
+     *
+     * \param[in] species_name the name of a species
+     * \return the PhysicalSpecies corresponding to species_name (if it exists)
+     */
+    std::optional<PhysicalSpecies> from_string (const std::string& species_name);
+
+    /**
+     * \brief Returns the charge associated to a PhysicalSpecies
+     *
+     * \param[in] ps the PhysicalSpecies
+     * \return the charge associated to the PhysicalSpecies
+     */
+    amrex::Real get_charge (const PhysicalSpecies& ps);
+
+    /**
+     * \brief Returns the mass associated to a PhysicalSpecies
+     *
+     * \param[in] ps the PhysicalSpecies
+     * \return the mass associated to the PhysicalSpecies
+     */
+    amrex::Real get_mass (const PhysicalSpecies& ps);
+
+    /**
+     * \brief Returns the name associated to a PhysicalSpecies
+     *
+     * \param[in] ps the PhysicalSpecies
+     * \return the name associated to the PhysicalSpecies
+     */
+    std::string get_name (const PhysicalSpecies& ps);
 }
 
 #endif // WARPX_SPECIESPHYSICALPROPERTIES_H_