From aa01df3e5db55cc26599f45df15f2069f30da17f Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 09:34:02 -0800 Subject: [PATCH 1/7] Fixing bug in hyper-resistivity calculation which had missing terms in vector laplacian evaluation. Additioanally fixing a staggering bug for density calculation in RZ. Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- .../HybridPICSolveE.cpp | 43 +++++++++++++++---- 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 47e45bbe753..0d844a5602a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -586,8 +586,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; + Real jr_val = 0._rt; if (solve_for_Faraday && resistivity_has_J_dependence) { - const Real jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); + jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); const Real jt_val = Interp(Jt, Jt_stag, Er_stag, coarsen, i, j, 0, 0); const Real jz_val = Interp(Jz, Jz_stag, Er_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); @@ -613,7 +614,8 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // 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); + 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,13 +635,14 @@ 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; + Real jt_val = 0._rt; if (solve_for_Faraday && resistivity_has_J_dependence) { const Real jr_val = Interp(Jr, Jr_stag, Et_stag, coarsen, i, j, 0, 0); - const Real jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); + jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); const Real jz_val = Interp(Jz, Jz_stag, Et_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); } @@ -659,7 +662,15 @@ 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) { + // Do not apply hyper-resistivity at r=0 + auto nabla2Jt = 0._rt; + if (r > 0.5_rt*dr) { + 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 +708,15 @@ 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) { - auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); + // r on nodal point (Jz is nodal in r) + Real const r = rmin + i*dr; + + auto nabla2Jz = 0._rt; + 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) + + T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); + } + Ez(i, j, 0) -= eta_h * nabla2Jz; } } @@ -918,7 +937,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 +979,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 +1021,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; } } From 829dd3126a2d48c0f2dec452f89e960c58196499 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 09:46:03 -0800 Subject: [PATCH 2/7] Removing trailing whitespace. Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 0d844a5602a..52fefd14318 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -716,7 +716,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( nabla2Jz = T_Algo::Dr_rDr_over_r(Jz, r, dr, coefs_r, n_coefs_r, i, j, 0, 0) + T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); } - + Ez(i, j, 0) -= eta_h * nabla2Jz; } } From d988e07e698c8947fd660d95f23a92cb18fdde6f Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 10:41:03 -0800 Subject: [PATCH 3/7] Changing location of current evaluation for hyper-resistivity to break dependencies on flags for standard resistivity. Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- .../HybridPICSolveE.cpp | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 52fefd14318..7788523076d 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -586,9 +586,8 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; - Real jr_val = 0._rt; if (solve_for_Faraday && resistivity_has_J_dependence) { - jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); + const Real jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); const Real jt_val = Interp(Jt, Jt_stag, Er_stag, coarsen, i, j, 0, 0); const Real jz_val = Interp(Jz, Jz_stag, Er_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); @@ -612,8 +611,8 @@ 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; - + 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; @@ -639,10 +638,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; - Real jt_val = 0._rt; if (solve_for_Faraday && resistivity_has_J_dependence) { const Real jr_val = Interp(Jr, Jr_stag, Et_stag, coarsen, i, j, 0, 0); - jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); + const Real jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); const Real jz_val = Interp(Jz, Jz_stag, Et_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); } @@ -663,12 +661,10 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( if (solve_for_Faraday) { Et(i, j, 0) += eta(rho_val, jtot_val) * Jt(i, j, 0); } if (include_hyper_resistivity_term) { - // Do not apply hyper-resistivity at r=0 - auto nabla2Jt = 0._rt; - if (r > 0.5_rt*dr) { - 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); - } + 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; } }, From c50a7720d2a3010814c52a0c025fc84731ac18a1 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 10:52:34 -0800 Subject: [PATCH 4/7] Update Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> --- .../FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 7788523076d..545c78e069f 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -707,10 +707,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // r on nodal point (Jz is nodal in r) Real const r = rmin + i*dr; - auto nabla2Jz = 0._rt; + 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) - + T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); + 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; From f6fe9dca8215a22f03f30f940c4a4251c0d76d56 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 10:55:58 -0800 Subject: [PATCH 5/7] adding missing semicolon Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 545c78e069f..2047e87b696 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -709,7 +709,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( 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) + 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; From 7e4022c2f173f3e625375521a87414038c10a307 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 11:49:47 -0800 Subject: [PATCH 6/7] Resetting Ohm's Law EM modes benchmarks for CI testing. Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- Examples/Tests/ohm_solver_em_modes/analysis_rz.py | 2 +- .../test_rz_ohm_solver_em_modes_picmi.json | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Examples/Tests/ohm_solver_em_modes/analysis_rz.py b/Examples/Tests/ohm_solver_em_modes/analysis_rz.py index 841e1177630..9c76265841c 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 } } From bb39e1b8a60d1ab09a0929655b6219e93fc42626 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 4 Feb 2025 11:53:35 -0800 Subject: [PATCH 7/7] Adding commas to array values for CI testing. Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com> --- Examples/Tests/ohm_solver_em_modes/analysis_rz.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Tests/ohm_solver_em_modes/analysis_rz.py b/Examples/Tests/ohm_solver_em_modes/analysis_rz.py index 9c76265841c..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([59.23850009 19.26746169 92.65794174 10.83627164]) + amps, np.array([59.23850009, 19.26746169, 92.65794174, 10.83627164]) )