Implement injection of particles from the embedded boundary

Examples/Tests/flux_injection/CMakeLists.txt

@@ -11,6 +11,16 @@ add_warpx_test(
     OFF  # dependency
 )
 
+add_warpx_test(
+    test_3d_flux_injection_from_eb  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_flux_injection_from_eb  # inputs
+    analysis_flux_injection_3d.py  # analysis
+    diags/diag1000002  # output
+    OFF  # dependency
+)
+
 add_warpx_test(
     test_rz_flux_injection  # name
     RZ  # dims

Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb

@@ -0,0 +1,52 @@
+# Maximum number of time steps
+max_step = 10
+
+# number of grid points
+amr.n_cell = 16 16 16
+
+# The lo and hi ends of grids are multipliers of blocking factor
+amr.blocking_factor = 8
+
+# Maximum allowable size of each subdomain in the problem domain;
+#    this is used to decompose the domain for parallel calculations.
+amr.max_grid_size = 16
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 3
+geometry.prob_lo     = -4 -4 -4
+geometry.prob_hi     =  4  4  4
+
+# Deactivate Maxwell solver
+algo.maxwell_solver = none
+warpx.const_dt = 7e-9
+
+# Boundary condition
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+warpx.eb_implicit_function = "-(x**2+y**2+z**2-0.1**2)"
+
+# particles
+particles.species_names = electron
+algo.particle_shape = 3
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = NFluxPerCell
+electron.num_particles_per_cell = 100
+electron.inject_from_embedded_boundary = 1
+electron.flux_profile = parse_flux_function
+electron.flux_function(x,y,z,t) = "1."
+electron.momentum_distribution_type = gaussianflux
+electron.ux_th = 0.1
+electron.uy_th = 0.1
+electron.uz_th = 0.1
+electron.uz_m = 0.07
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 1000
+diag1.diag_type = Full
+diag1.fields_to_plot = none

Source/Initialization/PlasmaInjector.H

@@ -131,6 +131,8 @@ public:
     int flux_normal_axis;
     int flux_direction; // -1 for left, +1 for right
 
+    bool m_inject_from_eb = false; // whether to inject from the embedded boundary
+
     bool radially_weighted = true;
 
     std::string str_flux_function;

Source/Initialization/PlasmaInjector.cpp

@@ -9,6 +9,7 @@
  */
 #include "PlasmaInjector.H"
 
+#include "EmbeddedBoundary/Enabled.H"
 #include "Initialization/GetTemperature.H"
 #include "Initialization/GetVelocity.H"
 #include "Initialization/InjectorDensity.H"
@@ -303,50 +304,65 @@ void PlasmaInjector::setupNFluxPerCell (amrex::ParmParse const& pp_species)
         "(Please visit PR#765 for more information.)");
     }
 #endif
-    utils::parser::getWithParser(pp_species, source_name, "surface_flux_pos", surface_flux_pos);
-    utils::parser::queryWithParser(pp_species, source_name, "flux_tmin", flux_tmin);
-    utils::parser::queryWithParser(pp_species, source_name, "flux_tmax", flux_tmax);
-    std::string flux_normal_axis_string;
-    utils::parser::get(pp_species, source_name, "flux_normal_axis", flux_normal_axis_string);
-    flux_normal_axis = -1;
+
+    // Check whether injection from the embedded boundary is requested
+    utils::parser::queryWithParser(pp_species, source_name, "inject_from_embedded_boundary", m_inject_from_eb);
+    if (m_inject_from_eb) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( EB::enabled(),
+            "Error: Embedded boundary injection is only available when "
+            "embedded boundaries are enabled.");
+        flux_normal_axis = 2; // Interpret z as the normal direction to the EB
+        flux_direction = 1;
+    } else {
+        // Injection is through a plane in this case.
+        // Parse the parameters of the plane (position, normal direction, etc.)
+
+        utils::parser::getWithParser(pp_species, source_name, "surface_flux_pos", surface_flux_pos);
+        utils::parser::queryWithParser(pp_species, source_name, "flux_tmin", flux_tmin);
+        utils::parser::queryWithParser(pp_species, source_name, "flux_tmax", flux_tmax);
+        std::string flux_normal_axis_string;
+        utils::parser::get(pp_species, source_name, "flux_normal_axis", flux_normal_axis_string);
+        flux_normal_axis = -1;
 #ifdef WARPX_DIM_RZ
-    if      (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") {
-        flux_normal_axis = 0;
-    }
-    if      (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") {
-        flux_normal_axis = 1;
-    }
+        if      (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") {
+            flux_normal_axis = 0;
+        }
+        if      (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") {
+            flux_normal_axis = 1;
+        }
 #else
 #    ifndef WARPX_DIM_1D_Z
-    if      (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") {
-        flux_normal_axis = 0;
-    }
+        if      (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") {
+            flux_normal_axis = 0;
+        }
 #    endif
 #endif
 #ifdef WARPX_DIM_3D
-    if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") {
-        flux_normal_axis = 1;
-    }
+        if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") {
+            flux_normal_axis = 1;
+        }
 #endif
