Fix: openPMD Weighting Loading (#3783)

* Fix: openPMD Weighting Loading

This fixes the loading of openPMD beams with per-particle
weighting. The previous implementation forgot to `flush()` the
data of the weighting attribute and assumed the weighting of the
first particle applies to all particles of the loaded beam.

This change loads the weighting for each particle.

* Updating the documentation: removed re-scaling feature while using the "external file" injection style.
* Issue a warning if a nonzero q_tot is found in the inputs file
* Formatting and comments

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

1 files changed, 13 insertions, 19 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 51f7569a6..af50b3cf7 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -586,7 +586,7 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot,
         openPMD::ParticleSpecies ps = it.particles.begin()->second;
 
         auto const npart = ps["position"]["x"].getExtent()[0];
-#if !defined(WARPX_DIM_1D_Z)
+#if !defined(WARPX_DIM_1D_Z)  // 2D, 3D, and RZ
         std::shared_ptr<ParticleReal> ptr_x = ps["position"]["x"].loadChunk<ParticleReal>();
         double const position_unit_x = ps["position"]["x"].unitSI();
 #endif
@@ -596,6 +596,9 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot,
         double const momentum_unit_x = ps["momentum"]["x"].unitSI();
         std::shared_ptr<ParticleReal> ptr_uz = ps["momentum"]["z"].loadChunk<ParticleReal>();
         double const momentum_unit_z = ps["momentum"]["z"].unitSI();
+        std::shared_ptr<ParticleReal> ptr_w = ps["weighting"][openPMD::RecordComponent::SCALAR].loadChunk<ParticleReal>();
+        double const w_unit = ps["weighting"][openPMD::RecordComponent::SCALAR].unitSI();
+
 #   if !(defined(WARPX_DIM_XZ) || defined(WARPX_DIM_1D_Z))
         std::shared_ptr<ParticleReal> ptr_y = ps["position"]["y"].loadChunk<ParticleReal>();
         double const position_unit_y = ps["position"]["y"].unitSI();
@@ -604,32 +607,23 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot,
         double momentum_unit_y = 1.0;
         if (ps["momentum"].contains("y")) {
             ptr_uy = ps["momentum"]["y"].loadChunk<ParticleReal>();
-             momentum_unit_y = ps["momentum"]["y"].unitSI();
+            momentum_unit_y = ps["momentum"]["y"].unitSI();
         }
         series->flush();  // shared_ptr data can be read now
 
-        ParticleReal weight = 1.0_prt;  // base standard: no info means "real" particles
         if (q_tot != 0.0) {
-            weight = std::abs(q_tot) / ( std::abs(charge) * ParticleReal(npart) );
-            if (ps.contains("weighting")) {
-                std::stringstream ss;
-                ss << "Both '" << ps_name << ".q_tot' and '"
-                        << ps_name << ".injection_file' specify a total charge.\n'"
-                        << ps_name << ".q_tot' will take precedence.";
-                ablastr::warn_manager::WMRecordWarning("Species", ss.str());
-            }
-        }
-        // ED-PIC extension?
-        else if (ps.contains("weighting")) {
-            // TODO: Add ASSERT_WITH_MESSAGE to test if weighting is a constant record
-            // TODO: Add ASSERT_WITH_MESSAGE for macroWeighted value in ED-PIC
-            ParticleReal w = ps["weighting"][openPMD::RecordComponent::SCALAR].loadChunk<ParticleReal>().get()[0];
-            double const w_unit = ps["weighting"][openPMD::RecordComponent::SCALAR].unitSI();
-            weight = w * w_unit;
+            std::stringstream warnMsg;
+            warnMsg << " Loading particle species from file. " << ps_name << ".q_tot is ignored.";
+            ablastr::warn_manager::WMRecordWarning("AddPlasmaFromFile",
+               warnMsg.str(), ablastr::warn_manager::WarnPriority::high);
         }
 
         for (auto i = decltype(npart){0}; i<npart; ++i){
+
+            ParticleReal const weight = ptr_w.get()[i]*w_unit;
+
 #if !defined(WARPX_DIM_1D_Z)
+
             ParticleReal const x = ptr_x.get()[i]*position_unit_x;
 #else
             ParticleReal const x = 0.0_prt;