Skip to content

Commit

Permalink
Refactored binary collision method. Density and temperature are only …
Browse files Browse the repository at this point in the history 
…computed once per cell, if need, rather than at each binary pairing (order Npcc^2 cost). Coulomb method now produces correct scattering physics with weighted particles.
  • Loading branch information
@JustinRayAngus
JustinRayAngus committed Jul 12, 2024
1 parent dc087bf commit 1696315
Show file tree
Hide file tree
Showing 6 changed files with 214 additions and 108 deletions.
115 changes: 115 additions & 0 deletions Source/Particles/Collision/BinaryCollision/BinaryCollision.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#define WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_

#include "Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H"
#include "Particles/Collision/BinaryCollision/Coulomb/ComputeTemperature.H"
#include "Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H"
#include "Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H"
#include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H"
Expand Down Expand Up @@ -360,6 +361,18 @@ public:
End of calculations only required when creating product particles
*/

// create vectors to store density and temperature on cell level
amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec;
amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec;
if (binary_collision_functor.m_computeSpeciesDensities) {
n1_vec.resize(n_cells);
}
if (binary_collision_functor.m_computeSpeciesTemperatures) {
T1_vec.resize(n_cells);
}
amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();

// Loop over cells
amrex::ParallelForRNG( n_cells,
[=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
Expand All @@ -372,6 +385,30 @@ public:
// Do not collide if there is only one particle in the cell
if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }

// compute local density [1/m^3]
if (binary_collision_functor.m_computeSpeciesDensities) {
amrex::ParticleReal wtot1 = 0.0;
amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
wtot1 += w1[ indices_1[i1] ];
}
#if defined WARPX_DIM_RZ
const int ri = (i_cell - i_cell%nz) / nz;
auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
#endif
n1_in_each_cell[i_cell] = wtot1/dV;
}

// compute local temperature [Joules]
if (binary_collision_functor.m_computeSpeciesTemperatures) {
amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
T1_in_each_cell[i_cell] = ComputeTemperature( cell_start_1, cell_stop_1, indices_1,
w1, u1x, u1y, u1z, m1 );
}

// shuffle
ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
}
Expand Down Expand Up @@ -410,6 +447,17 @@ public:
const int ri = (i_cell - i_cell%nz) / nz;
auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
#endif

// Get the local density and temperature for this cell
amrex::ParticleReal n1 = 0.0;
amrex::ParticleReal T1 = 0.0;
if (binary_collision_functor.m_computeSpeciesDensities) {
n1 = n1_in_each_cell[i_cell];
}
if (binary_collision_functor.m_computeSpeciesTemperatures) {
T1 = T1_in_each_cell[i_cell];
}

// Call the function in order to perform collisions
// If there are product species, mask, p_pair_indices_1/2, and
// p_pair_reaction_weight are filled here
Expand All @@ -418,6 +466,7 @@ public:
cell_half_1, cell_stop_1,
indices_1, indices_1,
soa_1, soa_1, get_position_1, get_position_1,
n1, n1, T1, T1,
q1, q1, m1, m1, dt, dV, coll_idx,
cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
p_pair_reaction_weight, engine);
Expand Down Expand Up @@ -562,6 +611,21 @@ public:
End of calculations only required when creating product particles
*/

// create vectors to store density and temperature on cell level
amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec, n2_vec;
amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec, T2_vec;
if (binary_collision_functor.m_computeSpeciesDensities) {
n1_vec.resize(n_cells);
n2_vec.resize(n_cells);
}
if (binary_collision_functor.m_computeSpeciesTemperatures) {
T1_vec.resize(n_cells);
T2_vec.resize(n_cells);
}
amrex::Real* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
amrex::Real* AMREX_RESTRICT n2_in_each_cell = n2_vec.dataPtr();
amrex::Real* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();
amrex::Real* AMREX_RESTRICT T2_in_each_cell = T2_vec.dataPtr();

// Loop over cells
amrex::ParallelForRNG( n_cells,
Expand All @@ -579,6 +643,43 @@ public:
// ux_1[ indices_1[i] ], where i is between
// cell_start_1 (inclusive) and cell_start_2 (exclusive)

// compute local densities [1/m^3]
if (binary_collision_functor.m_computeSpeciesDensities) {
amrex::ParticleReal w1tot = 0.0;
amrex::ParticleReal w2tot = 0.0;
amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
w1tot += w1[ indices_1[i1] ];
}
for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) {
w2tot += w2[ indices_2[i2] ];
}
#if defined WARPX_DIM_RZ
const int ri = (i_cell - i_cell%nz) / nz;
auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
#endif
n1_in_each_cell[i_cell] = w1tot/dV;
n2_in_each_cell[i_cell] = w2tot/dV;
}

