Ohm's law solver (hybrid kinetic-fluid extension) (#3665)

* Add "None" as an option for the Maxwell solver

* fixed some of the reasons for failing CI tests

* no longer pass `do_electrostatic` to `GuardCellManager`

* renamed `MaxwellSolverAlgo` to `ElectromagneticSolverAlgo`

* rename `do_electrostatic` to `electrostatic_solver_id`

* rename `maxwell_solver_id` to `electromagnetic_solver_id`

* start of hybrid solver logic

* changes requested during PR review

* remove `do_no_deposit` from tests without field evolution

* added `HybridSolveE.cpp`

* bulk of the hybrid solver implementation

* mostly reproduce 1d cold ion mirror results

* ion Bernstein modes reproduced with this version

* fix bug with reduced diagnostic FieldProbe in 1d

* added hybrid solver installation to PICMI and added example script generating ion-Bernstein modes

* enable the use of `FieldProbe` default parameter values

* use default field-probe values

* steady progress

* add `do_not_push` flag to picmi

* possibly use nodal fields

* added extra multifab for current calculated from curl B

* added `CalculateTotalCurrent` functions

* rewrote implementation to calculate J x curl B on a nodal grid first

* fill boundary for auxiliary MFs used during hybrid E solve

* properly handle nonzero resistivity

* updated hybrid model example

* clean up example scripts

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

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

* fixed invalid memory access for GPU and other code cleanups

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

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

* refinements on the example scripts

* added ion beam instability example

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

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

* added EM modes and ion beam examples to CI test suite

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

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

* started docs section on the hybrid model

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

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

* more progress on documentation

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

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

* added ion Landau damping verification test / example

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

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

* add checksum benchmark for Landau damping example

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

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

* added fields.py wrapper to access total current density in hybrid case

* refactored the charge deposition fix to be performed with the field data rather than individual particles

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

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

* also correct current density at PEC boundary

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

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

* made resistivity a parsed function of `rho`

* work on PEC boundary condition

* corrections pointed out during code review

* fix build fails due to unused variables

* fix issue with GPU builds

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

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

* actually apply rho boundary correction in EM case

* take one sided derivative at PEC boundary when calculating div Pe

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

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

* actually apply rho boundary correction in EM case

* removed specific treatment of E-field on PEC boundary for Ohm's law solver

* first round of CI fixes

* second round of CI fixes

* added description of deposition logic with PEC boundaries to documentation

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

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

* third round of CI fixing

* move J and rho boundary handling to after `SyncRho` and `SyncCurrent` calls

* properly order the application of boundary conditions on rho, for electrostatic simulations

* fourth round of CI fixing

