Stochastic particle reflection from absorbing domain boundaries (#2281)

* added tunable particle reflection from absorbing domain boundaries

* extended picmi.py to allow setting boundary reflection coefficients and added a CI test for the reflection implementation

* allow R(E) to be specified, except for embedded boundaries

* changed approach for particle reflection, now the ParticleBoundaries object will hold the reflection coefficient; reflection from EBs not implemented

* added functionality to reflect from EB; still needs to be tested for accuracy

* added support for energy dependent reflection models for domain boundaries

* fixed at least one issue causing CI fails - building reflection model parsers for not physical particle containers

* switched reflection coefficients to be functions of the velocity component perpendicular to the boundary rather than energy

* reverted initial work on reflecting from EBs

* changed naming convention for new CI test for particle reflection

* switched useMPI back to 1 in test

* breaking changes while trying to sort out GPU issue

* fixed issue with CUDA compilation - hopefully :)

* various code improvements from PR review suggestions

* fix of major issues

* no need to parse the reflection models at every step

* skip particles that are already flagged for removal in ApplyBoundaryConditions

1 files changed, 25 insertions, 6 deletions

diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 43b8e5fe6..5a76a24df 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -104,6 +104,18 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies)
     local_jx.resize(num_threads);
     local_jy.resize(num_threads);
     local_jz.resize(num_threads);
+
+    // The boundary conditions are read in in ReadBCParams but a child class
+    // can allow these value to be overwritten if different boundary
+    // conditions are desired for a specific species
+    m_boundary_conditions.SetBoundsX(WarpX::particle_boundary_lo[0], WarpX::particle_boundary_hi[0]);
+#ifdef WARPX_DIM_3D
+    m_boundary_conditions.SetBoundsY(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]);
+    m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[2], WarpX::particle_boundary_hi[2]);
+#else
+    m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]);
+#endif
+    m_boundary_conditions.BuildReflectionModelParsers();
 }
 
 void
@@ -1099,10 +1111,13 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p,
 }
 
 void
-WarpXParticleContainer::ApplyBoundaryConditions (ParticleBoundaries& boundary_conditions){
+WarpXParticleContainer::ApplyBoundaryConditions (){
     WARPX_PROFILE("WarpXParticleContainer::ApplyBoundaryConditions()");
 
-    if (boundary_conditions.CheckAll(ParticleBoundaryType::Periodic)) return;
+    // Periodic boundaries are handled in AMReX code
+    if (m_boundary_conditions.CheckAll(ParticleBoundaryType::Periodic)) return;
+
+    auto boundary_conditions = m_boundary_conditions.data;
 
     for (int lev = 0; lev <= finestLevel(); ++lev)
     {
@@ -1128,10 +1143,14 @@ WarpXParticleContainer::ApplyBoundaryConditions (ParticleBoundaries& boundary_co
             amrex::ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data();
 
             // Loop over particles and apply BC to each particle
-            amrex::ParallelFor(
+            amrex::ParallelForRNG(
                 pti.numParticles(),
-                [=] AMREX_GPU_DEVICE (long i) {
+                [=] AMREX_GPU_DEVICE (long i, amrex::RandomEngine const& engine) {
                     ParticleType& p = pp[i];
+
+                    // skip particles that are already flagged for removal
+                    if (p.id() < 0) return;
+
                     ParticleReal x, y, z;
                     GetPosition.AsStored(i, x, y, z);
                     // Note that for RZ, (x, y, z) is actually (r, theta, z).
@@ -1143,10 +1162,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (ParticleBoundaries& boundary_co
 #endif
                                                               z, zmin, zmax,
                                                               ux[i], uy[i], uz[i], particle_lost,
-                                                              boundary_conditions);
+                                                              boundary_conditions, engine);
 
                     if (particle_lost) {
-                        p.id() = -1;
+                        p.id() = -p.id();
                     } else {
                         SetPosition.AsStored(i, x, y, z);
                     }