-    if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") {
-        flux_normal_axis = 2;
-    }
+        if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") {
+            flux_normal_axis = 2;
+        }
 #ifdef WARPX_DIM_3D
-    const std::string flux_normal_axis_help = "'x', 'y', or 'z'.";
+        const std::string flux_normal_axis_help = "'x', 'y', or 'z'.";
 #else
 #    ifdef WARPX_DIM_RZ
-    const std::string flux_normal_axis_help = "'r' or 'z'.";
+        const std::string flux_normal_axis_help = "'r' or 'z'.";
 #    elif WARPX_DIM_XZ
-    const std::string flux_normal_axis_help = "'x' or 'z'.";
+        const std::string flux_normal_axis_help = "'x' or 'z'.";
 #    else
-    const std::string flux_normal_axis_help = "'z'.";
+        const std::string flux_normal_axis_help = "'z'.";
 #    endif
-#endif
-    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0,
-        "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help);
-    utils::parser::getWithParser(pp_species, source_name, "flux_direction", flux_direction);
-    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1,
-        "Error: flux_direction must be -1 or +1.");
+    #endif
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0,
+            "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help);
+        utils::parser::getWithParser(pp_species, source_name, "flux_direction", flux_direction);
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1,
+            "Error: flux_direction must be -1 or +1.");
+    }
+
     // Construct InjectorPosition with InjectorPositionRandom.
     h_flux_pos = std::make_unique<InjectorPosition>(
         (InjectorPositionRandomPlane*)nullptr,

Source/Particles/AddPlasmaUtilities.H

@@ -71,12 +71,46 @@ int compute_area_weights (const amrex::IntVect& iv, const int normal_axis) {
     return r;
 }
 
+
+#ifdef AMREX_USE_EB
+/*
+  This computes the scale_fac (used for setting the particle weights) on a on area basis
+  (used for flux injection from the embedded boundary).
+ */
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real compute_scale_fac_area_eb (
+    const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
+    const amrex::Real num_ppc_real,
+    amrex::Array4<const amrex::Real> const& eb_bnd_area_arr,
+    amrex::Array4<const amrex::Real> const& eb_bnd_normal_arr,
+    int i, int j, int k ) {
+    using namespace amrex::literals;
+    // Scale particle weight by the area of the emitting surface, within one cell
+    amrex::Real scale_fac = 0.0_rt;
+#if defined(WARPX_DIM_3D)
+    // By definition, eb_bnd_area_arr is normalized (unitless).
+    // Here we undo the normalization (i.e. multiply by the surface used for normalization).
+    const amrex::Real nx = eb_bnd_normal_arr(i,j,k,0);
+    const amrex::Real ny = eb_bnd_normal_arr(i,j,k,1);
+    const amrex::Real nz = eb_bnd_normal_arr(i,j,k,2);
+    scale_fac = eb_bnd_area_arr(i,j,k) *
+        std::sqrt(amrex::Math::powi<2>(nx*dx[1]*dx[2]) +
+                  amrex::Math::powi<2>(ny*dx[0]*dx[2]) +
+                  amrex::Math::powi<2>(nz*dx[0]*dx[1]));
+#else
+    // TODO Implement in RZ, 2D, and 1D
+#endif
+    scale_fac /= num_ppc_real;
+    return scale_fac;
+}
+#endif
+
 /*
   This computes the scale_fac (used for setting the particle weights) on a on area basis
-  (used for flux injection).
+  (used for flux injection from a plane).
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-amrex::Real compute_scale_fac_area (const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
+amrex::Real compute_scale_fac_area_plane (const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
                                     const amrex::Real num_ppc_real, const int flux_normal_axis) {
     using namespace amrex::literals;
     amrex::Real scale_fac = AMREX_D_TERM(dx[0],*dx[1],*dx[2])/num_ppc_real;