Hi,

I have solved the mystery of why the McLachlan test suites were failing on Kraken. The root cause of the problem was, again, the fact that the Exact thorn which sets up the gauge wave initial data calculates the extrinsic curvature imprecisely. The way the Exact thorn works is that the 4-metric is given in the code and is then finite differenced with a 2nd order accurate centered difference to compute the extrinsic curvature. It does not use the analytic derivatives. In order to make this accurate, a small step size, eps = 1e-6, is chosen for the finite differencing.

Unfortunately, this small value means that g(t+eps) and g(t-eps) are very similar to each other, and when they are subtracted from each other, the least significant bits in g become very significant in K. As Peter pointed out, this can be many orders of magnitude. So an error of 1e-15 caused by roundoff differences due to compilers, math libraries, etc in g can lead to an error of 1e-11 in K. This is exactly what was happening with the McLachlan test suite. kxx on the initial data, which has a maximum value of ~ 0.2, had differences between the test suite result and Kraken of O(1e-11). This led to a similar difference in the BSSN RHS computation, and from iteration 0 to iteration 1, the difference in gxx jumped 4 orders of magnitude.

To avoid this problem, Erik wrote a patch to implement fourth order accurate finite differencing in Exact. The idea is that you can get the same truncation error in the computation of K by using a higher order scheme with a larger step size. Since you are subtracting numbers which are not so similar to each other, the cancellation above does not occur, and roundoff errors are not amplified.

I have tested the patch with the gauge wave initial data by computing kxx with order = 4 and eps = 1e-1, 1e-2, 1e-3 and 1e-4. The convergence with eps is 4th order, and a Richardson extrapolation reveals that the 1e-4 result has an error which is O(1e-15); i.e. using order = 4 and eps = 1e-4 for this initial data yields a result for kxx which is correct up to roundoff error.

I attach a plot of the difference in gxx beween using order = 2, eps = 1e-6 and using order = 4, eps = 1e-4 between my laptop and Kraken. We see that after the initial data, the difference jumps from 1e-15 to 1e-12 with order = 2, eps = 1e-6, whereas with order = 4, eps = 1e-4, the difference remains small.

Since we are close to a release, Erik doesn't want to commit his patch at this point. We understand why the McLachlan tests fail on kraken, and I have verified that with the patch, they can be made to pass. After the release, we will commit Erik's patch and modify the test suites to use it.

We are not sure how to proceed in general with Exact. It is not clear that a single value of eps is suitable for all of the models in Exact. The current patch could be improved to perform a Richardson extrapolation to estimate the error, but it is then not clear how to reliably calculate the relative error that we need to minimise.