// compute local temperatures [Joules]
if (binary_collision_functor.m_computeSpeciesTemperatures) {
amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
T1_in_each_cell[i_cell] = ComputeTemperature( cell_start_1, cell_stop_1, indices_1,
w1, u1x, u1y, u1z, m1 );

amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
T2_in_each_cell[i_cell] = ComputeTemperature( cell_start_2, cell_stop_2, indices_2,
w2, u2x, u2y, u2z, m2 );
}

// Do not collide if one species is missing in the cell
if ( cell_stop_1 - cell_start_1 < 1 ||
cell_stop_2 - cell_start_2 < 1 ) { return; }
Expand Down Expand Up @@ -628,13 +729,27 @@ public:
const int ri = (i_cell - i_cell%nz) / nz;
auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
#endif

// Get the local densities and temperatures for this cell
amrex::ParticleReal n1 = 0.0, n2 = 0.0;
amrex::ParticleReal T1 = 0.0, T2 = 0.0;
if (binary_collision_functor.m_computeSpeciesDensities) {
n1 = n1_in_each_cell[i_cell];
n2 = n2_in_each_cell[i_cell];
}
if (binary_collision_functor.m_computeSpeciesTemperatures) {
T1 = T1_in_each_cell[i_cell];
T2 = T2_in_each_cell[i_cell];
}

// Call the function in order to perform collisions
// If there are product species, p_mask, p_pair_indices_1/2, and
// p_pair_reaction_weight are filled here
binary_collision_functor(
cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
indices_1, indices_2,
soa_1, soa_2, get_position_1, get_position_2,
n1, n2, T1, T2,
q1, q2, m1, m2, dt, dV, coll_idx,
cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
p_pair_reaction_weight, engine);
Expand Down
90 changes: 31 additions & 59 deletions Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_

#include "ComputeTemperature.H"
#include "UpdateMomentumPerezElastic.H"
#include "Particles/WarpXParticleContainer.H"
#include "Utils/WarpXConst.H"
Expand All @@ -25,12 +24,11 @@
* @param[in] I1e,I2e is the stop index for I1,I2 (exclusive).
* @param[in] I1,I2 the index arrays. They determine all elements that will be used.
* @param[in,out] soa_1,soa_2 the struct of array for species 1/2
* @param[in] n1,n2 density of species 1/2
* @param[in] T1,T2 temperature (Joules) of species 1/2
* only used if L <= 0
* @param[in] q1,q2 charge of species 1/2
* @param[in] m1,m2 mass of species 1/2
* @param[in] T1 temperature (Joule) of species 1
* and will be used if greater than zero,
* otherwise will be computed.
* @param[in] T2 temperature (Joule) of species 2, @see T1
* @param[in] dt is the time step length between two collision calls.
* @param[in] L is the Coulomb log and will be used if greater than zero,
* otherwise will be computed.
Expand All @@ -48,10 +46,11 @@ void ElasticCollisionPerez (
T_index const* AMREX_RESTRICT I1,
T_index const* AMREX_RESTRICT I2,
SoaData_type soa_1, SoaData_type soa_2,
T_PR const n1, T_PR const n2,
T_PR const T1, T_PR const T2,
T_PR const q1, T_PR const q2,
T_PR const m1, T_PR const m2,
T_PR const T1, T_PR const T2,
T_R const dt, T_PR const L, T_R const dV,
T_R const dt, T_PR const L, T_R const dV,
amrex::RandomEngine const& engine,
bool const isSameSpecies, T_index coll_idx)
{
Expand All @@ -70,57 +69,22 @@ void ElasticCollisionPerez (
T_PR * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
T_PR * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];

// get local T1t and T2t
T_PR T1t; T_PR T2t;
if ( T1 <= T_PR(0.0) && L <= T_PR(0.0) )
{
T1t = ComputeTemperature(I1s,I1e,I1,w1,u1x,u1y,u1z,m1);
}
else { T1t = T1; }
if ( T2 <= T_PR(0.0) && L <= T_PR(0.0) )
{
T2t = ComputeTemperature(I2s,I2e,I2,w2,u2x,u2y,u2z,m2);
}
else { T2t = T2; }

// local density
T_PR n1 = T_PR(0.0);
T_PR n2 = T_PR(0.0);
T_PR n12 = T_PR(0.0);
for (T_index i1=I1s; i1<I1e; ++i1) { n1 += w1[ I1[i1] ]; }
for (T_index i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; }
// Intra-species: the density is in fact the sum of the density of
// each sub-group of particles
if (isSameSpecies) {
n1 = n1 + n2;
n2 = n1;
}
n1 = n1 / dV; n2 = n2 / dV;
{
T_index i1 = I1s; T_index i2 = I2s;
for (T_index k = 0; k < max_N; ++k)
{
n12 += amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] );
++i1; if ( i1 == I1e ) { i1 = I1s; }
++i2; if ( i2 == I2e ) { i2 = I2s; }
}
n12 = n12 / dV;
// compute Debye length lmdD (if not using a fixed L = Coulomb Log)
T_PR lmdD = T_PR(-1.0);
if ( L < T_PR(0.0) ) {
lmdD = T_PR(1.0)/std::sqrt( n1*q1*q1/(T1*PhysConst::ep0) +
n2*q2*q2/(T2*PhysConst::ep0) );
}
// Intra-species: Apply correction in collision rate
if (isSameSpecies) { n12 *= T_PR(2.0); }

