Add: MultiFabRegister

Source/BoundaryConditions/PML.cpp

@@ -570,6 +570,8 @@ PML::PML (const int lev, const BoxArray& grid_ba,
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!eb_enabled, "PML: eb_enabled is true but was not compiled in.");
 #endif
 
+    using ablastr::fields::Direction;
+
     // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain or fine patch(es)
     // (instead of extending `ncell` outside of the physical domain or fine patch(es))
     // In order to implement this, we define a new reduced Box Array ensuring that it does not
@@ -697,12 +699,15 @@ PML::PML (const int lev, const BoxArray& grid_ba,
     const int ncompe = (m_dive_cleaning) ? 3 : 2;
     const int ncompb = (m_divb_cleaning) ? 3 : 2;
 
-    const amrex::BoxArray ba_Ex = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::Efield_fp, 0,0).ixType().toIntVect());
-    const amrex::BoxArray ba_Ey = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::Efield_fp, 0,1).ixType().toIntVect());
-    const amrex::BoxArray ba_Ez = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::Efield_fp, 0,2).ixType().toIntVect());
-    WarpX::AllocInitMultiFab(pml_E_fp[0], ba_Ex, dm, ncompe, nge, lev, "pml_E_fp[x]", 0.0_rt);
-    WarpX::AllocInitMultiFab(pml_E_fp[1], ba_Ey, dm, ncompe, nge, lev, "pml_E_fp[y]", 0.0_rt);
-    WarpX::AllocInitMultiFab(pml_E_fp[2], ba_Ez, dm, ncompe, nge, lev, "pml_E_fp[z]", 0.0_rt);
+    auto& warpx = WarpX::GetInstance();
+    using ablastr::fields::Direction;
+
+    const amrex::BoxArray ba_Ex = amrex::convert(ba, warpx.m_fields.get("Efield_fp",Direction{0},0)->ixType().toIntVect());
+    const amrex::BoxArray ba_Ey = amrex::convert(ba, warpx.m_fields.get("Efield_fp",Direction{1},0)->ixType().toIntVect());
+    const amrex::BoxArray ba_Ez = amrex::convert(ba, warpx.m_fields.get("Efield_fp",Direction{2},0)->ixType().toIntVect());
+    warpx.m_fields.alloc_init("pml_E_fp", Direction{0}, lev, ba_Ex, dm, ncompe, nge, 0.0_rt);
+    warpx.m_fields.alloc_init("pml_E_fp", Direction{1}, lev, ba_Ey, dm, ncompe, nge, 0.0_rt);
+    warpx.m_fields.alloc_init("pml_E_fp", Direction{2}, lev, ba_Ez, dm, ncompe, nge, 0.0_rt);
 
     const amrex::BoxArray ba_Bx = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::Bfield_fp, 0,0).ixType().toIntVect());
     const amrex::BoxArray ba_By = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::Bfield_fp, 0,1).ixType().toIntVect());
@@ -711,9 +716,9 @@ PML::PML (const int lev, const BoxArray& grid_ba,
     WarpX::AllocInitMultiFab(pml_B_fp[1], ba_By, dm, ncompb, ngb, lev, "pml_B_fp[y]", 0.0_rt);
     WarpX::AllocInitMultiFab(pml_B_fp[2], ba_Bz, dm, ncompb, ngb, lev, "pml_B_fp[z]", 0.0_rt);
 
-    const amrex::BoxArray ba_jx = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::current_fp, 0,0).ixType().toIntVect());
-    const amrex::BoxArray ba_jy = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::current_fp, 0,1).ixType().toIntVect());
-    const amrex::BoxArray ba_jz = amrex::convert(ba, WarpX::GetInstance().getField(FieldType::current_fp, 0,2).ixType().toIntVect());
+    const amrex::BoxArray ba_jx = amrex::convert(ba, WarpX::GetInstance().m_fields.get("current_fp", Direction{0}, 0)->ixType().toIntVect());
+    const amrex::BoxArray ba_jy = amrex::convert(ba, WarpX::GetInstance().m_fields.get("current_fp", Direction{1}, 0)->ixType().toIntVect());
+    const amrex::BoxArray ba_jz = amrex::convert(ba, WarpX::GetInstance().m_fields.get("current_fp", Direction{2}, 0)->ixType().toIntVect());
     WarpX::AllocInitMultiFab(pml_j_fp[0], ba_jx, dm, 1, ngb, lev, "pml_j_fp[x]", 0.0_rt);
     WarpX::AllocInitMultiFab(pml_j_fp[1], ba_jy, dm, 1, ngb, lev, "pml_j_fp[y]", 0.0_rt);
     WarpX::AllocInitMultiFab(pml_j_fp[2], ba_jz, dm, 1, ngb, lev, "pml_j_fp[z]", 0.0_rt);
