From 1c2ee5542fbcbffc24fb0c2df3ba2102cd931870 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Fri, 8 Nov 2019 15:12:14 -0800 Subject: Add CoulombCollisions files; Compile successfully. --- Source/Particles/MultiParticleContainer.cpp | 32 +++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index fca22daa2..6b3a34e38 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -617,6 +617,38 @@ MultiParticleContainer::doFieldIonization () } // pc_source } +void +MultiParticleContainer::doCoulombCollisions () +{ + BL_PROFILE("MPC::doCoulombCollisions"); + + // At this point, the code always collides the first and second species + // TODO: Read the user input to read the different types of collisions requested + // and loop over all types of collisions, selecting each time the + // two types of species that will be collided + auto& species1 = allcontainers[0]; + auto& species2 = allcontainers[1]; + + // Enable tiling + MFItInfo info; + if (Gpu::notInLaunchRegion()) info.EnableTiling(species1->tile_size); + + // Loop over refinement levels + for (int lev = 0; lev <= species1->finestLevel(); ++lev){ + + // Loop over all grids/tiles at this level +#ifdef _OPENMP + info.SetDynamic(true); + #pragma omp parallel +#endif + for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ + + doCoulombCollisionsWithinTile( lev, mfi, species1, species2 ); + + } + } +} + #ifdef WARPX_QED void MultiParticleContainer::InitQED () { -- cgit v1.2.3 From bd8fb800bcb8675429464292cdd96dae8ac9c761 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Tue, 12 Nov 2019 16:55:56 -0800 Subject: add CollisionType --- Source/Particles/Collision/CollisionType.H | 36 +++++ Source/Particles/Collision/CollisionType.cpp | 180 +++++++++++++++++++++ Source/Particles/Collision/CoulombCollisions.cpp | 2 +- Source/Particles/Collision/Make.package | 6 +- .../Collision/UpdateMomentumPerezElastic.H | 179 ++++++++++++++++++++ Source/Particles/MultiParticleContainer.H | 7 + Source/Particles/MultiParticleContainer.cpp | 60 ++++--- 7 files changed, 446 insertions(+), 24 deletions(-) create mode 100644 Source/Particles/Collision/CollisionType.H create mode 100644 Source/Particles/Collision/CollisionType.cpp create mode 100644 Source/Particles/Collision/UpdateMomentumPerezElastic.H (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H new file mode 100644 index 000000000..30058017f --- /dev/null +++ b/Source/Particles/Collision/CollisionType.H @@ -0,0 +1,36 @@ +#ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ +#define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ + +class CollisionType +{ +public: + int m_species1; + int m_species2; + + void CollisionType( + const std::vector& species_names, + const std::vector collision_name); + + void ElasticCollisionPerez( + amrex::DenseBins::index_type const I1s, + amrex::DenseBins::index_type const I1e, + amrex::DenseBins::index_type const I2s, + amrex::DenseBins::index_type const I2e, + amrex::DenseBins::index_type *I1, + amrex::DenseBins::index_type *I2, + amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, + amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, + const amrex::Real* w1, const amrex::Real* w2, + const amrex::Real q1, const amrex::Real q2, + const amrex::Real m1, const amrex::Real m2, + const amrex::Real T1, const amrex::Real T2, + const amrex::Real dt, const amrex::Real L, const amrex::Real V); + + void ShuffleFisherYates( + amrex::DenseBins::index_type *array, + amrex::DenseBins::index_type const is, + amrex::DenseBins::index_type const ie); + +}; + +#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..736223ae2 --- /dev/null +++ b/Source/Particles/Collision/CollisionType.cpp @@ -0,0 +1,180 @@ +#include "CollisionType.H" +#include "UpdateMomentumPerezElastic.H" + +#include "WarpXParticleContainer.H" +#include +#include +#include + +#include + +void CollisionType::CollisionType( + const std::vector& species_names, + const std::vector collision_name) +{ + + std::vector collision_species; + + ParmParse pp(collision_name); + pp.getarr("species", collision_species); + + for (int i=0; i::index_type const I1s, + amrex::DenseBins::index_type const I1e, + amrex::DenseBins::index_type const I2s, + amrex::DenseBins::index_type const I2e, + amrex::DenseBins::index_type *I1, + amrex::DenseBins::index_type *I2, + amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, + amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, + const amrex::Real* w1, const amrex::Real* w2, + const amrex::Real q1, const amrex::Real q2, + const amrex::Real m1, const amrex::Real m2, + const amrex::Real T1, const amrex::Real T2, + const amrex::Real dt, const amrex::Real L, const amrex::Real V) +{ + + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + int NI1 = I1e - I1s; + int NI2 = I2e - I2s; + + // shuffle I1 and I2 + ShuffleFisherYates(I1, I1s, I1e); + ShuffleFisherYates(I2, I2s, I2e); + + // get local T1t and T2t + amrex::Real T1t; amrex::Real T2t; + if ( T1 <= 0. ) + { + amrex::Real vx = 0.; amrex::Real vy = 0.; + amrex::Real vz = 0.; amrex::Real vs = 0.; + amrex::Real gm = 0.; amrex::Real us = 0.; + for (int i1 = I1s; i1 < I1e; ++i1) + { + us = ( u1x[ I1[i1] ] * u1x[ I1[i1] ] + + u1y[ I1[i1] ] * u1y[ I1[i1] ] + + u1z[ I1[i1] ] * u1z[ I1[i1] ] ); + gm = std::sqrt(1.+us*inv_c2); + vx += u1x[ I1[i1] ] / gm; + vy += u1y[ I1[i1] ] / gm; + vz += u1z[ I1[i1] ] / gm; + vs += us / gm / gm; + } + vx = vx / NI1; vy = vy / NI1; + vz = vz / NI1; vs = vs / NI1; + T1t = m1/3.*(vs-(vx*vx+vy*vy+vz*vz)); + } + else { T1t = T1; } + if ( T2 <= 0. ) + { + amrex::Real vx = 0.; amrex::Real vy = 0.; + amrex::Real vz = 0.; amrex::Real vs = 0.; + amrex::Real gm = 0.; amrex::Real us = 0.; + for (int i2 = I2s; i2 < I2e; ++i2) + { + us = ( u2x[ I2[i2] ] * u2x[ I2[i2] ] + + u2y[ I2[i2] ] * u2y[ I2[i2] ] + + u2z[ I2[i2] ] * u2z[ I2[i2] ] ); + gm = std::sqrt(1.+us*inv_c2); + vx += u2x[ I2[i2] ] / gm; + vy += u2y[ I2[i2] ] / gm; + vz += u2z[ I2[i2] ] / gm; + vs += us / gm / gm; + } + vx = vx / NI2; vy = vy / NI2; + vz = vz / NI2; vs = vs / NI2; + T2t = m2/3.*(vs-(vx*vx+vy*vy+vz*vz)); + } + else { T2t = T2; } + + // local density + amrex::Real n1 = 0.; + amrex::Real n2 = 0.; + amrex::Real n12 = 0.; + for (int i1=I1s; i1::index_type *array, + amrex::DenseBins::index_type const is, + amrex::DenseBins::index_type const ie) +{ + int j; int buf; + for (int i = ie-1; i >= is+1; --i) + { + // get random number is <= j <= i + while ( true ) { + j = (int) floor(amrex::Random()*(i-is+1)) + is; + if ( j >= is && j <= i ) { break; } + } + buf = array[i]; + array[i] = array[j]; + array[j] = buf; + } +} diff --git a/Source/Particles/Collision/CoulombCollisions.cpp b/Source/Particles/Collision/CoulombCollisions.cpp index a97dd5ed1..0ca13bca2 100644 --- a/Source/Particles/Collision/CoulombCollisions.cpp +++ b/Source/Particles/Collision/CoulombCollisions.cpp @@ -1,7 +1,7 @@ #include "WarpXParticleContainer.H" #include #include -#include "ElasticCollisionPerez.H" +#include "CollisionType.H" using namespace amrex; // Define shortcuts for frequently-used type names diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package index 12b7e6046..66d8859f4 100644 --- a/Source/Particles/Collision/Make.package +++ b/Source/Particles/Collision/Make.package @@ -1,6 +1,8 @@ -CEXE_headers += ElasticCollisionPerez.H -CEXE_headers += ShuffleFisherYates.H +CEXE_headers += CollisionType.H CEXE_headers += CoulombCollisions.H +CEXE_headers += ShuffleFisherYates.H + +CEXE_sources += CollisionType.cpp CEXE_sources += CoulombCollisions.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H new file mode 100644 index 000000000..d5c8620e3 --- /dev/null +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -0,0 +1,179 @@ +#ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ +#define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ + +#include "WarpXParticleContainer.H" +#include +#include +#include +#include "ShuffleFisherYates.H" + +/* \brief Update particle velocities according to + * F. Perez et al., Phys. Plasmas 19, 083104 (2012). + * LmdD is max(Debye length, minimal interparticle distance). + * L is the Coulomb log. A fixed L will be used if L > 0, + * otherwise L will be calculated based on the algorithm.*/ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumPerezElastic( + amrex::Real u1x, amrex::Real u1y, amrex::Real u1z, + amrex::Real u2x, amrex::Real u2y, amrex::Real u2z, + const amrex::Real n1, const amrex::Real n2, const amrex::Real n12, + const amrex::Real q1, const amrex::Real m1, const amrex::Real w1, + const amrex::Real q2, const amrex::Real m2, const amrex::Real w2, + const amrex::Real dt, const amrex::Real L, const amrex::Real lmdD) +{ + // Temporary variabls + amrex::Real buf1; + amrex::Real buf2; + + constexpr amrex::Real inv_c2 = 1./(PhysConst::c * PhysConst::c); + + // Compute Lorentz factor gamma + amrex::Real g1 = std::sqrt(1.+(u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2); + amrex::Real g2 = std::sqrt(1.+(u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2); + + // Compute momenta + amrex::Real p1x = u1x * m1; + amrex::Real p1y = u1y * m1; + amrex::Real p1z = u1z * m1; + amrex::Real p2x = u2x * m2; + amrex::Real p2y = u2y * m2; + amrex::Real p2z = u2z * m2; + + // Compute center-of-mass (COM) velocity and gamma + buf1 = m1 * g1; + buf2 = m2 * g2; + const amrex::Real vcx = (p1x+p2x) / (buf1+buf2); + const amrex::Real vcy = (p1y+p2y) / (buf1+buf2); + const amrex::Real vcz = (p1z+p2z) / (buf1+buf2); + const amrex::Real vcm = std::sqrt(vcx*vcx+vcy*vcy+vcz*vcz); + const amrex::Real gc = 1./std::sqrt(1.-vcm*vcm*inv_c2); + + // Compute vc dot v1 and v2 + const amrex::Real vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z)/g1; + const amrex::Real vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z)/g2; + + // Compute p1 star + buf2 = ( (gc-1.)/vcm/vcm*vcDv1 - gc )*buf1; + const amrex::Real p1sx = p1x + buf2*vcx; + const amrex::Real p1sy = p1y + buf2*vcy; + const amrex::Real p1sz = p1z + buf2*vcz; + const amrex::Real p1sm = std::sqrt(p1sx*p1sx+p1sy*p1sy+p1sz*p1sz); + + // Compute gamma star + const amrex::Real g1s = (1.-vcDv1*inv_c2)*gc*g1; + const amrex::Real g2s = (1.-vcDv2*inv_c2)*gc*g2; + + // Compute the Coulomb log lnLmd + amrex::Real lnLmd; + if ( L > 0. ) + { lnLmd = L; } + else + { + // Compute b0 (buf1) + buf1 = std::abs(q1*q2) * inv_c2 / + (4.*MathConst::pi*PhysConst::ep0) * gc/(m1*g1+m2*g2) * + ( m1*g1s*m2*g2s/p1sm/p1sm/inv_c2 + 1. ); + + // Compute the minimal impact parameter (buf2) + buf2 = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,buf1); + + // Compute the Coulomb log lnLmd + lnLmd = amrex::max( 2., 0.5*std::log(1.+lmdD*lmdD/buf2/buf2) ); + } + + // Compute s + amrex::Real s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / + ( 4. * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * + m1*g1*m2*g2/inv_c2/inv_c2 ) * gc*p1sm/(m1*g1+m2*g2) * + std::pow(m1*g1s*m2*g2s/inv_c2/p1sm/p1sm + 1., 2); + + // Compute s' + buf1 = (m1*g1+m2*g2)*p1sm/(m1*g1s*m2*g2s*gc); // vrel + const amrex::Real sp = std::pow(4.*MathConst::pi/3.,1./3.) * + n1*n2/n12 * dt * buf1 * (m1+m2) / + amrex::max(m1*std::pow(n1,2./3.),m2*std::pow(n2,2./3.)); + + // Determine s + s = amrex::min(s,sp); + + // Get random numbers + buf1 = amrex::Random(); + + // Compute scattering angle + amrex::Real cosXs; + amrex::Real sinXs; + if ( s <= 0.1 ) + { cosXs = 1. + s * std::log(buf1); } + else if ( s > 0.1 && s <= 3. ) + { + // buf2 is A inverse + buf2 = 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 = buf2*std::log(std::exp(-1./buf2)+2.*buf1*std::sinh(1.0/buf2)); + } + else if ( s > 3. && s <= 6. ) + { + // buf2 is A + buf2 = 3. * std::exp(-s); + cosXs = 1./buf2*std::log(std::exp(-buf2)+2.*buf1*std::sinh(buf2)); + } + else + { + cosXs = 2.*buf1 - 1.; + } + sinXs = std::sqrt(1. - cosXs*cosXs); + + // Get random azimuthal angle + buf1 = amrex::Random() * 2. * MathConst::pi; + + // Compute post-collision momenta in COM + // buf2 is the p1s perpendicular + // p1 is now p1f star + buf2 = std::sqrt(p1sx*p1sx + p1sy*p1sy); + p1x = ( p1sx*p1sx/buf2) * sinXs*std::cos(buf1) + + (-p1sy*p1sm/buf2) * sinXs*std::sin(buf1) + + ( p1sx ) * cosXs; + p1y = ( p1sy*p1sz/buf2) * sinXs*std::cos(buf1) + + ( p1sx*p1sm/buf2) * sinXs*std::sin(buf1) + + ( p1sy ) * cosXs; + p1z = (-buf2 ) * sinXs*std::cos(buf1) + + ( 0. ) * sinXs*std::sin(buf1) + + ( p1sz ) * cosXs; + p2x = -p1x; + p2y = -p1y; + p2z = -p1z; + + // Transform from COM to lab frame + // buf1 is vc dot p1 (p1fs) + // p1 is then updated to be p1f + buf1 = vcx*p1x + vcy*p1y + vcz*p1z; + buf2 = (gc-1.)/vcm/vcm*buf1 + m1*g1s*gc; + p1x = p1x + vcx * buf2; + p1y = p1y + vcy * buf2; + p1z = p1z + vcz * buf2; + // buf1 is vc dot p2 (p2fs) + // p2 is then updated to be p2f + buf1 = vcx*p2x + vcy*p2y + vcz*p2z; + buf2 = (gc-1.)/vcm/vcm*buf1 + m2*g2s*gc; + p2x = p2x + vcx * buf2; + p2y = p2y + vcy * buf2; + p2z = p2z + vcz * buf2; + + // Rejection method + buf1 = amrex::Random(); + if ( w2/amrex::max(w1,w2) > buf1 ) + { + u1x = p1x / m1; + u1y = p1y / m1; + u1z = p1z / m1; + } + buf2 = amrex::Random(); + if ( w1/amrex::max(w1,w2) > buf2 ) + { + u2x = p2x / m2; + u2y = p2y / m2; + u2z = p2z / m2; + } +} + +#endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 9f965a88f..edace2d4f 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -3,6 +3,7 @@ #include "ElementaryProcess.H" #include "CoulombCollisions.H" +#include "CollisionType.H" #include #include @@ -138,6 +139,7 @@ public: std::unique_ptr GetChargeDensity(int lev, bool local = false); void doFieldIonization (); + void doCoulombCollisions (); void Checkpoint (const std::string& dir) const; @@ -212,6 +214,10 @@ protected: std::vector lasers_names; + std::vector collision_names; + + std::vector allcollisions; + //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid std::vector m_deposit_on_main_grid; @@ -298,5 +304,6 @@ private: // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). + int ncollisions = 0; }; #endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6b3a34e38..c1c774cdc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -54,6 +54,12 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) } } ionization_process = IonizationProcess(); + + // collision + for (int i = 0; i < ncollisions; ++i) { + allcollisions[i].reset(new CollisionType(species_names, collision_names[i])); + } + } void @@ -120,6 +126,13 @@ MultiParticleContainer::ReadParameters () } } + // collisions + ParmParse pc("collisions"); + pc.query("ncollisions", ncollisions); + BL_ASSERT(ncollisions >= 0); + pc.getarr("collision_names", collision_names); + BL_ASSERT(collision_names.size() == ncollisions); + } pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr); @@ -622,29 +635,34 @@ MultiParticleContainer::doCoulombCollisions () { BL_PROFILE("MPC::doCoulombCollisions"); - // At this point, the code always collides the first and second species - // TODO: Read the user input to read the different types of collisions requested - // and loop over all types of collisions, selecting each time the - // two types of species that will be collided - auto& species1 = allcontainers[0]; - auto& species2 = allcontainers[1]; - - // Enable tiling - MFItInfo info; - if (Gpu::notInLaunchRegion()) info.EnableTiling(species1->tile_size); - - // Loop over refinement levels - for (int lev = 0; lev <= species1->finestLevel(); ++lev){ - - // Loop over all grids/tiles at this level + for (int i = 0; i < ncollisions; ++i) + { + // At this point, the code always collides the first and second species + // TODO: Read the user input to read the different types of collisions requested + // and loop over all types of collisions, selecting each time the + // two types of species that will be collided + // auto& species1 = allcontainers[0]; + // auto& species2 = allcontainers[1]; + auto& species1 = allcollisions[i].m_species1; + auto& species2 = allcollisions[i].m_species2; + + // Enable tiling + MFItInfo info; + if (Gpu::notInLaunchRegion()) info.EnableTiling(species1->tile_size); + + // Loop over refinement levels + for (int lev = 0; lev <= species1->finestLevel(); ++lev){ + + // Loop over all grids/tiles at this level #ifdef _OPENMP - info.SetDynamic(true); - #pragma omp parallel + info.SetDynamic(true); + #pragma omp parallel #endif - for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ - - doCoulombCollisionsWithinTile( lev, mfi, species1, species2 ); - + for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ + + doCoulombCollisionsWithinTile( lev, mfi, species1, species2 ); + + } } } } -- cgit v1.2.3 From 6e370ba4658d0be08f15de55f461f3c05ca0e3a9 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Tue, 12 Nov 2019 22:05:42 -0800 Subject: This version has compiling errors. --- Source/Particles/Collision/CollisionType.H | 7 ++++++- Source/Particles/Collision/CollisionType.cpp | 15 +++++---------- Source/Particles/Collision/CoulombCollisions.cpp | 2 +- Source/Particles/Collision/Make.package | 2 +- Source/Particles/Collision/UpdateMomentumPerezElastic.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 2 +- 6 files changed, 15 insertions(+), 15 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H index 30058017f..e8f1807a2 100644 --- a/Source/Particles/Collision/CollisionType.H +++ b/Source/Particles/Collision/CollisionType.H @@ -1,13 +1,18 @@ #ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ #define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ +#include "WarpXParticleContainer.H" +#include +#include +#include + class CollisionType { public: int m_species1; int m_species2; - void CollisionType( + CollisionType( const std::vector& species_names, const std::vector collision_name); diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index 736223ae2..3b89f7810 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -1,14 +1,7 @@ #include "CollisionType.H" #include "UpdateMomentumPerezElastic.H" -#include "WarpXParticleContainer.H" -#include -#include -#include - -#include - -void CollisionType::CollisionType( +CollisionType::CollisionType( const std::vector& species_names, const std::vector collision_name) { @@ -45,6 +38,7 @@ void CollisionType::CollisionType( * L is the Coulomb log and will be used if greater than zero, * otherwise will be computed. * V is the volume of the corresponding cell.*/ + void CollisionType::ElasticCollisionPerez( amrex::DenseBins::index_type const I1s, amrex::DenseBins::index_type const I1e, @@ -66,8 +60,8 @@ void CollisionType::ElasticCollisionPerez( int NI2 = I2e - I2s; // shuffle I1 and I2 - ShuffleFisherYates(I1, I1s, I1e); - ShuffleFisherYates(I2, I2s, I2e); + CollisionType::ShuffleFisherYates(I1, I1s, I1e); + CollisionType::ShuffleFisherYates(I2, I2s, I2e); // get local T1t and T2t amrex::Real T1t; amrex::Real T2t; @@ -160,6 +154,7 @@ void CollisionType::ElasticCollisionPerez( /* \brief Shuffle array according to Fisher-Yates algorithm. * Only shuffle the part between is <= i < ie, n = ie-is.*/ + void CollisionType::ShuffleFisherYates( amrex::DenseBins::index_type *array, amrex::DenseBins::index_type const is, diff --git a/Source/Particles/Collision/CoulombCollisions.cpp b/Source/Particles/Collision/CoulombCollisions.cpp index 0ca13bca2..3104d71fe 100644 --- a/Source/Particles/Collision/CoulombCollisions.cpp +++ b/Source/Particles/Collision/CoulombCollisions.cpp @@ -111,7 +111,7 @@ doCoulombCollisionsWithinTile ( int lev, MFIter const& mfi, // Call the function in order to perform collisions - ElasticCollisionPerez( + CollisionType::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, diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package index 66d8859f4..3134a29b8 100644 --- a/Source/Particles/Collision/Make.package +++ b/Source/Particles/Collision/Make.package @@ -1,6 +1,6 @@ CEXE_headers += CollisionType.H CEXE_headers += CoulombCollisions.H -CEXE_headers += ShuffleFisherYates.H +CEXE_headers += UpdateMomentumPerezElastic.H CEXE_sources += CollisionType.cpp CEXE_sources += CoulombCollisions.cpp diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H index d5c8620e3..156f03b19 100644 --- a/Source/Particles/Collision/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -5,13 +5,13 @@ #include #include #include -#include "ShuffleFisherYates.H" /* \brief Update particle velocities according to * F. Perez et al., Phys. Plasmas 19, 083104 (2012). * LmdD is max(Debye length, minimal interparticle distance). * L is the Coulomb log. A fixed L will be used if L > 0, * otherwise L will be calculated based on the algorithm.*/ + AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumPerezElastic( amrex::Real u1x, amrex::Real u1y, amrex::Real u1z, diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index c1c774cdc..823756447 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -126,7 +126,7 @@ MultiParticleContainer::ReadParameters () } } - // collisions + // collision ParmParse pc("collisions"); pc.query("ncollisions", ncollisions); BL_ASSERT(ncollisions >= 0); -- cgit v1.2.3 From 92484d9eade085bbe685bf33f3956ce9ec650ea8 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Wed, 13 Nov 2019 08:02:52 -0800 Subject: Use allcontainer. --- Source/Particles/MultiParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 823756447..bd67e198a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -643,8 +643,8 @@ MultiParticleContainer::doCoulombCollisions () // two types of species that will be collided // auto& species1 = allcontainers[0]; // auto& species2 = allcontainers[1]; - auto& species1 = allcollisions[i].m_species1; - auto& species2 = allcollisions[i].m_species2; + auto& species1 = allcontainers[ allcollisions[i].m_species1 ]; + auto& species2 = allcontainers[ allcollisions[i].m_species2 ]; // Enable tiling MFItInfo info; -- cgit v1.2.3 From cfcf3fb326ca57317cfd629cf6716a7d16c8a299 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Wed, 13 Nov 2019 09:08:15 -0800 Subject: Compile successfully! --- Source/Particles/Collision/CollisionType.H | 4 ++-- Source/Particles/Collision/CollisionType.cpp | 10 ++++++---- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 ++-- 4 files changed, 11 insertions(+), 9 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H index 3671296f3..327ac0280 100644 --- a/Source/Particles/Collision/CollisionType.H +++ b/Source/Particles/Collision/CollisionType.H @@ -2,7 +2,7 @@ #define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ #include "WarpXParticleContainer.H" -#include +//#include #include #include @@ -14,7 +14,7 @@ public: CollisionType( const std::vector& species_names, - const std::vector collision_name); + std::string collision_name); static void ElasticCollisionPerez( amrex::DenseBins::index_type const I1s, diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index 9e936b4fd..80f155a80 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -1,14 +1,16 @@ #include "CollisionType.H" #include "UpdateMomentumPerezElastic.H" +#include + CollisionType::CollisionType( const std::vector& species_names, - const std::vector collision_name) + std::string collision_name) { std::vector collision_species; - ParmParse pp(collision_name); + amrex::ParmParse pp(collision_name); pp.getarr("species", collision_species); for (int i=0; i::index_type const I1s, amrex::DenseBins::index_type const I1e, amrex::DenseBins::index_type const I2s, @@ -155,7 +157,7 @@ static void CollisionType::ElasticCollisionPerez( /* \brief Shuffle array according to Fisher-Yates algorithm. * Only shuffle the part between is <= i < ie, n = ie-is.*/ -static void CollisionType::ShuffleFisherYates( +void CollisionType::ShuffleFisherYates( amrex::DenseBins::index_type *array, amrex::DenseBins::index_type const is, amrex::DenseBins::index_type const ie) diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index edace2d4f..240b94273 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -216,7 +216,7 @@ protected: std::vector collision_names; - std::vector allcollisions; + amrex::Vector > allcollisions; //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid std::vector m_deposit_on_main_grid; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index bd67e198a..5158be996 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -643,8 +643,8 @@ MultiParticleContainer::doCoulombCollisions () // two types of species that will be collided // auto& species1 = allcontainers[0]; // auto& species2 = allcontainers[1]; - auto& species1 = allcontainers[ allcollisions[i].m_species1 ]; - auto& species2 = allcontainers[ allcollisions[i].m_species2 ]; + auto& species1 = allcontainers[ allcollisions[i]->m_species1 ]; + auto& species2 = allcontainers[ allcollisions[i]->m_species2 ]; // Enable tiling MFItInfo info; -- cgit v1.2.3 From 6e59c708665b63db88ffc890d679072506956e67 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Wed, 13 Nov 2019 09:42:28 -0800 Subject: Changed query. --- Source/Particles/MultiParticleContainer.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 5158be996..907dcea15 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -130,8 +130,10 @@ MultiParticleContainer::ReadParameters () ParmParse pc("collisions"); pc.query("ncollisions", ncollisions); BL_ASSERT(ncollisions >= 0); - pc.getarr("collision_names", collision_names); - BL_ASSERT(collision_names.size() == ncollisions); + if (ncollisions > 0) { + pc.getarr("collision_names", collision_names); + BL_ASSERT(collision_names.size() == ncollisions); + } } -- cgit v1.2.3 From e182ca4464c90f6c74f9c3d6b37d7f0c3a447a26 Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Wed, 13 Nov 2019 09:43:59 -0800 Subject: Changed format. --- Source/Particles/MultiParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 907dcea15..3d8dc0ddc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -130,7 +130,7 @@ MultiParticleContainer::ReadParameters () ParmParse pc("collisions"); pc.query("ncollisions", ncollisions); BL_ASSERT(ncollisions >= 0); - if (ncollisions > 0) { + if (ncollisions > 0) { pc.getarr("collision_names", collision_names); BL_ASSERT(collision_names.size() == ncollisions); } -- cgit v1.2.3 From 7c5126df0a02dbe72f8e9049be6d16c3600c723a Mon Sep 17 00:00:00 2001 From: Yin-YinjianZhao Date: Wed, 13 Nov 2019 16:11:46 -0800 Subject: First version that can run with collisions! --- Source/Evolve/WarpXEvolveEM.cpp | 1 + Source/Particles/Collision/CollisionType.H | 26 +- Source/Particles/Collision/CollisionType.cpp | 271 +++++++++------------ Source/Particles/Collision/ElasticCollisionPerez.H | 144 +++++++++++ Source/Particles/Collision/Make.package | 4 +- Source/Particles/Collision/ShuffleFisherYates.H | 33 +++ .../Collision/UpdateMomentumPerezElastic.H | 2 +- Source/Particles/MultiParticleContainer.H | 4 +- Source/Particles/MultiParticleContainer.cpp | 4 +- 9 files changed, 309 insertions(+), 180 deletions(-) create mode 100644 Source/Particles/Collision/ElasticCollisionPerez.H create mode 100644 Source/Particles/Collision/ShuffleFisherYates.H (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index b5fd52bdc..68ea807f1 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -280,6 +280,7 @@ WarpX::OneStep_nosub (Real cur_time) // Loop over species. For each ionizable species, create particles in // product species. mypc->doFieldIonization(); + mypc->doCoulombCollisions(); // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} // Deposit current j^{n+1/2} diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H index 327ac0280..c9ba643d8 100644 --- a/Source/Particles/Collision/CollisionType.H +++ b/Source/Particles/Collision/CollisionType.H @@ -2,9 +2,9 @@ #define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ #include "WarpXParticleContainer.H" -//#include -#include #include +#include +#include class CollisionType { @@ -16,25 +16,9 @@ public: const std::vector& species_names, std::string collision_name); - static void ElasticCollisionPerez( - amrex::DenseBins::index_type const I1s, - amrex::DenseBins::index_type const I1e, - amrex::DenseBins::index_type const I2s, - amrex::DenseBins::index_type const I2e, - amrex::DenseBins::index_type *I1, - amrex::DenseBins::index_type *I2, - amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, - amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, - const amrex::Real* w1, const amrex::Real* w2, - const amrex::Real q1, const amrex::Real q2, - const amrex::Real m1, const amrex::Real m2, - const amrex::Real T1, const amrex::Real T2, - const amrex::Real dt, const amrex::Real L, const amrex::Real V); - - static void ShuffleFisherYates( - amrex::DenseBins::index_type *array, - amrex::DenseBins::index_type const is, - amrex::DenseBins::index_type const ie); + static void doCoulombCollisionsWithinTile ( int lev, amrex::MFIter const& mfi, + std::unique_ptr& species1, + std::unique_ptr& species2 ); }; diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index 80f155a80..bd19c00cb 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -1,7 +1,123 @@ +#include #include "CollisionType.H" -#include "UpdateMomentumPerezElastic.H" +#include "ElasticCollisionPerez.H" + +using namespace amrex; +// Define shortcuts for frequently-used type names +using ParticleType = WarpXParticleContainer::ParticleType; +using ParticleTileType = WarpXParticleContainer::ParticleTileType; +using ParticleBins = DenseBins; +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 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; + } + +} + + +/* \brief Perform Coulomb collisions within one particle tile */ +void CollisionType::doCoulombCollisionsWithinTile + ( int lev, MFIter const& mfi, + std::unique_ptr& species_1, + std::unique_ptr& species_2 ) +{ -#include + // 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(); + Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); + Real* uy_1 = soa_1.GetRealData(PIdx::ux).data(); + Real* uz_1 = soa_1.GetRealData(PIdx::ux).data(); + Real* 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::ux).data(); + Real* uz_2 = soa_2.GetRealData(PIdx::ux).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) + + // 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, -1.0, -1.0, dt, -1.0, dV); + + } + ); + +} CollisionType::CollisionType( const std::vector& species_names, @@ -23,155 +139,4 @@ CollisionType::CollisionType( } -/* \brief Prepare information for and call - * UpdateMomentumPerezElastic(). - * i1s,i2s is the start index for I1,I2 (inclusive). - * i1e,i2e is the start index for I1,I2 (exclusive). - * I1 and I2 are the index arrays. - * u1 and u2 are the velocity arrays (u=v*gamma), - * they could be either different or the same, - * their lengths are not needed, - * I1 and I2 determine all elements that will be used. - * w1 and w2 are arrays of weights. - * q1 and q2 are charges. m1 and m2 are masses. - * T1 and T2 are temperatures (Joule) and will be used if greater than zero, - * otherwise will be computed. - * dt is the time step length between two collision calls. - * L is the Coulomb log and will be used if greater than zero, - * otherwise will be computed. - * V is the volume of the corresponding cell.*/ - -void CollisionType::ElasticCollisionPerez( - amrex::DenseBins::index_type const I1s, - amrex::DenseBins::index_type const I1e, - amrex::DenseBins::index_type const I2s, - amrex::DenseBins::index_type const I2e, - amrex::DenseBins::index_type *I1, - amrex::DenseBins::index_type *I2, - amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, - amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, - const amrex::Real* w1, const amrex::Real* w2, - const amrex::Real q1, const amrex::Real q2, - const amrex::Real m1, const amrex::Real m2, - const amrex::Real T1, const amrex::Real T2, - const amrex::Real dt, const amrex::Real L, const amrex::Real V) -{ - - constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); - int NI1 = I1e - I1s; - int NI2 = I2e - I2s; - // shuffle I1 and I2 - CollisionType::ShuffleFisherYates(I1, I1s, I1e); - CollisionType::ShuffleFisherYates(I2, I2s, I2e); - - // get local T1t and T2t - amrex::Real T1t; amrex::Real T2t; - if ( T1 <= 0. ) - { - amrex::Real vx = 0.; amrex::Real vy = 0.; - amrex::Real vz = 0.; amrex::Real vs = 0.; - amrex::Real gm = 0.; amrex::Real us = 0.; - for (int i1 = I1s; i1 < I1e; ++i1) - { - us = ( u1x[ I1[i1] ] * u1x[ I1[i1] ] + - u1y[ I1[i1] ] * u1y[ I1[i1] ] + - u1z[ I1[i1] ] * u1z[ I1[i1] ] ); - gm = std::sqrt(1.+us*inv_c2); - vx += u1x[ I1[i1] ] / gm; - vy += u1y[ I1[i1] ] / gm; - vz += u1z[ I1[i1] ] / gm; - vs += us / gm / gm; - } - vx = vx / NI1; vy = vy / NI1; - vz = vz / NI1; vs = vs / NI1; - T1t = m1/3.*(vs-(vx*vx+vy*vy+vz*vz)); - } - else { T1t = T1; } - if ( T2 <= 0. ) - { - amrex::Real vx = 0.; amrex::Real vy = 0.; - amrex::Real vz = 0.; amrex::Real vs = 0.; - amrex::Real gm = 0.; amrex::Real us = 0.; - for (int i2 = I2s; i2 < I2e; ++i2) - { - us = ( u2x[ I2[i2] ] * u2x[ I2[i2] ] + - u2y[ I2[i2] ] * u2y[ I2[i2] ] + - u2z[ I2[i2] ] * u2z[ I2[i2] ] ); - gm = std::sqrt(1.+us*inv_c2); - vx += u2x[ I2[i2] ] / gm; - vy += u2y[ I2[i2] ] / gm; - vz += u2z[ I2[i2] ] / gm; - vs += us / gm / gm; - } - vx = vx / NI2; vy = vy / NI2; - vz = vz / NI2; vs = vs / NI2; - T2t = m2/3.*(vs-(vx*vx+vy*vy+vz*vz)); - } - else { T2t = T2; } - - // local density - amrex::Real n1 = 0.; - amrex::Real n2 = 0.; - amrex::Real n12 = 0.; - for (int i1=I1s; i1::index_type *array, - amrex::DenseBins::index_type const is, - amrex::DenseBins::index_type const ie) -{ - int j; int buf; - for (int i = ie-1; i >= is+1; --i) - { - // get random number is <= j <= i - while ( true ) { - j = (int) floor(amrex::Random()*(i-is+1)) + is; - if ( j >= is && j <= i ) { break; } - } - buf = array[i]; - array[i] = array[j]; - array[j] = buf; - } -} diff --git a/Source/Particles/Collision/ElasticCollisionPerez.H b/Source/Particles/Collision/ElasticCollisionPerez.H new file mode 100644 index 000000000..53997e33c --- /dev/null +++ b/Source/Particles/Collision/ElasticCollisionPerez.H @@ -0,0 +1,144 @@ +#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ +#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ + +#include "WarpXParticleContainer.H" +//#include +#include +#include + +#include "ShuffleFisherYates.H" +#include "UpdateMomentumPerezElastic.H" + +/* \brief Prepare information for and call + * UpdateMomentumPerezElastic(). + * i1s,i2s is the start index for I1,I2 (inclusive). + * i1e,i2e is the start index for I1,I2 (exclusive). + * I1 and I2 are the index arrays. + * u1 and u2 are the velocity arrays (u=v*gamma), + * they could be either different or the same, + * their lengths are not needed, + * I1 and I2 determine all elements that will be used. + * w1 and w2 are arrays of weights. + * q1 and q2 are charges. m1 and m2 are masses. + * T1 and T2 are temperatures (Joule) and will be used if greater than zero, + * otherwise will be computed. + * dt is the time step length between two collision calls. + * L is the Coulomb log and will be used if greater than zero, + * otherwise will be computed. + * V is the volume of the corresponding cell.*/ + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ElasticCollisionPerez( + amrex::DenseBins::index_type const I1s, + amrex::DenseBins::index_type const I1e, + amrex::DenseBins::index_type const I2s, + amrex::DenseBins::index_type const I2e, + amrex::DenseBins::index_type *I1, + amrex::DenseBins::index_type *I2, + amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, + amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, + const amrex::Real* w1, const amrex::Real* w2, + const amrex::Real q1, const amrex::Real q2, + const amrex::Real m1, const amrex::Real m2, + const amrex::Real T1, const amrex::Real T2, + const amrex::Real dt, const amrex::Real L, const amrex::Real V) +{ + + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + int NI1 = I1e - I1s; + int NI2 = I2e - I2s; + + // shuffle I1 and I2 + ShuffleFisherYates(I1, I1s, I1e); + ShuffleFisherYates(I2, I2s, I2e); + + // get local T1t and T2t + amrex::Real T1t; amrex::Real T2t; + if ( T1 <= 0. ) + { + amrex::Real vx = 0.; amrex::Real vy = 0.; + amrex::Real vz = 0.; amrex::Real vs = 0.; + amrex::Real gm = 0.; amrex::Real us = 0.; + for (int i1 = I1s; i1 < I1e; ++i1) + { + us = ( u1x[ I1[i1] ] * u1x[ I1[i1] ] + + u1y[ I1[i1] ] * u1y[ I1[i1] ] + + u1z[ I1[i1] ] * u1z[ I1[i1] ] ); + gm = std::sqrt(1.+us*inv_c2); + vx += u1x[ I1[i1] ] / gm; + vy += u1y[ I1[i1] ] / gm; + vz += u1z[ I1[i1] ] / gm; + vs += us / gm / gm; + } + vx = vx / NI1; vy = vy / NI1; + vz = vz / NI1; vs = vs / NI1; + T1t = m1/3.*(vs-(vx*vx+vy*vy+vz*vz)); + } + else { T1t = T1; } + if ( T2 <= 0. ) + { + amrex::Real vx = 0.; amrex::Real vy = 0.; + amrex::Real vz = 0.; amrex::Real vs = 0.; + amrex::Real gm = 0.; amrex::Real us = 0.; + for (int i2 = I2s; i2 < I2e; ++i2) + { + us = ( u2x[ I2[i2] ] * u2x[ I2[i2] ] + + u2y[ I2[i2] ] * u2y[ I2[i2] ] + + u2z[ I2[i2] ] * u2z[ I2[i2] ] ); + gm = std::sqrt(1.+us*inv_c2); + vx += u2x[ I2[i2] ] / gm; + vy += u2y[ I2[i2] ] / gm; + vz += u2z[ I2[i2] ] / gm; + vs += us / gm / gm; + } + vx = vx / NI2; vy = vy / NI2; + vz = vz / NI2; vs = vs / NI2; + T2t = m2/3.*(vs-(vx*vx+vy*vy+vz*vz)); + } + else { T2t = T2; } + + // local density + amrex::Real n1 = 0.; + amrex::Real n2 = 0.; + amrex::Real n12 = 0.; + for (int i1=I1s; i1 +#include +#include + +/* \brief Shuffle array according to Fisher-Yates algorithm. + * Only shuffle the part between is <= i < ie, n = ie-is.*/ + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ShuffleFisherYates( + amrex::DenseBins::index_type *array, + amrex::DenseBins::index_type const is, + amrex::DenseBins::index_type const ie) +{ + int j; + amrex::DenseBins::index_type buf; + for (int i = (int) ie-1; i >= (int) is+1; --i) + { + // get random number is <= j <= i + while ( true ) { + j = (int) floor(amrex::Random()*(i-is+1)) + is; + if ( j >= is && j <= i ) { break; } + } + 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 index 156f03b19..015d9a3a1 100644 --- a/Source/Particles/Collision/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -2,9 +2,9 @@ #define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ #include "WarpXParticleContainer.H" -#include #include #include +#include /* \brief Update particle velocities according to * F. Perez et al., Phys. Plasmas 19, 083104 (2012). diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 240b94273..0efd60e91 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -2,8 +2,6 @@ #define WARPX_ParticleContainer_H_ #include "ElementaryProcess.H" -#include "CoulombCollisions.H" -#include "CollisionType.H" #include #include @@ -17,6 +15,8 @@ #include #endif +#include "CollisionType.H" + #include #include #include diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 3d8dc0ddc..55da8b938 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -56,6 +56,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) ionization_process = IonizationProcess(); // collision + allcollisions.resize(ncollisions); for (int i = 0; i < ncollisions; ++i) { allcollisions[i].reset(new CollisionType(species_names, collision_names[i])); } @@ -662,7 +663,8 @@ MultiParticleContainer::doCoulombCollisions () #endif for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ - doCoulombCollisionsWithinTile( lev, mfi, species1, species2 ); + CollisionType::doCoulombCollisionsWithinTile + ( lev, mfi, species1, species2 ); } } -- cgit v1.2.3 From e4d1eb5417129873b631cd248652191a8ff2d345 Mon Sep 17 00:00:00 2001 From: Yinjian Zhao Date: Mon, 18 Nov 2019 12:48:53 -0800 Subject: Add collisions between a single species. --- Source/Particles/Collision/CollisionType.H | 11 +- Source/Particles/Collision/CollisionType.cpp | 196 ++++++++++++++------- Source/Particles/Collision/ElasticCollisionPerez.H | 106 ++++++----- Source/Particles/MultiParticleContainer.cpp | 7 +- 4 files changed, 190 insertions(+), 130 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H index c9ba643d8..ef0d89464 100644 --- a/Source/Particles/Collision/CollisionType.H +++ b/Source/Particles/Collision/CollisionType.H @@ -9,16 +9,19 @@ class CollisionType { public: - int m_species1; - int m_species2; + int m_species1_index; + int m_species2_index; + bool m_isSameSpecies; CollisionType( const std::vector& species_names, std::string collision_name); - static void doCoulombCollisionsWithinTile ( int lev, amrex::MFIter const& mfi, + static void doCoulombCollisionsWithinTile ( + int lev, amrex::MFIter const& mfi, std::unique_ptr& species1, - std::unique_ptr& species2 ); + std::unique_ptr& species2, + bool isSameSpecies ); }; diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index 22cc3bb8a..fc14e9df6 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -1,6 +1,7 @@ -#include #include "CollisionType.H" +#include "ShuffleFisherYates.H" #include "ElasticCollisionPerez.H" +#include using namespace amrex; // Define shortcuts for frequently-used type names @@ -50,72 +51,128 @@ namespace { void CollisionType::doCoulombCollisionsWithinTile ( int lev, MFIter const& mfi, std::unique_ptr& species_1, - std::unique_ptr& species_2 ) + std::unique_ptr& species_2, + bool isSameSpecies ) { - // 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(); - Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); - Real* uy_1 = soa_1.GetRealData(PIdx::uy).data(); - Real* uz_1 = soa_1.GetRealData(PIdx::uz).data(); - Real* 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) - - // 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, 50.0, 50.0, dt, 2.0, dV); - - } - ); + 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(); + Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); + Real* uy_1 = soa_1.GetRealData(PIdx::uy).data(); + Real* uz_1 = soa_1.GetRealData(PIdx::uz).data(); + Real* 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]; + + // shuffle + ShuffleFisherYates( + indices_1, cell_start_1, (cell_start_1+cell_stop_1)/2+1 ); + + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start_1, (cell_start_1+cell_stop_1)/2, + (cell_start_1+cell_stop_1)/2, 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, -1.0, -1.0, dt, 2.0, 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(); + Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); + Real* uy_1 = soa_1.GetRealData(PIdx::uy).data(); + Real* uz_1 = soa_1.GetRealData(PIdx::uz).data(); + Real* 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) + + // 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, -1.0, -1.0, dt, 2.0, dV); + } + ); + } // end if ( isSameSpecies) } @@ -132,11 +189,16 @@ CollisionType::CollisionType( for (int i=0; i -#include -#include +//#include "WarpXParticleContainer.H" +//#include +//#include -#include "ShuffleFisherYates.H" #include "UpdateMomentumPerezElastic.H" +#include +#include /* \brief Prepare information for and call * UpdateMomentumPerezElastic(). @@ -20,55 +20,50 @@ * I1 and I2 determine all elements that will be used. * w1 and w2 are arrays of weights. * q1 and q2 are charges. m1 and m2 are masses. - * T1 and T2 are temperatures (Joule) and will be used if greater than zero, + * T1 and T2 are temperatures (Joule) + * and will be used if greater than zero, * otherwise will be computed. * dt is the time step length between two collision calls. * L is the Coulomb log and will be used if greater than zero, * otherwise will be computed. - * V is the volume of the corresponding cell.*/ + * V is the volume of the corresponding cell. +*/ +template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void ElasticCollisionPerez( - amrex::DenseBins::index_type const I1s, - amrex::DenseBins::index_type const I1e, - amrex::DenseBins::index_type const I2s, - amrex::DenseBins::index_type const I2e, - amrex::DenseBins::index_type *I1, - amrex::DenseBins::index_type *I2, - amrex::Real* u1x, amrex::Real* u1y, amrex::Real* u1z, - amrex::Real* u2x, amrex::Real* u2y, amrex::Real* u2z, - const amrex::Real* w1, const amrex::Real* w2, - const amrex::Real q1, const amrex::Real q2, - const amrex::Real m1, const amrex::Real m2, - const amrex::Real T1, const amrex::Real T2, - const amrex::Real dt, const amrex::Real L, const amrex::Real V) +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 V) { - constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); - int NI1 = (int) I1e - I1s; - int NI2 = (int) I2e - I2s; - - // shuffle I1 and I2 - ShuffleFisherYates(I1, I1s, I1e); - ShuffleFisherYates(I2, I2s, I2e); + T_R constexpr inv_c2 = 1./(PhysConst::c*PhysConst::c); + int NI1 = I1e - I1s; + int NI2 = I2e - I2s; // get local T1t and T2t - amrex::Real T1t; amrex::Real T2t; + T_R T1t; T_R T2t; if ( T1 <= 0. ) { - amrex::Real vx = 0.; amrex::Real vy = 0.; - amrex::Real vz = 0.; amrex::Real vs = 0.; - amrex::Real gm = 0.; amrex::Real us = 0.; - for (int i1 = (int) I1s; i1 < (int) I1e; ++i1) + T_R vx = 0.; T_R vy = 0.; + T_R vz = 0.; T_R vs = 0.; + T_R gm = 0.; T_R us = 0.; + for (int i1 = I1s; i1 < I1e; ++i1) { - int index = (int) I1[i1]; - us = ( u1x[ index ] * u1x[ index ] + - u1y[ index ] * u1y[ index ] + - u1z[ index ] * u1z[ index ] ); + us = ( u1x[ I1[i1] ] * u1x[ I1[i1] ] + + u1y[ I1[i1] ] * u1y[ I1[i1] ] + + u1z[ I1[i1] ] * u1z[ I1[i1] ] ); gm = std::sqrt(1.+us*inv_c2); - vx += u1x[ index ] / gm; - vy += u1y[ index ] / gm; - vz += u1z[ index ] / gm; + vx += u1x[ I1[i1] ] / gm; + vy += u1y[ I1[i1] ] / gm; + vz += u1z[ I1[i1] ] / gm; vs += us / gm / gm; } if ( NI1 == 0 ) { T1t = 0.; } @@ -82,19 +77,18 @@ void ElasticCollisionPerez( else { T1t = T1; } if ( T2 <= 0. ) { - amrex::Real vx = 0.; amrex::Real vy = 0.; - amrex::Real vz = 0.; amrex::Real vs = 0.; - amrex::Real gm = 0.; amrex::Real us = 0.; - for (int i2 = (int) I2s; i2 < (int) I2e; ++i2) + T_R vx = 0.; T_R vy = 0.; + T_R vz = 0.; T_R vs = 0.; + T_R gm = 0.; T_R us = 0.; + for (int i2 = I2s; i2 < I2e; ++i2) { - int index = (int) I2[i2]; - us = ( u2x[ index ] * u2x[ index ] + - u2y[ index ] * u2y[ index ] + - u2z[ index ] * u2z[ index ] ); + us = ( u2x[ I2[i2] ] * u2x[ I2[i2] ] + + u2y[ I2[i2] ] * u2y[ I2[i2] ] + + u2z[ I2[i2] ] * u2z[ I2[i2] ] ); gm = std::sqrt(1.+us*inv_c2); - vx += u2x[ index ] / gm; - vy += u2y[ index ] / gm; - vz += u2z[ index ] / gm; + vx += u2x[ I2[i2] ] / gm; + vy += u2y[ I2[i2] ] / gm; + vz += u2z[ I2[i2] ] / gm; vs += us / gm / gm; } if ( NI2 == 0 ) { T2t = 0.; } @@ -108,9 +102,9 @@ void ElasticCollisionPerez( else { T2t = T2; } // local density - amrex::Real n1 = 0.; - amrex::Real n2 = 0.; - amrex::Real n12 = 0.; + T_R n1 = 0.; + T_R n2 = 0.; + T_R n12 = 0.; for (int i1=I1s; i1m_species1 ]; - auto& species2 = allcontainers[ allcollisions[i]->m_species2 ]; + auto& species1 = allcontainers[ allcollisions[i]->m_species1_index ]; + auto& species2 = allcontainers[ allcollisions[i]->m_species2_index ]; // Enable tiling MFItInfo info; @@ -656,7 +656,8 @@ MultiParticleContainer::doCoulombCollisions () for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ CollisionType::doCoulombCollisionsWithinTile - ( lev, mfi, species1, species2 ); + ( lev, mfi, species1, species2, + allcollisions[i]->m_isSameSpecies ); } } -- cgit v1.2.3 From a9b158c647a20d030c98c18636d8bdb81d6dda35 Mon Sep 17 00:00:00 2001 From: Yinjian Zhao Date: Sun, 24 Nov 2019 13:27:37 -0700 Subject: Modified the collision part of MultiParticleContainer.cpp. --- Source/Particles/MultiParticleContainer.cpp | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index f6f9ab1a2..37b595996 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -58,7 +58,8 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // collision allcollisions.resize(ncollisions); for (int i = 0; i < ncollisions; ++i) { - allcollisions[i].reset(new CollisionType(species_names, collision_names[i])); + allcollisions[i].reset + (new CollisionType(species_names, collision_names[i])); } } @@ -632,33 +633,27 @@ MultiParticleContainer::doCoulombCollisions () for (int i = 0; i < ncollisions; ++i) { - // At this point, the code always collides the first and second species - // TODO: Read the user input to read the different types of collisions requested - // and loop over all types of collisions, selecting each time the - // two types of species that will be collided - // auto& species1 = allcontainers[0]; - // auto& species2 = allcontainers[1]; auto& species1 = allcontainers[ allcollisions[i]->m_species1_index ]; auto& species2 = allcontainers[ allcollisions[i]->m_species2_index ]; - + // Enable tiling MFItInfo info; if (Gpu::notInLaunchRegion()) info.EnableTiling(species1->tile_size); - + // Loop over refinement levels for (int lev = 0; lev <= species1->finestLevel(); ++lev){ - + // Loop over all grids/tiles at this level #ifdef _OPENMP info.SetDynamic(true); #pragma omp parallel #endif for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ - + CollisionType::doCoulombCollisionsWithinTile ( lev, mfi, species1, species2, allcollisions[i]->m_isSameSpecies ); - + } } } -- cgit v1.2.3 From bac7cea9ffe013c7ce84be8d8b5550d9b35b2d4f Mon Sep 17 00:00:00 2001 From: Yinjian Zhao Date: Mon, 25 Nov 2019 14:44:44 -0700 Subject: Add user inputs of Coulomb Log of each collision type. --- Source/Particles/Collision/CollisionType.H | 3 ++- Source/Particles/Collision/CollisionType.cpp | 12 ++++++++---- Source/Particles/MultiParticleContainer.cpp | 3 ++- 3 files changed, 12 insertions(+), 6 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H index ef0d89464..fc764c8c5 100644 --- a/Source/Particles/Collision/CollisionType.H +++ b/Source/Particles/Collision/CollisionType.H @@ -12,6 +12,7 @@ public: int m_species1_index; int m_species2_index; bool m_isSameSpecies; + amrex::Real m_CoulombLog; CollisionType( const std::vector& species_names, @@ -21,7 +22,7 @@ public: int lev, amrex::MFIter const& mfi, std::unique_ptr& species1, std::unique_ptr& species2, - bool isSameSpecies ); + bool isSameSpecies, amrex::Real CoulombLog ); }; diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index b0d6cd0d1..2a26e0108 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -13,6 +13,10 @@ CollisionType::CollisionType( amrex::ParmParse pp(collision_name); pp.getarr("species", collision_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_1, std::unique_ptr& species_2, - bool isSameSpecies ) + bool isSameSpecies, Real CoulombLog ) { - if ( isSameSpecies) // species_1 == species_2 + if ( isSameSpecies ) // species_1 == species_2 { // Extract particles in the tile that `mfi` points to ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); @@ -130,7 +134,7 @@ void CollisionType::doCoulombCollisionsWithinTile 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, -1.0, -1.0, dt, 2.0, dV); + q1, q1, m1, m1, -1.0, -1.0, dt, CoulombLog, dV); } } ); @@ -203,7 +207,7 @@ void CollisionType::doCoulombCollisionsWithinTile 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, -1.0, -1.0, dt, 2.0, dV); + q1, q2, m1, m2, -1.0, -1.0, dt, CoulombLog, dV); } } ); diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 37b595996..a090fdbec 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -652,7 +652,8 @@ MultiParticleContainer::doCoulombCollisions () CollisionType::doCoulombCollisionsWithinTile ( lev, mfi, species1, species2, - allcollisions[i]->m_isSameSpecies ); + allcollisions[i]->m_isSameSpecies, + allcollisions[i]->m_CoulombLog ); } } -- cgit v1.2.3 From 38b593dcb880fc292caa5fb3e2c75baa5d0b61d3 Mon Sep 17 00:00:00 2001 From: Yinjian Zhao Date: Wed, 18 Dec 2019 08:34:45 -0800 Subject: Modify based on suggestions. --- Source/Particles/Collision/CollisionType.cpp | 52 +++++++++++++--------- Source/Particles/Collision/ElasticCollisionPerez.H | 2 +- .../Collision/UpdateMomentumPerezElastic.H | 8 ++-- Source/Particles/MultiParticleContainer.cpp | 2 +- 4 files changed, 36 insertions(+), 28 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp index 6e21cc327..b8014579d 100644 --- a/Source/Particles/Collision/CollisionType.cpp +++ b/Source/Particles/Collision/CollisionType.cpp @@ -8,18 +8,11 @@ CollisionType::CollisionType( std::string const collision_name) { - // check the coordinate is Cartesian - amrex::ParmParse pp_check("geometry"); - int coord_sys; - pp_check.query("coord_sys",coord_sys); - if (coord_sys != 0) - { amrex::Abort("Collision only works for Cartesian coordinate for now."); } - - // check the dimension is 3 - std::vector prob_lo; - pp_check.getarr("prob_lo", prob_lo); - if (prob_lo.size() != 3) - { amrex::Abort("Collision only works for 3D for now."); } +#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 collision_species; @@ -90,8 +83,15 @@ namespace { } - -/* \brief Perform Coulomb collisions within one particle tile */ +/** 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& species_1, @@ -113,10 +113,14 @@ void CollisionType::doCoulombCollisionsWithinTile int const n_cells = bins_1.numBins(); // - Species 1 auto& soa_1 = ptile_1.GetStructOfArrays(); - Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); - Real* uy_1 = soa_1.GetRealData(PIdx::uy).data(); - Real* uz_1 = soa_1.GetRealData(PIdx::uz).data(); - Real* w_1 = soa_1.GetRealData(PIdx::w).data(); + 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(); @@ -171,10 +175,14 @@ void CollisionType::doCoulombCollisionsWithinTile int const n_cells = bins_1.numBins(); // - Species 1 auto& soa_1 = ptile_1.GetStructOfArrays(); - Real* ux_1 = soa_1.GetRealData(PIdx::ux).data(); - Real* uy_1 = soa_1.GetRealData(PIdx::uy).data(); - Real* uz_1 = soa_1.GetRealData(PIdx::uz).data(); - Real* w_1 = soa_1.GetRealData(PIdx::w).data(); + 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(); diff --git a/Source/Particles/Collision/ElasticCollisionPerez.H b/Source/Particles/Collision/ElasticCollisionPerez.H index d347802a0..8e16d95cc 100644 --- a/Source/Particles/Collision/ElasticCollisionPerez.H +++ b/Source/Particles/Collision/ElasticCollisionPerez.H @@ -6,7 +6,7 @@ #include #include -/* \brief Prepare information for and call +/** \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). diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H index efc9be465..948e8b075 100644 --- a/Source/Particles/Collision/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -9,8 +9,8 @@ /* \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). - * LmdD is max(Debye length, minimal interparticle distance). - * L is the Coulomb log. A fixed L will be used if L > 0, + * @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. @@ -229,7 +229,7 @@ void UpdateMomentumPerezElastic ( // Rejection method r = amrex::Random(); - if ( w2/amrex::max(w1,w2) > r ) + if ( w2 > r*amrex::max(w1, w2) ) { u1x = p1fx / m1; u1y = p1fy / m1; @@ -238,7 +238,7 @@ void UpdateMomentumPerezElastic ( AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z)); } r = amrex::Random(); - if ( w1/amrex::max(w1,w2) > r ) + if ( w1 > r*amrex::max(w1, w2) ) { u2x = p2fx / m2; u2y = p2fy / m2; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 29b5d2069..d84bc1afa 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -646,7 +646,7 @@ MultiParticleContainer::doCoulombCollisions () // Loop over all grids/tiles at this level #ifdef _OPENMP info.SetDynamic(true); - #pragma omp parallel + #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ -- cgit v1.2.3