domain: math info phys

text:
Philippe Leleux
email: pleleux@laas.fr

Christina Schwarz
email: christina.schwarz@uni-bayreuth.de

Martin Joachim Kühn
email: martin.kuehn@dlr.de

Carola Kruse
email: carola.kruse@cerfacs.fr

Ulrich Rüde
email: ulrich.ruede@fau.de

Complexity analysis and scalability of a matrix-free extrapolated geometric multigrid solver for curvilinear coordinates representations from fusion plasma applications

Keywords: Multigrid, complexity, curvilinear coordinates, parallelization, multithreading, Plasma fusion

Mathematics Subject Classification

(2020) 68Q25 • 65Y20 • 65Y05 • 65N55 • 65N06 • 65B99 1 Introduction
With accelerating climate change, plasma fusion is a promising alternative for future energy productions. Tokamak fusion reactors are actively studied objects to realize energy production from plasma fusion; see, e.g., [WC11] and the references therein. In these torus-shaped reactors, plasma is magnetically confined to obtain energy from physical processes. Empirical laws are important to design current reactors such as ITER or JET1 . With the advent of modern supercomputers, simulation and high performance computing (HPC) have opened up new possibilities to study potential reactor designs. Simulation codes can lead to improved designs by creating a better understanding for physical aspects such as the transport of plasma. However, the simulation of plasma in these reactors still involves high cost in terms of computation and memory, and thus requires the use of efficient codes. The gyrokinetic framework is of particular interest in this context since it is able to capture the turbulence in the plasma, and allows to decrease the dimensions of the system from 6 to 5 dimensions: three for the torus geometry and two for the velocity [BBG + 18, GAB + 16]. In gyrokinetic codes, a 5D Vlasov equation coupled with a 3D Poisson-like equation must be solved at each time step. While being of smaller dimension, the solution of the latter 3D system is still computationally expensive and its compute time as well as the needed resources are non negligible. In current implementations such as Gysela [GAB + 16] or Jorek [HHP + 21], the 3D equation is reduced to the solution of 2D Poisson-like equations on poloidal cross-sections of the tokamak. The description of the domain by curvilinear coordinates [BBG + 18, ZG19, Zon19], and the consideration of a varying density field [ZG19,Zon19], lead to anisotropic operators and nonuniform meshes. In this case, a standard derivation of a finite difference discretization may lead to a non-symmetric linear system, although the corresponding energy functional is symmetric. In [KKR21], symmetric finite difference schemes based on minimizing the energy functional were derived. A rigorous analysis showed that, when using nonstandard finite element discretizations, extrapolation formulas [JR96] applied on the discrete energy functional lead to higher order approximations. Numerical results showed that additionally, the likewise extrapolated finite difference schemes lead to an order of convergence of four.

Multigrid solvers are among the fastest numerical methods to solve elliptic partial differential equations and are often optimal in the sense that the number of arithmetic operations is proportional to the degrees of freedom of the discretized problem [TOS01]. Multigrid algorithms can be divided into the class of algebraic and geometric algorithms. The focus of this paper is on a tailored geometric multigrid solver for a Poisson-like equation on disk-like domains. The extrapolation techniques presented in [KKR21] can be incorporated naturally in a multigrid scheme, as already shown by [JR96,JR98]. In [KKR22], this implicitly extrapolated geometric multigrid method based on symmetric finite difference schemes, and named GMGPolar, was considered for a simplified version of the gyrokinetic Poisson equation. The algorithm was shown to deliver a higher order convergence in only a small number of solver iterations.

When it comes to computational expense, two different foci are memory and computation. Algorithms are called memory-bound if their execution speed is limited by the speed of data transfer (between different levels of memory). Algorithms are called compute-bound if their execution speed is limited by the speed of the computational units, e.g., the cores of the computing node. For distributed-memory-based applications also the network speed can become a limiting factor.

In the present article, a more realistic version of the Poisson-like equation from [KKR22] is considered. We provide an efficient C/C++ implementation of GMGPolar that is available opensource. The published code can be found GitHub under a permissive Apache-2.0 license. Its release version 1.0 [LKS + 23] has been used for this article. We demonstrate optimal linear complexity of the solver with a systematic analysis of computational and memory costs. The algorithm employs a task-based multithreaded parallelism which takes advantage of the disk-shaped geometry of the problem. Such parallel implementation is essential for large scale simulations of plasma in tokamak reactors. In [BLK + 23], it has recently been shown that GMGPolar is able to compete with other state-of-the-art solvers for the considered applications. Our rigorous computation and memory analysis is built in the sense that it can be transferred to other stencil-based operators.

The article is structured as follows. In Section 2, we introduce the model problem, and in Section 3 we introduce the discretization used for the model problem and geometry. In Section 4, we briefly rephrase the implicitly extrapolated geometric multigrid solver, then in Section 5 we provide a detailed complexity analysis. The corresponding details on the particular function evaluations and our rigorous framework for stencil-based operators is to be found in the Appendix. The parallelization is discussed in Section 6 and numerical results can be found in Section 7.

Model problem

We are interested in the solution of the Poisson-like equation

-∇ • ( ∇ ) + = in Ω, = on Ω, (1)
where Ω ⊂ R 2 is a disk-like domain, : Ω → R, ∈ C 0 (Ω), is the right hand side, , : Ω → R with ∈ C 1 (Ω) ∩ C 0 (Ω) and ∈ C(Ω) are coefficients corresponding to density profiles. Furthermore, we have the Dirichlet boundary conditions ∈ C 0 ( Ω). For more details on this equation, see [Zon19,eq. (4.36)]. In [KKR22], Equation (1) was considered with = 0. Particular choices for the coefficients and will be presented for numerical experiments in Section 7.

As the discretization method presented in the next section will be based on the energy functional corresponding to the partial differential equation (1), we also briefly introduce this energy functional. The problem (1) can be traced back to the minimization of the energy functional

( ) := ∫ Ω 1 2 |∇ | 2 + 1 2 2 -d( , ), (2)
over a suitable Sobolev space incorporating the boundary conditions .

As in [BBG + 18, ZG19, Zon19], we describe our domains Ω by curvilinear coordinates, i.e. based on a mapping from the Cartesian coordinates ( , ) to the (generalized) polar coordinates or curvilinear coordinates ( , θ) ∈ [ min , max ] × [0, 2 ) where is the (generalized) radius and θ the (generalized) angle. Figure 1 shows an example of such coordinates and mapping. This description introduces a periodicity constraint in the θ-direction. Figure 1 provides an example of a circular domain described by a grid in Cartesian and polar coordinates. Please note that the mapping between both systems is only bijective if min > 0. In [KKR22], we have considered different approaches to handle the area around the origin. Similarly to [GAB + 16, KKR22], we consider that min ≈ 0, e.g. min = 10 -6 , and prescribe (homogeneous) Dirichlet boundary conditions on the artificial inner boundary. As presented in Figure 1, we describe a circular geometry by simple polar coordinates and an invertible mapping from Ω := [ , ] × [0, 2 ) onto the circle. In order to model more realistic cross-sections, generalized polar or curvilinear coordinates need to be used. Particular geometries have been introduced in [BBG + 18, CH08, ZG19]. We now introduce two mappings and which describe geometries by ( Ω) and ( Ω), see Figure 2. The Shafranov geometry is a deformed ellipse defined by the mapping

( , θ) := ( , θ) ( , θ) = 0 + (1 -) cos θ -2 0 + (1 + ) sin θ , (3)
where is the elongation and is the Shafranov shift, see [BBG + 18, ZG19]. In [KKR21, KKR22], this geometry was simply denoted deformed geometry, and 0 = 0 = 0 were removed from the equation.

The Czarny geometry adds triangularity to the shape with the mapping

( , θ) := ( , θ) ( , θ) = 1 1 -1 + ( + 2 cos θ) 0 + sin θ 2 -1 + ( + 2 cos θ) , (4)
where is the inverse aspect ratio, the ellipticity, and

= 1/ 1 -2 /4, see [CH08, ZG19].
In Figure 2, we present the Shafranov geometry for 0 = 0 = 0, = 0.3, and = 0.2 as well as the Czarny geometry for 0 = 0, = 0.3, and = 1.4.

Discretization

In the following, we introduce the mesh as well as the finite difference discretization with respect to the logical domain Ω := [ , ] × [0, 2 ). The mesh Ω 1 will be in product (i.e., tensor) format and given by two separate subdivisions 1 , . . . , ∈ [ min , max ] and θ 1 , . . . , θ θ +1 ∈ [0, 2 ]. These subdivisions are arbitrary with definitions 1 = min , = max , θ 1 = 0, and θ θ +1 = 2 ;

+1 = + ℎ , ∈ {1, . . . , -1}, θ +1 := θ + , ∈ {1, . . . , θ }, (5)
ℎ 2 = ℎ 2 -1 , ∈ {1, . . . , ( -1)/2}, θ even, 2 = 2 -1 , ∈ {1, . . . , θ /2}. (6)
Remark 1 The assumptions in Equation ( 6) mean that there is a specific structure needed on pairs of rectangular elements. The assumption ( 6) is needed to prove the relations between quadratic and extrapolated linear finite elements in [KKR21] which increase the order of convergence of the resulting multigrid scheme. There is no proof for finite differences, but it was shown numerically on these grids that they behave similarly as their finite element counterpart. Note that this assumption is not needed if the discretization is used for geometric multigrid without the proposed extrapolation step. For more details, see [KKR21].

We further assume ℎ and to be uniformly bounded by

0 < ℎ min ≤ ℎ ≤ ℎ and 0 < min ≤ ≤ , (7)
as well as the existence of 0 < < ∞ such that ℎ = . Our finite difference discretization builds upon the localization (11) of the energy functional (2) first, the discretization of this localized form second, and the differentiation of the sum of the localized, discretized energies third.

In order to localize Equation (2), consider any rectangular element := [ , +ℎ ]×[θ , θ + ] of the logical domain Ω. For , , we have its curvilinear representation by ( ). Here, typically corresponds to the Shafranov or Czarny geometry but could be given by another invertible mapping.

The following paragraphs are essentially built upon the theory in [KKR21] and only differ due to the coefficients and . In [KKR21], we allowed for a more general diffusion tensor : R 2 → R 2 instead of : Ω → R which could, in principle, also be used here. Also, in [KKR21], we used = 0. In order to make the paper self-contained, we repeat some of the basic ideas.

To make the transformation between Cartesian and curvilinear coordinates clear, let us use

( , ) = ( , θ), ( , ) = ( , θ), ( , ) = ( , θ), ( , ) = ( , θ). ( 8
)
By transformation of the energy functional in Equation ( 2), we then obtain

, ( ) : = ∫ ( , ) 1 2 |∇ ( , ) | 2 + 1 2 2 - d( , ) = ∫ , 1 2 | -∇ ( ,θ) | 2 + 1 2 2 - | det | d( , θ), (9)
where is the Jacobian matrix of and -:= ( ) -1 . In principle, we have to distinguish the functions in Equation ( 8) as well as the operators with respect to the different coordinates throughout the paper but we will generally write or ∇ and assume that the variables are clear from the context.

In order to simplify the notation, we define

The symmetric finite difference stencils are now derived by discretizing , ( ). Here and in the following, we use for any function : Ω → R the notation

, := ( , θ ). ( 13
)
The discretization of

∫ , 2 + θ θ + θθ 2 θ -| det | d( , θ) (14)
can be obtained by [KKR21,(3.15) and (3.8)]. We only need to consider the remaining part