@@ -864,9 +869,9 @@ PML::PML (const int lev, const BoxArray& grid_ba,
             WarpX::AllocInitMultiFab( pml_G_cp, cba_G_nodal, cdm, 3, ngf, lev, "pml_G_cp", 0.0_rt);
         }
 
-        const amrex::BoxArray cba_jx = amrex::convert(cba, WarpX::GetInstance().getField(FieldType::current_cp, 1,0).ixType().toIntVect());
-        const amrex::BoxArray cba_jy = amrex::convert(cba, WarpX::GetInstance().getField(FieldType::current_cp, 1,1).ixType().toIntVect());
-        const amrex::BoxArray cba_jz = amrex::convert(cba, WarpX::GetInstance().getField(FieldType::current_cp, 1,2).ixType().toIntVect());
+        const amrex::BoxArray cba_jx = amrex::convert(cba, WarpX::GetInstance().m_fields.get("current_cp", Direction{0}, 0)->ixType().toIntVect());
+        const amrex::BoxArray cba_jy = amrex::convert(cba, WarpX::GetInstance().m_fields.get("current_cp", Direction{1}, 0)->ixType().toIntVect());
+        const amrex::BoxArray cba_jz = amrex::convert(cba, WarpX::GetInstance().m_fields.get("current_cp", Direction{2}, 0)->ixType().toIntVect());
         WarpX::AllocInitMultiFab(pml_j_cp[0], cba_jx, cdm, 1, ngb, lev, "pml_j_cp[x]", 0.0_rt);
         WarpX::AllocInitMultiFab(pml_j_cp[1], cba_jy, cdm, 1, ngb, lev, "pml_j_cp[y]", 0.0_rt);
         WarpX::AllocInitMultiFab(pml_j_cp[2], cba_jz, cdm, 1, ngb, lev, "pml_j_cp[z]", 0.0_rt);

Source/BoundaryConditions/PML_RZ.cpp

