Added new arrays with flags

Source/EmbeddedBoundary/WarpXInitEB.cpp

@@ -291,9 +291,60 @@ WarpX::ScaleAreas (ablastr::fields::VectorField& face_areas,
     }
 }
 
+void
+WarpX::MarkUpdateECells (
+    amrex::EBFArrayBoxFactory const & eb_fact,
+    int const lev )
+{
+
+    using ablastr::fields::Direction;
+    using warpx::fields::FieldType;
+
+    // Extract structures for embedded boundaries
+    amrex::FabArray<amrex::EBCellFlagFab> const& eb_flag = eb_fact.getMultiEBCellFlagFab();
+
+    for (int idim = 0; idim < 3; ++idim) {
+
+        // TODO: Handle guard cells
+
+        for (amrex::MFIter mfi(*m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)); mfi.isValid(); ++mfi) {
+
+            const amrex::Box& box = mfi.tilebox();
+            auto const & update_E_arr = m_eb_update_E[lev][idim]->array(mfi);
+
+            amrex::FabType const fab_type = eb_flag[mfi].getType(mfi.tilebox(amrex::IntVect::TheCellVector()));
+
+            if (fab_type == amrex::FabType::regular) {
+
+                // every cell in box is all regular
+                amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    update_E_arr(i, j, k) = 1;
+                });
+
+            } else if (fab_type == amrex::FabType::covered) {
+
+                // every cell in box is all covered
+                amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    update_E_arr(i, j, k) = 0;
+                });
+
+            } else {
+
+                auto const & flag = eb_flag[mfi].array();
+                amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    update_E_arr(i, j, k) = flag(i, j, k).isRegular();
+                });
+
+            }
+
+        }
+
+    }
+
+}
 
 void
-WarpX::MarkCells ()
+WarpX::MarkExtensionCells ()
 {
     using ablastr::fields::Direction;
     using warpx::fields::FieldType;
@@ -302,7 +353,7 @@ WarpX::MarkCells ()
     auto const &cell_size = CellSize(maxLevel());
 
 #if !defined(WARPX_DIM_3D) && !defined(WARPX_DIM_XZ)
-    WARPX_ABORT_WITH_MESSAGE("MarkCells only implemented in 2D and 3D");
+    WARPX_ABORT_WITH_MESSAGE("MarkExtensionCells only implemented in 2D and 3D");
 #endif
 
     for (int idim = 0; idim < 3; ++idim) {

Source/Initialization/WarpXInitData.cpp

@@ -1300,6 +1300,9 @@ void WarpX::InitializeEBGridData (int lev)
 
             auto const eb_fact = fieldEBFactory(lev);
 
+            MarkUpdateECells( eb_fact, lev );
+
+            // TODO: move inside if condition for ECT
             auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev);
             ComputeEdgeLengths(edge_lengths_lev, eb_fact);
             ScaleEdges(edge_lengths_lev, CellSize(lev));
@@ -1309,7 +1312,7 @@ void WarpX::InitializeEBGridData (int lev)
             ScaleAreas(face_areas_lev, CellSize(lev));
 
             if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) {
-                MarkCells();
+                MarkExtensionCells();
                 ComputeFaceExtensions();
             }
         }

Source/WarpX.H

@@ -1013,6 +1013,14 @@ public:
     void InitEB ();
 
 #ifdef AMREX_USE_EB
+    /** Set a flag to indicate which cells are in the valid domain and should be updated
+     * by the field solve, and which cells are in the EB and should not be updated.
+     * TODO: Add detail about criterion for Yee and ECT solver
+     */
+    void MarkUpdateECells (
+        amrex::EBFArrayBoxFactory const & eb_fact,
+        int const lev );
+
     /**
     * \brief Compute the length of the mesh edges. Here the length is a value in [0, 1].
     *        An edge of length 0 is fully covered.
@@ -1044,7 +1052,7 @@ public:
     *          - 2 for stable cells which have been intruded
     *        Here we cannot know if a cell is intruded or not so we initialize all stable cells with 1
     */
-    void MarkCells();
+    void MarkExtensionCells();
 #endif
 
     /**
@@ -1393,17 +1401,22 @@ private:
     mutable amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > Afield_dotMask;
     mutable amrex::Vector<            std::unique_ptr<amrex::iMultiFab>     > phi_dotMask;
 
+    /** EB: Flag to indicate whether the E location is inside the EB and therefore E should not be updated.
+     * (One array per level and per direction, due to staggering)
+    */
+    amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > > m_eb_update_E;
+
     /** EB: for every mesh face flag_info_face contains a:
      *          * 0 if the face needs to be extended
      *          * 1 if the face is large enough to lend area to other faces
      *          * 2 if the face is actually intruded by other face
-     * It is initialized in WarpX::MarkCells
+     * It is initialized in WarpX::MarkExtensionCells
      * This is only used for the ECT solver.*/
     amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>, 3 > > m_flag_info_face;
     /** EB: for every mesh face face flag_ext_face contains a:
      *          * 1 if the face needs to be extended
      *          * 0 otherwise
-     * It is initialized in WarpX::MarkCells and then modified in WarpX::ComputeOneWayExtensions
+     * It is initialized in WarpX::MarkExtensionCells and then modified in WarpX::ComputeOneWayExtensions
      * and in WarpX::ComputeEightWaysExtensions
      * This is only used for the ECT solver.*/
     amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>, 3 > > m_flag_ext_face;

Source/WarpX.cpp

@@ -323,6 +323,7 @@ WarpX::WarpX ()
     Afield_dotMask.resize(nlevs_max);
     phi_dotMask.resize(nlevs_max);
 
+    m_eb_update_E.resize(nlevs_max);
     m_flag_info_face.resize(nlevs_max);
     m_flag_ext_face.resize(nlevs_max);
     m_borrowing.resize(nlevs_max);
@@ -2301,6 +2302,15 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         // EB info are needed only at the finest level
         if (lev == maxLevel()) {
             if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) {
+
+               AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps,
+                                  guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]");
+                AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps,
+                                  guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]");
+                AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps,
+                                  guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]");
+
+                // TODO: do not allocate these arrays anymore with the Yee solver
                 //! EB: Lengths of the mesh edges
                 m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag),
                     dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt);