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

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
new file mode 100644
index 000000000..0793b3b27
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
@@ -0,0 +1,181 @@
+/* Copyright 2023 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Roelof Groenewald (TAE Technologies)
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_HYBRIDPICMODEL_H_
+#define WARPX_HYBRIDPICMODEL_H_
+
+#include "HybridPICModel_fwd.H"
+
+#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
+#include "Utils/Parser/ParserUtils.H"
+#include "Utils/WarpXConst.H"
+#include "Utils/WarpXProfilerWrapper.H"
+#include "WarpX.H"
+
+#include <AMReX_Array.H>
+#include <AMReX_REAL.H>
+
+
+/**
+ * \brief This class contains the parameters needed to evaluate hybrid field
+ * solutions (kinetic ions with fluid electrons).
+ */
+class HybridPICModel
+{
+public:
+    HybridPICModel (int nlevs_max); // constructor
+
+    /** Read user-defined model parameters. Called in constructor. */
+    void ReadParameters ();
+
+    /** 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,
+                           const 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);
+
+    /** Helper function to clear values from hybrid-PIC specific multifabs. */
+    void ClearLevel (int lev);
+
+    void InitData ();
+
+    /**
+     * \brief
+     * Function to calculate the total current based on Ampere's law while
+     * neglecting displacement current (J = curl x B). Used in the Ohm's law
+     * solver (kinetic-fluid hybrid model).
+     *
+     * \param[in] Bfield       Magnetic field from which the current is calculated.
+     * \param[in] edge_lengths Length of cell edges taking embedded boundaries into account
+     */
+    void CalculateCurrentAmpere (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths
+    );
+    void CalculateCurrentAmpere (
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        const int lev
+    );
+
+    /**
+     * \brief
+     * Function to update the E-field using Ohm's law (hybrid-PIC model).
+     */
+    void HybridPICSolveE (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
+        amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
+        DtType dt_type);
+    void HybridPICSolveE (
+        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& Bfield,
+        std::unique_ptr<amrex::MultiFab> const& rhofield,
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        const int lev, DtType dt_type);
+    void HybridPICSolveE (
+        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& Bfield,
+        std::unique_ptr<amrex::MultiFab> const& rhofield,
+        std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        const int lev, PatchType patch_type, DtType dt_type);
+
+    /**
+     * \brief
+     * Function to calculate the electron pressure at a given timestep type
+     * using the simulation charge density. Used in the Ohm's law solver
+     * (kinetic-fluid hybrid model).
+     */
+    void CalculateElectronPressure (                DtType a_dt_type);
+    void CalculateElectronPressure (const int lev,  DtType a_dt_type);
+
+    /**
+     * \brief Fill the electron pressure multifab given the kinetic particle
+     * charge density (and assumption of quasi-neutrality) using the user
+     * specified electron equation of state.
+     *
+     * \param[out] Pe   scalar electron pressure MultiFab at a given level
+     * \param[in] rhofield scalar ion chrge density Multifab at a given level
+     */
+    void FillElectronPressureMF (
+        std::unique_ptr<amrex::MultiFab> const& Pe,
+        amrex::MultiFab* const& rhofield );
+
+    // Declare variables to hold hybrid-PIC model parameters
+    /** Number of substeps to take when evolving B */
+    int m_substeps = 100;
+
+    /** Electron temperature in eV */
+    amrex::Real m_elec_temp;
+    /** Reference electron density */
+    amrex::Real m_n0_ref = 1.0;
+    /** Electron pressure scaling exponent */
+    amrex::Real m_gamma = 5.0/3.0;
+
+    /** Plasma density floor - if n < n_floor it will be set to n_floor */
+    amrex::Real m_n_floor = 1.0;
+
+    /** Plasma resistivity */
+    std::string m_eta_expression = "0.0";
+    std::unique_ptr<amrex::Parser> m_resistivity_parser;
+    amrex::ParserExecutor<1> m_eta;
+
+    // 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::unique_ptr<amrex::MultiFab>      > electron_pressure_fp;
+
+    // Helper functions to retrieve hybrid-PIC multifabs
+    amrex::MultiFab * get_pointer_current_fp_ampere  (int lev, int direction) const { return current_fp_ampere[lev][direction].get(); }
+    amrex::MultiFab * get_pointer_electron_pressure_fp  (int lev) const { return electron_pressure_fp[lev].get(); }
+
+    /** Gpu Vector with index type of the Jx multifab */
+    amrex::GpuArray<int, 3> Jx_IndexType;
+    /** Gpu Vector with index type of the Jy multifab */
+    amrex::GpuArray<int, 3> Jy_IndexType;
+    /** Gpu Vector with index type of the Jz multifab */
+    amrex::GpuArray<int, 3> Jz_IndexType;
+    /** Gpu Vector with index type of the Bx multifab */
+    amrex::GpuArray<int, 3> Bx_IndexType;
+    /** Gpu Vector with index type of the By multifab */
+    amrex::GpuArray<int, 3> By_IndexType;
+    /** Gpu Vector with index type of the Bz multifab */
+    amrex::GpuArray<int, 3> Bz_IndexType;
+    /** Gpu Vector with index type of the Ex multifab */
+    amrex::GpuArray<int, 3> Ex_IndexType;
+    /** Gpu Vector with index type of the Ey multifab */
+    amrex::GpuArray<int, 3> Ey_IndexType;
+    /** Gpu Vector with index type of the Ez multifab */
+    amrex::GpuArray<int, 3> Ez_IndexType;
+};
+
+/**
+ * \brief
+ * This struct contains only static functions to compute the electron pressure
+ * using the particle density at a given point and the user provided reference
+ * density and temperatures.
+ */
+struct ElectronPressure {
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real get_pressure (amrex::Real const n0,
+                                     amrex::Real const T0,
+                                     amrex::Real const gamma,
+                                     amrex::Real const rho) {
+        return n0 * T0 * pow((rho/PhysConst::q_e)/n0, gamma);
+    }
+};
+
+#endif // WARPX_HYBRIDPICMODEL_H_