ECT conformal solver (#1923)

* adding base implementation of the conformal solver

* adding some preprocessor directives

* qualifying the isnan's correctly

* getting rid of some unused fields

* computing S_stab on the fly

* using empty in length check

* Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp

* replaced a looponcpu with a parallelfor

* trying to not pass lending/borrowing info by reference in evolveB

* passing data using dataPtr

* converting inds into a device vector

* simplifying some +=

* extracting the inds BaseFabs as DataPtr

* Revert "extracting the inds BaseFabs as DataPtr"

This reverts commit bc709d2fcaa347c119514de79a3f57169a05af65.

* implementing new way ogf handling cell extensions (gpu compatible)

* fixing some white spaces

* removed a typo

* made a function static

* tidying up

* tidying up

* getting rid of the rho vector

* moving the rho update to evolveE

* refectoring the cell extension

* small chenge in input parameters

* updating WarpX.H

* bug fix

* fixing another bug

* rotating the cube

* updated gitignore

* using the conformal soler in the cube test

* fixing the signature of a function

* trying to fix the setVal problem

* Revert "trying to fix the setVal problem"

This reverts commit c7d0e5e3709b730ff570183b2a6df5f587ca4640.

* trying to fix the setVal problem

* cleaning some comments

* removing an abort in device code

* Including geometric information in the external field initializer

* cleaned up the test

* moving the geometric data initialization before the external fields initialization

* changing way of saving info about intruded cells

* fixed some white spaces

* Using S_mod instead of S_red and S_enl

* substituting a std::vector with amrex::Array1D

* bug fix in the uint8_t -> coords maps

* Condensed ect cell info into one single MultiFab

* fixing some includes

* fixing some more includes

* fixed a typo

* making some functions available in gpu code

* using tilebox instead of validbox for the geometry initialization

* communicating the geometric info in the guard cells

* fixing the guard cells initialization for ect

* fixing an unused parameter

* fixing an unused parameter

* ignored some unused variables when EB is not supported

* cleaning up

* cleaning up

* ignored some unused variables when EB is not supported

* evolving rho just for ect

* handling some more unused variables

* removing some white spaces

* adding the rotated cube test

* removed some white spaces

* change cells into faces

* small fix in unused parameters

* fixed a comment

* for safety for now just initialize the geometric data when lev==maxLevel

* adding some preprocessor directives to isolate EB code

* updating Makefile stuff

* Improving the rotated cube test

* bug fix in mesh refinement

* fixing boundary conditions in rotated cube test

* deactivating MPI in rotated cube test

* activating mpi in test

* improving outputs of cells extension

* improved the documentation

* Update Source/EmbeddedBoundary/WarpXInitEB.cpp

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Update Source/EmbeddedBoundary/WarpXInitEB.cpp

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Removed some obsolete isnan's

* undone change to gitignore

* adding some missing AMREX_GPU_DEVICEs

* Changing some amrex::Real(0) into 0.

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

* Making some variables const

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

* Adding amrex:: to some ignore_unused

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

* Enhancing readability of some parallelfor's

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

* Removed some unecessary includes

* Fixing some tags

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

* Fixing the function CountExtFaces

* added a comment for Rhofield

* Fixed a typo in CountExtCells

* trying to remove accesses to private members in ComputeEdgeLengths

* making some functions public

* undoing some useless changes

* undoing some useless changes

* adding some AMREX_GPU_DEVICEs

* adding some unused variables

* adding some AMREX_GPU_DEVICEs to fix warnings

* fixing relative error

* Making several variables const

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

* Fixed a comment

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

* changing some double to amrex::Real

* removing commented lines

* changing some double to amrex::Real

* removing commented lines

* inizialing nelems_x,y,z to avoid a warning

* Fixing number of cells

* Adding missing analysis routine in rotated cube test

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

* Making some variables constexpr

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Improving some reduce operations

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* reading time directly from the output

* fixed a few data types

* fixing another int

* Replacing or->||, and->&&, not->!

* Adding license info

* Adding a white space

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

* removed some unused imports

* Moving the doxygen documentation

* Including AMReX_LayoutData.H

* removing some useless parameters

* adding ect solver to the doc

* Removing some useless reductions

* Small change for consistency

* Handling the resizing of arrays

* Handling correctly the resing of arrays

* Fixing some whitespaces

* Fixing an indentation

* Improving a comment

* Removed a useless initialization

* Renamed Rhofield to ECTRhofield

* Renamed Rhofield to ECTRhofield

* Added some if conditions to isolate some EB related code

* Improved some comments

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Made some functions not member of WarpX anymore

* Isolated some EB-only code

* Isolated some EB-only code

* Using the square brackets to access vectors

* Ignoring some unused variables

* Bug fix in initialization

* Removed a redundant if

* Using enumeration for connectivity and improved names

* Added a comment

* added a few comments

* Removed some useless template parameters

* Printing some info only in verbose mode

* Revert "Removed some useless template parameters"

This reverts commit 2c527980e6872c0212767c27f00e2b53ecbcfd0a.

* Fixed a bug with the neighbours enumeration

* Initializing geom data alsoo in the ghost cells

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

6 files changed, 282 insertions, 7 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index f5bba9dea..493e91dc9 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -31,6 +31,7 @@
 #include <AMReX_LayoutData.H>
 #include <AMReX_MFIter.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_iMultiFab.H>
 #include <AMReX_REAL.H>
 #include <AMReX_Utility.H>
 
@@ -49,8 +50,17 @@ void FiniteDifferenceSolver::EvolveB (
     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,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
+    std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
+    std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
     int lev, amrex::Real const dt ) {
 
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing);
+#endif
+
    // Select algorithm (The choice of algorithm is a runtime option,
    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
@@ -69,7 +79,12 @@ void FiniteDifferenceSolver::EvolveB (
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
         EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
+#ifdef AMREX_USE_EB
+    } else if (m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
+        EvolveBCartesianECT(Bfield, face_areas, area_mod, ECTRhofield, Venl, flag_info_cell,
+                            borrowing, lev, dt);
+#endif
 #endif
     } else {
         amrex::Abort("EvolveB: Unknown algorithm");
@@ -197,6 +212,224 @@ 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
+    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;
+
+                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);
+
+            },
+
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (Sy(i, j, k) <= 0) return;
+
+                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);
+
+            },
+
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (Sz(i, j, k) <= 0) return;
+
+                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);
+
+            }
+        );
+
+        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);
+        }
+    }
+#else
+    amrex::ignore_unused(Efield, edge_lengths, face_areas, ECTRhofield, lev);
+#endif
+}
+
+void FiniteDifferenceSolver::EvolveBCartesianECT (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
+    std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
+    std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
+    const int lev, amrex::Real const dt ) {
+#ifdef AMREX_USE_EB
+    amrex::LayoutData<amrex::Real> *cost = WarpX::getCosts(lev);
+
+    Venl[0]->setVal(0.);
+    Venl[1]->setVal(0.);
+    Venl[2]->setVal(0.);
+
+    // 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(*Bfield[0]); mfi.isValid(); ++mfi) {
+
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
+
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            // Extract field data for this grid/tile
+            Array4<Real> const &B = Bfield[idim]->array(mfi);
+            Array4<Real> const &Rho = ECTRhofield[idim]->array(mfi);
+            Array4<Real> const &Venl_dim = Venl[idim]->array(mfi);
+
+            amrex::Array4<int> const &flag_info_cell_dim = flag_info_cell[idim]->array(mfi);
+            amrex::Array4<Real> const &S = face_areas[idim]->array(mfi);
+            amrex::Array4<Real> const &S_mod = area_mod[idim]->array(mfi);
+
+            auto &borrowing_dim = (*borrowing[idim])[mfi];
+            auto borrowing_dim_neigh_faces = borrowing_dim.neigh_faces.data();
+            auto borrowing_dim_area = borrowing_dim.area.data();
+
+            auto const &borrowing_inds = (*borrowing[idim])[mfi].inds.data();
+            auto const &borrowing_size = (*borrowing[idim])[mfi].size.array();
+            auto const &borrowing_inds_pointer = (*borrowing[idim])[mfi].inds_pointer.array();
+
+            // Extract tileboxes for which to loop
+            Box const &tb = mfi.tilebox(Bfield[idim]->ixType().toIntVect());
+
+            //Take care of the unstable cells
+            amrex::ParallelFor(tb, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+
+                if (S(i, j, k) <= 0) return;
+
+                if (!(flag_info_cell_dim(i, j, k) == 0))
+                    return;
+
+                Venl_dim(i, j, k) = Rho(i, j, k) * S(i, j, k);
+                amrex::Real rho_enl;
+
+                // First we compute the rho of the enlarged face
+                for (int offset = 0; offset<borrowing_size(i, j, k); offset++) {
+                    int ind = borrowing_inds[*borrowing_inds_pointer(i, j, k) + offset];
+                    auto vec = FaceInfoBox::uint8_to_inds(borrowing_dim_neigh_faces[ind]);
+                    int ip, jp, kp;
+                    if(idim == 0){
+                        ip = i;
+                        jp = j + vec(0);
+                        kp = k + vec(1);
+                    }else if(idim == 1){
+                        ip = i + vec(0);
+                        jp = j;
+                        kp = k + vec(1);
+                    }else{
+                        ip = i + vec(0);
+                        jp = j + vec(1);
+                        kp = k;
+                    }
+
+                    Venl_dim(i, j, k) += Rho(ip, jp, kp) * borrowing_dim_area[ind];
+
+                }
+
+                rho_enl = Venl_dim(i, j, k) / S_mod(i, j, k);
+
+                for (int offset = 0; offset < borrowing_size(i, j, k); offset++) {
+                    int ind = borrowing_inds[*borrowing_inds_pointer(i, j, k) + offset];
+                    auto vec = FaceInfoBox::uint8_to_inds(borrowing_dim_neigh_faces[ind]);
+                    int ip, jp, kp;
+                    if(idim == 0){
+                        ip = i;
+                        jp = j + vec(0);
+                        kp = k + vec(1);
+                    }else if(idim == 1){
+                        ip = i + vec(0);
+                        jp = j;
+                        kp = k + vec(1);
+                    }else{
+                        ip = i + vec(0);
+                        jp = j + vec(1);
+                        kp = k;
+                    }
+
+                    Venl_dim(ip, jp, kp) += rho_enl * borrowing_dim_area[ind];
+
+                }
+
+                B(i, j, k) = B(i, j, k) - dt * rho_enl;
+
+            });
+
+            //Take care of the stable cells
+            amrex::ParallelFor(tb, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+                if (S(i, j, k) <= 0) return;
+
+                if (flag_info_cell_dim(i, j, k) == 0) {
+                    return;
+                }
+                else if (flag_info_cell_dim(i, j, k) == 1) {
+                    //Stable cell which hasn't been intruded
+                    B(i, j, k) = B(i, j, k) - dt * Rho(i, j, k);
+                } else if (flag_info_cell_dim(i, j, k) == 2) {
+                    //Stable cell which has been intruded
+                    Venl_dim(i, j, k) += Rho(i, j, k) * S_mod(i, j, k);
+                    B(i, j, k) = B(i, j, k) - dt * Venl_dim(i, j, k) / S(i, j, k);
+                }
+
+            });
+
+        }
+        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);
+        }
+    }
+#else
+    amrex::ignore_unused(Bfield, face_areas, area_mod, ECTRhofield, Venl, flag_info_cell, borrowing,
+                         lev, dt);
+#endif
+}
+
 #else // corresponds to ifndef WARPX_DIM_RZ
 
 template<typename T_Algo>
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
index 299f92f2f..5b182812a 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
@@ -55,7 +55,7 @@ void FiniteDifferenceSolver::EvolveBPML (
 
         EvolveBPMLCartesian <CartesianNodalAlgorithm> (Bfield, Efield, dt, dive_cleaning);
 
-    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
         EvolveBPMLCartesian <CartesianYeeAlgorithm> (Bfield, Efield, dt, dive_cleaning);
 
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
index 241992539..ac19984cd 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -49,11 +49,21 @@ void FiniteDifferenceSolver::EvolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     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,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
     int lev, amrex::Real const dt ) {
 
-   // Select algorithm (The choice of algorithm is a runtime option,
-   // but we compile code for each algorithm, using templates)
+#ifdef AMREX_USE_EB
+    if (m_fdtd_algo != MaxwellSolverAlgo::ECT) {
+        amrex::ignore_unused(face_areas, ECTRhofield);
+    }
+#else
+    amrex::ignore_unused(face_areas, ECTRhofield);
+#endif
+
+    // Select algorithm (The choice of algorithm is a runtime option,
+    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
     if (m_fdtd_algo == MaxwellSolverAlgo::Yee){
         ignore_unused(edge_lengths);
@@ -63,10 +73,14 @@ void FiniteDifferenceSolver::EvolveE (
 
         EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
 
-    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+    } 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 );
@@ -157,13 +171,16 @@ void FiniteDifferenceSolver::EvolveECartesian (
                 // Skip field push if this cell is fully covered by embedded boundaries
                 if (ly(i,j,k) <= 0) return;
 #endif
+
                 Ey(i, j, k) += c2 * dt * (
                     - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k)
                     + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k)
                     - PhysConst::mu0 * jy(i, j, k) );
