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. The specific checks performed and the associated return codes include:

return code 1:: Whether the simulation ran for $0.1$ seconds without the projectile leaving the barrel, indicating that the projectile likely won’t ever leave the barrel.
return code 2:: Whether a control volume has negative mass.
return code 3:: Whether a control volume has negative temperature.
return code 4:: Whether the total system mass deviates from the starting mass by more than $0.001\%$ .
return code 5:: Whether the total system energy deviates from the starting energy by more than $0.01\%$ .
return code 6:: If automatic differentiation is used, whether any component of the gradient of the total system mass deviates by more than $10.0\cdot 10^{-9}$ (regardless of the units).
return code 7:: If automatic differentiation is used, whether any component of the gradient of the total system energy deviates by more than $10.0\cdot 10^{-6}$ (regardless of the units).
return code 8:: Whether the ideal gas equation of state has become inaccurate due to the pressure increasing above the critical pressure.
return code 9:: Whether the coordinates of any plunger become desynchronized between control volumes on either side of the mentioned plunger by more than $1.0\cdot 10^{-12}$ m.
return code 10:: Whether the coordinate of a projectile/plunger goes negative, indicating that the control volume has negative volume and that the projectile/plunger has not only impacted the end but gone through the end.
return code 11:: Whether the maximum number of iterations ( $50$ ) was reached when interpolating to a particular projectile/plunger location. This indicates that BlasterSim was unable to properly handle an event like projectile exit or a plunger impacting the end of the plunger tube.

The thresholds used for some of these checks are somewhat arbitrary and are based on what conservatively seems to never cause problems in practice.

For a discussion of return codes -1 and 0, see § 1.4 on the CSV output variable rc. These return codes indicate no errors were detected.

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 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 [11, 12]. 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.4 Verification test suite

The test suite mentioned at the start of § 3 has 366 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. At the moment, the tests listed here are incomplete.

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 [6, 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 [7, § 2.3].

Comparison with plunger impact exact solution

Figure 3.2: Plunger impact test case.

In figure 3.2, a constant velocity ( $\dot{x}$ ) 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\,A\,\dot{x}.

(3.9)

Note that $\dot{x}$ is negative here, so $\dot{m}$ is positive. Because $\dot{m}$ , $A$ , and $\dot{x}$ , are constant, $\rho$ is necessarily also constant.

Equation 3.9 will necessarily be the mass flow rate as $x$ approaches 0 (as the plunger approaches impact) 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.

The complete solution is:

$\displaystyle m(t)$	$\displaystyle=\rho\,A(x_{0}+\dot{x}\,t),$	(3.10)
$\displaystyle E(t)$	$\displaystyle=\rho\,Au(x_{0}+\dot{x}\,t),$	(3.11)
$\displaystyle x(t)$	$\displaystyle=x_{0}+\dot{x}\,t,$	(3.12)
$\displaystyle\dot{x}(t)$	$\displaystyle=\dot{x},$	(3.13)
$\displaystyle u(t)$	$\displaystyle=u,$	(3.14)
$\displaystyle\rho(t)$	$\displaystyle=\rho,$	(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_{0}}{\dot{x}}.

(3.16)

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.4 can not be applied as the numerical error is negligible. After the plunger has moved for $0.1$ s, 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.