Add convolutional PML (CPML) and fix PEC-mask handling at the PML interface

Docs/source/usage/parameters.rst

@@ -381,6 +381,51 @@ Additional PML parameters
     The characteristic depth, in number of cells, over which
     the absorption coefficients of the PML increases.
 
+* ``warpx.pml_type`` (`string`; default: ``split``)
+    PML formulation used with the FDTD solver.
+
+    * ``split``: split-field (Berenger) PML (original implementation).
+
+    * ``convolutional`` (or ``cpml``): convolutional PML (Roden & Gedney, 2000), using the
+      complex-frequency-shifted coordinate stretching
+      :math:`s_i = \kappa_i + \sigma_i/(\alpha_i + i\omega)`. The total (unsplit) field is
+      advanced in the PML with the *same* material update as in the interior (local
+      :math:`\epsilon`, :math:`\mu`, :math:`\sigma`, including spatially varying properties and
+      the ``laxwendroff``/``backwardeuler`` macroscopic schemes); only the spatial derivatives are
+      stretched, through auxiliary convolution variables. The PEC mask built from the GDS/numpy
+      layer (``algo.use_PEC_mask``) and its extension into the PML is applied to the CPML fields
+      in the same way as for the split-field PML. The :math:`\sigma` profile is the same as for the
+      split-field PML (set by ``warpx.pml_ncell`` and ``warpx.pml_delta``).
+      Only implemented for the (staggered) Yee solver, without divergence cleaning in the PML,
+      J-damping, or moving window.
+      With the default parameters (``kappa_max = 1``, ``alpha_max = 0``) its absorption is
+      essentially the same as the split-field PML. Using ``warpx.pml_delta`` around 15-20 together
+      with a small ``warpx.cpml_alpha_max`` (and optionally ``warpx.cpml_kappa_max = 2``)
+      typically lowers the reflection further; a too-large ``alpha_max`` degrades absorption
+      of the low frequencies.
+
+* ``algo.extend_PEC_mask_into_pml`` (`int`; default: 1)
+    Only used with ``algo.use_PEC_mask = 1``. If 1, the PEC mask built from the numpy (GDS)
+    layer is extended into the PML (the mask values on the domain edge are continued through
+    the PML), for both ``warpx.pml_type = split`` and ``convolutional``. If 0, no PEC mask is
+    applied in the PML.
+
+* ``warpx.cpml_kappa_max`` (`float`; default: 1)
+    CPML only. Maximum of the real stretching factor :math:`\kappa`, graded as
+    :math:`\kappa(d) = 1 + (\kappa_{max}-1)\,d^m` with :math:`d \in [0,1]` the normalized
+    depth into the PML. Values > 1 improve absorption of evanescent waves.
+
+* ``warpx.cpml_alpha_max`` (`float`, in 1/s; default: 0)
+    CPML only. Maximum of the complex-frequency-shift parameter, given as the rate
+    :math:`\alpha_{SI}/\epsilon_0`, graded as :math:`\alpha(d) = \alpha_{max} (1-d)^{m_a}`.
+    A typical choice is :math:`2\pi f_{min}`, with :math:`f_{min}` the lowest frequency of interest.
+
+* ``warpx.cpml_grading_order`` (`float`; default: 2)
+    CPML only. Polynomial order :math:`m` of the :math:`\kappa` grading.
+
+* ``warpx.cpml_alpha_grading_order`` (`float`; default: 1)
+    CPML only. Polynomial order :math:`m_a` of the :math:`\alpha` grading.
+
 * ``warpx.do_pml_in_domain`` (`int`; default: 0)
     Whether to create the PML inside the simulation area or outside. If inside,
     it allows the user to propagate particles in PML and to use extended PML

Source/BoundaryConditions/CPML_K.H

