Make buffer of scraped particles available to Python code (#2164)

* Added wrapper to get number of particle species tracked by the scraper

Not sure if this is going to be useful, but it demonstrates a method to get information
  from the ParticleBoundaryBuffer into Python.

* Stubbed out the main wrapper functions

* Added parameters to wrapper

* Added wrapper for getting the number of particles scraped of a species on a boundary

* added picmi arguments to scrape particles at the domain boundary

* Added wrapper to get the full particle buffer into python

* rearanged the getBuffer properties code a little

* Added docstrings +other suggested changes

* Added num_particles_impacted_boundary docstring

* fixed mistake in docstring

* Changed boundary parameter to be a string for clarity

* Fixed issue with the boundary parameter for scraping

* Fixed issue with the boundary input for scraping stats wrapper

* Added demonstration of particle scraping wrapper

* Added analysis.py file

* Fix typo in one of the dimension maps

Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>

* Added before esolve to warpx evolve

* added test for the scraped particle buffer wrappers

* Moved python PICMI particle boundary scrape test

* Renamed test file to the correct name

* Removed old test

* added special functionality to get the timestep at which particles were scraped

* removed debug print

* added python wrapper for the clearParticles() function of the scraper buffer

* added special wrapper function to get the timesteps at which the particles in the boundary buffer were scraped

* updated test to match the non-PICMI test for the particle scraper buffer

* Fix uncaught rebase mistake

* re-activated picmi test of accessing the scraped particle buffers via python

* added documentation for the new parameters involved in the scraped particle buffer and fixed remaining issue with picmi test

* changes requested during code review

Co-authored-by: mkieburtz <michaelkieburtz@gmail.com>
Co-authored-by: Roelof <roelof.groenewald@modernelectron.com>
Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>

11 files changed, 435 insertions, 3 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index b329928f6..8e5bc77de 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -759,6 +759,27 @@ Particle initialization
     If `1` is given, this species will not be pushed
     by any pusher during the simulation.
 
+* ``<species>.save_particles_at_xlo/ylo/zlo``,  ``<species>.save_particles_at_xhi/yhi/zhi`` and ``<species>.save_particles_at_eb`` (`0` or `1` optional, default `0`)
+    If `1` particles of this species will be copied to the scraped particle
+    buffer for the specified boundary if they leave the simulation domain in
+    the specified direction. **If USE_EB=TRUE** the ``save_particles_at_eb``
+    flag can be set to `1` to also save particle data for the particles of this
+    species that impact the embedded boundary.
+    The scraped particle buffer can be used to track particle fluxes out of the
+    simulation but is currently only accessible via the Python interface. The
+    ``pywarpx._libwarpx`` function ``get_particle_boundary_buffer()`` can be
+    used to access the scraped particle buffer. An entry is included for every
+    particle in the buffer of the timestep at which the particle was scraped.
+    This can be accessed by passing the argument ``comp_name="step_scraped"`` to
+    the above mentioned function.
+
+    .. note::
+        Currently the scraped particle buffer relies on the user to access the
+        data in the buffer for processing and periodically clear the buffer. The
+        buffer will grow unbounded as particles are scraped and therefore could
+        lead to memory issues if not periodically cleared. To clear the buffer
+        call ``warpx_clearParticleBoundaryBuffer()``.
+
 * ``<species>.do_back_transformed_diagnostics`` (`0` or `1` optional, default `1`)
     Only used when ``warpx.do_back_transformed_diagnostics=1``. When running in a
     boosted frame, whether or not to plot back-transformed diagnostics for
diff --git a/Examples/Modules/ParticleBoundaryScrape/PICMI_inputs_scrape.py b/Examples/Modules/ParticleBoundaryScrape/PICMI_inputs_scrape.py
new file mode 100644
index 000000000..79cf01a57
--- /dev/null
+++ b/Examples/Modules/ParticleBoundaryScrape/PICMI_inputs_scrape.py
@@ -0,0 +1,131 @@
+# --- Input file to test the particle scraper and the Python wrappers
+# --- to access the buffer of scraped particles.
+
+from pywarpx import picmi
+
+##########################
+# numerics parameters
+##########################
+
+# --- Number of time steps
+max_steps = 60
+diagnostic_intervals = 20
+
+# --- Grid
+nx = 64
+ny = 64
+nz = 128
+
+cfl = 0.99
+
+xmin = -125e-6
+ymin = -125e-6
+zmin = -149e-6
+xmax = 125e-6
+ymax = 125e-6
+zmax = 1e-6
+
+##########################
+# physics components
+##########################
+
+uniform_plasma_elec = picmi.UniformDistribution(
+    density = 1e23, # number of electrons per m^3
+    lower_bound = [-1e-5, -1e-5, -149e-6],
+    upper_bound = [1e-5, 1e-5, -129e-6],
+    directed_velocity = [0., 0., 2000.*picmi.constants.c] # uth the std of the (unitless) momentum
+)
+
+electrons = picmi.Species(
+    particle_type='electron', name='electrons',
+    initial_distribution=uniform_plasma_elec,
+    warpx_save_particles_at_xhi=1, warpx_save_particles_at_eb=1
+)
+
+##########################
+# numerics components
+##########################
+
+grid = picmi.Cartesian3DGrid(
+    number_of_cells = [nx, ny, nz],
+    lower_bound = [xmin, ymin, zmin],
+    upper_bound = [xmax, ymax, zmax],
+    lower_boundary_conditions=['none', 'none', 'none'],
+    upper_boundary_conditions=['none', 'none', 'none'],
+    lower_boundary_conditions_particles=['open', 'open', 'open'],
+    upper_boundary_conditions_particles=['open', 'open', 'open'],
+    warpx_max_grid_size = 128
+)
+
+solver = picmi.ElectromagneticSolver(
+    grid=grid, cfl=cfl
+)
+
+embedded_boundary = picmi.EmbeddedBoundary(
+    implicit_function="-max(max(max(x-12.5e-6,-12.5e-6-x),max(y-12.5e-6,-12.5e-6-y)),max(z-(-6.15e-5),-8.65e-5-z))"
+)
+
+##########################
+# diagnostics
+##########################
+
+field_diag = picmi.FieldDiagnostic(
+    name = 'diag1',
+    grid = grid,
+    period = diagnostic_intervals,
+    data_list = ['Ex', 'Ey', 'Ez', 'Bx', 'By', 'Bz'],
+    write_dir = '.',
+    warpx_file_prefix = 'Python_particle_scrape_plt'
+)
+
+##########################
+# simulation setup
+##########################
+
+sim = picmi.Simulation(
+    solver = solver,
+    max_steps = max_steps,
+    warpx_embedded_boundary=embedded_boundary,
+    verbose=True
+)
+
+sim.add_species(
+    electrons,
+    layout = picmi.GriddedLayout(
+        n_macroparticle_per_cell=[1, 1, 1], grid=grid
+    )
+)
+
+sim.add_diagnostic(field_diag)
+
+##########################
+# simulation run
+##########################
+
+# sim.write_input_file(file_name = 'inputs_from_PICMI')
+sim.step(max_steps)
+
+################################################
+# check that the wrappers to access the particle
+# buffer functions as intended
+################################################
+
+from pywarpx import _libwarpx
+
+n = _libwarpx.get_particle_boundary_buffer_size("electrons", 'eb')
+print("Number of electrons in buffer:", n)
+assert n == 612
+
+scraped_steps = _libwarpx.get_particle_boundary_buffer("electrons", 'eb', 'step_scraped', 0)
+for arr in scraped_steps:
+    assert all(arr > 40)
+
+weights = _libwarpx.get_particle_boundary_buffer("electrons", 'eb', 'w', 0)
+assert sum(len(arr) for arr in weights) == 612
+
+# clear the particle buffer
+_libwarpx.libwarpx.warpx_clearParticleBoundaryBuffer()
+# confirm that the buffer was cleared
+n = _libwarpx.get_particle_boundary_buffer_size("electrons", 'eb')
+print("Number of electrons in buffer:", n)
+assert n == 0
diff --git a/Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py b/Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py
index b970c4933..c325495db 100755
--- a/Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py
+++ b/Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py
@@ -1,6 +1,7 @@
 #! /usr/bin/env python
 
 import yt
+from pathlib import Path
 
 # This test shoots a beam of electrons at cubic embedded boundary geometry
 # At time step 40, none of the particles have hit the boundary yet. At time
@@ -9,11 +10,19 @@ import yt
 # the problem domain yet.
 
 # all particles are still there
-ds40 = yt.load("particle_scrape_plt00040")
+if Path("particle_scrape_plt00040").is_dir():
+    filename = "particle_scrape_plt00040"
+else:
+    filename = "Python_particle_scrape_plt00040"
+ds40 = yt.load(filename)
 np40 = ds40.index.particle_headers['electrons'].num_particles
 assert(np40 == 612)
 
 # all particles have been removed
-ds60 = yt.load("particle_scrape_plt00060")
+if Path("particle_scrape_plt00060").is_dir():
+    filename = "particle_scrape_plt00060"
+else:
+    filename = "Python_particle_scrape_plt00060"
+ds60 = yt.load(filename)
 np60 = ds60.index.particle_headers['electrons'].num_particles
 assert(np60 == 0)
diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py
index 110de5190..ca87f4f41 100755
--- a/Python/pywarpx/_libwarpx.py
+++ b/Python/pywarpx/_libwarpx.py
@@ -120,6 +120,7 @@ class Particle(ctypes.Structure):
 # some useful typenames
 _LP_particle_p = ctypes.POINTER(ctypes.POINTER(Particle))
 _LP_c_int = ctypes.POINTER(ctypes.c_int)
+_LP_LP_c_int = ctypes.POINTER(_LP_c_int)
 _LP_c_void_p = ctypes.POINTER(ctypes.c_void_p)
 _LP_c_real = ctypes.POINTER(c_real)
 _LP_LP_c_real = ctypes.POINTER(_LP_c_real)
@@ -185,6 +186,9 @@ libwarpx.warpx_getChargeDensityCP.restype = _LP_LP_c_real
 libwarpx.warpx_getChargeDensityCPLoVects.restype = _LP_c_int
 libwarpx.warpx_getChargeDensityFP.restype = _LP_LP_c_real
 libwarpx.warpx_getChargeDensityFPLoVects.restype = _LP_c_int
+libwarpx.warpx_getParticleBoundaryBufferSize.restype = ctypes.c_int
+libwarpx.warpx_getParticleBoundaryBuffer.restype = _LP_LP_c_particlereal
+libwarpx.warpx_getParticleBoundaryBufferScrapedSteps.restype = _LP_LP_c_int
 
 libwarpx.warpx_getEx_nodal_flag.restype = _LP_c_int
 libwarpx.warpx_getEy_nodal_flag.restype = _LP_c_int
@@ -763,6 +767,127 @@ def add_real_comp(species_name, pid_name, comm=True):
     )
 
 
