Fix BTD/Scrape Flush Count with Filters (#4657)

* Fix BTD/Scrape Flush Count with Filters

Move the counting of already flushed particles for writers that call the
I/O backends multiple time per data set, e.g., BTD and boundary scraping,
into the I/O backend.

Currently, filtering is done as the first step in I/O backends and thus
the previous count outside of the I/O backends was over-counting
particles that might still get filtered out.

Offset should be a `long`:
Overflow risk is very high for pure `int`.
Also, counter is `unsigned`, so `unsigned long` for now.

* Simplify: Remove `m_totalParticles_flushed_already`

Less state we can forget in checkpoint-restart and that we
have to transfer across API boundaries.

Author: ax3l

Source/Diagnostics/BTDiagnostics.H

@@ -399,8 +399,6 @@ private:
         lab-frame data. */
     void InitializeParticleFunctors () override;
 
-    /** Update total number of particles flushed for all species for ith snapshot */
-    void UpdateTotalParticlesFlushed(int i_buffer);
     /** Reset total number of particles in the particle buffer to 0 for ith snapshot */
     void ResetTotalParticlesInBuffer(int i_buffer);
     /** Clear particle data stored in the particle buffer */

Source/Diagnostics/BTDiagnostics.cpp

@@ -129,7 +129,6 @@ void BTDiagnostics::DerivedInitData ()
         }
     }
     m_particles_buffer.resize(m_num_buffers);
-    m_totalParticles_flushed_already.resize(m_num_buffers);
     m_totalParticles_in_buffer.resize(m_num_buffers);
 
     // check that simulation can fill all BTD snapshots
@@ -1065,12 +1064,12 @@ BTDiagnostics::Flush (int i_buffer, bool force_flush)
         }
     }
     m_flush_format->WriteToFile(
-        m_varnames, m_mf_output[i_buffer], m_geom_output[i_buffer], warpx.getistep(),
-        labtime, m_output_species[i_buffer], nlev_output, file_name, m_file_min_digits,
+        m_varnames, m_mf_output.at(i_buffer), m_geom_output.at(i_buffer), warpx.getistep(),
+        labtime,
+        m_output_species.at(i_buffer), nlev_output, file_name, m_file_min_digits,
         m_plot_raw_fields, m_plot_raw_fields_guards,
-        use_pinned_pc, isBTD, i_buffer, m_buffer_flush_counter[i_buffer],
-        m_max_buffer_multifabs[i_buffer], m_geom_snapshot[i_buffer][0], isLastBTDFlush,
-        m_totalParticles_flushed_already[i_buffer]);
+        use_pinned_pc, isBTD, i_buffer, m_buffer_flush_counter.at(i_buffer),
+        m_max_buffer_multifabs.at(i_buffer), m_geom_snapshot.at(i_buffer).at(0), isLastBTDFlush);
 
     // Rescaling the box for plotfile after WriteToFile. This is because, for plotfiles, when writing particles, amrex checks if the particles are within the bounds defined by the box. However, in BTD, particles can be (at max) 1 cell outside the bounds of the geometry. So we keep a one-cell bigger box for plotfile when writing out the particle data and rescale after.
     if (m_format == "plotfile") {
@@ -1104,7 +1103,6 @@ BTDiagnostics::Flush (int i_buffer, bool force_flush)
     NullifyFirstFlush(i_buffer);
     // if particles are selected for output then update and reset counters
     if (!m_output_species_names.empty()) {
-        UpdateTotalParticlesFlushed(i_buffer);
         ResetTotalParticlesInBuffer(i_buffer);
         ClearParticleBuffer(i_buffer);
     }
@@ -1271,10 +1269,10 @@ void BTDiagnostics::MergeBuffersForPlotfile (int i_snapshot)
                 InterleaveSpeciesHeader(recent_species_Header,snapshot_species_Header,
                                         m_output_species_names[i], m_buffer_flush_counter[i_snapshot]);
                 if (BufferSpeciesHeader.m_total_particles == 0) { continue; }
-                if (m_totalParticles_flushed_already[i_snapshot][i]==0) {
-                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                    std::rename(recent_ParticleHdrFilename.c_str(), snapshot_ParticleHdrFilename.c_str()) == 0,
-                    std::string("Renaming ").append(recent_ParticleHdrFilename).append(" to ").append(snapshot_ParticleHdrFilename).append(" has failed"));
+                if (!amrex::FileExists(snapshot_ParticleHdrFilename)) {
+                    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                        std::rename(recent_ParticleHdrFilename.c_str(), snapshot_ParticleHdrFilename.c_str()) == 0,
+                        std::string("Renaming ").append(recent_ParticleHdrFilename).append(" to ").append(snapshot_ParticleHdrFilename).append(" has failed"));
                 } else {
                     InterleaveParticleDataHeader(recent_ParticleHdrFilename,
                                                  snapshot_ParticleHdrFilename);
@@ -1435,10 +1433,8 @@ BTDiagnostics::InitializeParticleBuffer ()
     const MultiParticleContainer& mpc = warpx.GetPartContainer();
     for (int i = 0; i < m_num_buffers; ++i) {
         m_particles_buffer[i].resize(m_output_species_names.size());
-        m_totalParticles_flushed_already[i].resize(m_output_species_names.size());
         m_totalParticles_in_buffer[i].resize(m_output_species_names.size());
         for (int isp = 0; isp < m_particles_buffer[i].size(); ++isp) {
-            m_totalParticles_flushed_already[i][isp] = 0;
             m_totalParticles_in_buffer[i][isp] = 0;
             m_particles_buffer[i][isp] = std::make_unique<PinnedMemoryParticleContainer>(WarpX::GetInstance().GetParGDB());
             const int idx = mpc.getSpeciesID(m_output_species_names[isp]);
@@ -1489,15 +1485,6 @@ BTDiagnostics::PrepareParticleDataForOutput()
     }
 }
 
-void
-BTDiagnostics::UpdateTotalParticlesFlushed(int i_buffer)
-{
-    for (int isp = 0; isp < m_totalParticles_flushed_already[i_buffer].size(); ++isp) {
-        m_totalParticles_flushed_already[i_buffer][isp] += static_cast<int>(
-            m_particles_buffer[i_buffer][isp]->TotalNumberOfParticles());
-    }
-}
-
 void
 BTDiagnostics::ResetTotalParticlesInBuffer(int i_buffer)
 {

Source/Diagnostics/BoundaryScrapingDiagnostics.cpp

@@ -102,15 +102,6 @@ BoundaryScrapingDiagnostics::InitializeParticleBuffer ()
             m_output_species[i_buffer].push_back(ParticleDiag(m_diag_name, species_name, pc, bnd_buffer));
         }
     }
