Fixing bug in hyper-resistivity calculation which had missing terms i… (#5638)

…n vector laplacian evaluation. Additioanally fixing a staggering bug
for density calculation in RZ.

---------

Signed-off-by: S. Eric Clark <[email protected]>
Co-authored-by: Roelof Groenewald <[email protected]>

Author: clarkse

Examples/Tests/ohm_solver_em_modes/analysis_rz.py

@@ -179,5 +179,5 @@ def process(it):
     amps = np.abs(F_kw[2, 1, len(kz) // 2 - 2 : len(kz) // 2 + 2])
     print("Amplitude sample: ", amps)
     assert np.allclose(
-        amps, np.array([61.02377286, 19.80026021, 100.47687017, 10.83331295])
+        amps, np.array([59.23850009, 19.26746169, 92.65794174, 10.83627164])
     )

Regression/Checksum/benchmarks_json/test_rz_ohm_solver_em_modes_picmi.json

@@ -1,12 +1,12 @@
 {
   "lev=0": {},
   "ions": {
-    "particle_momentum_x": 5.0438993756415296e-17,
-    "particle_momentum_y": 5.0444406612873916e-17,
-    "particle_momentum_z": 5.0519292431385393e-17,
-    "particle_position_x": 143164.41713467025,
-    "particle_position_y": 143166.51845281923,
-    "particle_theta": 2573261.8729711357,
-    "particle_weight": 8.128680645366887e+18
+    "particle_momentum_x": 5.043784704795177e-17,
+    "particle_momentum_y": 5.0444695502620235e-17,
+    "particle_momentum_z": 5.05193106847111e-17,
+    "particle_position_x": 143164.53685279266,
+    "particle_position_y": 143166.5185853012,
+    "particle_theta": 2573262.446840369,
+    "particle_weight": 8.128680645366886e+18
   }
 }

Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp

@@ -611,9 +611,10 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
 
                 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);
+                    const Real r = rmin + (i + 0.5_rt)*dr;
+                    const Real jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0);
+                    auto nabla2Jr = T_Algo::Dr_rDr_over_r(Jr, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)
+                        + T_Algo::Dzz(Jr, coefs_z, n_coefs_z, i, j, 0, 0) - jr_val/(r*r);
                     Er(i, j, 0) -= eta_h * nabla2Jr;
                 }
             },
@@ -633,7 +634,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
                 }
 
                 // Interpolate to get the appropriate charge density in space
-                Real rho_val = Interp(rho, nodal, Er_stag, coarsen, i, j, 0, 0);
+                Real rho_val = Interp(rho, nodal, Et_stag, coarsen, i, j, 0, 0);
 
                 // Interpolate current to appropriate staggering to match E field
                 Real jtot_val = 0._rt;
@@ -659,7 +660,13 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
                 // Add resistivity only if E field value is used to update B
                 if (solve_for_Faraday) { 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
+                if (include_hyper_resistivity_term) {
+                    const Real jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0);
+                    auto nabla2Jt = T_Algo::Dr_rDr_over_r(Jt, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)
+                        + T_Algo::Dzz(Jt, coefs_z, n_coefs_z, i, j, 0, 0) - jt_val/(r*r);
+
+                    Et(i, j, 0) -= eta_h * nabla2Jt;
+                }
             },
 
             // Ez calculation
@@ -697,7 +704,14 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
                 if (solve_for_Faraday) { Ez(i, j, 0) += eta(rho_val, jtot_val) * Jz(i, j, 0); }
 
                 if (include_hyper_resistivity_term) {
+                    // r on nodal point (Jz is nodal in r)
+                    Real const r = rmin + i*dr;
+
                     auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0);
+                    if (r > 0.5_rt*dr) {
+                        nabla2Jz += T_Algo::Dr_rDr_over_r(Jz, r, dr, coefs_r, n_coefs_r, i, j, 0, 0);
+                    }
+
                     Ez(i, j, 0) -= eta_h * nabla2Jz;
                 }
             }
@@ -918,7 +932,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
                 if (solve_for_Faraday) { 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);
+                    auto nabla2Jx = T_Algo::Dxx(Jx, coefs_x, n_coefs_x, i, j, k)
+                        + T_Algo::Dyy(Jx, coefs_y, n_coefs_y, i, j, k)
+                        + T_Algo::Dzz(Jx, coefs_z, n_coefs_z, i, j, k);
                     Ex(i, j, k) -= eta_h * nabla2Jx;
                 }
             },
@@ -958,7 +974,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
                 if (solve_for_Faraday) { 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);
+                    auto nabla2Jy = T_Algo::Dxx(Jy, coefs_x, n_coefs_x, i, j, k)
+                        + T_Algo::Dyy(Jy, coefs_y, n_coefs_y, i, j, k)
+                        + T_Algo::Dzz(Jy, coefs_z, n_coefs_z, i, j, k);
                     Ey(i, j, k) -= eta_h * nabla2Jy;
                 }
             },
@@ -998,7 +1016,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
                 if (solve_for_Faraday) { 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);
+                    auto nabla2Jz = T_Algo::Dxx(Jz, coefs_x, n_coefs_x, i, j, k)
+                        + T_Algo::Dyy(Jz, coefs_y, n_coefs_y, i, j, k)
+                        + T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, k);
                     Ez(i, j, k) -= eta_h * nabla2Jz;
                 }
             }