Fixing the computation of ECT Rho Field (#2711)

* Fixed computation of ECT Rho

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Adding a missing preprocessor directive

* Updated contributors list

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Fixed an abort message

* Fix in Source/Initialization/WarpXInitData.cpp

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Added some comments

* Made EvolveRhoCartesianECT not static

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

Co-authored-by: lgiacome <lorenzo.giacome@cern.ch>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

6 files changed, 170 insertions, 96 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
index 2455d8907..850fe5b9b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
+++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
@@ -8,6 +8,7 @@ target_sources(WarpX
     EvolveF.cpp
     EvolveFPML.cpp
     EvolveG.cpp
+    EvolveECTRho.cpp
     FiniteDifferenceSolver.cpp
     MacroscopicEvolveE.cpp
     ApplySilverMuellerBoundary.cpp
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index 9de1e8f8b..5e0a944b2 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -197,97 +197,6 @@ void FiniteDifferenceSolver::EvolveBCartesian (
     }
 }
 
-void FiniteDifferenceSolver::EvolveRhoCartesianECT (
-    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
-    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
-    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
-    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, const int lev ) {
-#ifdef AMREX_USE_EB
-
-#if !(defined(WARPX_DIM_3D) || defined(WARPX_DIM_XZ))
-    amrex::Abort("EvolveRhoCartesianECT: Embedded Boundaries are only implemented in 2D3V and 3D3V");
-#endif
-
-    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
-
-    // Loop through the grids, and over the tiles within each grid
-#ifdef AMREX_USE_OMP
-#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
-#endif
-    for (MFIter mfi(*ECTRhofield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) {
-            amrex::Gpu::synchronize();
-        }
-        Real wt = amrex::second();
-
-        // Extract field data for this grid/tile
-        Array4<Real> const &Ex = Efield[0]->array(mfi);
-        Array4<Real> const &Ey = Efield[1]->array(mfi);
-        Array4<Real> const &Ez = Efield[2]->array(mfi);
-        Array4<Real> const &Rhox = ECTRhofield[0]->array(mfi);
-        Array4<Real> const &Rhoy = ECTRhofield[1]->array(mfi);
-        Array4<Real> const &Rhoz = ECTRhofield[2]->array(mfi);
-        amrex::Array4<amrex::Real> const &lx = edge_lengths[0]->array(mfi);
-        amrex::Array4<amrex::Real> const &ly = edge_lengths[1]->array(mfi);
-        amrex::Array4<amrex::Real> const &lz = edge_lengths[2]->array(mfi);
-        amrex::Array4<amrex::Real> const &Sx = face_areas[0]->array(mfi);
-        amrex::Array4<amrex::Real> const &Sy = face_areas[1]->array(mfi);
-        amrex::Array4<amrex::Real> const &Sz = face_areas[2]->array(mfi);
-
-        // Extract tileboxes for which to loop
-        Box const &trhox = mfi.tilebox(ECTRhofield[0]->ixType().toIntVect());
-        Box const &trhoy = mfi.tilebox(ECTRhofield[1]->ixType().toIntVect());
-        Box const &trhoz = mfi.tilebox(ECTRhofield[2]->ixType().toIntVect());
-
-        amrex::ParallelFor(trhox, trhoy, trhoz,
-
-            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-                if (Sx(i, j, k) <= 0) return;
-
-#ifndef WARPX_DIM_XZ
-                Rhox(i, j, k) = (Ey(i, j, k) * ly(i, j, k) - Ey(i, j, k + 1) * ly(i, j, k + 1) +
-                    Ez(i, j + 1, k) * lz(i, j + 1, k) - Ez(i, j, k) * lz(i, j, k)) / Sx(i, j, k);
-#endif
-            },
-
-            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-                if (Sy(i, j, k) <= 0) return;
-
-#ifdef WARPX_DIM_XZ
-                Rhoy(i, j, k) = (Ez(i, j, k) * lz(i, j, k) - Ez(i + 1, j, k) * lz(i + 1, j, k) +
-                    Ex(i, j + 1, k) * lx(i, j + 1, k) - Ex(i, j, k) * lx(i, j, k)) / Sy(i, j, k);
-#elif defined(WARPX_DIM_3D)
-                Rhoy(i, j, k) = (Ez(i, j, k) * lz(i, j, k) - Ez(i + 1, j, k) * lz(i + 1, j, k) +
-                    Ex(i, j, k + 1) * lx(i, j, k + 1) - Ex(i, j, k) * lx(i, j, k)) / Sy(i, j, k);
-#else
-                amrex::Abort("EvolveRhoCartesianECT: Embedded Boundaries are only implemented in 2D3V and 3D3V");
-#endif
-            },
-
-            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-                if (Sz(i, j, k) <= 0) return;
-
-#ifndef WARPX_DIM_XZ
-                Rhoz(i, j, k) =  (Ex(i, j, k) * lx(i, j, k) - Ex(i, j + 1, k) * lx(i, j + 1, k) +
-                    Ey(i + 1, j, k) * ly(i + 1, j, k) - Ey(i, j, k) * ly(i, j, k)) / Sz(i, j, k);
-#endif
-            }
-        );
-
-        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
-        {
-            amrex::Gpu::synchronize();
-            wt = amrex::second() - wt;
-            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
-        }
-#ifdef WARPX_DIM_XZ
-        amrex::ignore_unused(Ey, Rhox, Rhoz, ly);
-#endif
-    }
-#else
-    amrex::ignore_unused(Efield, edge_lengths, face_areas, ECTRhofield, lev);
-#endif
-}
 
 void FiniteDifferenceSolver::EvolveBCartesianECT (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
index 318081159..b044f3b70 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -76,11 +76,7 @@ void FiniteDifferenceSolver::EvolveE (
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
         EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
-#ifdef AMREX_USE_EB
-        if (m_fdtd_algo == MaxwellSolverAlgo::ECT) {
-            EvolveRhoCartesianECT(Efield, edge_lengths, face_areas, ECTRhofield, lev);
-        }
-#endif
+
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
         EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveECTRho.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveECTRho.cpp
new file mode 100644
index 000000000..fd70c0552
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveECTRho.cpp
@@ -0,0 +1,161 @@
+/* Copyright 2020 Remi Lehe, Lorenzo Giacomel
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "FiniteDifferenceSolver.H"
+
+#ifndef WARPX_DIM_RZ
+#   include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+#   include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+#   include "FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#else
+#   include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#endif
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "Utils/WarpXConst.H"
+#include "WarpX.H"
+
+#include <AMReX.H>
+#include <AMReX_Array4.H>
+#include <AMReX_Config.H>
+#include <AMReX_Extension.H>
+#include <AMReX_GpuAtomic.H>
+#include <AMReX_GpuContainers.H>
+#include <AMReX_GpuControl.H>
+#include <AMReX_GpuDevice.H>
+#include <AMReX_GpuLaunch.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_IndexType.H>
+#include <AMReX_LayoutData.H>
+#include <AMReX_MFIter.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_iMultiFab.H>
+#include <AMReX_REAL.H>
+#include <AMReX_Utility.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <array>
+#include <memory>
+
+using namespace amrex;
+
+/**
+ * \brief Update the B field, over one timestep
+ */
+void FiniteDifferenceSolver::EvolveECTRho (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+    const int lev) {
+
+#if !defined(WARPX_DIM_RZ) and defined(AMREX_USE_EB)
+    if (m_fdtd_algo == MaxwellSolverAlgo::ECT) {
+
+        EvolveRhoCartesianECT(Efield, edge_lengths, face_areas, ECTRhofield, lev);
+
+    }
+#else
+    amrex::ignore_unused(Efield, edge_lengths, face_areas, ECTRhofield, lev);
+#endif
+}
+
+// If we implement ECT in 1D we will need to take care of this #ifndef differently
+#ifndef WARPX_DIM_RZ
+void FiniteDifferenceSolver::EvolveRhoCartesianECT (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, const int lev ) {
+#ifdef AMREX_USE_EB
+
+#if !(defined(WARPX_DIM_3D) || defined(WARPX_DIM_XZ))
+    amrex::Abort("EvolveRhoCartesianECT: Embedded Boundaries are only implemented in 3D and XZ");
+#endif
+
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
+    // Loop through the grids, and over the tiles within each grid
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+    for (amrex::MFIter mfi(*ECTRhofield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) {
+            amrex::Gpu::synchronize();
+        }
+        amrex::Real wt = amrex::second();
+
+        // Extract field data for this grid/tile
+        amrex::Array4<amrex::Real> const &Ex = Efield[0]->array(mfi);
+        amrex::Array4<amrex::Real> const &Ey = Efield[1]->array(mfi);
+        amrex::Array4<amrex::Real> const &Ez = Efield[2]->array(mfi);
+        amrex::Array4<amrex::Real> const &Rhox = ECTRhofield[0]->array(mfi);
+        amrex::Array4<amrex::Real> const &Rhoy = ECTRhofield[1]->array(mfi);
+        amrex::Array4<amrex::Real> const &Rhoz = ECTRhofield[2]->array(mfi);
+        amrex::Array4<amrex::Real> const &lx = edge_lengths[0]->array(mfi);
+        amrex::Array4<amrex::Real> const &ly = edge_lengths[1]->array(mfi);
+        amrex::Array4<amrex::Real> const &lz = edge_lengths[2]->array(mfi);
+        amrex::Array4<amrex::Real> const &Sx = face_areas[0]->array(mfi);
+        amrex::Array4<amrex::Real> const &Sy = face_areas[1]->array(mfi);
+        amrex::Array4<amrex::Real> const &Sz = face_areas[2]->array(mfi);
+
+        // Extract tileboxes for which to loop
+        amrex::Box const &trhox = mfi.tilebox(ECTRhofield[0]->ixType().toIntVect());
+        amrex::Box const &trhoy = mfi.tilebox(ECTRhofield[1]->ixType().toIntVect());
+        amrex::Box const &trhoz = mfi.tilebox(ECTRhofield[2]->ixType().toIntVect());
+
+        amrex::ParallelFor(trhox, trhoy, trhoz,
+
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (Sx(i, j, k) <= 0) return;
+
+// If we implement ECT in 1D we will need to take care of this #ifndef differently
+#ifndef WARPX_DIM_XZ
+                Rhox(i, j, k) = (Ey(i, j, k) * ly(i, j, k) - Ey(i, j, k + 1) * ly(i, j, k + 1) +
+                    Ez(i, j + 1, k) * lz(i, j + 1, k) - Ez(i, j, k) * lz(i, j, k)) / Sx(i, j, k);
+#endif
+            },
+
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (Sy(i, j, k) <= 0) return;
+
+#ifdef WARPX_DIM_XZ
+                Rhoy(i, j, k) = (Ez(i, j, k) * lz(i, j, k) - Ez(i + 1, j, k) * lz(i + 1, j, k) +
+                    Ex(i, j + 1, k) * lx(i, j + 1, k) - Ex(i, j, k) * lx(i, j, k)) / Sy(i, j, k);
+#elif defined(WARPX_DIM_3D)
+                Rhoy(i, j, k) = (Ez(i, j, k) * lz(i, j, k) - Ez(i + 1, j, k) * lz(i + 1, j, k) +
+                    Ex(i, j, k + 1) * lx(i, j, k + 1) - Ex(i, j, k) * lx(i, j, k)) / Sy(i, j, k);
+#else
+                amrex::Abort("EvolveRhoCartesianECT: Embedded Boundaries are only implemented in 3D and XZ");
+#endif
+            },
+
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (Sz(i, j, k) <= 0) return;
+
+// If we implement ECT in 1D we will need to take care of this #ifndef differently
+#ifndef WARPX_DIM_XZ
+                Rhoz(i, j, k) =  (Ex(i, j, k) * lx(i, j, k) - Ex(i, j + 1, k) * lx(i, j + 1, k) +
+                    Ey(i + 1, j, k) * ly(i + 1, j, k) - Ey(i, j, k) * ly(i, j, k)) / Sz(i, j, k);
+#endif
+            }
+        );
+
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
+#ifdef WARPX_DIM_XZ
+        amrex::ignore_unused(Ey, Rhox, Rhoz, ly);
+#endif
+    }
+#else
+    amrex::ignore_unused(Efield, edge_lengths, face_areas, ECTRhofield, lev);
+#endif
+}
+#endif
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 178f81e8c..5b3ef3930 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -76,6 +76,12 @@ class FiniteDifferenceSolver
                       std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
                       amrex::Real const dt);
 
+        void EvolveECTRho ( std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+                            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+                            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
+                            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+                            const int lev );
+
         void ApplySilverMuellerBoundary(
             std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
index e1ddee9ac..898b54a2f 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package
+++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
@@ -3,6 +3,7 @@ CEXE_sources += EvolveB.cpp
 CEXE_sources += EvolveE.cpp
 CEXE_sources += EvolveF.cpp
 CEXE_sources += EvolveG.cpp
+CEXE_sources += EvolveECTRho.cpp
 CEXE_sources += ComputeDivE.cpp
 CEXE_sources += MacroscopicEvolveE.cpp
 CEXE_sources += EvolveBPML.cpp