Skip to content

Commit

Permalink
Fix bug with DSMC collisions in RZ (ECP-WarpX#5622)
Browse files Browse the repository at this point in the history
I realized there is a bug in the DSMC module for RZ geometry where the
velocity vectors for the colliding pair was not rotated so that a proper
center-of-momentum calculation could be done. This PR fixes the bug.

To check that the fix in this PR works, I compared the azimuthal
velocity distribution for energetic ions created from NBI with finite
impact parameter (such that a net rotation should be imparted on the
ions), when simulated in 3d versus RZ:

3d result:

![image](https://github.com/user-attachments/assets/458f7e9c-b7b2-46ec-b456-8733ce959f94)

RZ result:

![image](https://github.com/user-attachments/assets/ce37a51c-f127-44be-b9b9-4f9a1d7d2cbb)
  • Loading branch information
roelof-groenewald authored Feb 4, 2025
1 parent 57931b8 commit 57f6317
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 0 deletions.
1 change: 1 addition & 0 deletions Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
}
});

Expand Down

0 comments on commit 57f6317

Please sign in to comment.