Merge pull request #181 from lin-yang-ly/add_rmspe

Add rmspe on the V&V cases missing it

doc/content/verification_and_validation/val-2b.md

@@ -41,7 +41,7 @@ The model applies 13.3 kPa of D$_2$ for 50 hours and 15 seconds followed by evac
 ## Results
 
 
-[val-2b_comparison] shows the comparison of the TMAP8 calculation and the experimental data. There is good agreement between the TMAP predictions and the data.
+[val-2b_comparison] shows the comparison of the TMAP8 calculation and the experimental data. There is reasonable agreement between the TMAP predictions and the experimental data with a root mean square percentage error of RMSPE = 19.41 %. Note that the agreement could be improved by adjusting the model parameters. It is also possible to perform this optimization with [MOOSE's stochastic tools module](https://mooseframework.inl.gov/modules/stochastic_tools/index.html).
 
 !media comparison_val-2b.py
        image_name=val-2b_comparison.png

doc/content/verification_and_validation/ver-1a.md

@@ -62,7 +62,7 @@ The expression in Eq. (1) of [!cite](longhurst1992verification) writes $-\alpha_
 
 ### Results
 
-[ver-1a_comparison_analytical_TMAP4_release_fraction] shows the comparison of the TMAP8 calculation and the analytical solution provided in [!cite](longhurst1992verification). There is agreement between the two plots.
+[ver-1a_comparison_analytical_TMAP4_release_fraction] shows the comparison of the TMAP8 calculation and the analytical solution provided in [!cite](longhurst1992verification). There is agreement between the two plots with a root mean square percentage error of RMSPE = 0.19 % for $t \geq 1$ s.
 
 !media comparison_ver-1a.py
        image_name=ver-1a_comparison_analytical_TMAP4_release_fraction.png
@@ -110,7 +110,7 @@ which leads to
 
 ### Results
 
-[ver-1a_comparison_analytical_TMAP7_release_fraction] shows the comparison of the TMAP8 calculation and the analytical solution for release fraction provided in [!cite](ambrosek2008verification). There is agreement between the two plots.
+[ver-1a_comparison_analytical_TMAP7_release_fraction] shows the comparison of the TMAP8 calculation and the analytical solution for release fraction provided in [!cite](ambrosek2008verification). There is agreement between the two plots with a root mean square percentage error of RMSPE = 0.07 % for $t \geq 1$ s.
 
 !media comparison_ver-1a.py
        image_name=ver-1a_comparison_analytical_TMAP7_release_fraction.png
@@ -134,7 +134,7 @@ Again, be aware of the typos in [!cite](ambrosek2008verification) and the differ
 
 ### Results
 
-[ver-1a_comparison_analytical_TMAP7_flux] shows the comparison of the TMAP8 calculation and the analytical solution for flux at the outer surface of the SiC slab in [!cite](ambrosek2008verification). There is agreement between the two plots.
+[ver-1a_comparison_analytical_TMAP7_flux] shows the comparison of the TMAP8 calculation and the analytical solution for flux at the outer surface of the SiC slab in [!cite](ambrosek2008verification). There is agreement between the two plots with a root mean square percentage error of RMSPE = 0.26 % for $t \geq 1$ s.
 
 !media comparison_ver-1a.py
        image_name=ver-1a_comparison_analytical_TMAP7_flux.png

doc/content/verification_and_validation/ver-1b.md

@@ -48,7 +48,9 @@ The flux as given by [eq:flux] is compared with values calculated by TMAP8.
 The diffusivity, D, and the initial concentration, C$_o$, were both
 taken as unity, and the distance, x, was taken as 0.5 in this comparison.
 TMAP8 initially under predicts but the results match well subsequently. Comparison
-results are shown in [ver-1b_comparison_flux]
+results are shown in [ver-1b_comparison_flux] with a root mean square percentage
+error of RMSPE = 6.03 %. The error is calculated for $t \geq 10$ s due to infinite
+value at small $t$.
 
 !media comparison_ver-1b.py
        image_name=ver-1b_comparison_flux.png

doc/content/verification_and_validation/ver-1d.md

@@ -86,7 +86,7 @@ Using TMAP8 we examine these two different regimes, one where diffusion is the r
 
 ## Diffusion-limited
 
-For the effective diffusivity limit, we selected $\epsilon/k = 100 K$ to give $\zeta = 91.47 c/\rho$. The comparison results are presented in [ver-1d_comparison_diffusion].
+For the effective diffusivity limit, we selected $\epsilon/k = 100 K$ to give $\zeta = 91.47 c/\rho$. The comparison results are presented in [ver-1d_comparison_diffusion] with a root mean square percentage error of RMSPE = 0.96 % for $t \geq 0.4$ s..
 
 !media comparison_ver-1d.py
        image_name=ver-1d_comparison_diffusion.png

doc/content/verification_and_validation/ver-1dc.md

@@ -82,7 +82,7 @@ The $\zeta_i$ values of the three traps from [!cite](ambrosek2008verification) h
 
 ## Results and comparison against analytical solution
 
-The analytical solution for the permeation transient is compared with TMAP8 results in [ver-1dc_comparison_diffusion]. The graphs for the theoretical flux and the calculated flux are in good agreement, with root mean square percentage errors (RMSPE) of RMSPE = 0.41 % when time $t \geq 3$ s. The breakthrough time calculated from [eqn:tau_be] in analytical solution is 4.04 s, and the breakthrough time from TMAP8 is 4.12 s.
+The analytical solution for the permeation transient is compared with TMAP8 results in [ver-1dc_comparison_diffusion]. The graphs for the theoretical flux and the calculated flux are in good agreement, with root mean square percentage errors (RMSPE) of RMSPE = 0.41 % for $t \geq 3$ s. The breakthrough time calculated from [eqn:tau_be] in analytical solution is 4.04 s, and the breakthrough time from TMAP8 is 4.12 s.
 
 !media comparison_ver-1dc.py
        image_name=ver-1dc_comparison_diffusion.png

doc/content/verification_and_validation/ver-1fa.md

@@ -23,7 +23,7 @@ The slab is assumed to have a density of 1 kg/m$^3$ and a specific heat capabity
 
 Comparison of the temperature computed through TMAP8 and calculated analytically is shown in
 [ver-1fa_comparison_temperature]. The TMAP8 code predictions match very well with
-the analytical solution.
+the analytical solution with a root mean square percentage error of RMSPE = 0.05 %.
 
 !media comparison_ver-1fa.py
        image_name=ver-1fa_comparison_temperature.png

doc/content/verification_and_validation/ver-1fb.md

@@ -37,7 +37,7 @@ $k$ is the thermal conductivity, $\rho$ is the density and $C_p$ is the specific
 #
 
 
-Comparison of the temperature distribution in the slab, computed through TMAP8 and calculated analytically, is shown in [ver-1fb_comparison_temperature]. The TMAP8 code predictions match very well with the analytical solution.
+Comparison of the temperature distribution in the slab, computed through TMAP8 and calculated analytically, is shown in [ver-1fb_comparison_temperature]. The TMAP8 code predictions match very well with the analytical solution with the root mean square percentage errors of RMSPE = 0.09 % at $t = 0.1$ s, RMSPE = 0.03 % at $t = 0.5$ s, RMSPE = 0.02 % at $t = 1$ s, and RMSPE = 0.00 % at $t = 5$ s.
 
 !media comparison_ver-1fb.py
        image_name=ver-1fb_comparison_temperature.png

doc/content/verification_and_validation/ver-1g.md

@@ -42,7 +42,7 @@ C_{i_0} = 10^{-18} \frac{P_{i_0} N_a}{R T},
 
 where $P_{i_0}$ is the initial pressure, $N_a$ is Avogardro's constant, $R$ is the gas constant (from  [PhysicalConstants](source/utils/PhysicalConstants.md)), and $T$ = 298.15 K (25$\deg$ C) is the temperature. The factor $10^{-18}$ converts the concentration from atoms/m$^3$ to atoms/$\mu$m$^3$.
 
-A comparison of the concentration of AB as a function of time is plotted in [ver-1g_comparison_equal_conc] for case (a), and [ver-1g_comparison_diff_conc] for the cases (b) and (c), respectively. The TMAP8 calculations are found to be in good agreement with the analytical solution.
+A comparison of the concentration of AB as a function of time is plotted in [ver-1g_comparison_equal_conc] for case (a), and [ver-1g_comparison_diff_conc] for the cases (b) and (c), respectively. The TMAP8 calculations are found to be in good agreement with the analytical solution with the root mean square percentage errors of (a) RMSPE = 0.27 %, (b) RMSPE = 0.22 %, and (c) RMSPE = 0.24 %.
 
 !media comparison_ver-1g.py
        image_name=ver-1g_comparison_equal_conc.png

test/tests/val-2b/comparison_val-2b.py

@@ -9,6 +9,24 @@
 script_folder = os.path.dirname(__file__)
 os.chdir(script_folder)
 
+def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_output):
+    """Get new numerical solution based on the experimental input data points
+
+    Args:
+        experiment_input (float, ndarray): experimental input data points
+        tmap_input (float, ndarray): numerical input data points
+        tmap_output (float, ndarray): numerical output data points
+
+    Returns:
+        float, ndarray: updated tmap_output based on the data points in experiment_input
+    """
+    new_tmap_output = np.zeros(len(experiment_input))
+    for i in range(len(experiment_input)):
+        left_limit = np.argwhere((np.diff(tmap_input < experiment_input[i])))[0][0]
+        right_limit = left_limit + 1
+        new_tmap_output[i] = (experiment_input[i] - tmap_input[left_limit]) / (tmap_input[right_limit] - tmap_input[left_limit]) * (tmap_output[right_limit] - tmap_output[left_limit]) + tmap_output[left_limit]
+    return new_tmap_output
+
 fig = plt.figure(figsize=[6.5, 5.5])
 gs = gridspec.GridSpec(1, 1)
 ax = fig.add_subplot(gs[0])