// compute Debye length lmdD
T_PR lmdD;
if ( T1t < T_PR(0.0) || T2t < T_PR(0.0) ) {
lmdD = T_PR(0.0);
}
else {
lmdD = T_PR(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) +
n2*q2*q2/(T2t*PhysConst::ep0) );
}
T_PR rmin = std::pow( T_PR(4.0) * MathConst::pi / T_PR(3.0) *
amrex::max(n1,n2), T_PR(-1.0/3.0) );
lmdD = amrex::max(lmdD, rmin);
// compute atomic spacing
const T_PR maxn = amrex::max(n1,n2);
const auto rmin = static_cast<T_PR>( 1.0/std::cbrt(4.0*MathConst::pi/3.0*maxn) );

// bmax (i.e., screening length) cannot be smaller than atomic spacing
const T_PR bmax = amrex::max(lmdD, rmin);

// set max cross section based on mfp = atomic spacing
const T_PR sigma_max = T_PR(1.0) / (maxn * rmin);

#if (defined WARPX_DIM_RZ)
T_PR * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
Expand Down Expand Up @@ -151,13 +115,21 @@ void ElasticCollisionPerez (
u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
#endif

// compute the effective density used for s12 in UpdateMomentumPerezElastic()
// The effective density used for s12 is defined such that average value of s12
// after many scattering cycles is the same as that for the full NxN pairing method.
// The approach here is similar to that in JCP 413 (2020) by Higginson et al., but
// does not use a unique duplication factor for each macro-particle pair.
T_PR n12;
const T_PR wpmax = amrex::max(w1[ I1[i1] ],w2[ I2[i2] ]);
if (isSameSpecies) { n12 = wpmax*static_cast<T_PR>(min_N+max_N-1)/dV; }
else { n12 = wpmax*static_cast<T_PR>(min_N)/dV; }

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);
n12, sigma_max, L, bmax, dt, engine );

#if (defined WARPX_DIM_RZ)
T_PR const u1xbuf_new = u1x[I1[i1]];
Expand Down
11 changes: 9 additions & 2 deletions Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ public:
m_CoulombLog = CoulombLog;

m_exe.m_CoulombLog = m_CoulombLog;
if (m_CoulombLog<0.0) { m_exe.m_computeSpeciesTemperatures = true; }
m_exe.m_isSameSpecies = m_isSameSpecies;
}

Expand All @@ -72,6 +73,8 @@ public:
* @param[in] I1e,I2e is the stop index for I1,I2 (exclusive).
* @param[in] I1,I2 index arrays. They determine all elements that will be used.
* @param[in,out] soa_1,soa_2 contain the struct of array data of the two species.
* @param[in] n1,n2 are local densities.
* @param[in] T1,T2 are local temperatures.
* @param[in] q1,q2 are charges.
* @param[in] m1,m2 are masses.
* @param[in] dt is the time step length between two collision calls.
Expand All @@ -87,6 +90,8 @@ public:
index_type const* AMREX_RESTRICT I2,
const SoaData_type& soa_1, const SoaData_type& soa_2,
GetParticlePosition<PIdx> /*get_position_1*/, GetParticlePosition<PIdx> /*get_position_2*/,
amrex::ParticleReal const n1, amrex::ParticleReal const n2,
amrex::ParticleReal const T1, amrex::ParticleReal const T2,
amrex::ParticleReal const q1, amrex::ParticleReal const q2,
amrex::ParticleReal const m1, amrex::ParticleReal const m2,
amrex::Real const dt, amrex::Real const dV, index_type coll_idx,
Expand All @@ -99,12 +104,14 @@ public:

ElasticCollisionPerez(
I1s, I1e, I2s, I2e, I1, I2,
soa_1, soa_2,
q1, q2, m1, m2, -1.0_prt, -1.0_prt,
soa_1, soa_2, n1, n2, T1, T2,
q1, q2, m1, m2,
dt, m_CoulombLog, dV, engine, m_isSameSpecies, coll_idx);
}

amrex::ParticleReal m_CoulombLog;
bool m_computeSpeciesDensities = true;
bool m_computeSpeciesTemperatures = false;
bool m_isSameSpecies;
};

Expand Down
Loading

0 comments on commit 1696315

Please sign in to comment.