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. These checks are made in the check_sys subroutine of cva.f90. The specific checks performed include:

•

Whether the simulation ran for too long without the projectile leaving the barrel. (run subroutine return code 1)
•

Whether a control volume has negative mass. (run subroutine return code 2)
•

Whether a control volume has negative temperature. (run subroutine return code 3)
•

Whether the total system mass deviates from the starting mass by more than $0.001\%$ . (run subroutine return code 4)
•

Whether the total system energy deviates from the starting energy by more than $0.01\%$ . (run subroutine return code 5)
•

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). (run subroutine return code 6)
•

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). (run subroutine return code 7)
•

Whether the ideal gas equation of state has become inaccurate due to the pressure increasing above the critical pressure. (run subroutine return code 8)
•

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. (run subroutine return code 9)

The thresholds used for 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 [10, 11]. 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.

3.1.4 Comparison with 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: Test case which can be solved exactly.

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].