Skip to content

Commit

Permalink
update ImpactIonization.H for energy splitting
Browse files Browse the repository at this point in the history
  • Loading branch information
@k-wolfinger
k-wolfinger committed Sep 12, 2024
1 parent b8de621 commit fd62c4e
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 222 deletions.
2 changes: 0 additions & 2 deletions Python/pywarpx/picmi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1578,8 +1578,6 @@ def solver_initialize_inputs(self):

pywarpx.warpx.poisson_solver = self.method

pywarpx.warpx.poisson_solver = self.method


class GaussianLaser(picmistandard.PICMI_GaussianLaser):
def laser_initialize_inputs(self):
Expand Down
12 changes: 12 additions & 0 deletions Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,7 @@ BackgroundMCCCollision::BackgroundMCCCollision (std::string const& collision_nam
const std::string kw_energy = scattering_process + "_energy";
utils::parser::getWithParser(
pp_collision_name, kw_energy.c_str(), energy);
// std::cout << scattering_process << std::endl;
}

ScatteringProcess process(scattering_process, cross_section_file, energy);
Expand Down Expand Up @@ -371,6 +372,9 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile
amrex::ParallelForRNG(np,
[=] AMREX_GPU_HOST_DEVICE (long ip, amrex::RandomEngine const& engine)
{
// std::string s1 = std::to_string(total_collision_prob);
// std::cout << s1 << std::endl;

// determine if this particle should collide
if (amrex::Random(engine) > total_collision_prob) { return; }

Expand Down Expand Up @@ -417,9 +421,15 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile
// calculate normalized collision frequency
nu_i += n_a * sigma_E * v_coll / nu_max;

// std::cout << "looping" << std::endl;
// std::string s2 = std::to_string(nu_i);
// std::cout << s2 << std::endl;

// check if this collision should be performed
if (col_select > nu_i) { continue; }

// std::cout << "COLL" << std::endl;

// charge exchange is implemented as a simple swap of the projectile
// and target velocities which doesn't require any of the Lorentz
// transformations below; note that if the projectile and target
Expand Down Expand Up @@ -453,11 +463,13 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile
// transform to COM frame
ParticleUtils::doLorentzTransform(vx, vy, vz, uCOM_x, uCOM_y, uCOM_z);

// std::cout << (scattering_process.m_type == ScatteringProcessType::ELASTIC) << std::endl;
if ((scattering_process.m_type == ScatteringProcessType::ELASTIC)
|| (scattering_process.m_type == ScatteringProcessType::EXCITATION)) {
ParticleUtils::RandomizeVelocity(
vx, vy, vz, sqrt(vx*vx + vy*vy + vz*vz), engine
);
// std::cout << "elastic coll occurred" << std::endl;
}
else if (scattering_process.m_type == ScatteringProcessType::BACK) {
// elastic scattering with cos(chi) = -1 (i.e. 180 degrees)
Expand Down
243 changes: 23 additions & 220 deletions Source/Particles/Collision/BackgroundMCC/ImpactIonization.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,236 +204,40 @@ public:
// m_mass1 is the mass of the incident species
const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
ParticleUtils::getEnergy(u_coll2, m_mass1, E_coll);

ParticleUtils::getEnergy(u_coll2, m_mass1, E_coll);

// calculate unit vector of the incident species
const amrex::ParticleReal x_unit = ux / sqrt(u_coll2);
const amrex::ParticleReal y_unit = uy / sqrt(u_coll2);
const amrex::ParticleReal z_unit = uz / sqrt(u_coll2);

// the incident particle gets half of the available energy
// the ionized electron gets the other half of the available energy
// TODO: This should be changed to encode the correct physics
const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) / 2.0_prt * PhysConst::q_e);
// const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) / 2.0_prt * PhysConst::q_e);
const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) * PhysConst::q_e);
constexpr auto c2 = PhysConst::c * PhysConst::c;

// Corresponding momentum of the incident species
const auto mc2 = m_mass1 * c2;
const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m_mass1;
// isotropically scatter incident species
ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine);
// // isotropically scatter incident species
// ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine);
// Use same unit vector as before collision, essentially turning off scattering
ux = x_unit * up;
uy = y_unit * up;
uz = z_unit * up;

