May 27 2013

Abstract

This thorn implements a multigrid solver for systems of elliptic partial diﬀerential equations. It uses Carpet’s interface to handle loops over the grid hierarchy and pass information between diﬀerent components. The system of equations is passed to the solver via parameters pointing to the external grid functions that hold the equation coeﬃcients.

CT_MultiLevel is a Cactus thorn that implements a multigrid solver for elliptic partial diﬀerential equations (PDEs). The implementation is rather standard, and it also allows for local AMR grids (in which case it uses a rudimentary multilevel algorithm).

This thorn requires Carpet, which it uses to manage the access to the grid structure (via the BEGIN_*_LOOP and END_*_LOOP macros) and to pass information between the diﬀerent levels (via the restriction and prolongation operators).

Structurally, the solver can tackle any type of equation or system of equations, their coeﬃcients being passed to CT_MultiLevel via string parameters holding the names of the corresponding grid function. The thorn comes with a number of test cases which should be easy to customize and extend.

Multigrid schemes are designed to solve elliptic equations on a hierarchy of grids, with the goal of speeding up the convergence rate of relaxation algorithms (such as Gauss-Seidel) by eliminating diﬀerent frequency modes of the solution error on diﬀerent grids. It is known that relaxation techniques are very eﬀective at eliminating the high-frequency part of the error; when solving an elliptic PDE on a grid, it is then advantageous to relax it on a coarser grid beforehand and interpolate the result back onto the original grid. Complex schemes involving several grid levels, each with a diﬀerent PDE, can be designed to maximize the speed up. Please refer to [1] for the details.

CT_MultiLevel, along with the rest of the Cosmology arrangement, is released under the General Public License, version 2 and higher. The copyright of this thorn remains with myself.

CT_MultiLevel is publicly available in COSMOTOOLKIT’s git repository:

git clone https://eloisa@bitbucket.org/eloisa/cosmology.git

A helper thorn CT_Analytic is available under the same repository.

In order to solve an elliptic PDE with CT_MultiLevel, CT_MultiLevel can be simply compiled in, with the requirement that Carpet is compiled as well. The thorn also inherits from boundary and grid for boundary APIs and coordinate labels.

Currently, the equation is solved in a single schedule call at CCTK_INITIAL, in GLOBAL_LATE mode (to ensure that necessary objects such as the grid arrays have been populated).

As the equation coeﬃcients will be set by grid functions, a thorn allocating and setting these grid functions to the desired values has to be included. Another thorn in the Cosmology arrangement, CT_Analytic, serves this purpose. Grid functions of other origin are naturally just as good.

Each equation in the system to solve is parametrized as follows:

$$\begin{array}{rcll}{c}_{xx}{\partial}_{xx}\psi +{c}_{xy}{\partial}_{xy}\psi +{c}_{xz}{\partial}_{xz}\psi +{c}_{yy}{\partial}_{yy}\psi +{c}_{yz}{\partial}_{yz}\psi +{c}_{zz}{\partial}_{zz}\psi & +& & \text{(1)}\text{}\text{}\\ {c}_{x}{\partial}_{x}\psi +{c}_{y}{\partial}_{y}\psi +{c}_{z}{\partial}_{z}\psi +{c}_{0}{\psi}^{{n}_{0}}+{c}_{1}{\psi}^{{n}_{1}}+{c}_{2}{\psi}^{{n}_{2}}+{c}_{3}{\psi}^{{n}_{3}}+{c}_{4}{\psi}^{{n}_{4}}& =& 0& \text{(2)}\text{}\text{}\end{array}$$

All the coeﬃcients ${c}_{\ast}$ and the powers ${n}_{\ast}$ can be speciﬁed by setting the parameters *_gfname[*] to the desired grid function name. These parameters are arrays to allow for the solution of systems of equations.

The only other parameter which must be set to ensure correct operation is CT_MultiLevel::topMGlevel, which tells the solver which is the ﬁnest reﬁnement level that covers the entire domain in which the equation is to be solved. All levels below and including this will be used in the classical multigrid sense (e.g., in a V- or FMG-cycle). The ones above will be treated as local AMR boxes, and be solved via a multilevel prescription; presently, this involves simply interpolating the solution from the topMGlevel reﬁnement level at the end of the multigrid part and relaxing it progressively, from coarser to ﬁner, on all the local grids.

Other parameters control the stopping criteria (numbers of relaxation sweeps and residual tolerance), the ﬁnite-diﬀerencing order used to discretized the equation, the number of equations in the system, and whether the solution error should be calculated (this obviously requires the exact solution to be known and stored in a grid function, its name passed to CT_MultiLevel via the usual parameter mechanism).

Concrete examples on how to set these parameters are provided below.

When the equation coeﬃcients cannot be passed in via external grid functions (as is the case, for instance, when solving a coupled system of PDEs), then one can resort to auxiliary functions, which are recalculated after each relaxation sweep. The actual algorithm used to populate these coeﬃcient has to be provided by the user; some examples are provided under src/auxiliaries. These values are stored in CT_MultiLevel’s auxiliaries variable group.

