Solver methods

Basic principle of the solver

The multigrid elliptic solver we use aims to solve a set of coupled elliptic equations, each of the form

$$ {\bf a}(x) ~ {X} + {\bf b}(x) ~ \delta^{ij} \partial_i \partial_j {X} = {\bf S}(x)$$

Here ${\bf a}$ and ${\bf b}$ are coefficients that may be complicated functions of the matter input, and ${\bf S}$ is the source or "right hand side"/"rhs", which again may (and most likely will in GR) be a complicated function of the inputs.

Any equation given to the solver has to be linearised to cast it in the form above, if not already linear, since the coefficients and source cannot depend on the solver variables themselves. Once linearised, we start with a guess for the linearised variable, say $\psi_0$, and we solve for a correction $\delta \psi$ and update it over a series of non linear steps. This allows us to use the value of $\psi_0$ at the previous step in places where it appears in the coefficients and source. That is, the coefficients will depend on $\psi_0$ but not $\delta \psi$. If the solver converges, $\delta \psi$ should tend to zero over the non linear steps.

For example, in the Hybrid CTTK method, the equation for $\psi$ is linearised, but the equations for $V_i$ are already linear, so they are implemented as such. Whether the equation is linearised or not will affect how the corresponding variable is updated after each solver round, since for a linearised equation what is actually being solved is the change in the variable, e.g. $\delta \psi$. The variable updates are done in the function update_psi0 in Source/Core/Grids.cpp.

All the current methods require some non linear iterations, because although the individual equations are linear, they are coupled such that updates of the Hamiltonian constraint affect the momentum constraint at the next step and vice versa.

Please refer to the CTTK article for more detailed explanations of these points.

Current solver methods

The current version of the solver supports the use of two different methods for solving the constraint equations: CTTK and Hybrid CTTK. See this article for full details of each approach.

The basic difference between the two methods is that in the Hybrid approach the conformal factor $\psi$ varies in response to the $A_{ij}$ sources in the Hamiltonian constraint, whereas in the CTTK method it is simply set to a constant $\psi = 1$, and the variation in $K$ absorbs the $A_{ij}$ sources. CTTK Hybrid is usually the default for black hole spacetimes, whereas the CTTK method is mainly used in periodic/cosmological spacetimes. In principle, both can be used in either case, but often one provides better convergence in certain cases.

Setting the solver method

The desired method is chosen in the GNUMakefile, by leaving the line with the desired method uncommented and commenting out the other options. For example, for using the Hybrid CTTK method, this line should be used:

cxxcppflags = -DUSE_CTTKHybrid

Then, the command make realclean should be used on the terminal after which the whole solver should be recompiled so that the desired method is implemented in the compilation.

Creating a new solver method

To create a new solver method, one needs to create a class that provides methods for calculating the rhs, aCoef and bCoef values throughout the grid. These need to be defined for each of the variables to be solved for. Probably the easiest way to do this is to copy and amend one of the existing example methods (e.g., CTTK.hpp/CTTK.impl.hpp).

The methods themselves are coded in the Source/Method folder. Each method defines a class, with five corresponding functions. The functions and the variables defined through them are listed below:

• read_params : reads the appropriate parameters from the params.txt file

• solve_analytic : computes $K$, $\bar{A}{ij}$, and the ADM components of the EM tensor, $\rho$ $S_i$ and $S{ij}$.

• set_elliptic_terms : sets the RHS of the constraint equations and the aCoef (multiplying the constraint variable in the equation)

• initialise_method_vars: initialises the multigrid variables to 0 or values from params file

• initialise_constraint_vars : initialises the constraint variables to 0

In order to introduce a new method, the corresponding files should be added into the Source/Method folder (e.g. MyNewMethod.impl.hpp, MyNewMethod.hpp). In addition, the option has to be added to solver_main.cpp and the GNUMakefile, following the code for the existing examples (e.g., introduce USE_MYNEWMETHOD).

If new variables are added, they will need to be introduced in the rest of the code (for example, in MultigridUserVariables.hpp and BoundaryConditions.cpp), but the solver can solve any number of equations, provided that rhs, aCoef and bCoef are defined for all the variables to be solved.