Ohm's law solver (hybrid kinetic-fluid extension) (#3665)

* Add "None" as an option for the Maxwell solver

* fixed some of the reasons for failing CI tests

* no longer pass `do_electrostatic` to `GuardCellManager`

* renamed `MaxwellSolverAlgo` to `ElectromagneticSolverAlgo`

* rename `do_electrostatic` to `electrostatic_solver_id`

* rename `maxwell_solver_id` to `electromagnetic_solver_id`

* start of hybrid solver logic

* changes requested during PR review

* remove `do_no_deposit` from tests without field evolution

* added `HybridSolveE.cpp`

* bulk of the hybrid solver implementation

* mostly reproduce 1d cold ion mirror results

* ion Bernstein modes reproduced with this version

* fix bug with reduced diagnostic FieldProbe in 1d

* added hybrid solver installation to PICMI and added example script generating ion-Bernstein modes

* enable the use of `FieldProbe` default parameter values

* use default field-probe values

* steady progress

* add `do_not_push` flag to picmi

* possibly use nodal fields

* added extra multifab for current calculated from curl B

* added `CalculateTotalCurrent` functions

* rewrote implementation to calculate J x curl B on a nodal grid first

* fill boundary for auxiliary MFs used during hybrid E solve

* properly handle nonzero resistivity

* updated hybrid model example

* clean up example scripts

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fixed invalid memory access for GPU and other code cleanups

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* refinements on the example scripts

* added ion beam instability example

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* added EM modes and ion beam examples to CI test suite

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* started docs section on the hybrid model

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* more progress on documentation

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* added ion Landau damping verification test / example

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* add checksum benchmark for Landau damping example

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* added fields.py wrapper to access total current density in hybrid case

* refactored the charge deposition fix to be performed with the field data rather than individual particles

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* also correct current density at PEC boundary

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* made resistivity a parsed function of `rho`

* work on PEC boundary condition

* corrections pointed out during code review

* fix build fails due to unused variables

* fix issue with GPU builds

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* actually apply rho boundary correction in EM case

* take one sided derivative at PEC boundary when calculating div Pe

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* actually apply rho boundary correction in EM case

* removed specific treatment of E-field on PEC boundary for Ohm's law solver

* first round of CI fixes

* second round of CI fixes

* added description of deposition logic with PEC boundaries to documentation

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* third round of CI fixing

* move J and rho boundary handling to after `SyncRho` and `SyncCurrent` calls

* properly order the application of boundary conditions on rho, for electrostatic simulations

* fourth round of CI fixing

