Implement injection of particles from the embedded boundary

Docs/source/usage/parameters.rst

@@ -971,16 +971,21 @@ Particle initialization
       The ``external_file`` option is currently implemented for 2D, 3D and RZ geometries, with record components in the cartesian coordinates ``(x,y,z)`` for 3D and RZ, and ``(x,z)`` for 2D.
       For more information on the `openPMD format <https://github.com/openPMD>`__ and how to build WarpX with it, please visit :ref:`the install section <install-developers>`.
 
-    * ``NFluxPerCell``: Continuously inject a flux of macroparticles from a planar surface.
+    * ``NFluxPerCell``: Continuously inject a flux of macroparticles from a surface. The emitting surface can be chosen to be either a plane
+      defined by the user (using some of the parameters listed below), or the embedded boundary (see :ref:`Embedded Boundary Conditions <running-cpp-parameters-eb>`).
       This requires the additional parameters:
 
       * ``<species_name>.flux_profile`` (see the description of this parameter further below)
 
-      * ``<species_name>.surface_flux_pos`` (`double`, location of the injection plane [meter])
+      * ``<species_name>.inject_from_embedded_boundary`` (`0` or `1`, default `0` ; whether to inject from the embedded boundary or from a user-specified plane.
+        When injecting from the embedded boundary, the momentum distribution specified by the user along ``z`` (see e.g. ``uz_m``, ``uz_th`` below) is interpreted
+        as the momentum distribution along the local normal to the embedded boundary.)
 
-      * ``<species_name>.flux_normal_axis`` (`x`, `y`, or `z` for 3D, `x` or `z` for 2D, or `r`, `t`, or `z` for RZ. When `flux_normal_axis` is `r` or `t`, the `x` and `y` components of the user-specified momentum distribution are interpreted as the `r` and `t` components respectively)
+      * ``<species_name>.surface_flux_pos`` (only used when injecting from a plane, `double`, location of the injection plane [meter])
 
-      * ``<species_name>.flux_direction`` (`-1` or `+1`, direction of flux relative to the plane)
+      * ``<species_name>.flux_normal_axis`` (only used when injecting from a plane, `x`, `y`, or `z` for 3D, `x` or `z` for 2D, or `r`, `t`, or `z` for RZ. When `flux_normal_axis` is `r` or `t`, the `x` and `y` components of the user-specified momentum distribution are interpreted as the `r` and `t` components respectively)
+
+      * ``<species_name>.flux_direction`` (only used when injecting from a plane, `-1` or `+1`, direction of flux relative to the plane)
 
       * ``<species_name>.num_particles_per_cell`` (`double`)
 

Examples/Tests/flux_injection/CMakeLists.txt

@@ -20,3 +20,33 @@ add_warpx_test(
     diags/diag1000120  # output
     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_from_eb.py  # analysis
+    diags/diag1000010  # output
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_rz_flux_injection_from_eb  # name
+    RZ  # dims
+    2  # nprocs
+    inputs_test_rz_flux_injection_from_eb  # inputs
+    analysis_flux_injection_from_eb.py  # analysis
+    diags/diag1000010  # output
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_2d_flux_injection_from_eb  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_flux_injection_from_eb  # inputs
+    analysis_flux_injection_from_eb.py  # analysis
+    diags/diag1000010  # output
+    OFF  # dependency
+)

Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py

