Convergence testing

A numerical calculation done at only one resolution without studying how the solution behaves as the resolution is increased is meaningless – Miguel Alcubierre, Section 9.11 of "Introduction to 3+1 Numerical Relativity".

Convergence tests are a key step of any numerical simulation. They help validate the solution and provide an estimate of the numerical error. In the context of GRTresna, these tests are crucial not only for ensuring the existence of a solution but also to confirm that the solution satisfies the Hamiltonian and momentum constraints. It is important to remember that just seeing a "small" constraint violation is not sufficient to verify that the solver is working correctly - it needs to decrease at the right rate to be sure that the error comes from the finite derivative approximations and not other sources (e.g. code bugs).

When the exact solution is not known (the solution in the continuum limit), it is essential to perform simulations at three different resolutions. This approach, outlined in Eqns. 9.11.5 and 9.11.6 of Alcubierre's book, ensures that the solution converges correctly as the resolution improves. In numerical relativity, the Hamiltonian and momentum constraints should be exactly zero in the continuum limit. This allows convergence tests to be performed using only two resolutions: low resolution (LR) and high resolution (HR). For example, consider a spatial lineout of the Hamiltonian constraint at two resolutions $\Delta_\mathrm{LR}$ and $\Delta_\mathrm{HR}$, with the values $H_\mathrm{LR}$ and $H_\mathrm{HR}$ respectively. Then, assuming a finite-difference scheme of order $n$, the errors at two different resolutions should satisfy:

$\frac{H_\mathrm{LR}}{H_\mathrm{HR}} = \Big(\frac{\Delta_{\mathrm{LR}}}{\Delta_{\mathrm{HR}}}\Big)^n$

In the output file Ham_and_Mom_error.txt you can find the norm of the relative error of the Hamiltonian and momentum constraints for each iteration of the solver. Although this is indicative of the reduction in constraint errors as the solver progresses, it is highly important that you also check the convergence of the final output as above. That is, you should confirm that the solution shows the expected 2nd order convergence as you increase the resolution.

Convergence testing the ScalarFieldCosmo example

The easiest way to increase the resolution of a simulation is by doubling the number of cells N on the coarsest level. GRTresna is relatively cheap to run, so the code should produce relatively quickly a InitialDataFinal.3d.hdf5 file for each simulation. As explained in Using the outputs in GRChombo, we can start the NR evolution from the initial data by changing the path of the restart_file parameter (make sure the grid parameters are the same in both GRChombo and GRTresna!). For this example, we have prepared consistent parameter files in the convergence_test folder within GRTresna and in the GRTresna_restart_setup within the GRChombo evolution example.

As soon as the evolution restarts, GRChombo will write a constraints_lineout.dat file in the specified data folder , This file contains the values of the Hamiltonian and momentum constraints along the "x" direction of the grid. Repeat this process for both the low and high resolution cases. We have prepared a python script plot_restart.py that plots and compares the values of the constraints along these lineouts. You should get a plot like this one:

The red lines correspond to the Hamiltonian and and momentum constraints of the low-resolution simulation. The dashed black line is inferred from $H_\mathrm{HR} = (\Delta_{\mathrm{HR}}/\Delta_{\mathrm{LR}})^n H_\mathrm{LR}$, where $n=2$. The green lines are the obtained Hamiltonian and momentum constraints for the high-resolution simulation. Remember GRTresna uses 2nd order stencils and thus these lines must match!

Convergence testing the ScalarFieldBH example

One of the main differences between the ScalarFieldCosmo and ScalarFieldBH examples is the use of Adaptive Mesh Refinement (AMR), which makes convergence testing more difficult. Ideally, one would like the difference in resolution between the LR and HR simulations to be the same at every point on the grid, but this can be tricky with AMR. Normally the best you can do is toggle the tagging criteria such that each refinement level covers roughly the same area, as you do when convergence testing during the evolution (see convergence testing for the evolution code). For this case we have also prepared parameter files in both GRTresna and GRChombo that ensure the AMR grids are aligned. When repeating the same process as above, you should get the following convergence plot

Note that the spikes correspond to the AMR boundaries across the grid, which in this case correspond to having max_level = 3 refinement levels.

If you are trying to reproduce these tests using a modified version of GRTresna and the solution does not converge at the right rate (i.e. you find the lines do not match), you should increase the resolution, or check for bugs in the code.