Read species distribution from OPMD - part 2 (#859)

* Added <species>.profile=external_file and .profile_file

* Added description of input parameters to Docs

* Changed from profile to injection option for external file

* Fix typo in amrex abort message (due to copy paste)

* Added the OpenPMD use amrex abort message

* Minor fix - not sure how to remove EOL issue

* Tried to add AddExternalFileBeam functon to PhysicalParticleContainer

* Trued to fix EOL white space issue

* Added read/print species name from OPMD file

* Fixed OpenPMD charge and mass read

* Added number of particles in species

* Added nparts and converted charge/mass units to SI

* Fix to nr of particles

* Added q_tot parameter to determine part weight

* Added macroparticle's weight

* Fixed std::int typo - use only int

* Update Source/Initialization/PlasmaInjector.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* No need to include openPMD header yet

* Minor Fixes

* Fix EOL according to @ax3l's recommendation in #845

* Remove commented out AbortMessage

Co-Authored-By: MaxThevenet <[email protected]>

* Removed commented out part initialization (used only in branch for next PR)

Co-Authored-By: MaxThevenet <[email protected]>

* Added warning that this is WIP

Co-Authored-By: MaxThevenet <[email protected]>

* Changed function name to AddPlasmaFromFile

* Removed AMReX warning from loop

Co-Authored-By: MaxThevenet <[email protected]>

* Added fix suggested in PR# 847

* Changed to mevpc2_kg

* Fix EOL again

* Style fix adding spaces around ` = `

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Style fix added spaces in definitions and between comma

Co-Authored-By: Axel Huebl <[email protected]>

* Style fix

Co-Authored-By: Axel Huebl <[email protected]>

* Added Assert Message for #iterations and TODO const mass/charge

* Fix EOL again

* Update Docs/source/running_cpp/parameters.rst

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Corrected division by int

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <[email protected]>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-authored-by: Axel Huebl <[email protected]>
Co-authored-by: MaxThevenet <[email protected]>

Author: 3 people

Docs/source/running_cpp/parameters.rst

@@ -280,7 +280,7 @@ Particle initialization
       symmetrize the beam in the x and y directions).
 
     * ``external_file``: inject macroparticles with properties (charge, mass, position, and momentum) according to data in external file.
-      It requires the additional argument ``<species_name>.injection_file``, which is the string corresponding to the OpenPMD file name.
+      It requires the additional arguments ``<species_name>.injection_file`` and ``<species_name>.q_tot``, which are the string corresponding to the openPMD file name and the beam charge.
       When using this style, it is not necessary to add other ``<species_name>.(...)`` paramters, because they will be read directly from the file.
 
 * ``<species_name>.num_particles_per_cell_each_dim`` (`3 integers in 3D and RZ, 2 integers in 2D`)

Source/Initialization/PlasmaInjector.cpp

@@ -229,9 +229,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
 #ifdef WARPX_USE_OPENPMD
         external_file = true;
         pp.get("injection_file",str_injection_file);
+        pp.get("q_tot",q_tot);
 #else
-        amrex::Abort("WarpX has to be compiled with USE_OPENPMD=TRUE to be able to"
-                     " read the external openPMD file with species data");
+        amrex::Abort("WarpX has to be compiled with USE_OPENPMD=TRUE to be able"
+                     " to read the external openPMD file with species data");
 #endif
     } else {
         StringParseAbortMessage("Injection style", part_pos_s);

Source/Particles/PhysicalParticleContainer.H

@@ -206,8 +206,9 @@ public:
     /** Load a particle beam from an external file
      *
      * @param[in] s_f path to an external file in openPMD markup
+     * @param[in] q_tot total charge of the particle species to be initialized
      */
-    void AddPlasmaFromFile(const std::string s_f);
+    void AddPlasmaFromFile(const std::string s_f, amrex::Real q_tot);
 
     void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z,
                              std::array<amrex::Real, 3> u,

Source/Particles/PhysicalParticleContainer.cpp

@@ -307,17 +307,34 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
 }
 
 void
-PhysicalParticleContainer::AddPlasmaFromFile(const std::string s_f)
+PhysicalParticleContainer::AddPlasmaFromFile(const std::string s_f, amrex::Real q_tot)
 {
 #ifdef WARPX_USE_OPENPMD
     openPMD::Series series = openPMD::Series(s_f, openPMD::AccessType::READ_ONLY);
     amrex::Print() << "openPMD standard version " << series.openPMD() << "\n";
 
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(series.iterations.size() == 1u, "External "
+                                     "file should contain only one iteration\n");
     openPMD::Iteration& i = series.iterations[1];
-    amrex::Print()  << "File contains " << i.particles.size() << " specie(s):" << "\n";
-    for( auto const& ps : i.particles ) {
-        amrex::Print() << "\t" << ps.first << "\n";
-    }
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(i.particles.size() == 1u, "External file "
+                                     "should contain only one species\n");
+    std::pair<std::string,openPMD::ParticleSpecies> ps = *i.particles.begin();
+
+    //TODO: In future PRs will add AMREX_ALWAYS_ASSERT_WITH_MESSAGE to test if mass and charge are both const
+    amrex::Real p_m = ps.second["mass"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::Real>().get()[0];
+    amrex::Real p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::Real>().get()[0];
+    int npart = ps.second["position"]["x"].getExtent()[0];
+    series.flush();
+
+    mass = p_m*PhysConst::mevpc2_kg;
+    charge = p_q*PhysConst::q_e;
+    Real const weight = q_tot/(charge*amrex::Real(npart));
+
+    amrex::Print() << npart << " parts of species " << ps.first << "\nWith"
+    << " mass = " << mass << " and charge = " << charge << "\nTo initialize "
+    << npart << " macroparticles with weight " << weight << "\n";
+
     amrex::Print()<<"WARNING: this is WIP, no particle has been injected!!";
 #endif
     return;
@@ -380,7 +397,8 @@ PhysicalParticleContainer::AddParticles (int lev)
     }
 
     if (plasma_injector->external_file) {
-        AddPlasmaFromFile(plasma_injector->str_injection_file);
+        AddPlasmaFromFile(plasma_injector->str_injection_file,
+                          plasma_injector->q_tot);
         return;
     }
 

Source/Utils/WarpXConst.H

@@ -34,6 +34,7 @@ namespace PhysConst
                                          (45.*m_e*m_e*m_e*m_e*c*c*c*c*c); // SI value is 1.3050122.e-52
     static constexpr amrex::Real xi_c2 = xi * c * c; // This should be usable for single precision, though
     // very close to smallest number possible: smallest number = 1.2e-38, xi_c2 = 1.1e-35
+    static constexpr amrex::Real mevpc2_kg = 1.7826619216279e-30; // to convert MeV in kg
 }
 
 #endif