Add Dan Barnes' effective potential (semi-implicit) Poisson solver (#5079)

This PR adds a new electrostatic solver based on the semi-implicit
scheme developed by [Barnes, Journal of Comp. Phys., 424 (2021),
109852](https://www.sciencedirect.com/science/article/pii/S0021999120306264?via%3Dihub)
(see Appendix A of reference).

The implementation was tested through the simulation of a Langmuir probe
placed inside a uniform plasma. To this end a 2d simulation was
performed with an initial uniform plasma with a conducting disk inserted
in the center of the domain. The disk is kept at 0 V while Neumann
boundary conditions are used for the domain boundary. Electrons and
hydrogen ions were injected from all the domain boundaries using the
`UniformFluxDistribution` in order to simulate the particle flux from a
uniform plasma. Thermal boundaries were used for the particles with the
same temperatures as the initial plasma particles. The expected outcome
of this simulation is that a sheath develops around the "Langmuir probe"
with value
$V_s = 0.5T_e\ln\left(2\pi\frac{m_e}{2m_p}(1 + \frac{T_i}{T_e})\right)
$.
A pre-sheath of roughly $0.7T_e$ is also expected to form but the exact
value of this depends on the domain size since the pre-sheath extends
for many Debye lengths.
The figure below shows an example outcome from a simulation as described
above in which the semi-implicit factor was set to $C_{SI} = 10$.

![image](https://github.com/user-attachments/assets/de4c4f44-3eb1-4e4b-8297-4c047d1c0d98)

Required before merging:

- [x] Merging of AMReX-Codes/amrex#3968
- [x] Add example to CI tests
- [x] Update documentation

---------

Signed-off-by: roelof-groenewald <regroenewald@gmail.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Thomas Marks <marksta@umich.edu>

Author: 3 people

Docs/source/refs.bib

@@ -480,6 +480,18 @@ @article{VayFELB2009
   url = {https://doi.org/10.1063/1.3080930},
 }
 
+@article{Barnes2021,
+  title = {Improved C1 shape functions for simplex meshes},
+  author = {D.C. Barnes},
+  journal = {Journal of Computational Physics},
+  volume = {424},
+  pages = {109852},
+  year = {2021},
+  issn = {0021-9991},
+  doi = {https://doi.org/10.1016/j.jcp.2020.109852},
+  url = {https://www.sciencedirect.com/science/article/pii/S0021999120306264},
+}
+
 @article{Rhee1987,
     author = {Rhee, M. J. and Schneider, R. F. and Weidman, D. J.},
     title = "{Simple time‐resolving Thomson spectrometer}",

Docs/source/usage/parameters.rst

@@ -236,6 +236,23 @@ Overall simulation parameters
         \boldsymbol{\nabla}^2 \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi \\
         \boldsymbol{\nabla}^2 \boldsymbol{A} = - \mu_0 \boldsymbol{j} \qquad \boldsymbol{B} = \boldsymbol{\nabla}\times\boldsymbol{A}
 
+    * ``labframe-effective-potential``: Poisson's equation is solved with a modified dielectric function
+      (resulting in an "effective potential") to create a semi-implicit scheme which is robust to the
+      numerical instability seen in explicit electrostatic PIC when :math:`\Delta t \omega_{pe} > 2`.
+      If this option is used the additional parameter ``warpx.effective_potential_factor`` can also be
+      specified to set the value of :math:`C_{EP}` (default 4). The method is stable for :math:`C_{EP} \geq 1`
+      regardless of :math:`\Delta t`, however, the larger :math:`C_{EP}` is set, the lower the numerical plasma
+      frequency will be and therefore care must be taken to not set it so high that the plasma mode
+      hybridizes with other modes of interest.
+      Details of the method can be found in Appendix A of :cite:t:`param-Barnes2021` (note that in that paper
+      the method is referred to as "semi-implicit electrostatic" but here it has been renamed to "effective potential"
+      to avoid confusion with the semi-implicit method of Chen et al.).
+      In short, the code solves:
+
+      .. math::
+
+        \boldsymbol{\nabla}\cdot\left(1+\frac{C_{EP}}{4}\sum_{s \in \text{species}}(\omega_{ps}\Delta t)^2 \right)\boldsymbol{\nabla} \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi
+
     * ``relativistic``: Poisson's equation is solved **for each species**
       in their respective rest frame. The corresponding field
       is mapped back to the simulation frame and will produce both E and B

Examples/Tests/CMakeLists.txt

@@ -71,6 +71,7 @@ add_subdirectory(restart)
 add_subdirectory(restart_eb)
 add_subdirectory(rigid_injection)
 add_subdirectory(scraping)
+add_subdirectory(effective_potential_electrostatic)
 add_subdirectory(silver_mueller)
 add_subdirectory(single_particle)
 add_subdirectory(space_charge_initialization)

Examples/Tests/effective_potential_electrostatic/CMakeLists.txt

@@ -0,0 +1,12 @@
+# Add tests (alphabetical order) ##############################################
+#
+
+add_warpx_test(
+    test_3d_effective_potential_electrostatic_picmi  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_effective_potential_electrostatic_picmi.py  # inputs
+    analysis.py  # analysis
+    diags/field_diag/  # output
+    OFF  # dependency
+)

Examples/Tests/effective_potential_electrostatic/analysis.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+#
+# --- Analysis script for the effective potential Poisson solver test. This test
+# --- is based on the adiabatic plasma expansion benchmark from Connor et al. (2021)
+# --- doi.org/10.1109/TPS.2021.3072353.
+# --- The electron density distribution (as a function of radius) is compared
+# --- with the analytically calculated density based on the input parameters
+# --- of the test simulation at each output timestep.
+
+import os
+import sys
+
+import dill
+import matplotlib.pyplot as plt
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.interpolate import RegularGridInterpolator
+
+from pywarpx import picmi
+
+constants = picmi.constants
+
+# load simulation parameters
+with open("sim_parameters.dpkl", "rb") as f:
+    sim = dill.load(f)
+
+# characteristic expansion time
+tau = sim["sigma_0"] * np.sqrt(sim["M"] / (constants.kb * (sim["T_e"] + sim["T_i"])))
+
+
+def get_analytic_density(r, t):
+    expansion_factor = 1.0 + t**2 / tau**2
+    T = sim["T_e"] / expansion_factor
+    sigma = sim["sigma_0"] * np.sqrt(expansion_factor)
+    return (
+        sim["n_plasma"] * (T / sim["T_e"]) ** 1.5 * np.exp(-(r**2) / (2.0 * sigma**2))
+    )
+
+
+def get_radial_function(field, info):
+    """Helper function to transform Cartesian data to polar data and average
+    over theta and phi."""
+    r_grid = np.linspace(0, np.max(info.z) - 1e-3, 30)
+    theta_grid = np.linspace(0, 2 * np.pi, 40, endpoint=False)
+    phi_grid = np.linspace(0, np.pi, 20)
+
+    r, t, p = np.meshgrid(r_grid, theta_grid, phi_grid, indexing="ij")
+    x_sp = r * np.sin(p) * np.cos(t)
+    y_sp = r * np.sin(p) * np.sin(t)
+    z_sp = r * np.cos(p)
+
+    interp = RegularGridInterpolator((info.x, info.y, info.z), field)
+    field_sp = interp((x_sp, y_sp, z_sp))
+    return r_grid, np.mean(field_sp, axis=(1, 2))
+
+
+diag_dir = "diags/field_diag"
+ts = OpenPMDTimeSeries(diag_dir, check_all_files=True)
+
+rms_errors = np.zeros(len(ts.iterations))
+
+for ii, it in enumerate(ts.iterations):
+    rho_e, info = ts.get_field(field="rho_electrons", iteration=it)
+    r_grid, n_e = get_radial_function(-rho_e / constants.q_e, info)
+
+    n_e_analytic = get_analytic_density(r_grid, ts.t[ii])
+    rms_errors[ii] = (
+        np.sqrt(np.mean(np.sum((n_e - n_e_analytic) ** 2))) / n_e_analytic[0]
+    )
+
+    plt.plot(r_grid, n_e_analytic, "k--", alpha=0.6)
+    plt.plot(r_grid, n_e, label=f"t = {ts.t[ii]*1e6:.2f} $\mu$s")
+
+print("RMS error (%) in density: ", rms_errors)
+assert np.all(rms_errors < 0.05)
+
+plt.ylabel("$n_e$ (m$^{-3}$)")
+plt.xlabel("r (m)")
+plt.grid()
+plt.legend()
+plt.show()
+
+if len(sys.argv) > 1:
+    sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+    import checksumAPI
+
+    filename = sys.argv[1]
+
+    test_name = os.path.split(os.getcwd())[1]
+    checksumAPI.evaluate_checksum(test_name, filename, output_format="openpmd")

Examples/Tests/effective_potential_electrostatic/inputs_test_3d_effective_potential_electrostatic_picmi.py

@@ -0,0 +1,236 @@
+#!/usr/bin/env python3
+#
+# --- Test script for the effective potential Poisson solver. This test is based
+# --- on the adiabatic plasma expansion benchmark from Connor et al. (2021)
+# --- doi.org/10.1109/TPS.2021.3072353.
+# --- In the benchmark an expanding plasma ball with Gaussian density distribution
+# --- in the radial direction is simulated and the time evolution of the
+# --- density of the electron species is compared to an approximate analytic solution.
+# --- The example is modified slightly in the following ways:
+# --- 1) The original example used an electromagnetic solver with absorbing
+# ---    boundaries while the present case encloses the plasma in a conducting
+# ---    sphere.
+# --- 2) The domain and plasma parameters for this case has been modified to
+# ---    set the cell-size and time step such that the explicit electrostatic
+# ---    solver is unstable.
+
+import dill
+import numpy as np
+from mpi4py import MPI as mpi
+from scipy.special import erf
+
+from pywarpx import picmi
+
+constants = picmi.constants
+
+comm = mpi.COMM_WORLD
+
+simulation = picmi.Simulation(warpx_serialize_initial_conditions=True, verbose=0)
+
+m_ion = 25  # Ion mass (electron masses)
+
+# Plasma parameters
+n_plasma = 5e12  # Plasma density (m^-3)
+sigma_0 = 22  # Initial Gaussian distribution standard deviation (Debye lengths)
+T_e = 100.0  # Electron temperature (K)
+T_i = 10.0  # Ion temperature (K)
+
+# Spatial domain
+R = 86  # Radius of metallic sphere (Debye lengths)
+NZ = 72  # Number of cells in each direction
+
+# Temporal domain (if not run as a CI test)
+LT = 0.6e-6  # Simulation temporal length (s)
+
+# Numerical parameters
+NPARTS = 500000  # Seed number of particles
+DT = 0.8  # Time step (electron streaming)
+
+# Solver parameter
+C_EP = 1.0  # Effective potential factor
+
+#######################################################################
+# Calculate various plasma parameters based on the simulation input   #
+#######################################################################
+
+# Ion mass (kg)
+M = m_ion * constants.m_e
+
+# Electron plasma frequency (Hz)
+f_pe = np.sqrt(constants.q_e**2 * n_plasma / (constants.m_e * constants.ep0)) / (
+    2.0 * np.pi
+)
+
+# Debye length (m)
+lambda_e = np.sqrt(constants.ep0 * constants.kb * T_e / (n_plasma * constants.q_e**2))
+
+# Thermal velocities (m/s) from v_th = np.sqrt(kT / m)
+v_ti = np.sqrt(T_i * constants.kb / M)
+v_te = np.sqrt(T_e * constants.kb / constants.m_e)
+
+R *= lambda_e
+sigma_0 *= lambda_e
+
+dz = 2.0 * R / (NZ - 4)
+Lz = dz * NZ
+dt = DT * dz / v_te
+
+total_steps = int(LT / dt)
+diag_steps = total_steps // 3
+total_steps = diag_steps * 3
+
+# dump attributes needed for analysis to a dill pickle file
+if comm.rank == 0:
+    parameter_dict = {
+        "sigma_0": sigma_0,
+        "M": M,
+        "T_i": T_i,
+        "T_e": T_e,
+        "n_plasma": n_plasma,
+    }
+    with open("sim_parameters.dpkl", "wb") as f:
+        dill.dump(parameter_dict, f)
+
+# print out plasma parameters
+if comm.rank == 0:
+    print(
+        f"Initializing simulation with input parameters:\n"
+        f"\tT_e = {T_e:.1f} K\n"
+        f"\tT_i = {T_i:.1f} K\n"
+        f"\tn = {n_plasma:.1e} m^-3\n"
+    )
+    print(
+        f"Plasma parameters:\n"
+        f"\tlambda_e = {lambda_e:.1e} m\n"
+        f"\tt_pe = {1.0/f_pe:.1e} s\n"
+        f"\tv_ti = {v_ti:.1e} m/s\n"
+    )
+    print(
+        f"Numerical parameters:\n"
+        f"\tdz/lambda_e = {dz/lambda_e:.2f}\n"
+        f"\tdt*w_pe = {dt*f_pe*2.0*np.pi:.2f}\n"
+        f"\tdiag steps = {diag_steps:d}\n"
+        f"\ttotal steps = {total_steps:d}\n"
+    )
+
+
+#######################################################################
+# Set geometry and boundary conditions                                #
+#######################################################################
+
+grid = picmi.Cartesian3DGrid(
+    number_of_cells=[NZ] * 3,
+    lower_bound=[-Lz / 2.0] * 3,
+    upper_bound=[Lz / 2.0] * 3,
+    lower_boundary_conditions=["neumann"] * 3,
+    upper_boundary_conditions=["neumann"] * 3,
+    lower_boundary_conditions_particles=["absorbing"] * 3,
+    upper_boundary_conditions_particles=["absorbing"] * 3,
+    warpx_max_grid_size=NZ // 2,
+)
+simulation.time_step_size = dt
+simulation.max_steps = total_steps
+simulation.current_deposition_algo = "direct"
+simulation.particle_shape = 1
+simulation.verbose = 1
+
+#######################################################################
+# Insert spherical boundary as EB                                     #
+#######################################################################
+
+embedded_boundary = picmi.EmbeddedBoundary(
+    implicit_function=f"(x**2+y**2+z**2-{R**2})",
+    potential=0.0,
+)
+simulation.embedded_boundary = embedded_boundary
+
+#######################################################################
+# Field solver                                                        #
+#######################################################################
+
+solver = picmi.ElectrostaticSolver(
+    grid=grid,
+    method="Multigrid",
+    warpx_effective_potential=True,
+    warpx_effective_potential_factor=C_EP,
+    warpx_self_fields_verbosity=0,
+)
+simulation.solver = solver
+
+#######################################################################
+# Particle types setup                                                #
+#######################################################################
+
+total_parts = (
+    n_plasma
+    * sigma_0**2
+    * (
+        (2.0 * np.pi) ** 1.5 * sigma_0 * erf(R / (np.sqrt(2) * sigma_0))
+        + 4.0 * np.pi * R * np.exp(-(R**2) / (2.0 * sigma_0**2))
+    )
+)
+
+electrons = picmi.Species(
+    name="electrons",
+    particle_type="electron",
+    initial_distribution=picmi.GaussianBunchDistribution(
+        n_physical_particles=total_parts,
+        rms_bunch_size=[sigma_0] * 3,
+        rms_velocity=[v_te] * 3,
+    ),
+)
+simulation.add_species(
+    electrons,
+    layout=picmi.PseudoRandomLayout(grid=grid, n_macroparticles=NPARTS),
+)
+
+ions = picmi.Species(
+    name="ions",
+    charge="q_e",
+    mass=M,
+    initial_distribution=picmi.GaussianBunchDistribution(
+        n_physical_particles=total_parts,
+        rms_bunch_size=[sigma_0] * 3,
+        rms_velocity=[v_ti] * 3,
+    ),
+)
+simulation.add_species(
+    ions,
+    layout=picmi.PseudoRandomLayout(grid=grid, n_macroparticles=NPARTS),
+)
+
+#######################################################################
+# Add diagnostics                                                     #
+#######################################################################
+
+field_diag = picmi.FieldDiagnostic(
+    name="field_diag",
+    grid=grid,
+    period=diag_steps,
+    data_list=[
+        "E",
+        "J",
+        "T_electrons",
+        "T_ions",
+        "phi",
+        "rho_electrons",
+        "rho_ions",
+    ],
+    write_dir="diags",
+    warpx_format="openpmd",
+    warpx_openpmd_backend="h5",
+)
+simulation.add_diagnostic(field_diag)
+
+#######################################################################
+# Initialize simulation                                               #
+#######################################################################
+
+simulation.initialize_inputs()
+simulation.initialize_warpx()
+
+#######################################################################
+# Execute simulation                                                  #
+#######################################################################
+
+simulation.step()