@@ -50,7 +68,10 @@
 ax.legend(loc="best")
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+tmap_flux_for_rmspe = numerical_solution_on_experiment_input(expt_temp, tmap_temp, tmap_flux)
+RMSE = np.sqrt(np.mean((tmap_flux_for_rmspe-expt_flux)**2) )
+RMSPE = RMSE*100/np.mean(expt_flux)
+ax.text(870,3e16, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('val-2b_comparison.png', bbox_inches='tight')
 plt.close(fig)

test/tests/ver-1a/comparison_ver-1a.py

@@ -151,7 +151,21 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
     flux = 2 * S * P_0 * L * D * summation
     return flux
 
-time_analytical = np.linspace(0, 140, 1000)
+# Extract data from 'gold' TMAP8 run
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder = "../../../../test/tests/ver-1a/gold/ver-1a_csv.csv"
+else:                                  # if in test folder
+    csv_folder = "./gold/ver-1a_csv.csv"
+tmap8_prediction = pd.read_csv(csv_folder)
+tmap8_time = tmap8_prediction['time']
+tmap8_release_fraction_right = tmap8_prediction['released_fraction_right']
+tmap8_release_fraction_left = tmap8_prediction['released_fraction_left']
+tmap8_flux_right = tmap8_prediction['flux_surface_right'] # at/microns^2/s
+tmap8_flux_right = tmap8_flux_right*1e6*1e6 # at/m^2/s
+idx = np.where(tmap8_time >= 1.0)[0][0]
+
+# time_analytical = np.linspace(0, 140, 1000)
+time_analytical = np.array(tmap8_time)
 T=2373
 P_0=1e6
 D=1.58e-4*np.exp(-308000.0/(R*T))
@@ -190,18 +204,6 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
     l=l,
 )
 
