diff options
Diffstat (limited to 'Source')
33 files changed, 1019 insertions, 78 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 572030f73..51439430d 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -752,14 +752,14 @@ PML::CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp, void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) +PML::ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain) { ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) +PML::ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); @@ -844,11 +844,10 @@ PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, void PML::CopyToPML (MultiFab& pml, MultiFab& reg, const Geometry& geom) { - const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); const auto& period = geom.periodicity(); - pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period); + pml.ParallelCopy(reg, 0, 0, 1, IntVect(0), ngp, period); } void diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index b90df6945..c79a12066 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -119,6 +119,7 @@ private: /** This function saves the plot file * * @param[in] pc WarpX particle container + * @param[in] name species name * @param[in] iteration timestep * @param[in] write_real_comp The real attribute ids, from WarpX * @param[in] real_comp_names The real attribute names, from WarpX diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index b5fd52bdc..f5491ffe3 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} @@ -475,7 +476,7 @@ WarpX::OneStep_sub1 (Real curtime) } void -WarpX::PushParticlesandDepose (Real cur_time) +WarpX::PushParticlesandDepose (amrex::Real cur_time) { // Evolve particles to p^{n+1/2} and x^{n+1} // Depose current, j^{n+1/2} @@ -485,7 +486,7 @@ WarpX::PushParticlesandDepose (Real cur_time) } void -WarpX::PushParticlesandDepose (int lev, Real cur_time, DtType a_dt_type) +WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) { mypc->Evolve(lev, *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], @@ -583,7 +584,7 @@ WarpX::ComputeDt () * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void -WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ +WarpX::computeMaxStepBoostAccelerator(const amrex::Geometry& a_geom){ // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 8f0853484..edd4df34d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -3,10 +3,10 @@ using namespace amrex; /* \brief Initialize fields in spectral space, and FFT plans */ -SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, - const SpectralKSpace& k_space, - const DistributionMapping& dm, - const int n_field_required ) +SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const amrex::DistributionMapping& dm, + const int n_field_required ) { const BoxArray& spectralspace_ba = k_space.spectralspace_ba; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 9807665c6..4848b051e 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -90,7 +90,7 @@ WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) #endif void -WarpX::EvolveB (Real a_dt) +WarpX::EvolveB (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { EvolveB(lev, a_dt); @@ -98,7 +98,7 @@ WarpX::EvolveB (Real a_dt) } void -WarpX::EvolveB (int lev, Real a_dt) +WarpX::EvolveB (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); EvolveB(lev, PatchType::fine, a_dt); @@ -303,7 +303,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveE (Real a_dt) +WarpX::EvolveE (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -312,7 +312,7 @@ WarpX::EvolveE (Real a_dt) } void -WarpX::EvolveE (int lev, Real a_dt) +WarpX::EvolveE (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); EvolveE(lev, PatchType::fine, a_dt); @@ -611,7 +611,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveF (Real a_dt, DtType a_dt_type) +WarpX::EvolveF (amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -622,7 +622,7 @@ WarpX::EvolveF (Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -631,7 +631,7 @@ WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 5d3c14c14..93729a0ae 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -144,7 +144,7 @@ void Filter::DoFilter (const Box& tbx, * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil()"); ncomp = std::min(ncomp, srcmf.nComp()); @@ -179,8 +179,8 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, - const Box& tbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, + const amrex::Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index e32719adf..48aaac275 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -1,3 +1,6 @@ +#ifndef WARPX_F_H_ +#define WARPX_F_H_ + #include <AMReX_BLFort.H> #ifdef __cplusplus @@ -130,3 +133,5 @@ extern "C" #ifdef __cplusplus } #endif + +#endif //WARPX_F_H_ diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 300c41f7b..261a14fc5 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -21,7 +21,7 @@ namespace { amrex::Abort(string.c_str()); } - Real parseChargeName(const ParmParse pp, const std::string& name) { + Real parseChargeName(const ParmParse& pp, const std::string& name) { Real result; if (name == "q_e") { return PhysConst::q_e; @@ -33,13 +33,13 @@ namespace { } } - Real parseChargeString(const ParmParse pp, const std::string& name) { + Real parseChargeString(const ParmParse& pp, const std::string& name) { if(name.substr(0, 1) == "-") return -1.0 * parseChargeName(pp, name.substr(1, name.size() - 1)); return parseChargeName(pp, name); } - Real parseMassString(const ParmParse pp, const std::string& name) { + Real parseMassString(const ParmParse& pp, const std::string& name) { Real result; if (name == "m_e") { return PhysConst::m_e; @@ -109,6 +109,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) part_pos_s.end(), part_pos_s.begin(), ::tolower); + num_particles_per_cell_each_dim.assign(3, 0); if (part_pos_s == "python") { return; } else if (part_pos_s == "singleparticle") { diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 63ace31fb..1e83ca02f 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -10,6 +10,17 @@ #include <memory> #include <limits> +/** + * The main method to inject a laser pulse in WarpX is to use an artificial + * antenna: particles evenly distributed in a given plane (one particle per + * cell) move at each iteration and deposit a current J onto the grid, which + * in turns creates an electromagnetic field on the grid. The particles' + * displacements are prescribed to create the field requested by the user. + * + * These artificial particles are contained in the LaserParticleContainer. + * LaserParticleContainer derives directly from WarpXParticleContainer. It + * requires a DepositCurrent function, but no FieldGather function. + */ class LaserParticleContainer : public WarpXParticleContainer { diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H index 528309492..e0ec0dc28 100644 --- a/Source/Laser/LaserProfiles.H +++ b/Source/Laser/LaserProfiles.H @@ -45,7 +45,7 @@ public: /** Initialize Laser Profile * * Reads the section of the inputfile relative to the laser beam - * (e.g. <laser_name>.profile_t_peak, <laser_name>.profile_duration...) + * (e.g. laser_name.profile_t_peak, laser_name.profile_duration...) * and the "my_constants" section. It also receives some common * laser profile parameters. It uses these data to initialize the * member variables of the laser profile class. @@ -64,6 +64,7 @@ public: * * Xp, Yp and amplitude must be arrays with the same length * + * @param[in] np number of antenna particles * @param[in] Xp X coordinate of the particles of the antenna * @param[in] Yp Y coordinate of the particles of the antenna * @param[in] t time (seconds) diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H new file mode 100644 index 000000000..d020f47e8 --- /dev/null +++ b/Source/Particles/Collision/CollisionType.H @@ -0,0 +1,39 @@ +#ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ +#define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ + +#include "WarpXParticleContainer.H" +#include <AMReX_DenseBins.H> +#include <AMReX_REAL.H> +#include <AMReX_ParmParse.H> + +class CollisionType +{ +public: + int m_species1_index; + int m_species2_index; + bool m_isSameSpecies; + amrex::Real m_CoulombLog; + + CollisionType( + const std::vector<std::string>& species_names, + std::string const collision_name); + + /** Perform all binary collisions within a tile + * + * @param lev AMR level of the tile + * @param mfi iterator for multifab + * @param species1/2 pointer to species container + * @param isSameSpecies true if collision is between same species + * @param CoulombLog user input Coulomb logrithm + * + */ + + static void doCoulombCollisionsWithinTile ( + int const lev, amrex::MFIter const& mfi, + std::unique_ptr<WarpXParticleContainer>& species1, + std::unique_ptr<WarpXParticleContainer>& species2, + bool const isSameSpecies, amrex::Real const CoulombLog ); + +}; + +#endif // WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp new file mode 100644 index 000000000..b8014579d --- /dev/null +++ b/Source/Particles/Collision/CollisionType.cpp @@ -0,0 +1,241 @@ +#include "CollisionType.H" +#include "ShuffleFisherYates.H" +#include "ElasticCollisionPerez.H" +#include <WarpX.H> + +CollisionType::CollisionType( + const std::vector<std::string>& species_names, + std::string const collision_name) +{ + +#if defined WARPX_DIM_XZ + amrex::Abort("Collisions only work in 3D geometry for now."); +#elif defined WARPX_DIM_RZ + amrex::Abort("Collisions only work in Cartesian geometry for now."); +#endif + + // read collision species + std::vector<std::string> collision_species; + amrex::ParmParse pp(collision_name); + pp.getarr("species", collision_species); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(collision_species.size() == 2, + "Collision species must name exactly two species."); + + // default Coulomb log, if < 0, will be computed automatically + m_CoulombLog = -1.0; + pp.query("CoulombLog", m_CoulombLog); + + for (int i=0; i<species_names.size(); i++) + { + if (species_names[i] == collision_species[0]) + { m_species1_index = i; } + if (species_names[i] == collision_species[1]) + { m_species2_index = i; } + } + + if (collision_species[0] == collision_species[1]) + m_isSameSpecies = true; + else + m_isSameSpecies = false; + +} + +using namespace amrex; +// Define shortcuts for frequently-used type names +using ParticleType = WarpXParticleContainer::ParticleType; +using ParticleTileType = WarpXParticleContainer::ParticleTileType; +using ParticleBins = DenseBins<ParticleType>; +using index_type = ParticleBins::index_type; + +namespace { + + /* Find the particles and count the particles that are in each cell. + Note that this does *not* rearrange particle arrays */ + ParticleBins + findParticlesInEachCell( int const lev, MFIter const& mfi, + ParticleTileType const& ptile) { + + // Extract particle structures for this tile + int const np = ptile.numParticles(); + ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data(); + + // Extract box properties + Geometry const& geom = WarpX::GetInstance().Geom(lev); + Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box + const auto lo = lbound(cbx); + const auto dxi = geom.InvCellSizeArray(); + const auto plo = geom.ProbLoArray(); + + // Find particles that are in each cell ; + // results are stored in the object `bins`. + ParticleBins bins; + bins.build(np, particle_ptr, cbx, + // Pass lambda function that returns the cell index + [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect + { + return IntVect(AMREX_D_DECL((p.pos(0)-plo[0])*dxi[0] - lo.x, + (p.pos(1)-plo[1])*dxi[1] - lo.y, + (p.pos(2)-plo[2])*dxi[2] - lo.z)); + }); + + return bins; + } + +} + +/** Perform all binary collisions within a tile + * + * @param lev AMR level of the tile + * @param mfi iterator for multifab + * @param species1/2 pointer to species container + * @param isSameSpecies true if collision is between same species + * @param CoulombLog user input Coulomb logrithm + * + */ +void CollisionType::doCoulombCollisionsWithinTile + ( int const lev, MFIter const& mfi, + std::unique_ptr<WarpXParticleContainer>& species_1, + std::unique_ptr<WarpXParticleContainer>& species_2, + bool const isSameSpecies, Real const CoulombLog ) +{ + + if ( isSameSpecies ) // species_1 == species_2 + { + // Extract particles in the tile that `mfi` points to + ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); + + // Find the particles that are in each cell of this tile + ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 ); + + // Loop over cells, and collide the particles in each cell + + // Extract low-level data + int const n_cells = bins_1.numBins(); + // - Species 1 + auto& soa_1 = ptile_1.GetStructOfArrays(); + ParticleReal * const AMREX_RESTRICT ux_1 = + soa_1.GetRealData(PIdx::ux).data(); + ParticleReal * const AMREX_RESTRICT uy_1 = + soa_1.GetRealData(PIdx::uy).data(); + ParticleReal * const AMREX_RESTRICT uz_1 = + soa_1.GetRealData(PIdx::uz).data(); + ParticleReal const * const AMREX_RESTRICT w_1 = + soa_1.GetRealData(PIdx::w).data(); + index_type* indices_1 = bins_1.permutationPtr(); + index_type const* cell_offsets_1 = bins_1.offsetsPtr(); + Real q1 = species_1->getCharge(); + Real m1 = species_1->getMass(); + + const Real dt = WarpX::GetInstance().getdt(lev); + Geometry const& geom = WarpX::GetInstance().Geom(lev); + const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); + + // Loop over cells + amrex::ParallelFor( n_cells, + [=] AMREX_GPU_DEVICE (int i_cell) noexcept + { + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + index_type const cell_start_1 = cell_offsets_1[i_cell]; + index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; + index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2; + + // Do not collide if there is only one particle in the cell + if ( cell_stop_1 - cell_start_1 >= 2 ) + { + // shuffle + ShuffleFisherYates( + indices_1, cell_start_1, cell_half_1 ); + + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start_1, cell_half_1, + cell_half_1, cell_stop_1, + indices_1, indices_1, + ux_1, uy_1, uz_1, ux_1, uy_1, uz_1, w_1, w_1, + q1, q1, m1, m1, Real(-1.0), Real(-1.0), + dt, CoulombLog, dV ); + } + } + ); + } + else // species_1 != species_2 + { + // Extract particles in the tile that `mfi` points to + ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); + ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi); + + // Find the particles that are in each cell of this tile + ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 ); + ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 ); + + // Loop over cells, and collide the particles in each cell + + // Extract low-level data + int const n_cells = bins_1.numBins(); + // - Species 1 + auto& soa_1 = ptile_1.GetStructOfArrays(); + ParticleReal * const AMREX_RESTRICT ux_1 = + soa_1.GetRealData(PIdx::ux).data(); + ParticleReal * const AMREX_RESTRICT uy_1 = + soa_1.GetRealData(PIdx::uy).data(); + ParticleReal * const AMREX_RESTRICT uz_1 = + soa_1.GetRealData(PIdx::uz).data(); + ParticleReal const * const AMREX_RESTRICT w_1 = + soa_1.GetRealData(PIdx::w).data(); + index_type* indices_1 = bins_1.permutationPtr(); + index_type const* cell_offsets_1 = bins_1.offsetsPtr(); + Real q1 = species_1->getCharge(); + Real m1 = species_1->getMass(); + // - Species 2 + auto& soa_2 = ptile_2.GetStructOfArrays(); + Real* ux_2 = soa_2.GetRealData(PIdx::ux).data(); + Real* uy_2 = soa_2.GetRealData(PIdx::uy).data(); + Real* uz_2 = soa_2.GetRealData(PIdx::uz).data(); + Real* w_2 = soa_2.GetRealData(PIdx::w).data(); + index_type* indices_2 = bins_2.permutationPtr(); + index_type const* cell_offsets_2 = bins_2.offsetsPtr(); + Real q2 = species_2->getCharge(); + Real m2 = species_2->getMass(); + + const Real dt = WarpX::GetInstance().getdt(lev); + Geometry const& geom = WarpX::GetInstance().Geom(lev); + const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); + + // Loop over cells + amrex::ParallelFor( n_cells, + [=] AMREX_GPU_DEVICE (int i_cell) noexcept + { + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + index_type const cell_start_1 = cell_offsets_1[i_cell]; + index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; + // Same for species 2 + index_type const cell_start_2 = cell_offsets_2[i_cell]; + index_type const cell_stop_2 = cell_offsets_2[i_cell+1]; + + // ux from species1 can be accessed like this: + // ux_1[ indices_1[i] ], where i is between + // cell_start_1 (inclusive) and cell_start_2 (exclusive) + + // Do not collide if one species is missing in the cell + if ( cell_stop_1 - cell_start_1 >= 1 && + cell_stop_2 - cell_start_2 >= 1 ) + { + // shuffle + ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1); + ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2); + + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start_1, cell_stop_1, cell_start_2, cell_stop_2, + indices_1, indices_2, + ux_1, uy_1, uz_1, ux_2, uy_2, uz_2, w_1, w_2, + q1, q2, m1, m2, Real(-1.0), Real(-1.0), + dt, CoulombLog, dV ); + } + } + ); + } // end if ( isSameSpecies) + +} diff --git a/Source/Particles/Collision/ComputeTemperature.H b/Source/Particles/Collision/ComputeTemperature.H new file mode 100644 index 000000000..3cc96fb52 --- /dev/null +++ b/Source/Particles/Collision/ComputeTemperature.H @@ -0,0 +1,39 @@ +#ifndef WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ +#define WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ + +#include <WarpXConst.H> + +template <typename T_index, typename T_R> +T_R ComputeTemperature ( + T_index const Is, T_index const Ie, T_index const *I, + T_R const *ux, T_R const *uy, T_R const *uz, T_R const m ) +{ + + T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); + + int N = Ie - Is; + if ( N == 0 ) { return T_R(0.0); } + + T_R vx = T_R(0.0); T_R vy = T_R(0.0); + T_R vz = T_R(0.0); T_R vs = T_R(0.0); + T_R gm = T_R(0.0); T_R us = T_R(0.0); + + for (int i = Is; i < Ie; ++i) + { + us = ( ux[ I[i] ] * ux[ I[i] ] + + uy[ I[i] ] * uy[ I[i] ] + + uz[ I[i] ] * uz[ I[i] ] ); + gm = std::sqrt( T_R(1.0) + us*inv_c2 ); + vx += ux[ I[i] ] / gm; + vy += uy[ I[i] ] / gm; + vz += uz[ I[i] ] / gm; + vs += us / gm / gm; + } + + vx = vx / N; vy = vy / N; + vz = vz / N; vs = vs / N; + + return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz)); +} + +#endif // WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ diff --git a/Source/Particles/Collision/ElasticCollisionPerez.H b/Source/Particles/Collision/ElasticCollisionPerez.H new file mode 100644 index 000000000..8e16d95cc --- /dev/null +++ b/Source/Particles/Collision/ElasticCollisionPerez.H @@ -0,0 +1,105 @@ +#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ +#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ + +#include "UpdateMomentumPerezElastic.H" +#include "ComputeTemperature.H" +#include <WarpXConst.H> +#include <AMReX_Random.H> + +/** \brief Prepare information for and call + * UpdateMomentumPerezElastic(). + * @param[in] I1s,I2s is the start index for I1,I2 (inclusive). + * @param[in] I1e,I2e is the start index for I1,I2 (exclusive). + * @param[in] I1 and I2 are the index arrays. + * @param[in,out] u1 and u2 are the velocity arrays (u=v*gamma), + * they could be either different or the same, + * their lengths are not needed, + * @param[in] I1 and I2 determine all elements that will be used. + * @param[in] w1 and w2 are arrays of weights. + * @param[in] q1 and q2 are charges. m1 and m2 are masses. + * @param[in] T1 and T2 are temperatures (Joule) + * and will be used if greater than zero, + * otherwise will be computed. + * @param[in] dt is the time step length between two collision calls. + * @param[in] L is the Coulomb log and will be used if greater than zero, + * otherwise will be computed. + * @param[in] dV is the volume of the corresponding cell. +*/ + +template <typename T_index, typename T_R> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ElasticCollisionPerez ( + T_index const I1s, T_index const I1e, + T_index const I2s, T_index const I2e, + T_index *I1, T_index *I2, + T_R *u1x, T_R *u1y, T_R *u1z, + T_R *u2x, T_R *u2y, T_R *u2z, + T_R const *w1, T_R const *w2, + T_R const q1, T_R const q2, + T_R const m1, T_R const m2, + T_R const T1, T_R const T2, + T_R const dt, T_R const L, T_R const dV) +{ + + T_R constexpr inv_c2 = T_R(1.0)/(PhysConst::c*PhysConst::c); + int NI1 = I1e - I1s; + int NI2 = I2e - I2s; + + // get local T1t and T2t + T_R T1t; T_R T2t; + if ( T1 <= T_R(0.0) && L <= T_R(0.0) ) + { + T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1); + } + else { T1t = T1; } + if ( T2 <= T_R(0.0) && L <= T_R(0.0) ) + { + T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2); + } + else { T2t = T2; } + + // local density + T_R n1 = T_R(0.0); + T_R n2 = T_R(0.0); + T_R n12 = T_R(0.0); + for (int i1=I1s; i1<I1e; ++i1) { n1 += w1[ I1[i1] ]; } + for (int i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; } + n1 = n1 / dV; n2 = n2 / dV; + { + int i1 = I1s; int i2 = I2s; + for (int k = 0; k < amrex::max(NI1,NI2); ++k) + { + n12 += amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] ); + ++i1; if ( i1 == I1e ) { i1 = I1s; } + ++i2; if ( i2 == I2e ) { i2 = I2s; } + } + n12 = n12 / dV; + } + + // compute Debye length lmdD + T_R lmdD; + lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + + n2*q2*q2/(T2t*PhysConst::ep0) ); + T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) * + amrex::max(n1,n2), T_R(-1.0/3.0) ); + lmdD = amrex::max(lmdD, rmin); + + // call UpdateMomentumPerezElastic() + { + int i1 = I1s; int i2 = I2s; + for (int k = 0; k < amrex::max(NI1,NI2); ++k) + { + UpdateMomentumPerezElastic( + u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], + u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], + n1, n2, n12, + q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ], + dt, L, lmdD); + ++i1; if ( i1 == I1e ) { i1 = I1s; } + ++i2; if ( i2 == I2e ) { i2 = I2s; } + } + } + +} + +#endif // WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package new file mode 100644 index 000000000..163508fb7 --- /dev/null +++ b/Source/Particles/Collision/Make.package @@ -0,0 +1,10 @@ +CEXE_headers += CollisionType.H +CEXE_headers += ElasticCollisionPerez.H +CEXE_headers += ShuffleFisherYates.H +CEXE_headers += UpdateMomentumPerezElastic.H +CEXE_headers += ComputeTemperature.H + +CEXE_sources += CollisionType.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision diff --git a/Source/Particles/Collision/ShuffleFisherYates.H b/Source/Particles/Collision/ShuffleFisherYates.H new file mode 100644 index 000000000..621e654d6 --- /dev/null +++ b/Source/Particles/Collision/ShuffleFisherYates.H @@ -0,0 +1,29 @@ +#ifndef WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ +#define WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ + +#include <AMReX_Random.H> + +/* \brief Shuffle array according to Fisher-Yates algorithm. + * Only shuffle the part between is <= i < ie, n = ie-is. + * T_index shall be + * amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type +*/ + +template <typename T_index> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ShuffleFisherYates (T_index *array, T_index const is, T_index const ie) +{ + int j; + T_index buf; + for (int i = ie-1; i >= is+1; --i) + { + // get random number j: is <= j <= i + j = amrex::Random_int(i-is+1) + is; + // swop the ith array element with the jth + buf = array[i]; + array[i] = array[j]; + array[j] = buf; + } +} + +#endif // WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H new file mode 100644 index 000000000..948e8b075 --- /dev/null +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -0,0 +1,252 @@ +#ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ +#define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ + +#include <WarpXConst.H> +#include <AMReX_Random.H> +#include <cmath> // isnan() isinf() +#include <limits> // numeric_limits<float>::min() + +/* \brief Update particle velocities according to + * F. Perez et al., Phys.Plasmas.19.083104 (2012), + * which is based on Nanbu's method, PhysRevE.55.4642 (1997). + * @param[in] LmdD is max(Debye length, minimal interparticle distance). + * @param[in] L is the Coulomb log. A fixed L will be used if L > 0, + * otherwise L will be calculated based on the algorithm. + * To see if there are nan or inf updated velocities, + * compile with USE_ASSERTION=TRUE. +*/ + +template <typename T_R> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumPerezElastic ( + T_R& u1x, T_R& u1y, T_R& u1z, T_R& u2x, T_R& u2y, T_R& u2z, + T_R const n1, T_R const n2, T_R const n12, + T_R const q1, T_R const m1, T_R const w1, + T_R const q2, T_R const m2, T_R const w2, + T_R const dt, T_R const L, T_R const lmdD) +{ + + // If g = u1 - u2 = 0, do not collide. + if ( std::abs(u1x-u2x) < std::numeric_limits<T_R>::min() && + std::abs(u1y-u2y) < std::numeric_limits<T_R>::min() && + std::abs(u1z-u2z) < std::numeric_limits<T_R>::min() ) + { return; } + + T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); + + // Compute Lorentz factor gamma + T_R const g1 = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 ); + T_R const g2 = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 ); + + // Compute momenta + T_R const p1x = u1x * m1; + T_R const p1y = u1y * m1; + T_R const p1z = u1z * m1; + T_R const p2x = u2x * m2; + T_R const p2y = u2y * m2; + T_R const p2z = u2z * m2; + + // Compute center-of-mass (COM) velocity and gamma + T_R const mass_g = m1 * g1 + m2 * g2; + T_R const vcx = (p1x+p2x) / mass_g; + T_R const vcy = (p1y+p2y) / mass_g; + T_R const vcz = (p1z+p2z) / mass_g; + T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz; + T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2 ); + + // Compute vc dot v1 and v2 + T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1; + T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2; + + // Compute p1 star + T_R p1sx; + T_R p1sy; + T_R p1sz; + if ( vcms > std::numeric_limits<T_R>::min() ) + { + T_R const lorentz_tansform_factor = + ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1; + p1sx = p1x + vcx*lorentz_tansform_factor; + p1sy = p1y + vcy*lorentz_tansform_factor; + p1sz = p1z + vcz*lorentz_tansform_factor; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1sx = p1x; + p1sy = p1y; + p1sz = p1z; + } + T_R const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz ); + + // Compute gamma star + T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1; + T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2; + + // Compute the Coulomb log lnLmd + T_R lnLmd; + if ( L > T_R(0.0) ) { lnLmd = L; } + else + { + // Compute b0 + T_R const b0 = std::abs(q1*q2) * inv_c2 / + (T_R(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g * + ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) ); + + // Compute the minimal impact parameter + T_R bmin = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,b0); + + // Compute the Coulomb log lnLmd + lnLmd = amrex::max( T_R(2.0), + T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) ); + } + + // Compute s + T_R s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / + ( T_R(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * + m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * + std::pow(m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_R(1.0), 2.0); + + // Compute s' + T_R const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); + T_R const sp = std::pow(T_R(4.0)*MathConst::pi/T_R(3.0),T_R(1.0/3.0)) * + n1*n2/n12 * dt * vrel * (m1+m2) / + amrex::max( m1*std::pow(n1,T_R(2.0/3.0)), + m2*std::pow(n2,T_R(2.0/3.0)) ); + + // Determine s + s = amrex::min(s,sp); + + // Get random numbers + T_R r = amrex::Random(); + + // Compute scattering angle + T_R cosXs; + T_R sinXs; + if ( s <= T_R(0.1) ) + { + while ( true ) + { + cosXs = T_R(1.0) + s * std::log(r); + // Avoid the bug when r is too small such that cosXs < -1 + if ( cosXs >= T_R(-1.0) ) { break; } + r = amrex::Random(); + } + } + else if ( s > T_R(0.1) && s <= T_R(3.0) ) + { + T_R const Ainv = 0.0056958 + 0.9560202*s - 0.508139*s*s + + 0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s; + cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) + + T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) ); + } + else if ( s > T_R(3.0) && s <= T_R(6.0) ) + { + T_R const A = T_R(3.0) * std::exp(-s); + cosXs = T_R(1.0)/A * std::log( std::exp(-A) + + T_R(2.0) * r * std::sinh(A) ); + } + else + { + cosXs = T_R(2.0) * r - T_R(1.0); + } + sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs); + + // Get random azimuthal angle + T_R const phis = amrex::Random() * T_R(2.0) * MathConst::pi; + T_R const cosphis = std::cos(phis); + T_R const sinphis = std::sin(phis); + + // Compute post-collision momenta pfs in COM + T_R p1fsx; + T_R p1fsy; + T_R p1fsz; + // p1sp is the p1s perpendicular + T_R p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy ); + // Make sure p1sp is not almost zero + if ( p1sp > std::numeric_limits<T_R>::min() ) + { + p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis + + ( p1sy*p1sm/p1sp ) * sinXs*sinphis + + ( p1sx ) * cosXs; + p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis + + (-p1sx*p1sm/p1sp ) * sinXs*sinphis + + ( p1sy ) * cosXs; + p1fsz = (-p1sp ) * sinXs*cosphis + + ( T_R(0.0) ) * sinXs*sinphis + + ( p1sz ) * cosXs; + // Note a negative sign is different from + // Eq. (12) in Perez's paper, + // but they are the same due to the random nature of phis. + } + else + { + // If the previous p1sp is almost zero + // x->y y->z z->x + // This set is equivalent to the one in Nanbu's paper + p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz ); + p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis + + ( p1sz*p1sm/p1sp ) * sinXs*sinphis + + ( p1sy ) * cosXs; + p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis + + (-p1sy*p1sm/p1sp ) * sinXs*sinphis + + ( p1sz ) * cosXs; + p1fsx = (-p1sp ) * sinXs*cosphis + + ( T_R(0.0) ) * sinXs*sinphis + + ( p1sx ) * cosXs; + } + + T_R const p2fsx = -p1fsx; + T_R const p2fsy = -p1fsy; + T_R const p2fsz = -p1fsz; + + // Transform from COM to lab frame + T_R p1fx; T_R p2fx; + T_R p1fy; T_R p2fy; + T_R p1fz; T_R p2fz; + if ( vcms > std::numeric_limits<T_R>::min() ) + { + T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; + T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; + T_R const factor = (gc-T_R(1.0))/vcms; + T_R const factor1 = factor*vcDp1fs + m1*g1s*gc; + T_R const factor2 = factor*vcDp2fs + m2*g2s*gc; + p1fx = p1fsx + vcx * factor1; + p1fy = p1fsy + vcy * factor1; + p1fz = p1fsz + vcz * factor1; + p2fx = p2fsx + vcx * factor2; + p2fy = p2fsy + vcy * factor2; + p2fz = p2fsz + vcz * factor2; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1fx = p1fsx; + p1fy = p1fsy; + p1fz = p1fsz; + p2fx = p2fsx; + p2fy = p2fsy; + p2fz = p2fsz; + } + + // Rejection method + r = amrex::Random(); + if ( w2 > r*amrex::max(w1, w2) ) + { + u1x = p1fx / m1; + u1y = p1fy / m1; + u1z = p1fz / m1; + AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z)); + AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z)); + } + r = amrex::Random(); + if ( w1 > r*amrex::max(w1, w2) ) + { + u2x = p2fx / m2; + u2y = p2fy / m2; + u2z = p2fz / m2; + AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z)); + AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z)); + } + +} + +#endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 6d34fa0ef..8e336ba69 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -19,6 +19,7 @@ include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package include $(WARPX_HOME)/Source/Particles/Sorting/Make.package include $(WARPX_HOME)/Source/Particles/ParticleCreation/Make.package +include $(WARPX_HOME)/Source/Particles/Collision/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 72064ce8c..9db129b05 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -16,19 +16,27 @@ #include <QuantumSyncEngineWrapper.H> #endif +#include "CollisionType.H" + #include <memory> #include <map> #include <string> #include <algorithm> -// -// MultiParticleContainer holds multiple (nspecies or npsecies+1 when -// using laser) WarpXParticleContainer. WarpXParticleContainer is -// derived from amrex::ParticleContainer, and it is -// polymorphic. PhysicalParticleContainer and LaserParticleContainer are -// concrete class derived from WarpXParticlContainer. -// - +/** + * The class MultiParticleContainer holds multiple instances of the polymorphic + * class WarpXParticleContainer, stored in its member variable "allcontainers". + * The class WarpX typically has a single (pointer to an) instance of + * MultiParticleContainer. + * + * MultiParticleContainer typically has two types of functions: + * - Functions that loop over all instances of WarpXParticleContainer in + * allcontainers and calls the corresponding function (for instance, + * MultiParticleContainer::Evolve loops over all particles containers and + * calls the corresponding WarpXParticleContainer::Evolve function). + * - Functions that specifically handle multiple species (for instance + * ReadParameters or mapSpeciesProduct). + */ class MultiParticleContainer { @@ -142,6 +150,8 @@ public: void doFieldIonization (); + void doCoulombCollisions (); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -214,6 +224,10 @@ protected: std::vector<std::string> lasers_names; + std::vector<std::string> collision_names; + + amrex::Vector<std::unique_ptr<CollisionType> > allcollisions; + //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid std::vector<bool> m_deposit_on_main_grid; @@ -295,5 +309,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 bf1ffb74e..d84bc1afa 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -54,6 +54,14 @@ 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])); + } + } void @@ -120,6 +128,15 @@ MultiParticleContainer::ReadParameters () } } + // collision + ParmParse pc("collisions"); + pc.query("ncollisions", ncollisions); + BL_ASSERT(ncollisions >= 0); + if (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); @@ -609,6 +626,40 @@ MultiParticleContainer::doFieldIonization () } // pc_source } +void +MultiParticleContainer::doCoulombCollisions () +{ + BL_PROFILE("MPC::doCoulombCollisions"); + + for (int i = 0; i < ncollisions; ++i) + { + 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 if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ + + CollisionType::doCoulombCollisionsWithinTile + ( lev, mfi, species1, species2, + allcollisions[i]->m_isSameSpecies, + allcollisions[i]->m_CoulombLog ); + + } + } + } +} + #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/ParticleCreation/TransformParticle.H b/Source/Particles/ParticleCreation/TransformParticle.H index 1b71df723..c0158db78 100644 --- a/Source/Particles/ParticleCreation/TransformParticle.H +++ b/Source/Particles/ParticleCreation/TransformParticle.H @@ -21,7 +21,7 @@ void transformSourceParticle( /** * \brief small modifications on product particle * \param i index of particle - * \param particle particle struct + * \param v_particle pointer to vector of particles */ template < ElementaryProcessType ProcessT > AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 580f13f86..9b688cc59 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -4,6 +4,13 @@ #include <PhysicalParticleContainer.H> #include <AMReX_Vector.H> +/** + * Photon particles have no mass, they deposit no charge, and see specific QED + * effects. For these reasons, they are stored in the separate particle + * container PhotonParticleContainer, that inherts from + * PhysicalParticleContainer. The particle pusher and current deposition, in + * particular, are overriden in this container. + */ class PhotonParticleContainer : public PhysicalParticleContainer { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index e70b470b8..0192ffb37 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -15,6 +15,13 @@ #include <map> +/** + * PhysicalParticleContainer is the ParticleContainer class containing plasma + * particles (if a simulation has 2 plasma species, say "electrons" and + * "ions"), they will be two instances of PhysicalParticleContainer. + * + * PhysicalParticleContainer inherits from WarpXParticleContainer. + */ class PhysicalParticleContainer : public WarpXParticleContainer { @@ -68,6 +75,39 @@ public: int lev, int depos_lev); + /** + * \brief Evolve is the central function PhysicalParticleContainer that + * advances plasma particles for a time dt (typically one timestep). + * + * \param lev level on which particles are living + * \param Ex MultiFab from which field Ex is gathered + * \param Ey MultiFab from which field Ey is gathered + * \param Ez MultiFab from which field Ez is gathered + * \param Bx MultiFab from which field Bx is gathered + * \param By MultiFab from which field By is gathered + * \param Bz MultiFab from which field Bz is gathered + * \param jx MultiFab to which the particles' current jx is deposited + * \param jy MultiFab to which the particles' current jy is deposited + * \param jz MultiFab to which the particles' current jz is deposited + * \param cjx Same as jx (coarser, from lev-1), when using deposition buffers + * \param cjy Same as jy (coarser, from lev-1), when using deposition buffers + * \param cjz Same as jz (coarser, from lev-1), when using deposition buffers + * \param rho MultiFab to which the particles' charge is deposited + * \param crho Same as rho (coarser, from lev-1), when using deposition buffers + * \param cEx Same as Ex (coarser, from lev-1), when using gather buffers + * \param cEy Same as Ey (coarser, from lev-1), when using gather buffers + * \param cEz Same as Ez (coarser, from lev-1), when using gather buffers + * \param cBx Same as Bx (coarser, from lev-1), when using gather buffers + * \param cBy Same as By (coarser, from lev-1), when using gather buffers + * \param cBz Same as Bz (coarser, from lev-1), when using gather buffers + * \param t current physical time + * \param dt time step by which particles are advanced + * \param a_dt_type type of time step (used for sub-cycling) + * + * Evolve iterates over particle iterator (each box) and performs filtering, + * field gather, particle push and current deposition for all particles + * in the box. + */ virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -163,13 +203,39 @@ public: * \brief Apply NCI Godfrey filter to all components of E and B before gather * \param lev MR level * \param box box onto which the filter is applied - * \param exeli safeguard to avoid destructing arrays between ParIter iterations on GPU + * \param exeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ex + * \param eyeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ey + * \param ezeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ez + * \param bxeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Bx + * \param byeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field By + * \param bzeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Bz * \param filtered_Ex Array containing filtered value + * \param filtered_Ey Array containing filtered value + * \param filtered_Ez Array containing filtered value + * \param filtered_Bx Array containing filtered value + * \param filtered_By Array containing filtered value + * \param filtered_Bz Array containing filtered value * \param Ex Field array before filtering (not modified) - * \param ex_ptr pointer to the array used for field gather. + * \param Ey Field array before filtering (not modified) + * \param Ez Field array before filtering (not modified) + * \param Bx Field array before filtering (not modified) + * \param By Field array before filtering (not modified) + * \param Bz Field array before filtering (not modified) + * \param exfab pointer to the Ex field (modified) + * \param eyfab pointer to the Ey field (modified) + * \param ezfab pointer to the Ez field (modified) + * \param bxfab pointer to the Bx field (modified) + * \param byfab pointer to the By field (modified) + * \param bzfab pointer to the Bz field (modified) * * The NCI Godfrey filter is applied on Ex, the result is stored in filtered_Ex - * and the pointer is modified (before this function is called, it points to Ex + * and the pointer exfab is modified (before this function is called, it points to Ex * and after this function is called, it points to Ex_filtered) */ void applyNCIFilter ( diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 95e8ce4fe..94d9bc363 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -963,8 +963,12 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul void PhysicalParticleContainer::FieldGather (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) { const std::array<Real,3>& dx = WarpX::CellSize(lev); @@ -1078,11 +1082,11 @@ PhysicalParticleContainer::Evolve (int lev, for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const auto np = pti.numParticles(); - const auto lev = pti.GetLevel(); + const auto t_lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); tmp_particle_data.resize(finestLevel()+1); for (int i = 0; i < TmpIdx::nattribs; ++i) - tmp_particle_data[lev][index][i].resize(np); + tmp_particle_data[t_lev][index][i].resize(np); } } @@ -1427,15 +1431,19 @@ PhysicalParticleContainer::SplitParticles(int lev) { pti.GetPosition(xp, yp, zp); - // offset for split particles is computed as a function of cell size - // and number of particles per cell, so that a uniform distribution - // before splitting results in a uniform distribution after splitting const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); - amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0], - dx[1]/2._rt/ppc_nd[1], - dx[2]/2._rt/ppc_nd[2]}; - + amrex::Vector<Real> split_offset = {dx[0]/2._rt, + dx[1]/2._rt, + dx[2]/2._rt}; + if (ppc_nd[0] > 0){ + // offset for split particles is computed as a function of cell size + // and number of particles per cell, so that a uniform distribution + // before splitting results in a uniform distribution after splitting + split_offset[0] /= ppc_nd[0]; + split_offset[1] /= ppc_nd[1]; + split_offset[2] /= ppc_nd[2]; + } // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data @@ -2148,12 +2156,12 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, RealVector& Bxp, RealVector& Byp, RealVector& Bzp, - FArrayBox const * exfab, - FArrayBox const * eyfab, - FArrayBox const * ezfab, - FArrayBox const * bxfab, - FArrayBox const * byfab, - FArrayBox const * bzfab, + amrex::FArrayBox const * exfab, + amrex::FArrayBox const * eyfab, + amrex::FArrayBox const * ezfab, + amrex::FArrayBox const * bxfab, + amrex::FArrayBox const * byfab, + amrex::FArrayBox const * bzfab, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index a2473c5ad..fecb9c48e 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -4,6 +4,24 @@ #include <PhysicalParticleContainer.H> #include <AMReX_Vector.H> +/** + * When injecting a particle beam (typically for a plasma wakefield + * acceleration simulation), say propagating in the z direction, it can + * be necessary to make particles propagate in a straight line up to a given + * location z=z0. This is of particular importance when running in a boosted + * frame, where the beam may evolve due to its space charge fields before + * entering the plasma, causing the actual injected beam, and hence the whole + * simulation result, to depend on the Lorentz factor of the boost. + * + * This feature is implemented in RigidInjectedParticleContainer: At each + * iteration, for each particle, if z<z0 the particle moves in a straight line, + * and if z>z0 the particle evolves as a regular PhysicalParticleContainer. + * + * Note: This option is also useful to build self-consistent space charge + * fields for the particle beam. + * + * RigidInjectedParticleContainer derives from PhysicalParticleContainer. + */ class RigidInjectedParticleContainer : public PhysicalParticleContainer { diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 924c58382..fa5c586da 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -101,8 +101,32 @@ public: } }; +// Forward-declaration needed by WarpXParticleContainer below class MultiParticleContainer; +/** + * WarpXParticleContainer is the base polymorphic class from which all concrete + * particle container classes (that store a collection of particles) derive. Derived + * classes can be used for plasma particles, photon particles, or non-physical + * particles (e.g., for the laser antenna). + * It derives from amrex::ParticleContainer<0,0,PIdx::nattribs>, where the + * template arguments stand for the number of int and amrex::Real SoA and AoS + * data in amrex::Particle. + * - AoS amrex::Real: x, y, z (default), 0 additional (first template + * parameter) + * - AoS int: id, cpu (default), 0 additional (second template parameter) + * - SoA amrex::Real: PIdx::nattribs (third template parameter), see PIdx for + * the list. + * + * WarpXParticleContainer contains the main functions for initialization, + * interaction with the grid (field gather and current deposition) and particle + * push. + * + * Note: many functions are pure virtual (meaning they MUST be defined in + * derived classes, e.g., Evolve) or empty function (meaning they + * do not do anything, e.g., FieldGather, meant to be overriden by derived + * function) or actual functions (e.g. CurrentDeposition). + */ class WarpXParticleContainer : public amrex::ParticleContainer<0,0,PIdx::nattribs> { @@ -140,6 +164,11 @@ public: amrex::Real t, amrex::Real dt) = 0; #endif // WARPX_DO_ELECTROSTATIC + /** + * Evolve is the central WarpXParticleContainer function that advances + * particles for a time dt (typically one timestep). It is a pure virtual + * function for flexibility. + */ virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, @@ -288,6 +317,11 @@ public: std::map<std::string, int> getParticleComps () { return particle_comps;} std::map<std::string, int> getParticleiComps () { return particle_icomps;} + //amrex::Real getCharge () {return charge;} + amrex::ParticleReal getCharge () const {return charge;} + //amrex::Real getMass () {return mass;} + amrex::ParticleReal getMass () const {return mass;} + protected: std::map<std::string, int> particle_comps; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9f02da338..e6c6d31bd 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -420,7 +420,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const int * const ion_lev, - MultiFab* rho, int icomp, + amrex::MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) { @@ -507,7 +507,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, +WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, bool local, bool reset, bool do_rz_volume_scaling) { @@ -760,7 +760,7 @@ WarpXParticleContainer::PushXES (Real dt) } void -WarpXParticleContainer::PushX (Real dt) +WarpXParticleContainer::PushX (amrex::Real dt) { const int nLevels = finestLevel(); for (int lev = 0; lev <= nLevels; ++lev) { @@ -769,7 +769,7 @@ WarpXParticleContainer::PushX (Real dt) } void -WarpXParticleContainer::PushX (int lev, Real dt) +WarpXParticleContainer::PushX (int lev, amrex::Real dt) { BL_PROFILE("WPC::PushX()"); @@ -846,8 +846,9 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, { p.m_idata.id = DoSplitParticleID; } - // For the moment, do not do anything if particles goes - // to lower level. + if (pld.m_lev == lev-1){ + // For the moment, do not do anything if particles goes + // to lower level. } } diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 31315f64f..54dfe1955 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -14,8 +14,9 @@ namespace *num_boxes = mf.local_size(); int shapesize = AMREX_SPACEDIM; if (mf.nComp() > 1) shapesize += 1; - *shapes = (int*) malloc(shapesize * (*num_boxes) * sizeof(int)); - amrex::Real** data = (amrex::Real**) malloc((*num_boxes) * sizeof(amrex::Real*)); + *shapes = static_cast<int*>(malloc(sizeof(int)*shapesize * (*num_boxes))); + auto data = + static_cast<amrex::Real**>(malloc((*num_boxes) * sizeof(amrex::Real*))); #ifdef _OPENMP #pragma omp parallel diff --git a/Source/QED/QedTableParserHelperFunctions.H b/Source/QED/QedTableParserHelperFunctions.H index 9f9f37017..528613727 100644 --- a/Source/QED/QedTableParserHelperFunctions.H +++ b/Source/QED/QedTableParserHelperFunctions.H @@ -15,7 +15,7 @@ namespace QedUtils{ * This function safely extracts an amrex::Vector<T> from raw binary data. * T must be a simple datatype (e.g. an int, a float, a double...). * - * @param[in] p_char a pointer to the binary stream + * @param[in] p_data a pointer to the binary stream * @param[in] how_many how many T should be read from stream * @param[in] p_last a pointer to the last element of the char* array * @return {a tuple containing @@ -43,7 +43,7 @@ namespace QedUtils{ * This function safely extracts a T from raw binary data. * T must be a simple datatype (e.g. an int, a float, a double...). * - * @param[in] p_char a pointer to the binary stream + * @param[in] p_data a pointer to the binary stream * @param[in] p_last a pointer to the last element of the char* array * @return {a tuple containing * 1) flag (which is false if p_last is exceeded) diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index 7734d9011..70436cb72 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -13,7 +13,7 @@ namespace MathConst // https://physics.nist.gov/cuu/Constants/index.html namespace PhysConst { - static constexpr amrex::Real c = 299792458.; + static constexpr amrex::Real c = 299'792'458.; static constexpr amrex::Real ep0 = 8.8541878128e-12; static constexpr amrex::Real mu0 = 1.25663706212e-06; static constexpr amrex::Real q_e = 1.602176634e-19; diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index ab28c5446..195e309cf 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -1,3 +1,6 @@ +#ifndef WARPX_UTILS_H_ +#define WARPX_UTILS_H_ + #include <AMReX_REAL.H> #include <AMReX_Vector.H> #include <AMReX_MultiFab.H> @@ -21,3 +24,5 @@ namespace WarpXUtilIO{ */ bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); } + +#endif //WARPX_UTILS_H_ diff --git a/Source/WarpX.H b/Source/WarpX.H index caffa5e25..fe8b2840d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -171,7 +171,7 @@ public: void ComputeDt (); // Compute max_step automatically for simulations in a boosted frame. - void computeMaxStepBoostAccelerator(amrex::Geometry geom); + void computeMaxStepBoostAccelerator(const amrex::Geometry& geom); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c1eda30ba..377d103d1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1116,9 +1116,9 @@ WarpX::Evolve (int numsteps) { } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array<const MultiFab*, 3>& B, - const std::array<Real,3>& dx) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array<const amrex::MultiFab*, 3>& B, + const std::array<amrex::Real,3>& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1150,9 +1150,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array<const MultiFab*, 3>& B, - const std::array<Real,3>& dx, int ngrow) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array<const amrex::MultiFab*, 3>& B, + const std::array<amrex::Real,3>& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1184,9 +1184,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array<const MultiFab*, 3>& E, - const std::array<Real,3>& dx) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array<const amrex::MultiFab*, 3>& E, + const std::array<amrex::Real,3>& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1218,9 +1218,9 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array<const MultiFab*, 3>& E, - const std::array<Real,3>& dx, int ngrow) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array<const amrex::MultiFab*, 3>& E, + const std::array<amrex::Real,3>& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; |