parse and apply thermal boundary

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;

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;
+}

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 <AMReX_AmrCore.H>
 
 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;

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)