Adding hyper-resistivity to generalized ohms law hybrid solver. (#4772)

* Adding hyper-resistivity to generalized ohms law hybrid solver.

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

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

* Apply suggestions from code review

Style changes to clean up Laplacian operators.

Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>

* Removed extra arguments for triggering calculation of hyper_resistivity.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>

Python/pywarpx/picmi.py

@@ -1142,7 +1142,8 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
         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,
+                 n_floor=None, plasma_resistivity=None,
+                 plasma_hyper_resistivity=None, substeps=None,
                  Jx_external_function=None, Jy_external_function=None,
                  Jz_external_function=None, **kw):
         self.grid = grid
@@ -1153,6 +1154,7 @@ def __init__(self, grid, Te=None, n0=None, gamma=None,
         self.gamma = gamma
         self.n_floor = n_floor
         self.plasma_resistivity = plasma_resistivity
+        self.plasma_hyper_resistivity = plasma_hyper_resistivity
 
         self.substeps = substeps
 
@@ -1187,6 +1189,7 @@ def solver_initialize_inputs(self):
             'plasma_resistivity(rho,J)',
             pywarpx.my_constants.mangle_expression(self.plasma_resistivity, self.mangle_dict)
         )
+        pywarpx.hybridpicmodel.plasma_hyper_resistivity = self.plasma_hyper_resistivity
         pywarpx.hybridpicmodel.substeps = self.substeps
         pywarpx.hybridpicmodel.__setattr__(
             'Jx_external_grid_function(x,y,z,t)',

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H

@@ -100,6 +100,25 @@ struct CartesianYeeAlgorithm {
 #endif
     }
 
+    /**
+     * Perform second derivative along x on a cell-centered grid, from a cell-centered field `F`*/
+    template< typename T_Field>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Dxx (
+        T_Field const& F,
+        amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
+        int const i, int const j, int const k, int const ncomp=0 ) {
+
+        using namespace amrex;
+#if (defined WARPX_DIM_1D_Z)
+        amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
+        return 0._rt; // 1D Cartesian: derivative along x is 0
+#else
+        amrex::Real const inv_dx2 = coefs_x[0]*coefs_x[0];
+        return inv_dx2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) );
+#endif
+    }
+
     /**
      * Perform derivative along y on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -147,6 +166,25 @@ struct CartesianYeeAlgorithm {
 #endif
     }
 
+    /**
+     * Perform derivative along y on a nodal grid, from a cell-centered field `F`*/
+    template< typename T_Field>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Dyy (
+        T_Field const& F,
+        amrex::Real const * const coefs_y, int const /*n_coefs_y*/,
+        int const i, int const j, int const k, int const ncomp=0 ) {
+
+        using namespace amrex;
+#if defined WARPX_DIM_3D
+        Real const inv_dy2 = coefs_y[0]*coefs_y[0];
+        return inv_dy2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
+#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
+        amrex::ignore_unused(F, coefs_y, i, j, k, ncomp);
+        return 0._rt; // 1D and 2D Cartesian: derivative along y is 0
+#endif
+    }
+
     /**
      * Perform derivative along z on a cell-centered grid, from a nodal field `F`*/
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -186,6 +224,26 @@ struct CartesianYeeAlgorithm {
 #endif
     }
 
+    /**
+     * Perform second derivative along z on a cell-centered field `F`*/
+    template< typename T_Field>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Dzz (
+        T_Field const& F,
+        amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
+        int const i, int const j, int const k, int const ncomp=0 ) {
+
+        using namespace amrex;
+        Real const inv_dz2 = coefs_z[0]*coefs_z[0];
+#if defined WARPX_DIM_3D
+        return inv_dz2*( F(i,j,k-1,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j,k+1,ncomp) );
+#elif (defined WARPX_DIM_XZ)
+        return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
+#elif (defined WARPX_DIM_1D_Z)
+        return inv_dz2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) );
+#endif
+    }
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H

@@ -134,6 +134,24 @@ struct CylindricalYeeAlgorithm {
         return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) );
     }
 