All coeﬃcients, auxiliaries, and other useful grid functions are stored in the pointer structure coeffptr, to which the user can further append pointers to extra external grid functions, which are not used as PDE coeﬃcients (and thus do not need additional storage space and shouldn’t be introduced as auxiliaries), but are still needed by the solver. Examples can be found under src/extra.

Finally, if variable resetting is necessary after each relaxation step (for instance, when solving linear problems in periodic domains), a suitable calculation needs to be provided by the user, similar to the examples in src/integral.

CT_MultiLevel uses some of Carpet’s macros to access diﬀerent components of a grid hierarchy, and its prolongation and restriction operators to pass information between them.

Furthermore, CT_MultiLevel needs at least one external thorn to set the PDE coeﬃcients. The helper thorn CT_Analytic can be used to this purpose, but any set of external grid functions will serve the purpose.

The solver is described in [2], along with a number of examples. I rediscuss some of those below from the technical standpoint, highlighting which parameters need to be set in each case and whether other information is required from the user.

Poisson’s equation reads:

$$\Delta \psi +\rho =0$$ | (3) |

where $\rho $ is a known source function and $\psi $ is an unknown to solve for. Its Laplacian, $\Delta \psi $, simply represents the sum of its diagonal second-order derivatives:

$$\Delta \psi ={\partial}_{xx}\psi +{\partial}_{yy}\psi +{\partial}_{zz}\psi $$ | (4) |

so that cxx_gfname[0], cyy_gfname[0] and czz_gfname[0] all have to be set to one. The source term can be speciﬁed freely. CT_Analytic can be used to set all coeﬃcients, e.g. setting the CT_Analytic::free_data parameter to exact, where the coeﬃcients of the diagonal second-order derivatives are all one, the source term can be set to a linear combination of a gaussian, a sine, and a $1\u2215r$ term, and all the other coeﬃcients default to zero. To choose a gaussian source term, for instance, one can set:

CT_Analytic::free_data = "exact"

CT_Analytic::ampG = 1

CT_Analytic::sigma = 0.5

CT_Analytic::ampG = 1

CT_Analytic::sigma = 0.5

and then pass the relevant grid functions to CT_MultiLevel:

CT_MultiLevel::cxx_gfname[0] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[0] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[0] = "CT_Analytic::testczz"

CT_MultiLevel::c0_gfname[0] = "CT_Analytic::testc0"

CT_MultiLevel::n0[0] = 0

CT_MultiLevel::cyy_gfname[0] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[0] = "CT_Analytic::testczz"

CT_MultiLevel::c0_gfname[0] = "CT_Analytic::testc0"

CT_MultiLevel::n0[0] = 0

CT_Analytic can also be used to set the initial guess for the solver:

CT_MultiLevel::inipsi_gfname[0] = "CT_Analytic::testinipsi"

Any coeﬃcient that is not speciﬁed explicitly via a parameter will be initialized to zero.

Other parameters can be used to tune the solver’s behavior, including the cycling mode through the diﬀerent levels and the stopping criterion. The parameter ﬁle par/poisson.par provides a concrete example.

Let’s consider Einstein’s constraints in the Lichnerowicz-York form and in a conformally-ﬂat space, i.e. a coupled system of four PDEs for the four functions $\psi $ and ${X}_{i}$:

$$\begin{array}{rcll}\Delta \psi -\frac{{K}^{2}}{12}\phantom{\rule{0.3em}{0ex}}{\psi}^{5}+\frac{1}{8}{A}_{ij}{A}^{ij}{\psi}^{-7}=-2\pi \rho {\psi}^{5}& & & \text{(5)}\text{}\text{}\\ \Delta {X}_{i}+{\partial}_{i}{\partial}_{j}{X}^{j}-\frac{2}{3}{\psi}^{6}{\delta}^{ij}{\partial}_{i}K=8\pi {j}^{i}{\psi}^{10}& & & \text{(6)}\text{}\text{}\end{array}$$

where

$${A}_{ij}={\partial}_{i}{X}_{j}+{\partial}_{j}{X}_{i}+\frac{2}{3}\phantom{\rule{2.77695pt}{0ex}}{\delta}_{ij}{\partial}_{k}{X}^{k}$$ | (7) |

and $K$, $\rho $ and ${j}^{i}$ are freely speciﬁable functions.

Due to conformal ﬂatness, again $\Delta ={\partial}_{xx}+{\partial}_{yy}+{\partial}_{zz}$. However, in this case we have the additional complication that the coeﬃcient of the ${\psi}^{-7}$ term in the ﬁrst equation depends on ${X}^{i}$, and therefore has to be updated during the solution process. This is accomplished by storing this coeﬃcient in an auxiliary function, ﬁlled by the function CT_SetAuxiliaries initially and after each Gauss-Seidel iteration. The recipe to ﬁll this grid function is contained in src/auxiliaries/Lump.cc.

The parameter speciﬁcation for this case is then:

CT_Analytic::free_data = "Lump"

CT_Analytic::Kc = -0.1

CT_Analytic::ampG = 1

CT_Analytic::ampC = 1