+def _get_boundary_number(boundary):
+    '''
+
+    Utility function to find the boundary number given a boundary name.
+
+    Parameters
+    ----------
+
+        boundary       : the boundary from which to get the scraped particle data.
+                         In the form x/y/z_hi/lo or eb.
+
+    Returns
+    -------
+
+    Integer index in the boundary scraper buffer for the given boundary.
+    '''
+    if geometry_dim == '3d':
+        dimensions = {'x' : 0, 'y' : 1, 'z' : 2}
+    elif geometry_dim == '2d':
+        dimensions = {'x' : 0, 'z' : 1}
+    else:
+        raise NotImplementedError("RZ is not supported for particle scraping.")
+
+    if boundary != 'eb':
+        boundary_parts = boundary.split("_")
+        dim_num = dimensions[boundary_parts[0]]
+        if boundary_parts[1] == 'lo':
+            side = 0
+        elif boundary_parts[1] == 'hi':
+            side = 1
+        else:
+            raise RuntimeError(f'Unknown boundary specified: {boundary}')
+        boundary_num = 2 * dim_num + side
+    else:
+        boundary_num = 4 if geometry_dim == '2d' else 6
+
+    return boundary_num
+
+
+def get_particle_boundary_buffer_size(species_name, boundary):
+    '''
+
+    This returns the number of particles that have been scraped so far in the simulation
+    from the specified boundary and of the specified species.
+
+    Parameters
+    ----------
+
+        species_name   : return the number of scraped particles of this species
+        boundary       : the boundary from which to get the scraped particle data.
+                         In the form x/y/z_hi/lo
+
+    Returns
+    -------
+
+        The number of particles scraped so far from a boundary and of a species.
+
+    '''
+    return libwarpx.warpx_getParticleBoundaryBufferSize(
+        ctypes.c_char_p(species_name.encode('utf-8')),
+        _get_boundary_number(boundary)
+    )
+
+
+def get_particle_boundary_buffer(species_name, boundary, comp_name, level):
+    '''
+
+    This returns a list of numpy arrays containing the particle array data
+    for a species that has been scraped by a specific simulation boundary.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        species_name   : the species name that the data will be returned for.
+        boundary       : the boundary from which to get the scraped particle data.
+                         In the form x/y/z_hi/lo or eb.
+        comp_name      : the component of the array data that will be returned.
+                         If "step_scraped" the special attribute holding the
+                         timestep at which a particle was scraped will be
+                         returned.
+        level          : Which AMR level to retrieve scraped particle data from.
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+    particles_per_tile = _LP_c_int()
+    num_tiles = ctypes.c_int(0)
+    if comp_name == 'step_scraped':
+        data = libwarpx.warpx_getParticleBoundaryBufferScrapedSteps(
+            ctypes.c_char_p(species_name.encode('utf-8')),
+            _get_boundary_number(boundary), level,
+            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
+        )
+    else:
+        data = libwarpx.warpx_getParticleBoundaryBuffer(
+            ctypes.c_char_p(species_name.encode('utf-8')),
+            _get_boundary_number(boundary), level,
+            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile),
+            ctypes.c_char_p(comp_name.encode('utf-8'))
+        )
+
+    particle_data = []
+    for i in range(num_tiles.value):
+        arr = np.ctypeslib.as_array(data[i], (particles_per_tile[i],))
+        try:
+            # This fails on some versions of numpy
+            arr.setflags(write=1)
+        except ValueError:
+            pass
+        particle_data.append(arr)
+
+    _libc.free(particles_per_tile)
+    _libc.free(data)
+    return particle_data
+
+
 def _get_mesh_field_list(warpx_func, level, direction, include_ghosts):
     """
      Generic routine to fetch the list of field data arrays.
@@ -1677,6 +1802,7 @@ def get_mesh_charge_density_cp_lovects(level, include_ghosts=True):
     '''
     return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.warpx_getChargeDensityCPLoVects)
 