@@ -42,13 +42,15 @@ PML_RZ::PML_RZ (const int lev, const amrex::BoxArray& grid_ba, const amrex::Dist
       m_do_pml_in_domain(do_pml_in_domain),
       m_geom(geom)
 {
+    auto& warpx = WarpX::GetInstance();
+    using ablastr::fields::Direction;
 
-    const amrex::MultiFab & Er_fp = WarpX::GetInstance().getField(FieldType::Efield_fp, lev,0);
-    const amrex::MultiFab & Et_fp = WarpX::GetInstance().getField(FieldType::Efield_fp, lev,1);
+    const amrex::MultiFab & Er_fp = *warpx.m_fields.get("Efield_fp",Direction{0},lev);
+    const amrex::MultiFab & Et_fp = *warpx.m_fields.get("Efield_fp",Direction{1},lev);
     const amrex::BoxArray ba_Er = amrex::convert(grid_ba, Er_fp.ixType().toIntVect());
     const amrex::BoxArray ba_Et = amrex::convert(grid_ba, Et_fp.ixType().toIntVect());
-    WarpX::AllocInitMultiFab(pml_E_fp[0], ba_Er, grid_dm, Er_fp.nComp(), Er_fp.nGrowVect(), lev, "pml_E_fp[0]", 0.0_rt);
-    WarpX::AllocInitMultiFab(pml_E_fp[1], ba_Et, grid_dm, Et_fp.nComp(), Et_fp.nGrowVect(), lev, "pml_E_fp[1]", 0.0_rt);
+    warpx.m_fields.alloc_init("pml_E_fp", Direction{0}, lev, ba_Er, grid_dm, Er_fp.nComp(), Er_fp.nGrowVect(), 0.0_rt);
+    warpx.m_fields.alloc_init("pml_E_fp", Direction{1}, lev, ba_Et, grid_dm, Et_fp.nComp(), Et_fp.nGrowVect(), 0.0_rt);
 
     const amrex::MultiFab & Br_fp = WarpX::GetInstance().getField(FieldType::Bfield_fp, lev,0);
     const amrex::MultiFab & Bt_fp = WarpX::GetInstance().getField(FieldType::Bfield_fp, lev,1);

Source/BoundaryConditions/WarpXEvolvePML.cpp

@@ -63,9 +63,10 @@ WarpX::DampPML (const int lev, PatchType patch_type)
     WARPX_PROFILE("WarpX::DampPML()");
 #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT)
     if (pml_rz[lev]) {
-        pml_rz[lev]->ApplyDamping(Efield_fp[lev][1].get(), Efield_fp[lev][2].get(),
-                                  Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get(),
-                                  dt[lev]);
+        pml_rz[lev]->ApplyDamping( m_fields.get("Efield_fp",Direction{1},lev),
+                                   m_fields.get("Efield_fp",Direction{2},lev),
+                                   Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get(),
+                                   dt[lev]);
     }
 #endif
     if (pml[lev]) {
@@ -338,15 +339,21 @@ WarpX::DampJPML (int lev, PatchType patch_type)
 void
 WarpX::CopyJPML ()
 {
+    using ablastr::fields::Direction;
+
     for (int lev = 0; lev <= finest_level; ++lev)
     {
         if (pml[lev] && pml[lev]->ok()){
-            pml[lev]->CopyJtoPMLs({ current_fp[lev][0].get(),
-                                  current_fp[lev][1].get(),
-                                  current_fp[lev][2].get() },
-                                { current_cp[lev][0].get(),
-                                  current_cp[lev][1].get(),
-                                  current_cp[lev][2].get() });
+            pml[lev]->CopyJtoPMLs({
+                m_fields.get("current_fp", Direction{0}, lev),
+                m_fields.get("current_fp", Direction{1}, lev),
+                m_fields.get("current_fp", Direction{2}, lev)
+                }, {
+                m_fields.get("current_cp", Direction{0}, lev),
+                m_fields.get("current_cp", Direction{1}, lev),
+                m_fields.get("current_cp", Direction{2}, lev)
+                }
+            );
         }
     }
 }

Source/BoundaryConditions/WarpXFieldBoundaries.cpp

@@ -50,12 +50,14 @@ namespace
 
 void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
 {
+    using ablastr::fields::Direction;
+
     if (::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi)) {
         if (patch_type == PatchType::fine) {
             PEC::ApplyPECtoEfield(
-                    {getFieldPointer(FieldType::Efield_fp, lev, 0),
-                    getFieldPointer(FieldType::Efield_fp, lev, 1),
-                    getFieldPointer(FieldType::Efield_fp, lev, 2)},
+                    {m_fields.get("Efield_fp",Direction{0},lev),
+                     m_fields.get("Efield_fp",Direction{1},lev),
+                     m_fields.get("Efield_fp",Direction{2},lev)},
                     field_boundary_lo, field_boundary_hi,
                     get_ng_fieldgather(), Geom(lev),
                     lev, patch_type, ref_ratio);
@@ -92,9 +94,9 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type)
 
 #ifdef WARPX_DIM_RZ
     if (patch_type == PatchType::fine) {
-        ApplyFieldBoundaryOnAxis(getFieldPointer(FieldType::Efield_fp, lev, 0),
-                                 getFieldPointer(FieldType::Efield_fp, lev, 1),
-                                 getFieldPointer(FieldType::Efield_fp, lev, 2), lev);
+        ApplyFieldBoundaryOnAxis(m_fields.get("Efield_fp",Direction{0},lev),
+                                 m_fields.get("Efield_fp",Direction{1},lev),
+                                 m_fields.get("Efield_fp",Direction{2},lev), lev);
     } else {
         ApplyFieldBoundaryOnAxis(getFieldPointer(FieldType::Efield_cp, lev, 0),
                                  getFieldPointer(FieldType::Efield_cp, lev, 1),
@@ -131,8 +133,9 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d
     if (lev == 0) {
         if (a_dt_type == DtType::FirstHalf) {
             if(::isAnyBoundary<FieldBoundaryType::Absorbing_SilverMueller>(field_boundary_lo, field_boundary_hi)){
+                auto Efield_fp_new = m_fields.get_mr_levels_alldirs("Efield_fp",max_level); // JRA, new to prevent shadow
                 m_fdtd_solver_fp[0]->ApplySilverMuellerBoundary(
-                Efield_fp[lev], Bfield_fp[lev],
+                Efield_fp_new[lev], Bfield_fp[lev],
                 Geom(lev).Domain(), dt[lev],
                 field_boundary_lo, field_boundary_hi);
             }

Source/Diagnostics/BTDiagnostics.cpp

@@ -525,6 +525,8 @@ BTDiagnostics::DefineCellCenteredMultiFab(int lev)
 void
 BTDiagnostics::InitializeFieldFunctors (int lev)
 {
+    using ablastr::fields::Direction;
+
     // Initialize fields functors only if do_back_transformed_fields is selected
     if (!m_do_back_transformed_fields) { return; }
 
@@ -579,11 +581,11 @@ BTDiagnostics::InitializeFieldFunctors (int lev)
         } else if ( m_cellcenter_varnames[comp] == "Bz" ){
             m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Bfield_aux, lev, 2), lev, m_crse_ratio);
         } else if ( m_cellcenter_varnames[comp] == "jx" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 0), lev, m_crse_ratio);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp",Direction{0}, lev), lev, m_crse_ratio);
         } else if ( m_cellcenter_varnames[comp] == "jy" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 1), lev, m_crse_ratio);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp",Direction{1}, lev), lev, m_crse_ratio);
         } else if ( m_cellcenter_varnames[comp] == "jz" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 2), lev, m_crse_ratio);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp",Direction{2}, lev), lev, m_crse_ratio);
         } else if ( m_cellcenter_varnames[comp] == "rho" ){
             m_cell_center_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio);
         }
