AMReX_ParticleContainer.H

Go to the documentation of this file.

#ifndef AMREX_PARTICLE_CONTAINER_H_
#define AMREX_PARTICLE_CONTAINER_H_
#include <AMReX_Config.H>
 
#include <AMReX_ParticleContainerBase.H>
#include <AMReX_ParmParse.H>
#include <AMReX_ParGDB.H>
#include <AMReX_REAL.H>
#include <AMReX_IntVect.H>
#include <AMReX_Array.H>
#include <AMReX_Vector.H>
#include <AMReX_Utility.H>
#include <AMReX_Geometry.H>
#include <AMReX_VisMF.H>
#include <AMReX_RealBox.H>
#include <AMReX_Print.H>
#include <AMReX_MultiFabUtil.H>
#include <AMReX_NFiles.H>
#include <AMReX_VectorIO.H>
#include <AMReX_Particle_mod_K.H>
#include <AMReX_ParticleMPIUtil.H>
#include <AMReX_StructOfArrays.H>
#include <AMReX_ArrayOfStructs.H>
#include <AMReX_Particle.H>
#include <AMReX_ParticleTile.H>
#include <AMReX_TypeTraits.H>
#include <AMReX_GpuContainers.H>
#include <AMReX_ParticleUtil.H>
#include <AMReX_ParticleReduce.H>
#include <AMReX_ParticleBufferMap.H>
#include <AMReX_ParticleCommunication.H>
#include <AMReX_ParticleLocator.H>
#include <AMReX_Scan.H>
#include <AMReX_DenseBins.H>
#include <AMReX_SparseBins.H>
#include <AMReX_ParticleTransformation.H>
#include <AMReX_ParticleMesh.H>
#include <AMReX_OpenMP.H>
#include <AMReX_ParIter.H>
 
#ifdef AMREX_LAZY
#include <AMReX_Lazy.H>
#endif
 
#ifdef AMREX_USE_OMP
#include <omp.h>
#endif
 
#ifdef AMREX_USE_HDF5
#include <hdf5.h>
#endif
 
#include <algorithm>
#include <array>
#include <cstring>
#include <fstream>
#include <iostream>
#include <limits>
#include <map>
#include <memory>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace amrex {
 
#ifdef AMREX_USE_HDF5_ASYNC
    extern hid_t es_par_g;
#endif
struct ParticleCommData
{
    using RealType = ParticleReal;
    int     m_lev;
    int     m_grid;
    IntVect m_cell;
    RealType m_data[1 + AMREX_SPACEDIM];
};
 
struct ParticleLocData
{
    int     m_lev  = -1;
    int     m_grid = -1;
    int     m_tile = -1;
    IntVect m_cell {AMREX_D_DECL(-1,-1,-1)};
    Box     m_gridbox;
    Box     m_tilebox;
    Box     m_grown_gridbox;
};
 
template<int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
struct ParticleInitType
{
    std::array<double, NStructReal> real_struct_data;
    std::array<int,    NStructInt > int_struct_data;
    std::array<double, NArrayReal > real_array_data;
    std::array<int,    NArrayInt  > int_array_data;
};
 
template <bool is_const, typename T_ParticleType, int NArrayReal, int NArrayInt,
          template<class> class Allocator, class CellAssignor>
class ParIterBase_impl;
 
template <typename T_ParticleType, int T_NArrayReal, int T_NArrayInt,
          template<class> class Allocator=DefaultAllocator,
          class T_CellAssignor=DefaultAssignor>
 
class ParticleContainer_impl : public ParticleContainerBase
{
public:
    using ParticleType = T_ParticleType;
    using ConstParticleType = typename ParticleType::ConstType;
    using CellAssignor = T_CellAssignor;
 
    static constexpr int NStructReal = ParticleType::NReal;
    static constexpr int NStructInt = ParticleType::NInt;
    static constexpr int NArrayReal = T_NArrayReal;
    static constexpr int NArrayInt = T_NArrayInt;
 
private:
    friend class ParIterBase_impl<true,ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
    friend class ParIterBase_impl<false,ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
 
public:
    template <typename T>
    using AllocatorType = Allocator<T>;
    using SuperParticleType = Particle<NStructReal+NArrayReal, NStructInt+NArrayInt>;
    using RealType = typename Particle<NStructReal, NStructInt>::RealType;
 
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
    RealDescriptor ParticleRealDescriptor = FPC::Native32RealDescriptor();
#else
    RealDescriptor ParticleRealDescriptor = FPC::Native64RealDescriptor();
#endif
 
