Start implementation of spectral PML

6 files changed, 26 insertions, 13 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index c62c21f44..ee8376865 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,6 +1,8 @@
 CEXE_headers += SpectralBaseAlgorithm.H
 CEXE_headers += PsatdAlgorithm.H
 CEXE_sources += PsatdAlgorithm.cpp
+CEXE_headers += PMLPsatdAlgorithm.H
+CEXE_sources += PMLPsatdAlgorithm.cpp
 
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 0487e5226..52e587e7f 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -14,8 +14,11 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
                          const int norder_x, const int norder_y,
                          const int norder_z, const bool nodal,
                          const amrex::Real dt);
-        // Redefine update equation from base class
+        // Redefine functions from base class
         virtual void pushSpectralFields(SpectralFieldData& f) const override final;
+        virtual int getRequiredNumberOfFields() const override final {
+            return SpectralFieldIndex::n_fields;
+        }
 
     private:
         SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 602eb2473..5d5e376c1 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -14,9 +14,9 @@
 class SpectralBaseAlgorithm
 {
     public:
-        // Member function that updates the fields in spectral space ;
-        // meant to be overridden in subclasses
+        // Virtual member function ; meant to be overridden in subclasses
         virtual void pushSpectralFields(SpectralFieldData& f) const = 0;
+        virtual int getRequiredNumberOfFields() const = 0;
         // The destructor should also be a virtual function, so that
         // a pointer to subclass of `SpectraBaseAlgorithm` actually
         // calls the subclass's destructor.
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 8e58aa1d8..30cf3733b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -8,18 +8,24 @@
 // Declare type for spectral fields
 using SpectralField = amrex::FabArray< amrex::BaseFab <Complex> >;
 
-/* Index for the fields that will be stored in spectral space */
+/* Index for the regular fields, when stored in spectral space */
 struct SpectralFieldIndex {
   enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields };
   // n_fields is automatically the total number of fields
 };
 
+/* Index for the PML fields, when stored in spectral space */
+struct SpectralPMLIndex {
+  enum { Exy=0, Exz, Eyx, Eyz, Ezx, Ezy,
+           Bxy, Bxz, Byx, Byz, Bzx, Bzy, n_fields };
+  // n_fields is automatically the total number of fields
+};
+
 /* \brief Class that stores the fields in spectral space, and performs the
  *  Fourier transforms between real space and spectral space
  */
 class SpectralFieldData
 {
-    friend class PsatdAlgorithm;
 
     // Define the FFTplans type, which holds one fft plan per box
     // (plans are only initialized for the boxes that are owned by
@@ -32,8 +38,9 @@ class SpectralFieldData
 
     public:
         SpectralFieldData( const amrex::BoxArray& realspace_ba,
-                      const SpectralKSpace& k_space,
-                      const amrex::DistributionMapping& dm );
+                           const SpectralKSpace& k_space,
+                           const amrex::DistributionMapping& dm,
+                           const int n_field_required );
         SpectralFieldData() = default; // Default constructor
         SpectralFieldData& operator=(SpectralFieldData&& field_data) = default;
         ~SpectralFieldData();
@@ -41,10 +48,10 @@ class SpectralFieldData
                                const int field_index, const int i_comp);
         void BackwardTransform( amrex::MultiFab& mf,
                                const int field_index, const int i_comp);
-
-    private:
         // `fields` stores fields in spectral space, as multicomponent FabArray
         SpectralField fields;
+
+    private:
         // tmpRealField and tmpSpectralField store fields
         // right before/after the Fourier transform
         SpectralField tmpSpectralField; // contains Complexs
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index 02fa2015f..c45809dd5 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -5,14 +5,14 @@ using namespace amrex;
 /* \brief Initialize fields in spectral space, and FFT plans */
 SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba,
                             const SpectralKSpace& k_space,
-                            const DistributionMapping& dm )
+                            const DistributionMapping& dm,
+                            const int n_field_required )
 {
     const BoxArray& spectralspace_ba = k_space.spectralspace_ba;
 
     // Allocate the arrays that contain the fields in spectral space
     // (one component per field)
-    fields = SpectralField(spectralspace_ba, dm,
-                            SpectralFieldIndex::n_fields, 0);
+    fields = SpectralField(spectralspace_ba, dm, n_field_required, 0);
 
     // Allocate temporary arrays - in real space and spectral space
     // These arrays will store the data just before/after the FFT
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index c21c3cfb1..a91fcbc47 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -30,6 +30,7 @@ SpectralSolver::SpectralSolver(
             k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) );
 
     // - Initialize arrays for fields in spectral space + FFT plans
-    field_data = SpectralFieldData( realspace_ba, k_space, dm );
+    field_data = SpectralFieldData( realspace_ba, k_space, dm,
+            algorithm->getRequiredNumberOfFields() );
 
 };