Started implementing finite difference solver

Added Yee algorithm

7 files changed, 224 insertions, 4 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
new file mode 100644
index 000000000..e5514c697
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -0,0 +1,68 @@
+// TODO include statements
+
+FiniteDifferenceSolver::EvolveB ( VectorField Bfield,
+                                  ConstVectorField Efield,
+                                  amrex::Real dt ) {
+    // Select algorithm (The choice of algorithm is a runtime option,
+    // but we compile code for each algorithm, using templates)
+    if (fdtd_algo == MaxwellSolverAlgo::Yee){
+        EvolveB <YeeAlgorithm> ( Bfield, Efield, dt );
+    } else if (fdtd_algo == MaxwellSolverAlgo::CKC) {
+        EvolveB <CKCAlgorithm> ( Bfield, Efield, dt );
+    } else {
+        amrex::Abort("Unknown algorithm");
+    }
+)
+
+template<typename algo>
+FiniteDifferenceSolver::EvolveB ( VectorField Bfield,
+                                 ConstVectorField Efield,
+                                 amrex::Real dt ) {
+
+    // Loop through the grids, and over the tiles within each grid
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+        // Extract field data for this grid/tile
+        auto const& Bx = Bfield[0]->array(mfi);
+        auto const& By = Bfield[1]->array(mfi);
+        auto const& Bz = Bfield[2]->array(mfi);
+        auto const& Ex = Efield[0]->array(mfi);
+        auto const& Ey = Efield[1]->array(mfi);
+        auto const& Ez = Efield[2]->array(mfi);
+
+        // Extract stencil coefficients
+        Real const* AMREX_RESTRICT coefs_x = stencil_coefs_x.dataPtr();
+        Real const* AMREX_RESTRICT coefs_y = stencil_coefs_y.dataPtr();
+        Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr();
+
+        // Extract tileboxes for which to loop
+        const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
+        const Box& tby  = mfi.tilebox(By_nodal_flag);
+        const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
+
+        // Loop over the cells and update the fields
+        amrex::ParallelFor(tbx, tby, tbz,
+
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+                Bx(i, j, k) += dt * algo::UpwardDz( Ey, i, j, k, coefs_z)
+                             - dt * algo::UpwardDy( Ez, i, j, k, coefs_y);
+            },
+
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+                By(i, j, k) += dt * algo::UpwardDx( Ez, i, j, k, coefs_x)
+                             - dt * algo::UpwardDz( Ex, i, j, k, coefs_z);
+            },
+
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+                Bz(i, j, k) += dt * algo::UpwardDy( Ex, i, j, k, coefs_y)
+                             - dt * algo::UpwardDx( Ey, i, j, k, coefs_x);
+            }
+
+        );
+
+    }
+
+};
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H
new file mode 100644
index 000000000..7c9572298
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H
@@ -0,0 +1,80 @@
+#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_
+#define WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_
+
+struct YeeAlgorithm {
+
+    void InitializeStencilCoefficients( std::array<Real,3> cell_size,
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_x,
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_y,
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z ) {
+
+        // Store the inverse cell size along each direction in the coefficients
+        stencil_coefs_x.resize(1);
+        stencil_coefs_x[0] = 1./cell_size[0];
+        stencil_coefs_y.resize(1);
+        stencil_coefs_y[0] = 1./cell_size[1];
+        stencil_coefs_z.resize(1);
+        stencil_coefs_z[0] = 1./cell_size[2];
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real UpwardDx(
+        Array4 const F, int const i, int const j, int const k,
+        Real const* coefs_x, int const n_coefs_x ) {
+
+        amrex::Real inv_dx = coefs_x[0];
+        return inv_dx*( F(i+1,j,k) - F(i,j,k) );
+    };
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real DownwardDx(
+        Array4 const F, int const i, int const j, int const k,
+        Real const* coefs_x, int const n_coefs_x ) {
+
+        amrex::Real inv_dx = coefs_x[0];
+        return inv_dx*( F(i,j,k) - F(i-1,j,k) );
+    };
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real UpwardDy(
+        Array4 const F, int const i, int const j, int const k,
+        Real const* coefs_y, int const n_coefs_y ) {
+
+        #if defined WARPX_DIM_3D
+            amrex::Real inv_dy = coefs_y[0];
+            return inv_dy*( F(i,j+1,k) - F(i,j,k) );
+        #elif (defined WARPX_DIM_XZ)
+            return 0; // 2D Cartesian: derivative along y is 0
+        #endif
+    };
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real DownwardDy(
+        Array4 const F, int const i, int const j, int const k,
+        Real const* coefs_y, int const n_coefs_y ) {
+
+        #if defined WARPX_DIM_3D
+            amrex::Real inv_dy = coefs_y[0];
+            return inv_dy*( F(i,j,k) - F(i,j-1,k) );
+        #elif (defined WARPX_DIM_XZ)
+            return 0; // 2D Cartesian: derivative along y is 0
+        #endif
+    };
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real UpwardDz(
+        Array4 const F, int const i, int const j, int const k,
+        Real const* coefs_z, int const n_coefs_z ) {
+
+        amrex::Real inv_dz = coefs_z[0];
+        #if defined WARPX_DIM_3D
+            return inv_dz*( F(i,j,k+1) - F(i,j,k) );
+        #elif (defined WARPX_DIM_XZ)
+            return inv_dz*( F(i,j+1,k) - F(i,j,k) );
+        #endif
+    };
+
+
+}
+
+#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
new file mode 100644
index 000000000..c1fdef5dd
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -0,0 +1,38 @@
+#ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_
+#define WARPX_FINITE_DIFFERENCE_SOLVER_H_
+
+/**
+ * \brief Top-level class for the electromagnetic finite-difference solver
+ *
+ * TODO
+ */
+class FiniteDifferenceSolver
+{
+    public:
+
+        using VectorField = std::array< std::unique_ptr<amrex::MultiFab>, 3 >;
+        using ConstVectorField = std::array< std::unique_ptr<amrex::MultiFab const>, 3 >;
+
+        // Constructor
+        void FiniteDifferenceSolver( std::array<Real,3> dx );
+
+        void EvolveB( VectorField Bfield,
+                    ConstVectorField Efield,
+                    amrex::Real dt ) const;
+
+        };
+
+    private:
+
+        int fdtd_algo;
+
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_x;
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_y;
+        amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z;
+
+        template< typename fdtd_algo >
+        void EvolveB
+
+};
+
+#endif // WARPX_FINITE_DIFFERENCE_SOLVER_H_
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
new file mode 100644
index 000000000..a146b66b1
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
@@ -0,0 +1,15 @@
+
+
+void FiniteDifferenceSolver::FiniteDifferenceSolver ( std::array<Real,3> cell_size ) {
+    // Select algorithm (The choice of algorithm is a runtime option,
+    // but we compile code for each algorithm, using templates)
+    if (fdtd_algo == MaxwellSolverAlgo::Yee){
+        YeeAlgorithm::InitializeStencilCoefficients( cell_size,
+            stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+    } else if (fdtd_algo == MaxwellSolverAlgo::CKC) {
+        CKCAlgorithm::InitializeStencilCoefficients( cell_size,
+            stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+    } else {
+        amrex::Abort("Unknown algorithm");
+    }
+};
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 4848b051e..1137da741 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -111,6 +111,13 @@ WarpX::EvolveB (int lev, amrex::Real a_dt)
 void
 WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
 {
+    if (patch_type == PatchType::fine) {
+        fdtd_solver_fp->EvolveB( Bfield_fp[lev], Efield_fp[lev], a_dt );
+    } else {
+        fdtd_solver_cp->EvolveB( Bfield_cp[lev], Efield_cp[lev], a_dt );
+    }
+
+    // Goes into initializer
     const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
     const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
     const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2];
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 3bab73833..6775c982a 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -694,6 +694,9 @@ private:
 
     amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_fp;
     amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_cp;
+#else
+    amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> fdtd_solver_fp;
+    amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> fdtd_solver_cp;
 #endif
 
 #ifdef WARPX_USE_PSATD_HYBRID
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 48b4bbd55..50e73d4c1 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -223,9 +223,13 @@ WarpX::WarpX ()
 
     costs.resize(nlevs_max);
 
+    // Allocate field solver objects
 #ifdef WARPX_USE_PSATD
     spectral_solver_fp.resize(nlevs_max);
     spectral_solver_cp.resize(nlevs_max);
+#else
+    fdtd_solver_fp.resize(nlevs_max);
+    fdtd_solver_cp.resize(nlevs_max);
 #endif
 #ifdef WARPX_USE_PSATD_HYBRID
     Efield_fp_fft.resize(nlevs_max);
@@ -889,11 +893,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     if (fft_hybrid_mpi_decomposition == false){
         // Allocate and initialize the spectral solver
         std::array<Real,3> dx = CellSize(lev);
-#if (AMREX_SPACEDIM == 3)
+    #if (AMREX_SPACEDIM == 3)
         RealVect dx_vect(dx[0], dx[1], dx[2]);
-#elif (AMREX_SPACEDIM == 2)
+    #elif (AMREX_SPACEDIM == 2)
         RealVect dx_vect(dx[0], dx[2]);
-#endif
+    #endif
         // Get the cell-centered box, with guard cells
         BoxArray realspace_ba = ba;  // Copy box
         realspace_ba.enclosedCells().grow(ngE); // cell-centered + guard cells
@@ -901,8 +905,10 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         spectral_solver_fp[lev].reset( new SpectralSolver( realspace_ba, dm,
             nox_fft, noy_fft, noz_fft, do_nodal, dx_vect, dt[lev] ) );
     }
+#else
+        std::array<Real,3> dx = CellSize(lev);
+        fdtd_solver_fp[lev].reset( new FiniteDifferenceSolver(dx); )
 #endif
-
     //
     // The Aux patch (i.e., the full solution)
     //
@@ -986,6 +992,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
             spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm,
                 nox_fft, noy_fft, noz_fft, do_nodal, cdx_vect, dt[lev] ) );
         }
+#else
+    std::array<Real,3> dx = CellSize(lev-1);
+    fdtd_solver_cp[lev].reset( new FiniteDifferenceSolver(dx); )
 #endif
     }