-
Notifications
You must be signed in to change notification settings - Fork 17
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implement binary Coulomb collisions between plasma species (#676)
* add header files copy-pasted from WarpX, with modifications for the QSA PIC loop Co-authored-by: Axel Huebl <axelhuebl@lbl.gov> Co-authored-by: Andrew Myers <atmyers@lbl.gov> Co-authored-by: Yinjian Zhao <yinjianzhao@lbl.gov> Co-authored-by: Remi Lehe <rlehe@lbl.gov> * add rest of the implementation * fix error in CMake file * fix doxygen doc * add PhysConstSI object * update doc * Also expose SI constants as static constexpr * Update docs/source/run/parameters.rst * add CI for collisions with SI * add benchmark * update script to reset benchmarks * fix CI * fix * fix few typo and compiler warnings * fix benchmark by @SeverinDiederichs * Apply suggestions from code review Co-authored-by: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com> * check user input * compiler warning * Apply suggestions from code review * update to the new psi * fix typo * fix other typo, now code runs, let's see CI Co-authored-by: Axel Huebl <axelhuebl@lbl.gov> Co-authored-by: Andrew Myers <atmyers@lbl.gov> Co-authored-by: Yinjian Zhao <yinjianzhao@lbl.gov> Co-authored-by: Remi Lehe <rlehe@lbl.gov> Co-authored-by: Tools <tools@desy.de> Co-authored-by: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>
- Loading branch information
1 parent
5a75c20
commit 95c4f2c
Showing
19 changed files
with
929 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,49 @@ | ||
/* Copyright 2019-2020 Andrew Myers, Yinjian Zhao | ||
* | ||
* This file is part of WarpX. | ||
* | ||
* License: BSD-3-Clause-LBNL | ||
*/ | ||
|
||
#ifndef HIPACE_COMPUTE_TEMPERATURE_H_ | ||
#define HIPACE_COMPUTE_TEMPERATURE_H_ | ||
|
||
template <typename T_index, typename T_R> | ||
AMREX_GPU_HOST_DEVICE | ||
T_R ComputeTemperature ( | ||
T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I, | ||
T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT psi, | ||
T_R const m ) | ||
{ | ||
using namespace amrex::literals; | ||
T_R constexpr c2i = 1._rt/(PhysConstSI::c * PhysConstSI::c); | ||
|
||
const 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 < (int) Ie; ++i) | ||
{ | ||
// particle's Lorentz factor | ||
gm = (1.0_rt + ux[I[i]]*ux[I[i]]*c2i + uy[I[i]]*uy[I[i]]*c2i | ||
+ psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] ); | ||
const amrex::Real uz = PhysConstSI::c * (gm - psi[I[i]]); | ||
us = ( ux[ I[i] ] * ux[ I[i] ] + | ||
uy[ I[i] ] * uy[ I[i] ] + | ||
uz * uz); | ||
vx += ux[ I[i] ] / gm; | ||
vy += uy[ I[i] ] / gm; | ||
vz += uz / 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 // HIPACE_COMPUTE_TEMPERATURE_H_ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,45 @@ | ||
#ifndef HIPACE_COULOMB_COLLISION_H_ | ||
#define HIPACE_COULOMB_COLLISION_H_ | ||
|
||
#include "PlasmaParticleContainer.H" | ||
|
||
#include <AMReX_DenseBins.H> | ||
#include <AMReX_REAL.H> | ||
#include <AMReX_ParmParse.H> | ||
|
||
/** | ||
* \brief This class handles Coulomb collisions between 2 plasma species (can be the same). | ||
*/ | ||
class CoulombCollision | ||
{ | ||
public: | ||
int m_species1_index {-1}; | ||
int m_species2_index {-1}; | ||
bool m_isSameSpecies {false}; | ||
amrex::Real m_CoulombLog {-1.}; | ||
|
||
/** Constructor */ | ||
CoulombCollision( | ||
const std::vector<std::string>& species_names, | ||
std::string const collision_name); | ||
|
||
/** | ||
* \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] lev MR level | ||
* \param[in] bx transverse box (plasma particles will be sorted per-cell on this box) | ||
* \param[in] geom corresponding geometry object | ||
* \param[in,out] species1 first plasma species | ||
* \param[in,out] species2 second plasma species | ||
* \param[in] is_same_species whether both species are the same (intra-species collisions) | ||
* \param[in] CoulombLog Value of the Coulomb logarithm used for the collisions. If <0, the | ||
* Coulomb logarithm is deduced from the plasma temperature, measured in each cell. | ||
**/ | ||
static void doCoulombCollision ( | ||
int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, | ||
PlasmaParticleContainer& species2, const bool is_same_species, const amrex::Real CoulombLog); | ||
|
||
}; | ||
|
||
#endif // HIPACE_COULOMB_COLLISION_H_ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,184 @@ | ||
#include "CoulombCollision.H" | ||
#include "ShuffleFisherYates.H" | ||
#include "TileSort.H" | ||
#include "ElasticCollisionPerez.H" | ||
#include "utils/HipaceProfilerWrapper.H" | ||
|
||
CoulombCollision::CoulombCollision( | ||
const std::vector<std::string>& species_names, | ||
std::string const collision_name) | ||
{ | ||
using namespace amrex::literals; | ||
|
||
// TODO: ionization level | ||
// TODO: Fix dt | ||
|
||
// 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 is -1, if < 0 (e.g. not specified), will be computed automatically | ||
pp.query("CoulombLog", m_CoulombLog); | ||
|
||
for (int i=0; i<(int) 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; | ||
} | ||
m_isSameSpecies = collision_species[0] == collision_species[1] ? true : false; | ||
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( | ||
m_species1_index >= 0 && m_species2_index >= 0, | ||
"<collision name>.species must contain exactly the name of 2 species in plasma.names" | ||
); | ||
} | ||
|
||
void | ||
CoulombCollision::doCoulombCollision ( | ||
int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, | ||
PlasmaParticleContainer& species2, const bool is_same_species, const amrex::Real CoulombLog) | ||
{ | ||
HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); | ||
AMREX_ALWAYS_ASSERT(lev == 0); | ||
using namespace amrex::literals; | ||
|
||
if ( is_same_species ) // species_1 == species_2 | ||
{ | ||
// Logically particles per-cell, and return indices of particles in each cell | ||
PlasmaBins bins1 = findParticlesInEachTile(lev, bx, 1, species1, geom); | ||
int const n_cells = bins1.numBins(); | ||
|
||
// Counter to check there is only 1 box | ||
int count = 0; | ||
for (PlasmaParticleIterator pti(species1, lev); pti.isValid(); ++pti) { | ||
|
||
// Get particles SoA and AoS data | ||
auto& soa1 = pti.GetStructOfArrays(); | ||
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data(); | ||
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data(); | ||
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data(); | ||
const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data(); | ||
PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); | ||
PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); | ||
amrex::Real q1 = species1.GetCharge(); | ||
amrex::Real m1 = species1.GetMass(); | ||
|
||
const amrex::Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); | ||
const amrex::Real dt = geom.CellSize(2)/PhysConstSI::c; | ||
|
||
amrex::ParallelForRNG( | ||
n_cells, | ||
[=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept | ||
{ | ||
// The particles from species1 that are in the cell `i_cell` are | ||
// given by the `indices_1[cell_start_1:cell_stop_1]` | ||
PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; | ||
PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; | ||
PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2; | ||
|
||
if ( cell_stop1 - cell_start1 <= 1 ) return; | ||
// Do not collide if there is only one particle in the cell | ||
// shuffle | ||
ShuffleFisherYates( | ||
indices1, cell_start1, cell_half1, engine ); | ||
|
||
// TODO: FIX DT | ||
// Call the function in order to perform collisions | ||
ElasticCollisionPerez( | ||
cell_start1, cell_half1, | ||
cell_half1, cell_stop1, | ||
indices1, indices1, | ||
ux1, uy1, psi1, ux1, uy1, psi1, w1, w1, | ||
q1, q1, m1, m1, -1.0_rt, -1.0_rt, | ||
dt, CoulombLog, dV, engine ); | ||
} | ||
); | ||
count++; | ||
} | ||
AMREX_ALWAYS_ASSERT(count == 1); | ||
|
||
} else { | ||
|
||
// Logically particles per-cell, and return indices of particles in each cell | ||
PlasmaBins bins1 = findParticlesInEachTile(lev, bx, 1, species1, geom); | ||
PlasmaBins bins2 = findParticlesInEachTile(lev, bx, 1, species2, geom); | ||
|
||
int const n_cells = bins1.numBins(); | ||
|
||
// Counter to check there is only 1 box | ||
int count = 0; | ||
for (PlasmaParticleIterator pti(species1, lev); pti.isValid(); ++pti) { | ||
|
||
// Get particles SoA and AoS data for species 1 | ||
auto& soa1 = pti.GetStructOfArrays(); | ||
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data(); | ||
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data(); | ||
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data(); | ||
const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data(); | ||
PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); | ||
PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); | ||
amrex::Real q1 = species1.GetCharge(); | ||
amrex::Real m1 = species1.GetMass(); | ||
|
||
// Get particles SoA and AoS data for species 2 | ||
auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); | ||
auto& soa2 = ptile2.GetStructOfArrays(); | ||
amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux).data(); | ||
amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy).data(); | ||
amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi).data(); | ||
const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data(); | ||
PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); | ||
PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); | ||
amrex::Real q2 = species2.GetCharge(); | ||
amrex::Real m2 = species2.GetMass(); | ||
|
||
const amrex::Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); | ||
const amrex::Real dt = geom.CellSize(2)/PhysConstSI::c; | ||
|
||
// 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); | ||
// Loop over cells, and collide the particles in each cell | ||
|
||
// Loop over cells | ||
amrex::ParallelForRNG( | ||
n_cells, | ||
[=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept | ||
{ | ||
// The particles from species1 that are in the cell `i_cell` are | ||
// given by the `indices_1[cell_start_1:cell_stop_1]` | ||
PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; | ||
PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; | ||
// Same for species 2 | ||
PlasmaBins::index_type const cell_start2 = offsets2[i_cell]; | ||
PlasmaBins::index_type const cell_stop2 = offsets2[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_stop1 - cell_start1 < 1 || | ||
cell_stop2 - cell_start2 < 1 ) return; | ||
// shuffle | ||
ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine); | ||
ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine); | ||
|
||
// TODO: FIX DT. | ||
// Call the function in order to perform collisions | ||
ElasticCollisionPerez( | ||
cell_start1, cell_stop1, cell_start2, cell_stop2, | ||
indices1, indices2, | ||
ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, | ||
q1, q2, m1, m2, -1.0_rt, -1.0_rt, | ||
dt, CoulombLog, dV, engine ); | ||
} | ||
); | ||
count++; | ||
} | ||
AMREX_ALWAYS_ASSERT(count == 1); | ||
} | ||
} |
Oops, something went wrong.