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, 32 insertions, 0 deletions

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):
         '''