@@ -0,0 +1,161 @@
+#!/usr/bin/env python3
+#
+# Copyright 2024 Remi Lehe
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+"""
+This script tests the emission of particles from the embedded boundary.
+(In this case, the embedded boundary is a sphere in 3D and RZ, a cylinder in 2D.)
+We check that the embedded boundary emits the correct number of particles, and that
+the particle distributions are consistent with the expected distributions.
+"""
+
+import os
+import re
+import sys
+
+import matplotlib.pyplot as plt
+import numpy as np
+import yt
+from scipy.constants import c, m_e
+from scipy.special import erf
+
+sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+import checksumAPI
+
+yt.funcs.mylog.setLevel(0)
+
+# Open plotfile specified in command line
+fn = sys.argv[1]
+ds = yt.load(fn)
+ad = ds.all_data()
+t_max = ds.current_time.item()  # time of simulation
+
+# Extract the dimensionality of the simulation
+with open("./warpx_used_inputs", "r") as f:
+    warpx_used_inputs = f.read()
+if re.search("geometry.dims = 2", warpx_used_inputs):
+    dims = "2D"
+elif re.search("geometry.dims = RZ", warpx_used_inputs):
+    dims = "RZ"
+elif re.search("geometry.dims = 3", warpx_used_inputs):
+    dims = "3D"
+
+# Total number of electrons expected:
+# Simulation parameters determine the total number of particles emitted (Ntot)
+flux = 1.0  # in m^-2.s^-1, from the input script
+R = 2.0  # in m, radius of the sphere
+if dims == "3D" or dims == "RZ":
+    emission_surface = 4 * np.pi * R**2  # in m^2
+elif dims == "2D":
+    emission_surface = 2 * np.pi * R  # in m
+Ntot = flux * emission_surface * t_max
+
+# Parameters of the histogram
+hist_bins = 50
+hist_range = [-0.5, 0.5]
+
+
+# Define function that histograms and checks the data
+def gaussian_dist(u, u_th):
+    return 1.0 / ((2 * np.pi) ** 0.5 * u_th) * np.exp(-(u**2) / (2 * u_th**2))
+
+
+def gaussian_flux_dist(u, u_th, u_m):
+    normalization_factor = u_th**2 * np.exp(-(u_m**2) / (2 * u_th**2)) + (
+        np.pi / 2
+    ) ** 0.5 * u_m * u_th * (1 + erf(u_m / (2**0.5 * u_th)))
+    result = (
+        1.0
+        / normalization_factor
+        * np.where(u > 0, u * np.exp(-((u - u_m) ** 2) / (2 * u_th**2)), 0)
+    )
+    return result
+
+
+def compare_gaussian(u, w, u_th, label=""):
+    du = (hist_range[1] - hist_range[0]) / hist_bins
+    w_hist, u_hist = np.histogram(u, bins=hist_bins, weights=w / du, range=hist_range)
+    u_hist = 0.5 * (u_hist[1:] + u_hist[:-1])
+    w_th = Ntot * gaussian_dist(u_hist, u_th)
+    plt.plot(u_hist, w_hist, label=label + ": simulation")
+    plt.plot(u_hist, w_th, "--", label=label + ": theory")
+    assert np.allclose(w_hist, w_th, atol=0.07 * w_th.max())
+
+
+def compare_gaussian_flux(u, w, u_th, u_m, label=""):
+    du = (hist_range[1] - hist_range[0]) / hist_bins
+    w_hist, u_hist = np.histogram(u, bins=hist_bins, weights=w / du, range=hist_range)
+    u_hist = 0.5 * (u_hist[1:] + u_hist[:-1])
+    w_th = Ntot * gaussian_flux_dist(u_hist, u_th, u_m)
+    plt.plot(u_hist, w_hist, label=label + ": simulation")
+    plt.plot(u_hist, w_th, "--", label=label + ": theory")
+    assert np.allclose(w_hist, w_th, atol=0.05 * w_th.max())
+
+
+# Load data and perform check
+
+plt.figure()
+
+if dims == "3D":
+    x = ad["electron", "particle_position_x"].to_ndarray()
+    y = ad["electron", "particle_position_y"].to_ndarray()
+    z = ad["electron", "particle_position_z"].to_ndarray()
+elif dims == "2D":
+    x = ad["electron", "particle_position_x"].to_ndarray()
+    y = np.zeros_like(x)
+    z = ad["electron", "particle_position_y"].to_ndarray()
+elif dims == "RZ":
+    theta = ad["electron", "particle_theta"].to_ndarray()
+    r = ad["electron", "particle_position_x"].to_ndarray()
+    x = r * np.cos(theta)
+    y = r * np.sin(theta)
+    z = ad["electron", "particle_position_y"].to_ndarray()
+ux = ad["electron", "particle_momentum_x"].to_ndarray() / (m_e * c)
+uy = ad["electron", "particle_momentum_y"].to_ndarray() / (m_e * c)
+uz = ad["electron", "particle_momentum_z"].to_ndarray() / (m_e * c)
+w = ad["electron", "particle_weight"].to_ndarray()
+
+# Check that the total number of particles emitted is correct
+Ntot_sim = np.sum(w)
+print("Ntot_sim = ", Ntot_sim)
+print("Ntot = ", Ntot)
+assert np.isclose(Ntot_sim, Ntot, rtol=0.01)
+
+# Check that none of the particles are inside the EB
+# A factor 0.98 is applied to accomodate
+# the cut-cell approximation of the sphere
+assert np.all(x**2 + y**2 + z**2 > (0.98 * R) ** 2)
+
+# Check that the normal component of the velocity is consistent with the expected distribution
+r = np.sqrt(x**2 + y**2 + z**2)
+nx = x / r
+ny = y / r
+nz = z / r
+u_n = ux * nx + uy * ny + uz * nz  # normal component
+compare_gaussian_flux(u_n, w, u_th=0.1, u_m=0.07, label="u_n")
+
+# Pick a direction that is orthogonal to the normal direction, and check the distribution
+vx = ny / np.sqrt(nx**2 + ny**2)
+vy = -nx / np.sqrt(nx**2 + ny**2)
+vz = 0
+u_perp = ux * vx + uy * vy + uz * vz
+compare_gaussian(u_perp, w, u_th=0.01, label="u_perp")
+
+# Pick the other perpendicular direction, and check the distribution
+# The third direction is obtained by the cross product (n x v)
+wx = ny * vz - nz * vy
+wy = nz * vx - nx * vz
+wz = nx * vy - ny * vx
+u_perp2 = ux * wx + uy * wy + uz * wz
+compare_gaussian(u_perp2, w, u_th=0.01, label="u_perp")
+
+plt.tight_layout()
+plt.savefig("Distribution.png")
+
+# Verify checksum
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn)

