Merge pull request #103 from ECP-WarpX/factory

Use virtual function to select the spectral algorithm

8 files changed, 136 insertions, 62 deletions

diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package
index 50127914d..b0ee54095 100644
--- a/Source/FieldSolver/SpectralSolver/Make.package
+++ b/Source/FieldSolver/SpectralSolver/Make.package
@@ -1,11 +1,12 @@
 CEXE_headers += WarpX_Complex.H
 CEXE_headers += SpectralSolver.H
+CEXE_sources += SpectralSolver.cpp
 CEXE_headers += SpectralFieldData.H
 CEXE_sources += SpectralFieldData.cpp
-CEXE_headers += PsatdAlgorithm.H
-CEXE_sources += PsatdAlgorithm.cpp
 CEXE_headers += SpectralKSpace.H
 CEXE_sources += SpectralKSpace.cpp
 
+include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver
diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H
deleted file mode 100644
index acefcc466..000000000
--- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef WARPX_PSATD_ALGORITHM_H_
-#define WARPX_PSATD_ALGORITHM_H_
-
-#include <SpectralKSpace.H>
-#include <SpectralFieldData.H>
-
-/* \brief Class that updates the field in spectral space
- * and stores the coefficients of the corresponding update equation.
- */
-class PsatdAlgorithm
-{
-    using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >;
-
-    public:
-        PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
-                         const amrex::DistributionMapping& dm,
-                         const int norder_x, const int norder_y,
-                         const int norder_z, const bool nodal, const amrex::Real dt);
-        PsatdAlgorithm() = default; // Default constructor
-        PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default;
-        void pushSpectralFields(SpectralFieldData& f) const;
-
-    private:
-        // Modified finite-order vectors
-        KVectorComponent modified_kx_vec, modified_kz_vec;
-#if (AMREX_SPACEDIM==3)
-        KVectorComponent modified_ky_vec;
-#endif
-        SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
-};
-
-#endif // WARPX_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
new file mode 100644
index 000000000..c62c21f44
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -0,0 +1,6 @@
+CEXE_headers += SpectralBaseAlgorithm.H
+CEXE_headers += PsatdAlgorithm.H
+CEXE_sources += PsatdAlgorithm.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
new file mode 100644
index 000000000..0487e5226
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -0,0 +1,24 @@
+#ifndef WARPX_PSATD_ALGORITHM_H_
+#define WARPX_PSATD_ALGORITHM_H_
+
+#include <SpectralBaseAlgorithm.H>
+
+/* \brief Class that updates the field in spectral space
+ * and stores the coefficients of the corresponding update equation.
+ */
+class PsatdAlgorithm : public SpectralBaseAlgorithm
+{
+    public:
+        PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
+                         const amrex::DistributionMapping& dm,
+                         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
+        virtual void pushSpectralFields(SpectralFieldData& f) const override final;
+
+    private:
+        SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
+};
+
+#endif // WARPX_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index ada7506c3..37892d35a 100644
--- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -9,14 +9,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
                          const DistributionMapping& dm,
                          const int norder_x, const int norder_y,
                          const int norder_z, const bool nodal, const Real dt)