    using ParticleContainerType = ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
    using ParticleTileType = ParticleTile<ParticleType, NArrayReal, NArrayInt, Allocator>;
    using ParticleInitData = ParticleInitType<NStructReal, NStructInt, NArrayReal, NArrayInt>;
 
    using ParticleLevel = std::map<std::pair<int, int>, ParticleTileType>;
    using PTDType = typename ParticleTileType::ParticleTileDataType;
    using ConstPTDType = typename ParticleTileType::ConstParticleTileDataType;
    using AoS = typename ParticleTileType::AoS;
    using SoA = typename ParticleTileType::SoA;
 
    using RealVector       = typename SoA::RealVector;
    using IntVector        = typename SoA::IntVector;
    using ParticleVector   = typename AoS::ParticleVector;
    using CharVector       = Gpu::DeviceVector<char>;
    using ParIterType      = ParIter_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
    using ParConstIterType = ParConstIter_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>;
 
    static constexpr bool has_polymorphic_allocator =
        IsPolymorphicArenaAllocator<Allocator<RealType>>::value;
 
    static_assert(!ParticleType::is_soa_particle || T_NArrayReal >= AMREX_SPACEDIM,
                  "Pure SoA particle containers require at least AMREX_SPACEDIM real components for positions");
 
    ParticleContainer_impl ()
      :
      ParticleContainerBase()
    {
        Initialize ();
    }
 
    ParticleContainer_impl (ParGDBBase* gdb)
        :
        ParticleContainerBase(gdb)
    {
        Initialize ();
        ParticleContainer_impl::reserveData();
        ParticleContainer_impl::resizeData();
    }
 
    ParticleContainer_impl (const Geometry            & geom,
                            const DistributionMapping & dmap,
                            const BoxArray            & ba)
        :
        ParticleContainerBase(geom, dmap, ba)
    {
        Initialize ();
        ParticleContainer_impl::reserveData();
        ParticleContainer_impl::resizeData();
    }
 
    ParticleContainer_impl (const Vector<Geometry>            & geom,
                            const Vector<DistributionMapping> & dmap,
                            const Vector<BoxArray>            & ba,
                            const Vector<int>                 & rr)
        :
        ParticleContainerBase(geom, dmap, ba, rr)
    {
        Initialize ();
        ParticleContainer_impl::reserveData();
        ParticleContainer_impl::resizeData();
    }
 
    ParticleContainer_impl (const Vector<Geometry>            & geom,
                            const Vector<DistributionMapping> & dmap,
                            const Vector<BoxArray>            & ba,
                            const Vector<IntVect>             & rr)
        :
        ParticleContainerBase(geom, dmap, ba, rr)
    {
        Initialize ();
        ParticleContainer_impl::reserveData();
        ParticleContainer_impl::resizeData();
    }
 
    ~ParticleContainer_impl () override = default;
 
    ParticleContainer_impl ( const ParticleContainer_impl &) = delete;
    ParticleContainer_impl& operator= ( const ParticleContainer_impl & ) = delete;
 
    ParticleContainer_impl ( ParticleContainer_impl && ) noexcept = default;
    ParticleContainer_impl& operator= ( ParticleContainer_impl && ) noexcept = default;
 
 
    void Define (ParGDBBase* gdb)
    {
        this->ParticleContainerBase::Define(gdb);
        reserveData();
        resizeData();
    }
 
    void Define (const Geometry            & geom,
                 const DistributionMapping & dmap,
                 const BoxArray            & ba)
    {
        this->ParticleContainerBase::Define(geom, dmap, ba);
        reserveData();
        resizeData();
    }
 
