# Analytical testing/quality control for scientific software in professional setting

I am charged with maintaining a buildserver on Teamcity which is meant to test our in house FE software. Currently our test suite consists of a list of benchmarks which run every time a commit is made using svn as version control software. The benchmarks are fully meshed boundary value problems intended to test various features of the software. The features are tested by running the head revision of the kernel on each benchmark and then comparing an output file with a reference file. The reference file is a result for a previous run of the test which was deemed "correct" by an expert/academic. The tests are thus considered to be "broken" when the result of this output file is different from that generated by the initial run.

I would like to improve upon this system by having our test suite run comparisons based on physical/analytical solutions for each benchmark rather than an arbitrary output file. However, in this case, how does one judge when a test is "broken"? The analytical solution can and will never be exactly the same as the numerical one, so comparing two files is out of the question. Does there exist a general method/quality control process perhaps based on statistics which can ensure that our code remains scientifically viable? How do major software companies like ANSYS or Bentley ensure that their scientific software does not deviate from exact solutions?

Please note I am in a commercial software development office setting and am looking for things that can be implemented in a straightforward manner within the Teamcity/tortoisesvn/Visual Studio environment.

Has anyone had such experience before?

• Whether the baseline results are generated from a previous run or from analytical results, tests are almost never designed so that two floating point numbers must be exactly equal for the test to pass. The tests should be designed so that the new results compare with the baseline values within some tolerance; the tolerance is based on the judgement of the test designer. Most test frameworks, e.g. Google Test (github.com/google/googletest) have tolerance-based comparisons for floating point numbers. – Bill Greene Jul 18 '18 at 16:24
• For complex solid mechanics (e.g.., multiphase soil models with shocks) there are no analytical solutions and manufactured solutions are very difficult (if not impossible) to find. In those cases regression testing becomes the only option. I usually unit test individual parts of my code with googletests during compile time and do the regression tests before and after check in. – Biswajit Banerjee Jul 18 '18 at 21:13

There's a Sandia technical report discussing the Method of Manufactured Solutions (a generalisation of the concept of exact solutions to PDEs) which goes to some lengths to discuss the concept of acceptance criteria on verification tests for numerical codes. They propose 4 classes:

1. Expert Judgement,
2. Percentage Error,
3. Consistency,
4. Order of accuracy,

Here it appears you're currently using a version of expert judgment, tied to a regression test on the output data. Bill Greene's comment encompasses the idea of percent error, with an acceptable tolerance in the absolute difference between the numerically calculated value and the analytic one. A "good" choice of value of this tolerance is to some degree still a judgement call, and is a function of the discretisation, the implementation and the computer architecture being used.

Just testing values for a single discretisation of a given system is still less rigorous than convergence testing, in which runs are performed at multiple resolutions, and it is confirmed that the reduction in the error in the system scales with grid length in the manner theory suggests it should (of course this comparison of real numbers also requires a tolerance). The biggest annoyance with using this method under CI is that it can need moderately high resolution to sit in the region dominated by the leading order error term (rather than higher order terms above and machine precision below). This makes run times longer than comfortable at times when multiple commits are occuring.

I was once charged with solving this problem, in another setting. Here's what I did:

• I found conservation laws. These are pure gold, as is well known to physicists. Conservation of energy is present in many FE settings, but conservation laws abound in almost all physics models. If you can find a conservation law and compute the associated conserved quantity at (say) each time step, then you have vastly higher confidence in the associated solution.

• I manufactured solutions, as discussed in the other answer, and compared the numerical result to in in $L_{1}$, $L_2$ and $L_{\infty}$ norms (I like googletest for this purpose, but there are many good frameworks for the task). Of course these are never zero, but you can set some reasonable bound and when the tests break (and they will) you will have three good measures of whether that breakage is intolerable or not.

• For linear models, you have superposition and scale invariance. So you don't even need an analytic solution, you can just rescale your rhs and test that the new solution matches your scaled solution. You'd think this isn't a powerful test, but guess what . . . (Some will argue that this is the same as my first point.)

• I worked to make sure the code compiled on clang, not just MSVC. This allowed me to use AddressSanitizer, UndefinedBehaviorSanitizer, and ThreadSanitizer. These caught many memory leaks and race conditions, and provided output that allowed me to fix the associated problems without undue effort.

• I made the case that the output files which were being compared were forcing scientists to write parsers, which we know is exceptionally difficult. Once management became cognizant that most of the unit test failures were parser errors rather than code errors, my mandate and scope was expanded, and I had more authority to overhaul the code.