Add option for finite-order centering of currents (nodal to staggered) (#1763)

* Start adding centering of current

* Implement arbitrary order centering and split inputs

* No need to define a brand new interpolation function

* Update input file of hybrid CI tests

* Clean up

* Clean up more

* Fix bug and clean up

* Use current centering in two existing CI tests

* Update documentation

* Move Calls To UpdateCurrentNodalToStag Into SyncCurrent

* Add Doxygen For New Function UpdateCurrentNodalToStag

* Add Doxygen For New Functions Used For Stencil Coefficients

* Finite-Order Centering of Currents Not Implemented With MR

1 files changed, 142 insertions, 58 deletions

diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 92a2efb61..7417921ab 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -83,6 +83,8 @@ amrex::Vector<int> WarpX::field_boundary_hi(AMREX_SPACEDIM,0);
 amrex::Vector<int> WarpX::particle_boundary_lo(AMREX_SPACEDIM,0);
 amrex::Vector<int> WarpX::particle_boundary_hi(AMREX_SPACEDIM,0);
 
+bool WarpX::do_current_centering = false;
+
 int WarpX::n_rz_azimuthal_modes = 1;
 int WarpX::ncomps = 1;
 
@@ -90,11 +92,15 @@ long WarpX::nox = 1;
 long WarpX::noy = 1;
 long WarpX::noz = 1;
 
-// For momentum-conserving field gathering, order of interpolation from the
-// staggered positions to the grid nodes
-int WarpX::field_gathering_nox = 2;
-int WarpX::field_gathering_noy = 2;
-int WarpX::field_gathering_noz = 2;
+// Order of finite-order centering of fields (staggered to nodal)
+int WarpX::field_centering_nox = 2;
+int WarpX::field_centering_noy = 2;
+int WarpX::field_centering_noz = 2;
+
+// Order of finite-order centering of currents (nodal to staggered)
+int WarpX::current_centering_nox = 2;
+int WarpX::current_centering_noy = 2;
+int WarpX::current_centering_noz = 2;
 
 bool WarpX::use_fdtd_nci_corr = false;
 bool WarpX::galerkin_interpolation = true;
@@ -229,6 +235,11 @@ WarpX::WarpX ()
 
     current_store.resize(nlevs_max);
 
+    if (do_current_centering)
+    {
+        current_fp_nodal.resize(nlevs_max);
+    }
+
     F_cp.resize(nlevs_max);
     rho_cp.resize(nlevs_max);
     current_cp.resize(nlevs_max);
@@ -630,6 +641,21 @@ WarpX::ReadParameters ()
 
         // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1
         pp_warpx.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes);
+
+        // If true, the current is deposited on a nodal grid and then interpolated onto a Yee grid
+        pp_warpx.query("do_current_centering", do_current_centering);
+
+        // If do_nodal = 1, Maxwell's equations are solved on a nodal grid and
+        // the current should not be centered
+        if (do_nodal)
+        {
+            do_current_centering = false;
+        }
+
+        if ((maxLevel() > 0) && do_current_centering)
+        {
+            amrex::Abort("\nFinite-order centering of currents is not implemented with mesh refinement");
+        }
     }
 
     {
@@ -695,64 +721,50 @@ WarpX::ReadParameters ()
         pp_interpolation.query("noz", noz);
 
 #ifdef WARPX_USE_PSATD
+
         if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+
             // For momentum-conserving field gathering, read from input the order of
             // interpolation from the staggered positions to the grid nodes
-            if (field_gathering_algo == GatheringAlgo::MomentumConserving) {
-                pp_interpolation.query("field_gathering_nox", field_gathering_nox);
-                pp_interpolation.query("field_gathering_noy", field_gathering_noy);
-                pp_interpolation.query("field_gathering_noz", field_gathering_noz);
+            if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving) {
+                pp_interpolation.query("field_centering_nox", field_centering_nox);
+                pp_interpolation.query("field_centering_noy", field_centering_noy);
+                pp_interpolation.query("field_centering_noz", field_centering_noz);
+            }
+
+            // Read order of finite-order centering of currents (nodal to staggered)
+            if (WarpX::do_current_centering) {
+                pp_interpolation.query("current_centering_nox", current_centering_nox);
+                pp_interpolation.query("current_centering_noy", current_centering_noy);
+                pp_interpolation.query("current_centering_noz", current_centering_noz);
             }
 
-            if (maxLevel() > 0) {
+            if (maxLevel() > 0)
+            {
                 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-                    field_gathering_nox == 2 && field_gathering_noy == 2 && field_gathering_noz == 2,
-                    "High-order interpolation (order > 2) is not implemented with mesh refinement");
+                    field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2,
+                    "High-order centering of fields (order > 2) is not implemented with mesh refinement");
             }
 
