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

Docs/source/refs.bib

@@ -479,3 +479,15 @@ @article{VayFELB2009
   doi = {10.1063/1.3080930},
   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},
+}

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,81 @@
+#!/usr/bin/env python3
+
+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")