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>

1 files changed, 58 insertions, 0 deletions

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.