Adding grid based split field advance for Hybrid solver to improve pe…

…rformance.

Python/pywarpx/fields.py

@@ -696,6 +696,40 @@ def BzFPExternalWrapper(level=0, include_ghosts=False):
         mf_name="Bfield_fp_external[z]", level=level, include_ghosts=include_ghosts
     )
 
+def ExHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Efield_hyb_external[x]", level=level, include_ghosts=include_ghosts
+    )
+
+
+def EyHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Efield_hyb_external[y]", level=level, include_ghosts=include_ghosts
+    )
+
+
+def EzHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Efield_hyb_external[z]", level=level, include_ghosts=include_ghosts
+    )
+
+
+def BxHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Bfield_hyb_external[x]", level=level, include_ghosts=include_ghosts
+    )
+
+
+def ByHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Bfield_hyb_external[y]", level=level, include_ghosts=include_ghosts
+    )
+
+
+def BzHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="Bfield_hyb_external[z]", level=level, include_ghosts=include_ghosts
+    )
 
 def JxFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(

Python/pywarpx/picmi.py

@@ -1783,6 +1783,9 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
 
     Jx/y/z_external_function: str
         Function of space and time specifying external (non-plasma) currents.
+
+    Ex/y/z_external_function: str
+        Function of space and time specifying external (non-plasma) E-fields.
     """
 
     def __init__(
@@ -1798,6 +1801,12 @@ def __init__(
         Jx_external_function=None,
         Jy_external_function=None,
         Jz_external_function=None,
+        Ex_expression=None,
+        Ey_expression=None,
+        Ez_expression=None,
+        Bx_expression=None,
+        By_expression=None,
+        Bz_expression=None,
         **kw,
     ):
         self.grid = grid
@@ -1816,6 +1825,24 @@ def __init__(
         self.Jy_external_function = Jy_external_function
         self.Jz_external_function = Jz_external_function
 
+        self.add_external_fields = None
+
+        self.Ex_external_function = Ex_expression
+        self.Ey_external_function = Ey_expression
+        self.Ez_external_function = Ez_expression
+
+        self.Bx_external_function = Bx_expression
+        self.By_external_function = By_expression
+        self.Bz_external_function = Bz_expression
+
+        if (Ex_expression is not None
+            or Ey_expression is not None
+            or Ez_expression is not None
+            or Bx_expression is not None
+            or By_expression is not None
+            or Bz_expression is not None):
+            self.add_external_fields = True
+
         # Handle keyword arguments used in expressions
         self.user_defined_kw = {}
         for k in list(kw.keys()):
@@ -1864,6 +1891,43 @@ def solver_initialize_inputs(self):
                 self.Jz_external_function, self.mangle_dict
             ),
         )
+        pywarpx.hybridpicmodel.add_external_fields = self.add_external_fields
+        pywarpx.hybridpicmodel.__setattr__(
+            "Bx_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.Bx_external_function, self.mangle_dict
+            ),
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            "By_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.By_external_function, self.mangle_dict
+            ),
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            "Bz_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.Bz_external_function, self.mangle_dict
+            ),
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            "Ex_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.Ex_external_function, self.mangle_dict
+            ),
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            "Ey_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.Ey_external_function, self.mangle_dict
+            ),
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            "Ez_external_grid_function(x,y,z,t)",
+            pywarpx.my_constants.mangle_expression(
+                self.Ez_external_function, self.mangle_dict
+            ),
+        )
 
 
 class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H

@@ -21,6 +21,9 @@
 
 #include <AMReX_Array.H>
 #include <AMReX_REAL.H>
+#include <AMReX_BoxArray.H>
+#include <AMReX_IntVect.H>
+#include <AMReX_DistributionMapping.H>
 
 #include <optional>
 
@@ -38,10 +41,26 @@ public:
 
     /** Allocate hybrid-PIC specific multifabs. Called in constructor. */
     void AllocateMFs (int nlevs_max);
-    void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
-                           int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho,
-                           const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag,
-                           const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag);
+    void AllocateLevelMFs (
+        int lev,
+        const amrex::BoxArray& ba,
+        const amrex::DistributionMapping& dm,
+        const int ncomps,
+        const amrex::IntVect& ngJ,
+        const amrex::IntVect& ngRho,
+        const amrex::IntVect& ngE, 
+        const amrex::IntVect& ngB,
+        const amrex::IntVect& jx_nodal_flag,
+        const amrex::IntVect& jy_nodal_flag,
+        const amrex::IntVect& jz_nodal_flag,
+        const amrex::IntVect& rho_nodal_flag,
+        const amrex::IntVect& Ex_nodal_flag,
+        const amrex::IntVect& Ey_nodal_flag,
+        const amrex::IntVect& Ez_nodal_flag,
+        const amrex::IntVect& Bx_nodal_flag,
+        const amrex::IntVect& By_nodal_flag,
+        const amrex::IntVect& Bz_nodal_flag
+    );
 
     /** Helper function to clear values from hybrid-PIC specific multifabs. */
     void ClearLevel (int lev);
@@ -57,11 +76,39 @@ public:
     void GetCurrentExternal (
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths
     );
+
     void GetCurrentExternal (
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
         int lev
     );
 
+    void GetFieldsExternal (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths
+    );
+
+    void GetFieldsExternal (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths, 
+        amrex::Real t);
+
+    void GetFieldsExternal (
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        int lev
+    );
+
+    void GetExternalFieldFromExpression (
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
+        std::array< amrex::ParserExecutor<4>, 3> const& expression,
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        int lev
+    );
+
+    void GetExternalFieldFromExpression (
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& field,
+        std::array< amrex::ParserExecutor<4>, 3> const& expression,
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        int lev, amrex::Real t
+    );
+
     /**
      * \brief
      * Function to calculate the total current based on Ampere's law while
@@ -152,7 +199,7 @@ public:
      * charge density (and assumption of quasi-neutrality) using the user
      * specified electron equation of state.
      *
-     * \param[out] Pe_filed scalar electron pressure MultiFab at a given level
+     * \param[out] Pe_field scalar electron pressure MultiFab at a given level
      * \param[in] rho_field scalar ion charge density Multifab at a given level
      */
     void FillElectronPressureMF (
@@ -191,17 +238,17 @@ public:
     bool m_external_field_has_time_dependence = false;
 
     /** External E/B fields */
-    bool m_add_ext_particle_B_field = false;
-    std::string m_Bx_ext_part_function = "0.0";
-    std::string m_By_ext_part_function = "0.0";
-    std::string m_Bz_ext_part_function = "0.0";
+    bool m_add_external_fields = false;
+
+    std::string m_Bx_ext_grid_function = "0.0";
+    std::string m_By_ext_grid_function = "0.0";
+    std::string m_Bz_ext_grid_function = "0.0";
     std::array< std::unique_ptr<amrex::Parser>, 3> m_B_external_parser;
     std::array< amrex::ParserExecutor<4>, 3> m_B_external;
 
-    bool m_add_ext_particle_E_field = false;
-    std::string m_Ex_ext_part_function = "0.0";
-    std::string m_Ey_ext_part_function = "0.0";
-    std::string m_Ez_ext_part_function = "0.0";
+    std::string m_Ex_ext_grid_function = "0.0";
+    std::string m_Ey_ext_grid_function = "0.0";
+    std::string m_Ez_ext_grid_function = "0.0";
     std::array< std::unique_ptr<amrex::Parser>, 3> m_E_external_parser;
     std::array< amrex::ParserExecutor<4>, 3> m_E_external;
 
@@ -212,6 +259,11 @@ public:
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_external;
     amrex::Vector<            std::unique_ptr<amrex::MultiFab>      > electron_pressure_fp;
 
+    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_hyb_external;
+    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_hyb_external;
+    // amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_hyb_self;
+    // amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_hyb_self;
+
     // Helper functions to retrieve hybrid-PIC multifabs
     [[nodiscard]] amrex::MultiFab*
     get_pointer_current_fp_ampere  (int lev, int direction) const