Add thermalBC for particles

Docs/source/usage/parameters.rst

@@ -417,6 +417,12 @@ Domain Boundary Conditions
     * ``Reflecting``: Particles leaving the boundary are reflected from the boundary back into the domain.
       When ``boundary.reflect_all_velocities`` is false, the sign of only the normal velocity is changed, otherwise the sign of all velocities are changed.
 
+    * ``Thermal``: Particles leaving the boundary are reflected from the boundary back into the domain
+      and their velocities are thermalized. The tangential velocity components are sampled from ``gaussian`` distribution
+      and the component normal to the boundary is sampled from ``gaussian flux`` distribution.
+      The standard deviation for these distributions should be provided for each species using
+      ``boundary.<species>.u_th``. The same standard deviation is used to sample all components.
+
 * ``boundary.reflect_all_velocities`` (`bool`) optional (default `false`)
     For a reflecting boundary condition, this flags whether the sign of only the normal velocity is changed or all velocities.
 

Examples/Tests/particle_thermal_boundary/analysis_2d.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+#
+# Copyright 2023-2024 Revathi Jambunathan
+#
+# This file is part of WarpX
+#
+# License: BSD-3-Clause-LBNL
+
+"""
+The script checks to see the growth in total field energy and total particle energy
+
+The input file in 2D initializes uniform plasma (electrons,ions) with
+thermal boundary condition. We do not expect the particle energy to increase
+beyond 2% in the time that it takes all particles to cross the domain boundary
+"""
+
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(1,'../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+FE_rdiag = './diags/reducedfiles/EF.txt'
+init_Fenergy = np.loadtxt(FE_rdiag)[1,2]
+final_Fenergy = np.loadtxt(FE_rdiag)[-1,2]
+assert(final_Fenergy/init_Fenergy < 40)
+assert(final_Fenergy < 5.e-5)
+
+PE_rdiag = './diags/reducedfiles/EN.txt'
+init_Penergy = np.loadtxt(PE_rdiag)[0,2]
+final_Penergy = np.loadtxt(PE_rdiag)[-1,2]
+assert( abs(final_Penergy - init_Penergy)/init_Penergy < 0.02)
+filename = sys.argv[1]
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Tests/particle_thermal_boundary/inputs_2d

@@ -0,0 +1,77 @@
+max_step = 2000
+
+# number of grid points
+amr.n_cell =   16 16
+
+# Maximum allowable size of each subdomain in the problem domain;
+#    this is used to decompose the domain for parallel calculations.
+amr.max_grid_size = 128
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo     = 0.e-6   0.e-6    # physical domain
+geometry.prob_hi     = 2.5e-7   2.5e-7
+
+# Boundary condition
+boundary.field_lo = pml pml
+boundary.field_hi = pml pml
+boundary.particle_lo = thermal thermal
+boundary.particle_hi = thermal thermal
+boundary.electrons.u_th = uth_e
+boundary.C.u_th = uth_C
+warpx.do_dive_cleaning = 0
+
+# Verbosity
+warpx.verbose = 1
+
+# Order of particle shape factors
+algo.particle_shape = 2
+
+# CFL
+warpx.cfl = 0.98
+
+# Density
+my_constants.n0 = 4e26
+my_constants.uth_e = 0.06256112470898544
+my_constants.uth_C = 0.00042148059678527106
+
+# Particles
+particles.species_names = electrons C
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.num_particles_per_cell_each_dim = 8 8
+electrons.profile = constant
+electrons.density = n0
+electrons.momentum_distribution_type = gaussian
+electrons.ux_th = uth_e
+electrons.uy_th = uth_e
+electrons.uz_th = uth_e
+
+C.charge = 6*q_e
+C.mass = 12*m_p
+C.injection_style = "NUniformPerCell"
+C.num_particles_per_cell_each_dim = 8 8
+C.profile = constant
+C.density = n0/6
+C.momentum_distribution_type = gaussian
+C.ux_th = uth_C
+C.uy_th = uth_C
+C.uz_th = uth_C
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 3000
+diag1.write_species = 0
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz rho divE
+diag1.diag_type = Full
+
+warpx.reduced_diags_names = EN EF
+EN.intervals = 10
+EN.type = ParticleEnergy
+EF.intervals = 10
+EF.type = FieldEnergy

Regression/Checksum/benchmarks_json/particle_thermal_boundary.json

@@ -0,0 +1,12 @@
+{
+  "lev=0": {
+    "Bx": 229.28366410112736,
+    "By": 242.40221902487008,
+    "Bz": 164.70543954019666,
+    "Ex": 789234700854.9321,
+    "Ey": 15892271613.51144,
+    "Ez": 693631132513.1528,
+    "divE": 5.441439446237374e+19,
+    "rho": 331367728.98455656
+  }
+}

Regression/WarpX-tests.ini

@@ -4720,3 +4720,20 @@ compareParticles = 1
 particleTypes = electrons
 outputFile = particle_boundary_interaction_plt
 analysisRoutine = Examples/Tests/particle_boundary_interaction/analysis.py
+
+[particle_thermal_boundary]
+buildDir = .
+inputFile = Examples/Tests/particle_thermal_boundary/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_OPENPMD=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/particle_thermal_boundary/analysis_2d.py

Source/BoundaryConditions/WarpXFieldBoundaries.cpp

@@ -106,7 +106,8 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d
 void WarpX::ApplyRhofieldBoundary (const int lev, MultiFab* rho,
                                    PatchType patch_type)
 {
-    if (PEC::isAnyParticleBoundaryReflecting() || PEC::isAnyBoundaryPEC()) {
+    if (PEC::isAnyParticleBoundaryReflecting() || WarpX::isAnyParticleBoundaryThermal() || PEC::isAnyBoundaryPEC())
+    {
         PEC::ApplyReflectiveBoundarytoRhofield(rho, lev, patch_type);
     }
 }
@@ -115,7 +116,8 @@ void WarpX::ApplyJfieldBoundary (const int lev, amrex::MultiFab* Jx,
                                  amrex::MultiFab* Jy, amrex::MultiFab* Jz,
                                  PatchType patch_type)
 {
-    if (PEC::isAnyParticleBoundaryReflecting() || PEC::isAnyBoundaryPEC()) {
+    if (PEC::isAnyParticleBoundaryReflecting() || WarpX::isAnyParticleBoundaryThermal() || PEC::isAnyBoundaryPEC())
+    {
         PEC::ApplyReflectiveBoundarytoJfield(Jx, Jy, Jz, lev, patch_type);
     }
 }

Source/BoundaryConditions/WarpX_PEC.cpp

@@ -258,8 +258,10 @@ PEC::ApplyReflectiveBoundarytoRhofield (amrex::MultiFab* rho, const int lev, Pat
     amrex::GpuArray<GpuArray<int,2>, AMREX_SPACEDIM> mirrorfac;
     for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
         is_reflective[idim][0] = ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                              || ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal)
                               || ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC);
         is_reflective[idim][1] = ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                              || ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal)
                               || ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC);
         if (!is_reflective[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); }
         if (!is_reflective[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); }