    void Define (const Vector<Geometry>            & geom,
                 const Vector<DistributionMapping> & dmap,
                 const Vector<BoxArray>            & ba,
                 const Vector<int>                 & rr)
    {
        this->ParticleContainerBase::Define(geom, dmap, ba, rr);
        reserveData();
        resizeData();
    }
 
    void Define (const Vector<Geometry>            & geom,
                 const Vector<DistributionMapping> & dmap,
                 const Vector<BoxArray>            & ba,
                 const Vector<IntVect>             & rr)
    {
        this->ParticleContainerBase::Define(geom, dmap, ba, rr);
        reserveData();
        resizeData();
    }
 
    int numLocalTilesAtLevel (int lev) const {
        return (lev < m_particles.size()) ? m_particles[lev].size() : 0;
    }
 
    void reserveData () override;
 
    void resizeData () override;
 
    void InitFromAsciiFile (const std::string& file, int extradata,
                            const IntVect* Nrep = nullptr);
 
    void InitFromBinaryFile (const std::string& file, int extradata);
 
    void InitFromBinaryMetaFile (const std::string& file, int extradata);
 
    void InitRandom (Long icount, ULong iseed,
                     const ParticleInitData& pdata,
                     bool serialize = false, RealBox bx = RealBox());
 
 
    void InitRandomPerBox (Long icount, ULong iseed, const ParticleInitData& pdata);
 
 
    void InitOnePerCell (Real x_off, Real y_off, Real z_off,
                         const ParticleInitData& pdata);
 
 
    void InitNRandomPerCell (int n_per_cell, const ParticleInitData& pdata);
 
    void Increment (MultiFab& mf, int level);
 
    Long IncrementWithTotal (MultiFab& mf, int level, bool local = false);
 
    void Redistribute (int lev_min = 0, int lev_max = -1, int nGrow = 0, int local=0,
                       bool remove_negative=true);
 
 
    template <class index_type>
    void ReorderParticles (int lev, const MFIter& mfi, const index_type* permutations);
 
    void SortParticlesForDeposition (IntVect idx_type);
 
    void SortParticlesByCell ();
 
    void SortParticlesByBin (IntVect bin_size);
 
    bool OK (int lev_min = 0, int lev_max = -1, int nGrow = 0) const;
 
    std::array<Long, 3> ByteSpread () const;
 
    std::array<Long, 3> PrintCapacity () const;
 
    void ShrinkToFit ();
 
    Long NumberOfParticlesAtLevel (int level, bool only_valid = true, bool only_local = false) const;
 
    Vector<Long> NumberOfParticlesInGrid  (int level, bool only_valid = true, bool only_local = false) const;
 
    template <typename I, std::enable_if_t<std::is_integral_v<I> && (sizeof(I) >= sizeof(Long)), int> = 0>
    void CapacityOfParticlesInGrid (LayoutData<I>& mem, int lev) const;
 
    Long TotalNumberOfParticles (bool only_valid=true, bool only_local=false) const;
 
 
    void RemoveParticlesAtLevel (int level);
 
    void RemoveParticlesNotAtFinestLevel ();
 
    void CreateVirtualParticles (int level, AoS& virts) const;
 
    void CreateGhostParticles (int level, int ngrow, AoS& ghosts) const;
 
    void AddParticlesAtLevel (AoS& particles, int level, int nGrow=0);
 
    void CreateVirtualParticles (int level, ParticleTileType& virts) const;
 
    void CreateGhostParticles (int level, int ngrow, ParticleTileType& ghosts) const;
 
    void AddParticlesAtLevel (ParticleTileType& particles, int level, int nGrow=0);
 
 
    void clearParticles ();
 
 
    template <class PCType,
              std::enable_if_t<IsParticleContainer<PCType>::value, int> foo = 0>
    void copyParticles (const PCType& other, bool local=false);
 
    template <class PCType,
              std::enable_if_t<IsParticleContainer<PCType>::value, int> foo = 0>
    void addParticles (const PCType& other, bool local=false);
 
    template <class F, class PCType,
              std::enable_if_t<IsParticleContainer<PCType>::value, int> foo = 0,
              std::enable_if_t<! std::is_integral_v<F>, int> bar = 0>
    void copyParticles (const PCType& other, F&&f, bool local=false);
 
