Reflective particle boundary condition (#1728)

* Added ParticleBoundaries and reflecting boundary conditions

* Added ParticleBoundaries::AllNone

* Allowed different particle boundary conditions on each side of the domain

* Updated the documentation for particle boundaries

* Fix end of line space in Docs/source/running_cpp/parameters.rst

* Updated the reflecting BC to use boundary input group

* Fixes to reflective boundary conditions

* Bug fix in AsStored

* Added particle boundaries regression test particle_boundaries_3d

* Fixed particle_boundaries_3d.json

* Minor updates

* Added algo.particle_shape to test case

* Remove do_pml from test case

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

* Need to explicitly turn off pml in CI test

* Re-add include

* Fixed includes

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Revathi  Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

1 files changed, 25 insertions, 9 deletions

diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 8c5ad7f7f..b5ed32c94 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -1090,13 +1090,17 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p,
 }
 
 void
-WarpXParticleContainer::ApplyBoundaryConditions (ParticleBC boundary_conditions){
+WarpXParticleContainer::ApplyBoundaryConditions (ParticleBoundaries& boundary_conditions){
     WARPX_PROFILE("WarpXParticleContainer::ApplyBoundaryConditions()");
+
+    if (boundary_conditions.CheckAll(ParticleBoundaryType::Periodic)) return;
+
     for (int lev = 0; lev <= finestLevel(); ++lev)
     {
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
             auto GetPosition = GetParticlePosition(pti);
+            auto SetPosition = SetParticlePosition(pti);
             const Real xmin = Geom(lev).ProbLo(0);
             const Real xmax = Geom(lev).ProbHi(0);
 #ifdef WARPX_DIM_3D
@@ -1109,22 +1113,34 @@ WarpXParticleContainer::ApplyBoundaryConditions (ParticleBC boundary_conditions)
             ParticleTileType& ptile = ParticlesAt(lev, pti);
             ParticleType * const pp = ptile.GetArrayOfStructs()().data();
 
+            auto& soa = ptile.GetStructOfArrays();
+            amrex::ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data();
+            amrex::ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data();
+            amrex::ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data();
+
             // Loop over particles and apply BC to each particle
             amrex::ParallelFor(
                 pti.numParticles(),
                 [=] AMREX_GPU_DEVICE (long i) {
                     ParticleType& p = pp[i];
                     ParticleReal x, y, z;
-                    GetPosition(i, x, y, z);
+                    GetPosition.AsStored(i, x, y, z);
+                    // Note that for RZ, (x, y, z) is actually (r, theta, z).
+
+                    bool particle_lost = false;
+                    ParticleBoundaries::apply_boundaries(x, xmin, xmax,
 #ifdef WARPX_DIM_3D
-                    if (x < xmin || x > xmax || y < ymin || y > ymax || z < zmin || z > zmax){
-                        if (boundary_conditions == ParticleBC::absorbing) p.id() = -1;
-                    }
-#else
-                    if (x < xmin || x > xmax || z < zmin || z > zmax){
-                        if (boundary_conditions == ParticleBC::absorbing) p.id() = -1;
-                    }
+                                                         y, ymin, ymax,
 #endif
+                                                         z, zmin, zmax,
+                                                         ux[i], uy[i], uz[i], particle_lost,
+                                                         boundary_conditions);
+
+                    if (particle_lost) {
+                        p.id() = -1;
+                    } else {
+                        SetPosition.AsStored(i, x, y, z);
+                    }
                 }
             );
         }