diff --git a/Examples/Tests/ohm_solver_em_modes/analysis_rz.py b/Examples/Tests/ohm_solver_em_modes/analysis_rz.py index 841e1177630..7cd5086c408 100755 --- a/Examples/Tests/ohm_solver_em_modes/analysis_rz.py +++ b/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]) ) diff --git a/Regression/Checksum/benchmarks_json/test_rz_ohm_solver_em_modes_picmi.json b/Regression/Checksum/benchmarks_json/test_rz_ohm_solver_em_modes_picmi.json index ec1b6272092..feca88922e2 100644 --- a/Regression/Checksum/benchmarks_json/test_rz_ohm_solver_em_modes_picmi.json +++ b/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 } } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 47e45bbe753..2047e87b696 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/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; } }