Implement div(B) Cleaning With FDTD (#1829)

* Implement div(B) Cleaning With FDTD

* Add CI Test

* Clean Up

1 files changed, 31 insertions, 5 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index 3ddfc66b0..3d2dc42d8 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -25,6 +25,7 @@ using namespace amrex;
 void FiniteDifferenceSolver::EvolveB (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::unique_ptr<amrex::MultiFab> const& Gfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
     int lev, amrex::Real const dt ) {
 
@@ -38,21 +39,20 @@ void FiniteDifferenceSolver::EvolveB (
 #else
     if (m_do_nodal) {
 
-        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
 #endif
     } else {
         amrex::Abort("EvolveB: Unknown algorithm");
     }
-
 }
 
 
@@ -62,6 +62,7 @@ template<typename T_Algo>
 void FiniteDifferenceSolver::EvolveBCartesian (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::unique_ptr<amrex::MultiFab> const& Gfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
     int lev, amrex::Real const dt ) {
 
@@ -71,6 +72,8 @@ void FiniteDifferenceSolver::EvolveBCartesian (
 
     amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
+    Real constexpr c2 = PhysConst::c * PhysConst::c;
+
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
@@ -138,9 +141,32 @@ void FiniteDifferenceSolver::EvolveBCartesian (
                 Bz(i, j, k) += dt * T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k)
                              - dt * T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k);
             }
-
         );
 
+        // div(B) cleaning correction for errors in magnetic Gauss law (div(B) = 0)
+        if (Gfield)
+        {
+            // Extract field data for this grid/tile
+            Array4<Real> G = Gfield->array(mfi);
+
+            // Loop over cells and update G
+            amrex::ParallelFor(tbx, tby, tbz,
+
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    Bx(i,j,k) += c2 * dt * T_Algo::DownwardDx(G, coefs_x, n_coefs_x, i, j, k);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    By(i,j,k) += c2 * dt * T_Algo::DownwardDy(G, coefs_y, n_coefs_y, i, j, k);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    Bz(i,j,k) += c2 * dt * T_Algo::DownwardDz(G, coefs_z, n_coefs_z, i, j, k);
+                }
+            );
+        }
+
         if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
         {
             amrex::Gpu::synchronize();