diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 1 | ||||
| -rw-r--r-- | Source/Particles/Collision/CollisionType.H | 26 | ||||
| -rw-r--r-- | Source/Particles/Collision/CollisionType.cpp | 271 | ||||
| -rw-r--r-- | Source/Particles/Collision/ElasticCollisionPerez.H | 144 | ||||
| -rw-r--r-- | Source/Particles/Collision/Make.package | 4 | ||||
| -rw-r--r-- | Source/Particles/Collision/ShuffleFisherYates.H | 33 | ||||
| -rw-r--r-- | Source/Particles/Collision/UpdateMomentumPerezElastic.H | 2 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 4 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 4 |
9 files changed, 309 insertions, 180 deletions
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 <WarpX.H> -#include <AMReX_REAL.H> #include <AMReX_DenseBins.H> +#include <AMReX_REAL.H> +#include <AMReX_ParmParse.H> class CollisionType { @@ -16,25 +16,9 @@ public: const std::vector<std::string>& species_names, std::string collision_name); - static void ElasticCollisionPerez( - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I1s, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I1e, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2s, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2e, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type *I1, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::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<WarpXParticleContainer::ParticleType>::index_type *array, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const is, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const ie); + static void doCoulombCollisionsWithinTile ( int lev, amrex::MFIter const& mfi, + std::unique_ptr<WarpXParticleContainer>& species1, + std::unique_ptr<WarpXParticleContainer>& 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 <WarpX.H> #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<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 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<WarpXParticleContainer>& species_1, + std::unique_ptr<WarpXParticleContainer>& species_2 ) +{ -#include <AMReX_ParmParse.H> + // 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<std::string>& 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<WarpXParticleContainer::ParticleType>::index_type const I1s, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I1e, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2s, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2e, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type *I1, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::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<I1e; ++i1) { n1 += w1[ I1[i1] ]; } - for (int i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; } - n1 = n1 / V; n2 = n2 / V; - { - 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 / V; - } - - // compute Debye length lmdD - amrex::Real lmdD; - lmdD = 1./std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + - n2*q2*q2/(T2t*PhysConst::ep0) ); - amrex::Real rmin = - std::pow(4.*MathConst::pi/3.*amrex::max(n1,n2),-1./3.); - 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; } - } - } - -} - -/* \brief Shuffle array according to Fisher-Yates algorithm. - * Only shuffle the part between is <= i < ie, n = ie-is.*/ - -void CollisionType::ShuffleFisherYates( - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type *array, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const is, - amrex::DenseBins<WarpXParticleContainer::ParticleType>::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 <WarpX.H> +#include <AMReX_DenseBins.H> +#include <AMReX_REAL.H> + +#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<WarpXParticleContainer::ParticleType>::index_type const I1s, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I1e, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2s, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const I2e, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type *I1, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::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<I1e; ++i1) { n1 += w1[ I1[i1] ]; } + for (int i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; } + n1 = n1 / V; n2 = n2 / V; + { + 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 / V; + } + + // compute Debye length lmdD + amrex::Real lmdD; + lmdD = 1./std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + + n2*q2*q2/(T2t*PhysConst::ep0) ); + amrex::Real rmin = + std::pow(4.*MathConst::pi/3.*amrex::max(n1,n2),-1./3.); + 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 index 3134a29b8..a302f0c02 100644 --- a/Source/Particles/Collision/Make.package +++ b/Source/Particles/Collision/Make.package @@ -1,9 +1,9 @@ CEXE_headers += CollisionType.H -CEXE_headers += CoulombCollisions.H +CEXE_headers += ElasticCollisionPerez.H +CEXE_headers += ShuffleFisherYates.H CEXE_headers += UpdateMomentumPerezElastic.H CEXE_sources += CollisionType.cpp -CEXE_sources += CoulombCollisions.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..a3249fca8 --- /dev/null +++ b/Source/Particles/Collision/ShuffleFisherYates.H @@ -0,0 +1,33 @@ +#ifndef WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ +#define WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ + +#include "WarpXParticleContainer.H" +//#include <WarpX.H> +#include <AMReX_DenseBins.H> +#include <AMReX_REAL.H> + +/* \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<WarpXParticleContainer::ParticleType>::index_type *array, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const is, + amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type const ie) +{ + int j; + amrex::DenseBins<WarpXParticleContainer::ParticleType>::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 <WarpX.H> #include <AMReX_REAL.H> #include <AMReX_DenseBins.H> +#include <WarpXConst.H> /* \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 <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> @@ -17,6 +15,8 @@ #include <QuantumSyncEngineWrapper.H> #endif +#include "CollisionType.H" + #include <memory> #include <map> #include <string> 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 ); } } |