openPMD: Handle Zero Particles Well (#2980)

* openPMD: Handle Zero Particles Well

When a time step for output encounters zero particles in a species,
then we still want to dump them as "empty" species in openPMD. That
simplifies post-processing a lot and we have the mechanisms in
openPMD for it :)

* openPMD: Emtpy Particle Writes

Write empty records for iterations (steps or lab steps for BTD,
respectively) without particles in a species.

* Re-order: ED-PIC & Constant Particle Records

- set attributes once
- set constant records once
- clean up into appropriate functions

* Enable BTD for ADIOS :)

Works now as well :tada:

* Fix lingo in comments (Reva)

Thank you!! :)

Co-authored-by: Revathi  Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com>

Co-authored-by: Revathi  Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com>

4 files changed, 262 insertions, 215 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index cdefa2259..383781868 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -1908,7 +1908,6 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
     ``bp`` is the `ADIOS I/O library <https://csmd.ornl.gov/adios>`_, ``h5`` is the `HDF5 format <https://www.hdfgroup.org/solutions/hdf5/>`_, and ``json`` is a `simple text format <https://en.wikipedia.org/wiki/JSON>`_.
     ``json`` only works with serial/single-rank jobs.
     When WarpX is compiled with openPMD support, the first available backend in the order given above is taken.
-    Note that when using ``BackTransformed`` diagnostic type, the openpmd format supports only ``h5`` backend for both species and fields, while ``bp`` backend can be used for visualizing fields, but not particles. The code will abort if ``bp`` is selected for particle output.
 
 * ``<diag_name>.openpmd_encoding`` (optional, ``v`` (variable based), ``f`` (file based) or ``g`` (group based) ) only read if ``<diag_name>.format = openpmd``.
      openPMD `file output encoding <https://openpmd-api.readthedocs.io/en/0.14.0/usage/concepts.html#iteration-and-series>`__.
diff --git a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp
index 951013555..4c31dbecd 100644
--- a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp
+++ b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp
@@ -100,26 +100,6 @@ FlushFormatOpenPMD::FlushFormatOpenPMD (const std::string& diag_name)
     engine_type, engine_parameters,
     warpx.getPMLdirections()
   );
-
-  // Temporarily adding Abort for adios filetype if species is selected for BTD output
-  bool species_output = true;
-  int write_species = 1;
-  std::vector< std::string > output_species_names;
-  bool species_specified = pp_diag_name.queryarr("species", output_species_names);
-  if (species_specified == true and output_species_names.size() > 0) {
-      species_output = true;
-  } else {
-      // By default species output is computed for all diagnostics, if write_species is not set to 0
-      species_output = true;
-  }
-  // Check user-defined option to turn off species output
-  pp_diag_name.query("write_species", write_species);
-  if (write_species == 0) species_output = false;
-  if (diag_type_str == "BackTransformed" and species_output == true) {
-      if (m_OpenPMDPlotWriter->OpenPMDFileType() == "bp") {
-          amrex::Abort(" Currently BackTransformed diagnostics type does not support species output for ADIOS backend. Please select h5 as openpmd backend");
-      }
-  }
 }
 
 void
@@ -155,7 +135,7 @@ FlushFormatOpenPMD::WriteToFile (
         varnames, mf, geom, output_levels, output_iteration, time, isBTD, full_BTD_snapshot);
 
     // particles: all (reside only on locally finest level)
-    m_OpenPMDPlotWriter->WriteOpenPMDParticles(particle_diags, isBTD, totalParticlesFlushedAlready);
+    m_OpenPMDPlotWriter->WriteOpenPMDParticles(particle_diags, isBTD, isLastBTDFlush, totalParticlesFlushedAlready);
 
     // signal that no further updates will be written to this iteration
     m_OpenPMDPlotWriter->CloseStep(isBTD, isLastBTDFlush);
diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H
index 819035360..19d0aec5e 100644
--- a/Source/Diagnostics/WarpXOpenPMD.H
+++ b/Source/Diagnostics/WarpXOpenPMD.H
@@ -135,6 +135,7 @@ public:
   void WriteOpenPMDParticles (
               const amrex::Vector<ParticleDiag>& particle_diags,
               const bool isBTD = false,
+              const bool isLastBTDFlush = false,
               const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>());
 
   /** Write out all openPMD fields for all active MR levels
@@ -209,18 +210,31 @@ private:
       std::string& comp_name
   ) const;
 
-  /** This function sets up the entries for storing the particle positions, global IDs, and constant records (charge, mass)
+  /** This function sets up the entries for storing the particle positions and global IDs
   *
   * @param[in] currSpecies Corresponding openPMD species
   * @param[in] np          Number of particles
-  * @param[in] charge      Charge of the particles (note: fix for ions)
-  * @param[in] mass        Mass of the particles
+  * @param[in] isBTD       Is this a back-transformed diagnostics output?
   */
   void SetupPos (
         openPMD::ParticleSpecies& currSpecies,
         const unsigned long long& np,
+        bool const isBTD = false);
+
+  /** This function sets constant particle records and ED-PIC attributes.
+   *
+   * Sets the entries for storing particle position offset, constant records (charge, mass) and ED-PIC attributes.
+   *
+   * @param[in] currSpecies Corresponding openPMD species
+   * @param[in] np          Number of particles
+   * @param[in] charge      Charge of the particles (note: fix for ions)
+   * @param[in] mass        Mass of the particles
+   */
+  void SetConstParticleRecordsEDPIC (
+        openPMD::ParticleSpecies& currSpecies,
+        const unsigned long long& np,
         amrex::ParticleReal const charge,
-        amrex::ParticleReal const mass, bool const isBTD = false);
+        amrex::ParticleReal const mass);
 
   /** This function sets up the entries for particle properties
    *
@@ -269,7 +283,9 @@ private:
    * @param[in] int_comp_names The int attribute names, from WarpX
    * @param[in] charge         Charge of the particles (note: fix for ions)
    * @param[in] mass           Mass of the particles
-   * @param[in] isBTD is this a backtransformed diagnostics write?
+   * @param[in] isBTD is this a backtransformed diagnostics (BTD) write?
+   * @param[in] isLastBTDFlush is this the last time we will flush this BTD station?
+   * @param[in] ParticleFlushOffset previously flushed number of particles in BTD
    */
   void DumpToFile (ParticleContainer* pc,
             const std::string& name,
@@ -279,7 +295,9 @@ private:
             const amrex::Vector<std::string>& real_comp_names,
             const amrex::Vector<std::string>&  int_comp_names,
             amrex::ParticleReal const charge,
-            amrex::ParticleReal const mass, const bool isBTD = false,
+            amrex::ParticleReal const mass,
+            const bool isBTD = false,
+            const bool isLastBTDFlush = false,
             int ParticleFlushOffset = 0);
 
   /** Get the openPMD-api filename for openPMD::Series
@@ -294,9 +312,6 @@ private:
 
   std::unique_ptr<openPMD::Series> m_Series;
 
-  /** To set particle meta-data for the openpmd dump. */
-  bool m_doParticleSetUp = true;
-
   /** This is the output directory
    *
    * This usually does not yet end in a `/`.
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index a08227779..f94325a72 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -541,7 +541,7 @@ WarpXOpenPMDPlot::Init (openPMD::Access access, bool isBTD)
 
 void
 WarpXOpenPMDPlot::WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& particle_diags,
-                  const bool isBTD, const amrex::Vector<int>& totalParticlesFlushedAlready)
+                  const bool isBTD, const bool isLastBTDFlush, const amrex::Vector<int>& totalParticlesFlushedAlready)
 {
   WARPX_PROFILE("WarpXOpenPMDPlot::WriteOpenPMDParticles()");
 
@@ -633,7 +633,8 @@ WarpXOpenPMDPlot::WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& part
              int_flags,
              real_names, int_names,
              pc->getCharge(), pc->getMass(),
-             isBTD, totalParticlesFlushedAlready[i]
+             isBTD, isLastBTDFlush,
+             totalParticlesFlushedAlready[i]
           );
       } else {
           DumpToFile(&tmp,
@@ -643,7 +644,8 @@ WarpXOpenPMDPlot::WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& part
              int_flags,
              real_names, int_names,
              pc->getCharge(), pc->getMass(),
-             isBTD, 0
+             isBTD, isLastBTDFlush,
+             0
           );
       }
     }
@@ -662,161 +664,156 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc,
                     const amrex::Vector<std::string>& real_comp_names,
                     const amrex::Vector<std::string>&  int_comp_names,
                     amrex::ParticleReal const charge,
-                    amrex::ParticleReal const mass, const bool isBTD,
-                    int ParticleFlushOffset)
-{
-  WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Series != nullptr, "openPMD: series must be initialized");
-
-  AMREX_ALWAYS_ASSERT(write_real_comp.size() == pc->NumRealComps());
-  AMREX_ALWAYS_ASSERT( write_int_comp.size() == pc->NumIntComps());
-  AMREX_ALWAYS_ASSERT(real_comp_names.size() == pc->NumRealComps());
-  AMREX_ALWAYS_ASSERT( int_comp_names.size() == pc->NumIntComps());
-
-  WarpXParticleCounter counter(pc);
-  if (counter.GetTotalNumParticles() == 0) return;
-  openPMD::Iteration currIteration = GetIteration(iteration, isBTD);
-
-  openPMD::ParticleSpecies currSpecies = currIteration.particles[name];
-  // meta data for ED-PIC extension
-  currSpecies.setAttribute( "particleShape", double( WarpX::noz ) );
-  // TODO allow this per direction in the openPMD standard, ED-PIC extension?
-  currSpecies.setAttribute( "particleShapes", [](){
-      return std::vector< double >{
-#if AMREX_SPACEDIM>=2
-          double(WarpX::nox),
-#endif
-#if defined(WARPX_DIM_3D)
-          double(WarpX::noy),
-#endif
-          double(WarpX::noz)
-      };
-  }() );
-  currSpecies.setAttribute( "particlePush", [](){
-      switch( WarpX::particle_pusher_algo ) {
-          case ParticlePusherAlgo::Boris : return "Boris";
-          case ParticlePusherAlgo::Vay : return "Vay";
-          case ParticlePusherAlgo::HigueraCary : return "HigueraCary";
-          default: return "other";
-      }
-  }() );
-  currSpecies.setAttribute( "particleInterpolation", [](){
-      switch( WarpX::field_gathering_algo ) {
-          case GatheringAlgo::EnergyConserving : return "energyConserving";
-          case GatheringAlgo::MomentumConserving : return "momentumConserving";
-          default: return "other";
-      }
-  }() );
-  currSpecies.setAttribute( "particleSmoothing", "none" );
-  currSpecies.setAttribute( "currentDeposition", [](){
-      switch( WarpX::current_deposition_algo ) {
-          case CurrentDepositionAlgo::Esirkepov : return "Esirkepov";
-          case CurrentDepositionAlgo::Vay : return "Vay";
-          default: return "directMorseNielson";
-      }
-  }() );
+                    amrex::ParticleReal const mass,
+                    const bool isBTD,
+                    const bool isLastBTDFlush,
+                    int ParticleFlushOffset) {
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Series != nullptr, "openPMD: series must be initialized");
+
+    AMREX_ALWAYS_ASSERT(write_real_comp.size() == pc->NumRealComps());
+    AMREX_ALWAYS_ASSERT(write_int_comp.size() == pc->NumIntComps());
+    AMREX_ALWAYS_ASSERT(real_comp_names.size() == pc->NumRealComps());
+    AMREX_ALWAYS_ASSERT(int_comp_names.size() == pc->NumIntComps());
+
+    WarpXParticleCounter counter(pc);
+    auto const num_dump_particles = counter.GetTotalNumParticles();
+
+    openPMD::Iteration currIteration = GetIteration(iteration, isBTD);
+    openPMD::ParticleSpecies currSpecies = currIteration.particles[name];
+
+    // prepare data structures the first time BTD has non-zero particles
+    //   we set some of them to zero extent, so we need to time that well
+    bool const is_first_flush_with_particles = num_dump_particles > 0 && ParticleFlushOffset == 0;
+    // BTD: we flush multiple times to the same lab step and thus need to resize
+    //   our declared particle output sizes
+    bool const is_resizing_flush = num_dump_particles > 0 && ParticleFlushOffset > 0;
+    // write structure & declare particles in this (lab) step empty:
+    //   if not BTD, then this is the only (and last) time we flush to this step
+    //   if BTD, then we may do this multiple times until it is the last BTD flush
+    bool const is_last_flush_to_step = !isBTD || (isBTD && isLastBTDFlush);
+    // well, even in BTD we have to recognize that some lab stations may have no
+    //   particles - so we mark them empty at the end of station reconstruction
+    bool const is_last_flush_and_never_particles =
+            is_last_flush_to_step && num_dump_particles == 0 && ParticleFlushOffset == 0;
 
-  //
-  // define positions & offsets
-  //
-  const unsigned long long NewParticleVectorSize = counter.GetTotalNumParticles() + ParticleFlushOffset;
-  m_doParticleSetUp = false;
-  if (counter.GetTotalNumParticles() > 0 and ParticleFlushOffset == 0) {
-      // This will trigger meta-data flush for particles the first-time non-zero number of particles are flushed.
-      m_doParticleSetUp = true;
-  }
-  SetupPos(currSpecies, NewParticleVectorSize, charge, mass, isBTD);
-  SetupRealProperties(pc, currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, NewParticleVectorSize, isBTD);
-  // open files from all processors, in case some will not contribute below
-  m_Series->flush();
-  for (auto currentLevel = 0; currentLevel <= pc->finestLevel(); currentLevel++)
-    {
-      uint64_t offset = static_cast<uint64_t>( counter.m_ParticleOffsetAtRank[currentLevel] );
-      // For BTD, the offset include the number of particles already flushed
-      if (isBTD) offset += ParticleFlushOffset;
-      for (ParticleIter pti(*pc, currentLevel); pti.isValid(); ++pti) {
-         auto const numParticleOnTile = pti.numParticles();
-         uint64_t const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile );
-
-         // Do not call storeChunk() with zero-sized particle tiles:
-         //   https://github.com/openPMD/openPMD-api/issues/1147
-         if (numParticleOnTile == 0) continue;
-
-         // get position and particle ID from aos
-         // note: this implementation iterates the AoS 4x...
-         // if we flush late as we do now, we can also copy out the data in one go
-         const auto& aos = pti.GetArrayOfStructs();  // size =  numParticlesOnTile
-         {
-           // Save positions
-           auto const positionComponents = detail::getParticlePositionComponentLabels();
+    //
+    // prepare structure and meta-data
+    //
+
+    // define positions & offset structure
+    const unsigned long long NewParticleVectorSize = num_dump_particles + ParticleFlushOffset;
+    // we will set up empty particles unless it's BTD, where we might add some in a following buffer dump
+    //   during this setup, we mark some particle properties as constant and potentially zero-sized
+    bool doParticleSetup = true;
+    if (isBTD)
+        doParticleSetup = is_first_flush_with_particles || is_last_flush_and_never_particles;
+
+    // this setup stage also implicitly calls "makeEmpty" if needed (i.e., is_last_flush_and_never_particles)
+    //   for BTD, we call this multiple times as we may resize in subsequent dumps if number of particles in the buffer > 0
+    if (doParticleSetup || is_resizing_flush) {
+        SetupPos(currSpecies, NewParticleVectorSize, isBTD);
+        SetupRealProperties(pc, currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names,
+                            NewParticleVectorSize, isBTD);
+    }
+
+    if (is_last_flush_to_step) {
+        SetConstParticleRecordsEDPIC(currSpecies, NewParticleVectorSize, charge, mass);
+    }
+
+    // open files from all processors, in case some will not contribute below
+    m_Series->flush();
+
+    // dump individual particles
+    for (auto currentLevel = 0; currentLevel <= pc->finestLevel(); currentLevel++) {
+        uint64_t offset = static_cast<uint64_t>( counter.m_ParticleOffsetAtRank[currentLevel] );
+        // For BTD, the offset include the number of particles already flushed
+        if (isBTD) offset += ParticleFlushOffset;
+        for (ParticleIter pti(*pc, currentLevel); pti.isValid(); ++pti) {
+            auto const numParticleOnTile = pti.numParticles();
+            uint64_t const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile );
+
+            // Do not call storeChunk() with zero-sized particle tiles:
+            //   https://github.com/openPMD/openPMD-api/issues/1147
+            if (numParticleOnTile == 0) continue;
+
+            // get position and particle ID from aos
+            // note: this implementation iterates the AoS 4x...
+            // if we flush late as we do now, we can also copy out the data in one go
+            const auto &aos = pti.GetArrayOfStructs();  // size =  numParticlesOnTile
+            {
+                // Save positions
+                auto const positionComponents = detail::getParticlePositionComponentLabels();
 #if defined(WARPX_DIM_RZ)
-           {
-              std::shared_ptr<amrex::ParticleReal> z(
-                      new amrex::ParticleReal[numParticleOnTile],
-                      [](amrex::ParticleReal const *p) { delete[] p; }
-              );
-              for (auto i = 0; i < numParticleOnTile; i++)
-                  z.get()[i] = aos[i].pos(1);  // {0: "r", 1: "z"}
-              std::string const positionComponent = "z";
-              currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64});
-           }
+                {
+                   std::shared_ptr<amrex::ParticleReal> z(
+                           new amrex::ParticleReal[numParticleOnTile],
+                           [](amrex::ParticleReal const *p) { delete[] p; }
+                   );
+                   for (auto i = 0; i < numParticleOnTile; i++)
+                       z.get()[i] = aos[i].pos(1);  // {0: "r", 1: "z"}
+                   std::string const positionComponent = "z";
+                   currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64});
+                }
 
-           //   reconstruct x and y from polar coordinates r, theta
-           auto const& soa = pti.GetStructOfArrays();
-           amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr();
-           WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer.");
-           WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile,
-                                            "openPMD: theta and tile size do not match");
-           {
-               std::shared_ptr< amrex::ParticleReal > x(
-                       new amrex::ParticleReal[numParticleOnTile],
-                       [](amrex::ParticleReal const *p){ delete[] p; }
-               );
-               std::shared_ptr< amrex::ParticleReal > y(
-                       new amrex::ParticleReal[numParticleOnTile],
-                       [](amrex::ParticleReal const *p){ delete[] p; }
-               );
-               for (auto i=0; i<numParticleOnTile; i++) {
-                   auto const r = aos[i].pos(0);  // {0: "r", 1: "z"}
-                   x.get()[i] = r * std::cos(theta[i]);
-                   y.get()[i] = r * std::sin(theta[i]);
-               }
-               currSpecies["position"]["x"].storeChunk(x, {offset}, {numParticleOnTile64});
-               currSpecies["position"]["y"].storeChunk(y, {offset}, {numParticleOnTile64});
-           }
+                //   reconstruct x and y from polar coordinates r, theta
+                auto const& soa = pti.GetStructOfArrays();
+                amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr();
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer.");
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile,
+                                                 "openPMD: theta and tile size do not match");
+                {
+                    std::shared_ptr< amrex::ParticleReal > x(
+                            new amrex::ParticleReal[numParticleOnTile],
+                            [](amrex::ParticleReal const *p){ delete[] p; }
+                    );
+                    std::shared_ptr< amrex::ParticleReal > y(
+                            new amrex::ParticleReal[numParticleOnTile],
+                            [](amrex::ParticleReal const *p){ delete[] p; }
+                    );
+                    for (auto i=0; i<numParticleOnTile; i++) {
+                        auto const r = aos[i].pos(0);  // {0: "r", 1: "z"}
+                        x.get()[i] = r * std::cos(theta[i]);
+                        y.get()[i] = r * std::sin(theta[i]);
+                    }
+                    currSpecies["position"]["x"].storeChunk(x, {offset}, {numParticleOnTile64});
+                    currSpecies["position"]["y"].storeChunk(y, {offset}, {numParticleOnTile64});
+                }
 #else
-           for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) {
-                std::shared_ptr< amrex::ParticleReal > curr(
-                    new amrex::ParticleReal[numParticleOnTile],
-                    [](amrex::ParticleReal const *p){ delete[] p; }
-                );
-                for (auto i=0; i<numParticleOnTile; i++) {
-                     curr.get()[i] = aos[i].pos(currDim);
+                for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) {
+                    std::shared_ptr<amrex::ParticleReal> curr(
+                            new amrex::ParticleReal[numParticleOnTile],
+                            [](amrex::ParticleReal const *p) { delete[] p; }
+                    );
+                    for (auto i = 0; i < numParticleOnTile; i++) {
+                        curr.get()[i] = aos[i].pos(currDim);
+                    }
+                    std::string const positionComponent = positionComponents[currDim];
+                    currSpecies["position"][positionComponent].storeChunk(curr, {offset},
+                                                                          {numParticleOnTile64});
                 }
-                std::string const positionComponent = positionComponents[currDim];
-                currSpecies["position"][positionComponent].storeChunk(curr, {offset}, {numParticleOnTile64});
-           }
 #endif
 
-           // save particle ID after converting it to a globally unique ID
-           std::shared_ptr< uint64_t > ids(
-               new uint64_t[numParticleOnTile],
-               [](uint64_t const *p){ delete[] p; }
-           );
-           for (auto i=0; i<numParticleOnTile; i++) {
-               ids.get()[i] = WarpXUtilIO::localIDtoGlobal( aos[i].id(), aos[i].cpu() );
-           }
-           auto const scalar = openPMD::RecordComponent::SCALAR;
-           currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64});
-         }
-         //  save "extra" particle properties in AoS and SoA
-         SaveRealProperty(pti,
-             currSpecies,
-             offset,
-             write_real_comp, real_comp_names,
-             write_int_comp, int_comp_names);
-
-         offset += numParticleOnTile64;
-      }
+                // save particle ID after converting it to a globally unique ID
+                std::shared_ptr<uint64_t> ids(
+                        new uint64_t[numParticleOnTile],
+                        [](uint64_t const *p) { delete[] p; }
+                );
+                for (auto i = 0; i < numParticleOnTile; i++) {
+                    ids.get()[i] = WarpXUtilIO::localIDtoGlobal(aos[i].id(), aos[i].cpu());
+                }
+                auto const scalar = openPMD::RecordComponent::SCALAR;
+                currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64});
+
+            }
+            //  save "extra" particle properties in AoS and SoA
+            SaveRealProperty(pti,
+                             currSpecies,
+                             offset,
+                             write_real_comp, real_comp_names,
+                             write_int_comp, int_comp_names);
+
+            offset += numParticleOnTile64;
+        }
     }
     m_Series->flush();
 }
@@ -856,9 +853,6 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc,
         }
     }
 
-    // attributes need to be set only the first time BTD flush is called for a snapshot
-    if (isBTD and m_doParticleSetUp == false) return;
-
     std::set< std::string > addedRecords; // add meta-data per record only once
     for (auto idx=0; idx<pc->NumRealComps(); idx++) {
         auto ii = ParticleContainer::NStructReal + idx; // jump over extra AoS names
@@ -976,51 +970,110 @@ void
 WarpXOpenPMDPlot::SetupPos (
     openPMD::ParticleSpecies& currSpecies,
     const unsigned long long& np,
-    amrex::ParticleReal const charge,
-    amrex::ParticleReal const mass,
     bool const isBTD)
 {
   std::string options = "{}";
   if (isBTD) options = "{ \"resizable\": true }";
-  auto realType= openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}, options);
+  auto realType = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}, options);
   auto idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}, options);
 
   auto const positionComponents = detail::getParticlePositionComponentLabels();
   for( auto const& comp : positionComponents ) {
-      currSpecies["positionOffset"][comp].resetDataset( realType );
       currSpecies["position"][comp].resetDataset( realType );
   }
 
   auto const scalar = openPMD::RecordComponent::SCALAR;
   currSpecies["id"][scalar].resetDataset( idType );
-  currSpecies["charge"][scalar].resetDataset( realType );
-  currSpecies["mass"][scalar].resetDataset( realType );
+}
 
-  if (isBTD and m_doParticleSetUp == false) return;
-  // make constant
-  for( auto const& comp : positionComponents ) {
-      currSpecies["positionOffset"][comp].makeConstant( 0. );
-  }
-  currSpecies["charge"][scalar].makeConstant( charge );
-  currSpecies["mass"][scalar].makeConstant( mass );
+void
+WarpXOpenPMDPlot::SetConstParticleRecordsEDPIC (
+        openPMD::ParticleSpecies& currSpecies,
+        const unsigned long long& np,
+        amrex::ParticleReal const charge,
+        amrex::ParticleReal const mass)
+{
+    auto realType = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np});
 
-  // meta data
-  currSpecies["position"].setUnitDimension( detail::getUnitDimension("position") );
-  currSpecies["positionOffset"].setUnitDimension( detail::getUnitDimension("positionOffset") );
-  currSpecies["charge"].setUnitDimension( detail::getUnitDimension("charge") );
-  currSpecies["mass"].setUnitDimension( detail::getUnitDimension("mass") );
-
-  // meta data for ED-PIC extension
-  currSpecies["position"].setAttribute( "macroWeighted", 0u );
-  currSpecies["position"].setAttribute( "weightingPower", 0.0 );
-  currSpecies["positionOffset"].setAttribute( "macroWeighted", 0u );
-  currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 );
-  currSpecies["id"].setAttribute( "macroWeighted", 0u );
-  currSpecies["id"].setAttribute( "weightingPower", 0.0 );
-  currSpecies["charge"].setAttribute( "macroWeighted", 0u );
-  currSpecies["charge"].setAttribute( "weightingPower", 1.0 );
-  currSpecies["mass"].setAttribute( "macroWeighted", 0u );
-  currSpecies["mass"].setAttribute( "weightingPower", 1.0 );
+    auto const positionComponents = detail::getParticlePositionComponentLabels();
+    for( auto const& comp : positionComponents ) {
+        currSpecies["positionOffset"][comp].resetDataset( realType );
+    }
+
+    // make constant
+    using namespace amrex::literals;
+    auto const scalar = openPMD::RecordComponent::SCALAR;
+    for( auto const& comp : positionComponents ) {
+        currSpecies["positionOffset"][comp].makeConstant( 0._prt );
+    }
+    currSpecies["charge"][scalar].makeConstant( charge );
+    currSpecies["mass"][scalar].makeConstant( mass );
+
+    // meta data
+    currSpecies["position"].setUnitDimension( detail::getUnitDimension("position") );
+    currSpecies["positionOffset"].setUnitDimension( detail::getUnitDimension("positionOffset") );
+    currSpecies["charge"].setUnitDimension( detail::getUnitDimension("charge") );
+    currSpecies["mass"].setUnitDimension( detail::getUnitDimension("mass") );
+
+    // meta data for ED-PIC extension
+    currSpecies["position"].setAttribute( "macroWeighted", 0u );
+    currSpecies["position"].setAttribute( "weightingPower", 0.0 );
+    currSpecies["positionOffset"].setAttribute( "macroWeighted", 0u );
+    currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 );
+    currSpecies["id"].setAttribute( "macroWeighted", 0u );
+    currSpecies["id"].setAttribute( "weightingPower", 0.0 );
+    currSpecies["charge"].setAttribute( "macroWeighted", 0u );
+    currSpecies["charge"].setAttribute( "weightingPower", 1.0 );
+    currSpecies["mass"].setAttribute( "macroWeighted", 0u );
+    currSpecies["mass"].setAttribute( "weightingPower", 1.0 );
+
+    // more ED-PIC attributes
+    currSpecies.setAttribute("particleShape", double(WarpX::noz));
+    // TODO allow this per direction in the openPMD standard, ED-PIC extension?
+    currSpecies.setAttribute("particleShapes", []() {
+        return std::vector<double>{
+#if AMREX_SPACEDIM >= 2
+                double(WarpX::nox),
+#endif
+#if defined(WARPX_DIM_3D)
+                double(WarpX::noy),
+#endif
+                double(WarpX::noz)
+        };
+    }());
+    currSpecies.setAttribute("particlePush", []() {
+        switch (WarpX::particle_pusher_algo) {
+            case ParticlePusherAlgo::Boris :
+                return "Boris";
+            case ParticlePusherAlgo::Vay :
+                return "Vay";
+            case ParticlePusherAlgo::HigueraCary :
+                return "HigueraCary";
+            default:
+                return "other";
+        }
+    }());
+    currSpecies.setAttribute("particleInterpolation", []() {
+        switch (WarpX::field_gathering_algo) {
+            case GatheringAlgo::EnergyConserving :
+                return "energyConserving";
+            case GatheringAlgo::MomentumConserving :
+                return "momentumConserving";
+            default:
+                return "other";
+        }
+    }());
+    currSpecies.setAttribute("particleSmoothing", "none");
+    currSpecies.setAttribute("currentDeposition", []() {
+        switch (WarpX::current_deposition_algo) {
+            case CurrentDepositionAlgo::Esirkepov :
+                return "Esirkepov";
+            case CurrentDepositionAlgo::Vay :
+                return "Vay";
+            default:
+                return "directMorseNielson";
+        }
+    }());
 }