* moved calculation of total current (Ampere's law) to seperate function

* add random seed specification to `picmi`

* code clean-up -> renamed hybrid model to hybrid-PIC model

* added magnetic reconnection example

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

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

* code cleanup & benchmark updates

* update PICMI class name for hybrid solver to `HybridPICSolver`

* don't apply J field boundary in

* don't apply J field boundary in `MultiParticleContainer::DepositCurrent`

* apply changes requested during code review

* Apply suggestions from code review

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

* Loosen tolerance on failing CI test

* Removed unused variable

* code cleanup: make use of `MultiParticleContainer::DepositCurrent` in `AddSpaceChargeFieldLabFrame`

* switch to using a rho_fp_temp multifab for old and averaged charge density field, also no longer require particles to move only one cell per step

* use `ablastr::coarsen::sample` namespace in `HybridPICSolveECartesian`

* switched to using `MultiFab::LinComb` instead of self written GPU kernels to calculated averaged or extrapolated current density

* add verbosity flag for the Ohm solver tests

* deal with fine versus coarse patches

* add theoretical instability growth / damping rates to hybrid-PIC examples

* update ion-Bernstein mode plot in documentation

* move the `ApplyRhofieldBoundary` call to after `SumBoundary`

* use a uniform calculation for the number of cells a given index is from the boundary

* remove unused variable

* limit number of ghost cells updated during PEC BC application

* the number of ghost cells to consider depends on whether the field is nodal or not

* attempt 1 to fix failing CI tests

* attempt 2 to fix failing CI tests and code cleanup

* attempt 3 to fix failing CI tests

* attempt 4 to fix failing CI tests and docs cleanup

* switched to using bibtex citations

* move hybrid solver input parameters documentation to its own section

* clean up ion beam instability analysis script

* Apply suggestions from code review

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

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

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

* Apply suggestions from code review

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

* add inline comments describing the meaning of each argument for the `amrex::MultiFab::LinComb` calls used

* make `HybridPICSolver` a child class of `picmistandard.base._ClassWithInit`

* apply changes requested during code review

* add warning about using hybrid-PIC solver with Esirkepov current deposition

* add Stanier 2020 reference to recommend linear particles with hybrid-PIC

* add call to FillBoundary for `current_fp` - needed for accurate interpolation to nodal grid

* changes requested from code review

* Apply suggestions from code review

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

* include physics accuracy check for ion beam instability; switch CI tests to use direct current deposition

* reset benchmark values after switching to direct current deposition

* update ion beam instability benchmark

* minor changes requested during code review

* remove guard cells for `enE_nodal_mf` as well as corresponding `FillBoundary` call

* refactor: moved hybrid-PIC specific multifabs and `CalculateElectronPressure()` to `HybridPICModel`

* add assert that load balancing is not used with the hybrid-PIC solver since the new multifabs are not yet properly redistributed

* move `CalculateCurrentAmpere` to `HybridPICModel`

* refactor: moved `HybridPICSolveE` to `HybridPICModel` class

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 511 insertions, 0 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
new file mode 100644
index 000000000..24e17d040
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
@@ -0,0 +1,511 @@
+/* Copyright 2023 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Roelof Groenewald (TAE Technologies)
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "FiniteDifferenceSolver.H"
+
+#ifdef WARPX_DIM_RZ
+#   include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+#   include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+#endif
+#include "HybridPICModel/HybridPICModel.H"
+#include "Utils/TextMsg.H"
+#include "WarpX.H"
+
+#include <ablastr/coarsen/sample.H>
+
+using namespace amrex;
+
+void FiniteDifferenceSolver::CalculateCurrentAmpere (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3>& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev )
+{
+   // Select algorithm (The choice of algorithm is a runtime option,
+   // but we compile code for each algorithm, using templates)
+    if (m_fdtd_algo == ElectromagneticSolverAlgo::HybridPIC) {
+#ifdef WARPX_DIM_RZ
+        CalculateCurrentAmpereCylindrical <CylindricalYeeAlgorithm> (
+            Jfield, Bfield, edge_lengths, lev
+        );
+
+#else
+        CalculateCurrentAmpereCartesian <CartesianYeeAlgorithm> (
+            Jfield, Bfield, edge_lengths, lev
+        );
+
+#endif
+    } else {
+        amrex::Abort(Utils::TextMsg::Err(
+            "CalculateCurrentAmpere: Unknown algorithm choice."));
+    }
+}
+
+// /**
+//   * \brief Calculate electron current from Ampere's law without displacement
+//   * current and the kinetically tracked ion currents i.e. J_e = curl B.
+//   *
+//   * \param[out] Jfield  vector of electron current MultiFabs at a given level
+//   * \param[in] Bfield   vector of magnetic field MultiFabs at a given level
+//   */
+#ifdef WARPX_DIM_RZ
+template<typename T_Algo>
+void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev
+)
+{
+    amrex::ignore_unused(Jfield, Bfield, edge_lengths, lev);
+    amrex::Abort(Utils::TextMsg::Err(
+        "currently hybrid E-solve does not work for RZ"));
+}
+#else
+template<typename T_Algo>
+void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev
+)
+{
+    // for the profiler
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
+
+    // reset Jfield
+    Jfield[0]->setVal(0);
+    Jfield[1]->setVal(0);
+    Jfield[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(*Jfield[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& Jx = Jfield[0]->array(mfi);
+        Array4<Real> const& Jy = Jfield[1]->array(mfi);
+        Array4<Real> const& Jz = Jfield[2]->array(mfi);
+        Array4<Real const> const& Bx = Bfield[0]->const_array(mfi);
+        Array4<Real const> const& By = Bfield[1]->const_array(mfi);
+        Array4<Real const> const& Bz = Bfield[2]->const_array(mfi);
+
+#ifdef AMREX_USE_EB
+        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);
+#endif
+
+        // Extract stencil coefficients
+        Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
+        int const n_coefs_x = m_stencil_coefs_x.size();
+        Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
+        int const n_coefs_y = m_stencil_coefs_y.size();
+        Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+        int const n_coefs_z = m_stencil_coefs_z.size();
+
+        // Extract tileboxes for which to loop
+        Box const& tjx  = mfi.tilebox(Jfield[0]->ixType().toIntVect());
+        Box const& tjy  = mfi.tilebox(Jfield[1]->ixType().toIntVect());
+        Box const& tjz  = mfi.tilebox(Jfield[2]->ixType().toIntVect());
+
+        Real const one_over_mu0 = 1._rt / PhysConst::mu0;
+
+        // First calculate the total current using Ampere's law on the
+        // same grid as the E-field
+        amrex::ParallelFor(tjx, tjy, tjz,
+
+            // Jx calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip if this cell is fully covered by embedded boundaries
+                if (lx(i, j, k) <= 0) return;
+#endif
+                Jx(i, j, k) = one_over_mu0 * (
+                    - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k)
+                    + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k)
+                );
+            },
+
+            // Jy calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip if this cell is fully covered by embedded boundaries
+#ifdef WARPX_DIM_3D
+                if (ly(i,j,k) <= 0) return;
+#elif defined(WARPX_DIM_XZ)
+                // In XZ Jy is associated with a mesh node, so we need to check if the mesh node is covered
+                amrex::ignore_unused(ly);
+                if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0) return;
+#endif
+#endif
+                Jy(i, j, k) = one_over_mu0 * (
+                    - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k)
+                    + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k)
+                );
+            },
+
+            // Jz calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip if this cell is fully covered by embedded boundaries
+                if (lz(i,j,k) <= 0) return;
+#endif
+                Jz(i, j, k) = one_over_mu0 * (
+                    - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k)
+                    + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, 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);
+        }
+    }
+}
+#endif
+
+
+void FiniteDifferenceSolver::HybridPICSolveE (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jifield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
+    std::unique_ptr<amrex::MultiFab> const& rhofield,
+    std::unique_ptr<amrex::MultiFab> const& Pefield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev, HybridPICModel const* hybrid_model,
+    DtType a_dt_type )
+{
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        WarpX::grid_type == GridType::Staggered,
+        "Ohm's law E-solve only works with a staggered (Yee) grid.");
+
+
+   // Select algorithm (The choice of algorithm is a runtime option,
+   // but we compile code for each algorithm, using templates)
+    if (m_fdtd_algo == ElectromagneticSolverAlgo::HybridPIC) {
+#ifdef WARPX_DIM_RZ
+
+        HybridPICSolveECylindrical <CylindricalYeeAlgorithm> (
+            Efield, Jfield, Jifield, Bfield, rhofield, Pefield,
+            edge_lengths, lev, hybrid_model, a_dt_type
+        );
+
+#else
+
+        HybridPICSolveECartesian <CartesianYeeAlgorithm> (
+            Efield, Jfield, Jifield, Bfield, rhofield, Pefield,
+            edge_lengths, lev, hybrid_model, a_dt_type
+        );
+
+#endif
+    } else {
+        amrex::Abort(Utils::TextMsg::Err(
+            "HybridSolveE: The hybrid-PIC electromagnetic solver algorithm must be used"));
+    }
+}
+
+#ifdef WARPX_DIM_RZ
+template<typename T_Algo>
+void FiniteDifferenceSolver::HybridPICSolveECylindrical (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jifield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
+    std::unique_ptr<amrex::MultiFab> const& rhofield,
+    std::unique_ptr<amrex::MultiFab> const& Pefield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev, HybridPICModel const* hybrid_model,
+    DtType a_dt_type )
+{
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
+    amrex::ignore_unused(
+        Efield, Jfield, Jifield, Bfield, rhofield, Pefield, edge_lengths,
+        lev, hybrid_model, a_dt_type
+    );
+    amrex::Abort(Utils::TextMsg::Err(
+        "currently hybrid E-solve does not work for RZ"));
+}
+
+#else
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::HybridPICSolveECartesian (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jifield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
+    std::unique_ptr<amrex::MultiFab> const& rhofield,
+    std::unique_ptr<amrex::MultiFab> const& Pefield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    int lev, HybridPICModel const* hybrid_model,
+    DtType a_dt_type )
+{
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
+
+    // for the profiler
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
+    using namespace ablastr::coarsen::sample;
+
+    // get hybrid model parameters
+    const auto eta = hybrid_model->m_eta;
+    const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e;
+
+    // Index type required for interpolating fields from their respective
+    // staggering to the Ex, Ey, Ez locations
+    amrex::GpuArray<int, 3> const& Ex_stag = hybrid_model->Ex_IndexType;
+    amrex::GpuArray<int, 3> const& Ey_stag = hybrid_model->Ey_IndexType;
+    amrex::GpuArray<int, 3> const& Ez_stag = hybrid_model->Ez_IndexType;
+    amrex::GpuArray<int, 3> const& Jx_stag = hybrid_model->Jx_IndexType;
+    amrex::GpuArray<int, 3> const& Jy_stag = hybrid_model->Jy_IndexType;
+    amrex::GpuArray<int, 3> const& Jz_stag = hybrid_model->Jz_IndexType;
+    amrex::GpuArray<int, 3> const& Bx_stag = hybrid_model->Bx_IndexType;
+    amrex::GpuArray<int, 3> const& By_stag = hybrid_model->By_IndexType;
+    amrex::GpuArray<int, 3> const& Bz_stag = hybrid_model->Bz_IndexType;
+
+    // Parameters for `interp` that maps from Yee to nodal mesh and back
+    amrex::GpuArray<int, 3> const& nodal = {1, 1, 1};
+    // The "coarsening is just 1 i.e. no coarsening"
+    amrex::GpuArray<int, 3> const& coarsen = {1, 1, 1};
+
+    // The E-field calculation is done in 2 steps:
+    // 1) The J x B term is calculated on a nodal mesh in order to ensure
+    //    energy conservation.
+    // 2) The nodal E-field values are averaged onto the Yee grid and the
+    //    electron pressure & resistivity terms are added (these terms are
+    //    naturally located on the Yee grid).
+
+    // Create a temporary multifab to hold the nodal E-field values
+    // Note the multifab has 3 values for Ex, Ey and Ez which we can do here
+    // since all three components will be calculated on the same grid.
+    // Also note that enE_nodal_mf does not need to have any guard cells since
+    // these values will be interpolated to the Yee mesh which is contained
+    // by the nodal mesh.
+    auto const& ba = convert(rhofield->boxArray(), IntVect::TheNodeVector());
+    MultiFab enE_nodal_mf(ba, rhofield->DistributionMap(), 3, IntVect::TheZeroVector());
+
+    // Loop through the grids, and over the tiles within each grid for the
+    // initial, nodal calculation of E
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+    for ( MFIter mfi(enE_nodal_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
+
+        Array4<Real> const& enE_nodal = enE_nodal_mf.array(mfi);
+        Array4<Real const> const& Jx = Jfield[0]->const_array(mfi);
+        Array4<Real const> const& Jy = Jfield[1]->const_array(mfi);
+        Array4<Real const> const& Jz = Jfield[2]->const_array(mfi);
+        Array4<Real const> const& Jix = Jifield[0]->const_array(mfi);
+        Array4<Real const> const& Jiy = Jifield[1]->const_array(mfi);
+        Array4<Real const> const& Jiz = Jifield[2]->const_array(mfi);
+        Array4<Real const> const& Bx = Bfield[0]->const_array(mfi);
+        Array4<Real const> const& By = Bfield[1]->const_array(mfi);
+        Array4<Real const> const& Bz = Bfield[2]->const_array(mfi);
+
+        // Loop over the cells and update the nodal E field
+        amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+            // interpolate the total current to a nodal grid
+            auto const jx_interp = Interp(Jx, Jx_stag, nodal, coarsen, i, j, k, 0);
+            auto const jy_interp = Interp(Jy, Jy_stag, nodal, coarsen, i, j, k, 0);
+            auto const jz_interp = Interp(Jz, Jz_stag, nodal, coarsen, i, j, k, 0);
+
+            // interpolate the ion current to a nodal grid
+            auto const jix_interp = Interp(Jix, Jx_stag, nodal, coarsen, i, j, k, 0);
+            auto const jiy_interp = Interp(Jiy, Jy_stag, nodal, coarsen, i, j, k, 0);
+            auto const jiz_interp = Interp(Jiz, Jz_stag, nodal, coarsen, i, j, k, 0);
+
+            // interpolate the B field to a nodal grid
+            auto const Bx_interp = Interp(Bx, Bx_stag, nodal, coarsen, i, j, k, 0);
+            auto const By_interp = Interp(By, By_stag, nodal, coarsen, i, j, k, 0);
+            auto const Bz_interp = Interp(Bz, Bz_stag, nodal, coarsen, i, j, k, 0);
+
+            // calculate enE = (J - Ji) x B
+            enE_nodal(i, j, k, 0) = (
+                (jy_interp - jiy_interp) * Bz_interp
+                - (jz_interp - jiz_interp) * By_interp
+            );
+            enE_nodal(i, j, k, 1) = (
+                (jz_interp - jiz_interp) * Bx_interp
+                - (jx_interp - jix_interp) * Bz_interp
+            );
+            enE_nodal(i, j, k, 2) = (
+                (jx_interp - jix_interp) * By_interp
+                - (jy_interp - jiy_interp) * Bx_interp
+            );
+        });
+
+        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);
+        }
+    }
+
+    // Loop through the grids, and over the tiles within each grid again
+    // for the Yee grid calculation of the E field
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+    for ( MFIter mfi(*Efield[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> const& Jx = Jfield[0]->const_array(mfi);
+        Array4<Real const> const& Jy = Jfield[1]->const_array(mfi);
+        Array4<Real const> const& Jz = Jfield[2]->const_array(mfi);
+        Array4<Real const> const& enE = enE_nodal_mf.const_array(mfi);
+        Array4<Real const> const& rho = rhofield->const_array(mfi);
+        Array4<Real> const& Pe = Pefield->array(mfi);
+
+#ifdef AMREX_USE_EB
+        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);
+#endif
+
+        // Extract stencil coefficients
+        Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
+        int const n_coefs_x = m_stencil_coefs_x.size();
+        Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
+        int const n_coefs_y = m_stencil_coefs_y.size();
+        Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+        int const n_coefs_z = m_stencil_coefs_z.size();
+
+        Box const& tex  = mfi.tilebox(Efield[0]->ixType().toIntVect());
+        Box const& tey  = mfi.tilebox(Efield[1]->ixType().toIntVect());
+        Box const& tez  = mfi.tilebox(Efield[2]->ixType().toIntVect());
+
+        // Loop over the cells and update the E field
+        amrex::ParallelFor(tex, tey, tez,
+
+            // Ex calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip if this cell is fully covered by embedded boundaries
+                if (lx(i, j, k) <= 0) return;
+#endif
+                // Interpolate to get the appropriate charge density in space
+                Real rho_val = Interp(rho, nodal, Ex_stag, coarsen, i, j, k, 0);
+
+                // safety condition since we divide by rho_val later
+                if (rho_val < rho_floor) rho_val = rho_floor;
+
+                // Get the gradient of the electron pressure
+                auto grad_Pe = T_Algo::UpwardDx(Pe, coefs_x, n_coefs_x, i, j, k);
+
+                // interpolate the nodal neE values to the Yee grid
+                auto enE_x = Interp(enE, nodal, Ex_stag, coarsen, i, j, k, 0);
+
+                Ex(i, j, k) = (enE_x - grad_Pe) / rho_val;
+
+                // Add resistivity only if E field value is used to update B
+                if (a_dt_type != DtType::Full) Ex(i, j, k) += eta(rho_val) * Jx(i, j, k);
+            },
+
+            // Ey calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip field solve if this cell is fully covered by embedded boundaries
+#ifdef WARPX_DIM_3D
+                if (ly(i,j,k) <= 0) return;
+#elif defined(WARPX_DIM_XZ)
+                //In XZ Ey is associated with a mesh node, so we need to check if the mesh node is covered
+                amrex::ignore_unused(ly);
+                if (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0) return;
+#endif
+#endif
+                // Interpolate to get the appropriate charge density in space
+                Real rho_val = Interp(rho, nodal, Ey_stag, coarsen, i, j, k, 0);
+
+                // safety condition since we divide by rho_val later
+                if (rho_val < rho_floor) rho_val = rho_floor;
+
+                // Get the gradient of the electron pressure
+                auto grad_Pe = T_Algo::UpwardDy(Pe, coefs_y, n_coefs_y, i, j, k);
+
+                // interpolate the nodal neE values to the Yee grid
+                auto enE_y = Interp(enE, nodal, Ey_stag, coarsen, i, j, k, 1);
+
+                Ey(i, j, k) = (enE_y - grad_Pe) / rho_val;
+
+                // Add resistivity only if E field value is used to update B
+                if (a_dt_type != DtType::Full) Ey(i, j, k) += eta(rho_val) * Jy(i, j, k);
+            },
+
+            // Ez calculation
+            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+#ifdef AMREX_USE_EB
+                // Skip field solve if this cell is fully covered by embedded boundaries
+                if (lz(i,j,k) <= 0) return;
+#endif
+                // Interpolate to get the appropriate charge density in space
+                Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);
+
+                // safety condition since we divide by rho_val later
+                if (rho_val < rho_floor) rho_val = rho_floor;
+
+                // Get the gradient of the electron pressure
+                auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, k);
+
+                // interpolate the nodal neE values to the Yee grid
+                auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, k, 2);
+
+                Ez(i, j, k) = (enE_z - grad_Pe) / rho_val;
+
+                // Add resistivity only if E field value is used to update B
+                if (a_dt_type != DtType::Full) Ez(i, j, k) += eta(rho_val) * Jz(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);
+        }
+    }
+}
+#endif