Skip to content

Commit

Permalink
Add framework for SEE
Browse files Browse the repository at this point in the history 
* Determine number of SEE electrons at each boundary
* Begin allocation of new particle parameters
  • Loading branch information
@budjensen
budjensen committed Feb 27, 2025
1 parent b15bda2 commit f4308a6
Show file tree
Hide file tree
Showing 2 changed files with 161 additions and 1 deletion.
147 changes: 146 additions & 1 deletion Source/Particles/ParticleBoundaries_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,10 @@ namespace ApplyParticleBoundaries {
void
apply_boundary (amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax,
bool& change_sign_ux, bool& rethermalize_x, bool& particle_lost,
bool& SEE_particle_added, int& SEE_num_particles_added,
ParticleBoundaryType xmin_bc, ParticleBoundaryType xmax_bc,
amrex::Real refl_probability_xmin, amrex::Real refl_probability_xmax,
amrex::Real SEE_probability_xmin, amrex::Real SEE_probability_xmax,
amrex::RandomEngine const& engine )
{
if (x < xmin) {
Expand All @@ -38,6 +40,20 @@ namespace ApplyParticleBoundaries {
x = 2*xmin - x;
change_sign_ux = true;
}
if (SEE_probability_xmin > 0) {
// Configure SEE boundary emission
amrex::Real SEE_probability = SEE_probability_xmin;
if (SEE_probability_xmin >= 1) {
SEE_particle_added = true;
SEE_num_particles_added = static_cast<int>(SEE_probability_xmin);

SEE_probability = SEE_probability_xmin - static_cast<amrex::Real>(SEE_num_particles_added);
}
if (amrex::Random(engine) < SEE_probability) {
SEE_particle_added = true;
SEE_num_particles_added += 1;
}
}
}
else if (xmin_bc == ParticleBoundaryType::Reflecting) {
x = 2*xmin - x;
Expand All @@ -61,6 +77,20 @@ namespace ApplyParticleBoundaries {
x = 2*xmax - x;
change_sign_ux = true;
}
if (SEE_probability_xmax > 0) {
// Configure SEE boundary emission
amrex::Real SEE_probability = SEE_probability_xmax;
if (SEE_probability_xmax >= 1) {
SEE_particle_added = true;
SEE_num_particles_added = static_cast<int>(SEE_probability_xmax);

SEE_probability = SEE_probability_xmax - static_cast<amrex::Real>(SEE_num_particles_added);
}
if (amrex::Random(engine) < SEE_probability) {
SEE_particle_added = true;
SEE_num_particles_added += 1;
}
}
}
else if (xmax_bc == ParticleBoundaryType::Reflecting) {
x = 2*xmax - x;
Expand Down Expand Up @@ -98,13 +128,15 @@ namespace ApplyParticleBoundaries {
* coefficient is zero by default.
* For thermal boundaries, the particle is first reflected and the position of the particle
* is changed appropriately.
* For SEE boundaries, the particle is removed and a new SEE particle is added.
* Note that periodic boundaries are handled in AMReX code.
*
* \param x, xmin, xmax: particle x position, location of x boundary
* \param y, ymin, ymax: particle y position, location of y boundary (3D only)
* \param z, zmin, zmax: particle z position, location of z boundary
* \param ux, uy, uz: particle momenta
* \param particle_lost: output, flags whether the particle was lost
* \param SEE_particle_added: output, flags whether a SEE particle was added
* \param boundaries: object with boundary condition settings
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -114,7 +146,7 @@ namespace ApplyParticleBoundaries {
[[maybe_unused]] amrex::ParticleReal& z,
amrex::XDim3 gridmin, amrex::XDim3 gridmax,
amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
bool& particle_lost,
bool& particle_lost, bool& SEE_particle_added, int& SEE_num_particles_added,
ParticleBoundaries::ParticleBoundariesData const& boundaries,
amrex::RandomEngine const& engine)
{
Expand All @@ -125,8 +157,10 @@ namespace ApplyParticleBoundaries {
#ifndef WARPX_DIM_1D_Z
bool rethermalize_x = false; // stores if particle crosses x boundary and needs to be thermalized
apply_boundary(x, gridmin.x, gridmax.x, change_sign_ux, rethermalize_x, particle_lost,
SEE_particle_added, SEE_num_particles_added,
boundaries.xmin_bc, boundaries.xmax_bc,
boundaries.reflection_model_xlo(-ux), boundaries.reflection_model_xhi(ux),
boundaries.SEE_probability_xlo, boundaries.SEE_probability_xhi,
engine);
if (rethermalize_x) {
thermalize_boundary_particle(ux, uy, uz, boundaries.m_uth, engine);
Expand All @@ -135,17 +169,21 @@ namespace ApplyParticleBoundaries {
#ifdef WARPX_DIM_3D
bool rethermalize_y = false; // stores if particle crosses y boundary and needs to be thermalized
apply_boundary(y, gridmin.y, gridmax.y, change_sign_uy, rethermalize_y, particle_lost,
SEE_particle_added, SEE_num_particles_added,
boundaries.ymin_bc, boundaries.ymax_bc,
boundaries.reflection_model_ylo(-uy), boundaries.reflection_model_yhi(uy),
boundaries.SEE_probability_ylo, boundaries.SEE_probability_yhi,
engine);
if (rethermalize_y) {
thermalize_boundary_particle(uy, uz, ux, boundaries.m_uth, engine);
}
#endif
bool rethermalize_z = false; // stores if particle crosses z boundary and needs to be thermalized
apply_boundary(z, gridmin.z, gridmax.z, change_sign_uz, rethermalize_z, particle_lost,
SEE_particle_added, SEE_num_particles_added,
boundaries.zmin_bc, boundaries.zmax_bc,
boundaries.reflection_model_zlo(-uz), boundaries.reflection_model_zhi(uz),
boundaries.SEE_probability_zlo, boundaries.SEE_probability_zhi,
engine);
if (rethermalize_z) {
thermalize_boundary_particle(uz, ux, uy, boundaries.m_uth, engine);
Expand Down Expand Up @@ -178,5 +216,112 @@ namespace ApplyParticleBoundaries {
if (change_sign_uz) { uz = -uz; }
}


/* \brief Applies absorbing, reflecting or thermal boundary condition to the input particles, along all axis.
* For reflecting boundaries, the position of the particle is changed appropriately and
* the sign of the velocity is changed (depending on the reflect_all_velocities flag).
* For absorbing, a flag is set whether the particle has been lost (it is up to the calling
* code to take appropriate action to remove any lost particles). Absorbing boundaries can
* be given a reflection coefficient for stochastic reflection of particles, this
* coefficient is zero by default.
* For thermal boundaries, the particle is first reflected and the position of the particle
* is changed appropriately.
* For SEE boundaries, the particle is removed and a new SEE particle is added.
* Note that periodic boundaries are handled in AMReX code.
*
* \param xpos, ypos, zpos: particle x, y, z position before boundary conditions were applied
* \param gridmin, gridmax: location of boundaries
* \param max_SEE_probability: maximum probability of adding a SEE particle at any boundary
* \param SEE_num_particles_added: number of SEE particles added
* \param engine: random number generator
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
do_SEE(amrex::ParticleReal& xpos, amrex::ParticleReal& ypos, amrex::ParticleReal& zpos,
amrex::XDim3 gridmin, amrex::XDim3 gridmax,
amrex::Real const& max_SEE_probability, amrex::Real const& v_SEE,
int& SEE_num_particles_added, amrex::RandomEngine const& engine)
{
// Make sure there are not more particles added (in case two boundaries were crossed by one particle)
if (SEE_num_particles_added > static_cast<int>(max_SEE_probability) + 1) {
SEE_num_particles_added = 0;
amrex::Real SEE_probability = max_SEE_probability;

// Redraw a random number to determine if SEE particles are added
if (max_SEE_probability >= 1) {
SEE_num_particles_added = static_cast<int>(max_SEE_probability);

SEE_probability = max_SEE_probability - static_cast<amrex::Real>(SEE_num_particles_added);
}
if (amrex::Random(engine) < SEE_probability) {
SEE_num_particles_added += 1;
}
}

// Determine normal vector pointing inward from the boundary that was hit
amrex::ParticleReal nx = 0.0;
amrex::ParticleReal ny = 0.0;
amrex::ParticleReal nz = 0.0;

// Determine which boundaries were crossed
#ifndef WARPX_DIM_1D_Z
if (xpos < gridmin.x) {
xpos = gridmin.x;
nx = 1.0;
}
else if (xpos > gridmax.x) {
xpos = gridmax.x;
nx = -1.0;
}
#endif
#ifdef WARPX_DIM_3D
if (ypos < gridmin.y) {
ypos = gridmin.y;
ny = 1.0;
}
else if (ypos > gridmax.y) {
ypos = gridmax.y;
ny = -1.0;
}
#endif
if (zpos < gridmin.z) {
zpos = gridmin.z;
nz = 1.0;
}
else if (zpos > gridmax.z) {
zpos = gridmax.z;
nz = -1.0;
}

// Do we need to do anything special to account for RZ simulations?

// Initialize arrays to hold new particle data
amrex::Vector<amrex::ParticleReal> xp_new(SEE_num_particles_added, xpos);
amrex::Vector<amrex::ParticleReal> yp_new(SEE_num_particles_added, ypos);
amrex::Vector<amrex::ParticleReal> zp_new(SEE_num_particles_added, zpos);
amrex::Vector<amrex::ParticleReal> uxp_new(SEE_num_particles_added);
amrex::Vector<amrex::ParticleReal> uyp_new(SEE_num_particles_added);
amrex::Vector<amrex::ParticleReal> uzp_new(SEE_num_particles_added);
amrex::Vector<amrex::ParticleReal> wp_new(SEE_num_particles_added, wp[i]); // Use same weight as parent

// Generate velocities pointing inward from the boundary
for (int n = 0; n < SEE_num_particles_added; n++) {
// Generate random angles in the hemisphere pointing inward
const amrex::ParticleReal cos_theta = amrex::Random(engine); // Range [0,1] for hemisphere
const amrex::ParticleReal sin_theta = std::sqrt(1.0 - cos_theta*cos_theta);
const amrex::ParticleReal phi = 2.0 * MathConst::pi * amrex::Random(engine);

// Calculate velocity components in local hemisphere coordinate system
const amrex::ParticleReal vr = v_SEE * sin_theta;
const amrex::ParticleReal vz_local = v_SEE * cos_theta;
const amrex::ParticleReal vx_local = vr * std::cos(phi);
const amrex::ParticleReal vy_local = vr * std::sin(phi);

// TODO: Align z velocity component to boundary normal and then confirm only one boundary was crossed. If
// multiple boundaries were crossed, the particle needs to ensure it is directed inwards.
}

// TODO: Add the SEE particle, maybe with WarpXParticleContainer::AddNParticles
}
}
#endif //WARPX_PARTICLEBOUNDARIES_K_H_
15 changes: 15 additions & 0 deletions Source/Particles/WarpXParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1648,11 +1648,26 @@ WarpXParticleContainer::ApplyBoundaryConditions (){
GetPosition.AsStored(i, x, y, z);
// Note that for RZ, (x, y, z) is actually (r, theta, z).

// Save the particle position before applying the boundary conditions
ParticleReal xpos = x;
ParticleReal ypos = y;
ParticleReal zpos = z;

bool particle_lost = false;
bool SEE_particle_added = false;
int SEE_num_particles_added = 0;
ApplyParticleBoundaries::apply_boundaries(x, y, z, gridmin, gridmax,
ux[i], uy[i], uz[i], particle_lost,
SEE_particle_added, SEE_num_particles_added,
boundary_conditions, engine);

if (SEE_particle_added) {
ApplyParticleBoundaries::do_SEE(xpos, ypos, zpos, gridmin, gridmax,
boundary_conditions.max_SEE_probability,
boundary_conditions.v_SEE,
SEE_num_particles_added, engine);
}

if (particle_lost) {
pidw.make_invalid();
} else {
Expand Down

0 comments on commit f4308a6

Please sign in to comment.