Add external current handling to Ohm's law solver (#4405)

* add `WarpX::getedgelengths` and WarpX::getfaceareas` functions to return pointers to those multifabs

* add external current support to the hybrid-PIC solver

* add external current specification (for hybrid-PIC scheme) to picmi

* add RZ support for `FiniteDifferenceSolver::CalculateCurrentAmpere`

* allow an initial Bz field to be set in RZ

* code cleanup and addition of CI test

* avoid lambda capture issue

* update documentation to show external current implementation

* revert unwanted changes

* restore RZ support that was lost during rebase

* fix segfault when EB support is OFF

* fix codeQL issue

* add some details to docs

* the external current only needs to be calculated once per field solve step

* add description of `hybrid_pic_model.J[x/y/z]_external_grid_function(x, y, z, t)` to the documentation

Docs/source/theory/kinetic_fluid_hybrid_model.rst

@@ -35,25 +35,51 @@ neglecting the displacement current term :cite:p:`c-NIELSON1976`, giving,
 
         \mu_0\vec{J} = \vec{\nabla}\times\vec{B},
 
-where :math:`\vec{J} = \vec{J}_i - \vec{J}_e` is the total electrical current, i.e.
-the sum of electron and ion currents. Since ions are treated in the regular
-PIC manner, the ion current, :math:`\vec{J}_i`, is known during a simulation. Therefore,
+where :math:`\vec{J} = \sum_{s\neq e}\vec{J}_s + \vec{J}_e + \vec{J}_{ext}` is the total electrical current,
+i.e. the sum of electron and ion currents as well as any external current (not captured through plasma
+particles). Since ions are treated in the regular
+PIC manner, the ion current, :math:`\sum_{s\neq e}\vec{J}_s`, is known during a simulation. Therefore,
 given the magnetic field, the electron current can be calculated.
 
-If we now further assume electrons are inertialess, the electron momentum
-equation yields,
+The electron momentum transport equation (obtained from multiplying the Vlasov equation by mass and
+integrating over velocity), also called the generalized Ohm's law, is given by:
 
     .. math::
 
-        \frac{d(n_em_e\vec{V}_e)}{dt} = 0 = -en_e\vec{E}-\vec{J}_e\times\vec{B}-\nabla\cdot\vec{P}_e+en_e\vec{\eta}\cdot\vec{J},
+        en_e\vec{E} = \frac{m}{e}\frac{\partial \vec{J}_e}{\partial t} + \frac{m}{e^2}\left( \vec{U}_e\cdot\nabla \right) \vec{J}_e - \nabla\cdot {\overleftrightarrow P}_e - \vec{J}_e\times\vec{B}+\vec{R}_e
 
-where :math:`\vec{V_e}=\vec{J}_e/(en_e)`, :math:`\vec{P}_e` is the electron pressure
-tensor and :math:`\vec{\eta}` is the resistivity tensor. An expression for the electric field
-(generalized Ohm's law) can be obtained from the above as:
+where :math:`\vec{U}_e = \vec{J}_e/(en_e)` is the electron fluid velocity,
+:math:`{\overleftrightarrow P}_e` is the electron pressure tensor and
+:math:`\vec{R}_e` is the drag force due to collisions between electrons and ions.
+Applying the above momentum equation to the Maxwell-Faraday equation (:math:`\frac{\partial\vec{B}}{\partial t} = -\nabla\times\vec{E}`)
+and substituting in :math:`\vec{J}` calculated from the Maxwell-Ampere equation, gives,
 
     .. math::
 
-        \vec{E} = -\frac{1}{en_e}\left( \vec{J}_e\times\vec{B} + \nabla\cdot\vec{P}_e \right)+\vec{\eta}\cdot\vec{J}.
+        \frac{\partial\vec{J}_e}{\partial t} = -\frac{1}{\mu_0}\nabla\times\left(\nabla\times\vec{E}\right) - \frac{\partial\vec{J}_{ext}}{\partial t} - \sum_{s\neq e}\frac{\partial\vec{J}_s}{\partial t}.
+
+Plugging this back into the generalized Ohm' law gives:
+
+    .. math::
+
+        \left(en_e +\frac{m}{e\mu_0}\nabla\times\nabla\times\right)\vec{E} =&
+        - \frac{m}{e}\left( \frac{\partial\vec{J}_{ext}}{\partial t} + \sum_{s\neq e}\frac{\partial\vec{J}_s}{\partial t} \right) \\
+        &+ \frac{m}{e^2}\left( \vec{U}_e\cdot\nabla \right) \vec{J}_e - \nabla\cdot {\overleftrightarrow P}_e - \vec{J}_e\times\vec{B}+\vec{R}_e.
+
+If we now further assume electrons are inertialess (i.e. :math:`m=0`), the above equation simplifies to,
+
+    .. math::
+
+        en_e\vec{E} = -\vec{J}_e\times\vec{B}-\nabla\cdot{\overleftrightarrow P}_e+\vec{R}_e.
+
+Making the further simplifying assumptions that the electron pressure is isotropic and that
+the electron drag term can be written as a simple resistance
+i.e. :math:`\vec{R}_e = en_e\vec{\eta}\cdot\vec{J}`, brings us to the implemented form of
+Ohm's law:
+
+    .. math::
+
+        \vec{E} = -\frac{1}{en_e}\left( \vec{J}_e\times\vec{B} + \nabla P_e \right)+\vec{\eta}\cdot\vec{J}.
 
 Lastly, if an electron temperature is given from which the electron pressure can
 be calculated, the model is fully constrained and can be evolved given initial

Docs/source/usage/parameters.rst

@@ -2126,25 +2126,28 @@ Maxwell solver: kinetic-fluid hybrid
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 * ``hybrid_pic_model.elec_temp`` (`float`)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the electron temperature, in eV, used to calculate
-     the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the electron temperature, in eV, used to calculate
+    the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
 
 * ``hybrid_pic_model.n0_ref`` (`float`)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the reference density, in :math:`m^{-3}`, used to calculate
-     the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the reference density, in :math:`m^{-3}`, used to calculate
+    the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
 
 * ``hybrid_pic_model.gamma`` (`float`) optional (default ``5/3``)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the exponent used to calculate
-     the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the exponent used to calculate
+    the electron pressure (see :ref:`here <theory-hybrid-model-elec-temp>`).
 
 * ``hybrid_pic_model.plasma_resistivity(rho)`` (`float` or `str`) optional (default ``0``)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the plasma resistivity in :math:`\Omega m`.
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the plasma resistivity in :math:`\Omega m`.
+
+* ``hybrid_pic_model.J[x/y/z]_external_grid_function(x, y, z, t)`` (`float` or `str`) optional (default ``0``)
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the external current (on the grid) in :math:`A/m^2`.
 
 * ``hybrid_pic_model.n_floor`` (`float`) optional (default ``1``)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the plasma density floor, in :math:`m^{-3}`, which is useful since the generalized Ohm's law used to calculate the E-field includes a :math:`1/n` term.
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the plasma density floor, in :math:`m^{-3}`, which is useful since the generalized Ohm's law used to calculate the E-field includes a :math:`1/n` term.
 
 * ``hybrid_pic_model.substeps`` (`int`) optional (default ``100``)
-     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the number of sub-steps to take during the B-field update.
+    If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the number of sub-steps to take during the B-field update.
 
 .. note::
 

Python/pywarpx/picmi.py

@@ -1136,9 +1136,14 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
 
     substeps: int, default=100
         Number of substeps to take when updating the B-field.
+
+    Jx/y/z_external_function: str
+        Function of space and time specifying external (non-plasma) currents.
     """
     def __init__(self, grid, Te=None, n0=None, gamma=None,
-                 n_floor=None, plasma_resistivity=None, substeps=None, **kw):
+                 n_floor=None, plasma_resistivity=None, substeps=None,
+                 Jx_external_function=None, Jy_external_function=None,
+                 Jz_external_function=None, **kw):
         self.grid = grid
         self.method = "hybrid"
 
@@ -1150,6 +1155,10 @@ def __init__(self, grid, Te=None, n0=None, gamma=None,
 
         self.substeps = substeps
 
+        self.Jx_external_function = Jx_external_function
+        self.Jy_external_function = Jy_external_function
+        self.Jz_external_function = Jz_external_function
+
         self.handle_init(kw)
 
     def initialize_inputs(self):
@@ -1166,6 +1175,15 @@ def initialize_inputs(self):
             'plasma_resistivity(rho)', self.plasma_resistivity
         )
         pywarpx.hybridpicmodel.substeps = self.substeps
+        pywarpx.hybridpicmodel.__setattr__(
+            'Jx_external_grid_function(x,y,z,t)', self.Jx_external_function
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            'Jy_external_grid_function(x,y,z,t)', self.Jy_external_function
+        )
+        pywarpx.hybridpicmodel.__setattr__(
+            'Jz_external_grid_function(x,y,z,t)', self.Jz_external_function
+        )
 
 
 class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H

@@ -139,6 +139,7 @@ class FiniteDifferenceSolver
           * \param[out] Efield  vector of electric field MultiFabs updated at a given level
           * \param[in] Jfield   vector of total current MultiFabs at a given level
           * \param[in] Jifield  vector of ion current density MultiFabs at a given level
+          * \param[in] Jextfield  vector of external current density MultiFabs at a given level
           * \param[in] Bfield   vector of magnetic field MultiFabs at a given level
           * \param[in] rhofield scalar ion charge density Multifab at a given level
           * \param[in] Pefield  scalar electron pressure MultiFab at a given level
@@ -150,6 +151,7 @@ class FiniteDifferenceSolver
         void HybridPICSolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3>& Jfield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jifield,
+                      std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jextfield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
                       std::unique_ptr<amrex::MultiFab> const& rhofield,
                       std::unique_ptr<amrex::MultiFab> const& Pefield,
@@ -233,6 +235,7 @@ class FiniteDifferenceSolver
             std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jifield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jextfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
             std::unique_ptr<amrex::MultiFab> const& rhofield,
             std::unique_ptr<amrex::MultiFab> const& Pefield,
@@ -337,6 +340,7 @@ class FiniteDifferenceSolver
             std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jifield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jextfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
             std::unique_ptr<amrex::MultiFab> const& rhofield,
             std::unique_ptr<amrex::MultiFab> const& Pefield,

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H

@@ -46,6 +46,20 @@ public:
 
     void InitData ();
 
+    /**
+     * \brief
+     * Function to evaluate the external current expressions and populate the
+     * 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 (
+        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
+    );
+
     /**
      * \brief
      * Function to calculate the total current based on Ampere's law while
@@ -133,10 +147,19 @@ public:
     std::unique_ptr<amrex::Parser> m_resistivity_parser;
     amrex::ParserExecutor<1> m_eta;
 
+    /** External current */
+    std::string m_Jx_ext_grid_function = "0.0";
+    std::string m_Jy_ext_grid_function = "0.0";
+    std::string m_Jz_ext_grid_function = "0.0";
+    std::array< std::unique_ptr<amrex::Parser>, 3> m_J_external_parser;
+    std::array< amrex::ParserExecutor<4>, 3> m_J_external;
+    bool m_external_field_has_time_dependence = false;
+
     // Declare multifabs specifically needed for the hybrid-PIC model
     amrex::Vector<            std::unique_ptr<amrex::MultiFab>      > rho_fp_temp;
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_ampere;
+    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_external;
     amrex::Vector<            std::unique_ptr<amrex::MultiFab>      > electron_pressure_fp;
 
     // Helper functions to retrieve hybrid-PIC multifabs