CT_Analytic::ampV = 0.1

CT_Analytic::sigma = 0.2

CT_Analytic::vecA = 0.6

CT_Analytic::Kc = -0.1

CT_Analytic::ampG = 1

CT_Analytic::ampC = 1

CT_Analytic::ampV = 0.1

CT_Analytic::sigma = 0.2

CT_Analytic::vecA = 0.6

which speciﬁes the equation coeﬃcients and the free data $K$, $\rho $ and ${X}^{i}$, and:

CT_MultiLevel::inipsi_gfname[0] = "CT_Analytic::testinipsi"

CT_MultiLevel::cxx_gfname[0] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[0] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[0] = "CT_Analytic::testczz"

CT_MultiLevel::n0[0] = 5

CT_MultiLevel::c0_gfname[0] = "CT_Analytic::testc0"

CT_MultiLevel::n1[0] = 0

CT_MultiLevel::c1_gfname[0] = "CT_Analytic::testc1"

CT_MultiLevel::n2[0] = -7

CT_MultiLevel::c2_gfname[0] = "CT_MultiLevel::ct_auxiliary[0]"

CT_MultiLevel::inipsi_gfname[1] = "CT_Analytic::testinixx"

CT_MultiLevel::cxx_gfname[1] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[1] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[1] = "CT_Analytic::testczz"

CT_MultiLevel::n0[1] = 0

CT_MultiLevel::c0_gfname[1] = "CT_Analytic::testc2"

CT_MultiLevel::n1[1] = 0

CT_MultiLevel::c1_gfname[1] = "CT_MultiLevel::ct_auxiliary[1]"

CT_MultiLevel::inipsi_gfname[2] = "CT_Analytic::testinixy"

CT_MultiLevel::cxx_gfname[2] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[2] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[2] = "CT_Analytic::testczz"

CT_MultiLevel::n0[2] = 0

CT_MultiLevel::c0_gfname[2] = "CT_Analytic::testc3"

CT_MultiLevel::n1[2] = 0

CT_MultiLevel::c1_gfname[2] = "CT_MultiLevel::ct_auxiliary[2]"

CT_MultiLevel::inipsi_gfname[3] = "CT_Analytic::testinixz"

CT_MultiLevel::cxx_gfname[3] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[3] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[3] = "CT_Analytic::testczz"

CT_MultiLevel::n0[3] = 0

CT_MultiLevel::c0_gfname[3] = "CT_Analytic::testc4"

CT_MultiLevel::n1[3] = 0

CT_MultiLevel::c1_gfname[3] = "CT_MultiLevel::ct_auxiliary[3]"

CT_MultiLevel::cxx_gfname[0] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[0] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[0] = "CT_Analytic::testczz"

CT_MultiLevel::n0[0] = 5

CT_MultiLevel::c0_gfname[0] = "CT_Analytic::testc0"

CT_MultiLevel::n1[0] = 0

CT_MultiLevel::c1_gfname[0] = "CT_Analytic::testc1"

CT_MultiLevel::n2[0] = -7

CT_MultiLevel::c2_gfname[0] = "CT_MultiLevel::ct_auxiliary[0]"

CT_MultiLevel::inipsi_gfname[1] = "CT_Analytic::testinixx"

CT_MultiLevel::cxx_gfname[1] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[1] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[1] = "CT_Analytic::testczz"

CT_MultiLevel::n0[1] = 0

CT_MultiLevel::c0_gfname[1] = "CT_Analytic::testc2"

CT_MultiLevel::n1[1] = 0

CT_MultiLevel::c1_gfname[1] = "CT_MultiLevel::ct_auxiliary[1]"

CT_MultiLevel::inipsi_gfname[2] = "CT_Analytic::testinixy"

CT_MultiLevel::cxx_gfname[2] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[2] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[2] = "CT_Analytic::testczz"

CT_MultiLevel::n0[2] = 0

CT_MultiLevel::c0_gfname[2] = "CT_Analytic::testc3"

CT_MultiLevel::n1[2] = 0

CT_MultiLevel::c1_gfname[2] = "CT_MultiLevel::ct_auxiliary[2]"

CT_MultiLevel::inipsi_gfname[3] = "CT_Analytic::testinixz"

CT_MultiLevel::cxx_gfname[3] = "CT_Analytic::testcxx"

CT_MultiLevel::cyy_gfname[3] = "CT_Analytic::testcyy"

CT_MultiLevel::czz_gfname[3] = "CT_Analytic::testczz"

CT_MultiLevel::n0[3] = 0

CT_MultiLevel::c0_gfname[3] = "CT_Analytic::testc4"

CT_MultiLevel::n1[3] = 0

CT_MultiLevel::c1_gfname[3] = "CT_MultiLevel::ct_auxiliary[3]"

which informs CT_MultiLevel on the grid functions holding the coeﬃcients and on the existence of auxiliary functions.

Further parameters control other aspects of the solution, and can be found in par/constraints_spherical.par. In particular, the parameter CT_MultiLevel::number_of_equations has to be set explicitly to four to accomodate the full system of equations.