diff options
Diffstat (limited to 'Source')
3 files changed, 38 insertions, 4 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H b/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H index fd74d3e87..5e94a0f93 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<int>(I1e) ) { i1 = I1s; } ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; } } diff --git a/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H index b5acd85a7..38ddc5481 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<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 ); diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp index 073092574..dd78f160f 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<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"; |