+    /** Applies the differential operator `1/r * d(r * dF/dr)/dr`,
+     * where `F` is on a *cell-centered* or a nodal grid in `r`
+     * The input parameter `r` is given at the cell-centered position */
+    template< typename T_Field>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Dr_rDr_over_r (
+        T_Field const& F,
+        amrex::Real const r, amrex::Real const dr,
+        amrex::Real const * const coefs_r, int const /*n_coefs_r*/,
+        int const i, int const j, int const k, int const comp ) {
+
+        using namespace amrex;
+
+        Real const inv_dr2 = coefs_r[0]*coefs_r[0];
+        return 1._rt/r * inv_dr2*( (r+0.5_rt*dr)*(F(i+1,j,k,comp) - F(i,j,k,comp))
+                                 - (r-0.5_rt*dr)*(F(i,j,k,comp) - F(i-1,j,k,comp)) );
+    }
+
     /**
      * Perform derivative along z on a cell-centered grid, from a nodal field `F` */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -162,6 +180,21 @@ struct CylindricalYeeAlgorithm {
         return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) );
     }
 
+    /**
+     * Perform second derivative along z on a cell-centered field `F`*/
+    template< typename T_Field>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Dzz (
+        T_Field const& F,
+        amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
+        int const i, int const j, int const k, int const ncomp=0 ) {
+
+        using namespace amrex;
+        Real const inv_dz2 = coefs_z[0]*coefs_z[0];
+
+        return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
+    }
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H

@@ -175,6 +175,9 @@ public:
     amrex::ParserExecutor<2> m_eta;
     bool m_resistivity_has_J_dependence = false;
 
+    /** Plasma hyper-resisitivity */
+    amrex::Real m_eta_h = 0.0;
+
     /** External current */
     std::string m_Jx_ext_grid_function = "0.0";
     std::string m_Jy_ext_grid_function = "0.0";

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp

@@ -40,6 +40,8 @@ void HybridPICModel::ReadParameters ()
     pp_hybrid.query("plasma_resistivity(rho,J)", m_eta_expression);
     utils::parser::queryWithParser(pp_hybrid, "n_floor", m_n_floor);
 
+    utils::parser::queryWithParser(pp_hybrid, "plasma_hyper_resistivity", m_eta_h);
+
     // convert electron temperature from eV to J
     m_elec_temp *= PhysConst::q_e;
 

Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp

@@ -375,7 +375,7 @@ void FiniteDifferenceSolver::HybridPICSolveE (
     std::unique_ptr<amrex::MultiFab> const& Pefield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
     int lev, HybridPICModel const* hybrid_model,
-    const bool include_resistivity_term )
+    const bool include_resistivity_term)
 {
     // Select algorithm (The choice of algorithm is a runtime option,
     // but we compile code for each algorithm, using templates)
@@ -432,9 +432,12 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
 
     // get hybrid model parameters
     const auto eta = hybrid_model->m_eta;
+    const auto eta_h = hybrid_model->m_eta_h;
     const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e;
     const auto resistivity_has_J_dependence = hybrid_model->m_resistivity_has_J_dependence;
 
+    const bool include_hyper_resistivity_term = (eta_h > 0.0) && include_resistivity_term;
+
     // Index type required for interpolating fields from their respective
     // staggering to the Ex, Ey, Ez locations
     amrex::GpuArray<int, 3> const& Er_stag = hybrid_model->Ex_IndexType;
@@ -612,6 +615,14 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
 
                 // Add resistivity only if E field value is used to update B
                 if (include_resistivity_term) { Er(i, j, 0) += eta(rho_val, jtot_val) * Jr(i, j, 0); }
+
+                if (include_hyper_resistivity_term) {
+                    // r on cell-centered point (Jr is cell-centered in r)
+                    Real const r = rmin + (i + 0.5_rt)*dr;
+
+                    auto nabla2Jr = T_Algo::Dr_rDr_over_r(Jr, r, dr, coefs_r, n_coefs_r, i, j, 0, 0);
+                    Er(i, j, 0) -= eta_h * nabla2Jr;
+                }
             },
 
             // Et calculation
