Implement injection of particles from the embedded boundary (ECP-Warp…

…X#5208)

# Overview

This PR implements flux injection of particles from the embedded
boundary.

It also adds a test that emits particles from a sphere in 3D as
represented here:

![movie](https://github.com/user-attachments/assets/1e76cf87-fd7d-4fa3-8c83-363956226a42)
as well as RZ and 2D versions of this test. (In 2D, the particles are
emitted from a cylinder.)

As can be seen in the above movie, particles are emitted from a single
point within each cell (the centroid of the EB), instead of being
emitted uniformly on the surface of the EB within the cell. This could
be improved in a future PR.

The implementation as well as the user interface largely re-use the
infrastructure for the flux injection from a plane. However, as a
result, the user interface is perhaps not very intuitive. In particular,
when specify the velocity distribution, `uz` represents the direction
normal to the EB while `ux`, `uy` represent the tangential directions.
This again will be improved in follow-up PR.

# Follow-up PRs

- [ ] Change the interface of `gaussianflux` so as to specify the
tangential and normal distribution. In other words, instead of:
```
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
```
we would do:
```
electron.momentum_distribution_type = gaussianflux
electron.u_tangential_th = 0.01  # Tangential to the emitting surface
electron.u_normal_th = 0.1   # Normal to the emitting surface
electron.u_normal_m = 0.07
```

- [ ] Change the interface so that the user does not need to specify the
number of macroparticles per cell (which is problematic for EB, since
difference cell contain different EB surface, and should in general emit
different numbers of macroparticles). Instead, we would specify the
weight of macroparticles, i.e. instead of
```
electron.injection_style = NFluxPerCell
electron.num_particles_per_cell = 100   
electron.flux_function(x,y,z,t) = "1.”
```
we would do
```
electron.injection_style = NFluxPerCell
electron.flux_macroweight = 200   # Number of physical particles per macroparticle
electron.flux_function(x,y,z,t) = "4e12” # Number of physical particles emitted per unit time and surface  
```

- [ ] Add a way for the user to specify the total flux across the whole
emitting surface
Example:
```
electron.flux_function(x,y,z,t) = "(x>-1)*(x<1)"
electron.total_flux = 4e12 # physical particle / second (not per unit area)
```
(In that case, `flux_function` would be rescaled internally by WarpX so
as to emit the right number of particles.)

- [ ] Add PICMI interface
- [ ] Emit the particles uniformly from the surface of the EB within one
cell

Docs/source/usage/parameters.rst

@@ -1002,16 +1002,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;