// Corresponding momentum of the electron species
constexpr auto mc2_e = PhysConst::m_e * c2;
const amrex::ParticleReal up_e = sqrt(E_out * (E_out + 2.0_prt*mc2_e) / c2) / PhysConst::m_e;
// isotropically scatter incident species
ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine);

// get velocities for the ion from a Maxwellian distribution
i_ux = ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine);
i_uy = ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine);
i_uz = ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine);
}

private:
amrex::ParticleReal m_energy_cost;
double m_mass1;
amrex::ParticleReal m_sqrt_kb_m;
amrex::ParserExecutor<4> m_T_a_func;
amrex::Real m_t;
};

// Personal functions

/**
* \brief Filter functor for ION impact ionization
*/
class IonImpactIonizationFilterFunc
{
public:

/**
* \brief Constructor of the IonImpactIonizationFilterFunc functor.
*
* This function sample a random number and compares it to the total
* collision probability to see if the given particle should be considered
* for an ionization event. If the particle passes this stage the collision
* cross-section is calculated given its energy and another random number
* is used to determine whether it actually collides.
* Note that the mass and energy quantities have to be stored as doubles to
* ensure accurate energy calculations (otherwise errors occur with single
* or mixed precision builds of WarpX).
*
* @param[in] mcc_process an ScatteringProcess object associated with the ionization
* @param[in] mass colliding particle's mass (could also assume electron)
* @param[in] total_collision_prob total probability for a collision to occur
* @param[in] nu_max maximum collision frequency
* @param[in] n_a_func ParserExecutor<4> function to get the background
density in m^-3 as a function of space and time
* @param[in] t the current simulation time
*/
IonImpactIonizationFilterFunc(
ScatteringProcess const& mcc_process,
double const mass1,
double const mass2,
amrex::ParticleReal const total_collision_prob,
amrex::ParticleReal const nu_max,
amrex::ParserExecutor<4> const& n_a_func,
amrex::Real t
) : m_mcc_process(mcc_process.executor()), m_mass(mass1),
m_total_collision_prob(total_collision_prob),
m_nu_max(nu_max), m_n_a_func(n_a_func), m_t(t) { }

/**
* \brief Functor call. This method determines if a given (electron) particle
* should undergo an ionization collision.
*
* @param[in] ptd particle tile data
* @param[in] i particle index
* @param[in] engine the random number state and factory
* @return true if a collision occurs, false otherwise
*/
template <typename PData>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator() (
const PData& ptd, int const i, amrex::RandomEngine const& engine
) const noexcept
{
using namespace amrex;
using std::sqrt;

// determine if this particle should collide
if (Random(engine) > m_total_collision_prob) { return false; }

// get references to the particle to get its position
const auto& p = ptd.getSuperParticle(i);
ParticleReal x, y, z;
double E_coll;
get_particle_position(p, x, y, z);

// calculate neutral density at particle location
const ParticleReal n_a = m_n_a_func(x, y, z, m_t);

// get the particle velocity
const ParticleReal ux = ptd.m_rdata[PIdx::ux][i];
const ParticleReal uy = ptd.m_rdata[PIdx::uy][i];
const ParticleReal uz = ptd.m_rdata[PIdx::uz][i];

// calculate kinetic energy
const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
ParticleUtils::getEnergy(u_coll2, m_mass, E_coll);

// get collision cross-section
const ParticleReal sigma_E = m_mcc_process.getCrossSection(static_cast<amrex::ParticleReal>(E_coll));

// calculate normalized collision frequency
const ParticleReal nu_i = n_a * sigma_E * sqrt(u_coll2) / m_nu_max;

// check if this collision should be performed
return (Random(engine) <= nu_i);
}

private:
ScatteringProcess::Executor m_mcc_process;
double m_mass;
amrex::ParticleReal m_total_collision_prob = 0;
amrex::ParticleReal m_nu_max;
amrex::ParserExecutor<4> m_n_a_func;
amrex::Real m_t;
};


