3 Verification and validation 3 Verification and validation 3.2 Validation

3.1 Verification

Verification tests check that the simulation model has been implemented consistently. In other words, verification tests check that the math of the simulator is consistent with the intended governing equations. These tests do not test that the governing equations are appropriate in the first place, which is done through BlasterSim’s validation tests.

3.1.1 Run-time consistency checks

Every time step, BlasterSim performs internal consistency checks to alert the user if the simulation is becoming inaccurate in a detectable way. If an error is detected, BlasterSim terminates the simulation and prints out the return code. These checks are made in the check_sys subroutine of cva.f90. See § 1.5 for details.

3.1.2 Debugging run-time assertions

BlasterSim contains hundreds of run-time assertions which perform deeper internal consistency checks. These additional checks are disabled in the released versions of BlasterSim for speed but are enabled during all developer testing runs for debugging purposes.

3.1.3 Fuzz testing

BlasterSim is periodically tested with randomly generated input files, which is called fuzz testing. Fuzz testing has mostly improved BlasterSim’s robustness and had no impact on the correctness of the result. Before fuzz testing, BlasterSim would crash on certain inputs either due to mass or energy going negative in a control volume, the energy conservation tolerance being violated, or an invalid floating point operation, for example. These errors were mostly solved by adding stricter input validation to catch physically impossible inputs (that the fuzz tester itself can not know are physically impossible), better internal error detection and reporting, adding adaptive time stepping in response to mass and energy conservation errors reaching a certain threshold, and improving the numerical stability of floating point operations. These robustness improvements will be useful when design optimization is added as an optimizer is similar to a fuzz tester in that it will try many strange inputs, and BlasterSim not crashing is necessary for the process to continue.

3.1.4 Compile-time unit checking

Variables in BlasterSim are assigned units like meters or seconds, allowing the Fortran compiler to check for inconsistencies. This is achieved through a system called genunits which generates a custom Fortran module containing the required Fortran types [15, 16]. This system has proved useful to locate bugs as the compiler will often know precisely where the unit inconsistency is. A failing test will rarely provide such precision.

Some examples will illustrate how this works from the programmer’s perspective. Consider the following Fortran code:

⬇

program test_units_pass

use units

implicit none

type(si_length) :: x

type(si_time) :: t

type(si_velocity) :: v

x%v = 1.0

t%v = 1.0

v = x / t

print "(f3.1)", v

end program test_units_pass

The units module contains Fortran derived types and associated operations which allow for unit checking. When test_units_pass.f90 is compiled and run, the following output is produced:

⬇

1.0

However, consider what happens the velocity v is given the wrong unit:

⬇

program test_units_fail

use units

implicit none

type(si_length) :: x, v

type(si_time) :: t

x%v = 1.0

t%v = 1.0

v = x / t

print "(f3.1)", v

end program test_units_fail

test_units_fail.f90 will not compile, returning the following when compiled with gfortran :

⬇

test_units_fail.f90:12:4:

12 | v = x / t

| 1

Error: Cannot convert TYPE(si_velocity) to TYPE(si_length) at (1)

The compiler error message quickly tells the programmer that this code has a bug. In this case, the bug is in the variable definition and not in the highlighted line, however, often the compiler will identify precisely which term in the line has the bug.

The specific units used in cva.f90 and BlasterSim’s core simulation procedures otherwise are SI units. BlasterSim input and output may use other units, but the genunits types themselves represent SI units.

3.1.5 Verification test suite

The test suite mentioned at the start of § 3 has 390 the vast majority of which are verification tests checking that various functions and subroutines of BlasterSim produce the expected results. This section details some of the more important verification tests.

Comparison with single control volume exact solution

To test that the time integration in BlasterSim is working correctly, an exact solution for a simple case was constructed. This case has only two control volumes, one for the barrel, and one for the atmosphere. See figure 3.1 for an illustration of this test case.

Figure 3.1: Single control volume test case.

The barrel is initially filled with pressurized gas and the projectile has both dynamic and static friction pressure set to $p_{\text{f}}$ . The tube has a constant cross-sectional area of $A$ , so the volume of the barrel is $V=A\,x$ . The gas is ideal with a constant ratio of specific heats, $\gamma$ , and the process is adiabatic, so $p\,V^{\gamma}=\text{constant}$ in the barrel [8, p. 129, eq. 3.53]. So for any time $t$ , the pressure in the barrel is

p_{0}{(A\,x_{0})}^{\gamma}=p(t){(A\,x(t))}^{\gamma},

(3.1)

which can be rearranged and simplified to

p(t)=p_{0}{\left(\frac{x_{0}}{x(t)}\right)}^{\gamma},

(3.2)

where $P_{0}$ and $x_{0}$ are the initial pressure and projectile position, respectively.

The equation of motion of the projectile is

	$\displaystyle m_{\text{p}}\derivative{\dot{x}}{t}$	$\displaystyle=A\left(p(t)-p_{\text{atm}}-p_{\text{f}}\right),$		(3.3)
		$\displaystyle=A\left[p_{0}{\left(\frac{x_{0}}{x(t)}\right)}^{\gamma}-p_{\text{% atm}}-p_{\text{f}}\right].$		(3.4)

The trick $\displaystyle\derivative{\dot{x}}{t}\,\differential x=\dot{x}\,\differential% \dot{x}$ can be applied to solve this nonlinear differential equation for $\dot{x}$ as a function of $x$ rather than $t$ . Substituting in the trick returns

