1 files changed, 68 insertions, 0 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);
+            }
+
+        );
+
+    }
+
+};