* moved calculation of total current (Ampere's law) to seperate function

* add random seed specification to `picmi`

* code clean-up -> renamed hybrid model to hybrid-PIC model

* added magnetic reconnection example

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* code cleanup & benchmark updates

* update PICMI class name for hybrid solver to `HybridPICSolver`

* don't apply J field boundary in

* don't apply J field boundary in `MultiParticleContainer::DepositCurrent`

* apply changes requested during code review

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Loosen tolerance on failing CI test

* Removed unused variable

* code cleanup: make use of `MultiParticleContainer::DepositCurrent` in `AddSpaceChargeFieldLabFrame`

* switch to using a rho_fp_temp multifab for old and averaged charge density field, also no longer require particles to move only one cell per step

* use `ablastr::coarsen::sample` namespace in `HybridPICSolveECartesian`

* switched to using `MultiFab::LinComb` instead of self written GPU kernels to calculated averaged or extrapolated current density

* add verbosity flag for the Ohm solver tests

* deal with fine versus coarse patches

* add theoretical instability growth / damping rates to hybrid-PIC examples

* update ion-Bernstein mode plot in documentation

* move the `ApplyRhofieldBoundary` call to after `SumBoundary`

* use a uniform calculation for the number of cells a given index is from the boundary

* remove unused variable

* limit number of ghost cells updated during PEC BC application

* the number of ghost cells to consider depends on whether the field is nodal or not

* attempt 1 to fix failing CI tests

* attempt 2 to fix failing CI tests and code cleanup

* attempt 3 to fix failing CI tests

* attempt 4 to fix failing CI tests and docs cleanup

* switched to using bibtex citations

* move hybrid solver input parameters documentation to its own section

* clean up ion beam instability analysis script

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* add inline comments describing the meaning of each argument for the `amrex::MultiFab::LinComb` calls used

* make `HybridPICSolver` a child class of `picmistandard.base._ClassWithInit`

* apply changes requested during code review

* add warning about using hybrid-PIC solver with Esirkepov current deposition

* add Stanier 2020 reference to recommend linear particles with hybrid-PIC

* add call to FillBoundary for `current_fp` - needed for accurate interpolation to nodal grid

* changes requested from code review

* Apply suggestions from code review

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* include physics accuracy check for ion beam instability; switch CI tests to use direct current deposition

* reset benchmark values after switching to direct current deposition

* update ion beam instability benchmark

* minor changes requested during code review

* remove guard cells for `enE_nodal_mf` as well as corresponding `FillBoundary` call

* refactor: moved hybrid-PIC specific multifabs and `CalculateElectronPressure()` to `HybridPICModel`

* add assert that load balancing is not used with the hybrid-PIC solver since the new multifabs are not yet properly redistributed

* move `CalculateCurrentAmpere` to `HybridPICModel`

* refactor: moved `HybridPICSolveE` to `HybridPICModel` class

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

6 files changed, 125 insertions, 0 deletions

diff --git a/Python/pywarpx/HybridPICModel.py b/Python/pywarpx/HybridPICModel.py
new file mode 100644
index 000000000..e21bba4a2
--- /dev/null
+++ b/Python/pywarpx/HybridPICModel.py
@@ -0,0 +1,11 @@
+# Copyright 2023 The WarpX Community
+#
+# This file is part of WarpX.
+#
+# Authors: Roelof Groenewald (TAE Technologies)
+#
+# License: BSD-3-Clause-LBNL
+
+from .Bucket import Bucket
+
+hybridpicmodel = Bucket('hybrid_pic_model')
diff --git a/Python/pywarpx/WarpX.py b/Python/pywarpx/WarpX.py
index e8c170a76..79b248617 100644
--- a/Python/pywarpx/WarpX.py
+++ b/Python/pywarpx/WarpX.py
@@ -18,6 +18,7 @@ from .Constants import my_constants
 from .Diagnostics import diagnostics, reduced_diagnostics
 from .EB2 import eb2
 from .Geometry import geometry
+from .HybridPICModel import hybridpicmodel
 from .Interpolation import interpolation
 from .Langmuirwave import langmuirwave
 from .Lasers import lasers, lasers_list
@@ -43,6 +44,7 @@ class WarpX(Bucket):
         argv += amr.attrlist()
         argv += amrex.attrlist()
         argv += geometry.attrlist()
+        argv += hybridpicmodel.attrlist()
         argv += boundary.attrlist()
         argv += algo.attrlist()
         argv += langmuirwave.attrlist()
diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
index 5000180f3..1e80ccfdf 100644
--- a/Python/pywarpx/__init__.py
+++ b/Python/pywarpx/__init__.py
@@ -13,6 +13,7 @@ from .Constants import my_constants
 from .Diagnostics import diagnostics, reduced_diagnostics
 from .EB2 import eb2
 from .Geometry import geometry
+from .HybridPICModel import hybridpicmodel
 from .Interpolation import interpolation
 from .Langmuirwave import langmuirwave
 from .Lasers import lasers
diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py
index 02ca22afb..ebd18e5ee 100755
--- a/Python/pywarpx/_libwarpx.py
+++ b/Python/pywarpx/_libwarpx.py
@@ -200,6 +200,7 @@ class LibWarpX():
         self.libwarpx_so.warpx_getCurrentDensityCPLoVects_PML.restype = _LP_c_int
         self.libwarpx_so.warpx_getCurrentDensityFP_PML.restype = _LP_LP_c_real
         self.libwarpx_so.warpx_getCurrentDensityFPLoVects_PML.restype = _LP_c_int
+        self.libwarpx_so.warpx_getCurrentDensityFP_Ampere.restype = _LP_LP_c_real
         self.libwarpx_so.warpx_getChargeDensityCP.restype = _LP_LP_c_real
         self.libwarpx_so.warpx_getChargeDensityCPLoVects.restype = _LP_c_int
         self.libwarpx_so.warpx_getChargeDensityFP.restype = _LP_LP_c_real
@@ -1711,6 +1712,37 @@ class LibWarpX():
         except ValueError:
             raise Exception('PML not initialized')
 
+    def get_mesh_current_density_fp_ampere(self, level, direction, include_ghosts=True):
+        '''
+
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process calculated from the curl of B. This
+        quantity is calculated in the kinetic-fluid hybrid model to get the
+        electron current. This function returns the current density on the fine
+        patch for the given level.
+
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
+
+        Parameters
+        ----------
+
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
+
+        Returns
+        -------
+
+            A List of numpy arrays.
+
+        '''
+
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensityFP_Ampere, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('Current multifab not allocated.')
+
     def get_mesh_charge_density_cp(self, level, include_ghosts=True):
         '''
 
diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py
index 2577902c4..40ce6aa4a 100644
--- a/Python/pywarpx/fields.py
+++ b/Python/pywarpx/fields.py
@@ -897,6 +897,27 @@ def JzFPPMLWrapper(level=0, include_ghosts=False):
                             get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
+def JxFPAmpereWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=0,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_ampere,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
+                            level=level, include_ghosts=include_ghosts)
+
+def JyFPAmpereWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=1,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_ampere,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
+                            level=level, include_ghosts=include_ghosts)
+
+def JzFPAmpereWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=2,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_ampere,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
+                            level=level, include_ghosts=include_ghosts)
+
 def FFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
                             get_lovects=libwarpx.get_mesh_F_fp_lovects_pml,
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index 8e4860078..621757132 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -889,6 +889,7 @@ class Cartesian3DGrid(picmistandard.PICMI_Cartesian3DGrid):
         else:
             pywarpx.amr.max_level = 0
 
+
 class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver):
     """
     See `Input Parameters <https://warpx.readthedocs.io/en/latest/usage/parameters.html>`_ for more information.
@@ -998,6 +999,63 @@ class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver):
         pywarpx.warpx.pml_has_particles = self.pml_has_particles
         pywarpx.warpx.do_pml_j_damping = self.do_pml_j_damping
 
+
+class HybridPICSolver(picmistandard.base._ClassWithInit):
+    """
+    Hybrid-PIC solver based on Ohm's law.
+    See `Theory Section <https://warpx.readthedocs.io/en/latest/theory/kinetic_fluid_hybrid_model.html>`_ for more information.
+
+    Parameters
+    ----------
+    Te: float
+        Electron temperature in eV.
+
+    n0: float
+        Reference plasma density in m^-3.
+
+    gamma: float, default=3/2
+        Exponent in calculation of electron pressure.
+
+    n_floor: float, optional
+        Minimum density used in Ohm's law calculation.
+
+    plasma_resistivity: float or str
+        Value or expression to use for the plasma resistivity.
+
+    substeps: int, default=100
+        Number of substeps to take when updating the B-field.
+    """
+    def __init__(self, grid, Te=None, n0=None, gamma=None,
+                 n_floor=None, plasma_resistivity=None, substeps=None, **kw):
+        self.grid = grid
+        self.method = "hybrid"
+
+        self.Te = Te
+        self.n0 = n0
+        self.gamma = gamma
+        self.n_floor = n_floor
+        self.plasma_resistivity = plasma_resistivity
+
+        self.substeps = substeps
+
+        self.handle_init(kw)
+
+    def initialize_inputs(self):
+
+        self.grid.initialize_inputs()
+
+        pywarpx.algo.maxwell_solver = self.method
+
+        pywarpx.hybridpicmodel.elec_temp = self.Te
+        pywarpx.hybridpicmodel.n0_ref = self.n0
+        pywarpx.hybridpicmodel.gamma = self.gamma
+        pywarpx.hybridpicmodel.n_floor = self.n_floor
+        pywarpx.hybridpicmodel.__setattr__(
+            'plasma_resistivity(rho)', self.plasma_resistivity
+        )
+        pywarpx.hybridpicmodel.substeps = self.substeps
+
+
 class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):
     """
     See `Input Parameters <https://warpx.readthedocs.io/en/latest/usage/parameters.html>`_ for more information.