calculate J_plasma without taking account of EBs

roelof-groenewald · Sep 20, 2024 · 89de956 · 89de956
1 parent 7ce9538
commit 89de956
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 19 deletions.
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
@@ -103,11 +103,11 @@ void HybridPICModel::AllocateLevelMFs (int lev, const BoxArray& ba, const Distri
     // the external current density multifab matches the current staggering but no
     // ghost cells are needed since this multifab is only used to update local values
     WarpX::AllocInitMultiFab(current_fp_external[lev][0], amrex::convert(ba, jx_nodal_flag),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[x]", 0.0_rt);
+        dm, ncomps, IntVect::TheUnitVector(), lev, "current_fp_external[x]", 0.0_rt);
     WarpX::AllocInitMultiFab(current_fp_external[lev][1], amrex::convert(ba, jy_nodal_flag),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[y]", 0.0_rt);
+        dm, ncomps, IntVect::TheUnitVector(), lev, "current_fp_external[y]", 0.0_rt);
     WarpX::AllocInitMultiFab(current_fp_external[lev][2], amrex::convert(ba, jz_nodal_flag),
-        dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[z]", 0.0_rt);
+        dm, ncomps, IntVect::TheUnitVector(), lev, "current_fp_external[z]", 0.0_rt);
 
 #ifdef WARPX_DIM_RZ
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
@@ -418,13 +418,11 @@ void HybridPICModel::CalculatePlasmaCurrent (
     // ApplyJfieldBoundary(lev, Jfield[0].get(), Jfield[1].get(), Jfield[2].get());
     for (int i=0; i<3; i++) { current_fp_plasma[lev][i]->FillBoundary(warpx.Geom(lev).periodicity()); }
 
-    // Subtract external current from "Ampere" current calculated above
+    // 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.
     for (int i=0; i<3; i++) {
-        MultiFab::LinComb(
-            *current_fp_plasma[lev][i],
-            1._rt, *current_fp_plasma[lev][i], 0,
-            -1._rt, *current_fp_external[lev][i], 0, 0, 1, 0
-        );
+        current_fp_plasma[lev][i]->minus(*current_fp_external[lev][i], 0, 1, 1);
     }
 }
 

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
@@ -297,13 +297,17 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian (
         Real const one_over_mu0 = 1._rt / PhysConst::mu0;
 
         // Calculate the total current, using Ampere's law, on the same grid
-        // as the E-field
+        // as the E-field.
+        // Note that embedded boundary restrictions are not applied here since
+        // it is assumed that the magnetic field was initialized to respect
+        // the embedded boundary and since the B-field is not advanced inside
+        // the EB that boundary should remain accurate.
         amrex::ParallelFor(tjx, tjy, tjz,
 
             // Jx calculation
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 // Skip if this cell is fully covered by embedded boundaries
-                if (lx && lx(i, j, k) <= 0) { return; }
+                // if (lx && lx(i, j, k) <= 0) { return; }
 
                 Jx(i, j, k) = one_over_mu0 * (
                     - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k)
@@ -314,13 +318,13 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian (
             // Jy calculation
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 // Skip if this cell is fully covered by embedded boundaries
-#ifdef WARPX_DIM_3D
-                if (ly && 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 && (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0)) { return; }
-#endif
+// #ifdef WARPX_DIM_3D
+//                 if (ly && 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 && (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0)) { return; }
+// #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)
@@ -330,7 +334,7 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian (
             // Jz calculation
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 // Skip if this cell is fully covered by embedded boundaries
-                if (lz && lz(i,j,k) <= 0) { return; }
+                // if (lz && lz(i,j,k) <= 0) { return; }
 
                 Jz(i, j, k) = one_over_mu0 * (
                     - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k)