∫ , 1 2 2 | det | d( , θ). (15)
Defining := 1 2 2 | det |, the discretization of Equation ( 15) can be obtained by using the trapezoidal rule:

∫ , d( , θ) = ℎ 4 ( +1, +1 + +1, + , +1 + , ) + O (ℎ 3 + ℎ 3 ). ( 16
)
Remark 2 Alternatively, we could use the midpoint rule and operators (3.3) and (3.4) of [KKR21] to obtain the same result.

The discretization [START_REF] Grandgirard | A 5D gyrokinetic full-f global semi-Lagrangian code for flux-driven ion turbulence simulations[END_REF] as well as the discretizations of the two summands in Equation ( 14) are linear and quadratic expressions in := ( 1,1 , . . . , , θ+1 ) with a local discretization error of O (ℎ 3 + ℎ 3 ) (see [KKR21,(3.15) and (3.8)]). The sum over all elements then implicitly yields a discrete energy operator as well as a vector for | det | such that

( ) = -1 =1 θ =1 , ( ) + O (ℎ 2 + 2 ) = 1 2 - + O (ℎ 2 + 2 ), (17)
where , ( ) stands for the discretization of , ( ) in Equation ( 11). Instead of searching the minimum of the continuous energy functional ( ), we now compute the derivative of the discretized formulation 1 2 -with respect to the nodal values of and set it to zero, i.e., ,

1 2 - = -1 =1 θ =1 , , ( ) ! = 0. ( 18
)
Due to the difference, averaging, and integration formulas used for discretization, we have

, , ( ) = 0, for ( , ) ∉ , := {( , ), ( -1, ), ( , -1), ( -1, -1)}.
Computing the derivatives for ( , ) ∈ , and solving Equation ( 18) yields the nine point stencil with the horizontal and vertical connections

+1, : ( * 9) +1, = - + -1 ℎ ( , ) + ( + 1, ) 2 , -1, : ( * 9) -1, = - + -1 ℎ -1 ( -1, ) + ( , ) 2 ,
, +1 : ( * 9) , +1 = -

ℎ + ℎ -1 θθ ( , ) + θθ ( , + 1) 2 , , -1 : ( * 9) , -1 = - ℎ + ℎ -1 -1 θθ ( , -1) + θθ ( , ) 2 , (19)
the central update

, : ( * 9) , = -( * 9) +1, + ( * 9) -1, + ( * 9) , +1 + ( * 9) , -1 + (ℎ + ℎ -1 ) ( + -1 ) 4 , | det , |, (20)
and the diagonal connections

+1, +1 : ( * 9) +1, +1 = θ ( + 1, ) + θ ( , + 1) 4 , +1, -1 : ( * 9) +1, -1 = θ ( + 1, ) + θ ( , -1) 4 , -1, +1 : ( * 9) -1, +1 = θ ( -1, ) + θ ( , + 1) 4 , -1, -1 : ( * 9) -1, -1 = θ ( -1, ) + θ ( , -1) 4 . (21)
Finally, the right hand side is given by

(ℎ + ℎ -1 ) ( + -1 ) 4 , | det , |. (22)
4 Geometric multigrid solver Multigrid solvers are powerful numerical methods for solving partial differential equations; see, e.g., [TOS01]. In our paper, we will focus on a tailored geometric multigrid solver for Equation (1). In [KKR22], such a solver was proposed using two ingredients: a suited smoothing procedure for curvilinear discretizations, and an implicit extrapolation scheme. The smoothing procedure is motivated by [Bar88] while the implicit extrapolation scheme was introduced in [JR96, JR98]. The resulting solver is named GMGPolar.

To make this paper self-contained, we briefly introduce the solver used in the following numerical simulations. For more details on the derivation of the developed multigrid scheme, see [KKR22].

Multigrid hierarchy

In order to build a multigrid scheme, we need to construct a hierarchy of grids. Based on Equation (5), the finest grid is defined as a rectangular mesh Ω 1 of the logical domain, see also Figure 1. We then use standard coarsening to obtain coarser meshes. Basically, starting from the first node, we take every second value in { 1 , . . . , } and {θ 1 , . . . , θ θ +1 } to build the second grid level Ω 2 , and so on and so forth. Note that for all levels, the last values are naturally kept in each direction ( and θ θ +1 ).

Assuming that there exist constants , θ ∈ R such that = 2 + 1 and θ = θ 2 for ∈ N, we can construct a hierarchy of grids Ω 1 , . . . , Ω from finest to coarsest. On each level ∈ {1, . . . , }, we then use the discretizations ( 19)-( 22) to construct the system

= , (23)
where is a matrix of size × and , are vectors of size . On level , with and θ + 1 the number of divisions in the and θ directions, we have × ( θ + 1) nodes in the grid. In the θ-direction, we add the periodicity condition θ 1 = θ θ+1 and reduce the size of the matrix to = • θ . Since standard coarsening is used on the 2D grid, we have that

4 ( -1) ≈ 1 . ( 24
)
The levels are linked by two transfer operators, the prolongation operator +1 transferring the information from the coarser level + 1 to the finer level , and the restriction operator +1 , which is often defined by +1 := +1 . The latter definition ensures that the operator

+1 +1 +1 (25)
is symmetric when is symmetric. In [KKR22], the prolongation is defined as the bilinear interpolation operator on an anisotropic mesh. In Table 1, we provide an illustration of the stencil used for the interpolation +1 . For more details, see Section 4.3. Except for the implicit extrapolation, which will act only between the finest two levels, these stencils are qualitatively identical on any level, only depending on the prolongated node. Quantitatively, the values of the operators depend on the distance to the neighboring nodes.

In order to avoid a proliferation of indices, we skip the index for a particular grid where our remarks hold for any level of the hierarchical grid Ω 1 , . . . , Ω .

Prolongation Stencil Reference

Coarse Vertical Horizontal Diagonal central node ( , ) odd, odd even, odd odd, even even, even

+1 (interpolation) (injection) None None None (triangular)
Table 1: The stencils used for the prolongation operators depending on which node is prolongated.

The four types of central nodes depend on their neighboring coarse nodes, and are labeled "coarse", "vertical", "horizontal", and "diagonal". White circles are coarse nodes, black dots are fine nodes. A line corresponds to a relation with the central node in the prolongation stencil. Also, we give the reference element for each prolongation type: finite difference for +1 and , P 2 finite element for .

Smoothing procedure

Pointwise smoothing leads to slow convergence of the multigrid scheme when polar or curvilinear coordinate formulations are considered with standard coarsening [TOS01, Chap. 5.1.1]. [Bar88] proposed the use of circle-and radial-line smoothing operations each using two alternating colors (black and white) for every line on the full unit disk or annulus. On an annulus, for compact stencils (i.e., stencils of a maximum length of one in each direction), all lines of one color of a particular smoother (i.e., circle or radial) can be considered in parallel. On full disks, special care has to be taken for radial smoothing as its colors get connected by the origin; see also Figure 4 for a motivation. Due to the anisotropy introduced by the coordinate transformation, circle and radial smoothing behave differently on different parts of the computational domain. This means that smoothing each node with both smoothers, i.e. first all nodes with the first smoother and then all nodes with the second one, may be inefficient in terms of computational cost. Based on the results of [Bar88], in [KKR22], it was numerically shown that one smoothing step per grid node (either radial or circle) can be sufficient to obtain fast convergence. Therefore, we partition the domain into two subdomains where either circle or radial smoothing is applied. The criterion to change from circle (starting with the innermost circle line) to radial smoothing is to find a minimum such that it exists a with ℎ > 1;

(26) see Figure 3 for an illustration. We define as the index of the outermost radius of the interior subdomain. In Figures 4 and3, we have = 4 and the corresponding circle would be given by the nodes with global numbers 25 to 32. Then the number of nodes in a black (or white) radial line in the exterior subdomain is -= 4, e.g., nodes {36, 44, 52, 60} in Figure 4. The line smoother used is a simple block Gauss-Seidel algorithm and these results also hold for more general geometries. The sequential execution of the two partial smoothers (circle and radial), considering the update from the first one in the second one, reduces the iteration numbers. Moreover, a partially parallel implementation can substantially reduce the computation time without increasing the iteration number. For more details, see [KKR22, Sec. 4.1] and Section 6.

We further introduce the double index to refer to a subset of nodes used in partial smoothing operations. We use ∈ {C, R} with for circle and for radial as well as ∈ {B, W} with for black and for white. For instance, with = R and = B, we have the square submatrix

(27)
of , corresponding to a restriction on the black colored nodes where radial smoothing is applied.

In Figure 4, we provide a small example with the nodes highlighted with different background colors, and the corresponding nonzero pattern of the submatrix . Note that can easily be reordered as a block diagonal matrix when using, e.g., compact five, seven, or nine point stencils. In this case, all lines can be handled in parallel. Similarly to , the part of the solution vector (resp. right hand side) corresponding to the nodes on the lines is denoted by (resp. ). Finally, we define

:= / and := / , (28)
the ratios of unknowns touched by the partial smoothing against the total number of unknowns on the current grid. Here, represents the number of unknowns treated by the partial smoothing procedure (both colors for either circle or radial smoothing steps), while refers to the subset of unknowns treated by only one color . Similarly, we can define := / as the number of unknowns of color (either black or white), and we have ≈ ≈ /2. For the consideration of the complexity of the overall algorithm, we break all the ingredients down to the smallest possible unit. Therefore, let us further consider a submatrix of the previously introduced submatrix . For both partial smoothers, a grid line refers to one line that is smoothed together. In Figure 4 (left), we have four radial black grid lines, highlighted by the different background colors yellow, blue, red, and green. Let us denote the submatrix corresponding to the -th grid line

( ) . (29)
For instance, we obtain the matrix (0) by extracting all yellow squares from Figure 4 (right).

Similarly, we write ( ) for the corresponding right hand side part. The matrix ( ) ⊥ then consists of all remaining elements, which are the rows corresponding to the -th grid line, and all columns that do not belong to ; that is, either they belong to the other color of the same partial smoother or any color of the other partial smoother. One relaxation step for this line is then expressed as

( ) = ( ) -1 ( ) -( ) ⊥ ( ) ⊥ . ( 30
)
Table 2 gives the stencils of and ⊥ for each type of partial smoother. Finally, a complete relaxation step in GMGPolar applies both partial smoothers ∈ { , } with both colors, black ( ) and white ( ). Different colors are treated sequentially while all lines of one color are smoothed in parallel. That means, e.g., that the scheme of the second radial color uses the update from the Gauss-Seidel scheme of the first radial color.

Multigrid iterations and implicit extrapolation

Using the previously introduced hierarchy of systems and the tailored smoother on each level, GMGPolar follows a traditional -cycle to solve the discretized system. Algo. 1 gives the whole -cycle process, where S ( , , ) returns the result of relaxation steps applied onto on level .

Algorithm 1 MG( , ): Multigrid V-cycle with 1 pre-smoothing and 2 post-smoothing steps on level ( levels with = 1 the finest)

1: if coarsest level then 2:
return -1 (exact solve) 3: else 4:

← S 1 ( , , ) (pre-smoothing) 5:

+1 ← +1 ( -) (restriction of the residual) 6:

+1 ← MG( +1 , + 1) (coarse grid correction) 7:

← +1 +1 (interpolation of the error) 8:

← + (update of the solution) 9:

return S 2 ( , , ) (post-smoothing) 10: end if

