diff --git a/src/Hipace.H b/src/Hipace.H index 4a38fb33ff..12f7728761 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -147,8 +147,9 @@ public: /** * \brief does Coulomb collisions between plasmas and beams + * \param[in] islice slice number */ - void doCoulombCollision (); + void doCoulombCollision (int islice); /** * \brief add external fields to the field grid diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 4c81cdf2fc..4608d9c24f 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -858,7 +858,7 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st m_multi_beam.shiftSlippedParticles(islice, m_3D_geom[0]); // collisions for plasmas and beams - doCoulombCollision(); + doCoulombCollision(islice); // get minimum beam uz after push m_adaptive_time_step.GatherMinUzSlice(m_multi_beam, false); @@ -1286,7 +1286,7 @@ Hipace::AddGridExternalFields (const int lev, const int islice) } void -Hipace::doCoulombCollision () +Hipace::doCoulombCollision (int islice) { // collisions for all particles calculated on level 0 @@ -1301,7 +1301,7 @@ Hipace::doCoulombCollision () // TODO: enable tiling - CoulombCollision::doBeamPlasmaCoulombCollision( + CoulombCollision::doBeamPlasmaCoulombCollision( islice, m_all_collisions[i].m_collide_every, lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_CoulombLog, m_background_density_SI); } else { @@ -1311,7 +1311,7 @@ Hipace::doCoulombCollision () // TODO: enable tiling - CoulombCollision::doPlasmaPlasmaCoulombCollision( + CoulombCollision::doPlasmaPlasmaCoulombCollision( islice, m_all_collisions[i].m_collide_every, lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_isSameSpecies, m_all_collisions[i].m_CoulombLog, m_background_density_SI); } diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H index 4d201e3a88..c2597c8a07 100644 --- a/src/particles/collisions/CoulombCollision.H +++ b/src/particles/collisions/CoulombCollision.H @@ -20,6 +20,7 @@ public: int m_nbeams {0}; bool m_isSameSpecies {false}; amrex::Real m_CoulombLog {-1.}; + int m_collide_every {1}; /** Read parameters from the input file */ void ReadParameters ( @@ -31,6 +32,8 @@ public: * \brief Perform Coulomb collisions of plasma species over longitudinal push by 1 cell. * Particles of both species are sorted per cell, paired, and collided pairwise. * + * \param[in] islice current slice number + * \param[in] collide_every interval between collision calculations * \param[in] lev MR level * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box) * \param[in] geom corresponding geometry object @@ -41,7 +44,7 @@ public: * Coulomb logarithm is deduced from the plasma temperature, measured in each cell. * \param[in] background_density_SI background plasma density (only needed for normalized units) **/ - static void doPlasmaPlasmaCoulombCollision ( + static void doPlasmaPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog, amrex::Real background_density_SI); @@ -50,6 +53,8 @@ public: * \brief Perform Coulomb collisions of a beam with a plasma species over a push by one beam time step * Particles of both species are sorted per cell, paired, and collided pairwise. * + * \param[in] islice current slice number + * \param[in] collide_every interval between collision operations * \param[in] lev MR level * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box) * \param[in] geom corresponding geometry object @@ -59,7 +64,7 @@ public: * Coulomb logarithm is deduced from the plasma temperature, measured in each cell. * \param[in] background_density_SI background plasma density (only needed for normalized units) **/ - static void doBeamPlasmaCoulombCollision ( + static void doBeamPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI); diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index 5d1466f73a..e29c7dd3b4 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -25,6 +25,8 @@ CoulombCollision::ReadParameters( // default Coulomb log is -1, if < 0 (e.g. not specified), will be computed automatically pp.query("CoulombLog", m_CoulombLog); + // how often the collision operator should be applied - every m_collide_every'th slice + pp.query("collide_every", m_collide_every); for (int i=0; i<(int) beam_species_names.size(); i++) { if (beam_species_names[i] == collision_species[0]) m_nbeams += 1; @@ -58,11 +60,13 @@ CoulombCollision::ReadParameters( } void -CoulombCollision::doPlasmaPlasmaCoulombCollision ( +CoulombCollision::doPlasmaPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog, amrex::Real background_density_SI) { + if (islice%collide_every != 0) {return;} + HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); AMREX_ALWAYS_ASSERT(lev == 0); @@ -104,10 +108,11 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); // static_cast to avoid precision problems in FP32 const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * - PhysConstSI::q_e*PhysConstSI::q_e / - (PhysConstSI::ep0*PhysConstSI::m_e)); - const amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + PhysConstSI::q_e*PhysConstSI::q_e / + (PhysConstSI::ep0*PhysConstSI::m_e)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; amrex::ParallelForRNG( n_cells, @@ -183,10 +188,11 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); // static_cast to avoid precision problems in FP32 const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * - PhysConstSI::q_e*PhysConstSI::q_e / - (PhysConstSI::ep0*PhysConstSI::m_e)); - const amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + PhysConstSI::q_e*PhysConstSI::q_e / + (PhysConstSI::ep0*PhysConstSI::m_e)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; // 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); @@ -211,7 +217,7 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( // Do not collide if one species is missing in the cell if ( cell_stop1 - cell_start1 < 1 || - cell_stop2 - cell_start2 < 1 ) return; + cell_stop2 - cell_start2 < 1 ) return; // shuffle ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine); ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine); @@ -234,11 +240,13 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( } void -CoulombCollision::doBeamPlasmaCoulombCollision ( +CoulombCollision::doBeamPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI) { + if (islice%collide_every != 0) {return;} + HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()"); AMREX_ALWAYS_ASSERT(lev == 0);