From 1696315b79cd141560fcdcf4fe08e70183b78023 Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Fri, 12 Jul 2024 11:11:51 -0700 Subject: [PATCH] Refactored binary collision method. Density and temperature are only 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. --- .../BinaryCollision/BinaryCollision.H | 115 ++++++++++++++++++ .../Coulomb/ElasticCollisionPerez.H | 90 +++++--------- .../Coulomb/PairWiseCoulombCollisionFunc.H | 11 +- .../Coulomb/UpdateMomentumPerezElastic.H | 98 ++++++++------- .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 4 + .../NuclearFusion/NuclearFusionFunc.H | 4 + 6 files changed, 214 insertions(+), 108 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index 84e0fef8c9c..5960eba3562 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -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" @@ -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 n1_vec; + amrex::Gpu::DeviceVector 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 @@ -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 n1_vec, n2_vec; + amrex::Gpu::DeviceVector 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, @@ -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( 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]; @@ -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(min_N+max_N-1)/dV; } + else { n12 = wpmax*static_cast(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]]; diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index 26782c2e42a..3322c19ed2d 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -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; } @@ -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. @@ -87,6 +90,8 @@ public: index_type const* AMREX_RESTRICT I2, const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*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, @@ -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; }; diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H index 3101047e211..c48153b6b92 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H @@ -18,9 +18,12 @@ /* \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] bmax 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. + * @param[in] n12 = max(w1,w2)*min(N1,N2)/dV is the effective density used for s12 + * @param[in] sigma_max is the maximum cross section based on mfp = atomic spacing + * used for the normalized scattering length s12 * To see if there are nan or inf updated velocities, * compile with USE_ASSERTION=TRUE. * @@ -33,11 +36,11 @@ template AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumPerezElastic ( T_PR& u1x, T_PR& u1y, T_PR& u1z, T_PR& u2x, T_PR& u2y, T_PR& u2z, - T_PR const n1, T_PR const n2, T_PR const n12, T_PR const q1, T_PR const m1, T_PR const w1, T_PR const q2, T_PR const m2, T_PR const w2, - T_R const dt, T_PR const L, T_PR const lmdD, - amrex::RandomEngine const& engine) + T_PR const n12, T_PR const sigma_max, + T_PR const L, T_PR const bmax, + T_R const dt, amrex::RandomEngine const& engine ) { T_PR const diffx = amrex::Math::abs(u1x-u2x); @@ -102,58 +105,59 @@ void UpdateMomentumPerezElastic ( T_PR const g1s = ( T_PR(1.0) - vcDv1*inv_c2 )*gc*g1; T_PR const g2s = ( T_PR(1.0) - vcDv2*inv_c2 )*gc*g2; - // Compute s - T_PR s = 0; + // Compute variant relative velocity in cm frame + T_PR const muRst = g1s*m1*g2s*m2/(g1s*m1 + g2s*m2); + T_PR const vrelst = p1sm/muRst; + + // Compute invariant relative velocity (frame independent) + T_PR const denom = T_PR(1.0) + p1sm*p1sm/(m1*g1s*m2*g2s)*inv_c2; // (1.0 - v1s*v2s/c^2) + T_PR const vrelst_invar = vrelst/denom; + + // Compute s12 + T_PR s12 = 0; if (p1sm > std::numeric_limits::min()) { - // s is non-zero (i.e. particles scatter) only if the relative - // motion between particles is not negligible (p1sm non-zero) + // Writing b0 in a form that is directly analagous to the well-known non-relativistic form. + // See Eq. 3.3.2 in Principles of Plasma Discharges and Material Processing by + // M. A. Lieberman and A. J. Lichtenberg. + // Note that b0 on Eq. 22 of Perez POP 19 (2012) is bmin = b0/2, + // Note: there is a typo in Eq 22 of Perez, the last square is incorrect! + // See the SMILEI documentation: https://smileipic.github.io/Smilei/Understand/collisions.html + // and https://github.com/ECP-WarpX/WarpX/files/3799803/main.pdf from GitHub #429 + T_PR const b0 = amrex::Math::abs(q1*q2) / + (T_PR(2.0)*MathConst::pi*PhysConst::ep0*muRst*vrelst*vrelst_invar); + // Compute the Coulomb log lnLmd first T_PR lnLmd; if ( L > T_PR(0.0) ) { lnLmd = L; } else { - // Compute b0 according to eq (22) from Perez et al., Phys.Plasmas.19.083104 (2012) - // Note: there is a typo in the equation, the last square is incorrect! - // See the SMILEI documentation: https://smileipic.github.io/Smilei/Understand/collisions.html - // and https://github.com/ECP-WarpX/WarpX/files/3799803/main.pdf from GitHub #429 - T_PR const b0 = amrex::Math::abs(q1*q2) * inv_c2 / - (T_PR(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g * - ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_PR(1.0) ); - - // Compute the minimal impact parameter - constexpr T_PR hbar_pi = static_cast(PhysConst::hbar*MathConst::pi); - const T_PR bmin = amrex::max(hbar_pi/p1sm, b0); + + // Compute the minimum impact parameter from quantum: bqm = lDB/(4*pi) = hbar/(2*p1sm) + // See NRL formulary. Also see "An introduction to the physics of the Coulomb logarithm, + // with emphasis on quantum-mechanical effects", J. Plasma Phys. vol. 85 (2019). by J.A. Krommes. + // Note: The formula in Perez 2012 and in Lee and More 1984 uses h rather than + // hbar for bqm. If this is used, then the transition energy where bmin goes from classical + // to quantum is only 2.5 eV for electrons; compared to 100 eV when using hbar. + const T_PR bmin_qm = static_cast(PhysConst::hbar*0.5/p1sm); + + // Set the minimum impact parameter + const T_PR bmin = amrex::max(bmin_qm, T_PR(0.5)*b0); // Compute the Coulomb log lnLmd lnLmd = amrex::max( T_PR(2.0), - T_PR(0.5)*std::log(T_PR(1.0)+lmdD*lmdD/(bmin*bmin)) ); + T_PR(0.5)*std::log(T_PR(1.0) + bmax*bmax/(bmin*bmin)) ); } - // Compute s - const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_PR(1.0); - const auto tts2 = tts*tts; - s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / - ( T_PR(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * - m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2; - - // Compute s' - const auto cbrt_n1 = std::cbrt(n1); - const auto cbrt_n2 = std::cbrt(n2); - const auto coeff = static_cast( - std::pow(4.0*MathConst::pi/3.0,1.0/3.0)); - T_PR const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); - T_PR const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) / - amrex::max( m1*cbrt_n1*cbrt_n1, - m2*cbrt_n2*cbrt_n2); - - // Determine s - s = amrex::min(s,sp); + // Compute s12 with sigma limited by sigma_max where mfp = atomic spacing) + const T_PR sigma_eff = amrex::min(MathConst::pi*b0*b0*lnLmd,sigma_max); + s12 = sigma_eff * n12 * dt * vrelst * g1s*g2s/(g1*g2); + } // Only modify momenta if is s is non-zero - if (s > std::numeric_limits::min()) { + if (s12 > std::numeric_limits::min()) { // Get random numbers T_PR r = amrex::Random(engine); @@ -161,27 +165,27 @@ void UpdateMomentumPerezElastic ( // Compute scattering angle T_PR cosXs; T_PR sinXs; - if ( s <= T_PR(0.1) ) + if ( s12 <= T_PR(0.1) ) { while ( true ) { - cosXs = T_PR(1.0) + s * std::log(r); + cosXs = T_PR(1.0) + s12 * std::log(r); // Avoid the bug when r is too small such that cosXs < -1 if ( cosXs >= T_PR(-1.0) ) { break; } r = amrex::Random(engine); } } - else if ( s > T_PR(0.1) && s <= T_PR(3.0) ) + else if ( s12 > T_PR(0.1) && s12 <= T_PR(3.0) ) { T_PR const Ainv = static_cast( - 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); + 0.0056958 + 0.9560202*s12 - 0.508139*s12*s12 + + 0.47913906*s12*s12*s12 - 0.12788975*s12*s12*s12*s12 + 0.02389567*s12*s12*s12*s12*s12); cosXs = Ainv * std::log( std::exp(T_PR(-1.0)/Ainv) + T_PR(2.0) * r * std::sinh(T_PR(1.0)/Ainv) ); } - else if ( s > T_PR(3.0) && s <= T_PR(6.0) ) + else if ( s12 > T_PR(3.0) && s12 <= T_PR(6.0) ) { - T_PR const A = T_PR(3.0) * std::exp(-s); + T_PR const A = T_PR(3.0) * std::exp(-s12); cosXs = T_PR(1.0)/A * std::log( std::exp(-A) + T_PR(2.0) * r * std::sinh(A) ); } diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index f28f3fb984c..30b466e2ec2 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -96,6 +96,8 @@ public: index_type const* AMREX_RESTRICT I2, const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*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, @@ -189,6 +191,8 @@ public: } int m_process_count; + bool m_computeSpeciesDensities = false; + bool m_computeSpeciesTemperatures = false; ScatteringProcess::Executor* m_scattering_processes_data; }; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index b6d6fe171de..3113dc69839 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -131,6 +131,8 @@ public: index_type const* AMREX_RESTRICT I2, const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*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, @@ -232,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; };