@@ -656,6 +658,8 @@ void
 BTDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev)
 {
 #ifdef WARPX_DIM_RZ
+    using ablastr::fields::Direction;
+
     auto & warpx = WarpX::GetInstance();
     const int ncomp_multimodefab = warpx.getFieldPointer(FieldType::Efield_aux, 0,0)->nComp();
     const int ncomp = ncomp_multimodefab;
@@ -695,11 +699,11 @@ BTDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev)
         } else if ( m_cellcenter_varnames_fields[comp] == "Bz" ){
             m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Bfield_aux, lev, 2), lev, m_crse_ratio, false, ncomp);
         } else if ( m_cellcenter_varnames_fields[comp] == "jr" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 0), lev, m_crse_ratio, false, ncomp);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp", Direction{0}, lev), lev, m_crse_ratio, false, ncomp);
         } else if ( m_cellcenter_varnames_fields[comp] == "jt" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 1), lev, m_crse_ratio, false, ncomp);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp", Direction{1}, lev), lev, m_crse_ratio, false, ncomp);
         } else if ( m_cellcenter_varnames_fields[comp] == "jz" ){
-            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::current_fp, lev, 2), lev, m_crse_ratio, false, ncomp);
+            m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.m_fields.get("current_fp", Direction{2}, lev), lev, m_crse_ratio, false, ncomp);
         } else if ( m_cellcenter_varnames_fields[comp] == "rho" ){
             m_cell_center_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio, false, -1, false, ncomp);
         }

Source/Diagnostics/ComputeDiagFunctors/JFunctor.cpp

@@ -29,9 +29,11 @@ JFunctor::JFunctor (const int dir, int lev,
 void
 JFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/) const
 {
+    using ablastr::fields::Direction;
+
     auto& warpx = WarpX::GetInstance();
     /** pointer to source multifab (can be multi-component) */
-    amrex::MultiFab* m_mf_src = warpx.getFieldPointer(FieldType::current_fp, m_lev, m_dir);
+    amrex::MultiFab* m_mf_src = warpx.m_fields.get("current_fp",Direction{m_dir},m_lev);
 
     // Deposit current if no solver or the electrostatic solver is being used
     if (m_deposit_current)
@@ -40,16 +42,16 @@ JFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
         current_fp_temp.resize(1);
 
-        const auto& current_fp_x = warpx.getField(FieldType::current_fp, m_lev,0);
+        const auto& current_fp_x = *warpx.m_fields.get("current_fp",Direction{0},m_lev);
         current_fp_temp[0][0] = std::make_unique<amrex::MultiFab>(
             current_fp_x, amrex::make_alias, 0, current_fp_x.nComp()
         );
 
-        const auto& current_fp_y = warpx.getField(FieldType::current_fp, m_lev,1);
+        const auto& current_fp_y = *warpx.m_fields.get("current_fp",Direction{1},m_lev);
         current_fp_temp[0][1] = std::make_unique<amrex::MultiFab>(
             current_fp_y, amrex::make_alias, 0, current_fp_y.nComp()
         );
-        const auto& current_fp_z = warpx.getField(FieldType::current_fp, m_lev,2);
+        const auto& current_fp_z = *warpx.m_fields.get("current_fp",Direction{2},m_lev);
         current_fp_temp[0][2] = std::make_unique<amrex::MultiFab>(
             current_fp_z, amrex::make_alias, 0, current_fp_z.nComp()
         );

Source/Diagnostics/ComputeDiagFunctors/JdispFunctor.cpp

@@ -27,11 +27,13 @@ JdispFunctor::JdispFunctor (int dir, int lev,
 void
 JdispFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/) const
 {
+    using ablastr::fields::Direction;
+
     auto& warpx = WarpX::GetInstance();
     auto* hybrid_pic_model = warpx.get_pointer_HybridPICModel();
 
     /** pointer to total simulation current (J) multifab */
-    amrex::MultiFab* mf_j = warpx.getFieldPointer(FieldType::current_fp, m_lev, m_dir);
+    amrex::MultiFab* mf_j = warpx.m_fields.get("current_fp",Direction{m_dir},m_lev);
 
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(hybrid_pic_model,
         "Displacement current diagnostic is only implemented for the HybridPICModel.");