Correct collision algorithm in RZ.

Docs/source/usage/parameters.rst

@@ -1303,6 +1303,8 @@ Collision initialization
 
 WarpX provides a relativistic elastic Monte Carlo binary collision model,
 following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`_.
+When the RZ mode is used, `warpx.n_rz_azimuthal_modes` must be set to 1 at the moment,
+since the current implementation of the collision module assumes axisymmetry.
 A non-relativistic Monte Carlo treatment for particles colliding
 with a neutral, uniform background gas is also available. The implementation follows the so-called
 null collision strategy discussed for example in `Birdsall (IEEE Transactions on

Examples/Tests/collision/analysis_collision_rz.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+
+# Copyright 2019-2021 Yinjian Zhao
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This script tests the collision module in RZ.
+# Electrons are set with the same vr,
+# and thus particles are traveling locally in the same direction.
+# Physically, the collision rate should therefore be very low.
+# Thus the initial electron momentum will be compared with the final momentum.
+
+# Possible errors:
+# tolerance: 1.0e-30
+# Possible running time: ~ 1.0 s
+
+import os
+import sys
+import yt
+import numpy as np
+from glob import glob
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+tolerance = 1.0e-15
+
+last_fn = sys.argv[1]
+if (last_fn[-1] == "/"): last_fn = last_fn[:-1]
+fn_list = glob(last_fn[:-5] + "?????")
+
+for fn in fn_list:
+    # get time index j
+    j = int(fn[-5:])
+    if j==0:
+        # load file
+        ds = yt.load( fn )
+        ad = ds.all_data()
+        px1 = ad['particle_momentum_x'].to_ndarray()
+        py1 = ad['particle_momentum_y'].to_ndarray()
+    if j==150:
+        # load file
+        ds = yt.load( fn )
+        ad = ds.all_data()
+        px2 = ad['particle_momentum_x'].to_ndarray()
+        py2 = ad['particle_momentum_y'].to_ndarray()
+
+error = np.max( abs(px1-px2)+abs(py1-py2) )
+
+print('error = ', error)
+print('tolerance = ', tolerance)
+assert(error < tolerance)
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn, do_particles=False)

Examples/Tests/collision/inputs_rz

@@ -0,0 +1,56 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 150
+amr.n_cell = 8 8
+amr.max_grid_size = 8
+amr.blocking_factor = 8
+amr.max_level = 0
+geometry.coord_sys   = 1
+geometry.prob_lo     = 0.    0.
+geometry.prob_hi     = 4.154046151855669e2  4.154046151855669e2
+
+#################################
+###### Boundary Condition #######
+#################################
+boundary.field_lo = none periodic
+boundary.field_hi = none periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_ics = 1
+warpx.verbose = 1
+warpx.cfl = 1.0
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = electron
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "NUniformPerCell"
+electron.num_particles_per_cell_each_dim = 1 10 1
+electron.profile = constant
+electron.density = 1.0e21
+electron.momentum_distribution_type = parse_momentum_function
+electron.momentum_function_ux(x,y,z) = "if(x*x+y*y>0.0, 1.0*x/sqrt(x*x+y*y), 0.0)"
+electron.momentum_function_uy(x,y,z) = "if(x*x+y*y>0.0, 1.0*y/sqrt(x*x+y*y), 0.0)"
+electron.momentum_function_uz(x,y,z) = "0"
+electron.do_not_deposit = 1
+electron.do_not_push = 1
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = collision
+collision.species = electron electron
+collision.CoulombLog = 15.9
+
+diagnostics.diags_names = diag1
+diag1.intervals = 10
+diag1.diag_type = Full

Regression/Checksum/benchmarks_json/collisionRZ.json

@@ -0,0 +1,13 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 0.0,
+    "Ey": 0.0,
+    "Ez": 0.0,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  }
+}

Regression/WarpX-tests.ini

@@ -1788,6 +1788,24 @@ analysisRoutine = Examples/Tests/collision/analysis_collision_2d.py
 aux1File = Regression/PostProcessingUtils/post_processing_utils.py
 tolerance = 1.e-14
 
+[collisionRZ]
+buildDir = .
+inputFile = Examples/Tests/collision/inputs_rz
+runtime_params =
+dim = 2
+addToCompileString = USE_RZ=TRUE USE_PSATD=FALSE
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/collision/analysis_collision_rz.py
+aux1File = Regression/PostProcessingUtils/post_processing_utils.py
+tolerance = 1.e-15
+
 [Maxwell_Hybrid_QED_solver]
 buildDir = .
 inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d

Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H

@@ -107,18 +107,46 @@ void ElasticCollisionPerez (
                amrex::max(n1,n2), T_R(-1.0/3.0) );
     lmdD = amrex::max(lmdD, rmin);
 
+#if (defined WARPX_DIM_RZ)
+    T_R * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
+    T_R * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
+#endif
+
     // call UpdateMomentumPerezElastic()
     {
       int i1 = I1s; int i2 = I2s;
       for (int k = 0; k < amrex::max(NI1,NI2); ++k)
       {
+
+#if (defined WARPX_DIM_RZ)
+          /* In RZ geometry, macroparticles can collide with other macroparticles
+           * in the same *cylindrical* cell. For this reason, collisions between macroparticles
+           * are actually not local in space. In this case, the underlying assumption is that
+           * particles within the same cylindrical cell represent a cylindrically-symmetry
+           * momentum distribution function. Therefore, here, we temporarily rotate the
+           * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
+           * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
+           * there is a corresponding assert statement at initialization.) */
+          T_R const theta = theta2[I2[i2]]-theta1[I1[i1]];
+          T_R const u1xbuf = u1x[I1[i1]];
+          u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
+          u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
+#endif
+
           UpdateMomentumPerezElastic(
               u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
               u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
               n1, n2, n12,
               q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ],
               dt, L, lmdD,
               engine);
+
+#if (defined WARPX_DIM_RZ)
+          T_R const u1xbuf_new = u1x[I1[i1]];
+          u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
+          u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
+#endif
+
           ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; }
           ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; }
       }

Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H

@@ -36,11 +36,14 @@ void UpdateMomentumPerezElastic (
     amrex::RandomEngine const& engine)
 {
 
+    T_R const diffx = amrex::Math::abs(u1x-u2x);
+    T_R const diffy = amrex::Math::abs(u1y-u2y);
+    T_R const diffz = amrex::Math::abs(u1z-u2z);
+    T_R const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz);
+    T_R const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z);
     // If g = u1 - u2 = 0, do not collide.
-    if ( amrex::Math::abs(u1x-u2x) < std::numeric_limits<T_R>::min() &&
-         amrex::Math::abs(u1y-u2y) < std::numeric_limits<T_R>::min() &&
-         amrex::Math::abs(u1z-u2z) < std::numeric_limits<T_R>::min() )
-    { return; }
+    // Or if the relative difference is less than 1.0e-10.
+    if ( diffm < std::numeric_limits<T_R>::min() || diffm/summm < 1.0e-10 ) { return; }
 
     T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
 

Source/Particles/Collision/CollisionHandler.cpp

@@ -30,6 +30,9 @@ CollisionHandler::CollisionHandler(MultiParticleContainer const * const mypc)
     for (int i = 0; i < static_cast<int>(ncollisions); ++i) {
         amrex::ParmParse pp_collision_name(collision_names[i]);
 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::n_rz_azimuthal_modes==1,
+        "RZ mode `warpx.n_rz_azimuthal_modes` must be 1 when using the binary collision module.");
+
         // For legacy, pairwisecoulomb is the default
         std::string type = "pairwisecoulomb";