-# Extract data from 'gold' TMAP8 run
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
-    csv_folder = "../../../../test/tests/ver-1a/gold/ver-1a_csv.csv"
-else:                                  # if in test folder
-    csv_folder = "./gold/ver-1a_csv.csv"
-tmap8_prediction = pd.read_csv(csv_folder)
-tmap8_time = tmap8_prediction['time']
-tmap8_release_fraction_right = tmap8_prediction['released_fraction_right']
-tmap8_release_fraction_left = tmap8_prediction['released_fraction_left']
-tmap8_flux_right = tmap8_prediction['flux_surface_right'] # at/microns^2/s
-tmap8_flux_right = tmap8_flux_right*1e6*1e6 # at/m^2/s
-
 # Plot figure for verification of release fraction as determined in TMAP4 (SiC outer layer)
 fig = plt.figure(figsize=[6.5,5.5])
 gs = gridspec.GridSpec(1,1)
@@ -215,6 +217,9 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
 ax.set_xlim(right=140)
 ax.set_ylim(bottom=0)
 plt.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((tmap8_release_fraction_right-analytical_release_fraction_TMAP4)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_release_fraction_TMAP4[idx:])
+ax.text(60,0.6, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('ver-1a_comparison_analytical_TMAP4_release_fraction.png', bbox_inches='tight');
 plt.close(fig)
@@ -232,6 +237,9 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
 ax.set_xlim(right=140)
 ax.set_ylim(bottom=0)
 plt.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((tmap8_release_fraction_left-analytical_release_fraction_TMAP7)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_release_fraction_TMAP7[idx:])
+ax.text(40,0.6, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('ver-1a_comparison_analytical_TMAP7_release_fraction.png', bbox_inches='tight');
 plt.close(fig)
@@ -249,6 +257,9 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
 ax.set_xlim(right=140)
 ax.set_ylim(bottom=0)
 plt.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((tmap8_flux_right-analytical_flux_TMAP7)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_flux_TMAP7 [idx:])
+ax.text(60,0.6e19, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('ver-1a_comparison_analytical_TMAP7_flux.png', bbox_inches='tight');
 plt.close(fig)

test/tests/ver-1b/comparison_ver-1b.py

@@ -22,6 +22,7 @@
 tmap_sol = pd.read_csv(csv_folder)
 tmap_time = tmap_sol['time']
 tmap_conc = tmap_sol['conc_point1']
+idx = np.where(tmap_time >= 10.0)[0][0]
 
 analytical_time = tmap_time
 x = 0.2
@@ -39,9 +40,11 @@
 ax.set_xlim(right=45)
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+RMSE = np.sqrt(np.mean((tmap_conc-analytical_conc)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_conc[idx:])
+ax.text(5,0.9, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
-plt.savefig('ver-1b_comparison_time.png', bbox_inches='tight');
+plt.savefig('ver-1b_comparison_time.png', bbox_inches='tight')
 plt.close(fig)
 
 
@@ -75,9 +78,11 @@
 ax.set_xlim(right=50)
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+RMSE = np.sqrt(np.mean((tmap_conc-analytical_conc)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_conc[idx:])
+ax.text(10,0.4, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
-plt.savefig('ver-1b_comparison_dist.png', bbox_inches='tight');
+plt.savefig('ver-1b_comparison_dist.png', bbox_inches='tight')
 plt.close(fig)
 #================== Comparison of flux as a function of time ===================
 
@@ -110,7 +115,9 @@
 ax.set_xlim(right=45)
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+RMSE = np.sqrt(np.mean((tmap_flux-analytical_flux)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_flux[idx:])
+ax.text(10,0.25, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
-plt.savefig('ver-1b_comparison_flux.png', bbox_inches='tight');
+plt.savefig('ver-1b_comparison_flux.png', bbox_inches='tight')
 plt.close(fig)

test/tests/ver-1d/comparison_ver-1d.py

@@ -55,9 +55,17 @@ def summation_term(num_terms, time):
         sum += (-1)**m * np.exp(-1 * m**2 * time / (2*tau_be))
     return sum
 
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder = "../../../../test/tests/ver-1d/gold/ver-1d-diffusion_out.csv"
+else:                                  # if in test folder
+    csv_folder = "./gold/ver-1d-diffusion_out.csv"
+tmap_sol = pd.read_csv(csv_folder)
+tmap_time = tmap_sol['time']
+tmap_perm = tmap_sol['scaled_outflux']
+idx = np.where(tmap_time >= 0.4)[0][0]
 
 c_o = c
-analytical_time = np.linspace(0.0001, 3.0, 500)
+analytical_time = tmap_time
 Jp = N_o * (c_o * D / l) * \
     (1 + 2 * summation_term(num_summation_terms, analytical_time))
 
@@ -69,13 +77,6 @@ def summation_term(num_terms, time):
 ax.plot(analytical_time, analytical_permeation,
         label=r"Analytical", c='k', linestyle='--')
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
-    csv_folder = "../../../../test/tests/ver-1d/gold/ver-1d-diffusion_out.csv"
-else:                                  # if in test folder
-    csv_folder = "./gold/ver-1d-diffusion_out.csv"
-tmap_sol = pd.read_csv(csv_folder)
-tmap_time = tmap_sol['time']
-tmap_perm = tmap_sol['scaled_outflux']
 ax.plot(tmap_time, tmap_perm, label=r"TMAP8", c='tab:gray')
 
 ax.set_xlabel(u'Time(s)')
@@ -84,7 +85,9 @@ def summation_term(num_terms, time):
 ax.set_xlim(left=0)
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+RMSE = np.sqrt(np.mean((tmap_perm-analytical_permeation)[idx:]**2) )
+RMSPE = RMSE*100/np.mean(analytical_permeation[idx:])
+ax.text(1.0,0.5e18, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('ver-1d_comparison_diffusion.png', bbox_inches='tight')
 plt.close(fig)

test/tests/ver-1fa/comparison_ver-1fa.py

@@ -9,25 +9,25 @@
 script_folder = os.path.dirname(__file__)
 os.chdir(script_folder)
 
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder = "../../../../test/tests/ver-1fa/gold/ver-1fa_csv_line_0011.csv"
+else:                                  # if in test folder
+    csv_folder = "./gold/ver-1fa_csv_line_0011.csv"
+tmap_sol = pd.read_csv(csv_folder)
+tmap_x = tmap_sol['id']
+tmap_temp = tmap_sol['temp']
+
 fig = plt.figure(figsize=[6.5,5.5])
 gs = gridspec.GridSpec(1,1)
 ax = fig.add_subplot(gs[0])
 
-analytical_x = np.linspace(0.0, 1.6, 40)
+analytical_x = tmap_x
 Ts = 300
 k = 10
 L = 1.6
 Q = 10000
 analytical_temp = Ts + Q*L**2 * (1- analytical_x**2/L**2) / (2*k)
 ax.scatter(analytical_x,analytical_temp,label=r"Analytical",c='k', marker='^')
-
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
-    csv_folder = "../../../../test/tests/ver-1fa/gold/ver-1fa_csv_line_0011.csv"
-else:                                  # if in test folder
-    csv_folder = "./gold/ver-1fa_csv_line_0011.csv"
-tmap_sol = pd.read_csv(csv_folder)
-tmap_x = tmap_sol['id']
-tmap_temp = tmap_sol['temp']
 ax.plot(tmap_x,tmap_temp,label=r"TMAP8",c='tab:gray')
 
 ax.set_xlabel(u'Distance along slab (m)')
@@ -36,7 +36,9 @@
 #ax.set_xlim(left=0)
 ax.set_ylim(bottom=0)
 plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3)
-
+RMSE = np.sqrt(np.mean((tmap_temp-analytical_temp)**2) )
+RMSPE = RMSE*100/np.mean(analytical_temp)
+ax.text(0.5,1000, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold')
 ax.minorticks_on()
 plt.savefig('ver-1fa_comparison_temperature.png', bbox_inches='tight');
 plt.close(fig)