@@ -655,16 +666,18 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
 
                 // Add resistivity only if E field value is used to update B
                 if (include_resistivity_term) { Et(i, j, 0) += eta(rho_val, jtot_val) * Jt(i, j, 0); }
+
+                // Note: Hyper-resisitivity should be revisited here when modal decomposition is implemented
             },
 
             // Ez calculation
-            [=] AMREX_GPU_DEVICE (int i, int j, int k){
+            [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
 #ifdef AMREX_USE_EB
                 // Skip field solve if this cell is fully covered by embedded boundaries
                 if (lz(i,j,0) <= 0) { return; }
 #endif
                 // Interpolate to get the appropriate charge density in space
-                Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0);
+                Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, 0, 0);
 
                 // Interpolate current to appropriate staggering to match E field
                 Real jtot_val = 0._rt;
@@ -679,15 +692,20 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
                 if (rho_val < rho_floor) { rho_val = rho_floor; }
 
                 // Get the gradient of the electron pressure
-                auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, k, 0);
+                auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, 0, 0);
 
                 // interpolate the nodal neE values to the Yee grid
-                auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, k, 2);
+                auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, 0, 2);
 
-                Ez(i, j, k) = (enE_z - grad_Pe) / rho_val;
+                Ez(i, j, 0) = (enE_z - grad_Pe) / rho_val;
 
                 // Add resistivity only if E field value is used to update B
-                if (include_resistivity_term) { Ez(i, j, k) += eta(rho_val, jtot_val) * Jz(i, j, k); }
+                if (include_resistivity_term) { Ez(i, j, 0) += eta(rho_val, jtot_val) * Jz(i, j, 0); }
+
+                if (include_hyper_resistivity_term) {
+                    auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0);
+                    Ez(i, j, 0) -= eta_h * nabla2Jz;
+                }
             }
         );
 
@@ -726,9 +744,12 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
 
     // get hybrid model parameters
     const auto eta = hybrid_model->m_eta;
+    const auto eta_h = hybrid_model->m_eta_h;
     const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e;
     const auto resistivity_has_J_dependence = hybrid_model->m_resistivity_has_J_dependence;
 
+    const bool include_hyper_resistivity_term = (eta_h > 0.) && include_resistivity_term;
+
     // Index type required for interpolating fields from their respective
     // staggering to the Ex, Ey, Ez locations
     amrex::GpuArray<int, 3> const& Ex_stag = hybrid_model->Ex_IndexType;
@@ -904,6 +925,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
 
                 // Add resistivity only if E field value is used to update B
                 if (include_resistivity_term) { Ex(i, j, k) += eta(rho_val, jtot_val) * Jx(i, j, k); }
+
+                if (include_hyper_resistivity_term) {
+                    auto nabla2Jx = T_Algo::Dxx(Jx, coefs_x, n_coefs_x, i, j, k);
+                    Ex(i, j, k) -= eta_h * nabla2Jx;
+                }
             },
 
             // Ey calculation
@@ -943,6 +969,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
 
                 // Add resistivity only if E field value is used to update B
                 if (include_resistivity_term) { Ey(i, j, k) += eta(rho_val, jtot_val) * Jy(i, j, k); }
+
+                if (include_hyper_resistivity_term) {
+                    auto nabla2Jy = T_Algo::Dyy(Jy, coefs_y, n_coefs_y, i, j, k);
+                    Ey(i, j, k) -= eta_h * nabla2Jy;
+                }
             },
 
             // Ez calculation
@@ -976,6 +1007,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
 
                 // Add resistivity only if E field value is used to update B
                 if (include_resistivity_term) { Ez(i, j, k) += eta(rho_val, jtot_val) * Jz(i, j, k); }
+
+                if (include_hyper_resistivity_term) {
+                    auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, k);
+                    Ez(i, j, k) -= eta_h * nabla2Jz;
+                }
             }
         );
 

Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp

@@ -154,8 +154,7 @@ void WarpX::HybridPICEvolveFields ()
     // Update the E field to t=n+1 using the extrapolated J_i^n+1 value
     m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
     m_hybrid_pic_model->HybridPICSolveE(
-        Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths,
-        false
+        Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, false
     );
     FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);