In [KKR21], the previously introduced variant of GMGPolar showed an approximation error of degree 2 (compared to a manufactured solution of the Poisson-like problem (1) with = 0) in a modest number of iteration steps. In order to increase the approximation order of the algorithm from order 2 to order 3 or 4, the authors in [KKR22] presented a mechanism called implicit extrapolation [Rüd91, JR96] which can naturally be implemented in the multigrid scheme. The implicitly extrapolated multigrid is based on the -extrapolation technique [Ber97]. For more details, see the motivation in [Sch21, Sec. 4.5] and the references therein. This approach allows to increase the approximation order without the explicit construction of a higher order discretization, and thus, keeps the modest computational cost of the original multigrid method GMGPolar.

The implicit extrapolation uses additional information from the next coarser grid to estimate the local error, and to eliminate certain dominating error terms on the finest level. Due to the suitable combination of lower order systems when computing the new residual, we can implicitly achieve a higher order solution. Note that with implicit extrapolation, we only relax the fine nodes on the finest level, ensuring that the multigrid iteration converges to a unique higher order solution. For more details, see [JR96].

For this paper to be self-contained, we briefly recall how the components of the multigrid cycle are affected by this method. The extrapolation technique is only applied between the finest two grids. Without loss of generality, let us assume that the unknowns are ordered such that the coarse nodes are ordered first and the remaining fine (or fine-only) nodes second.

As the extrapolation only affects the first two levels, let us drop the index of = or explicitly set it to = 1 for (31)-(33). The system = then writes

⊥ ⊥ ⊥ ⊥ ⊥ = ⊥ , (31)
where the indices • and • ⊥ , respectively correspond to fine nodes which do not belong to the underlying coarse mesh, and the coarse nodes. The multigrid cycle with implicit extrapolation then differs in two aspects. First, and contrary to the -extrapolation [Ber97], the smoothing is only performed on the unknowns on the fine grid excluding the coarse unknowns ⊥ ; for more details see [JR96]. Thus, a relaxation step follows Equation (30) but, for each line, the matrix is restricted to the fine unknowns , and ⊥ also affects the coarse unknowns ⊥ of the lines , see Table 2. We denote this smoothing procedure by

S 1 ( , , ) := ⊥ S 1 ( , , - ⊥ ⊥ ). (32)
Secondly, we use a different intergrid transfer operator between the two finest levels. The restricted residual on the first level ( = 1) of the implicitly extrapolated scheme is given by

:= 4 3 ( 1 -1 1 ) - 1 3 ( 2 -2 1 ). ( 33
)
Here, is the simple injection operator such that extracts the coarse nodes from 1 . The operator is specific to the extrapolation approach. is based on triangle finite element theory and for more details, see [KKR22, (4.8)] as well as [KKR21]. It is the identity operator on the coarse nodes. For the remaining fine nodes, it returns 1/2 of the values of the adjacent coarse nodes in the triangle finite element. Table 1 gives the stencil corresponding to and .

Configuration Stencil

⊥

Circle smoother

Radial smoother

Extrap. fine nodes *

Table 2: The stencils used for the smoother matrices and ⊥ . A black dot corresponds to a relation with the central node in the stencil. * : "Extrap. fine nodes" refers to the stencil used when smoothing fine nodes (of the same color as the coarse nodes) with the implicitly extrapolated scheme, identically for both circle and radial smoothing.

The complete implicitly extrapolated multigrid cycle is given in Algo. 2. In the context of finite elements, it can be shown that the discretization, with linear elements and extrapolation, is equivalent to a quadratic discretization when nonstandard integration rules are applied. For more details, see [Rüd91, JR98, KKR21]. In practice, approximation orders between 3 and 4 are observed for the finite difference scheme introduced in [KKR21, KKR22], where the authors use = 0 in Equation (1). In Section 7, we demonstrate similar approximation orders using a general value for the coefficient , and the discretization from Section 3.

Algorithm 2 iexMG( 1 ): Implicit extrapolation applied on the finest level = 1 with 1 presmoothing and 2 post-smoothing steps

1: 1 ← S 1 ( 1 , 1 , 1 )
(pre-smoothing fine nodes)

2: ← 4 3 ( 1 -1 1 ) -1 3 2 -2 1 (extrapolated residual) 3: 2 ← MG( , 2) (coarse grid correction) 4: 1 ← 2 (interpolation of the error) 5: 1 ← 1 + 1 (update of the solution) 6: return S 2 ( 1 , 1 , 1 )
(post-smoothing fine nodes)

Computation and complexity

Numerical experiments are an essential tool to allow considerations which are difficult or impossible to realize by physical experiments. However, when it comes to the application of numerical algorithms on supercomputers, their cost in terms of computations and memory can also become important. In order to keep the cost reasonable and to make best use of the resources, these algorithms should be designed with care. This means that potential bottlenecks should be investigated and possibly removed or reduced. Here, we solve the related problem (1) with a target size between millions and a billion of dofs, mainly with a shared memory approach. However, many hundreds or thousands instances of these problems have to be solved in each time step of a gyrokinetic code such as Gysela [GAB + 16].

Multigrid methods can be asymptotically optimal, in the sense that the complexity to solve a linear system with sufficient accuracy grows only linearly with the number of unknowns . This has been shown for the full multigrid method [BL11]. Let us recall that GMGPolar follows a matrix-free implementation in order to decrease the memory-footprint, thus a problem of a size up to a billion can be solved in a single computing node, at the expense of a higher computational cost. In this section, we provide the theoretical complexity of GMGPolar in terms of computation and memory consumption. In particular, we show that GMGPolar inherits an optimal complexity. The computation of such complexity is rarely done for complete solvers. This section, completed by the appendix Appendix 8.4.1, can serve as a general approach for other frameworks.

Remark 3 The computational complexity of our method is expressed here in terms of floating point operations (flops). Although we are aware that different operations do not correspond to the same number of cycles on an actual machine, e.g. the addition or the division, we consider here that they are all equivalent. This number of flops still gives a fair idea of how expensive an algorithm is.

General considerations

Let us consider = × θ , the number of rows and columns of the matrix, as well as the number of nodes (i.e., degrees of freedom) in the grid, on the finest level. Our goal here is to obtain an estimation of the asymptotic complexity with respect to , in terms of both memory and computation. We denote by (X) the computational cost of the operation X, and by (V) the memory cost of the structure (or the function) V. Since we are interested in the asymptotic complexity of GMGPolar, we simplify the calculus by neglecting the boundaries (interior radius , exterior radius , as well as the boundary between radial and circle-line smoothers). We also neglect the manipulation of integers, to only focus on floating point numbers and their operations of the order of O ( ), ∈ N * .

In order to compute the complexity on the whole grid, we often consider the complexity per node (X) and (V). If the complexity applies to all nodes likewise (except the effect of boundaries), we have for the overall complexity

(X) = O (( (X) ) and (V) = O ( (V) ) . ( 34
)
Let us start with an example: As input, we need to have a right hand side (RHS) vector of size defined on the finest level. We consider it as an input vector of size . Then, ( ) = 0, and the corresponding memory consumption is (RHS) = . Also, we need to store the mesh information, i.e. the grid, on all levels. As we use a tensor-format mesh, the grid is represented by two arrays of respective size and θ containing the divisions in the -and θ-directions.

For a non-degenerated grid, we assume ≈ θ ≈ √ . Because the size of the -and θ-arrays are of order √ , we choose to neglect them to estimate the asymptotic complexity and memory consumption. In order to decrease the computations, we also consider that we have precomputed the values of cos(θ) and sin(θ) in the setup phase (only dependent on θ) as well as the values of the coefficients and from Equation (1) (only dependent on ). The corresponding number of flops and memory consumption is of the order √ , and thus neglected.

Let us now consider a matrix , which entries can be written as the sum of several factors, i.e. = ( ) . The computational cost to build , denoted by ( ), then includes the cost to build each factor ( ) , as well as an addition to add this factors into . Furthermore, considering a vector and the vector = , we have for each entry = = , ( ( ) ) . The computational cost to build includes one addition and one multiplication for each factor ( ) composing the matrix. In a matrix-free approach, each entry of the matrix must be recomputed on-the-fly at each application, since the matrix is not kept in memory. As a result, the computational cost for the matrix-free application of can be decomposed in the cost to construct the matrix (including one addition per factor), plus a cost to apply the matrix which corresponds only to the multiplication by the entries of (the addition is already accounted for). We write that

( - ) = ( ) + ( ), (35)
where ( ) is simply the number of multiplications by , i.e. the total number of different factors ( ) composing the matrix. This result is true for the matrix , as well as the prolongation and smoothing operators.

We now study the complexity of the three main components performed in the multigrid cycle detailed in Algo. 1 and 2: the computation of the residual (including the restriction), the smoothing procedure, and the coarsest grid solve. We first consider that the implicit extrapolation is turned off, and will detail the differences in Section 5.5. In order to lighten this very technical section, we only provide the main steps in order to compute the final complexity here. The whole computation is added as an appendix. The final complexity estimates are given in Table 3.

Computation of the residual

The computation of the residual mainly involves the matrix-free application of the matrix and the restriction .

Matrix-free application of

As introduced in Section 3, the matrix is based on a finite difference scheme with a nine point stencil in GMGPolar. When applying , nine entries must be computed for each node in the grid. These entries imply the computation of the three functions , θθ , and θ on the current node and its neighboring nodes, from Equation ( 19), which we generically denote by . In order to apply , we distinguish two approaches, which are illustrated in Figure 5 for the computation of a stencil. In a naive approach, here called -, the algorithm loops over all the grid nodes and computes the entries of the stencil directly, including the values for the neighboring nodes. Overall, the values for and θθ must be computed three times for each Table 3: Cost and number of applications for the functions used in the multigrid cycle. For each group of functions (surrounded by horizontal lines), we provide a summary statistics in bold face.

If memory gets reused, the summarized memory is less than the sum of the previous rows. For computations, the summary statistics is the sum of the previous rows. 1 and 2 are the number of pre-and post-smoothing sweeps. We denote the number of iterations for the convergence of the multigrid cycle.

Function Memory Computation

Applications Levels

right hand side input - = 1 matrix-free O ( [ ( ) + 99] ) = 1, • • • , -1 matrix-free O ( [7.25] ) residual 3 O ( [ ( ) + 107.25] ) = 1, • • • , -1 build O (9 ) O ( ( ) + 74) 1 = factorization O 2 /2 O 3 /3 1 coarse setup O 2 /2 O 3 /3 + ( ) + 74 1 = compute -1 O 2 2 - = coarse apply O 2 2 - = build 3 O ( [2 ( ) + 67] ) 1 = 1, • • • , -1 factorization O ( [3 + 2 ] ) O ( [3 + 7 ] ) 1 smoother setup O ( [3 + 2 ] ) O ( [2 ( ) + 70 + 7 ] ) 1 = 1, • • • , -1 compute -1 O ( [5 + 4 ] ) ( 1 + 2 ) = 1, • • • , -1 matrix-free ⊥ ⊥ O ( [2 ( ) + 74] ) ( 1 + 2 ) smoother apply 2 O ( [2 ( ) + 80 + 4 ] ) ( 1 + 2 ) = 1, • • • , -1
node, i.e., for the stencil of the current node and the two neighbors respectively from the same and θ. The values for the θ must be computed four times. In another approach, called -, we loop over all nodes to compute their local values only once, and update the entries in the matrix corresponding to the neighboring nodes. Thus, we directly divide by around three the number of computations for the functions . The functions are composed of similar elements, and we further decrease the computational cost by computing them simultaneously for each node.

Remark 4 Based on the stencil definition from Equations ( 19)-( 21), the θ functions used for the diagonal updates originate from the neighbors with the same or θ as the current node, thus the curved arrows in Figure 5. In other words, we never need θ function evaluations at the diagonal neighbors' positions. In the rest of the paper, we simplify the representation by drawing diagonal arrows for these relations.

Remark 5 In a finite element context, the A-give approach would correspond to looping over the grid nodes and adding in the matrix the contributions of surrounding elements. The A-take approach amounts to looping over the elements to update at once all the entries of the stiffness matrix corresponding to its vertices.

Remember that the values of cos(θ), sin(θ), , and are precomputed and their cost is asymptotically negligible as they only depend on either θ or . The core element of the functions is the computation of the Jacobian of the transformation , whose fixed cost depends on the chosen geometry. The computation of the three functions can then be obtained for ( ) = 19 + ( ) flops per node, see Appendix 8.2. If a subset of the is not needed, its specific cost is subtracted from the total sum (4 flops for or θθ , and 5 flops for θ ) but there remains an irreducible common cost of 6 flops corresponding to .) Additionally, by computing at once the of a whole line, with a fixed radius, it is possible to benefit from the vectorization of modern computing architectures in order to increase the computing intensity.

Adding the small cost from computing the scaling coefficients in Equations ( 19)-( 21) (see Appendix 8.3) gives the complexity of building the matrix , i.e.,

( ) = ( ) + 74 flops. ( 36
)
As explained in Section 5.1, in order to apply to a vector in a matrix-free fashion, each factor composing the entries of is multiplied with the corresponding vector value resulting in ( ) = 25 flops. Finally, with (35) and its surrounding explanation, we get

( - ) = ( ) + ( ) = ( ) + 99 flops. ( 37
)

Prolongation operator

The prolongation operator used in GMGPolar is a bilinear interpolation handling non-uniform grids. The stencil of this operator differs depending on the type of the prolongated node, i.e. depending on the neighboring coarse nodes. As defined in Table 1, there are four types of nodes which updates are defined as

( , θ ) =                      +1 ( , θ ) (coarse node injection), if odd and odd, 1 + -1 +1 ( , θ -1 ) + -1 +1 ( , θ +1 )
, if odd and even,

1 ℎ +ℎ -1 ℎ +1 ( -1 , θ ) + ℎ -1 +1 ( +1 , θ ) , if even and odd, 1 (ℎ +ℎ -1 ) ( + -1 ) ℎ +1 ( -1 , θ -1 ) +ℎ -1 +1 ( -1 , θ +1 ) + ℎ -1 +1 ( +1 , θ -1 ) +ℎ -1 -1 +1 ( +1 , θ +1 ) ,
if even and even.

(

) 38
By simple calculus, for a fine node ( , θ ), the number of flops in order to build and apply are respectively

( ) =            0, 4, 4, 12, and
( ) =           
1, if odd and odd, 2, if odd and even, 2, if even and odd, 4, if even and even.

(39)

Since each type of update concerns around /4 nodes, we finally obtain that

( ) = 5 and ( ) = 2.25, ( - ) = ( ) + ( ) = 7.25. ( 40
)
As for the restriction operator, i.e. = , we have the same complexity. The only difference is that we loop over the fine nodes of the grid to apply the prolongation, and we loop over the coarse nodes to apply the restriction. This way, the update of a specific node is completely independent from the others, which gives a natural and efficient parallelization.

Full residual expression

Finally, the restricted residual is expressed as ( -). In terms of memory consumption, only the right-hand side, the solution and the residual vectors are needed (the restricted residual becomes the right-hand side on the next coarser level). We then have

( ) = 3 . (41)
And, by taking into account the difference between vectors and in the residual (1 flop per node), we obtain the final complexity ( ) = ( -) +1+ ( -) for the restricted residual ( -), as given in Table 3.

Smoothing

We now turn to the application of the smoother, i.e., determining how expensive is the application of Equation (30) to each line. The matrices and ⊥ are submatrices of , up to some reordering of the rows and columns. We construct and/or apply them following an Agive approach. Thus, we obtain similar computational complexities for their construction and application as for , see Appendix 8.4.1 for the full details.

Based on the stencils from Table 2, not all of the functions need to be computed for the operators and ⊥ . The matrices are fully assembled in a setup phase in order to apply -1 during the smoothing steps. In this case, only the and θθ functions need to be computed. As for the application of ⊥ , the θ functions are needed for all nodes, however, either the computation of the or the θθ functions are required for the circle and the radial partial smoothers respectively.

Again, it is possible to group the computation of the functions and use vectorization. Although for the matrix , whole circle lines were computed at once, here we group circle lines or radial lines respectively for each partial line smoother.

In order to build the matrix (resp. apply ⊥ ), we need to loop over the nodes corresponding to the partial smoother but also over their neighbors, i.e. the nodes of the same partial smoother but with the other color. For example, in an A-give approach, to apply the smoother , we also loop over the nodes . Then, after applying all the partial smoothers, we have looped over all the grid nodes twice, implying to compute the Jacobian DF of the mapping twice. Additionally, while the stencil of the is greatly reduced, a major part of the computations from is kept because the central update of the stencil is the sum of the left, right, top, and bottom updates. Thus, the computational cost corresponding to the application of ⊥ is similar to the one for the construction of , with only the additional diagonal updates. In the end, setting aside the computation of the Jacobian which must be performed twice per node, the computational cost to build and to apply ⊥ represent respectively 2/3 and 3/4 of the operations required to build , see Table 3 andAppendix 8.3. Since the computational cost related to and θ is the same, the cost of the circle and radial-line smoothers are identical. In fact, in terms of computational cost, the stencil of the radial partial smoother is the same as a 90 • rotation of the stencil for the circle partial smoother with.

Once has been built, the relaxation step is applied at each iteration, for each type of partial smoother. Considering each line of the smoothers separately as in Equation ( 30), the structure of the corresponding matrix ( ) is specific to the partial smoother used:

for the circle-line smoother, ( ) is a tridiagonal matrix with periodicity condition, see Figure 6a, -for the radial-line smoother, ( ) is a tridiagonal matrix, see Figure 6b. In the actual GMG-Polar, the Dirichlet boundary condition would also appear in the form of a 1 as first entry in the diagonal but we neglect this aspect for the computation of the asymptotic complexity. Using a modified Thomas' algorithm, we then compute and store an LU decomposition of the matrices. For one line of unknowns, the computational cost of the factorization is then respectively O (10 ) and O (3 ) for the circle and the radial-line smoothers. See Appendix 8.4.2 for the algorithms. The total computational cost to factorize the matrices is obtained by summing on all lines and all partial smoothers. Considering and the ratio of unknowns associated with each type of smoother with respect to , we obtain

( . ) = 10 + 3 = 3 + 7 , (42)
since + = 1. The latter cost is an abstraction giving the average cost per node to compute the LU factorization. While there is no fill-in associated to the factorization of the matrix for the radial smoother, the last row and column of the matrix will be filled for the LU decomposition of the circle-line smoother. The LU decomposition of the matrices are stored in a sparse matrix format, replacing the original in memory. Accounting for the fill-in, the associated memory consumption is ( .

) = 3 + 2 . ( 43
)
At each relaxation, simple forward and backward substitutions are then applied to solve Equation (30), based on the computed LU factorizations, for a cost of O (9 ) and O ( 5) respectively for the circle and radial smoothers. The result is

( -1 ) = 9 + 5 = 5 + 4 , (44)
see Appendix 8.4.2 for the algorithms. Thus, we obtain the final complexities of the smoothing steps: ( ℎ ) (computed once), ( ℎ ) (at each iteration). As for the corresponding memory consumption, we have

( ℎ ) = ( . ), then ( ℎ ) = 2 corresponding to the two vectors -1 , ⊥ ⊥ . Overall, the smoothing procedure costs ( ℎ ) = ( ℎ ) + ( ℎ )
, as shown in Table 3.

Solving on the coarsest grid

On the coarsest grid of the multigrid hierarchy, GMGPolar uses a direct solver at each iteration to compute the coarse grid correction through the solution of the equation

= , (45)
where is the coarsest grid operator of size , and is the restricted residual transferred from the finer grid level -1. For this purpose, during the solver setup phase, the matrix is first constructed with the computational cost ( ) as detailed in Equation ( 36). The corresponding memory consumption is ( ) = 9 due to the nine point stencil. Since the operator is symmetric, we can then compute its factorization = ΓΣΓ , where Γ is lower triangular and Σ is diagonal, for the known cost

( . ) = O ( 3 /3). ( 46
)
This factorization is also computed only once during setup and the Γ factor is stored in a sparse matrix format. Due to the structure of the matrix , resulting from the use of the nine point finite differences stencil, the factorization fills-in, and the storage of the Γ and Σ factors is ( . ) = O ( 2 /2). These factors replace in memory. At each iteration, simple forward substitution (solve Γ = ) and backward substitution (solve Γ = Σ -1 ) are then applied to solve Equation ( 45) for ( -1 ) = O (2 2 -) flops, which can be neglected compared to the ΓΣΓ factorization above.

Although the computational cost associated with the coarsest grid correction has a non-linear complexity, this is contrasted by the relatively low dimension of the coarsest grid, i.e., ≈ /4 . In order for the cost (46) to become lower than O ( ), we then need the number of grid levels to respect

3 3 < ⇐⇒ > ln 2 3 3 ln(4) . ( 47
)
For example, considering a finest system of size = 10 6 , the required number of levels is then > 6. There is a trade-off to find between a mild cost to solve on the coarsest grid, and a faster convergence.

we can write that ≈ , ≈ ≈ 1/2. Then, as detailed in Appendix 8.4.1, we obtain

Δ ( ℎ ) = -9.75 -2.5 < -9.75, Δ ( ℎ ) = -2.5 -2 < -2.5, Δ ( ℎ ) = -1.5 - < -1.5, (49)
which means that the smoothing procedure is less expensive with the implicit extrapolation. Furthermore, the implicit extrapolation introduces two prolongation operators. Considering that +1 = /4, the operator simply corresponds to an injection operator and thus

( ) = 0 and ( ) = 1, if coarse node, 0, = 1/4. ( 50
)
The operator is very similar to the classical interpolation , except when the neighboring coarse nodes are on the diagonals. In this case (for a quarter of the nodes), we only keep the top-left and bottom-right updates as seen in Table 1. Based on Equation ( 38), these updates cost 8 flops to construct , and 2 flops to apply. We then obtain

( ) = 4, and
( ) = 1.75, ( - ) = ( ) + ( ) = 5.75. ( 51
)
Finally, the computation of the residual on the coarsest grid changes to Equation ( 33), and thus

Δ (X) ( ) = [ ( - ) + 1 + ( - ) + 1/4( ( - ) + 1) + ( - )] -[ ( - ) + 1 + ( - )] ≈ 1/4( - ). (52)
As a summary, the only differences in terms of computational cost from using the extrapolation is mainly a less expensive construction of once during setup (minus at least 10 flops per node; depending on the size of the circle smoother domain), and an overhead of a quarter of an application of per iteration (i.e. the application of from the coarser grid level) to compute the modified residual. This overhead does not change the overall complexity, in particular when compared to what would be obtained from using an explicit higher degree discretization.

Complete and asymptotic complexities

Now, we summarize the complexity of the whole multigrid solver, separating the setup from the actual multigrid cycle. We consider that we apply 1 and 2 pre-and post-smoothing steps and that the convergence is obtained in iterations. With and without extrapolation, we get

( ) = ( ℎ ) -1 =1 + ( ), ( ) = ( [( 1 + 2 ) ( ℎ ) + ( )] -1 =1 + ( )), ( ) = [ ( ) + ( ℎ )] -1 =1 + ( ). (53)
Due to standard coarsening, we have ≈ /4 . As we are only considering orders, not exact equality, we use the equal sign and obtain a geometric series expression