+
 def get_mesh_charge_density_fp_lovects(level, include_ghosts=True):
     '''
 
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index e68636b89..71c18769b 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -85,6 +85,15 @@ class Species(picmistandard.PICMI_Species):
         self.self_fields_verbosity = kw.pop('warpx_self_fields_verbosity', None)
         self.save_previous_position = kw.pop('warpx_save_previous_position', None)
 
+        # For the scraper buffer
+        self.save_particles_at_xlo = kw.pop('warpx_save_particles_at_xlo', None)
+        self.save_particles_at_xhi = kw.pop('warpx_save_particles_at_xhi', None)
+        self.save_particles_at_ylo = kw.pop('warpx_save_particles_at_ylo', None)
+        self.save_particles_at_yhi = kw.pop('warpx_save_particles_at_yhi', None)
+        self.save_particles_at_zlo = kw.pop('warpx_save_particles_at_zlo', None)
+        self.save_particles_at_zhi = kw.pop('warpx_save_particles_at_zhi', None)
+        self.save_particles_at_eb = kw.pop('warpx_save_particles_at_eb', None)
+
     def initialize_inputs(self, layout,
                           initialize_self_fields = False,
                           injection_plane_position = None,
@@ -108,6 +117,13 @@ class Species(picmistandard.PICMI_Species):
                                              self_fields_required_precision = self.self_fields_required_precision,
                                              self_fields_max_iters = self.self_fields_max_iters,
                                              self_fields_verbosity = self.self_fields_verbosity,
+                                             save_particles_at_xlo = self.save_particles_at_xlo,
+                                             save_particles_at_xhi = self.save_particles_at_xhi,
+                                             save_particles_at_ylo = self.save_particles_at_ylo,
+                                             save_particles_at_yhi = self.save_particles_at_yhi,
+                                             save_particles_at_zlo = self.save_particles_at_zlo,
+                                             save_particles_at_zhi = self.save_particles_at_zhi,
+                                             save_particles_at_eb = self.save_particles_at_eb,
                                              save_previous_position = self.save_previous_position)
         pywarpx.Particles.particles_list.append(self.species)
 
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index d6e4611af..bcab916f0 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -2348,6 +2348,25 @@ particleTypes = electrons
 analysisRoutine = Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py
 tolerance = 1.0e-4
 
+[Python_particle_scrape]
+buildDir = .
+inputFile = Examples/Modules/ParticleBoundaryScrape/PICMI_inputs_scrape.py
+runtime_params =
+customRunCmd = python PICMI_inputs_scrape.py
+dim = 3
+addToCompileString = USE_EB=TRUE USE_PYTHON_MAIN=TRUE PYINSTALLOPTIONS="--user --prefix="
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 0
+numthreads = 0
+compileTest = 0
+doVis = 0
+compareParticles = 1
+particleTypes = electrons
+analysisRoutine = Examples/Modules/ParticleBoundaryScrape/analysis_scrape.py
+tolerance = 1.0e-4
+
 [Python_particle_attr_access]
 buildDir = .
 inputFile = Examples/Tests/ParticleDataPython/PICMI_inputs_2d.py
diff --git a/Source/Particles/ParticleBoundaryBuffer.H b/Source/Particles/ParticleBoundaryBuffer.H
index 452722d56..592dc574f 100644
--- a/Source/Particles/ParticleBoundaryBuffer.H
+++ b/Source/Particles/ParticleBoundaryBuffer.H
@@ -40,6 +40,10 @@ public:
 
     void printNumParticles () const;
 
+    int getNumParticlesInContainer(const std::string species_name, int boundary);
+
+    ParticleBuffer::BufferType<amrex::PinnedArenaAllocator>& getParticleBuffer(const std::string species_name, int boundary);
+
     constexpr int numBoundaries () const {
         return AMREX_SPACEDIM*2
 #ifdef AMREX_USE_EB
diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp
index f84d05b28..16465ac76 100644
--- a/Source/Particles/ParticleBoundaryBuffer.cpp
+++ b/Source/Particles/ParticleBoundaryBuffer.cpp
@@ -97,7 +97,7 @@ void ParticleBoundaryBuffer::printNumParticles () const {
                 int np = buffer[i].isDefined() ? buffer[i].TotalNumberOfParticles(false) : 0;
                 amrex::Print() << "Species " << getSpeciesNames()[i] << " has "
                                << np << " particles in the boundary buffer "
-                               << " for side " << iside << " of dim " << idim << "\n";
+                               << "for side " << iside << " of dim " << idim << "\n";
             }
         }
     }
@@ -232,3 +232,28 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
     amrex::ignore_unused(distance_to_eb, dxi);
 #endif
 }
+
+int ParticleBoundaryBuffer::getNumParticlesInContainer(
+        const std::string species_name, int boundary) {
+
+    auto& buffer = m_particle_containers[boundary];
+    auto index = WarpX::GetInstance().GetPartContainer().getSpeciesID(species_name);
+
+    if (buffer[index].isDefined()) return buffer[index].TotalNumberOfParticles(false);
+    else return 0;
+}
+
+ParticleBuffer::BufferType<amrex::PinnedArenaAllocator>&
+ParticleBoundaryBuffer::getParticleBuffer(const std::string species_name, int boundary) {
+
+    auto& buffer = m_particle_containers[boundary];
+    auto index = WarpX::GetInstance().GetPartContainer().getSpeciesID(species_name);
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_boundary_buffer[boundary][index],
+                                     "Attempted to get particle buffer for boundary "
+                                     + boundary + ", which is not used!");
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(buffer[index].isDefined(),
+                                     "Tried to get a buffer that is not defined!");
+
+    return buffer[index];
+}
diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp
index 1450011f9..dd4cf0e75 100644
--- a/Source/Python/WarpXWrappers.cpp
+++ b/Source/Python/WarpXWrappers.cpp
@@ -9,6 +9,7 @@
 #include "BoundaryConditions/PML.H"
 #include "Initialization/WarpXAMReXInit.H"
 #include "Particles/MultiParticleContainer.H"
+#include "Particles/ParticleBoundaryBuffer.H"
 #include "Particles/WarpXParticleContainer.H"
 #include "Utils/WarpXUtil.H"
 #include "WarpX.H"
@@ -471,6 +472,72 @@ extern "C"
         mypc.defineAllParticleTiles();
     }
 
+    int warpx_getParticleBoundaryBufferSize(const char* species_name, int boundary)
+    {
+        const std::string name(species_name);
+        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
+        return particle_buffers.getNumParticlesInContainer(species_name, boundary);
+    }
+
+    int** warpx_getParticleBoundaryBufferScrapedSteps(const char* species_name, int boundary, int lev,
+                     int* num_tiles, int** particles_per_tile)
+    {
+        const std::string name(species_name);
+        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
+        auto& particle_buffer = particle_buffers.getParticleBuffer(species_name, boundary);
+
+        const int comp = particle_buffer.NumIntComps() - 1;
+
+        int i = 0;
+        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {}
+
+        // *num_tiles = myspc.numLocalTilesAtLevel(lev);
+        *num_tiles = i;
+        *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int));
+
+        int** data = (int**) malloc(*num_tiles*sizeof(int*));
+        i = 0;
+        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {
+            auto& soa = pti.GetStructOfArrays();
+            data[i] = (int*) soa.GetIntData(comp).dataPtr();
+            (*particles_per_tile)[i] = pti.numParticles();
+        }
+
+        return data;
+    }
+
+    amrex::ParticleReal** warpx_getParticleBoundaryBuffer(const char* species_name, int boundary, int lev,
+                     int* num_tiles, int** particles_per_tile, const char* comp_name)
+    {
+        const std::string name(species_name);
+        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
+        auto& particle_buffer = particle_buffers.getParticleBuffer(species_name, boundary);
+
+        const int comp = warpx_getParticleCompIndex(species_name, comp_name);
+
+        int i = 0;
+        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {}
+
+        // *num_tiles = myspc.numLocalTilesAtLevel(lev);
+        *num_tiles = i;
+        *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int));
+
+        amrex::ParticleReal** data = (amrex::ParticleReal**) malloc(*num_tiles*sizeof(amrex::ParticleReal*));
+        i = 0;
+        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {
+            auto& soa = pti.GetStructOfArrays();
+            data[i] = (amrex::ParticleReal*) soa.GetRealData(comp).dataPtr();
+            (*particles_per_tile)[i] = pti.numParticles();
+        }
+
+        return data;
+    }
+
+    void warpx_clearParticleBoundaryBuffer () {
+        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
+        particle_buffers.clearParticles();
+    }
+
     void warpx_ComputeDt () {
         WarpX& warpx = WarpX::GetInstance();
         warpx.ComputeDt ();
diff --git a/Source/Python/WarpXWrappers.h b/Source/Python/WarpXWrappers.h
index 53d461e65..71eb5b4fd 100644
--- a/Source/Python/WarpXWrappers.h
+++ b/Source/Python/WarpXWrappers.h
@@ -104,6 +104,18 @@ extern "C" {
     void warpx_addRealComp(
         const char* char_species_name, const char* char_comp_name, bool comm);
 
+    int warpx_getParticleBoundaryBufferSize(const char* species_name, int boundary);
+
+    int** warpx_getParticleBoundaryBufferScrapedSteps(
+        const char* species_name, int boundary, int lev,
+        int* num_tiles, int** particles_per_tile);
+
+    amrex::ParticleReal** warpx_getParticleBoundaryBuffer(
+        const char* species_name, int boundary, int lev,
+        int* num_tiles, int** particles_per_tile, const char* comp_name);
+
+    void warpx_clearParticleBoundaryBuffer ();
+
   void warpx_ComputeDt ();
   void warpx_MoveWindow (int step, bool move_j);
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index dcc473ff8..3b3df6790 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -98,6 +98,8 @@ public:
 
     MultiParticleContainer& GetPartContainer () { return *mypc; }
 
+    ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; }
+
     static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom,
                          int num_shift, int dir, amrex::Real external_field=0.0,
                          bool useparser = false, amrex::ParserExecutor<3> const& field_parser={});