Use only plasma current in HybridPICSolveE (ECP-WarpX#5273)

* use plasma current rather than total current in `HybridPICSolveE` * remove logic to subtract J_ext from plasma current in `JdispFunctor` * add one ghost cell to the hybrid-pic external current since we interpolate to a nodal grid * Fix Doxygen Signed-off-by: roelof-groenewald <[email protected]>
dpgrote · Oct 1, 2024 · 2d61720 · 2d61720
1 parent 617d7ba
commit 2d61720
Show file tree

Hide file tree

Showing 9 changed files with 92 additions and 141 deletions.
diff --git a/...ohm_solver_magnetic_reconnection/inputs_test_2d_ohm_solver_magnetic_reconnection_picmi.py b/...ohm_solver_magnetic_reconnection/inputs_test_2d_ohm_solver_magnetic_reconnection_picmi.py
@@ -303,7 +303,7 @@ def check_fields(self):
 
         rho = fields.RhoFPWrapper(include_ghosts=False)[:, :]
         Jiy = fields.JyFPWrapper(include_ghosts=False)[...] / self.J0
-        Jy = fields.JyFPAmpereWrapper(include_ghosts=False)[...] / self.J0
+        Jy = fields.JyFPPlasmaWrapper(include_ghosts=False)[...] / self.J0
         Bx = fields.BxFPWrapper(include_ghosts=False)[...] / self.B0
         By = fields.ByFPWrapper(include_ghosts=False)[...] / self.B0
         Bz = fields.BzFPWrapper(include_ghosts=False)[...] / self.B0

diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py
@@ -33,7 +33,7 @@
 ExFPPMLWrapper, EyFPPMLWrapper, EzFPPMLWrapper
 BxFPPMLWrapper, ByFPPMLWrapper, BzFPPMLWrapper
 JxFPPMLWrapper, JyFPPMLWrapper, JzFPPMLWrapper
-JxFPAmpereWrapper, JyFPAmpereWrapper, JzFPAmpereWrapper
+JxFPPlasmaWrapper, JyFPPlasmaWrapper, JzFPPlasmaWrapper
 FFPPMLWrapper, GFPPMLWrapper
 
 ExCPPMLWrapper, EyCPPMLWrapper, EzCPPMLWrapper
@@ -873,27 +873,27 @@ def FaceAreaszWrapper(level=0, include_ghosts=False):
     )
 
 
-def JxFPAmpereWrapper(level=0, include_ghosts=False):
+def JxFPPlasmaWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="hybrid_current_fp_ampere",
+        mf_name="hybrid_current_fp_plasma",
         idir=0,
         level=level,
         include_ghosts=include_ghosts,
     )
 
 
-def JyFPAmpereWrapper(level=0, include_ghosts=False):
+def JyFPPlasmaWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="hybrid_current_fp_ampere",
+        mf_name="hybrid_current_fp_plasma",
         idir=1,
         level=level,
         include_ghosts=include_ghosts,
     )
 
 
-def JzFPAmpereWrapper(level=0, include_ghosts=False):
+def JzFPPlasmaWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="hybrid_current_fp_ampere",
+        mf_name="hybrid_current_fp_plasma",
         idir=2,
         level=level,
         include_ghosts=include_ghosts,

diff --git a/Source/Diagnostics/ComputeDiagFunctors/JdispFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/JdispFunctor.cpp
@@ -1,8 +1,11 @@
-/* This file is part of Warpx.
+/* Copyright 2023-2024 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Avigdor Veksler (TAE Technologies)
  *
- * Authors: Avigdor Veksler
  * License: BSD-3-Clause-LBNL
-*/
+ */
 #include "JdispFunctor.H"
 
 #include "WarpX.H"
@@ -40,7 +43,7 @@ JdispFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buff
     AMREX_ASSUME(hybrid_pic_model != nullptr);
 
      /** pointer to current calculated from Ampere's Law (Jamp) multifab */
