Skip to content

Commit

Permalink
Fixing bug in hyper-resistivity calculation which had missing terms i… (
Browse files Browse the repository at this point in the history 
#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]>
  • Loading branch information
@clarkse @roelof-groenewald
clarkse and roelof-groenewald authored Feb 4, 2025
1 parent 93466dd commit 12269a0
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 16 deletions.
2 changes: 1 addition & 1 deletion Examples/Tests/ohm_solver_em_modes/analysis_rz.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
)
14 changes: 7 additions & 7 deletions Regression/Checksum/benchmarks_json/test_rz_ohm_solver_em_modes_picmi.json
View file
Original file line number Diff line number Diff line change
@@ -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
}
}
36 changes: 28 additions & 8 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
},
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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;
}
}
Expand Down Expand Up @@ -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;
}
},
Expand Down Expand Up @@ -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;
}
},
Expand Down Expand Up @@ -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;
}
}
Expand Down

0 comments on commit 12269a0

Please sign in to comment.