-            // Host vectors for Fornberg stencil coefficients
-            amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_x;
-            amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_y;
-            amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_z;
-
-            host_Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(field_gathering_nox, false);
-            host_Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(field_gathering_noy, false);
-            host_Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(field_gathering_noz, false);
-
-            // Host vectors for ordered Fornberg stencil coefficients used for finite-order centering
-            host_centering_stencil_coeffs_x.resize(field_gathering_nox);
-            host_centering_stencil_coeffs_y.resize(field_gathering_noy);
-            host_centering_stencil_coeffs_z.resize(field_gathering_noz);
-
-            // Re-order Fornberg stencil coefficients
-            // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2)
-            ReorderFornbergCoefficients(host_centering_stencil_coeffs_x,
-                                        host_Fornberg_stencil_coeffs_x, field_gathering_nox);
-            ReorderFornbergCoefficients(host_centering_stencil_coeffs_y,
-                                        host_Fornberg_stencil_coeffs_y, field_gathering_noy);
-            ReorderFornbergCoefficients(host_centering_stencil_coeffs_z,
-                                        host_Fornberg_stencil_coeffs_z, field_gathering_noz);
-
-            // Device vectors for ordered Fornberg stencil coefficients used for finite-order centering
-            device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size());
-            amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
-                                  host_centering_stencil_coeffs_x.begin(),
-                                  host_centering_stencil_coeffs_x.end(),
-                                  device_centering_stencil_coeffs_x.begin());
-
-            device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size());
-            amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
-                                  host_centering_stencil_coeffs_y.begin(),
-                                  host_centering_stencil_coeffs_y.end(),
-                                  device_centering_stencil_coeffs_y.begin());
-
-            device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size());
-            amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
-                                  host_centering_stencil_coeffs_z.begin(),
-                                  host_centering_stencil_coeffs_z.end(),
-                                  device_centering_stencil_coeffs_z.begin());
-
-            amrex::Gpu::synchronize();
+            if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving)
+            {
+                AllocateCenteringCoefficients(device_field_centering_stencil_coeffs_x,
+                                              device_field_centering_stencil_coeffs_y,
+                                              device_field_centering_stencil_coeffs_z,
+                                              field_centering_nox,
+                                              field_centering_noy,
+                                              field_centering_noz);
+            }
+
+            if (WarpX::do_current_centering)
+            {
+                AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x,
+                                              device_current_centering_stencil_coeffs_y,
+                                              device_current_centering_stencil_coeffs_z,
+                                              current_centering_nox,
+                                              current_centering_noy,
+                                              current_centering_noz);
+            }
         }
 #endif
 
@@ -1031,6 +1043,11 @@ WarpX::ClearLevel (int lev)
 
         current_store[lev][i].reset();
 
+        if (do_current_centering)
+        {
+            current_fp_nodal[lev][i].reset();
+        }
+
         current_cp[lev][i].reset();
         Efield_cp [lev][i].reset();
         Bfield_cp [lev][i].reset();
@@ -1209,6 +1226,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     current_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ,tag("current_fp[y]"));
     current_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ,tag("current_fp[z]"));
 
+    if (do_current_centering)
+    {
+        amrex::BoxArray const& nodal_ba = amrex::convert(ba, amrex::IntVect::TheNodeVector());
+        current_fp_nodal[lev][0] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+        current_fp_nodal[lev][1] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+        current_fp_nodal[lev][2] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ);
+    }
+
     Bfield_avg_fp[lev][0] = std::make_unique<MultiFab>(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[x]"));
     Bfield_avg_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[y]"));
     Bfield_avg_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[z]"));
@@ -1768,9 +1793,10 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma
 #endif
 }
 
-void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coeffs,
-                                        amrex::Vector<amrex::Real>& unordered_coeffs,
-                                        const int order)
+#ifdef WARPX_USE_PSATD
+void WarpX::ReorderFornbergCoefficients (amrex::Vector<amrex::Real>& ordered_coeffs,
+                                         amrex::Vector<amrex::Real>& unordered_coeffs,
+                                         const int order)
 {
     const int n = order / 2;
     for (int i = 0; i < n; i++) {
@@ -1781,6 +1807,64 @@ void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coef
     }
 }
 
+void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_x,
+                                           amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_y,
+                                           amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z,
+                                           const int centering_nox,
+                                           const int centering_noy,
+                                           const int centering_noz)
+{
+    // Vectors of Fornberg stencil coefficients
+    amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_x;
+    amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_y;
+    amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_z;
+
+    // Host vectors of stencil coefficients used for finite-order centering
+    amrex::Vector<amrex::Real> host_centering_stencil_coeffs_x;
+    amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y;
+    amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z;
+
+    Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(centering_nox, false);
+    Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(centering_noy, false);
+    Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(centering_noz, false);
+
+    host_centering_stencil_coeffs_x.resize(centering_nox);
+    host_centering_stencil_coeffs_y.resize(centering_noy);
+    host_centering_stencil_coeffs_z.resize(centering_noz);
+
+    // Re-order Fornberg stencil coefficients:
+    // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2)
+    ReorderFornbergCoefficients(host_centering_stencil_coeffs_x,
+                                Fornberg_stencil_coeffs_x, centering_nox);
+    ReorderFornbergCoefficients(host_centering_stencil_coeffs_y,
+                                Fornberg_stencil_coeffs_y, centering_noy);
+    ReorderFornbergCoefficients(host_centering_stencil_coeffs_z,
+                                Fornberg_stencil_coeffs_z, centering_noz);
+
+    // Device vectors of stencil coefficients used for finite-order centering
+
+    device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                          host_centering_stencil_coeffs_x.begin(),
+                          host_centering_stencil_coeffs_x.end(),
+                          device_centering_stencil_coeffs_x.begin());
+
+    device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                          host_centering_stencil_coeffs_y.begin(),
+                          host_centering_stencil_coeffs_y.end(),
+                          device_centering_stencil_coeffs_y.begin());
+
+    device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                          host_centering_stencil_coeffs_z.begin(),
+                          host_centering_stencil_coeffs_z.end(),
+                          device_centering_stencil_coeffs_z.begin());
+
+    amrex::Gpu::synchronize();
+}
+#endif
+
 const iMultiFab*
 WarpX::CurrentBufferMasks (int lev)
 {