@@ -268,9 +270,11 @@ PEC::ApplyReflectiveBoundarytoRhofield (amrex::MultiFab* rho, const int lev, Pat
         // components of the current density
         is_tangent_to_bndy[idim] = true;
 
-        psign[idim][0] = (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+        psign[idim][0] = ( ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                        || ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) )
                          ? 1._rt : -1._rt;
-        psign[idim][1] = (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+        psign[idim][1] = ( (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                        || ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) )
                          ? 1._rt : -1._rt;
         mirrorfac[idim][0] = 2*domain_lo[idim] - (1 - rho_nodal[idim]);
         mirrorfac[idim][1] = 2*domain_hi[idim] + (1 - rho_nodal[idim]);
@@ -352,8 +356,10 @@ PEC::ApplyReflectiveBoundarytoJfield(amrex::MultiFab* Jx, amrex::MultiFab* Jy,
     amrex::GpuArray<GpuArray<GpuArray<int, 2>, AMREX_SPACEDIM>, 3> mirrorfac;
     for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
         is_reflective[idim][0] = ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                              || ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal)
                               || ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC);
         is_reflective[idim][1] = ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                              || ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal)
                               || ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC);
         if (!is_reflective[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); }
         if (!is_reflective[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); }
@@ -375,15 +381,19 @@ PEC::ApplyReflectiveBoundarytoJfield(amrex::MultiFab* Jx, amrex::MultiFab* Jy,
 #endif
 
             if (is_tangent_to_bndy[icomp][idim]){
-                psign[icomp][idim][0] = (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                psign[icomp][idim][0] = ( (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                                     || ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) )
                                         ? 1._rt : -1._rt;
-                psign[icomp][idim][1] = (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                psign[icomp][idim][1] = ( (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                                     || ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) )
                                         ? 1._rt : -1._rt;
             }
             else {
-                psign[icomp][idim][0] = (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                psign[icomp][idim][0] = ( (WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting)
+                                     || ( WarpX::particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) )
                                         ? -1._rt : 1._rt;
-                psign[icomp][idim][1] = (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                psign[icomp][idim][1] = ( (WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting)
+                                     || ( WarpX::particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) )
                                         ? -1._rt : 1._rt;
             }
         }

Source/Initialization/InjectorMomentum.H

@@ -12,6 +12,7 @@
 #include "GetVelocity.H"
 #include "TemperatureProperties.H"
 #include "VelocityProperties.H"
+#include "SampleGaussianFluxDistribution.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXConst.H"
 
@@ -90,75 +91,6 @@ private:
     amrex::Real m_ux_th, m_uy_th, m_uz_th;
 };
 
-namespace {
-    /** Return u sampled according to the probability distribution:
-      * p(u) \propto u \exp(-(u-u_m)^2/2u_th^2)
-      *
-      * @param u_m Central momentum
-      * @param u_th Momentum spread
-      * @param engine Object used to generate random numbers
-      */
-    [[nodiscard]]
-    AMREX_FORCE_INLINE
-    AMREX_GPU_HOST_DEVICE
-    amrex::Real
-    generateGaussianFluxDist( amrex::Real u_m, amrex::Real u_th, amrex::RandomEngine const& engine ) {
-
-        using namespace amrex::literals;
-
-        // Momentum to be returned at the end of this function
-        amrex::Real u = 0._rt;
-
-        const amrex::Real abs_u_m = std::abs(u_m);
-
-        if (u_th == 0._rt) {
-            u = u_m; // Trivial case ; avoids division by 0 in the rest of the code below
-        } else if (abs_u_m < 0.6*u_th) {
-            // Mean velocity magnitude is less than thermal velocity
-            // Use the distribution u*exp(-u**2*(1-abs(u_m)/u_th)/(2*u_th**2)) as an approximation
-            // and then use the rejection method to correct it
-            // ( stop rejecting with probability exp(-abs(u_m)/(2*u_th**3)*(u-sign(u_m)*u_th)**2) )
-            // Note that this is the method that is used in the common case u_m=0
-            const amrex::Real umsign = std::copysign(1._rt, u_m);
-            const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - abs_u_m/u_th );
-            const amrex::Real reject_prefactor = (abs_u_m/u_th)/(2._rt*u_th*u_th); // To save computation
-            bool reject = true;
-            while (reject) {
-                // Generates u according to u*exp(-u**2/(2*approx_u_th**2)),
-                // using the method of the inverse cumulative function
-                amrex::Real xrand = 1._rt - amrex::Random(engine); // ensures urand > 0
-                u = approx_u_th * std::sqrt(2._rt*std::log(1._rt/xrand));
-                // Rejection method
-                xrand = amrex::Random(engine);
-                if (xrand < std::exp(-reject_prefactor*(u - umsign*u_th)*(u - umsign*u_th))) { reject = false; }
-            }
-        } else {
-            // Mean velocity magnitude is greater than thermal velocity
-            // Use the distribution exp(-(u-u_m-u_th**2/abs(u_m))**2/(2*u_th**2)) as an approximation
-            // and then use the rejection method to correct it
-            // ( stop rejecting with probability (u/abs(u_m))*exp(1-(u/abs(u_m))) ; note
-            // that this number is always between 0 and 1 )
-            // Note that in the common case `u_m = 0`, this rejection method
-            // is not used, and the above rejection method is used instead.
-            bool reject = true;
-            const amrex::Real approx_u_m = u_m + u_th*u_th/abs_u_m;
-            const amrex::Real inv_um = 1._rt/abs_u_m; // To save computation
-            while (reject) {
-                // Approximate distribution: normal distribution, where we only retain positive u
-                u = -1._rt;
-                while (u < 0) {
-                    u = amrex::RandomNormal(approx_u_m, u_th, engine);
-                }
-                // Rejection method
-                const amrex::Real xrand = amrex::Random(engine);
-                if (xrand < u*inv_um* std::exp(1._rt - u*inv_um)) { reject = false; }
-            }
-        }
-
-        return u;
-    }
-}
-
 // struct whose getMomentum returns momentum for 1 particle, from random
 // gaussian flux distribution in the specified direction.
 // Along the normal axis, the distribution is v*Gaussian,