aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/Collision
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/Collision')
-rw-r--r--Source/Particles/Collision/CollisionType.H39
-rw-r--r--Source/Particles/Collision/CollisionType.cpp241
-rw-r--r--Source/Particles/Collision/ComputeTemperature.H39
-rw-r--r--Source/Particles/Collision/ElasticCollisionPerez.H105
-rw-r--r--Source/Particles/Collision/Make.package10
-rw-r--r--Source/Particles/Collision/ShuffleFisherYates.H29
-rw-r--r--Source/Particles/Collision/UpdateMomentumPerezElastic.H252
7 files changed, 715 insertions, 0 deletions
diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H
new file mode 100644
index 000000000..d020f47e8
--- /dev/null
+++ b/Source/Particles/Collision/CollisionType.H
@@ -0,0 +1,39 @@
+#ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
+#define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
+
+#include "WarpXParticleContainer.H"
+#include <AMReX_DenseBins.H>
+#include <AMReX_REAL.H>
+#include <AMReX_ParmParse.H>
+
+class CollisionType
+{
+public:
+ int m_species1_index;
+ int m_species2_index;
+ bool m_isSameSpecies;
+ amrex::Real m_CoulombLog;
+
+ CollisionType(
+ const std::vector<std::string>& species_names,
+ std::string const collision_name);
+
+ /** Perform all binary collisions within a tile
+ *
+ * @param lev AMR level of the tile
+ * @param mfi iterator for multifab
+ * @param species1/2 pointer to species container
+ * @param isSameSpecies true if collision is between same species
+ * @param CoulombLog user input Coulomb logrithm
+ *
+ */
+
+ static void doCoulombCollisionsWithinTile (
+ int const lev, amrex::MFIter const& mfi,
+ std::unique_ptr<WarpXParticleContainer>& species1,
+ std::unique_ptr<WarpXParticleContainer>& species2,
+ bool const isSameSpecies, amrex::Real const CoulombLog );
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp
new file mode 100644
index 000000000..b8014579d
--- /dev/null
+++ b/Source/Particles/Collision/CollisionType.cpp
@@ -0,0 +1,241 @@
+#include "CollisionType.H"
+#include "ShuffleFisherYates.H"
+#include "ElasticCollisionPerez.H"
+#include <WarpX.H>
+
+CollisionType::CollisionType(
+ const std::vector<std::string>& species_names,
+ std::string const collision_name)
+{
+
+#if defined WARPX_DIM_XZ
+ amrex::Abort("Collisions only work in 3D geometry for now.");
+#elif defined WARPX_DIM_RZ
+ amrex::Abort("Collisions only work in Cartesian geometry for now.");
+#endif
+
+ // read collision species
+ std::vector<std::string> collision_species;
+ amrex::ParmParse pp(collision_name);
+ pp.getarr("species", collision_species);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(collision_species.size() == 2,
+ "Collision species must name exactly two species.");
+
+ // default Coulomb log, if < 0, will be computed automatically
+ m_CoulombLog = -1.0;
+ pp.query("CoulombLog", m_CoulombLog);
+
+ for (int i=0; i<species_names.size(); i++)
+ {
+ if (species_names[i] == collision_species[0])
+ { m_species1_index = i; }
+ if (species_names[i] == collision_species[1])
+ { m_species2_index = i; }
+ }
+
+ if (collision_species[0] == collision_species[1])
+ m_isSameSpecies = true;
+ else
+ m_isSameSpecies = false;
+
+}
+
+using namespace amrex;
+// Define shortcuts for frequently-used type names
+using ParticleType = WarpXParticleContainer::ParticleType;
+using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+using ParticleBins = DenseBins<ParticleType>;
+using index_type = ParticleBins::index_type;
+
+namespace {
+
+ /* Find the particles and count the particles that are in each cell.
+ Note that this does *not* rearrange particle arrays */
+ ParticleBins
+ findParticlesInEachCell( int const lev, MFIter const& mfi,
+ ParticleTileType const& ptile) {
+
+ // Extract particle structures for this tile
+ int const np = ptile.numParticles();
+ ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data();
+
+ // Extract box properties
+ Geometry const& geom = WarpX::GetInstance().Geom(lev);
+ Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
+ const auto lo = lbound(cbx);
+ const auto dxi = geom.InvCellSizeArray();
+ const auto plo = geom.ProbLoArray();
+
+ // Find particles that are in each cell ;
+ // results are stored in the object `bins`.
+ ParticleBins bins;
+ bins.build(np, particle_ptr, cbx,
+ // Pass lambda function that returns the cell index
+ [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect
+ {
+ return IntVect(AMREX_D_DECL((p.pos(0)-plo[0])*dxi[0] - lo.x,
+ (p.pos(1)-plo[1])*dxi[1] - lo.y,
+ (p.pos(2)-plo[2])*dxi[2] - lo.z));
+ });
+
+ return bins;
+ }
+
+}
+
+/** Perform all binary collisions within a tile
+ *
+ * @param lev AMR level of the tile
+ * @param mfi iterator for multifab
+ * @param species1/2 pointer to species container
+ * @param isSameSpecies true if collision is between same species
+ * @param CoulombLog user input Coulomb logrithm
+ *
+ */
+void CollisionType::doCoulombCollisionsWithinTile
+ ( int const lev, MFIter const& mfi,
+ std::unique_ptr<WarpXParticleContainer>& species_1,
+ std::unique_ptr<WarpXParticleContainer>& species_2,
+ bool const isSameSpecies, Real const CoulombLog )
+{
+
+ if ( isSameSpecies ) // species_1 == species_2
+ {
+ // Extract particles in the tile that `mfi` points to
+ ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi);
+
+ // Find the particles that are in each cell of this tile
+ ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
+
+ // Loop over cells, and collide the particles in each cell
+
+ // Extract low-level data
+ int const n_cells = bins_1.numBins();
+ // - Species 1
+ auto& soa_1 = ptile_1.GetStructOfArrays();
+ ParticleReal * const AMREX_RESTRICT ux_1 =
+ soa_1.GetRealData(PIdx::ux).data();
+ ParticleReal * const AMREX_RESTRICT uy_1 =
+ soa_1.GetRealData(PIdx::uy).data();
+ ParticleReal * const AMREX_RESTRICT uz_1 =
+ soa_1.GetRealData(PIdx::uz).data();
+ ParticleReal const * const AMREX_RESTRICT w_1 =
+ soa_1.GetRealData(PIdx::w).data();
+ index_type* indices_1 = bins_1.permutationPtr();
+ index_type const* cell_offsets_1 = bins_1.offsetsPtr();
+ Real q1 = species_1->getCharge();
+ Real m1 = species_1->getMass();
+
+ const Real dt = WarpX::GetInstance().getdt(lev);
+ Geometry const& geom = WarpX::GetInstance().Geom(lev);
+ const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2);
+
+ // Loop over cells
+ amrex::ParallelFor( n_cells,
+ [=] AMREX_GPU_DEVICE (int i_cell) noexcept
+ {
+ // The particles from species1 that are in the cell `i_cell` are
+ // given by the `indices_1[cell_start_1:cell_stop_1]`
+ index_type const cell_start_1 = cell_offsets_1[i_cell];
+ index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
+ index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
+
+ // Do not collide if there is only one particle in the cell
+ if ( cell_stop_1 - cell_start_1 >= 2 )
+ {
+ // shuffle
+ ShuffleFisherYates(
+ indices_1, cell_start_1, cell_half_1 );
+
+ // Call the function in order to perform collisions
+ ElasticCollisionPerez(
+ cell_start_1, cell_half_1,
+ cell_half_1, cell_stop_1,
+ indices_1, indices_1,
+ ux_1, uy_1, uz_1, ux_1, uy_1, uz_1, w_1, w_1,
+ q1, q1, m1, m1, Real(-1.0), Real(-1.0),
+ dt, CoulombLog, dV );
+ }
+ }
+ );
+ }
+ else // species_1 != species_2
+ {
+ // Extract particles in the tile that `mfi` points to
+ ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi);
+ ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi);
+
+ // Find the particles that are in each cell of this tile
+ ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
+ ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 );
+
+ // Loop over cells, and collide the particles in each cell
+
+ // Extract low-level data
+ int const n_cells = bins_1.numBins();
+ // - Species 1
+ auto& soa_1 = ptile_1.GetStructOfArrays();
+ ParticleReal * const AMREX_RESTRICT ux_1 =
+ soa_1.GetRealData(PIdx::ux).data();
+ ParticleReal * const AMREX_RESTRICT uy_1 =
+ soa_1.GetRealData(PIdx::uy).data();
+ ParticleReal * const AMREX_RESTRICT uz_1 =
+ soa_1.GetRealData(PIdx::uz).data();
+ ParticleReal const * const AMREX_RESTRICT w_1 =
+ soa_1.GetRealData(PIdx::w).data();
+ index_type* indices_1 = bins_1.permutationPtr();
+ index_type const* cell_offsets_1 = bins_1.offsetsPtr();
+ Real q1 = species_1->getCharge();
+ Real m1 = species_1->getMass();
+ // - Species 2
+ auto& soa_2 = ptile_2.GetStructOfArrays();
+ Real* ux_2 = soa_2.GetRealData(PIdx::ux).data();
+ Real* uy_2 = soa_2.GetRealData(PIdx::uy).data();
+ Real* uz_2 = soa_2.GetRealData(PIdx::uz).data();
+ Real* w_2 = soa_2.GetRealData(PIdx::w).data();
+ index_type* indices_2 = bins_2.permutationPtr();
+ index_type const* cell_offsets_2 = bins_2.offsetsPtr();
+ Real q2 = species_2->getCharge();
+ Real m2 = species_2->getMass();
+
+ const Real dt = WarpX::GetInstance().getdt(lev);
+ Geometry const& geom = WarpX::GetInstance().Geom(lev);
+ const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2);
+
+ // Loop over cells
+ amrex::ParallelFor( n_cells,
+ [=] AMREX_GPU_DEVICE (int i_cell) noexcept
+ {
+ // The particles from species1 that are in the cell `i_cell` are
+ // given by the `indices_1[cell_start_1:cell_stop_1]`
+ index_type const cell_start_1 = cell_offsets_1[i_cell];
+ index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
+ // Same for species 2
+ index_type const cell_start_2 = cell_offsets_2[i_cell];
+ index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
+
+ // ux from species1 can be accessed like this:
+ // ux_1[ indices_1[i] ], where i is between
+ // cell_start_1 (inclusive) and cell_start_2 (exclusive)
+
+ // 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 )
+ {
+ // shuffle
+ ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1);
+ ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2);
+
+ // Call the function in order to perform collisions
+ ElasticCollisionPerez(
+ cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
+ indices_1, indices_2,
+ ux_1, uy_1, uz_1, ux_2, uy_2, uz_2, w_1, w_2,
+ q1, q2, m1, m2, Real(-1.0), Real(-1.0),
+ dt, CoulombLog, dV );
+ }
+ }
+ );
+ } // end if ( isSameSpecies)
+
+}
diff --git a/Source/Particles/Collision/ComputeTemperature.H b/Source/Particles/Collision/ComputeTemperature.H
new file mode 100644
index 000000000..3cc96fb52
--- /dev/null
+++ b/Source/Particles/Collision/ComputeTemperature.H
@@ -0,0 +1,39 @@
+#ifndef WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+#define WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+
+#include <WarpXConst.H>
+
+template <typename T_index, typename T_R>
+T_R ComputeTemperature (
+ T_index const Is, T_index const Ie, T_index const *I,
+ T_R const *ux, T_R const *uy, T_R const *uz, T_R const m )
+{
+
+ T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
+
+ int N = Ie - Is;
+ if ( N == 0 ) { return T_R(0.0); }
+
+ 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);
+
+ for (int i = Is; i < 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;
+ }
+
+ vx = vx / N; vy = vy / N;
+ vz = vz / N; vs = vs / N;
+
+ return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
+}
+
+#endif // WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
diff --git a/Source/Particles/Collision/ElasticCollisionPerez.H b/Source/Particles/Collision/ElasticCollisionPerez.H
new file mode 100644
index 000000000..8e16d95cc
--- /dev/null
+++ b/Source/Particles/Collision/ElasticCollisionPerez.H
@@ -0,0 +1,105 @@
+#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
+#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
+
+#include "UpdateMomentumPerezElastic.H"
+#include "ComputeTemperature.H"
+#include <WarpXConst.H>
+#include <AMReX_Random.H>
+
+/** \brief Prepare information for and call
+ * UpdateMomentumPerezElastic().
+ * @param[in] I1s,I2s is the start index for I1,I2 (inclusive).
+ * @param[in] I1e,I2e is the start index for I1,I2 (exclusive).
+ * @param[in] I1 and I2 are the index arrays.
+ * @param[in,out] u1 and u2 are the velocity arrays (u=v*gamma),
+ * they could be either different or the same,
+ * their lengths are not needed,
+ * @param[in] I1 and I2 determine all elements that will be used.
+ * @param[in] w1 and w2 are arrays of weights.
+ * @param[in] q1 and q2 are charges. m1 and m2 are masses.
+ * @param[in] T1 and T2 are temperatures (Joule)
+ * and will be used if greater than zero,
+ * otherwise will be computed.
+ * @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.
+ * @param[in] dV is the volume of the corresponding cell.
+*/
+
+template <typename T_index, typename T_R>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void ElasticCollisionPerez (
+ T_index const I1s, T_index const I1e,
+ T_index const I2s, T_index const I2e,
+ T_index *I1, T_index *I2,
+ T_R *u1x, T_R *u1y, T_R *u1z,
+ T_R *u2x, T_R *u2y, T_R *u2z,
+ T_R const *w1, T_R const *w2,
+ T_R const q1, T_R const q2,
+ T_R const m1, T_R const m2,
+ T_R const T1, T_R const T2,
+ T_R const dt, T_R const L, T_R const dV)
+{
+
+ T_R constexpr inv_c2 = T_R(1.0)/(PhysConst::c*PhysConst::c);
+ int NI1 = I1e - I1s;
+ int NI2 = I2e - I2s;
+
+ // get local T1t and T2t
+ T_R T1t; T_R T2t;
+ if ( T1 <= T_R(0.0) && L <= T_R(0.0) )
+ {
+ T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1);
+ }
+ else { T1t = T1; }
+ if ( T2 <= T_R(0.0) && L <= T_R(0.0) )
+ {
+ T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2);
+ }
+ else { T2t = T2; }
+
+ // local density
+ T_R n1 = T_R(0.0);
+ T_R n2 = T_R(0.0);
+ T_R n12 = T_R(0.0);
+ for (int i1=I1s; i1<I1e; ++i1) { n1 += w1[ I1[i1] ]; }
+ for (int i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; }
+ n1 = n1 / dV; n2 = n2 / dV;
+ {
+ int i1 = I1s; int i2 = I2s;
+ for (int k = 0; k < amrex::max(NI1,NI2); ++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
+ T_R lmdD;
+ lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) +
+ n2*q2*q2/(T2t*PhysConst::ep0) );
+ T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) *
+ amrex::max(n1,n2), T_R(-1.0/3.0) );
+ lmdD = amrex::max(lmdD, rmin);
+
+ // call UpdateMomentumPerezElastic()
+ {
+ int i1 = I1s; int i2 = I2s;
+ for (int k = 0; k < amrex::max(NI1,NI2); ++k)
+ {
+ 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);
+ ++i1; if ( i1 == I1e ) { i1 = I1s; }
+ ++i2; if ( i2 == I2e ) { i2 = I2s; }
+ }
+ }
+
+}
+
+#endif // WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package
new file mode 100644
index 000000000..163508fb7
--- /dev/null
+++ b/Source/Particles/Collision/Make.package
@@ -0,0 +1,10 @@
+CEXE_headers += CollisionType.H
+CEXE_headers += ElasticCollisionPerez.H
+CEXE_headers += ShuffleFisherYates.H
+CEXE_headers += UpdateMomentumPerezElastic.H
+CEXE_headers += ComputeTemperature.H
+
+CEXE_sources += CollisionType.cpp
+
+INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision
+VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision
diff --git a/Source/Particles/Collision/ShuffleFisherYates.H b/Source/Particles/Collision/ShuffleFisherYates.H
new file mode 100644
index 000000000..621e654d6
--- /dev/null
+++ b/Source/Particles/Collision/ShuffleFisherYates.H
@@ -0,0 +1,29 @@
+#ifndef WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
+#define WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
+
+#include <AMReX_Random.H>
+
+/* \brief Shuffle array according to Fisher-Yates algorithm.
+ * Only shuffle the part between is <= i < ie, n = ie-is.
+ * T_index shall be
+ * amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type
+*/
+
+template <typename T_index>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void ShuffleFisherYates (T_index *array, T_index const is, T_index const ie)
+{
+ int j;
+ T_index buf;
+ for (int i = ie-1; i >= is+1; --i)
+ {
+ // get random number j: is <= j <= i
+ j = amrex::Random_int(i-is+1) + is;
+ // swop the ith array element with the jth
+ buf = array[i];
+ array[i] = array[j];
+ array[j] = buf;
+ }
+}
+
+#endif // WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H
new file mode 100644
index 000000000..948e8b075
--- /dev/null
+++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H
@@ -0,0 +1,252 @@
+#ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
+#define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
+
+#include <WarpXConst.H>
+#include <AMReX_Random.H>
+#include <cmath> // isnan() isinf()
+#include <limits> // numeric_limits<float>::min()
+
+/* \brief Update particle velocities according to
+ * F. Perez et al., Phys.Plasmas.19.083104 (2012),
+ * which is based on Nanbu's method, PhysRevE.55.4642 (1997).
+ * @param[in] LmdD is max(Debye length, minimal interparticle distance).
+ * @param[in] L is the Coulomb log. A fixed L will be used if L > 0,
+ * otherwise L will be calculated based on the algorithm.
+ * To see if there are nan or inf updated velocities,
+ * compile with USE_ASSERTION=TRUE.
+*/
+
+template <typename T_R>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumPerezElastic (
+ T_R& u1x, T_R& u1y, T_R& u1z, T_R& u2x, T_R& u2y, T_R& u2z,
+ T_R const n1, T_R const n2, T_R const n12,
+ T_R const q1, T_R const m1, T_R const w1,
+ T_R const q2, T_R const m2, T_R const w2,
+ T_R const dt, T_R const L, T_R const lmdD)
+{
+
+ // If g = u1 - u2 = 0, do not collide.
+ if ( std::abs(u1x-u2x) < std::numeric_limits<T_R>::min() &&
+ std::abs(u1y-u2y) < std::numeric_limits<T_R>::min() &&
+ std::abs(u1z-u2z) < std::numeric_limits<T_R>::min() )
+ { return; }
+
+ T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
+
+ // Compute Lorentz factor gamma
+ T_R const g1 = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 );
+ T_R const g2 = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 );
+
+ // Compute momenta
+ T_R const p1x = u1x * m1;
+ T_R const p1y = u1y * m1;
+ T_R const p1z = u1z * m1;
+ T_R const p2x = u2x * m2;
+ T_R const p2y = u2y * m2;
+ T_R const p2z = u2z * m2;
+
+ // Compute center-of-mass (COM) velocity and gamma
+ T_R const mass_g = m1 * g1 + m2 * g2;
+ T_R const vcx = (p1x+p2x) / mass_g;
+ T_R const vcy = (p1y+p2y) / mass_g;
+ T_R const vcz = (p1z+p2z) / mass_g;
+ T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz;
+ T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2 );
+
+ // Compute vc dot v1 and v2
+ T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
+ T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
+
+ // Compute p1 star
+ T_R p1sx;
+ T_R p1sy;
+ T_R p1sz;
+ if ( vcms > std::numeric_limits<T_R>::min() )
+ {
+ T_R const lorentz_tansform_factor =
+ ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1;
+ p1sx = p1x + vcx*lorentz_tansform_factor;
+ p1sy = p1y + vcy*lorentz_tansform_factor;
+ p1sz = p1z + vcz*lorentz_tansform_factor;
+ }
+ else // If vcms = 0, don't do Lorentz-transform.
+ {
+ p1sx = p1x;
+ p1sy = p1y;
+ p1sz = p1z;
+ }
+ T_R const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz );
+
+ // Compute gamma star
+ T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1;
+ T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2;
+
+ // Compute the Coulomb log lnLmd
+ T_R lnLmd;
+ if ( L > T_R(0.0) ) { lnLmd = L; }
+ else
+ {
+ // Compute b0
+ T_R const b0 = std::abs(q1*q2) * inv_c2 /
+ (T_R(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g *
+ ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) );
+
+ // Compute the minimal impact parameter
+ T_R bmin = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,b0);
+
+ // Compute the Coulomb log lnLmd
+ lnLmd = amrex::max( T_R(2.0),
+ T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) );
+ }
+
+ // Compute s
+ T_R s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 /
+ ( T_R(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 *
+ m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g *
+ std::pow(m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_R(1.0), 2.0);
+
+ // Compute s'
+ T_R const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc);
+ T_R const sp = std::pow(T_R(4.0)*MathConst::pi/T_R(3.0),T_R(1.0/3.0)) *
+ n1*n2/n12 * dt * vrel * (m1+m2) /
+ amrex::max( m1*std::pow(n1,T_R(2.0/3.0)),
+ m2*std::pow(n2,T_R(2.0/3.0)) );
+
+ // Determine s
+ s = amrex::min(s,sp);
+
+ // Get random numbers
+ T_R r = amrex::Random();
+
+ // Compute scattering angle
+ T_R cosXs;
+ T_R sinXs;
+ if ( s <= T_R(0.1) )
+ {
+ while ( true )
+ {
+ cosXs = T_R(1.0) + s * std::log(r);
+ // Avoid the bug when r is too small such that cosXs < -1
+ if ( cosXs >= T_R(-1.0) ) { break; }
+ r = amrex::Random();
+ }
+ }
+ else if ( s > T_R(0.1) && s <= T_R(3.0) )
+ {
+ T_R const Ainv = 0.0056958 + 0.9560202*s - 0.508139*s*s +
+ 0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s;
+ cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) +
+ T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) );
+ }
+ else if ( s > T_R(3.0) && s <= T_R(6.0) )
+ {
+ T_R const A = T_R(3.0) * std::exp(-s);
+ cosXs = T_R(1.0)/A * std::log( std::exp(-A) +
+ T_R(2.0) * r * std::sinh(A) );
+ }
+ else
+ {
+ cosXs = T_R(2.0) * r - T_R(1.0);
+ }
+ sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs);
+
+ // Get random azimuthal angle
+ T_R const phis = amrex::Random() * T_R(2.0) * MathConst::pi;
+ T_R const cosphis = std::cos(phis);
+ T_R const sinphis = std::sin(phis);
+
+ // Compute post-collision momenta pfs in COM
+ T_R p1fsx;
+ T_R p1fsy;
+ T_R p1fsz;
+ // p1sp is the p1s perpendicular
+ T_R p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy );
+ // Make sure p1sp is not almost zero
+ if ( p1sp > std::numeric_limits<T_R>::min() )
+ {
+ p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
+ ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
+ ( p1sx ) * cosXs;
+ p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
+ (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
+ ( p1sy ) * cosXs;
+ p1fsz = (-p1sp ) * sinXs*cosphis +
+ ( T_R(0.0) ) * sinXs*sinphis +
+ ( p1sz ) * cosXs;
+ // Note a negative sign is different from
+ // Eq. (12) in Perez's paper,
+ // but they are the same due to the random nature of phis.
+ }
+ else
+ {
+ // If the previous p1sp is almost zero
+ // x->y y->z z->x
+ // This set is equivalent to the one in Nanbu's paper
+ p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz );
+ p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
+ ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
+ ( p1sy ) * cosXs;
+ p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
+ (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
+ ( p1sz ) * cosXs;
+ p1fsx = (-p1sp ) * sinXs*cosphis +
+ ( T_R(0.0) ) * sinXs*sinphis +
+ ( p1sx ) * cosXs;
+ }
+
+ T_R const p2fsx = -p1fsx;
+ T_R const p2fsy = -p1fsy;
+ T_R const p2fsz = -p1fsz;
+
+ // Transform from COM to lab frame
+ T_R p1fx; T_R p2fx;
+ T_R p1fy; T_R p2fy;
+ T_R p1fz; T_R p2fz;
+ if ( vcms > std::numeric_limits<T_R>::min() )
+ {
+ T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz;
+ T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz;
+ T_R const factor = (gc-T_R(1.0))/vcms;
+ T_R const factor1 = factor*vcDp1fs + m1*g1s*gc;
+ T_R const factor2 = factor*vcDp2fs + m2*g2s*gc;
+ p1fx = p1fsx + vcx * factor1;
+ p1fy = p1fsy + vcy * factor1;
+ p1fz = p1fsz + vcz * factor1;
+ p2fx = p2fsx + vcx * factor2;
+ p2fy = p2fsy + vcy * factor2;
+ p2fz = p2fsz + vcz * factor2;
+ }
+ else // If vcms = 0, don't do Lorentz-transform.
+ {
+ p1fx = p1fsx;
+ p1fy = p1fsy;
+ p1fz = p1fsz;
+ p2fx = p2fsx;
+ p2fy = p2fsy;
+ p2fz = p2fsz;
+ }
+
+ // Rejection method
+ r = amrex::Random();
+ if ( w2 > r*amrex::max(w1, w2) )
+ {
+ u1x = p1fx / m1;
+ u1y = p1fy / m1;
+ u1z = p1fz / m1;
+ AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z));
+ AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
+ }
+ r = amrex::Random();
+ if ( w1 > r*amrex::max(w1, w2) )
+ {
+ u2x = p2fx / m2;
+ u2y = p2fy / m2;
+ u2z = p2fz / m2;
+ AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z));
+ AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
+ }
+
+}
+
+#endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-13 16:06:40 +0000