Implement PML for the outer RZ boundary with PSATD

Examples/Tests/PML/analysis_pml_psatd_rz.py

@@ -0,0 +1,53 @@
+#! /usr/bin/env python
+
+# Copyright 2021 David Grote
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+"""
+This script tests the absorption of fields in the PML in RZ geometry.
+
+The input file inputs_particle_rz is used: it features an electron
+moving radially that launches a pulse. This scripts runs until
+most of the pulse escapes the radial boundary. If the PML fails,
+the pulse will remain with in the domain.
+"""
+import sys
+import os
+import numpy as np
+import yt
+yt.funcs.mylog.setLevel(0)
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# Open plotfile specified in command line
+filename = sys.argv[1]
+ds = yt.load( filename )
+
+# yt 4.0+ has rounding issues with our domain data:
+# RuntimeError: yt attempted to read outside the boundaries
+#               of a non-periodic domain along dimension 0.
+if 'force_periodicity' in dir(ds): ds.force_periodicity()
+
+# Check that the field is low enough
+ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex_array = ad0['boxlib', 'Ex'].to_ndarray()
+Ez_array = ad0['boxlib', 'Ez'].to_ndarray()
+max_Ex = np.abs(Ex_array).max()
+max_Ez = np.abs(Ez_array).max()
+print( f'max Ex = {max_Ex}' )
+print( f'max Ez = {max_Ez}' )
+max_Efield = max(max_Ex, max_Ez)
+
+# This tolerance was obtained empirically. As the simulation progresses, the field energy is leaking
+# out through PML so that the max field diminishes with time. When the PML is working properly,
+# the field level falls below 2 at the end of the simulation.
+tolerance_abs = 2.
+print('tolerance_abs: ' + str(tolerance_abs))
+assert max_Efield < tolerance_abs
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Tests/PML/inputs_rz

@@ -0,0 +1,48 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 500
+amr.n_cell = 32 128
+amr.max_level = 0
+geometry.coord_sys   = 1  # 1: Cylindrical
+geometry.prob_lo     =    0.     -100.e-6
+geometry.prob_hi     =  50.e-6   +100.e-6
+
+warpx.n_rz_azimuthal_modes = 2
+
+#################################
+######## Boundary condition #####
+#################################
+boundary.field_lo = none periodic
+boundary.field_hi = pml  periodic
+
+# PML
+warpx.pml_ncell = 10
+warpx.do_pml_in_domain = 0
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.maxwell_solver = psatd
+warpx.use_filter = 0
+algo.particle_shape = 1
+
+#################################
+############ PARTICLE ###########
+#################################
+particles.species_names = electron
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "singleparticle"
+electron.single_particle_pos = 0. 0. 0.
+electron.single_particle_vel = 10. 0. 0.
+electron.single_particle_weight = 1.
+
+#################################
+########## DIAGNOSTICS ##########
+#################################
+diagnostics.diags_names = diag1
+diag1.intervals = 500
+diag1.fields_to_plot = Bx By Bz Ex Ey Ez jx jy jz rho
+diag1.diag_type = Full

Regression/Checksum/benchmarks_json/pml_psatd_rz.json

@@ -0,0 +1,14 @@
+{
+  "lev=0": {
+    "Ex": 1317.2337418101674,
+    "Ey": 7.056477765087081e-12,
+    "Ez": 1156.1377910176022,
+    "Bx": 2.594232124519114e-20,
+    "By": 3.2688462949057505e-06,
+    "Bz": 4.292128402662449e-21,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0,
+    "rho": 0.0
+  }
+}

Regression/WarpX-tests.ini

@@ -126,6 +126,22 @@ doVis = 0
 analysisRoutine = Examples/analysis_default_regression.py
 tolerance = 1.e-14
 
+[pml_psatd_rz]
+buildDir = .
+inputFile = Examples/Tests/PML/inputs_rz
+runtime_params =
+dim = 2
+addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+analysisRoutine = Examples/Tests/PML/analysis_pml_psatd_rz.py
+tolerance = 1.e-14
+
 [silver_mueller_2d_x]
 buildDir = .
 inputFile = Examples/Tests/SilverMueller/inputs_2d_x

Source/BoundaryConditions/CMakeLists.txt

@@ -5,3 +5,10 @@ target_sources(WarpX
     WarpXFieldBoundaries.cpp
     WarpX_PEC.cpp
 )
+
+if(WarpX_DIMS STREQUAL RZ)
+    target_sources(WarpX
+      PRIVATE
+        PML_RZ.cpp
+    )
+endif()

Source/BoundaryConditions/Make.package

@@ -1,5 +1,8 @@
 CEXE_sources += PML.cpp WarpXEvolvePML.cpp
 CEXE_sources += WarpXFieldBoundaries.cpp WarpX_PEC.cpp
 
+ifeq ($(USE_RZ),TRUE)
+    CEXE_sources += PML_RZ.cpp
+endif
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/BoundaryConditions

Source/BoundaryConditions/PML_RZ.H

@@ -0,0 +1,75 @@
+/* Copyright 2021 David Grote
+ *
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PML_RZ_H_
+#define WARPX_PML_RZ_H_
+
+#include "PML_RZ_fwd.H"
+
+#ifdef WARPX_USE_PSATD
+#   include "FieldSolver/SpectralSolver/SpectralSolverRZ.H"
+#endif
+
+#include <AMReX_MultiFab.H>
+#include <AMReX_BoxArray.H>
+#include <AMReX_Config.H>
+#include <AMReX_REAL.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <array>
+#include <string>
+
+enum struct PatchType : int;
+
+class PML_RZ
+{
+public:
+    PML_RZ (const int lev, const amrex::BoxArray& grid_ba, const amrex::DistributionMapping& grid_dm,
+            const amrex::Geometry* geom, const int ncell, const int do_pml_in_domain);
+
+    void ApplyDamping(amrex::MultiFab* Et_fp, amrex::MultiFab* Ez_fp,
+                      amrex::MultiFab* Bt_fp, amrex::MultiFab* Bz_fp,
+                      amrex::Real dt);
+
+    std::array<amrex::MultiFab*,2> GetE_fp ();
+    std::array<amrex::MultiFab*,2> GetB_fp ();
+
+#ifdef WARPX_USE_PSATD
+    void PushPSATD (const int lev);
+#endif
+
+    void FillBoundaryE ();
+    void FillBoundaryB ();
+    void FillBoundaryE (PatchType patch_type);
+    void FillBoundaryB (PatchType patch_type);
+
+    void CheckPoint (const std::string& dir) const;
+    void Restart (const std::string& dir);
+
+    ~PML_RZ () = default;
+
+private:
+
+    const int m_ncell;
+    const int m_do_pml_in_domain;
+    const amrex::Geometry* m_geom;
+
+    // Only contains Er and Et, and Br and Bt
+    std::array<std::unique_ptr<amrex::MultiFab>,2> pml_E_fp;
+    std::array<std::unique_ptr<amrex::MultiFab>,2> pml_B_fp;
+
+#ifdef WARPX_USE_PSATD
+    void PushPMLPSATDSinglePatchRZ ( const int lev,
+                SpectralSolverRZ& solver,
+                std::array<std::unique_ptr<amrex::MultiFab>,2>& pml_E,
+                std::array<std::unique_ptr<amrex::MultiFab>,2>& pml_B);
+#endif
+
+};
+
+#endif