From dbd252d9b83ce3bbdb66021dc0ed50f5b94b7830 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Mon, 16 Nov 2020 15:33:15 -0800 Subject: 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 --- Examples/Tests/PML/analysis_pml_psatd.py | 42 ++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 18 deletions(-) (limited to 'Examples/Tests/PML/analysis_pml_psatd.py') 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) -- cgit v1.2.3