@@ -0,0 +1,138 @@
+/* Convolutional PML (CPML) helpers.
+ *
+ * The CPML uses the complex-frequency-shifted coordinate stretching
+ *
+ *     s_d = kappa_d + sigma_d / (alpha_d + i omega)       (d = x, y, z)
+ *
+ * (sigma_d and alpha_d are rates, i.e. sigma_SI/eps0 and alpha_SI/eps0), which turns
+ * each spatial derivative into
+ *
+ *     d/dd  ->  (1/kappa_d) d/dd + psi_d
+ *
+ * where psi_d is the recursive convolution (Roden & Gedney, 2000)
+ *
+ *     psi_d^{new} = b_d psi_d^{old} + a_d (dF/dd),
+ *     b_d = exp(-(sigma_d/kappa_d + alpha_d) dt),
+ *     a_d = sigma_d / (kappa_d (sigma_d + kappa_d alpha_d)) (b_d - 1).
+ *
+ * Because the stretching only acts on the spatial derivatives, the material
+ * response (eps, mu, sigma of the medium) is untouched: the E update keeps the
+ * same alpha/beta (Lax-Wendroff or backward-Euler) coefficients as in the
+ * interior, evaluated with the local eps(x) and sigma(x). This is what makes the
+ * CPML reflectionless for arbitrary (and spatially varying) materials continued
+ * into the PML, in contrast to formulations whose damping must be matched to eps.
+ */
+#ifndef WARPX_CPML_K_H_
+#define WARPX_CPML_K_H_
+
+#include "BoundaryConditions/PML.H"
+
+#include <AMReX.H>
+#include <AMReX_Array4.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_IndexType.H>
+#include <AMReX_REAL.H>
+
+#include <cmath>
+
+/** Device-copyable view of the CPML 1D profiles of one SigmaBox, indexed by the
+ *  physical direction (0=x, 1=y, 2=z). Directions that are not resolved by the
+ *  simulation (y in 2D XZ) have null pointers and act as "no stretching". */
+struct CPMLCoefs
+{
+    const amrex::Real* sig [3] = {nullptr, nullptr, nullptr}; // nodal positions
+    const amrex::Real* sigs[3] = {nullptr, nullptr, nullptr}; // cell-centered positions
+    const amrex::Real* kap [3] = {nullptr, nullptr, nullptr};
+    const amrex::Real* kaps[3] = {nullptr, nullptr, nullptr};
+    const amrex::Real* alp [3] = {nullptr, nullptr, nullptr};
+    const amrex::Real* alps[3] = {nullptr, nullptr, nullptr};
+    int lo[3] = {0, 0, 0};
+
+    CPMLCoefs () = default;
+
+    explicit CPMLCoefs (SigmaBox const& sb)
+    {
+#if defined(WARPX_DIM_3D)
+        constexpr int dirs[AMREX_SPACEDIM] = {0, 1, 2};
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+        constexpr int dirs[AMREX_SPACEDIM] = {0, 2};
+#else
+        constexpr int dirs[AMREX_SPACEDIM] = {2};
+#endif
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            const int d = dirs[idim];
+            sig [d] = sb.sigma[idim].data();
+            sigs[d] = sb.sigma_star[idim].data();
+            kap [d] = sb.kappa[idim].data();
+            kaps[d] = sb.kappa_star[idim].data();
+            alp [d] = sb.alpha[idim].data();
+            alps[d] = sb.alpha_star[idim].data();
+            lo  [d] = sb.sigma[idim].lo();
+        }
+    }
+};
+
+/** Staggering of a field component, indexed by physical direction (1 = nodal) */
+inline amrex::GpuArray<int,3> cpml_stag (amrex::IndexType const& ixt)
+{
+    amrex::IntVect const iv = ixt.toIntVect();
+#if defined(WARPX_DIM_3D)
+    return {iv[0], iv[1], iv[2]};
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+    return {iv[0], 0, iv[1]};
+#else
+    return {0, 0, iv[0]};
+#endif
+}
+
+/** Grid index along physical direction dir */
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int cpml_index (int dir, int i, int j, int k) noexcept
+{
+#if defined(WARPX_DIM_3D)
+    return (dir == 0) ? i : ((dir == 1) ? j : k);
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+    amrex::ignore_unused(k);
+    return (dir == 0) ? i : j;
+#else
+    amrex::ignore_unused(dir, j, k);
+    return i;
+#endif
+}
+
+/** Coefficients (b, a, 1/kappa) of the CPML recursive convolution along direction
+ *  dir, at grid index ii, for a field with nodal (1) or cell-centered (0) staggering
+ *  along dir, for a time step dt. */
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cpml_coef (CPMLCoefs const& c, int dir, int nodal, int ii, amrex::Real dt,
+                amrex::Real& b, amrex::Real& a, amrex::Real& inv_kappa) noexcept
+{
+    using namespace amrex::literals;
+    if (c.sig[dir] == nullptr) {
+        b = 1._rt; a = 0._rt; inv_kappa = 1._rt;
+        return;
+    }
+    const int n = ii - c.lo[dir];
+    const amrex::Real s  = nodal ? c.sig[dir][n] : c.sigs[dir][n];
+    const amrex::Real kp = nodal ? c.kap[dir][n] : c.kaps[dir][n];
+    const amrex::Real al = nodal ? c.alp[dir][n] : c.alps[dir][n];
+    inv_kappa = 1._rt/kp;
+    b = std::exp(-(s*inv_kappa + al)*dt);
+    const amrex::Real denom = kp*(s + kp*al);
+    a = (denom > 0._rt) ? s/denom*(b - 1._rt) : 0._rt;
+}
+
+/** Stretched derivative: if update_psi, first advances the memory variable psi
+ *  (in place) over dt_psi; returns (1/kappa) dF + psi */
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real cpml_stretched_derivative (CPMLCoefs const& c, int dir, int nodal,
+                                       int i, int j, int k, amrex::Real dt_psi, bool update_psi,
+                                       amrex::Real dF, amrex::Real& psi) noexcept
+{
+    amrex::Real b, a, ik;
+    cpml_coef(c, dir, nodal, cpml_index(dir, i, j, k), dt_psi, b, a, ik);
+    if (update_psi) { psi = b*psi + a*dF; }
+    return ik*dF + psi;
+}
+
+#endif // WARPX_CPML_K_H_