    template <class F, class PCType,
              std::enable_if_t<IsParticleContainer<PCType>::value, int> foo = 0,
              std::enable_if_t<! std::is_integral_v<F>, int> bar = 0>
    void addParticles (const PCType& other, F const& f, bool local=false);
 
    void WriteParticleRealData (void* data, size_t size, std::ostream& os) const;
 
    void ReadParticleRealData (void* data, size_t size, std::istream& is);
 
    void Checkpoint (const std::string& dir, const std::string& name,
                     const Vector<std::string>& real_comp_names = Vector<std::string>(),
                     const Vector<std::string>& int_comp_names = Vector<std::string>()) const
    {
        Checkpoint(dir, name, true, real_comp_names, int_comp_names);
    }
 
    void Checkpoint (const std::string& dir, const std::string& name, bool is_checkpoint,
                     const Vector<std::string>& real_comp_names = Vector<std::string>(),
                     const Vector<std::string>& int_comp_names = Vector<std::string>()) const;
 
    void Checkpoint (const std::string& dir, const std::string& name,
                     const Vector<int>& write_real_comp,
                     const Vector<int>& write_int_comp,
                     const Vector<std::string>& real_comp_names,
                     const Vector<std::string>& int_comp_names) const;
 
    template <class F>
    void WriteBinaryParticleData (const std::string& dir,
                                  const std::string& name,
                                  const Vector<int>& write_real_comp,
                                  const Vector<int>& write_int_comp,
                                  const Vector<std::string>& real_comp_names,
                                  const Vector<std::string>&  int_comp_names,
                                  F&& f, bool is_checkpoint=false) const;
 
    void CheckpointPre ();
 
    void CheckpointPost ();
 
    void Restart (const std::string& dir, const std::string& file);
 
    void Restart (const std::string& dir, const std::string& file, bool is_checkpoint);
 
    void WritePlotFile (const std::string& dir, const std::string& name) const;
 
    template <class F, std::enable_if_t<!std::is_same_v<std::decay_t<F>, Vector<std::string>>>* = nullptr>
    void WritePlotFile (const std::string& dir, const std::string& name, F&& f) const;
 
    void WritePlotFile (const std::string& dir, const std::string& name,
                        const Vector<std::string>& real_comp_names,
                        const Vector<std::string>&  int_comp_names) const;
 
    template <class F>
    void WritePlotFile (const std::string& dir, const std::string& name,
                        const Vector<std::string>& real_comp_names,
                        const Vector<std::string>&  int_comp_names, F&& f) const;
 
    void WritePlotFile (const std::string& dir, const std::string& name,
                        const Vector<std::string>& real_comp_names) const;
 
    template <class F, std::enable_if_t<!std::is_same_v<std::decay_t<F>, Vector<std::string>>>* = nullptr>
    void WritePlotFile (const std::string& dir, const std::string& name,
                        const Vector<std::string>& real_comp_names, F&& f) const;
 
    void WritePlotFile (const std::string& dir,
                        const std::string& name,
                        const Vector<int>& write_real_comp,
                        const Vector<int>& write_int_comp) const;
 
    template <class F>
    void WritePlotFile (const std::string& dir,
                        const std::string& name,
                        const Vector<int>& write_real_comp,
                        const Vector<int>& write_int_comp, F&& f) const;
 
    void WritePlotFile (const std::string& dir,
                        const std::string& name,
                        const Vector<int>& write_real_comp,
                        const Vector<int>& write_int_comp,
                        const Vector<std::string>& real_comp_names,
                        const Vector<std::string>&  int_comp_names) const;
 
    template <class F>
    void WritePlotFile (const std::string& dir,
                        const std::string& name,
                        const Vector<int>& write_real_comp,
                        const Vector<int>& write_int_comp,
                        const Vector<std::string>& real_comp_names,
                        const Vector<std::string>&  int_comp_names,
                        F&& f) const;
 
    void WritePlotFilePre ();
 
    void WritePlotFilePost ();
 
    void WriteAsciiFile (const std::string& file);
 
    const Vector<ParticleLevel>& GetParticles () const { return m_particles; }
 
