Read species distribution from OPMD - part 3 (#883)

* 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

* Added read x,y and z + ifdef 3D

* Added velocity of particles and stored in container

* Converted velocities to SI

* No need to specfy momentum distribution if external file is used

* Update Source/Initialization/PlasmaInjector.cpp

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.H

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/Particles/PhysicalParticleContainer.cpp

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* No need to include openPMD header yet

* Fixed const in AddExternalFileBeam

* Fix compatibility with read_opmd and read_opmd_2

* Corrected position and charge units/sign

* Added Doc note and abort message for RZ

* Fix typo and EOL

* Minor fixes

* Fixed details. Added fix to EOL - testing

* Changed to physical_charge and explained it in Docs

* Fix header extra openPMD include files

* Removed additional debugging comments

* Try doxygen again

* Fix plasmainjector physical_q_tot

* Change apart to auto (not long as in gaussian beam injection style)

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Correct for new type of npart

Co-Authored-By: Axel Huebl <axel.huebl@plasma.ninja>

* Removed Geant4 renormalization

* Read openPMD file and checked its units

* Trying to correct momentum information in plotfiles

* Compilation error - > SegFault 11

* Added unitSI() to each position/momentum direction. Compilation SegFault

* Re-structured code to use only once series.flush()

* Commented out amrex::Print() s

* Path to fix issue

* Update after review

* Chanegd <physical_q_tot> to <q_tot> and made it optional

* Fixed documentation typo and re-organized it

* Added const to npart as per reviewer recommendation

* Fixed issue with velocity - which became momentum

* Fixed 2D and 3D options + documentation

* Implemented fixes from past reviewsand fixed duplicate entry in Docs

* Fix auto in iterator in for loop

* Added fixes according to reviewer

* Revert "Added fixes according to reviewer"

This reverts commit 0a485d83c28014c8b6f53a30f26719f23f89bad5.

* Fixed IO block and reordered logic with Axel

* Fix spaces

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

1 files changed, 93 insertions, 28 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 0bfbeea16..59b6c5ca3 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -254,7 +254,7 @@ PhysicalParticleContainer::AddGaussianBeam (
     std::normal_distribution<double> disty(y_m, y_rms);
     std::normal_distribution<double> distz(z_m, z_rms);
 
-    // Allocate temporary vectors on the CPU
+    // Declare temporary vectors on the CPU
     Gpu::HostVector<ParticleReal> particle_x;
     Gpu::HostVector<ParticleReal> particle_y;
     Gpu::HostVector<ParticleReal> particle_z;
@@ -326,36 +326,101 @@ PhysicalParticleContainer::AddGaussianBeam (
 }
 
 void
-PhysicalParticleContainer::AddPlasmaFromFile (const std::string s_f, amrex::Real q_tot)
+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_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));
+    // Declare temporary vectors on the CPU
+    Gpu::HostVector<ParticleReal> particle_x;
+    Gpu::HostVector<ParticleReal> particle_z;
+    Gpu::HostVector<ParticleReal> particle_ux;
+    Gpu::HostVector<ParticleReal> particle_uz;
+    Gpu::HostVector<ParticleReal> particle_w;
+    Gpu::HostVector<ParticleReal> particle_y;
+    Gpu::HostVector<ParticleReal> particle_uy;
 