+
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
 #ifdef AMREX_USE_EB
                 // Skip field push if this cell is fully covered by embedded boundaries
                 if (lz(i,j,k) <= 0) return;
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
index ea5ecb381..b182a836b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
@@ -61,7 +61,7 @@ void FiniteDifferenceSolver::EvolveEPML (
         EvolveEPMLCartesian <CartesianNodalAlgorithm> (
             Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
 
-    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+    } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
         EvolveEPMLCartesian <CartesianYeeAlgorithm> (
             Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 820f2b876..2a96ae44d 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -8,6 +8,7 @@
 #ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_
 #define WARPX_FINITE_DIFFERENCE_SOLVER_H_
 
+#include "EmbeddedBoundary/WarpXFaceInfoBox.H"
 #include "FiniteDifferenceSolver_fwd.H"
 
 #include "BoundaryConditions/PML_fwd.H"
@@ -49,12 +50,19 @@ class FiniteDifferenceSolver
                        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,
+                       std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
+                       std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+                       std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
+                       std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
+                       std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
                        int lev, amrex::Real const dt );
 
         void EvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
                        std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
                        std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
                        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,
                        std::unique_ptr<amrex::MultiFab> const& Ffield,
                        int lev, amrex::Real const dt );
 
@@ -200,6 +208,23 @@ class FiniteDifferenceSolver
             std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
             amrex::Real const dt);
 
+        void 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, int lev);
+
+        void EvolveBCartesianECT (
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl,
+            std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell,
+            std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
+            int lev, amrex::Real const dt
+        );
+
         template< typename T_Algo >
         void ComputeDivECartesian (
             const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
index aab4a69c5..cf6a603c7 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
@@ -65,7 +65,7 @@ FiniteDifferenceSolver::FiniteDifferenceSolver (
         CartesianNodalAlgorithm::InitializeStencilCoefficients( cell_size,
             m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z );
 
-    } else if (fdtd_algo == MaxwellSolverAlgo::Yee) {
+    } else if (fdtd_algo == MaxwellSolverAlgo::Yee || fdtd_algo == MaxwellSolverAlgo::ECT) {
 
         CartesianYeeAlgorithm::InitializeStencilCoefficients( cell_size,
             m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z );