    Vector      <ParticleLevel>& GetParticles ()       { return m_particles; }
 
    const ParticleLevel& GetParticles (int lev) const { return m_particles[lev]; }
 
    ParticleLevel      & GetParticles (int lev)       { return m_particles[lev]; }
 
    const ParticleTileType& ParticlesAt (int lev, int grid, int tile) const
    { return m_particles[lev].at(std::make_pair(grid, tile)); }
 
    ParticleTileType&       ParticlesAt (int lev, int grid, int tile)
    { return m_particles[lev].at(std::make_pair(grid, tile)); }
 
    template <class Iterator>
    const ParticleTileType& ParticlesAt (int lev, const Iterator& iter) const
        { return ParticlesAt(lev, iter.index(), iter.LocalTileIndex()); }
 
    template <class Iterator>
    ParticleTileType&       ParticlesAt (int lev, const Iterator& iter)
        { return ParticlesAt(lev, iter.index(), iter.LocalTileIndex()); }
 
    ParticleTileType& DefineAndReturnParticleTile (int lev, int grid, int tile)
    {
        m_particles[lev][std::make_pair(grid, tile)].define(
            NumRuntimeRealComps(), NumRuntimeIntComps(),
            &m_soa_rdata_names, &m_soa_idata_names,
            arena()
        );
        return ParticlesAt(lev, grid, tile);
    }
 
    template <class Iterator>
    ParticleTileType& DefineAndReturnParticleTile (int lev, const Iterator& iter)
    {
        auto index = std::make_pair(iter.index(), iter.LocalTileIndex());
        m_particles[lev][index].define(
            NumRuntimeRealComps(), NumRuntimeIntComps(),
            &m_soa_rdata_names, &m_soa_idata_names,
            arena()
        );
        return ParticlesAt(lev, iter);
    }
 
    void AssignDensity (int rho_index,
                        Vector<std::unique_ptr<MultiFab> >& mf_to_be_filled,
                        int lev_min, int ncomp, int finest_level, int ngrow=2) const;
 
    void AssignCellDensitySingleLevel (int rho_index, MultiFab& mf, int level,
                                       int ncomp=1, int particle_lvl_offset = 0) const;
 
    template <typename P, typename Assignor=CellAssignor>
    IntVect Index (const P& p, int lev) const;
 
    ParticleLocData Reset (ParticleType& prt, bool update, bool verbose=true,
                           ParticleLocData pld = ParticleLocData()) const;
 
    template <typename P>
    bool PeriodicShift (P& p) const;
 
    void SetLevelDirectoriesCreated (bool tf) { levelDirectoriesCreated = tf; }
 
    bool GetLevelDirectoriesCreated () const { return levelDirectoriesCreated; }
 
    void SetUsePrePost (bool tf) const {
      usePrePost = tf;
    }
    bool GetUsePrePost () const {
      return usePrePost;
    }
 
    Long GetMaxNextIDPrePost () const { return maxnextidPrePost; }
    Long GetNParticlesPrePost () const { return nparticlesPrePost; }
 
    void SetUseUnlink (bool tf) const {
      doUnlink = tf;
    }
 
    bool GetUseUnlink () const {
      return doUnlink;
    }
 
    void RedistributeCPU (int lev_min = 0, int lev_max = -1, int nGrow = 0, int local=0,
                          bool remove_negative=true);
 
    void RedistributeGPU (int lev_min = 0, int lev_max = -1, int nGrow = 0, int local=0,
                          bool remove_negative=true);
 
    void ReserveForRedistribute (ParticleCopyPlan const& plan);
 
    Long superParticleSize() const { return superparticle_size; }
 
    void AddRealComp (std::string const & name, int communicate=1)
    {
        // names must be unique
        auto const it = std::find(m_soa_rdata_names.begin(), m_soa_rdata_names.end(), name);
        if (it != m_soa_rdata_names.end()) {
            throw std::runtime_error("AddRealComp: name '" + name + "' is already present in the SoA.");
        }
        m_soa_rdata_names.push_back(name);
 
        m_runtime_comps_defined = true;
        m_num_runtime_real++;
        h_redistribute_real_comp.push_back(communicate);
        SetParticleSize();
        this->resizeData();
 
        // resize runtime SoA
        for (int lev = 0; lev < numLevels(); ++lev) {
            for (ParIterType pti(*this,lev); pti.isValid(); ++pti) {
                auto& tile = DefineAndReturnParticleTile(lev, pti);
                auto np = tile.numParticles();
                if (np > 0) {
                    auto& soa = tile.GetStructOfArrays();
                    soa.resize(np);
                }
            }
        }
    }
 