-    amrex::Print() << npart << " parts of species " << ps.first << "\nWith"
-    << " mass = " << mass << " and charge = " << charge << "\nTo initialize "
-    << npart << " macroparticles with weight " << weight << "\n";
+#ifdef WARPX_USE_OPENPMD
+    //TODO: Make changes for read/write in multiple MPI ranks
+    if (ParallelDescriptor::IOProcessor()) {
+        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 1 iteration\n");
+        openPMD::Iteration& i = series.iterations[1];
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(i.particles.size() == 1u, "External file "
+                                         "should contain only 1 species\n");
+        std::pair<std::string,openPMD::ParticleSpecies> ps = *i.particles.begin();
+
+        //TODO: Add ASSERT_WITH_MESSAGE to test if mass and charge are both const
+        amrex::ParticleReal p_m = ps.second["mass"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
+        double const mass_unit = ps.second["mass"][openPMD::RecordComponent::SCALAR].unitSI();
+        amrex::ParticleReal p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
+        double const charge_unit = ps.second["charge"][openPMD::RecordComponent::SCALAR].unitSI();
+#   if (defined WARPX_DIM_3D) || (defined WARPX_DIM_2D)
+        auto const npart = ps.second["position"]["x"].getExtent()[0];
+        std::shared_ptr<amrex::ParticleReal> ptr_x = ps.second["position"]["x"].loadChunk<amrex::ParticleReal>();
+        double const position_unit_x = ps.second["position"]["x"].unitSI();
+        std::shared_ptr<amrex::ParticleReal> ptr_z = ps.second["position"]["z"].loadChunk<amrex::ParticleReal>();
+        double const position_unit_z = ps.second["position"]["z"].unitSI();
+        std::shared_ptr<amrex::ParticleReal> ptr_ux = ps.second["momentum"]["x"].loadChunk<amrex::ParticleReal>();
+        double const momentum_unit_x = ps.second["momentum"]["x"].unitSI();
+        std::shared_ptr<amrex::ParticleReal> ptr_uz = ps.second["momentum"]["z"].loadChunk<amrex::ParticleReal>();
+        double const momentum_unit_z = ps.second["momentum"]["z"].unitSI();
+#   else
+        amrex::Abort("AddPlasmaFromFile is only implemented for 2D and 3D\n")
+#   endif
+#   if (defined WARPX_DIM_3D)
+        std::shared_ptr<amrex::ParticleReal> ptr_y = ps.second["position"]["y"].loadChunk<amrex::ParticleReal>();
+        double const position_unit_y = ps.second["position"]["y"].unitSI();
+        std::shared_ptr<amrex::ParticleReal> ptr_uy = ps.second["momentum"]["y"].loadChunk<amrex::ParticleReal>();
+        double const momentum_unit_y = ps.second["momentum"]["y"].unitSI();
+#   endif
+        series.flush();
+
+        mass=p_m*mass_unit;
+        charge=p_q*charge_unit;
+
+        amrex::Real weight;
+        if (q_tot != 0.0){
+            weight = q_tot/(p_q*amrex::Real(npart));
+        }
+        else {
+            weight = charge;
+        }
 
-    amrex::Print()<<"WARNING: this is WIP, no particle has been injected!!";
-#endif
+        for (auto i = decltype(npart){0}; i<npart; ++i){
+            amrex::ParticleReal const x = ptr_x.get()[i]*position_unit_x;
+            amrex::ParticleReal const z = ptr_z.get()[i]*position_unit_z;
+#   if (defined WARPX_DIM_2D)
+            amrex::Real const y = 0.0;
+#   elif (defined WARPX_DIM_3D)
+            amrex::ParticleReal const y = ptr_y.get()[i]*position_unit_y;
+#   endif
+            if (plasma_injector->insideBounds(x, y, z)) {
+                amrex::ParticleReal const ux = ptr_ux.get()[i]*momentum_unit_x/PhysConst::m_e;
+                amrex::ParticleReal const uz = ptr_uz.get()[i]*momentum_unit_z/PhysConst::m_e;
+#   if (defined WARPX_DIM_2D)
+                amrex::ParticleReal const uy = 0.0;
+#   elif (defined WARPX_DIM_3D)
+                amrex::ParticleReal const uy = ptr_uy.get()[i]*momentum_unit_y/PhysConst::m_e;
+#   endif
+                CheckAndAddParticle(x, y, z, { ux, uy, uz}, weight,
+                particle_x,  particle_y,  particle_z,
+                particle_ux, particle_uy, particle_uz,
+                particle_w);
+            }
+        }
+        auto const np = particle_z.size();
+        if (np < npart) {
+            amrex::Print()<<"WARNING: Simulation box doesn't cover all particles\n";
+        }
+    } //IO Processor
+    auto const np = particle_z.size();
+    AddNParticles(0,np,
+                  particle_x.dataPtr(),  particle_y.dataPtr(),  particle_z.dataPtr(),
+                  particle_ux.dataPtr(), particle_uy.dataPtr(), particle_uz.dataPtr(),
+                  1, particle_w.dataPtr(),1);
+#endif //OPENPMD
     return;
 }