Examples/Tests/flux_injection/inputs_base_from_eb

@@ -0,0 +1,42 @@
+# Maximum number of time steps
+max_step = 10
+
+# 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 = 8
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Deactivate Maxwell solver
+algo.maxwell_solver = none
+warpx.const_dt = 1e-9
+
+# Embedded boundary
+warpx.eb_implicit_function = "-(x**2+y**2+z**2-2**2)"
+
+# particles
+particles.species_names = electron
+algo.particle_shape = 3
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = NFluxPerCell
+electron.inject_from_embedded_boundary = 1
+electron.num_particles_per_cell = 100
+electron.flux_profile = parse_flux_function
+electron.flux_function(x,y,z,t) = "1."
+electron.momentum_distribution_type = gaussianflux
+electron.ux_th = 0.01
+electron.uy_th = 0.01
+electron.uz_th = 0.1
+electron.uz_m = 0.07
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 10
+diag1.diag_type = Full
+diag1.fields_to_plot = none

Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb

@@ -0,0 +1,13 @@
+FILE = inputs_base_from_eb
+
+# number of grid points
+amr.n_cell = 16 16
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo     = -4 -4
+geometry.prob_hi     =  4  4
+
+# Boundary condition
+boundary.field_lo = periodic periodic
+boundary.field_hi = periodic periodic

Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb

@@ -0,0 +1,13 @@
+FILE = inputs_base_from_eb
+
+# number of grid points
+amr.n_cell = 16 16 16
+
+# Geometry
+geometry.dims = 3
+geometry.prob_lo     = -4 -4 -4
+geometry.prob_hi     =  4  4  4
+
+# Boundary condition
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic

Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb

@@ -0,0 +1,15 @@
+FILE = inputs_base_from_eb
+
+# number of grid points
+amr.n_cell = 8 16
+
+# Geometry
+geometry.dims = RZ
+geometry.prob_lo     =  0 -4
+geometry.prob_hi     =  4  4
+
+# Boundary condition
+boundary.field_lo = none periodic
+boundary.field_hi = pec periodic
+
+electron.num_particles_per_cell = 300

Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json

@@ -0,0 +1,11 @@
+{
+  "lev=0": {},
+  "electron": {
+    "particle_momentum_x": 6.990772711451971e-19,
+    "particle_momentum_y": 5.4131306169803364e-20,
+    "particle_momentum_z": 6.997294931789925e-19,
+    "particle_position_x": 35518.95120597846,
+    "particle_position_y": 35517.855675902414,
+    "particle_weight": 1.25355e-07
+  }
+}

Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json

@@ -0,0 +1,12 @@
+{
+  "lev=0": {},
+  "electron": {
+    "particle_momentum_x": 4.371688233196277e-18,
+    "particle_momentum_y": 4.368885079657374e-18,
+    "particle_momentum_z": 4.367429424105371e-18,
+    "particle_position_x": 219746.94401890738,
+    "particle_position_y": 219690.7015248918,
+    "particle_position_z": 219689.45580938633,
+    "particle_weight": 4.954974999999999e-07
+  }
+}

Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json

@@ -0,0 +1,12 @@
+{
+  "lev=0": {},
+  "electron": {
+    "particle_momentum_x": 6.734984863106283e-19,
+    "particle_momentum_y": 6.786279785869023e-19,
+    "particle_momentum_z": 1.0527983828124758e-18,
+    "particle_position_x": 53309.270966506396,
+    "particle_position_y": 53302.3776094842,
+    "particle_theta": 58707.74469425615,
+    "particle_weight": 4.991396867417661e-07
+  }
+}

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;