Skip to content

Option to coarsen the calling of elastic collision operator #1403

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ibranajmudin wants to merge 4 commits into
base: development
Choose a base branch
from
Open

Option to coarsen the calling of elastic collision operator #1403

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
167f7e3
Added option collide_every, as an interval between applying the colli…
ibranajmudin Aug 26, 2026
4c07e75
Merge branch 'Hi-PACE:development' into collide_every
ibranajmudin Sep 10, 2026
2fba767
Finished addition of a collide_every option for coarsening collisions…
ibranajmudin Sep 10, 2026
0d8facb
Removed trailing whitespaces
ibranajmudin Sep 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/Hipace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand All @@ -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 {
Expand All @@ -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);
}
Expand Down
9 changes: 7 additions & 2 deletions src/particles/collisions/CoulombCollision.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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
Expand All @@ -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);
Expand Down
26 changes: 17 additions & 9 deletions src/particles/collisions/CoulombCollision.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

@AlexanderSinn AlexanderSinn Sep 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you change the name to "collision_period" to make it consistent with the output_period parameter name we use in the diagnostics? Additionally, can you add an entry in the documentation for it here https://github.com/Hi-PACE/hipace/blob/92d3da222f62481e9eacf3207b98aa22fea47b62/docs/source/run/parameters.rst?plain=1#L1457 

@AlexanderSinn AlexanderSinn Sep 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also I would be good to add an assert that it is >= 1.


for (int i=0; i<(int) beam_species_names.size(); i++) {
if (beam_species_names[i] == collision_species[0]) m_nbeams += 1;
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -104,10 +108,11 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2);
// static_cast<double> to avoid precision problems in FP32
const amrex::Real wp = std::sqrt(static_cast<double>(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,
Expand Down Expand Up @@ -183,10 +188,11 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2);
// static_cast<double> to avoid precision problems in FP32
const amrex::Real wp = std::sqrt(static_cast<double>(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);
Expand All @@ -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);
Expand All @@ -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;}

@AlexanderSinn AlexanderSinn Sep 22, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For Beam Plasma collisions this is a bit tricky. Beam particles usually stay on the same slice between time steps, so with this some won't get any collisions while others would get them every time step. Maybe we should just not allow the use of the parameter for this case. Or do something like (islice + step) % collision_period == 0 or step % collision_period == 0 . dt would also also need to be rescaled for that.


HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()");
AMREX_ALWAYS_ASSERT(lev == 0);

Expand Down
Loading