diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index f5c8e67640f..3ff582e536f 100644 --- a/Docs/source/usage/parameters.rst +++ b/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) `_. +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 diff --git a/Examples/Tests/collision/analysis_collision_rz.py b/Examples/Tests/collision/analysis_collision_rz.py new file mode 100755 index 00000000000..82b2e6da7e5 --- /dev/null +++ b/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) diff --git a/Examples/Tests/collision/inputs_rz b/Examples/Tests/collision/inputs_rz new file mode 100644 index 00000000000..687b24bd6c2 --- /dev/null +++ b/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 diff --git a/Regression/Checksum/benchmarks_json/collisionRZ.json b/Regression/Checksum/benchmarks_json/collisionRZ.json new file mode 100644 index 00000000000..d2493847113 --- /dev/null +++ b/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 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index f0403f38f89..a0ffe344a39 100644 --- a/Regression/WarpX-tests.ini +++ b/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 diff --git a/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H b/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H index fd74d3e87e8..5e94a0f93a9 100644 --- a/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H +++ b/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H @@ -107,11 +107,32 @@ 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] ], @@ -119,6 +140,13 @@ void ElasticCollisionPerez ( 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(I1e) ) { i1 = I1s; } ++i2; if ( i2 == static_cast(I2e) ) { i2 = I2s; } } diff --git a/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H index b5acd85a706..38ddc548182 100644 --- a/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H +++ b/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::min() && - amrex::Math::abs(u1y-u2y) < std::numeric_limits::min() && - amrex::Math::abs(u1z-u2z) < std::numeric_limits::min() ) - { return; } + // Or if the relative difference is less than 1.0e-10. + if ( diffm < std::numeric_limits::min() || diffm/summm < 1.0e-10 ) { return; } T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp index 07309257449..dd78f160fc7 100644 --- a/Source/Particles/Collision/CollisionHandler.cpp +++ b/Source/Particles/Collision/CollisionHandler.cpp @@ -30,6 +30,9 @@ CollisionHandler::CollisionHandler(MultiParticleContainer const * const mypc) for (int i = 0; i < static_cast(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";