-1 =1 = -1 =1 1 4 = 4 -1 4 -1 3 →∞ → 4 3 . ( 54
)
By taking a large number of grid levels and thus considering that the coarsest grid cost can be neglected, we obtain the linear asymptotic complexities

lim →∞ ( ) = 4 3 ( ℎ ) × , lim →∞ ( ) = 4 3 [( 1 + 2 ) ( ℎ ) + ( )] × , lim →∞ ( ) = 4 3 [ ( ) + ( ℎ )] × ≈ 12 . ( 55
)
The number of iterations was empirically shown to be independent of the problem size in [KKR22], and the complexity of GMGPolar is thus optimal in both computations and memory.

Parallel implementation

In the previous section, we have demonstrated optimal linear complexity for GMGPolar, in both memory and computations, using suited implementations. Also, most of the operations can be vectorized as explained in the previous section. The implementation has laid the base for parallelization as, e.g., the prolongation operator updates all of its entries independently and the matrix operator updates all neighboring nodes by avoiding duplicate computations. In order to make use of the shared memory infrastructure of modern computers, we will detail a multihreading approach for GMGPolar. This parallelization is implemented with the OpenMP-4.5 API.

Application of the intergrid transfer operators

The prolongation operator is the simplest operator to parallelize in our case. Its application to a vector is performed by updating all entries of the output vector independently. In fact, all entries in the output vector correspond to a fine node, and the values are computed as a weighted sum of the neighboring coarse nodes values; see Equation ( 38). However, updating all nodes in parallel can end up with each thread having too little computation to perform, we say that the granularity may be too low. We can then improve the parallelization by grouping the fine nodes of each row (with a fixed radius), which also has the advantage of vectorizing the computation by grouping memory accesses.

As for the restriction operator, the idea is the same but the loop is performed over the coarse nodes instead of the fine nodes, thus updating again the entries of the output vector independently.

Matrix-free application of

Concerning the application of the matrix to a vector , we could implement the same kind of parallelization by updating all entries of the output vector independently, i.e. using theapproach with duplicate computations of the functions , as introduced in Section 5.2.1. Using theapproach, we avoid duplicate computations of the functions , by looping over all nodes to update the neighboring nodes in the nine point stencil. Then, each entry of the output vector, corresponding to a single grid node, is updated by all of its 8 neighbors. In terms of parallelization, this introduces concurrent memory accesses between threads, and means that we need some mechanism to prevent incoherence from such conflicts. We solve this issue by using task-based parallelism and dependencies, respectively introduced in the revisions 3.0 and 4.0 of OpenMP.

The principle of task-based parallelism is to separate parts of a complete code as tasks. These tasks are put in a pool of tasks by the master thread. At runtime, the OpenMP threads can pick and perform tasks as they become available, and so on until all the tasks have been accomplished. Additionally, we can explicitly define dependency relations between tasks in order to avoid conflicts from parallel updates. OpenMP interprets these dependencies by creating a direct acyclic graph (DAG), named dependency graph where each vertex is a task, and the edges define dependencies between these tasks. The threads can only run tasks which dependencies have been resolved before. We obtain an asynchronous mechanism, which is essential to keep a high level of parallelism and efficiency of our application on modern supercomputers.

As before, to prevent issues of granularity and use vectorization, we consider grouping the computations corresponding to one circle line of the grid, i.e., a line with constant radius. We then define as task ( ) the sequential loop over all nodes in the line with radius . In this loop, we update the parts of corresponding to the neighbors of the nodes in line , i.e. all nodes in the lines , -1, and + 1. Since the nodes are considered sequentially inside a row, there can be no conflict of update between two nodes of the same row. However, row is also updated by the rows -1 and + 1, and thus, there can be conflicts from concurrent updates between these three rows. We say that the tasks ( ) , ( -1) , and ( +1) are dependent from each other. Consequently, we have that ( ) is independent from ( +3) .

We can resolve these conflicts in OpenMP through the depend clauses with the following relations for all ∈ {1, •

• • , } ( ) -→ ∅, if ∃ ∈ N 0 s.t. = 3 + 1, ( -1) , ( +2) , else. (56)
Through the dependency mechanism in OpenMP, these relations define a directed graph which automatically prevents conflicts when applying the matrix . Taking the example of a grid with nine divisions in the direction, Figure 7 illustrates the relation (56), i.e. the dependency graph between tasks. In summary, the tasks (3 +1) are can be run right from the start because they are independent from all other tasks. The execution continues with the tasks (3 +2) , then with the tasks (3 +3) . We then have a theoretical maximum speedup of around /3 which ensures a good scaling for large problems.

Smoothing

As in the Section 5.5, we define the first line for the circle and radial partial smoother as white.

Given the change from circle to radial smoothing of (26), the color of the outermost circle smoother line will be defined. However, for large grids, the difference in smoothing between two 1 2 3 4 5 6 7 8 9

Fig. 7: Dependency graph for the application of . Each circle/square is a task, i.e., the handling of one line of nodes. Tasks (3 +1) , = 0, 1, 2, . . . are within a square and are independent from all other tasks. An arrow from task ( ) to task ( ) , indicates that task ( ) must wait for task ( ) to finish.

lines is minimal -so, to simplify the notation, we consider that the outermost line of the circle partial smoother is white too. Based on the parallelism implemented for the application of , we can also parallelize the construction of the matrices during the setup. Here, we focus on the parallelization of the complete smoother, i.e., applying all four of the partial smoothers ( ∈ { , }, ∈ { , }) in a block Gauss-Seidel approach.

Let's remind Equation (30), i.e.

( ) ( ) = ( ) -( ) ⊥ ( ) ⊥ ,
showing that the smoothing procedure implies for each line the application of the operator ( ) ⊥ and the solution of the system itself. For these, we define the tasks ( ) ⊥ and ( ) , respectively. As before, in the case of the circle partial smoother, a line is a circle with constant radius, and for the radial partial smoother, a line consists of the nodes with a constant angle. For reasons outlined below, we number the lines in increasing order of the angle for the radial partial smoother, i.e. indices = 1, . . . , θ correspond to angles θ 1 , . . . , θ θ , but number the lines in decreasing order for the circle partial smoother, i.e. indices = 1, . . . , correspond to radii , . . . , 1 . As detailed in Section 5.3, to use an A-give approach for the application of the ⊥ matrices in the smoother, it is required to loop twice over all nodes, once per color of each smoother. Also, the first circle in the domain corresponding to the radial partial smoother, i.e., the line with radius +1 , must be handled in the circle partial smoother in order to fully update the line. This line is classically called a halo, and linked to the task (0) ⊥ . Thus, there are + 1 tasks of type ( ) ⊥ , tasks of type ( ) ⊥ , and finally θ of type ( ) ⊥ (resp. ( ) ⊥ ). However, the application of the specific ( ( ) ) -1 matrix is only required for the line and not for the neighboring ones. Thus, there are approximately /2 and θ /2 tasks of type ( ) for the circle and radial partial smoothers. For these tasks, we keep the numbering of the complete circle (or radial) partial smoothers, e.g. task ( ) corresponds to the -th line of the radial partial smoother which is a line with constant angle θ .

At each relaxation iteration, it is then possible to run in parallel the application of all tasks, for all smoothing types, as long as we respect the following dependencies:

1. Similarly to the application of the matrix , there are conflicts from concurrent updates between lines of the same type of smoother. These conflicts can be prevented using dependencies between tasks with the rules

( ) ⊥ -→ ∅, if ∃ s.t. = 3 + 1, ( -1)
⊥ , ( +2) ⊥ , else. Tasks:

tasks

(3 +1) ⊥ , ( 3 +2)
⊥ , ( 3 +3) ⊥ tasks ( ) * halo Dependencies: 2 3 4
Fig. 8: Dependency graph for the application of the smoother. Each node is a task, i.e., the handling of a circle line of nodes. The tasks ( ) ⊥ , = 2, 3, 5 and = 2, 3, 5, 6, are represented by circles and tasks ( ) ⊥ , = 1, 4 are squares. Tasks ( ) , = 1, 2, . . ., are represented by triangles. The different type of dependencies are colored: 2=blue, 3=green, and 4=red. For the sake of a simple presentation, we do not display the dependencies from (57), as they are identical to what is shown in Figure 7.

For the coefficients and , we set The solver GMGPolar was initially developed for the Equation (1) with = 0. Here, we use a more general form to show that GMGPolar naturally extends to a nontrivial coefficient, exhibiting similar properties of convergence and accuracy as in the simpler case.

( ) = exp       -tanh ( , ) -       , (62)
( )= exp       tanh ( , ) -       , (63)
To handle the singularity introduced by the polar geometry at 1 = 0, GMGPolar has demonstrated in [KKR22] the possibility through special considerations of the finite difference discretizations. However, we simplify these considerations by simply setting a minimum radius Our experiments are performed using an OpenMP parallelization on a single node with 2sockets having 256 GB memory. Each socket is an AMD EPYC 7702 "Rome" processor at 2GHz containing 64 cores, constituting a separate NUMA (non-uniform memory access) domain.

Here, we fix the CPU frequency to 1.8 Ghz and only use one thread per core. We use LIKWID 5.2.2 [HWT10, GEHW22] for measuring the amount of computations in flops, and the memory consumption in Mbytes. Finally, to solve the problem on the coarsest grid of GMGPolar, we use MUMPS v5.4.1 [AAB + 15].

In this section, we show the different types of scalability for GMGPolar from numerical through weak and finally strong scaling.

Numerical scalability of GMGPolar

In this section, we first demonstrate numerical scalability for the model problem and the geometries considered. Note thatn numerical scalability of GMGPolar for other configurations and smaller model problems has already been shown in [KKR22].

Here, we consider the right hand side given by the manufactured solution (61), called the Cartesian solution. The problem is defined on the Shafranov geometry, see Figure 9a, and the Czarny geometry, see Figure 9b. We construct grids of increasing refinement to get problems of size from 8 5 to 5 7 DOFs.

Tables 4 and5 provide the errors in (weighted) ℓ 2 and ∞ norm with respect to the manufactured solution, and provide the error reduction order as ord. for both norms. The order is calculated as the error reduction from one row to its predecessor, i.e., ord = log err err -1 /log gridsize -1 gridsize .

For both geometries, we obtain an approximation order around 3.5 (which is roughly below 4, which we obtained for the polar use case) in ℓ 2 -norm as expected from the use of the implicit extrapolation scheme. In fact, we obtain similar approximation orders as in [KKR21] and [KKR22] but solving here the equation ( 1) with a non-zero coefficient and also the Czarny geometry. Note that the number of iterations in all cases stays relatively stable with respect to the problem size. In fact, the number of iterations even tends to decrease for larger problems, e.g. from 56 to 49 iterations from the smallest to the largest problems on the Czarny geometry.

Weak scalability of GMGPolar

We now provide the weak scaling experiment for the exact same problems of the previous sections, starting from 5 7 degrees of freedom (DOFs). We demonstrate the overall parallel efficiency of GMGPolar in shared memory parallelism.

In Figure 11 we give the amount of computation in flops for the multigrid iterations of GMGPolar depending on the problem size, as measured using likwid. We include the asymptotic theoretical number of flops. The theoretical asymptotic computational complexity is computed based on Eq. ( 55) with the actual number of iterations, completed by the results in Table 3 and the specifities of the implicit extrapolation detailed in Section 5.5.