Source/BoundaryConditions/PML.H

@@ -70,6 +70,20 @@ struct SigmaBox
     SigmaVect sigma_star_cumsum_fac;
     amrex::Real v_sigma;
 
+    // Convolutional PML (CPML) coordinate-stretching parameters,
+    // s_i = kappa_i + sigma_i/(alpha_i + j*omega), using the same staggering as
+    // sigma (nodal) and sigma_star (cell-centered). sigma and alpha are rates (1/s),
+    // i.e. sigma_SI/eps0 and alpha_SI/eps0.
+    SigmaVect kappa;
+    SigmaVect kappa_star;
+    SigmaVect alpha;
+    SigmaVect alpha_star;
+    amrex::IntVect m_ncell;
+    amrex::Real m_cpml_kappa_max = 1.0;
+    amrex::Real m_cpml_alpha_max = 0.0;
+    amrex::Real m_cpml_grading_order = 2.0;
+    amrex::Real m_cpml_alpha_grading_order = 1.0;
+
 };
 
 class SigmaBoxFactory
@@ -163,6 +177,20 @@ public:
     amrex::MultiFab* GetE_fp (int comp);
     amrex::MultiFab* GetPEC_fp (int comp);
 
+    /** Convolutional PML (warpx.pml_type = convolutional) auxiliary fields.
+     *  psi_E[d] lives at the location of E_d and stores, in its two components,
+     *  the recursive convolutions of the two curl(H) derivatives entering the
+     *  update of E_d (same component ordering as PMLComp, e.g. Ex: xy=0 (d/dy), xz=1 (d/dz)).
+     *  psi_B[d] is the analog for the B (or H, with LLG) update, using curl(E).
+     *  These return nullptr with the split-field PML. */
+    std::array<amrex::MultiFab*,3> GetPsiE_fp ();
+    std::array<amrex::MultiFab*,3> GetPsiB_fp ();
+    std::array<amrex::MultiFab*,3> GetPsiE_cp ();
+    std::array<amrex::MultiFab*,3> GetPsiB_cp ();
+
+    /** True if this PML uses the convolutional (CPML) formulation */
+    bool isConvolutional () const { return m_convolutional; }
+
     // Used when WarpX::do_pml_dive_cleaning = true
     amrex::MultiFab* GetF_fp ();
     amrex::MultiFab* GetF_cp ();
@@ -252,6 +280,13 @@ private:
 
     std::array<std::unique_ptr<amrex::MultiFab>,3> pml_PEC_fp;
 
+    // CPML auxiliary (memory) variables
+    bool m_convolutional = false;
+    std::array<std::unique_ptr<amrex::MultiFab>,3> pml_psiE_fp;
+    std::array<std::unique_ptr<amrex::MultiFab>,3> pml_psiB_fp;
+    std::array<std::unique_ptr<amrex::MultiFab>,3> pml_psiE_cp;
+    std::array<std::unique_ptr<amrex::MultiFab>,3> pml_psiB_cp;
+
     std::array<std::unique_ptr<amrex::MultiFab>,3> pml_edge_lengths;
 
     std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_cp;