Option to specify the type of physical species (#746)

* first implementation, doesnt link

* fix species type and use it in some examples

* eol

* typo in input file

* Apply suggestions from code review

Co-Authored-By: Luca Fedeli <[email protected]>
Co-Authored-By: Yinjian Zhao <[email protected]>

* changes suggested by review

* put species functions into a namespace

Co-authored-by: Luca Fedeli <[email protected]>
Co-authored-by: Yinjian Zhao <[email protected]>

Author: 3 people

Docs/source/running_cpp/parameters.rst

@@ -210,11 +210,17 @@ Particle initialization
     particles are pushed in a standard way, using the specified pusher.
     (see the parameter ``<species_name>.zinject_plane`` below)
 
-* ``<species_name>.charge`` (`float`)
+* ``<species_name>.species_type`` (`string`) optional (default `unspecified`)
+    Type of physical species, ``"electron"``, ``"positron"``, ``"photon"``, ``"hydrogen"``.
+    Either this 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.
+    If ``species_type`` is specified, the charge will be set to the physical value and ``charge`` is optional.
 
-* ``<species_name>.mass`` (`float`)
+* ``<species_name>.mass`` (`float`) optional (default `NaN`)
     The mass of one `physical` particle of this species.
+    If ``species_type`` is specified, the mass will be set to the physical value and ``mass`` is optional.
 
 * ``<species_name>.injection_style`` (`string`)
     Determines how the particles will be injected in the simulation.

Examples/Physics_applications/laser_acceleration/inputs_2d

@@ -39,8 +39,7 @@ warpx.moving_window_v = 1.0 # units of speed of light
 particles.nspecies = 2 # number of species
 particles.species_names = electrons beam
 
-electrons.charge = -q_e
-electrons.mass = m_e
+electrons.species_type = electron
 electrons.injection_style = "NUniformPerCell"
 electrons.num_particles_per_cell_each_dim = 1 1 1
 electrons.xmin = -20.e-6
@@ -51,8 +50,7 @@ electrons.density = 2.e23  # number of electrons per m^3
 electrons.momentum_distribution_type = "constant"
 electrons.do_continuous_injection = 1
 
-beam.charge = -q_e
-beam.mass = m_e
+beam.species_type = electron
 beam.injection_style = "gaussian_beam"
 beam.x_rms = .5e-6
 beam.y_rms = .5e-6

Examples/Physics_applications/plasma_acceleration/inputs_2d_boost

@@ -51,8 +51,7 @@ particles.species_names = driver plasma_e plasma_p beam
 particles.use_fdtd_nci_corr = 1
 particles.rigid_injected_species = driver beam
 
-driver.charge = -q_e
-driver.mass = m_e
+driver.species_type = electron
 driver.injection_style = "gaussian_beam"
 driver.x_rms = 2.e-6
 driver.y_rms = 2.e-6
@@ -74,8 +73,7 @@ driver.rigid_advance = true
 driver.projected = true
 driver.focused = false
 
-plasma_e.charge = -q_e
-plasma_e.mass = m_e
+plasma_e.species_type = electron
 plasma_e.injection_style = "NUniformPerCell"
 plasma_e.zmin = 0.e-6
 plasma_e.zmax = 0.2
@@ -90,8 +88,7 @@ plasma_e.uy = 0.0
 plasma_e.uz = 0.0
 plasma_e.do_continuous_injection = 1
 
-plasma_p.charge = q_e
-plasma_p.mass = m_p
+plasma_p.species_type = hydrogen
 plasma_p.injection_style = "NUniformPerCell"
 plasma_p.zmin = 0.e-6
 plasma_p.zmax = 0.2
@@ -106,8 +103,7 @@ plasma_p.uy = 0.0
 plasma_p.uz = 0.0
 plasma_p.do_continuous_injection = 1
 
-beam.charge = -q_e
-beam.mass = m_e
+beam.species_type = electron
 beam.injection_style = "gaussian_beam"
 beam.x_rms = .5e-6
 beam.y_rms = .5e-6

Examples/Physics_applications/uniform_plasma/inputs_3d

@@ -25,8 +25,7 @@ amr.plot_int = 10
 particles.nspecies = 1
 particles.species_names = electrons
 
-electrons.charge = -q_e
-electrons.mass = m_e
+electrons.species_type = electron
 electrons.injection_style = "NUniformPerCell"
 electrons.num_particles_per_cell_each_dim = 2 2 2
 electrons.profile = constant

Source/Initialization/PlasmaInjector.H

@@ -9,6 +9,7 @@
 #ifndef PLASMA_INJECTOR_H_
 #define PLASMA_INJECTOR_H_
 
+#include "SpeciesPhysicalProperties.H"
 #include "InjectorPosition.H"
 #include "InjectorDensity.H"
 #include "InjectorMomentum.H"
@@ -47,6 +48,7 @@ public:
 
     amrex::Real getCharge () {return charge;}
     amrex::Real getMass () {return mass;}
+    PhysicalSpecies getPhysicalSpecies() const {return physical_species;}
 
     bool doInjection () const noexcept { return inj_pos != NULL;}
 
@@ -87,6 +89,8 @@ protected:
 
     amrex::Real mass, charge;
 
+    PhysicalSpecies physical_species;
+
     amrex::Real density;
 
     int species_id;

Source/Initialization/PlasmaInjector.cpp

@@ -8,7 +8,7 @@
  * License: BSD-3-Clause-LBNL
  */
 #include "PlasmaInjector.H"
-
+#include "SpeciesPhysicalProperties.H"
 #include "Utils/WarpXConst.H"
 #include "Utils/WarpXUtil.H"
 #include "WarpX.H"
@@ -95,22 +95,50 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
     pp.query("density_min", density_min);
     pp.query("density_max", density_max);
 
+    std::string physical_species_s;
+    bool species_is_specified = pp.query("species_type", physical_species_s);
+    if (species_is_specified){
+        physical_species = species::from_string( physical_species_s );
+        // charge = SpeciesCharge[physical_species];
+        charge = species::get_charge( physical_species );
+        // mass = SpeciesMass[physical_species];
+        mass = species::get_mass( physical_species );
+    }
+
     // parse charge and mass
     std::string charge_s;
-    pp.get("charge", charge_s);
-    std::transform(charge_s.begin(),
-                   charge_s.end(),
-                   charge_s.begin(),
-                   ::tolower);
-    charge = parseChargeString(pp, charge_s);
+    bool charge_is_specified = pp.query("charge", charge_s);
+    if (charge_is_specified){
+        std::transform(charge_s.begin(),
+                       charge_s.end(),
+                       charge_s.begin(),
+                       ::tolower);
+        charge = parseChargeString(pp, charge_s);
+    }
+    if ( charge_is_specified && species_is_specified ){
+        Print()<<"WARNING: Both <species>.charge and <species>species_type specified\n";
+        Print()<<"         The charge in <species>.mass overwrite the one from <species>.species_type\n";
+    }
+    if (!charge_is_specified && !species_is_specified){
+        amrex::Abort("Need to specify at least one of species_type or charge");
+    }
 
     std::string mass_s;
-    pp.get("mass", mass_s);
-    std::transform(mass_s.begin(),
-                   mass_s.end(),
-                   mass_s.begin(),
-                   ::tolower);
-    mass = parseMassString(pp, mass_s);
+    bool mass_is_specified = pp.query("mass", mass_s);
+    if (mass_is_specified){
+        std::transform(mass_s.begin(),
+                       mass_s.end(),
+                       mass_s.begin(),
+                       ::tolower);
+        mass = parseMassString(pp, mass_s);
+    }
+    if ( mass_is_specified && species_is_specified ){
+        Print()<<"WARNING: Both <species>.mass and <species>species_type specified\n";
+        Print()<<"         The mass in <species>.mass overwrite the one from <species>.species_type\n";
+    }
+    if (!mass_is_specified && !species_is_specified){
+        amrex::Abort("Need to specify at least one of species_type or mass");
+    }
 
     // parse injection style
     std::string part_pos_s;
@@ -362,4 +390,3 @@ PlasmaInjector::getInjectorMomentum ()
 {
     return inj_mom.get();
 }
-

Source/Particles/PhysicalParticleContainer.cpp

@@ -40,6 +40,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
       species_name(name)
 {
     plasma_injector.reset(new PlasmaInjector(species_id, species_name));
+    physical_species = plasma_injector->getPhysicalSpecies();
     charge = plasma_injector->getCharge();
     mass = plasma_injector->getMass();
 

Source/Particles/SpeciesPhysicalProperties.H

@@ -0,0 +1,78 @@
+/* Copyright 2020 Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_SPECIESPHYSICALPROPERTIES_H_
+#define WARPX_SPECIESPHYSICALPROPERTIES_H_
+
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_REAL.H>
+#include <AMReX_AmrCore.H>
+
+#include <map>
+#include <limits>
+
+enum struct PhysicalSpecies{unspecified=0, electron, positron, photon, hydrogen};
+
+namespace species
+{
+    static 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;
+        amrex::Abort("unknown PhysicalSpecies");
+        return PhysicalSpecies::unspecified;
+    }
+
+    static 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;
+        default:
+            amrex::Abort("unknown PhysicalSpecies");
+            return 0.;
+        }
+    }
+
+    static 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;
+        default:
+            amrex::Abort("unknown PhysicalSpecies");
+            return 0.;
+        }
+    }
+}
+
+#endif // WARPX_SPECIESPHYSICALPROPERTIES_H_

Source/Particles/WarpXParticleContainer.H

@@ -10,6 +10,8 @@
 #ifndef WARPX_WarpXParticleContainer_H_
 #define WARPX_WarpXParticleContainer_H_
 
+#include "Utils/WarpXConst.H"
+#include "SpeciesPhysicalProperties.H"
 #include "Evolve/WarpXDtType.H"
 
 #include <AMReX_Particles.H>
@@ -319,6 +321,7 @@ protected:
 
     amrex::Real charge;
     amrex::Real mass;
+    PhysicalSpecies physical_species;
 
     //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid
     bool m_deposit_on_main_grid = false;