-    amrex::MultiFab* mf_curlB = warpx.m_fields.get(FieldType::hybrid_current_fp_ampere, Direction{m_dir}, m_lev);
+    amrex::MultiFab* mf_curlB = warpx.m_fields.get(FieldType::hybrid_current_fp_plasma, Direction{m_dir}, m_lev);
 
     //if (!hybrid_pic_model) {
         // To finish this implementation, we need to implement a method to
@@ -63,51 +66,6 @@ JdispFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buff
         -1, *mf_j, 0, 0, 1, Jdisp.nGrowVect()
     );
 
-    if (hybrid_pic_model) {
-        // Subtract the interpolated j_external value from j_displacement.
-        /** pointer to external currents (Jext) multifab */
-        amrex::MultiFab* mf_j_external = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{m_dir}, m_lev);
-
-        // Index type required for interpolating Jext from their respective
-        // staggering (nodal) to the Jx_displacement, Jy_displacement, Jz_displacement
-        // locations. The staggering of J_displacement is the same as the
-        // staggering for J, so we use J_stag as the interpolation map.
-        // For interp to work below, the indices of the undefined dimensions
-        // must match. We set them as (1,1,1).
-        amrex::GpuArray<int, 3> Jext_IndexType = {1, 1, 1};
-        amrex::GpuArray<int, 3> J_IndexType = {1, 1, 1};
-        amrex::IntVect Jext_stag = mf_j_external->ixType().toIntVect();
-        amrex::IntVect J_stag = mf_j->ixType().toIntVect();
-
-        // Index types for the dimensions simulated are overwritten.
-        for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-            Jext_IndexType[idim] = Jext_stag[idim];
-            J_IndexType[idim] = J_stag[idim];
-        }
-
-        // Parameters for `interp` that maps from Jext to J.
-        // The "coarsening is just 1 i.e. no coarsening"
-        amrex::GpuArray<int, 3> const& coarsen = {1, 1, 1};
-
-        // Loop through the grids, and over the tiles within each grid to
-        // subtract the interpolated Jext from J_displacement.
-#ifdef AMREX_USE_OMP
-#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
-#endif
-        for ( MFIter mfi(Jdisp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-
-            Array4<Real> const& Jdisp_arr = Jdisp.array(mfi);
-            Array4<Real const> const& Jext = mf_j_external->const_array(mfi);
-
-            //  Loop over cells and update the Jdisp MultiFab
-            amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                // Interpolate Jext to the staggering of J
-                auto const jext_interp = ablastr::coarsen::sample::Interp(Jext, Jext_IndexType, J_IndexType, coarsen, i, j, k, 0);
-                Jdisp_arr(i, j, k, 0) -= jext_interp;
-            });
-        }
-    }
-
     InterpolateMFForDiag(mf_dst, Jdisp, dcomp, warpx.DistributionMap(m_lev),
                          m_convertRZmodes2cartesian);
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -141,9 +141,8 @@ class FiniteDifferenceSolver
           * https://link.springer.com/chapter/10.1007/3-540-36530-3_8
           *
           * \param[out] Efield  vector of electric field MultiFabs updated at a given level
-          * \param[in] Jfield   vector of total current MultiFabs at a given level
+          * \param[in] Jfield   vector of total plasma current MultiFabs at a given level
           * \param[in] Jifield  vector of ion current density MultiFabs at a given level
-          * \param[in] Jextfield  vector of external current density MultiFabs at a given level
           * \param[in] Bfield   vector of magnetic field MultiFabs at a given level
           * \param[in] rhofield scalar ion charge density Multifab at a given level
           * \param[in] Pefield  scalar electron pressure MultiFab at a given level
