diff --git a/Source/Particles/ParticleBoundaries.H b/Source/Particles/ParticleBoundaries.H index 78b090056b1..7e578f190d7 100644 --- a/Source/Particles/ParticleBoundaries.H +++ b/Source/Particles/ParticleBoundaries.H @@ -22,6 +22,8 @@ struct ParticleBoundaries void Set_reflect_all_velocities (bool flag); void SetAll (ParticleBoundaryType bc); + /** Sets thermal velocity in ParticleBoundariesData 'data.m_uth' to u_th */ + void SetThermalVelocity (amrex::Real u_th); void SetBoundsX (ParticleBoundaryType bc_lo, ParticleBoundaryType bc_hi); void SetBoundsY (ParticleBoundaryType bc_lo, ParticleBoundaryType bc_hi); @@ -31,6 +33,9 @@ struct ParticleBoundaries void BuildReflectionModelParsers (); + /** True if any of the particle boundary condition type is Thermal */ + bool isAnyParticleBoundaryThermal (); + // reflection models for absorbing boundaries std::string reflection_model_xlo_str = "0.0"; std::string reflection_model_xhi_str = "0.0"; @@ -54,6 +59,7 @@ struct ParticleBoundaries ParticleBoundaryType ymax_bc; ParticleBoundaryType zmin_bc; ParticleBoundaryType zmax_bc; + amrex::Real m_uth = 0.; amrex::ParserExecutor<1> reflection_model_xlo; amrex::ParserExecutor<1> reflection_model_xhi; diff --git a/Source/Particles/ParticleBoundaries.cpp b/Source/Particles/ParticleBoundaries.cpp index 5b51fa3fd25..93412b30a46 100644 --- a/Source/Particles/ParticleBoundaries.cpp +++ b/Source/Particles/ParticleBoundaries.cpp @@ -6,6 +6,7 @@ */ #include "ParticleBoundaries.H" +#include "WarpX.H" #include "Utils/Parser/ParserUtils.H" @@ -32,6 +33,12 @@ ParticleBoundaries::SetAll (ParticleBoundaryType bc) data.zmax_bc = bc; } +void +ParticleBoundaries::SetThermalVelocity (amrex::Real u_th) +{ + data.m_uth = u_th; +} + void ParticleBoundaries::SetBoundsX (ParticleBoundaryType bc_lo, ParticleBoundaryType bc_hi) { @@ -87,3 +94,13 @@ ParticleBoundaries::BuildReflectionModelParsers () utils::parser::makeParser(reflection_model_zhi_str, {"v"})); data.reflection_model_zhi = reflection_model_zhi_parser->compile<1>(); } + +bool +ParticleBoundaries::isAnyParticleBoundaryThermal () +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) {return true;} + if (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) {return true;} + } + return false; +} diff --git a/Source/Particles/ParticleBoundaries_K.H b/Source/Particles/ParticleBoundaries_K.H index bbec1f54e01..69d7dc25eb1 100644 --- a/Source/Particles/ParticleBoundaries_K.H +++ b/Source/Particles/ParticleBoundaries_K.H @@ -1,4 +1,4 @@ -/* Copyright 2021 David Grote +/* Copyright 2021 David Grote, Remi Lehe, Revathi Jambunathan * * This file is part of WarpX. * @@ -8,18 +8,19 @@ #define PARTICLEBOUNDARIES_K_H_ #include "ParticleBoundaries.H" +#include "Initialization/SampleGaussianDistribution.H" #include namespace ApplyParticleBoundaries { - + using namespace amrex::literals; /* \brief Applies the boundary condition on a specific axis * This is called by apply_boundaries. */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void apply_boundary (amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax, - bool& change_sign_ux, bool& particle_lost, + bool& change_sign_ux, bool& rethermalize_x, bool& particle_lost, ParticleBoundaryType xmin_bc, ParticleBoundaryType xmax_bc, amrex::Real refl_probability_xmin, amrex::Real refl_probability_xmax, amrex::RandomEngine const& engine ) @@ -42,6 +43,10 @@ namespace ApplyParticleBoundaries { x = 2*xmin - x; change_sign_ux = true; } + else if (xmin_bc == ParticleBoundaryType::Thermal) { + x = 2*xmin - x; + rethermalize_x = true; + } } else if (x > xmax) { if (xmax_bc == ParticleBoundaryType::Open) { @@ -61,16 +66,38 @@ namespace ApplyParticleBoundaries { x = 2*xmax - x; change_sign_ux = true; } + else if (xmax_bc == ParticleBoundaryType::Thermal) { + x = 2*xmax - x; + rethermalize_x = true; + } } } - /* \brief Applies absorbing or reflecting boundary condition to the input particles, along all axis. + /* \brief Thermalize particles that have been identified to cross the boundary. + * The normal component samples from a half-Maxwellian, + * and the two tangential components will sample from full Maxwellian disbutions + * with thermal velocity uth + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void thermalize_boundary_particle (amrex::ParticleReal& u_norm, amrex::ParticleReal& u_tang1, + amrex::ParticleReal& u_tang2, amrex::Real uth, + amrex::RandomEngine const& engine) + { + u_tang1 = (uth > 0._rt) ? PhysConst::c * amrex::RandomNormal(0._rt, uth, engine) : 0._rt; + u_tang2 = (uth > 0._rt) ? PhysConst::c * amrex::RandomNormal(0._rt, uth, engine) : 0._rt; + u_norm = -1._rt * PhysConst::c * generateGaussianFluxDist(0._rt, uth, engine) * ((u_norm>0)-(u_norm<0)); + } + + + /* \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 positrion of the particle + * is change appropriately. The sign of the ve * Note that periodic boundaries are handled in AMReX code. * * \param x, xmin, xmax: particle x position, location of x boundary @@ -101,24 +128,39 @@ namespace ApplyParticleBoundaries { bool change_sign_ux = false; bool change_sign_uy = false; bool change_sign_uz = false; + bool rethermalize_x = false; // stores if particle crosses x boundary and needs to be thermalized + bool rethermalize_y = false; // stores if particle crosses y boundary and needs to be thermalized + bool rethermalize_z = false; // stores if particle crosses z boundary and needs to be thermalized #ifndef WARPX_DIM_1D_Z - apply_boundary(x, xmin, xmax, change_sign_ux, particle_lost, + apply_boundary(x, xmin, xmax, change_sign_ux, rethermalize_x, particle_lost, boundaries.xmin_bc, boundaries.xmax_bc, boundaries.reflection_model_xlo(-ux), boundaries.reflection_model_xhi(ux), engine); #endif #ifdef WARPX_DIM_3D - apply_boundary(y, ymin, ymax, change_sign_uy, particle_lost, + apply_boundary(y, ymin, ymax, change_sign_uy, rethermalize_y, particle_lost, boundaries.ymin_bc, boundaries.ymax_bc, boundaries.reflection_model_ylo(-uy), boundaries.reflection_model_yhi(uy), engine); #endif - apply_boundary(z, zmin, zmax, change_sign_uz, particle_lost, + apply_boundary(z, zmin, zmax, change_sign_uz, rethermalize_z, particle_lost, boundaries.zmin_bc, boundaries.zmax_bc, boundaries.reflection_model_zlo(-uz), boundaries.reflection_model_zhi(uz), engine); +#ifndef WARPX_DIM_RZ + if (rethermalize_x) { + thermalize_boundary_particle(ux, uy, uz, boundaries.m_uth, engine); + } + if (rethermalize_y) { + thermalize_boundary_particle(uy, uz, ux, boundaries.m_uth, engine); + } + if (rethermalize_z) { + thermalize_boundary_particle(uz, ux, uy, boundaries.m_uth, engine); + } +#endif + if (boundaries.reflect_all_velocities && (change_sign_ux | change_sign_uy | change_sign_uz)) { change_sign_ux = true; change_sign_uy = true; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 160eac0d19c..b9b2d10156d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -444,6 +444,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp_boundary.query("reflect_all_velocities", flag); m_boundary_conditions.Set_reflect_all_velocities(flag); + // currently supports only isotropic thermal distribution + // same distribution is applied to all boundaries + if (m_boundary_conditions.isAnyParticleBoundaryThermal()) { + amrex::Real boundary_uth; + utils::parser::getWithParser(pp_species_name,"boundary_uth",boundary_uth); + amrex::Print() << " thermal vel : " << boundary_uth << "\n"; + m_boundary_conditions.SetThermalVelocity(boundary_uth); + } } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core)