m_{\text{p}}\dot{x}\,\differential\dot{x}=A\left[p_{0}{\left(\frac{x_{0}}{x}% \right)}^{\gamma}-p_{\text{atm}}-p_{\text{f}}\right]\,\differential x.

(3.5)

Now integrate from an initial velocity ( $\dot{x}_{0}$ ):

\int_{\dot{x}_{0}}^{\dot{x}}m_{\text{p}}\hat{\dot{x}}\,\differential\hat{\dot{% x}}=\int_{x_{0}}^{x}A\left[p_{0}{\left(\frac{x_{0}}{\hat{x}}\right)}^{\gamma}-% p_{\text{atm}}-p_{\text{f}}\right]\,\differential\hat{x}.

(3.6)

Evaluating the integrals returns the projectile kinetic energy as a function of $x$ :

\frac{1}{2}m_{\text{p}}\dot{x}^{2}=\frac{1}{2}m_{\text{p}}\dot{x}_{0}^{2}+A% \left[\frac{p_{0}\left(x_{0}^{\gamma}\,x^{1-\gamma}-x_{0}\right)}{1-\gamma}-(p% _{\text{atm}}+p_{\text{f}})(x-x_{0})\right].

(3.7)

In terms of the projectile velocity, the result is

\dot{x}=\sqrt{\dot{x}_{0}^{2}+\frac{2A}{m_{\text{p}}}\left[\frac{p_{0}\left(x_% {0}^{\gamma}\,x^{1-\gamma}-x_{0}\right)}{1-\gamma}-(p_{\text{atm}}+p_{\text{f}% })(x-x_{0})\right]}.

(3.8)

This case was implemented in test_cva.f90 in the subroutine test_exact. For $\Delta t=27.3\cdot 10^{-6}$ s, the numerical error was measured at $0.484\cdot 10^{-9}$ m/s, which is negligible. The numerical order-of-accuracy was measured at $4.012$ , which is close to the theoretically expected order-of-accuracy of $4$ for the RK4 integrator. In other words, not only does BlasterSim closely match the exact solution, but BlasterSim’s numerical error decreases at the theoretically expected rate as the time step decreases. This is the most rigorous test that can be made with an exact solution. For more information on order-of-accuracy verification, see Roy [10, § 2.3].

Comparison with plunger impact exact solution

Figure 3.2: Plunger impact test case.

In figure 3.2, a constant velocity ( $\dot{x}_{1}$ ) plunger is approaching the end of the plunger tube. In BlasterSim, a constant velocity plunger can be created by making the plunger have infinite mass, so no force will be able to change the plunger’s velocity. The mass flow rate out of the plunger tube is set to a constant $\dot{m}$ defined so that it perfectly balances the plunger mass flow rate:

\dot{m}=-\rho_{1}\,A\,\dot{x}_{1}.

(3.9)

Note that $\dot{x}_{1}$ is negative here, so $\dot{m}$ is positive. Because $\dot{m}$ , $A$ , and $\dot{x}_{1}$ , are constant, $\rho_{1}$ is necessarily also constant. The internal energy $u_{1}$ and temperature $T_{1}$ also end up being constant due to the balance of the mass flow rate, but the derivation is omitted for brevity.

Equation 3.9 will necessarily be the mass flow rate as $x$ approaches 0 because there will be no volume to absorb the difference between the plunger mass flow rate and the mass flow rate out of the plunger tube, so this test case should approximate conditions as the plunger approaches impact. Note that for numerical reasons, the minimum $x$ position of the plunger will be $x_{\text{min}}$ , which is greater than zero, as discussed in § 2.1.9, so $x$ will not actually approach zero.

The solution before impact is:

$\displaystyle m(t)$	$\displaystyle=\rho_{1}\,A(x_{0}+\dot{x}_{1}\,t),$	(3.10)
$\displaystyle E(t)$	$\displaystyle=\rho_{1}\,A\,u(x_{0}+\dot{x}_{1}\,t),$	(3.11)
$\displaystyle x(t)$	$\displaystyle=x_{0}+\dot{x}_{1}\,t,$	(3.12)
$\displaystyle\dot{x}(t)$	$\displaystyle=\dot{x}_{1},$	(3.13)
$\displaystyle u(t)$	$\displaystyle=u_{1},$	(3.14)
$\displaystyle\rho(t)$	$\displaystyle=\rho_{1},$	(3.15)

where $x_{0}$ is the initial plunger position.

The plunger impacts the end of the plunger tube at time

t_{\text{impact}}=\frac{x_{\text{min}}-x_{0}}{\dot{x}_{1}}.

(3.16)

The solution after plunger impact might not be possible to express in closed form and is not of particular interest aside from the plunger velocity after impact, $\dot{x}_{2}=-e_{\text{p}}\dot{x}_{1}$ .

The solution in this case is simply a series of linear algebraic equations, and BlasterSim solves these equations extremely accurately even for large time steps. Consequently, the order-of-accuracy test used in § 3.1.5 can not be applied as the numerical error is negligible. After the plunger has moved for $0.10$ s (which is well before plunger impact), the error in $m$ is $3.109\cdot 10^{-15}$ kg, the error in $E$ is $2.298$ J (negligible compared against the expected energy of $1.441\cdot 10^{6}$ J), the error in $x$ is $0.0$ m, and the error in $\dot{x}$ is $0.0$ m/s. The error in $t_{\text{impact}}$ is $2.916\cdot 10^{-12}$ s. Immediately after impact, the error in $\dot{x}$ is $0.0$ m/s.