@@ -153,15 +152,14 @@ class FiniteDifferenceSolver
           * \param[in] solve_for_Faraday boolean flag for whether the E-field is solved to be used in Faraday's equation
           */
         void HybridPICSolveE ( ablastr::fields::VectorField const& Efield,
-                      ablastr::fields::VectorField & Jfield,
-                      ablastr::fields::VectorField const& Jifield,
-                      ablastr::fields::VectorField const& Jextfield,
-                      ablastr::fields::VectorField const& Bfield,
-                      amrex::MultiFab const& rhofield,
-                      amrex::MultiFab const& Pefield,
-                      ablastr::fields::VectorField const& edge_lengths,
-                      int lev, HybridPICModel const* hybrid_model,
-                      bool solve_for_Faraday );
+                               ablastr::fields::VectorField & Jfield,
+                               ablastr::fields::VectorField const& Jifield,
+                               ablastr::fields::VectorField const& Bfield,
+                               amrex::MultiFab const& rhofield,
+                               amrex::MultiFab const& Pefield,
+                               ablastr::fields::VectorField const& edge_lengths,
+                               int lev, HybridPICModel const* hybrid_model,
+                               bool solve_for_Faraday );
 
         /**
           * \brief Calculation of total current using Ampere's law (without
@@ -241,7 +239,6 @@ class FiniteDifferenceSolver
             ablastr::fields::VectorField const& Efield,
             ablastr::fields::VectorField const& Jfield,
             ablastr::fields::VectorField const& Jifield,
-            ablastr::fields::VectorField const& Jextfield,
             ablastr::fields::VectorField const& Bfield,
             amrex::MultiFab const& rhofield,
             amrex::MultiFab const& Pefield,
@@ -346,7 +343,6 @@ class FiniteDifferenceSolver
             ablastr::fields::VectorField const& Efield,
             ablastr::fields::VectorField const& Jfield,
             ablastr::fields::VectorField const& Jifield,
-            ablastr::fields::VectorField const& Jextfield,
             ablastr::fields::VectorField const& Bfield,
             amrex::MultiFab const& rhofield,
             amrex::MultiFab const& Pefield,

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
@@ -63,18 +63,18 @@ public:
 
     /**
      * \brief
-     * Function to calculate the total current based on Ampere's law while
-     * neglecting displacement current (J = curl x B). Used in the Ohm's law
-     * solver (kinetic-fluid hybrid model).
+     * Function to calculate the total plasma current based on Ampere's law while
+     * neglecting displacement current (J = curl x B). Any external current is
+     * subtracted as well. Used in the Ohm's law solver (kinetic-fluid hybrid model).
      *
      * \param[in] Bfield       Magnetic field from which the current is calculated.
      * \param[in] edge_lengths Length of cell edges taking embedded boundaries into account
      */
