diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 5a3c925e9bd..6051aab1b59 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -176,6 +176,7 @@ public: m_process_count, m_scattering_processes_data, engine); #if (defined WARPX_DIM_RZ) + /* Undo the earlier velocity rotation. */ amrex::ParticleReal 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); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 473199a6b21..239a76c50c7 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -154,6 +154,25 @@ public: auto& uy2 = soa_products_data[1].m_rdata[PIdx::uy][product2_index]; auto& uz2 = soa_products_data[1].m_rdata[PIdx::uz][product2_index]; +#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.) + */ + amrex::ParticleReal const theta = ( + soa_products_data[1].m_rdata[PIdx::theta][product2_index] + - soa_products_data[0].m_rdata[PIdx::theta][product1_index] + ); + amrex::ParticleReal const ux1buf = ux1; + ux1 = ux1buf*std::cos(theta) - uy1*std::sin(theta); + uy1 = ux1buf*std::sin(theta) + uy1*std::cos(theta); +#endif + // for simplicity (for now) we assume non-relativistic particles // and simply calculate the center-of-momentum velocity from the // rest masses @@ -213,6 +232,13 @@ public: ux2 += uCOM_x; uy2 += uCOM_y; uz2 += uCOM_z; + +#if (defined WARPX_DIM_RZ) + /* Undo the earlier velocity rotation. */ + amrex::ParticleReal const ux1buf_new = ux1; + ux1 = ux1buf_new*std::cos(-theta) - uy1*std::sin(-theta); + uy1 = ux1buf_new*std::sin(-theta) + uy1*std::cos(-theta); +#endif } });