    void AddRealComp (int communicate=1)
    {
        AddRealComp(getDefaultCompNameReal<ParticleType>(NArrayReal+m_num_runtime_real), communicate);
    }
 
    void AddIntComp (std::string const & name, int communicate=1)
    {
        // names must be unique
        auto const it = std::find(m_soa_idata_names.begin(), m_soa_idata_names.end(), name);
        if (it != m_soa_idata_names.end()) {
            throw std::runtime_error("AddIntComp: name '" + name + "' is already present in the SoA.");
        }
        m_soa_idata_names.push_back(name);
 
        m_runtime_comps_defined = true;
        m_num_runtime_int++;
        h_redistribute_int_comp.push_back(communicate);
        SetParticleSize();
        this->resizeData();
 
        // resize runtime SoA
        for (int lev = 0; lev < numLevels(); ++lev) {
            for (ParIterType pti(*this,lev); pti.isValid(); ++pti) {
                auto& tile = DefineAndReturnParticleTile(lev, pti);
                auto np = tile.numParticles();
                if (np > 0) {
                    auto& soa = tile.GetStructOfArrays();
                    soa.resize(np);
                }
            }
        }
    }
 
    void AddIntComp (int communicate=1)
    {
        AddIntComp(getDefaultCompNameInt<ParticleType>(NArrayInt+m_num_runtime_int), communicate);
    }
 
    int NumRuntimeRealComps () const { return m_num_runtime_real; }
    int NumRuntimeIntComps  () const { return m_num_runtime_int;  }
 
    int NumRealComps () const { return NArrayReal + NumRuntimeRealComps(); }
    int NumIntComps  () const { return NArrayInt  + NumRuntimeIntComps() ; }
 
    void ResizeRuntimeRealComp (int new_size, bool communicate);
 
    void ResizeRuntimeIntComp (int new_size, bool communicate);
 
    template <template<class> class NewAllocator=amrex::DefaultAllocator>
    using ContainerLike = amrex::ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, NewAllocator>;
 
    template <template<class> class NewAllocator=Allocator>
    ContainerLike<NewAllocator>
    make_alike () const
    {
        ContainerLike<NewAllocator> tmp(m_gdb);
 
        std::vector<std::string> const real_names = this->GetRealSoANames();
        std::vector<std::string> real_ct_names(NArrayReal);
        for (int ic = 0; ic < NArrayReal; ++ic) { real_ct_names.at(ic) = real_names[ic]; }
 
        std::vector<std::string> const int_names = this->GetIntSoANames();
        std::vector<std::string> int_ct_names(NArrayInt);
        for (int ic = 0; ic < NArrayInt; ++ic) { int_ct_names.at(ic) = int_names[ic]; }
 
        tmp.SetSoACompileTimeNames(real_ct_names, int_ct_names);
 
        // add runtime real comps to tmp
        for (int ic = 0; ic < this->NumRuntimeRealComps(); ++ic) {
            tmp.AddRealComp(real_names.at(ic + NArrayReal));
        }
 
        // add runtime int comps to tmp
        for (int ic = 0; ic < this->NumRuntimeIntComps(); ++ic) {
            tmp.AddIntComp(int_names.at(ic + NArrayInt));
        }
 
        tmp.h_redistribute_real_comp = h_redistribute_real_comp;
        tmp.h_redistribute_int_comp = h_redistribute_int_comp;
 
        return tmp;
    }
 
    Vector<int> h_redistribute_real_comp;
    Vector<int> h_redistribute_int_comp;
 