-    // Initialize total number of particles flushed
-    m_totalParticles_flushed_already.resize(m_num_buffers);
-    for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) {
-        int const n_species = static_cast<int>(m_output_species_names.size());
-        m_totalParticles_flushed_already[i_buffer].resize(n_species);
-        for (int i_species=0; i_species<n_species; i_species++) {
-            m_totalParticles_flushed_already[i_buffer][i_species] = 0;
-        }
-    }
 }
 
 bool
@@ -157,11 +148,13 @@ BoundaryScrapingDiagnostics::Flush (int i_buffer, bool /* force_flush */)
     const std::string file_prefix = m_file_prefix + "/particles_at_" + particle_buffer.boundaryName(i_buffer);
 
     m_flush_format->WriteToFile(
-        m_varnames, m_mf_output[i_buffer], m_geom_output[i_buffer], warpx.getistep(),
-        warpx.gett_new(0), m_output_species[i_buffer], nlev_output, file_prefix,
+        m_varnames, m_mf_output.at(i_buffer), m_geom_output.at(i_buffer), warpx.getistep(),
+        warpx.gett_new(0),
+        m_output_species.at(i_buffer),
+        nlev_output, file_prefix,
         m_file_min_digits, false, false, use_pinned_pc, isBTD,
         warpx.getistep(0), bufferID, numBTDBuffers, geom,
-        isLastBTD, m_totalParticles_flushed_already[i_buffer]);
+        isLastBTD);
 
     // Now that the data has been written out, clear out the buffer
     particle_buffer.clearParticles(i_buffer);

Source/Diagnostics/Diagnostics.H

@@ -309,11 +309,6 @@ protected:
     /** Vector of pointers to functors to compute particle output per species*/
     amrex::Vector< std::unique_ptr<ComputeParticleDiagFunctor> > m_all_particle_functors;
 
-    /** Vector of total number of particles previously flushed, per species, per snapshot.
-     *  The first vector is for total number of snapshots and second vector loops
-     *  over the total number of species selected for diagnostics.
-     */
-    amrex::Vector< amrex::Vector <int> > m_totalParticles_flushed_already;
     /** Vector of total number of particles in the buffer, per species, per snapshot.
      *  The first vector is for total number of snapshots and second vector loops
      *  over the total number of species selected for diagnostics.

Source/Diagnostics/FlushFormats/FlushFormat.H

@@ -24,8 +24,7 @@ public:
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const = 0;
+        bool isLastBTDFlush = false) const = 0;
 
     FlushFormat () = default;
     virtual ~FlushFormat() = default;

Source/Diagnostics/FlushFormats/FlushFormatAscent.H

@@ -41,8 +41,7 @@ public:
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override;
+        bool isLastBTDFlush = false ) const override;
 
 #ifdef AMREX_USE_ASCENT
     /** \brief Do in-situ visualization for particle data.

Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp

@@ -21,7 +21,7 @@ FlushFormatAscent::WriteToFile (
     const bool /*use_pinned_pc*/,
     bool isBTD, int /*snapshotID*/, int /*bufferID*/, int /*numBuffers*/,
     const amrex::Geometry& /*full_BTD_snapshot*/,
