General parser function for external fields (#5349)

This feature replaces the InitializeExternalFieldsOnGridUsingParser with
a more general function ComputeExternalFieldOnGridUsingParser. The new
function takes parsed functions with four dimensions (x,y,z,t), so that
it can be used to evaluate functions throughout the simulation, for
example time-dependent external fields.

The function ComputeExternalFieldOnGridUsingParser has optional edge
length and face areas.

---------

Co-authored-by: Kristoffer Lindvall <[email protected]>

Author: kli-jfp

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H

@@ -53,13 +53,7 @@ public:
      * external current multifab. Note the external current can be a function
      * of time and therefore this should be re-evaluated at every step.
      */
-    void GetCurrentExternal (
-        ablastr::fields::MultiLevelVectorField const& edge_lengths
-    );
-    void GetCurrentExternal (
-        ablastr::fields::VectorField const& edge_lengths,
-        int lev
-    );
+    void GetCurrentExternal ();
 
     /**
      * \brief

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp

@@ -221,167 +221,32 @@ void HybridPICModel::InitData ()
     // if the current is time dependent which is what needs to be done to
     // write time independent fields on the first step.
     for (int lev = 0; lev <= warpx.finestLevel(); ++lev) {
-        auto edge_lengths = std::array<std::unique_ptr<amrex::MultiFab>, 3>();
-#ifdef AMREX_USE_EB
-        if (EB::enabled()) {
-            using ablastr::fields::Direction;
-            auto const & edge_lengths_x = *warpx.m_fields.get(FieldType::edge_lengths, Direction{0}, lev);
-            auto const & edge_lengths_y = *warpx.m_fields.get(FieldType::edge_lengths, Direction{1}, lev);
-            auto const & edge_lengths_z = *warpx.m_fields.get(FieldType::edge_lengths, Direction{2}, lev);
-
-            edge_lengths = std::array< std::unique_ptr<amrex::MultiFab>, 3 >{
-                std::make_unique<amrex::MultiFab>(
-                    edge_lengths_x, amrex::make_alias, 0, edge_lengths_x.nComp()),
-                std::make_unique<amrex::MultiFab>(
-                    edge_lengths_y, amrex::make_alias, 0, edge_lengths_y.nComp()),
-                std::make_unique<amrex::MultiFab>(
-                    edge_lengths_z, amrex::make_alias, 0, edge_lengths_z.nComp())
-            };
-        }
-#endif
-        GetCurrentExternal(ablastr::fields::a2m(edge_lengths), lev);
+        warpx.ComputeExternalFieldOnGridUsingParser(
+            FieldType::hybrid_current_fp_external,
+            m_J_external[0],
+            m_J_external[1],
+            m_J_external[2],
+            lev, PatchType::fine, 'e',
+            warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev),
+            warpx.m_fields.get_alldirs(FieldType::face_areas, lev));
     }
 }
 
-void HybridPICModel::GetCurrentExternal (
-    ablastr::fields::MultiLevelVectorField const& edge_lengths)
+void HybridPICModel::GetCurrentExternal ()
 {
     if (!m_external_field_has_time_dependence) { return; }
 
     auto& warpx = WarpX::GetInstance();
     for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
     {
-        GetCurrentExternal(edge_lengths[lev], lev);
-    }
-}
-
-
-void HybridPICModel::GetCurrentExternal (
-    ablastr::fields::VectorField const& edge_lengths,
-    int lev)
-{
-    // This logic matches closely to WarpX::InitializeExternalFieldsOnGridUsingParser
-    // except that the parsers include time dependence.
-    auto & warpx = WarpX::GetInstance();
-
-    auto t = warpx.gett_new(lev);
-
-    auto dx_lev = warpx.Geom(lev).CellSizeArray();
-    const RealBox& real_box = warpx.Geom(lev).ProbDomain();
-
-    using ablastr::fields::Direction;
-    amrex::MultiFab * mfx = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{0}, lev);
-    amrex::MultiFab * mfy = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{1}, lev);
-    amrex::MultiFab * mfz = warpx.m_fields.get(FieldType::hybrid_current_fp_external, Direction{2}, lev);
-
-    const amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect();
-    const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect();
-    const amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect();
-
-    // avoid implicit lambda capture
-    auto Jx_external = m_J_external[0];
-    auto Jy_external = m_J_external[1];
-    auto Jz_external = m_J_external[2];
-
-    for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-    {
-        const amrex::Box& tbx = mfi.tilebox( x_nodal_flag, mfx->nGrowVect() );
-        const amrex::Box& tby = mfi.tilebox( y_nodal_flag, mfy->nGrowVect() );
-        const amrex::Box& tbz = mfi.tilebox( z_nodal_flag, mfz->nGrowVect() );
-
-        auto const& mfxfab = mfx->array(mfi);
-        auto const& mfyfab = mfy->array(mfi);
-        auto const& mfzfab = mfz->array(mfi);
-
-        amrex::Array4<amrex::Real> lx, ly, lz;
-        if (EB::enabled()) {
-            lx = edge_lengths[0]->array(mfi);
-            ly = edge_lengths[1]->array(mfi);
-            lz = edge_lengths[2]->array(mfi);
-        }
-
-        amrex::ParallelFor (tbx, tby, tbz,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                // skip if node is covered by an embedded boundary
-                if (lx && lx(i, j, k) <= 0) { return; }
-
-                // Shift required in the x-, y-, or z- position
-                // depending on the index type of the multifab
-#if defined(WARPX_DIM_1D_Z)
-                const amrex::Real x = 0._rt;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                const amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
-#else
-                const amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real fac_y = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
-                const amrex::Real fac_z = (1._rt - x_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
-                const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
-#endif
-                // Initialize the x-component of the field.
-                mfxfab(i,j,k) = Jx_external(x,y,z,t);
-            },
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                // skip if node is covered by an embedded boundary
-                if (ly && ly(i, j, k) <= 0) { return; }
-
-#if defined(WARPX_DIM_1D_Z)
-                const amrex::Real x = 0._rt;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                const amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
-#elif defined(WARPX_DIM_3D)
-                const amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real fac_y = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
-                const amrex::Real fac_z = (1._rt - y_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
-                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);
-            },
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                // skip if node is covered by an embedded boundary
-                if (lz && lz(i, j, k) <= 0) { return; }
-
-#if defined(WARPX_DIM_1D_Z)
-                const amrex::Real x = 0._rt;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[0] + real_box.lo(0) + fac_z;
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                const amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real y = 0._rt;
-                const amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z;
-#elif defined(WARPX_DIM_3D)
-                const amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt;
-                const amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x;
-                const amrex::Real fac_y = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt;
-                const amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y;
-                const amrex::Real fac_z = (1._rt - z_nodal_flag[2]) * dx_lev[2] * 0.5_rt;
-                const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z;
-#endif
-                // Initialize the z-component of the field.
-                mfzfab(i,j,k) = Jz_external(x,y,z,t);
-            }
-        );
+        warpx.ComputeExternalFieldOnGridUsingParser(
+            FieldType::hybrid_current_fp_external,
+            m_J_external[0],
+            m_J_external[1],
+            m_J_external[2],
+            lev, PatchType::fine, 'e',
+            warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev),
+            warpx.m_fields.get_alldirs(FieldType::face_areas, lev));
     }
 }
 

Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp

@@ -68,8 +68,7 @@ void WarpX::HybridPICEvolveFields ()
     const int sub_steps = m_hybrid_pic_model->m_substeps;
 
     // Get the external current
-    m_hybrid_pic_model->GetCurrentExternal(
-        m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level));
+    m_hybrid_pic_model->GetCurrentExternal();
 
     // Reference hybrid-PIC multifabs
     ablastr::fields::MultiLevelScalarField rho_fp_temp = m_fields.get_mr_levels(FieldType::hybrid_rho_fp_temp, finest_level);

Source/Fluids/WarpXFluidContainer.cpp

@@ -1010,24 +1010,23 @@ void WarpXFluidContainer::GatherAndPush (
     // External field parsers
     external_e_fields = (m_E_ext_s == "parse_e_ext_function");
     external_b_fields = (m_B_ext_s == "parse_b_ext_function");
+
     amrex::ParserExecutor<4> Exfield_parser;
     amrex::ParserExecutor<4> Eyfield_parser;
     amrex::ParserExecutor<4> Ezfield_parser;
     amrex::ParserExecutor<4> Bxfield_parser;
     amrex::ParserExecutor<4> Byfield_parser;
     amrex::ParserExecutor<4> Bzfield_parser;
+
     if (external_e_fields){
-        constexpr int num_arguments = 4; //x,y,z,t
-        Exfield_parser = m_Ex_parser->compile<num_arguments>();
-        Eyfield_parser = m_Ey_parser->compile<num_arguments>();
-        Ezfield_parser = m_Ez_parser->compile<num_arguments>();
+        Exfield_parser = m_Ex_parser->compile<4>();
+        Eyfield_parser = m_Ey_parser->compile<4>();
+        Ezfield_parser = m_Ez_parser->compile<4>();
     }
-
     if (external_b_fields){
-        constexpr int num_arguments = 4; //x,y,z,t
-        Bxfield_parser = m_Bx_parser->compile<num_arguments>();
-        Byfield_parser = m_By_parser->compile<num_arguments>();
-        Bzfield_parser = m_Bz_parser->compile<num_arguments>();
+        Bxfield_parser = m_Bx_parser->compile<4>();
+        Byfield_parser = m_By_parser->compile<4>();
+        Bzfield_parser = m_Bz_parser->compile<4>();
     }
 
 

Source/Initialization/ExternalField.cpp

@@ -127,11 +127,11 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx)
             str_Bz_ext_grid_function);
 
         Bxfield_parser = std::make_unique<amrex::Parser>(
-            utils::parser::makeParser(str_Bx_ext_grid_function,{"x","y","z"}));
+            utils::parser::makeParser(str_Bx_ext_grid_function,{"x","y","z","t"}));
         Byfield_parser = std::make_unique<amrex::Parser>(
-            utils::parser::makeParser(str_By_ext_grid_function,{"x","y","z"}));
+            utils::parser::makeParser(str_By_ext_grid_function,{"x","y","z","t"}));
         Bzfield_parser = std::make_unique<amrex::Parser>(
-            utils::parser::makeParser(str_Bz_ext_grid_function,{"x","y","z"}));
+            utils::parser::makeParser(str_Bz_ext_grid_function,{"x","y","z","t"}));
     }
     //___________________________________________________________________________
 
@@ -163,11 +163,11 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx)
            str_Ez_ext_grid_function);
 
         Exfield_parser = std::make_unique<amrex::Parser>(
-           utils::parser::makeParser(str_Ex_ext_grid_function,{"x","y","z"}));
+           utils::parser::makeParser(str_Ex_ext_grid_function,{"x","y","z","t"}));
         Eyfield_parser = std::make_unique<amrex::Parser>(
-           utils::parser::makeParser(str_Ey_ext_grid_function,{"x","y","z"}));
+           utils::parser::makeParser(str_Ey_ext_grid_function,{"x","y","z","t"}));
         Ezfield_parser = std::make_unique<amrex::Parser>(
-           utils::parser::makeParser(str_Ez_ext_grid_function,{"x","y","z"}));
+           utils::parser::makeParser(str_Ez_ext_grid_function,{"x","y","z","t"}));
     }
     //___________________________________________________________________________