Fix particle reflecting boundaries in RZ (#3837)

* Fix the reflection of velocities in RZ * Add warpx_reflect_all_velocities to the PICMI interface
ECP-WarpX · Apr 17, 2023 · 296c381 · 296c381
1 parent 71566b7
commit 296c381
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 4 deletions.
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
@@ -537,6 +537,10 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid):
 
     warpx_potential_hi_z: float, default=0.
        Electrostatic potential on the upper longitudinal boundary
+
+    warpx_reflect_all_velocities: bool default=False
+        Whether the sign of all of the particle velocities are changed upon
+        reflection on a boundary, or only the velocity normal to the surface
     """
     def init(self, kw):
         self.max_grid_size = kw.pop('warpx_max_grid_size', 32)
@@ -552,6 +556,7 @@ def init(self, kw):
         self.potential_ymax = None
         self.potential_zmin = kw.pop('warpx_potential_lo_z', None)
         self.potential_zmax = kw.pop('warpx_potential_hi_z', None)
+        self.reflect_all_velocities = kw.pop('warpx_reflect_all_velocities', None)
 
         # Geometry
         # Set these as soon as the information is available
@@ -582,6 +587,7 @@ def initialize_inputs(self):
         pywarpx.boundary.field_hi = [BC_map[bc] for bc in self.upper_boundary_conditions]
         pywarpx.boundary.particle_lo = self.lower_boundary_conditions_particles
         pywarpx.boundary.particle_hi = self.upper_boundary_conditions_particles
+        pywarpx.boundary.reflect_all_velocities = self.reflect_all_velocities
 
         if self.moving_window_velocity is not None and np.any(np.not_equal(self.moving_window_velocity, 0.)):
             pywarpx.warpx.do_moving_window = 1

diff --git a/Source/Particles/ParticleBoundaries_K.H b/Source/Particles/ParticleBoundaries_K.H
@@ -86,8 +86,11 @@ namespace ApplyParticleBoundaries {
 #ifndef WARPX_DIM_1D_Z
               amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax,
 #endif
-#ifdef WARPX_DIM_3D
-                      amrex::ParticleReal& y, amrex::Real ymin, amrex::Real ymax,
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
+                      amrex::ParticleReal& y,
+#endif
+#if (defined WARPX_DIM_3D)
+                      amrex::Real ymin, amrex::Real ymax,
 #endif
                       amrex::ParticleReal& z, amrex::Real zmin, amrex::Real zmax,
                       amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
@@ -121,8 +124,25 @@ namespace ApplyParticleBoundaries {
             change_sign_uy = true;
             change_sign_uz = true;
         }
+#ifdef WARPX_DIM_RZ
+        // Note that the reflection of the position does "r = 2*rmax - r", but this is only approximate.
+        // The exact calculation requires the position at the start of the step.
+        if (change_sign_ux && change_sign_uy) {
+            ux = -ux;
+            uy = -uy;
+        } else if (change_sign_ux) {
+            // Reflect only ur
+            // Note that y is theta
+            amrex::Real ur = ux*std::cos(y) + uy*std::sin(y);
+            amrex::Real ut = -ux*std::sin(y) + uy*std::cos(y);
+            ur = -ur;
+            ux = ur*std::cos(y) - ut*std::sin(y);
+            uy = ur*std::sin(y) + ut*std::cos(y);
+        }
+#else
         if (change_sign_ux) ux = -ux;
         if (change_sign_uy) uy = -uy;
+#endif
         if (change_sign_uz) uz = -uz;
     }
 

diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
@@ -1100,8 +1100,11 @@ WarpXParticleContainer::ApplyBoundaryConditions (){
 #ifndef WARPX_DIM_1D_Z
                                                               x, xmin, xmax,
 #endif
-#ifdef WARPX_DIM_3D
-                                                              y, ymin, ymax,
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
+                                                              y,
+#endif
+#if (defined WARPX_DIM_3D)
+                                                              ymin, ymax,
 #endif
                                                               z, zmin, zmax,
                                                               ux[i], uy[i], uz[i], particle_lost,