-// Compute and assign the modified k vectors
-: modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)),
-#if (AMREX_SPACEDIM==3)
-  modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)),
-  modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal))
-#else
-  modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal))
-#endif
+     // Initialize members of base class
+     : SpectralBaseAlgorithm( spectral_kspace, dm,
+                              norder_x, norder_y, norder_z, nodal )
 {
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
new file mode 100644
index 000000000..602eb2473
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -0,0 +1,51 @@
+#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_
+#define WARPX_SPECTRAL_BASE_ALGORITHM_H_
+
+#include <SpectralKSpace.H>
+#include <SpectralFieldData.H>
+
+/* \brief Class that updates the field in spectral space
+ * and stores the coefficients of the corresponding update equation.
+ *
+ * `SpectralBaseAlgorithm` is only a base class and cannot be used directly.
+ * Instead use its subclasses, which implement the specific field update
+ * equations for a given spectral algorithm.
+ */
+class SpectralBaseAlgorithm
+{
+    public:
+        // Member function that updates the fields in spectral space ;
+        // meant to be overridden in subclasses
+        virtual void pushSpectralFields(SpectralFieldData& f) const = 0;
+        // The destructor should also be a virtual function, so that
+        // a pointer to subclass of `SpectraBaseAlgorithm` actually
+        // calls the subclass's destructor.
+        virtual ~SpectralBaseAlgorithm() {};
+
+    protected: // Meant to be used in the subclasses
+
+        using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >;
+
+        // Constructor
+        SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace,
+                              const amrex::DistributionMapping& dm,
+                              const int norder_x, const int norder_y,
+                              const int norder_z, const bool nodal)
+          // Compute and assign the modified k vectors
+          : modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)),
+#if (AMREX_SPACEDIM==3)
+            modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)),
+            modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal))
+#else
+            modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal))
+#endif
+          {};
+
+        // Modified finite-order vectors
+        KVectorComponent modified_kx_vec, modified_kz_vec;
+#if (AMREX_SPACEDIM==3)
+        KVectorComponent modified_ky_vec;
+#endif
+};
+
+#endif // WARPX_SPECTRAL_BASE_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 7444452af..d4019a9a3 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -1,8 +1,7 @@
 #ifndef WARPX_SPECTRAL_SOLVER_H_
 #define WARPX_SPECTRAL_SOLVER_H_
 
-#include <SpectralKSpace.H>
-#include <PsatdAlgorithm.H>
+#include <SpectralBaseAlgorithm.H>
 #include <SpectralFieldData.H>
 
 /* \brief Top-level class for the electromagnetic spectral solver
@@ -14,28 +13,17 @@
 class SpectralSolver
 {
     public:
-        // Inline definition of the member functions of `SpectralSolver`
+        // Inline definition of the member functions of `SpectralSolver`,
+        // except the constructor (see `SpectralSolver.cpp`)
         // The body of these functions is short, since the work is done in the
         // underlying classes `SpectralFieldData` and `PsatdAlgorithm`
 
-        /* \brief Initialize the spectral solver */
+        // Constructor
         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 ) {
-            // Initialize all structures using the same distribution mapping dm
-
-            // - Initialize k space object (Contains info about the size of
-            // the spectral space corresponding to each box in `realspace_ba`,
-            // as well as the value of the corresponding k coordinates)
-            const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx);
-            // - Initialize the algorithm (coefficients) over this space
-            algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y,
-                                        norder_z, nodal, dt );
-            // - Initialize arrays for fields in Fourier space + FFT plans
-            field_data = SpectralFieldData( realspace_ba, k_space, dm );
-        };
+                        const amrex::RealVect dx, const amrex::Real dt );
 
         /* \brief Transform the component `i_comp` of MultiFab `mf`
          *  to spectral space, and store the corresponding result internally
@@ -59,14 +47,20 @@ class SpectralSolver
         /* \brief Update the fields in spectral space, over one timestep */
         void pushSpectralFields(){
             BL_PROFILE("SpectralSolver::pushSpectralFields");
-            algorithm.pushSpectralFields( field_data );
+            // Virtual function: the actual function used here depends
+            // on the sub-class of `SpectralBaseAlgorithm` that was
+            // initialized in the constructor of `SpectralSolver`
+            algorithm->pushSpectralFields( field_data );
         };
 
     private:
         SpectralFieldData field_data; // Store field in spectral space
                                       // and perform the Fourier transforms
-        PsatdAlgorithm algorithm; // Contains Psatd coefficients
-                                  // and field update equation
+        std::unique_ptr<SpectralBaseAlgorithm> algorithm;
+        // Defines field update equation in spectral space,
+        // and the associated coefficients.
+        // SpectralBaseAlgorithm is a base class ; this pointer is meant
+        // to point an instance of a *sub-class* defining a specific algorithm
 };
 
 #endif // WARPX_SPECTRAL_SOLVER_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
new file mode 100644
index 000000000..c21c3cfb1
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -0,0 +1,35 @@
+#include <SpectralKSpace.H>
+#include <SpectralSolver.H>
+#include <PsatdAlgorithm.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.
+ */
+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 ) {
+
+    // Initialize all structures using the same distribution mapping dm
+
+    // - Initialize k space object (Contains info about the size of
+    // the spectral space corresponding to each box in `realspace_ba`,
+    // as well as the value of the corresponding k coordinates)
+    const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx);
+
+    // - 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(
+            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 );
+
+};