    mutable bool levelDirectoriesCreated;
    mutable bool usePrePost;
    mutable bool doUnlink;
    Long maxnextidPrePost;
    mutable int nOutFilesPrePost;
    Long nparticlesPrePost;
    Vector<Long> nParticlesAtLevelPrePost;
    mutable Vector<Vector<int>>  whichPrePost;
    mutable Vector<Vector<int>>  countPrePost;
    mutable Vector<Vector<Long>> wherePrePost;
    mutable std::string HdrFileNamePrePost;
    mutable Vector<std::string> filePrefixPrePost;
 
protected:
 
    template <typename P>
    bool Where (const P& prt, ParticleLocData& pld,
                int lev_min = 0, int lev_max = -1, int nGrow=0, int local_grid=-1) const;
 
 
    template <typename P>
    bool EnforcePeriodicWhere (P& prt, ParticleLocData& pld,
                               int lev_min = 0, int lev_max = -1, int local_grid=-1) const;
 
public:
    void
    WriteParticles (int level, std::ofstream& ofs, int fnum,
                    Vector<int>& which, Vector<int>& count, Vector<Long>& where,
                    const Vector<int>& write_real_comp, const Vector<int>& write_int_comp,
                    const Vector<std::map<std::pair<int, int>,IntVector>>& particle_io_flags, bool is_checkpoint) const;
 
#ifdef AMREX_USE_HDF5
#include "AMReX_ParticlesHDF5.H"
#endif
 
public: // NOLINT(readability-redundant-access-specifiers)
    void SetSoACompileTimeNames (std::vector<std::string> const & rdata_name, std::vector<std::string> const & idata_name);
 
    std::vector<std::string> GetRealSoANames () const {return m_soa_rdata_names;}
 
    std::vector<std::string> GetIntSoANames () const {return m_soa_idata_names;}
 
    bool HasRealComp (std::string const & name);
 
    bool HasIntComp (std::string const & name);
 
    int GetRealCompIndex (std::string const & name);
 
    int GetIntCompIndex (std::string const & name);
 
protected:
 
    template <class RTYPE>
    void ReadParticles (int cnt, int grd, int lev, std::ifstream& ifs, int finest_level_in_file, bool convert_ids);
 
    void SetParticleSize ();
 
    DenseBins<typename ParticleTileType::ParticleTileDataType> m_bins;
 
private:
    virtual void particlePostLocate (ParticleType& /*p*/, const ParticleLocData& /*pld*/,
                                     const int /*lev*/) {}
 
    virtual void correctCellVectors (int /*old_index*/, int /*new_index*/,
                                     int /*grid*/, const ParticleType& /*p*/) {}
 
    void RedistributeMPI (std::map<int, Vector<char> >& not_ours,
                          int lev_min = 0, int lev_max = 0, int nGrow = 0, int local=0);
 
    template <typename P>
    void locateParticle (P& p, ParticleLocData& pld,
                         int lev_min, int lev_max, int nGrow, int local_grid=-1) const;
 
    void Initialize ();
 
    bool m_runtime_comps_defined{false};
    int m_num_runtime_real{0};
    int m_num_runtime_int{0};
 
    size_t particle_size, superparticle_size;
    int num_real_comm_comps, num_int_comm_comps;
    Vector<ParticleLevel> m_particles;
 
    // names of both compile-time and runtime Real and Int SoA data
    std::vector<std::string> m_soa_rdata_names;
    std::vector<std::string> m_soa_idata_names;
};
 
template <int T_NStructReal, int T_NStructInt, int T_NArrayReal, int T_NArrayInt, template<class> class Allocator, class CellAssignor>
using ParticleContainer = ParticleContainer_impl<Particle<T_NStructReal, T_NStructInt>, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor>;
 
template <int T_NArrayReal, int T_NArrayInt, template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
using ParticleContainerPureSoA = ParticleContainer_impl<SoAParticle<T_NArrayReal, T_NArrayInt>, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor>;
 
}
 
#include "AMReX_ParticleInit.H"
#include "AMReX_ParticleContainerI.H"
#include "AMReX_ParticleIO.H"
 
#ifdef AMREX_USE_HDF5
#include "AMReX_ParticleHDF5.H"
#endif
 
#endif /*_PARTICLES_H_*/