As expected, we observe that the complexity of GMGPolar is increasing linearly with the size of the problem. Proving this is one of the main results from this article. The computational cost measured in practice is twice as large as the theoretical one, which leaves space for optimization. However, note that the expected theoretical numbers are asymptotic and that lower order computations are not considered in the theoretical sum. Finally, we perform a weak scalability experiment with the parallel application of GMGPolar on these same problems of size up to 5 7 DOFs, from 1 to 64 OpenMP threads. Here, each thread handles a constant number of 5 7 DOFs, and the ideal behavior is obtained when the execution time stays identical for all problem sizes. The setup phase of GMGPolar is (currently) sequential, thus we focus again on the parallel execution of the multigrid iterations. Figure 12 gives the parallel efficiency of the multigrid cycle phase depending on the problem size. We observe an excellent efficiency from 1 to 16 with 71% and 84% respectively for the Shafranov and Czarny geometries. This parallel efficiency is obtained thanks to the efficient parallelization scheme proposed in Section 6, using task-based parallelism. Using the full socket of 64 cores still provides additional gain overall, with an efficiency staying around 50% in the case of Czarny, although there is also potential for future improvement. In fact, the total computation time only doubles while the problem size is multiplied by 64. Note that GMGPolar seems slightly more scalable in the case of the Czarny geometry, it is likely that with the more expensive and repetitive computation of the mapping of the geometry in the case of Czarny, the granularity of the problem increases and this in turn improves the parallel efficiency of the implementation. In the next section, we will detail the scalability of each parts of the multigrid cycle to identify the ones currently limiting our parallel implementation.

Strong scalability of GMGPolar

In this section, we focus on a finer discretizations of Eq. ( 1) for which we consider again the two different geometries and the Cartesian solution (61). First, we perform a strong scaling experiment by applying GMGPolar with 1 up to 64 OpenMP threads, on a grid of size 6145 × 8192, i.e. a grid with about 50 million DOFs. Figure 13 presents the resulting execution times for the multigrid cycle phase respectively for the Shafranov and Czarny geometries. Again, the setup phase is excluded since it is sequential in our implementation. Additionally to the complete multigrid iterations, we detail the efficiency of several inner steps: the smoothing, the computation of the residual, the computation of the finest residual, the restriction and prolongation, the solution of the coarse problem.

From Figure 13, we observe a very good overall scaling from 1 to 16 cores, where we get a good speed-up of around 13. Although using the full socket of 64 cores still decreases the execution time, the scalability decreases and we see there is potential for future improvement for the other phases. The most efficient step is the smoothing which is often viewed as the critical part in a multigrid method. In fact, up to 32 threads, the smoothing procedure shows a great scalability. Then, the performance of the smoothing decreases but we still get a speed-up of around 35 when using 64 threads. This good scalability comes from the implemented task-based parallelism, making use of the decomposition in lines of fixed radius and the limited interactions between these.

Most of the other operations in GMGPolar have a minor effect on the total run time. For example, the matrix-free application of the residual and prolongation operators do not scale but have a low impact on the total execution time. Only the computation of the finest residual significantly worsens the impact scalability and needs to be improved for future applications. In fact, with 1 thread, computing the residual amounts to around 10% of the total execution time, and increases to 30% with 64 threads.

Conclusion

In this paper, we have presented a matrix-free approach of GMGPolar which uses a nine-point stencil inside a taylored, implicitly extrapolated multigrid scheme with zebra line smoothing adapted to different parts of the domain. We have newly presented a rigorous analysis of flops and memory usage and demonstrated optimal, linear asymptotic, complexity of the considered solver in both computations and memory. Although the solution of the coarsest grid problem does not have a linear complexity, considering a sufficiently large number of grid levels, it can be neglected for the overall cost analysis of the GMGPolar solver. Our rigorous analysis is kept as general as possible such that its framework ideas can be transferred to matrix-free approaches of other stencil-based operators.

Furthermore, we have presented a task-based multithreaded parallelism for GMGPolar which takes advantage of the disk-shaped geometry of the problem and defines dependencies accordingly. We have demonstrated numerical, weak, and strong scalability for use cases from fusion plasma applications. GMGPolar can now sufficiently handle disk-like geometries with tens and hundreds of millions of degrees of freedom. Hereby, the implemented node-level parallelism scales very well from 1 to a 16, 32, or 64-fold of cores. Together with the results of [BLK + 23], the solver is ready to be coupled with gyrokinetic codes such as Gysela [GAB + 16].

Future developments of GMGPolar will involve further optimizations of the run-time and the use of accelerators for computation.

Algorithm 4 ℎ : (8 mult., 5 add, 5 div. = 18 flops) Input: ( , θ), (cos θ, sin θ), ( , ), . 1:

= (-2 cos θ + cos θ)/ (3 mult., 2 add, 2 div.) 2: θ = ( sin θsin θ) (2 mult., 1 add, 1 div.) 3: θ = ( + 1) sin θ/

(1 mult., 1 add, 1 div.) 4: θθ = ( cos θ + cos θ) (2 mult., 1 add, 1 div.)

Algorithm 5 : (52 mult., 25 add, 23 div., 12 sqrt = 112 flops)

Input: ( , θ), (cos θ, sin θ), ( , ), .

1:

= -cos θ ( +2 cos θ) +1 / ( 3
mult., 2 add, 3 div., 1 sqrt)

2: θ = sin θ ( +2 cos θ) +1
/ (4 mult., 2 add, 3 div., 1 sqrt)

3: θ = sin θ cos θ √ 1-(1/4) 2 2-( +2 cos θ) +1 2 ( +2 cos θ) +1 + sin θ √ 1-(1/4) 2 2-( +2 cos θ) +1
/ (22 mult., 11 add, 9 div., 5 sqrt)

4: θθ = - sin 2 θ √ 1-(1/4) 2 2-( +2 cos θ) +1 2 ( +2 cos θ) +1 + cos θ √ 1-(1/4) 2 2-( +2 cos θ) +1 / ( 23
mult., 11 add, 9 div., 5 sqrt)

functions

Using Equation ( 10), the functions are defined as

1 2 θ 1 2 θ θθ = 1 2 -1 -| det( )| ⇐⇒            = 1 2 2 θθ + 2 θ /| det( )|, θθ = 1 2 2 θ + 2 /| det( )|, θ = ( θθ θ + θ ) /| det( )|, (67)
with det( ) = θθθ θ . Since the functions are built from the same elements, it is more efficient to compute all three functions simultaneously with the cost ( ) = 19 + ( ), as detailed in Algo. 6. In order to build/apply or the smoother matrices and ⊥ , different combinations of the functions must be computed with a corresponding computational cost. In fact, considering Algo. 6, there is a common cost of 6 flops for all three functions (lines 1 to 4) plus 4 flops, 4 flops, and 5 flops, respectively to compute , θθ , and θ .

Matrix-free application of

We now detail the computation of the complexity for the construction and the matrix-free application of the matrix . We recall that the matrix is based on the nine point finite-difference discretization introduced in Section 3 with the entries from Equations ( 19)- [START_REF] Hoelzl | The JOREK nonlinear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF].

