RZ spectral current correction and Galilean (#1216)

* Added stub for current correction in RZ spectral solver

* Fixed comments in RZ spectral for current correction stub

* Modified automated test for Galilean PSATD (#1033)

* Impemented current correction in RZ spectral

* Implementation Galilean version of RZ spectral solver

* For RZ spectral, do forward and backward transform with current correction

* Big fix in DivEFunctor.cpp for RZ spectral

* Added RZ rho diagnostic for saving the modes

* Implemented fft_periodic_single_box for RZ spectral

* Moved routines from SpectralSolverRZ.H to .cpp

* Added hook for VayDeposition in GalileanPsatdAlgorithmRZ

* Bug fix in DivEFunctor

* Fixes and cleanup for GalileanPsatdAlgorithmRZ

* Fix line spacing in SpectralSolverRZ.H

* Fix factor 1/2 in update of Ep_m

* Fix factor 1/2 in update of Em_m

* Fix sign error in current correction in GalileanPsatdAlgorithmRZ.cpp

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Add Langmuir RZ PSATD test with current correction

* Add Galilean tests with/without current correction

* For RZ psatd, simplified copy for forward transform

* Added GalileanPsatdAlgorithmRZ.cpp to CMakeLists

* Minor cleanup in RZ spectral solver

* In GalileanPsatdAlgorithmRZ.cpp use member initialization for m_v_galilean

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Added some _rt to GalileanPsatdAlgorithmRZ.cpp

Co-authored-by: Olga Shapoval <30510597+oshapoval@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>

1 files changed, 18 insertions, 8 deletions

diff --git a/Examples/Tests/galilean/analysis_2d.py b/Examples/Tests/galilean/analysis_2d.py
index e51689969..ac4ea2efd 100755
--- a/Examples/Tests/galilean/analysis_2d.py
+++ b/Examples/Tests/galilean/analysis_2d.py
@@ -27,6 +27,7 @@ filename = sys.argv[1]
 # Parse test name
 averaged = True if re.search( 'averaged', filename ) else False
 current_correction = True if re.search( 'current_correction', filename ) else False
+dims_RZ  = True if re.search('rz', filename) else False
 
 ds = yt.load( filename )
 
@@ -38,15 +39,24 @@ if (averaged):
     # energyE_ref was calculated with Galilean PSATD method (v_galilean = (0,0,0.99498743710662))
     energyE_ref = 26913.546573259937
     tolerance_rel = 1e-5
-elif (current_correction):
+elif (not dims_RZ and not current_correction):
+    # energyE_ref was calculated with standard PSATD method (v_galilean = (0.,0.,0.))
+    energyE_ref = 38362.88743899688
+    tolerance_rel = 1e-8
+elif (not dims_RZ and current_correction):
     # energyE_ref was calculated with standard PSATD method (v_galilean = (0.,0.,0.)):
-    # difference with respect to reference energy below due to absence of filter
+    # difference with respect to reference energy above due to absence of real-space filter
     energyE_ref = 745973.5742103161
-    tolerance_rel = 1e-8;
-else:
+    tolerance_rel = 1e-8
+elif (dims_RZ and not current_correction):
     # energyE_ref was calculated with standard PSATD method (v_galilean = (0.,0.,0.))
-    energyE_ref = 38362.88743899688
-    tolerance_rel = 1e-8;
+    energyE_ref = 178013.54481470847
+    tolerance_rel = 1e-8
+elif (dims_RZ and current_correction):
+    # energyE_ref was calculated with standard PSATD method (v_galilean = (0.,0.,0.))
+    # difference with respect to reference energy above due to absence of k-space filter
+    energyE_ref = 10955626.277865639
+    tolerance_rel = 1e-8
 
 energyE = np.sum(scc.epsilon_0/2*(Ex**2+Ey**2+Ez**2))
 
@@ -59,9 +69,9 @@ assert( error_rel < tolerance_rel )
 
 # Check charge conservation (relative L-infinity norm of error) with current correction
 if current_correction:
-    rho  = ds.index.grids[0]['boxlib', 'rho' ].squeeze().v
     divE = ds.index.grids[0]['boxlib', 'divE'].squeeze().v
-    error_rel = np.amax( np.abs( divE - rho/scc.epsilon_0 ) ) / np.amax( np.abs( rho/scc.epsilon_0 ) )
+    rho  = ds.index.grids[0]['boxlib', 'rho' ].squeeze().v / scc.epsilon_0
+    error_rel = np.amax(np.abs(divE - rho)) / max(np.amax(divE), np.amax(rho))
     tolerance = 1e-9
     print("Check charge conservation:")
     print("error_rel = {}".format(error_rel))