/**
* \brief Transform functor for ION impact ionization
*/
class IonImpactIonizationTransformFunc
{
public:

/**
* \brief Constructor of the IonImpactIonizationTransformFunc functor.
*
* The transform is responsible for appropriately decreasing the kinetic
* energy of the colliding particle and assigning appropriate velocities
* to the two newly created particles. To this end the energy cost of
* ionization is passed to the constructor as well as the mass of the
* colliding species and the standard deviation of the ion velocity
* (normalized temperature).
* Note that the mass and energy quantities have to be stored as doubles to
* ensure accurate energy calculations (otherwise errors occur with single
* or mixed precision builds of WarpX).
*
* @param[in] energy_cost energy cost of ionization
* @param[in] mass1 mass of the colliding species
* @param[in] sqrt_kb_m value of sqrt(kB/m), where kB is Boltzmann's constant
and m is the background neutral mass
* @param[in] T_a_func ParserExecutor<4> function to get the background
temperature in Kelvin as a function of space and time
* @param[in] t the current simulation time
*/
IonImpactIonizationTransformFunc(
amrex::ParticleReal energy_cost, double mass1, amrex::ParticleReal sqrt_kb_m,
amrex::ParserExecutor<4> const& T_a_func, amrex::Real t
) : m_energy_cost(energy_cost), m_mass1(mass1),
m_sqrt_kb_m(sqrt_kb_m), m_T_a_func(T_a_func), m_t(t) { }

/**
* \brief Functor call. It determines the properties of the generated pair
* and decreases the kinetic energy of the colliding particle. Inputs are
* standard from the FilterCopyTransfrom::filterCopyTransformParticles
* function.
*
* @param[in,out] dst1 target species 1 (electrons)
* @param[in,out] dst2 target species 2 (ions)
* @param[in] src source species (electrons)
* @param[in] i_src particle index of the source species
* @param[in] i_dst1 particle index of target species 1
* @param[in] i_dst2 particle index of target species 2
* @param[in] engine random number generator engine
*/
template <typename DstData, typename SrcData>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void operator() (DstData& dst1, DstData& dst2, SrcData& src,
int const i_src, int const i_dst1, int const i_dst2,
amrex::RandomEngine const& engine) const noexcept
{
using namespace amrex;
using std::sqrt;

// get references to the particle to get its position
const auto& p = src.getSuperParticle(i_src);
ParticleReal x, y, z;
double E_coll;
get_particle_position(p, x, y, z);

// calculate standard deviation in neutral velocity distribution using
// the local temperature
const ParticleReal ion_vel_std = m_sqrt_kb_m * std::sqrt(m_T_a_func(x, y, z, m_t));

// get references to the original particle's velocity
auto& ux = src.m_rdata[PIdx::ux][i_src];
auto& uy = src.m_rdata[PIdx::uy][i_src];
auto& uz = src.m_rdata[PIdx::uz][i_src];

// get references to the new particles' velocities
auto& e_ux = dst1.m_rdata[PIdx::ux][i_dst1];
auto& e_uy = dst1.m_rdata[PIdx::uy][i_dst1];
auto& e_uz = dst1.m_rdata[PIdx::uz][i_dst1];
auto& i_ux = dst2.m_rdata[PIdx::ux][i_dst2];
auto& i_uy = dst2.m_rdata[PIdx::uy][i_dst2];
auto& i_uz = dst2.m_rdata[PIdx::uz][i_dst2];

// calculate kinetic energy
const ParticleReal u_coll2 = ux*ux + uy*uy + uz*uz;
ParticleUtils::getEnergy(u_coll2, m_mass1, E_coll);

// each electron gets half the energy (could change this later)
const auto E_out = static_cast<amrex::ParticleReal>((E_coll - m_energy_cost) / 2.0_prt * PhysConst::q_e);

// precalculate often used value
constexpr auto c2 = PhysConst::c * PhysConst::c;
const auto mc2 = m_mass1*c2;

const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m_mass1;

// isotropically scatter electrons
ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine);
ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine);

const amrex::ParticleReal up_e = 0.0; // sqrt(E_out * (E_out + 2.0_prt*mc2_e) / c2) / PhysConst::m_e;
// // isotropically scatter incident species
// ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up_e, engine);
// Use same unit vector as before collision, essentially turning off scattering
e_ux = x_unit * up_e;
e_uy = y_unit * up_e;
e_uz = z_unit * up_e;

// get velocities for the ion from a Maxwellian distribution
i_ux = ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine);
i_uy = ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine);
Expand All @@ -448,5 +252,4 @@ private:
amrex::Real m_t;
};


#endif // WARPX_PARTICLES_COLLISION_IMPACT_IONIZATION_H_

0 comments on commit fd62c4e

Please sign in to comment.