diff --git a/Source/Particles/ParticleBoundaries_K.H b/Source/Particles/ParticleBoundaries_K.H index 5dd340d9053..55729990d7f 100644 --- a/Source/Particles/ParticleBoundaries_K.H +++ b/Source/Particles/ParticleBoundaries_K.H @@ -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) { @@ -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(SEE_probability_xmin); + + SEE_probability = SEE_probability_xmin - static_cast(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; @@ -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(SEE_probability_xmax); + + SEE_probability = SEE_probability_xmax - static_cast(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; @@ -98,6 +128,7 @@ 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 @@ -105,6 +136,7 @@ namespace ApplyParticleBoundaries { * \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 @@ -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) { @@ -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); @@ -135,8 +169,10 @@ 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); @@ -144,8 +180,10 @@ namespace ApplyParticleBoundaries { #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); @@ -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(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(max_SEE_probability); + + SEE_probability = max_SEE_probability - static_cast(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 xp_new(SEE_num_particles_added, xpos); + amrex::Vector yp_new(SEE_num_particles_added, ypos); + amrex::Vector zp_new(SEE_num_particles_added, zpos); + amrex::Vector uxp_new(SEE_num_particles_added); + amrex::Vector uyp_new(SEE_num_particles_added); + amrex::Vector uzp_new(SEE_num_particles_added); + amrex::Vector 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_ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 8e91093d95b..1757f91689a 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -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 {