In this section, we establish some tools in order to compute the computational cost per node required for the construction of a matrix in the context of discretized PDEs. We believe that our method can be generalized to most approaches (using e.g. other finite difference stencils or finite Algorithm 6 computation (DF, 13 mult., 4 add, 1 div., 1 abs. = 19 + ( ))

1: Compute = θ θ θθ ( (
= 2coeff × ( θθ θ + θ ) (4 mult., 1 add) = 6 + ( ) = 4 = 4 = 5
elements) in order to get a complete knowledge of how expensive constructing the discretized operator is. As detailed below, our method is based on decomposing the computation of a stencil in multiple contributions, and visualizing the corresponding costs as a diagram.

Construction of A:

Upon implementation, one should take care of precomputing as much variables as possible, and computing only once factors which are used several times. How this can be done depends on the code. The general idea is to loop over all nodes and compute the entries of the matrix corresponding to each node independently. As discussed in Section 5.2.1, we can consider two approaches. A-take: In a first approach, the loop is performed over each node ( , ) for which all entries in the corresponding row are computed. In this approach, we precompute all the functions required for the stencil, i.e. all three functions ( , θθ , θ ) for the current node and its neighbors (top, left, bottom, right) except θθ for the top and bottom nodes, and for the left and right nodes. Based on Algo. 6, we get that ( ) = 5 ( ) + 79 flops.

Using the notation above, we then consider that each appearance of the functions or the coefficient in the stencil (68) defines a factor, resulting in 25 factors. The factors are numbered in their order of appearance for each update in the stencil. It is easy for each of these factors to get the corresponding number of flops. In fact, all of the 8 diagonal factors require 1 flop each and appear only once, e.g. ⊗ appears in the top and central updates. Finally, there is a single factor corresponding to the coefficient which requires 7 flops.

We give a representation of this approach in the first line of Table 6. We draw one diagram, called complexity, displaying the computational complexity corresponding to all factors, under the form (#flops + #occurences) × #factors. For example, on the top place of the stencil, we write the complexity (4 + 2) × 2 corresponding to each factor (1) ⊗↑ and (2) ⊗↑ , both needed in the top and central update. A second diagram, called update, specifies for each complexity, which factors of the stencil are computed. By summing all of these complexities, we then obtain the final complexity per node A-give: Alternatively, the A-give approach ensures that each function is only computed once. In this approach, the loop is performed over each node, where the local functions are computed. The corresponding complexity is ( ) = ( ) + 19. Then, all the factors requiring these functions are computed to update the matrix. These factors are part of the matrix entries in the row and column corresponding to the current node.

It is then possible, as before, to simply compute factors for the stencil of the current node, e.g. (1) ⊗↑ = -( + -1 ) 2ℎ ⊗ is involved in the top and central updates. However, this same factor also appears in the bottom and, again, in the central update of the neighboring node placed above in the grid, i.e. in ↑⊗ and ↑↑ . The corresponding complexity for this factor is then (4 + 4) × 1. Compared to the A-take approach, less factors are computed per node but their number of occurrence is increased.

We give again the representation of this approach in Table 6. In the update column, we additionally indicate with a " * " symbol which factors are computed to update the stencil of a neighboring grid node, with an arrow showing which neighbor. Factors without a " * " symbol simply correspond to updates for the current node. Note that the factors for the diagonal updates of the stencil require only the functions from the neighboring nodes. For example, the topright update depends on the functions θ of the top and right neighboring nodes. Thus, there never is a purple dot in the corresponding representation. By combinating all the factors, we then obtain the final complexity per node The computational cost for the construction of is thus much lower when using the A-give approach. This is especially true for complex geometries such as Czarny (see Algo. 5), where computing the Jacobian of the mapping would become the dominating part when building the matrix.

In terms of memory, the matrix contains nine entries per row and thus we have ( ) = 9. This cost only appears on the coarsest grid level where a direct solver is used to solve a linear system based on . On the other levels, is applied with a matrix-free approach.

Table 7: Stencil representations for the construction of and ⊥ for the circle white partial smoother ( , ) in the case without extrapolation. This representation follows the same principles as Table 6. We distinguish two cases for each matrix: when the current node is of the same color as the partial smoother (white), or when it is of the other color (black). Updates of the local stencil are only present when considering a node of the same color as the partial smoother. Note that all other partial smoothers behave in the same way. (1) ↑→

(1) ↑⊗ ( 1) ↑←

(1) ↓→ ( 2) ↓⊗

(1) ↓←

Finally, in order to apply ⊥ in a matrix-free fashion, as detailed in Section 5.1, we only need the additional application cost ( ⊥ ) = 12 which corresponds to the number of Table 8: Stencil representations for the construction of and ⊥ for white nodes when using the circle white partial smoother ( , ), in the case where implicit extrapolation is used. This representation follows the same principles as Table 7. Additionally, we distinguish two cases for each matrix: when the current node is fine or when it is coarse. (2) →↑ which can be done with the combination of an LU factorization, using a modified Thomas' algorithm, and the corresponding backward/forward substitution. Starting from the particular sparsity pattern of ( ) (see Section 5), we consider a separate line of size .

Fig. 1 :

 1 Fig. 1: A circular geometry described (Left) in Cartesian coordinates and (Right) in polar coordinates ( , θ) ∈ [ , ] × [0, 2 ) with mappings and -1 between the different descriptions.

  Fig. 2: Two geometries are considered.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Left: Schematic subdivision of a circular domain into an interior subdomain where circleline smoothing is used and an exterior subdomain where radial smoothing is used. Right: Circular mesh with 64 nodes. Circular smoothing applied on the four interior circles (dark gray background color), radial smoothing applied on the four exterior circles (light gray background color). Nodes colored according to smoother color.

Fig. 5 :

 5 Fig. 5: Schematic representation of the two approaches to apply the matrix : take or give the evaluations of the functions evaluations from or to the neighboring nodes, respectively.

Fig. 6 :

 6 Fig. 6: Sparsity pattern of the matrix for a line of eight nodes.

where = 0

 0 .05 and = 0.7, as in [BLK + 23], and = 1.3 by our own choice as in [KKR22]. Note that was set to 1 in [BLK + 23]. See Figure 10 for an illustration of the coefficients and .

Fig. 9 :

 9 Fig. 9: Illustration of the manufactured solution (61) for the Shafranov and Czarny geometries.

Fig. 10 :

 10 Fig. 10: Coefficients (62) and (63) used in Equation (1).

Fig. 11 :

 11 Fig.11: Computational cost in flops when applying GMGPolar for the manufactured Cartesian solution (61) on the geometries Shafranov and Czarny. We include the number of flops as measured by LIKWID as well as the theoretical (asymptotic) expectations.

Fig. 12 :

 12 Fig. 12: Weak scaling of GMGPolar. Parallel efficiency compared to the sequential run from 1 to 64 cores on grids of size 769 × 1024 to 6145 × 8192, manufactured Cartesian solution (61).

Fig. 13 :

 13 Fig. 13: Strong scaling of GMGPolar. From 1 to 64 cores on a grid of size 6145 × 8192, Shafranov geometry Figure 2a (left) and Czarny geometry Figure 2b (right) and manufactured solution (61).

  appears only in ⊗ . The top, bottom, left, and right factors require 4 flops each and appear twice, e.g. (1) ⊗↑ = -( + -1 ) 2ℎ

Table 4 :

 4 Numerical scalability of GMGPolar from 787k to 50m degrees of freedom on the Shafranov geometry and Cartesian solution (61).

	×	its	ℓ 2	ord.	∞	ord.
	769 × 1 024	70	3.18e-08	-	9.40e-07	-
	1 537 × 2 048	67	2.86e-09 3.48 1.18e-07 3.00
	3 073 × 4 096	64	2.55e-10 3.49 1.47e-08 3.00
	6 145 × 8 192	61	2.26e-11 3.49 1.84e-09 3.00

Table 5 :

 5 Numerical scalability of GMGPolar from 787k to 50m degrees of freedom on the Czarny geometry and the Cartesian solution (61).

	×	its	ℓ 2	ord.	∞	ord.
	769 × 1 024	56	2.31e-08	-	7.05e-07	-
	1 537 × 2 048	54	2.01e-09 3.51 8.81e-08 3.00
	3 073 × 4 096	51	1.79e-10 3.50 1.10e-08 3.00
	6 145 × 8 192	49	1.58e-11 3.50 1.38e-09 3.00

Table 6 :

 6 Stencil representations of the construction of in terms of computations for a single grid node ( , ). Two implementations are considered: A-take and A-give. Functions to precompute for the current node are given. The complexity column gives for each point in the stencil what is the computational cost expressed as (#flops + #occurrences)×#factors, where #flops is the complexity of computing the scaling coefficient of each factor, #occurrences is the number of times this factor appears in the stencil updates (= 1 addition in the matrix entry), and #factors is the number of different factors to compute with the same complexity. The colored lines indicate which function is required for each computation: green, blue, and red respectively for , θθ , and θ . The update column gives the factors in the stencil corresponding to each complexity. The " * " symbols indicate which factors are computed for the stencil of a neighboring grid node, indicated by the arrow direction. Other factors are updated in the stencil of the current node. Purple dots show the updated stencil points, and a rectangle stands for the central update.

https://www.iter.org/fr/sci/beyonditer

= 1 -8 and constructing Dirichlet boundary conditions on this innermost circle of the domain from the manufactured solution. When considering the asymptotic computational and memory complexities, we did not consider the boundary conditions, so this choice for 1 does not affect this study.In our experiments, the iterations are stopped based on a relative reduction threshold of 1 -11 for the residual, or after a maximum of 300 iterations of GMGPolar. At each level of the multigrid cycle, we perform one pre-and one post-smoothing step. We always use implicit extrapolation and the matrix-free implementation of GMGPolar.

https://www.sympy.org/en/index.html

Availability of data and material

No particular data was used.

Code availability

The code is available on GitHub https://github.com/SciCompMod/GMGPolar and with tagged version 1.0.2 and persistent identifier on zenodo https://doi.org/10.5281/zenodo.1040929. Grids have been created with version 1.0.0, simulations have been performed with version 1.0.2.

Implicit extrapolation

The implicit extrapolation technique allows to increase the order of approximation of the solution computed in GMGPolar without the explicit use of a higher order discretization. We now detail all changes in complexity involved in the use of the implicit extrapolation, and show in particular that the complexity stays mostly unchanged. We denote by (X) and (V) the computational cost of operation X and memory requirement of data structure V when using the extrapolation.

Note again that the extrapolation only acts between the two finest levels and that the smoother on the finest level is only applied on the fine nodes. Let's remind the reader that the coarse grids are defined using standard coarsening as detailed in Section 4.1. Thus, starting from the first grid node, the coarse nodes are chosen with a distance of two from each other. Assuming without loss of generality that the first circle line and the first radial line are white, half of their points are coarse nodes, i.e. kept on the next coarser grid level. As defined in (32), these coarse nodes are kept out of the white smoothers. Written in terms of a linear equation, instead of solving (30), we have ( )

, ⊥ ( )

for each line of a white-colored smoother. Here, , ⊥ is the identity matrix acting on the coarse nodes and likewise we have denoted subvectors and matrices of and . Consequently, we only solve the second line of the previous equation set. In order to simplify the notation, we will in the following generically use ( ) for the nontrivial smoothing matrix ( ) , , .

-Using implicit extrapolation, the ( ) matrices are diagonal. In fact, the ( ) matrices correspond to the submatrix of restricted to the white fine nodes of the line as also shown in the previous equation. These fine nodes are independent from each other when using the nine-point stencil, i.e. a stencil of length one. Consequently, only the central point of the stencil is needed for each node of this type, i.e. the diagonal entry in the matrix, see Table 2. Thus, the smoother matrices change from tridiagonal to diagonal. Then, setting aside the cost to compute the Jacobian , the resulting computational cost represents around 2/3 of the cost without extrapolation. The matrices ( ( ) ) -1 are then trivial to apply, thus ( -1 ) = . We only need to keep the diagonals in memory, i.e. ( ) = .

-For the coarse white nodes corresponding to each partial smoother, the ( ) ⊥ matrices now apply the nine point stencil except for the central update. In fact, ( ) ⊥ contains around /2 additional columns corresponding to the coarse white nodes. The result is an increased cost per line, due to the additional nodes where all functions must be computed, and an increased number of corresponding additional updates. However, since the partial smoother is only applied on the fine nodes, ( ) ⊥ contains around /2 less rows when using extrapolation. As a result, ( -⊥ ) stays unchanged compared to the case without extrapolation, as detailed in Appendix 8.4.1.

As for the black points, the extrapolation does not change the partial smoothers. Let's introduce Δ (X) = (X) -(X), and Δ (V) = (V) -(V), the overheads from using the extrapolation scheme. Since the number of black and white lines only differs by one at the maximum, where ∈ {1, • • • , } for the circle smoother, and ∈ {1, • • • , θ } for the radial smoother. 2. Naturally, ( ( ) ) -1 is applied on ( ) -( ) ⊥ ( ) ⊥ , and thus depends on the update of ( ) ⊥ for this line and the two neighboring lines. Thus, we impose the rule

where ∈ {1, • • • , } for the circle smoother, and ∈ {1, • • • , θ } for the radial smoother. 3. Then, the application of a black smoother for a line must wait for the application of the corresponding white smoother on the neighboring lines, i.e.

( )

, ( +1) , else.

(

where ∈ {1, • • • , } for the circle smoother, and ∈ {1, • • • , θ } for the radial smoother. 4. Finally, there are also conflicts from the updates between the circle and radial smoothers. In fact, to start the smoothing of radial lines, the outermost line of the circle smoother must be finished. Thus, considering that the line is white and corresponds to task (1) , we impose the additional dependency

In order to start the radial partial smoother as soon as possible, in parallel with the circle partial smoother, we need to finish task (1) first which is dependent on (1) ⊥ , (2) ⊥ and

⊥ . Thus, we choose to start the relaxation by sequentially performing the three latter tasks, in parallel with all independent tasks (3 +1) ⊥ , which explains the reverse numbering of tasks for the circle partial smoother.

These dependencies are illustrated in Figure 8 on an example with five circle lines, and six radial lines. Note that we do not include the dependencies corresponding to the first rule above, since they follow the same rule as for the application of A but in a descending order, see Figure 7.

By accumulating all of these rules through the OpenMP dependencies, we obtain a directed graph of tasks which prevents any conflicts from happening. Similarly to the application of , we then have an asymptotic maximum speed-up of /3 + θ /3 , also ensuring a good scaling for large problem sizes.

Numerical experiments

In the following, we will present numerical results for the model problem (1). The coefficients and used here as well as the manufactured solution are obtained from the use cases in [Zon19], [KKR22] and, precisely, [BLK + 23, Eq. ( 18)-( 19) and ( 23)]. These cases are inspired by the simulation of plasma in tokamak fusion reactors. Here, we consider a Cartesian solution with oscillations aligned with the Cartesian grid:

Figure 9 illustrate the solution for the Shafranov and Czarny geometries.

Declarations

Funding and Acknowledgements

We acknowledge the computing time granted on the HPC cluster CARO. Furthermore, we would like to thank Nils Kanning for support on CARO and Melven Röhrig-Zöllner and Johannes Wendler for fruitful hints and discussions on LIKWID and performance monitoring issues.

Conflicts of interest/Competing interests

None.

Authors' contributions MJK, CK, and UR developed the numerical algorithm, the implementation was done by PL, CS, MJK and CK. MJK and CK did the initial implementation in Matlab, PL and CS ported the code to C/C++ and did the task-based parallelization. Simulations were conducted by MJK. All authors wrote and revised the script.

Appendix A: Detailed computational complexity and analysis framework

In this appendix, we detail how to compute the complexity of the components constituting GMGPolar. The main steps and results of this appendix are given in Section 5, and we use the same notation as introduced in Section 5.1. We remind the reader that cos θ, sin θ, , and are precomputed with a negligible cost compared to a linear complexity. While Appendix 8.1 provides exact number of flops for the geometries considered, Appendix 8.2 is kept general for applications with mappings on curvilinear coordinates.

Jacobian of the mapping

Note that the algorithms we provide to compute the Jacobian are far from optimal. To stay general with respect to the application, the corresponding codes are generated from the analytical formula of the mapping function using the python library SymPy 2 for symbolic mathematics. Appendix 8.3 is written in general framework form so that it can be reused for other stencilbased operators and Appendix 8.4 is applicable for circle-radial line smoothers (also providing the extrapolation case).

In order to compute the entries of the matrix using the nine point stencil as introduced in Section 3, the Jacobian of the geometry mapping

is required. This Jacobian naturally depends on the choice of the geometry. Based on Section 2, we can use:

1. The polar geometry with the Jacobian

which can be computed with ( ) = 6 flops based on Algo. 3.

Algorithm 3 : (2 mult., 4 div. = 6 flops)

Input: ( , θ), (cos θ, sin θ), . 1:

= cos θ/ (1 div.) 2: θ =sin θ/

(1 mult., 1 div.) 3: θ = sin θ/

(1 div.) 4: θθ = cos θ/

(1 mult., 1 div.)

2. The Shafranov geometry (3) with the Jacobian

which can be computed with ( ℎ ) = 18 flops based on Algo. 4. 3. The Czarny geometry (4) with the Jacobian

with = ( + 2 cos θ) + 1 and = 1/ 1 -2 /4, which can be computed with ( ) = 112 flops based on Algo. 5.

Construction of a generic matrix

Let us consider a square matrix of size such that its entries can be expressed as the sum of a few contributions, i.e. = ∈ ( ) , where depends on the entry ( , ). Note that the contributions are not uniquely defined and may vary depending on what shall be precomputed and reused.

If we consider that the matrix is built from a discretization using a specific stencil, then the non-zero entries of the matrix are defined by the form of the stencil, and the sets of contributions can defined by the expression of each update in the stencil. To explain our general approach, let us now consider the example of the nine-point stencil ( 19)-( 21) used in GMGPolar and a node with coordinate indices ( , ).

In this example, we directly use the coordinate indices ( , ) of each node instead of or and we use instead of to distinguish the example from the general consideration. We assume to be ordered with the coordinate indices. This way, the matrix is defined up to a certain renumbering which can be chosen freely upon implementation. Generally, we use = * θ + , i.e. numbering nodes following the circles lines. Finally, to make the notation more visual, we replace the coordinates ( , ) by ⊗, and coordinates ( ± 1, ± 1) by an arrow indicating the direction of the corresponding node with respect to the node ( , ) (e.g.

for the coordinates ( -1, + 1)). We have the expression of each update

Continuing the above example, we get the set of contributions ⊗↑ = {1, 2} with

From Equation ( 68), we understand that some elements which are common to different updates should be precomputed to save flops. In our example, the functions are to be computed only once per node ( , ) while the sums and quotients of ℎ and are cheap to compute. In a general sense, we further write

where we have certain variables = ( ) ∈ {1,..., }, ∈ to precompute and corresponding scalars ( ) ∈ {1,..., }, ∈ that will be recomputed each time. The distinction between and the is given by the individual costs of each value. In order to keep the memory usage small, should only contain a limited number of costly evaluations. Multiple appearances of the same value in are to be understood as a single evaluation. The computational complexity to build a matrix , can then be expressed as

where ( ) and ( ( ) ) are the computational complexity respectively for the precomputation of the variables and the coefficients ( ) (including one flop for the multiplication ( )

) plus one flop corresponding to summing the factor ( ) into the matrix entry .

Additionally there can be factors ( ) which appear multiple times, i.e. ∃( 1 , 1 , 1 ), ( 2 , 2 , 2 )

2 2 . An obvious example from GMGPolar is the fact that the central update of the nine-point stencil is the sum of the left, right, bottom, and top updates, as seen in Equation (68). We can then write:

where defines the set of factors ( ) such that factors ( ) were factorized, and is the number of occurrences of the factor ( ) , i.e., which can be reused -1 times. Finally, if there are factors which are numerically different but correspond to the same computational cost and number of occurences, we can write

where was reduced further to factorize factors with the same cost, and ( ) gives the number of factors considered to have the same cost and number of occurences as ( ) for the final calculus of the computational complexity. To simplify, setting aside the precomputation of the elements, for each aggregated factor in the cost is decomposed as

Matrix-free application (with A-give approach):

As we advocate to use the favorable A-give approach, we only present the final result for Agive. The result for A-take can be computed easily. As detailed in Section 5.1, the computational complexity from applying a matrix in a matrix-free fashion can be decomposed in its construction cost ( ), and the cost for its application ( ). The latter corresponds to the number of factors in the stencil. With (75), we thus have We now detail the complexity from the construction of the matrices and the matrix-free application of the ⊥ . For the same reasons as for the matrix in section 8.3, we follow an A-give type of approach. We focus on the case without extrapolation first, then detail the differences when using extrapolation.

Smoothing without extrapolation:

When constructing the operators, either or ⊥ , we need to loop over all nodes connected to the nodes of the current partial smoother, i.e., the nodes with the same color (from ), and the nodes with the other color (from ⊥ ). For example, when handling the (circle white) partial smoother, we also need to loop over the nodes corresponding to the ⊥ = partial smoother since these nodes are linked through the nine point stencil. Also, in order to update the nodes of the outermost circle line of the circle partial smoother, we need to loop over the nodes in the first circle line of the radial partial smoother. This line is classically called a halo. This is explained in Section 5.3 in the main body of the article.

To get the corresponding computational complexities, we use the same approach and notation as for the matrix . We thus detail the complexity stencils in Table 7 in the case of the circle white smoother ( ), distinguishing when the current node is from the current partial smoother, or has the other color. The complete complexity for the smoother is then obtained by summing these two cases.

Remark 6 Note that if we merge all the stencils in Table 7 ( and ⊥ , for a node and a node ⊥ ), we reconstruct the original 9-point stencil for , see Table 6.

The complexities per node for the radial partial smoother are obtained with a 90°rotation, with identical computational cost per node. Also, there is no difference in computational cost per node depending on which partial smoother is applied. Thus, all partial smoothers behave in the same way. Using Table 7 similarly as for the matrix , we obtain the final complexities per node

factors to compute ⊥ . The latter number can be obtained by counting the factors from both and ⊥ cases in Table 6. Finally, we get

Smoothing with implicit extrapolation: When using the implicit extrapolation, the smoother is applied only on the fine nodes. Let's consider without loss of generality that the first circle line of the grid ( = 1 ) is white. Let's also remind that the coarse grids are defined using standard coarsening as detailed in Section 4.1. Thus, starting from the first grid node ( 1 , θ 1 ), the coarse nodes are chosen with a distance of 2 from each other. With respect to and ⊥ , nothing changes in the partial smoothers for black nodes when using extrapolation because they are all fine nodes. However, concerning the white nodes, we distinguish the two cases, whether the current node in the loop is a fine node or a coarse node.

We provide Table 8, again representing the complexity per node as well as the computed factors for the construction of and ⊥ for white nodes of the circle partial smoother. This table does not include the case of black nodes, which can still be found in Table 7 (case without extrapolation). As before, the complexities for the radial smoother are identical to those of the circle smoother.

Let us remember that the ratio of nodes for each color is expressed as = / ≈ 1/2, where is the total number of grid points, and the number of grid points of color ∈ { , }. Then, the ratio of white fine and coarse nodes are = = /2.

Remark 7 The number of black lines is identical to the number of white lines plus or minus 1.

Asymptotically, the previous approximations can be considered exact.

The complexity per node for black nodes is unchanged compared to the one without extrapolation. Following the same principles as before, using

Thus, when activating the extrapolation, the cost of applying ⊥ stays the same, and the cost from building is decreased.

LU factorization for the partial smoothing

To apply the partial smoothers, it is necessary to solve the equation

-Circle smoother: tridiagonal matrix with periodicity conditions on the first and last rows. Using Algo. 7, it is possible to compute an LU factorization of the matrix with 10 flops.

Using this factorization, it is possible to solve Equation (80) with 9 flops using Algo. 8. Due to the last row and column filling-in, the resulting memory consumption is (LU( )) = 5. In order to save some memory footprint, the LU decomposition is written directly inside the matrix. L is in the lower triangular part of the modified ( ) with an omitted diagonal of ones, U is the diagonal and upper triangular part of the modified ( ) .

-Radial smoother: tridiagonal matrix. Using Algo. 9, it is possible to compute an LU factorization of the matrix with 3 flops. Using this factorization, it is possible to solve Equation ( 80) with 5 flops using Algo. 10. There is no fill-in, so the resulting memory consumption is (LU( )) = ( ) = 3. Again, the LU decomposition is directly written over .

When using the extrapolation, the difference is that only the fine nodes are smoothed. Then, in the case of the fine nodes in the white smoother, the matrix becomes diagonal and can be inverted trivially for 1 flop per entry. We remind the reader that the white nodes represent half of all nodes, i.e., asymptotically we have = and = = 1/2, and that all nodes are either white or black so + = 1. By summing all of the lines, we finally get:

without extrapolation:

with extrapolation, since only the fine nodes are smoothed, i.e., half of the white colored nodes:

(LU( )) = 0 /2 + 10 /2 + 3 /2 = 1.5 + 3.5

Algorithm 7 Compute ( ) = ( ) ( ) : O (10 ) (here is the length of a circle line)

Input:

( )

Output:

( )

// update row + 1 3:

( + 1, ) = ( + 1, )/ ( , ) (1 div.) 4:

= ( + 1, ) 5:

( + 1, + 1) = ( + 1, + 1) -( , + 1) (1 add, 1 mult.) 6:

( + 1, ) = ( + 1, ) -( , ) (1 add, 1 mult.) 7:

// update last row 8:

( , ) = ( , )/ ( , ) (1 div.) 9:

= ( , ) 10:

( , + 1) = ( , + 1) -( , + 1) (1 add, 1 mult.) 11:

( , ) = ( , ) -( , ) (1 add, 1 mult.) 12: end for

Complete smoothing

The only step left is to perform the difference between and ( ) = ( ) 4:

( + 1) = ( + 1) -( + 1, ) ( ) (1 add, 1 mult.) 5:

( ) = ( ) -( , ) ( ) (1 add, 1 mult.) 6: end for 7: ( -1) = ( -1) 8: ( ) = ( ) -( , -1) ( -1) (1 add, 1 mult.) 9: ( ) = ( ) 10: // Backward substitution: = 11: ( ) = ( )/ ( , )(1 div.) 12: ( -1) = ( -1) -( -1, ) ( ) (1 add, 1 mult.) 13: for = -1, • • • , 2 do 14:

( ) = ( )/ ( , ) (1 div.) 15:

( -1) = ( -1) -( -1, ) ( ) -( -1, ) ( ) (2 add, 2 mult.) 16: end for 17: (1) = (1)/ (1, 1) (1 div.)

Algorithm 9 Compute ( ) = ( ) ( ) : O (3 ) (here is the length of a radial line)

Input:

( )

Output:

( ) contains the LU decomposition 1: for = 1, • • • , -1 do 2:

( + 1, ) = ( + 1, )/ ( , ) (1 div.) 3:

( + 1, + 1) = ( + 1, + 1) -( , + 1) ( + 1, ) (1 add, 1 mult.) 4: end for Algorithm 10 Solve ( ) = ( )/ ( , ) (1 div.) 10:

( -1) = ( -1) -( -1, ) ( ) (1 add, 1 mult.) 11: end for 12: (1) = (1)/ (1, 1) (1 div.)

during the multigrid cycles, we obtain the final complexity as