-    void CalculateCurrentAmpere (
+    void CalculatePlasmaCurrent (
         ablastr::fields::MultiLevelVectorField const& Bfield,
         ablastr::fields::MultiLevelVectorField const& edge_lengths
     );
-    void CalculateCurrentAmpere (
+    void CalculatePlasmaCurrent (
         ablastr::fields::VectorField const& Bfield,
         ablastr::fields::VectorField const& edge_lengths,
         int lev

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
@@ -67,18 +67,18 @@ void HybridPICModel::AllocateLevelMFs (ablastr::fields::MultiFabRegister & field
 
     // The "hybrid_electron_pressure_fp" multifab stores the electron pressure calculated
     // from the specified equation of state.
-    // The "hybrid_rho_fp_temp" multifab is used to store the ion charge density
-    // interpolated or extrapolated to appropriate timesteps.
-    // The "hybrid_current_fp_temp" multifab is used to store the ion current density
-    // interpolated or extrapolated to appropriate timesteps.
-    // The "hybrid_current_fp_ampere" multifab stores the total current calculated as
-    // the curl of B.
     fields.alloc_init(FieldType::hybrid_electron_pressure_fp,
         lev, amrex::convert(ba, rho_nodal_flag),
         dm, ncomps, ngRho, 0.0_rt);
+
+    // The "hybrid_rho_fp_temp" multifab is used to store the ion charge density
+    // interpolated or extrapolated to appropriate timesteps.
     fields.alloc_init(FieldType::hybrid_rho_fp_temp,
         lev, amrex::convert(ba, rho_nodal_flag),
         dm, ncomps, ngRho, 0.0_rt);
+
+    // The "hybrid_current_fp_temp" multifab is used to store the ion current density
+    // interpolated or extrapolated to appropriate timesteps.
     fields.alloc_init(FieldType::hybrid_current_fp_temp, Direction{0},
         lev, amrex::convert(ba, jx_nodal_flag),
         dm, ncomps, ngJ, 0.0_rt);
@@ -89,28 +89,29 @@ void HybridPICModel::AllocateLevelMFs (ablastr::fields::MultiFabRegister & field
         lev, amrex::convert(ba, jz_nodal_flag),
         dm, ncomps, ngJ, 0.0_rt);
 
-    fields.alloc_init(FieldType::hybrid_current_fp_ampere, Direction{0},
+    // The "hybrid_current_fp_plasma" multifab stores the total plasma current calculated
+    // as the curl of B minus any external current.
+    fields.alloc_init(FieldType::hybrid_current_fp_plasma, Direction{0},
         lev, amrex::convert(ba, jx_nodal_flag),
         dm, ncomps, ngJ, 0.0_rt);
-    fields.alloc_init(FieldType::hybrid_current_fp_ampere, Direction{1},
+    fields.alloc_init(FieldType::hybrid_current_fp_plasma, Direction{1},
         lev, amrex::convert(ba, jy_nodal_flag),
         dm, ncomps, ngJ, 0.0_rt);
-    fields.alloc_init(FieldType::hybrid_current_fp_ampere, Direction{2},
+    fields.alloc_init(FieldType::hybrid_current_fp_plasma, Direction{2},
         lev, amrex::convert(ba, jz_nodal_flag),
         dm, ncomps, ngJ, 0.0_rt);
 
-    // the external current density multifab is made nodal to avoid needing to interpolate
-    // to a nodal grid as has to be done for the ion and total current density multifabs
-    // this also allows the external current multifab to not have any ghost cells
+    // the external current density multifab matches the current staggering and
+    // one ghost cell is used since we interpolate the current to a nodal grid
     fields.alloc_init(FieldType::hybrid_current_fp_external, Direction{0},
-        lev, amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), 0.0_rt);
+        lev, amrex::convert(ba, jx_nodal_flag),
+        dm, ncomps, IntVect(1), 0.0_rt);
     fields.alloc_init(FieldType::hybrid_current_fp_external, Direction{1},
-        lev, amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), 0.0_rt);
+        lev, amrex::convert(ba, jy_nodal_flag),
+        dm, ncomps, IntVect(1), 0.0_rt);
     fields.alloc_init(FieldType::hybrid_current_fp_external, Direction{2},
-        lev, amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), 0.0_rt);
+        lev, amrex::convert(ba, jz_nodal_flag),
+        dm, ncomps, IntVect(1), 0.0_rt);
 
 #ifdef WARPX_DIM_RZ
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
@@ -352,7 +353,7 @@ void HybridPICModel::GetCurrentExternal (
                 const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
 #endif
                 // Initialize the y-component of the field.
-                mfyfab(i,j,k)  = Jy_external(x,y,z,t);
+                mfyfab(i,j,k) = Jy_external(x,y,z,t);
             },
             [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                 // skip if node is covered by an embedded boundary
@@ -384,35 +385,44 @@ void HybridPICModel::GetCurrentExternal (
     }
 }
 
-void HybridPICModel::CalculateCurrentAmpere (
+void HybridPICModel::CalculatePlasmaCurrent (
     ablastr::fields::MultiLevelVectorField const& Bfield,
     ablastr::fields::MultiLevelVectorField const& edge_lengths)
 {
     auto& warpx = WarpX::GetInstance();
     for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
     {
-        CalculateCurrentAmpere(Bfield[lev], edge_lengths[lev], lev);
+        CalculatePlasmaCurrent(Bfield[lev], edge_lengths[lev], lev);
     }
 }
 
-void HybridPICModel::CalculateCurrentAmpere (
+void HybridPICModel::CalculatePlasmaCurrent (
     ablastr::fields::VectorField const& Bfield,
     ablastr::fields::VectorField const& edge_lengths,
     const int lev)
 {
-    WARPX_PROFILE("WarpX::CalculateCurrentAmpere()");
+    WARPX_PROFILE("HybridPICModel::CalculatePlasmaCurrent()");
 
     auto& warpx = WarpX::GetInstance();
-    ablastr::fields::VectorField current_fp_ampere = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_ampere, lev);
+    ablastr::fields::VectorField current_fp_plasma = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_plasma, lev);
     warpx.get_pointer_fdtd_solver_fp(lev)->CalculateCurrentAmpere(
-        current_fp_ampere, Bfield, edge_lengths, lev
+        current_fp_plasma, Bfield, edge_lengths, lev
     );
 
     // we shouldn't apply the boundary condition to J since J = J_i - J_e but
     // the boundary correction was already applied to J_i and the B-field
     // boundary ensures that J itself complies with the boundary conditions, right?
     // ApplyJfieldBoundary(lev, Jfield[0].get(), Jfield[1].get(), Jfield[2].get());
-    for (int i=0; i<3; i++) { current_fp_ampere[i]->FillBoundary(warpx.Geom(lev).periodicity()); }
+    for (int i=0; i<3; i++) { current_fp_plasma[i]->FillBoundary(warpx.Geom(lev).periodicity()); }
+
+    // Subtract external current from "Ampere" current calculated above. Note
+    // we need to include 1 ghost cell since later we will interpolate the
+    // plasma current to a nodal grid.
+    ablastr::fields::VectorField current_fp_external = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_external, lev);
+    for (int i=0; i<3; i++) {
+        current_fp_plasma[i]->minus(*current_fp_external[i], 0, 1, 1);
+    }
+
 }
 
 void HybridPICModel::HybridPICSolveE (
@@ -463,19 +473,15 @@ void HybridPICModel::HybridPICSolveE (
     const int lev, PatchType patch_type,
     const bool solve_for_Faraday) const
 {
-
     auto& warpx = WarpX::GetInstance();
 
-    ablastr::fields::VectorField current_fp_ampere = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_ampere, lev);
-    const ablastr::fields::VectorField current_fp_external = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_external, lev);
+    ablastr::fields::VectorField current_fp_plasma = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_plasma, lev);
     const ablastr::fields::ScalarField electron_pressure_fp = warpx.m_fields.get(FieldType::hybrid_electron_pressure_fp, lev);
 
     // Solve E field in regular cells
     warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE(
-        Efield, current_fp_ampere, Jfield, current_fp_external,
-        Bfield, rhofield,
-        *electron_pressure_fp,
-        edge_lengths, lev, this, solve_for_Faraday
+        Efield, current_fp_plasma, Jfield, Bfield, rhofield,
+        *electron_pressure_fp, edge_lengths, lev, this, solve_for_Faraday
     );
     warpx.ApplyEfieldBoundary(lev, patch_type);
 }
@@ -679,8 +685,8 @@ void HybridPICModel::FieldPush (
 {
     auto& warpx = WarpX::GetInstance();
 
-    // Calculate J = curl x B / mu0
-    CalculateCurrentAmpere(Bfield, edge_lengths);
+    // Calculate J = curl x B / mu0 - J_ext
+    CalculatePlasmaCurrent(Bfield, edge_lengths);
     // Calculate the E-field from Ohm's law
     HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
     warpx.FillBoundaryE(ng, nodal_sync);