Fix PSATD equations used with PML (#1513)

* Implement new PML PSATD equations

* Update CI test and benchmark

* Compute coefficients C1,...,C22 on the fly

* Add check on initial energy from diagnostics

1 files changed, 24 insertions, 18 deletions

diff --git a/Examples/Tests/PML/analysis_pml_psatd.py b/Examples/Tests/PML/analysis_pml_psatd.py
index a28541c53..d219fdd1d 100755
--- a/Examples/Tests/PML/analysis_pml_psatd.py
+++ b/Examples/Tests/PML/analysis_pml_psatd.py
@@ -7,7 +7,6 @@
 #
 # License: BSD-3-Clause-LBNL
 
-
 import sys
 import yt ; yt.funcs.mylog.setLevel(0)
 import numpy as np
@@ -20,7 +19,23 @@ filename = sys.argv[1]
 ############################
 ### INITIAL LASER ENERGY ###
 ############################
-energy_start = 9.1301289517e-08
+energy_start = 7.297301362346945e-08 # electromagnetic energy at iteration 50
+
+# Check consistency of field energy diagnostics with initial energy above
+ds = yt.load('pml_x_psatd_plt00050')
+all_data_level_0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Bx = all_data_level_0['boxlib', 'Bx'].v.squeeze()
+By = all_data_level_0['boxlib', 'By'].v.squeeze()
+Bz = all_data_level_0['boxlib', 'Bz'].v.squeeze()
+Ex = all_data_level_0['boxlib', 'Ex'].v.squeeze()
+Ey = all_data_level_0['boxlib', 'Ey'].v.squeeze()
+Ez = all_data_level_0['boxlib', 'Ez'].v.squeeze()
+energyE = np.sum(0.5 * scc.epsilon_0 * (Ex**2 + Ey**2 + Ez**2))
+energyB = np.sum(0.5 / scc.mu_0 * (Bx**2 + By**2 + Bz**2))
+energy_start_diags = energyE + energyB
+error = abs(energy_start - energy_start_diags) / energy_start
+tolerance = 1e-14
+assert (error < tolerance)
 
 ##########################
 ### FINAL LASER ENERGY ###
@@ -33,26 +48,17 @@ Bz = all_data_level_0['boxlib', 'Bz'].v.squeeze()
 Ex = all_data_level_0['boxlib', 'Ex'].v.squeeze()
 Ey = all_data_level_0['boxlib', 'Ey'].v.squeeze()
 Ez = all_data_level_0['boxlib', 'Ez'].v.squeeze()
-rho = all_data_level_0['boxlib','rho'].v.squeeze()
-divE = all_data_level_0['boxlib','divE'].v.squeeze()
-energyE = np.sum(scc.epsilon_0/2*(Ex**2+Ey**2+Ez**2))
-energyB = np.sum(1./scc.mu_0/2*(Bx**2+By**2+Bz**2))
+energyE = np.sum(0.5 * scc.epsilon_0 * (Ex**2 + Ey**2 + Ez**2))
+energyB = np.sum(0.5 / scc.mu_0 * (Bx**2 + By**2 + Bz**2))
 energy_end = energyE + energyB
 
-Reflectivity = energy_end/energy_start
-Reflectivity_theory = 1.3806831258153887e-06
-
-error_rel = abs(Reflectivity-Reflectivity_theory) / Reflectivity_theory
-tolerance_rel = 5./100
-
-print("error_rel    : " + str(error_rel))
-print("tolerance_rel: " + str(tolerance_rel))
+reflectivity = energy_end / energy_start
+reflectivity_max = 1e-06
 
-assert( error_rel < tolerance_rel )
+print("reflectivity     = " + str(reflectivity))
+print("reflectivity_max = " + str(reflectivity_max))
 
-# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E)
-Linf_norm = np.amax( np.abs( rho/scc.epsilon_0 - divE ) ) / np.amax( np.abs( rho/scc.epsilon_0 ) )
-assert( Linf_norm < 2.e-2 )
+assert(reflectivity < reflectivity_max)
 
 test_name = filename[:-9] # Could also be os.path.split(os.getcwd())[1]
 checksumAPI.evaluate_checksum(test_name, filename)