1 files changed, 19 insertions, 5 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index c21c3cfb1..4b9def013 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -1,19 +1,29 @@
 #include <SpectralKSpace.H>
 #include <SpectralSolver.H>
 #include <PsatdAlgorithm.H>
+#include <PMLPsatdAlgorithm.H>
 
 /* \brief Initialize the spectral Maxwell solver
  *
  * This function selects the spectral algorithm to be used, allocates the
  * corresponding coefficients for the discretized field update equation,
  * and prepares the structures that store the fields in spectral space.
+ *
+ * \param norder_x Order of accuracy of the spatial derivatives along x
+ * \param norder_y Order of accuracy of the spatial derivatives along y
+ * \param norder_z Order of accuracy of the spatial derivatives along z
+ * \param nodal    Whether the solver is applied to a nodal or staggered grid
+ * \param dx       Cell size along each dimension
+ * \param dt       Time step
+ * \param pml      Whether the boxes in which the solver is applied are PML boxes
  */
 SpectralSolver::SpectralSolver(
                 const amrex::BoxArray& realspace_ba,
                 const amrex::DistributionMapping& dm,
                 const int norder_x, const int norder_y,
                 const int norder_z, const bool nodal,
-                const amrex::RealVect dx, const amrex::Real dt ) {
+                const amrex::RealVect dx, const amrex::Real dt,
+                const bool pml ) {
 
     // Initialize all structures using the same distribution mapping dm
 
@@ -24,12 +34,16 @@ SpectralSolver::SpectralSolver(
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
-    // TODO: Add more algorithms + selection depending on input parameters
-    //       For the moment, this only uses the standard PsatdAlgorithm
-    algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm(
+    if (pml) {
+        algorithm = std::unique_ptr<PMLPsatdAlgorithm>( new PMLPsatdAlgorithm(
+            k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) );
+    } else {
+        algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm(
             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() );
 
 };