-    bool /*isLastBTDFlush*/, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const
+    bool /*isLastBTDFlush*/) const
 {
 #ifdef AMREX_USE_ASCENT
     WARPX_PROFILE("FlushFormatAscent::WriteToFile()");

Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H

@@ -28,8 +28,7 @@ class FlushFormatCheckpoint final : public FlushFormatPlotfile
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const final;
+        bool isLastBTDFlush = false) const final;
 
     void CheckpointParticles (const std::string& dir,
                               const amrex::Vector<ParticleDiag>& particle_diags) const;

Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp

@@ -39,7 +39,7 @@ FlushFormatCheckpoint::WriteToFile (
         bool /*isBTD*/, int /*snapshotID*/,
         int /*bufferID*/, int /*numBuffers*/,
         const amrex::Geometry& /*full_BTD_snapshot*/,
-        bool /*isLastBTDFlush*/, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const
+        bool /*isLastBTDFlush*/) const
 {
     WARPX_PROFILE("FlushFormatCheckpoint::WriteToFile()");
 

Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.H

@@ -40,8 +40,7 @@ public:
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override;
+        bool isLastBTDFlush = false ) const override;
 
     ~FlushFormatOpenPMD () override = default;
 

Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp

@@ -126,7 +126,7 @@ FlushFormatOpenPMD::WriteToFile (
     const bool use_pinned_pc,
     bool isBTD, int snapshotID, int bufferID, int numBuffers,
     const amrex::Geometry& full_BTD_snapshot,
-    bool isLastBTDFlush, const amrex::Vector<int>& totalParticlesFlushedAlready) const
+    bool isLastBTDFlush) const
 {
     WARPX_PROFILE("FlushFormatOpenPMD::WriteToFile()");
     const std::string& filename = amrex::Concatenate(prefix, iteration[0], file_min_digits);
@@ -164,7 +164,7 @@ FlushFormatOpenPMD::WriteToFile (
 
     // particles: all (reside only on locally finest level)
     m_OpenPMDPlotWriter->WriteOpenPMDParticles(
-        particle_diags, static_cast<amrex::Real>(time), use_pinned_pc, isBTD, isLastBTDFlush, totalParticlesFlushedAlready);
+        particle_diags, static_cast<amrex::Real>(time), use_pinned_pc, isBTD, isLastBTDFlush);
 
     // signal that no further updates will be written to this iteration
     m_OpenPMDPlotWriter->CloseStep(isBTD, isLastBTDFlush);

Source/Diagnostics/FlushFormats/FlushFormatPlotfile.H

@@ -35,8 +35,7 @@ public:
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override;
+        bool isLastBTDFlush = false) const override;
 
     /** Write general info of the run into the plotfile */
     void WriteJobInfo(const std::string& dir) const;

Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp

@@ -65,7 +65,7 @@ FlushFormatPlotfile::WriteToFile (
     const bool /*use_pinned_pc*/,
     bool isBTD, int snapshotID,  int bufferID, int numBuffers,
     const amrex::Geometry& /*full_BTD_snapshot*/,
-    bool isLastBTDFlush, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const
+    bool isLastBTDFlush) const
 {
     WARPX_PROFILE("FlushFormatPlotfile::WriteToFile()");
     auto & warpx = WarpX::GetInstance();
@@ -340,9 +340,9 @@ FlushFormatPlotfile::WriteWarpXHeader(
 void
 FlushFormatPlotfile::WriteParticles(const std::string& dir,
                                     const amrex::Vector<ParticleDiag>& particle_diags,
-                                    const amrex::Real time, bool isBTD) const
+                                    const amrex::Real time,
+                                    bool isBTD) const
 {
-
     for (const auto& part_diag : particle_diags) {
         WarpXParticleContainer* pc = part_diag.getParticleContainer();
         PinnedMemoryParticleContainer* pinned_pc = part_diag.getPinnedParticleContainer();

Source/Diagnostics/FlushFormats/FlushFormatSensei.H

@@ -61,8 +61,7 @@ public:
         bool isBTD = false, int snapshotID = -1,
         int bufferID = 1, int numBuffers = 1,
         const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(),
-        bool isLastBTDFlush = false,
-        const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override;
+        bool isLastBTDFlush = false) const override;
 
     /** \brief Do in-situ visualization for particle data.
      * \param[in] particle_diags Each element of this vector handles output of 1 species.

Source/Diagnostics/FlushFormats/FlushFormatSensei.cpp

@@ -53,14 +53,13 @@ FlushFormatSensei::WriteToFile (
     bool plot_raw_fields, bool plot_raw_fields_guards,
     const bool use_pinned_pc,
     bool isBTD, int /*snapshotID*/, int /*bufferID*/, int /*numBuffers*/,
-    const amrex::Geometry& /*full_BTD_snapshot*/, bool /*isLastBTDFlush*/,
-    const amrex::Vector<int>& totalParticlesFlushedAlready) const
+    const amrex::Geometry& /*full_BTD_snapshot*/, bool /*isLastBTDFlush*/) const
 {
     amrex::ignore_unused(
         geom, nlev, prefix, file_min_digits,
         plot_raw_fields, plot_raw_fields_guards,
-        use_pinned_pc,
-        totalParticlesFlushedAlready);
+        use_pinned_pc
+    );
 
 #ifndef AMREX_USE_SENSEI_INSITU
     amrex::ignore_unused(varnames, mf, iteration, time, particle_diags,

Source/Diagnostics/FullDiagnostics.cpp

@@ -133,8 +133,9 @@ FullDiagnostics::Flush ( int i_buffer, bool /* force_flush */ )
     auto & warpx = WarpX::GetInstance();
 
     m_flush_format->WriteToFile(
-        m_varnames, m_mf_output[i_buffer], m_geom_output[i_buffer], warpx.getistep(),
-        warpx.gett_new(0), m_output_species[i_buffer], nlev_output, m_file_prefix,
+        m_varnames, m_mf_output.at(i_buffer), m_geom_output.at(i_buffer), warpx.getistep(),
+        warpx.gett_new(0),
+        m_output_species.at(i_buffer), nlev_output, m_file_prefix,
         m_file_min_digits, m_plot_raw_fields, m_plot_raw_fields_guards);
 
     FlushRaw();

Source/Diagnostics/OpenPMDHelpFunction.H

@@ -14,7 +14,25 @@
 #include <string>
 
 
+/** Determine the preferred file ending if unspecified
+ *
+ * @return file ending without the "."
+ */
 std::string
 WarpXOpenPMDFileType ();
 
+#ifdef WARPX_USE_OPENPMD
+/** Determine how many particles were already written in this species and step
+ *
+ * This checks for a particle species the current size of the id attribute, if it exists,
+ * and if it does it takes its extent as the number of particles already on disk.
+ *
+ * Note that this checks declared size, not necessarily written size.
+ *
+ * @return exisitng extent of the "id" attribute or zero.
+ */
+unsigned long
+num_already_flushed (openPMD::ParticleSpecies & currSpecies);
+#endif
+
 #endif  // WARPX_OPENPMDHELPFUNCTION_H_

Source/Diagnostics/OpenPMDHelpFunction.cpp

@@ -27,3 +27,23 @@ WarpXOpenPMDFileType ()
 #endif // WARPX_USE_OPENPMD
     return openPMDFileType;
 }
+
+#ifdef WARPX_USE_OPENPMD
+unsigned long
+num_already_flushed (openPMD::ParticleSpecies & currSpecies)
+{
+    const auto *const scalar = openPMD::RecordComponent::SCALAR;
+
+    unsigned long ParticleFlushOffset = 0;
+
+    if (currSpecies.contains("id")) {
+        if (currSpecies["id"].contains(scalar)) {
+            if (!currSpecies["id"][scalar].empty()) {
+                ParticleFlushOffset = currSpecies["id"][scalar].getExtent().at(0);
+            }
+        }
+    }
+
+    return ParticleFlushOffset;
+}
+#endif

Source/Diagnostics/WarpXOpenPMD.H

@@ -125,8 +125,7 @@ public:
               amrex::Real time,
               bool use_pinned_pc = false,
               bool isBTD = false,
-              bool isLastBTDFlush = false,
-              const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>());
+              bool isLastBTDFlush = false);
 
   /** Write out all openPMD fields for all active MR levels
    *
@@ -290,9 +289,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[inout] ParticleFlushOffset previously flushed number of particles in BTD
    * @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,
@@ -304,8 +303,7 @@ private:
             amrex::ParticleReal charge,
             amrex::ParticleReal mass,
             bool isBTD = false,
-            bool isLastBTDFlush = false,
-            int ParticleFlushOffset = 0);
+            bool isLastBTDFlush = false);
 
   /** Get the openPMD-api filename for openPMD::Series
    *