Verify numerical stability

[spco]

Related to #47 .

As a step towards verifying the output to a suitable standard, I think we need a suite of tests with known solutions. Currently our tests only really cover 'did I break anything from the last commit?', but in reality we also want to have confidence that our solutions are right. Typically this is done with the 'method of manufactured solutions', and I think it would be good to have a few tests touching on that. These typically can be very uncomplicated and short to run, so they shouldn't add considerably to the tests running time. I'm thinking of a few classes of tests:

1. Constant solution. Create mechanism files that, coupled with the initial conditions, give a constant solution. e.g.
   A + B -> C, rate k
   A + C -> B, rate k
   B + C -> A, rate k.
   A(0) = B(0) = C(0) = 1/3. I think this shouldn't matter whether we fix the zenith angles etc, as the mechanism is decoupled from the temperature and lighting. I'd probably actually have 2 versions - one with light varying and one without. Both should be identical.

2. Steady-state solution. Species concentrations vary from the start, but tend toward a steady-state solution (all concentrations hit a constant value and stay there). This can be used to verify that (a) the solution converges to a steady state at all, and (b) that steady state is the correct one according to the known solution.

3. 'Manufactured solution'. We pick the RHS of the differential equations so that the solution is known. Check that the calculated solution sticks to the true solution. These can be made arbitrarily complicated, which is useful for testing the limits of what the model can do.

[rs028]

Agree. Besides the current suite of tests has a lot of redundancy. You don't need a full MCM model to test a few switches.

[rs028]

Marked as high priority because the new tests (see PR #268) seem to indicate some minor calculation errors. Possibily related to the numerical issues we have been worried about.

[spco]

Just putting this here as it seems the most appropriate place: when looking at numerical stability, difference in the versions of BLAS/LAPACK should be considered as a possible source of issues. More rigorous unit testing of our own routines will help narrow this all down.

[spco]

#340 adds tests for static, first-order and second-order reactions. The results are good, but the first-order test shows that there are instabilities at the smallest concentrations. This is likely because CVODE is using the same (not the latest) concentration values throughout its iteration before it stops for the next timestep. This can result in overshoots. Perhaps it can be refactored a little to use the values from the latest sub-iteration rather than the beginning of the latest outer iteration.

[rs028]

The accuracy of the solver (i.e., the default values in solver.parameters) needs to be investigated. See notes in (#384 (comment)) for more information.

[rs028]

Additional notes about accuracy and stability of the solver:

• need to rationalize and expand the test suite - baseline numerical tests are needed (see Extend tests to cover more .fac file examples #47, Un-useful tests #326)
• need to upgrade CVODE (see Upgrade CVODE #257)
• numerical differences between Linux and macOS cause several tests to fail - currently the tests are set to high tolerance to avoid this problem and openlibm is used
• currently LAPACK/BLAS have been disabled (PR Disable LAPACK support in CVODE #460). Can LAPACK support in CVODE help with numerical differences/instabilities?