Skip to content

Commit

Permalink
moved Temp calc from inside ElasticCollisionPerez to cell level in Bi…
Browse files Browse the repository at this point in the history 
…naryCollision.H. Cost is now order N rather than order N^2. Also generalized ComputeTemperature.H to work for weighted particles.
  • Loading branch information
@JustinRayAngus
JustinRayAngus committed Jul 9, 2024
1 parent 372cf91 commit a9a9e2d
Show file tree
Hide file tree
Showing 6 changed files with 81 additions and 50 deletions.
69 changes: 52 additions & 17 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 @@ -373,15 +374,27 @@ 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
// compute local density [1/m^3]
n1 = 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) { n1 += w1[ indices_1[i1] ]; }
if (binary_collision_functor.m_computeSpeciesDensities) {
amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) { n1 += 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;
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 /= dV;
n1 /= dV;
}

// compute local temperature [Joules]
T1 = 0.0;
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 = ComputeTemperature(cell_start_1,cell_stop_1,indices_1,w1,u1x,u1y,u1z,m1);
}

// shuffle
ShuffleFisherYates(
Expand Down Expand Up @@ -430,7 +443,8 @@ public:
cell_half_1, cell_stop_1,
indices_1, indices_1,
soa_1, soa_1, get_position_1, get_position_1,
n1, n1, q1, q1, m1, m1, dt, dV, coll_idx,
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 @@ -591,19 +605,38 @@ public:
// ux_1[ indices_1[i] ], where i is between
// cell_start_1 (inclusive) and cell_start_2 (exclusive)

// compute local densities
// compute local densities [1/m^3]
n1 = 0.0;
n2 = 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) { n1 += w1[ indices_1[i1] ]; }
for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) { n2 += w2[ indices_2[i2] ]; }
if (binary_collision_functor.m_computeSpeciesDensities) {
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) { n1 += w1[ indices_1[i1] ]; }
for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) { n2 += 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;
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 /= dV;
n2 /= dV;
n1 /= dV;
n2 /= dV;
}

// compute local temperatures [Joules]
T1 = 0.0;
T2 = 0.0;
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 = 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 = 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 ||
Expand Down Expand Up @@ -661,7 +694,8 @@ public:
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, q1, q2, m1, m2, dt, dV, coll_idx,
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 Expand Up @@ -693,6 +727,7 @@ private:
bool m_isSameSpecies;
bool m_have_product_species;
mutable amrex::ParticleReal n1, n2; // local densities
mutable amrex::ParticleReal T1, T2; // local temperatures
amrex::Vector<std::string> m_product_species;
// functor that performs collisions within a cell
CollisionFunctor m_binary_collision_functor;
Expand Down
15 changes: 9 additions & 6 deletions Source/Particles/Collision/BinaryCollision/Coulomb/ComputeTemperature.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ template <typename T_index, typename T_R>
AMREX_GPU_HOST_DEVICE
T_R ComputeTemperature (
T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I,
T_R const * AMREX_RESTRICT w,
T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT uz,
T_R const m )
{
Expand All @@ -26,21 +27,23 @@ T_R ComputeTemperature (
T_R vx = T_R(0.0); T_R vy = T_R(0.0);
T_R vz = T_R(0.0); T_R vs = T_R(0.0);
T_R gm = T_R(0.0); T_R us = T_R(0.0);
T_R wtot = T_R(0.0);

for (int i = Is; i < static_cast<int>(Ie); ++i)
{
us = ( ux[ I[i] ] * ux[ I[i] ] +
uy[ I[i] ] * uy[ I[i] ] +
uz[ I[i] ] * uz[ I[i] ] );
gm = std::sqrt( T_R(1.0) + us*inv_c2 );
vx += ux[ I[i] ] / gm;
vy += uy[ I[i] ] / gm;
vz += uz[ I[i] ] / gm;
vs += us / gm / gm;
wtot += w[ I[i] ];
vx += w[ I[i] ]*ux[ I[i] ] / gm;
vy += w[ I[i] ]*uy[ I[i] ] / gm;
vz += w[ I[i] ]*uz[ I[i] ] / gm;
vs += w[ I[i] ]*us / gm / gm;
}

vx = vx / N; vy = vy / N;
vz = vz / N; vs = vs / N;
vx = vx / wtot; vy = vy / wtot;
vz = vz / wtot; vs = vs / wtot;

return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
}
Expand Down
32 changes: 7 additions & 25 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 @@ -29,8 +28,7 @@
* @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.
* and will be used if L is less than zero,
* @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,
Expand All @@ -50,9 +48,9 @@ void ElasticCollisionPerez (
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,
amrex::RandomEngine const& engine,
bool const isSameSpecies, T_index coll_idx)
Expand All @@ -72,27 +70,11 @@ 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,u1x,u1y,u1z,m1);
}
else { T1t = T1; }
if ( T2 <= T_PR(0.0) && L <= T_PR(0.0) )
{
T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2);
}
else { T2t = T2; }

// 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) );
// 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) );
}

// compute atomic spacing
Expand Down
9 changes: 7 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 @@ -73,6 +74,7 @@ public:
* @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 @@ -89,6 +91,7 @@ public:
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 @@ -101,12 +104,14 @@ public:

ElasticCollisionPerez(
I1s, I1e, I2s, I2e, I1, I2,
soa_1, soa_2,
n1, n2, 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
3 changes: 3 additions & 0 deletions Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ public:
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 Down Expand Up @@ -190,6 +191,8 @@ public:
}

int m_process_count;
bool m_computeSpeciesDensities = false;
bool m_computeSpeciesTemperatures = false;
ScatteringProcess::Executor* m_scattering_processes_data;
};

Expand Down
3 changes: 3 additions & 0 deletions Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,7 @@ public:
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 Down Expand Up @@ -233,6 +234,8 @@ public:
amrex::ParticleReal m_probability_threshold;
amrex::ParticleReal m_probability_target_value;
NuclearFusionType m_fusion_type;
bool m_computeSpeciesDensities = false;
bool m_computeSpeciesTemperatures = false;
bool m_isSameSpecies;
};

Expand Down

0 comments on commit a9a9e2d

Please sign in to comment.