{'fraction_non_alphanumeric': 0.04140974256204717, 'fraction_numerical': 0.014028056112224449, 'mean_word_length': 4.311873080859775, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A ghost-point smoothing strategy for geometric multigrid on curved boundaries

We present a Boundary Local Fourier Analysis (BLFA) to optimize the relaxation parameters of boundary conditions in a multigrid framework. The method is implemented to solve elliptic equations on curved domains embedded in a uniform Cartesian mesh, although it is designed to be extended for general PDEs in curved domains, wherever a multigrid technique can be implemented. The boundary is implicitly defined by a level-set function and a ghost-point technique is employed to treat the boundary conditions. Existing strategies in literature adopt a constant relaxation parameter on the whole boundary. In this paper, the relaxation parameters are optimized in terms of the distance between ghost points and boundary, with the goal of smoothing the residual along the tangential direction. Theoretical results are confirmed by several numerical tests in 1D, 2D and 3D, showing that the convergence factor associated with the smoothing on internal equations is not degraded by boundary effects.


Introduction
Multigrid (MG) method is among the most efficient iterative solver for discrete Partial Differential Equations (PDEs). Initially designed for elliptic equations, MG methods have been widely adopted to solve countless problems thanks to its optimal computation cost that scales linearly with respect to the number of computational nodes for sparse matrices, outperforming many other numerical solvers. Although initially MG approaches were proposed for simple domains such as rectangles in 2D and cuboids in 3D, the application to more complexshaped domains gained an increasing interest in the scientific community in the last years in several contexts [GFO18].
For example, in the context of energy storage problems one of the recent challenges in computer simulations is the representation of 3D microstructure effects of lithium-ion batteries. It is common practice to represent the complex geometry of this class of problems by adopting simplified approaches, such as union of spheres [WGB + 13]. In this way, the complexity of the geometry is more manageable, though it comes at the expense of accuracy. Representing complex-shaped objects is a bottleneck to countless applications in scientific computing, ranging from micro to macro scale problems, such as fluid dynamics past complex capillary structure in cardiovascular systems or hydrothermal circulations in subsurface reservoirs.
In the context of MG for complex-shaped domains, it is well known that boundary effects can propagate from the boundary to the overall domain and then degrade the efficiency of the entire MG unless a proper treatment is adopted. To this purpose, it is common practice in literature to add extra relaxations on the boundary conditions in order to improve the accuracy near the boundary and maintain the overall performance acceptable [BL11]. This approach has been adopted, for example, in [Coc20,CR18,CR13]. The rationale is that the extra computational effort required to perform the local boundary relaxations is negligible with respect to the overall internal relaxations and vanishes when h → 0, where h is the spatial step [TOS00]. This rationale is well justified for computational grids that have an almost uniform resolution in the whole domain. However, several real-world problems require a more refined mesh close to the boundary to capture the curvature of complex-shaped domains and to maintain the overall numerical error almost uniform everywhere. In these cases, the computational cost of extra relaxations along the boundary might dominate the overall cost, since the majority of the mesh is concentrated in the surrounding area of the boundary, resulting in a sub-optimal solver.
A way to overcome this difficulty is to adopt Algebraic Multigrid (AMG) methods [RS87], in which the MG operators are defined merely from the linear system. We refer the reader to [XZ17] for a complete review on the a general theory of AMG methods. The two main approaches used to setup the hierarchy of prolongation matrices in a completely algebraic setting are classical coarsening and coarsening by aggregation; see e.g. [BVEM00,Not10]. The main advantage of AMG methods is that they only requires matrices and can thus be easily integrated into existing code (see [DDF21] for a recent implementation of parallel coarsening algorithms based on aggregation of unknowns). On the other hand, in AMG methods one cannot easily incorporate a priori knowledge about the model problem, and more efficient solvers could be designed if geometrical information is included into the definition of the operators. In the context of MG methods for curved boundaries, this was highlighted in [LC05], where a comparison between geometric and algebraic approaches is performed.
Complex-shaped domains are usually modelled by fitted-boundary methods such as the Finite Element Method (FEM) [LOV18,Sha13,dPVvB + 20,DLOV17], that has been by far the dominant method to represent the geometry of domains, thanks to the use of isoparametric elements that allow a greater flexibility in the domain discretization. However, there are critical limitations in FEM when applied to highly complex-shaped domains: (i) the mesh generation process may become cumbersome when generating elements to conform highly varying curvatures of the boundary and (ii) the design and implementation of parallel solvers in FEM require impractical implementation efforts to partition the mesh for a balanced distribution of the computational load between cores. For this reason, unfitted-boundary methods such as those based on Finite Difference Methods (FDM) and level-set approaches [GFO18] have become more popular in recent years, since they allow a natural design of parallel solvers in the context of High Performance Computing and at the same time they do not require any mesh generation effort, since the domain is immersed into a fixed grid.
FDM in non-rectangular domains on unfitted grid is usually approached by the Ghost-Point method. The idea appeared for the first time in the Immersed Boundary Method [Pes77] to model blood flows in the heart, subsequently extended by LeVeque and Li in [LL94] to higher accuracy. More recent development can be found in [FAMO99,GF02,GF05,NMG09].
To capture the multi-scale nature of the problem, a very effective and popular approach is the Adaptive Mesh Refinement (AMR), where a uniform Cartesian grid is refined according to a recursive criterion and the entire grid structure is effectively stored in a quad-tree in 2D and oct-tree in 3D. MG methods for AMR (and its parallel implementation) have been elected by the community in scientific computing to be very efficient solvers for PDEs in highly complex-shaped geometries [DAB + 14, TE18,RAZA21]. However, the treatment of boundary conditions in MG methods has been so far undervalued, leading to suboptimal solvers.
In this paper, we improve the technique proposed in [CR13] and propose an efficient and accurate technique to embed boundary conditions in multigrid techniques on curved domains without the necessity of extra-relaxations on the boundary equations. The study is proposed for uniform Cartesian grids, providing a building-block to future implementations of MG for AMR, where the main benefit in terms of computational cost will be observed.
Following the approach proposed by one of the authors in [CR13] in the framework of finite-difference unfitted-boundary methods, a ghost-point technique will be adopted to enforce the boundary conditions, together with a relaxation technique for boundary equations which is based on adopting a fictitious time approach. For example, a Dirichlet boundary condition u = g is relaxed by the fictitious time evolution ∂u ∂τ = g − u. In [CR13], the fictitious time step was determined by a CFL stability condition to guarantee the convergence of the smoother. However, after a few iterations the boundary relaxations spoiled the smoothness of the solution and degraded the overall performance, needing some extra-relaxations along the boundary to maintain an appropriate overall performance. To overcome this issue, an improved technique will be proposed in this paper by adopting a local fictitious time step that varies along the boundary for each ghost point and that depends on geometrical properties like the distance of the ghost point from the actual boundary. The local fictitious time step will be determined by optimising the boundary relaxation parameters through a Boundary Local Fourier Analysis (BLFA) in order to maximise the smoothing property along the tangential direction. In this way, the overall performance is not degraded and then the MG will maintain an efficiency similar to that of classical MG on rectangular domains. The approach will be described in 1D, 2D and 3D. In the latter case, the implementation will be performed in C++ using PETSc [BAA + 22, BGMS97] to allow for a more natural parallel implementation. Our MG Poisson solver will be compared with the AGgregation-based algebraic MultiGrid (AGMG) of [Not,Not10].
The paper is organized as follows. In §2 we pose the problem in the 1D setting. Such case, that can indeed be easily managed by eliminated boundary conditions, is intended as starting point for the presentation and discussion of our approach in 2D and 3D settings given in §3 and §4, respectively. We end with §5 where we draw conclusions.

One-dimensional case
Consider a one-dimensional uniform grid G h = {x 0 , x 1 , . . . , x N } with spatial step h = x i − x i−1 , for i = 1, . . . , N . The elliptic boundary-value problem is where Bu(a) = g a is a generic boundary condition on x = a, with x N −1 < a ≤ x N (see Fig. 1). We call ϑ = (x N − a)/h. Then 0 ≤ ϑ < 1. In this paper we will discuss Dirichlet (Bu(a) = u(a)) and Neumann (Bu(a) = u (a)) boundary conditions. The extension to more general boundary conditions, such as Robin boundary conditions, is straightforward. The elliptic equation −∆u = f is discretized by central differences on x i ∈ Ω = (x 0 , a) while the boundary condition is discretized by Bu(a) = Bq(a), where q(x) is the polynomial of degree s − 1 that interpolates u on nodes u N , u N −1 , . . . , u N −s+1 , where s is the stencil size. In general, the discretization of the boundary condition for s ≤ 3 is where the coefficients c −k are the values (respectively derivatives) of the Lagrange polynomials associated to the nodes −s + 1, . . . , 0. In particular: (1, −2, 1)) for Neumann b.c. and s = 3.
The grid point x N is called a ghost point and u N is the ghost value. Although the boundary condition can be easily eliminated in 1D by solving (3) for u N and substituting into the internal equations (as it has been done for the boundary condition on x 0 , see Eq. 2), we maintain a non eliminated approach since we aim at the extension to higher dimensions, where the eliminated boundary condition approach is not possible (or extremely impractical). We note that a possible approach based on eliminated boundary conditions in higher dimensions for Dirichlet boundary conditions is proposed in [GF02], although the extension to Neumann boundary conditions is not straightforward.
The discretized problem is then a linear system

Multigrid method
We quickly recall the main idea of the multigrid method, referring the reader to more comprehensive literature for more details (see, for example, [TOS00]). First, a relaxation scheme is applied to the linear system M h u h = f h for a few iterations, say ν 1 pre-smoothing iterations, obtaining an approximated solution u * . Commonly, the relaxation scheme is of Richardson type, namely: for a suitable preconditioner matrix P . The residual r h = f h − M h u * h is evaluated on the grid G h and transferred to a coarser grid G 2h through a restriction operator r 2h = I h 2h r h , defined below. Then, the residual equation M 2h e 2h = r 2h is solved exactly on the coarse grid. The algebraic error e 2h is transferred back to the fine grid G h through an interpolation operator e h = I 2h h e 2h . The approximation u * h is then updated by u * h := u * h + e h and another few, say ν 2 , post-smoothing iterations of (4) are applied to the fine grid. The approximation u * represents the output of a single two-grid iteration (cycle). The cycle is iterated until the residual on the fine grid falls under a suitable tolerance.
The above approach is also called Two Grid Correction Scheme (TGCS), since only two grids are involved in the process, G h and G 2h . More grids can be adopted if the residual equation M 2h e 2h = r 2h is not solved exactly, but is instead approximated by γ cycles of the TGCS on grids G 2h and G 4h , and so on recursively, until a very coarse grid G H is reached and where the residual equation can be solved exactly at a low computational cost. Depending on the value of γ, we have the V-cycle (γ = 1) or W-cycle (γ = 2). It is uncommon to use higher values of γ since the it would not be convenient in terms computational effort [TOS00].
We now describe the relaxation scheme and the transfer operators that we adopt for the multigrid strategy in 1D.

Relaxation scheme
In order to have an effective multigrid strategy, the relaxation scheme (4) should be able to quickly reduce the high-oscillatory modes of the error e h (smoothing property), regardless of the the reduction of low-oscillatory modes. A relaxation scheme with this property is called a smoother. It is well known that the Jacobi scheme (when P = diag(M h ) is the diagonal part of M h ) is not a proper smoother. Relaxation schemes that have the smoothing property are weighted-Jacobi P = ωdiag(M h ) + (1 − ω)I with ω = 2/3 in 1D, where I is the identity matrix, and the Gauss-Seidel scheme, obtained when the preconditioner matrix P = L h is the lower triangular part of M h (including the diagonal). It can be proved that the convergence of both Jacobi and Gauss-Seidel schemes is not guaranteed for the presence of the last equation in the linear system, especially when ϑ > 1/2. For example, when ϑ is close to 1 then the diagonal entry c 0 of M h is in magnitude much larger than off-diagonal entries c −1 , c −2 and we observe numerically that the iteration scheme is not convergent.
To maintain the convergence, the ghost value u N is not updated by the classical Gauss-Seidel scheme. It is instead obtained by a relaxation of the boundary conditions on x = a: for some fictitious positive relaxation parameter ∆τ > 0 to be chosen in order to ensure the convergence of the entire iteration scheme. Observe that (5) reverts to the classical Gauss-Seidel scheme if ∆τ = 1/c 0 . With (5), the preconditioner matrix P is the lower triangular part of M h except for P N N that is replaced by We have used the symbol ∆τ because the iteration (5) can be also seen as the discretization of the evolutionary differential equation: in a fictitious time τ (see [CR13]), where the time derivative in the right hand side is discretized on x = x N by the explicit Euler scheme and the left hand side on x = a. The steady-state solution of (6) satisfies the boundary condition u(a) = g a .

Restriction operator
Following the approach of [CR13], the transfer operators (restriction I h 2h and interpolation I 2h h ) are defined as follows. The residual r h contains the residual of internal equations for x 1 ≤ x i ≤ x N −1 and that of the boundary condition for the ghost point x N . These two sub-residuals scale at different powers of h: internal equations scale at h 2 , while the relaxation of the boundary condition scales at h 0 or h for Dirichlet and Neumann boundary conditions, respectively. The restriction is then discontinuous across the boundary x = a and for this reason, the ghost value r N +1 is not used in the restriction of internal equations. In detail, if x + h < a the restriction is Full-Weighting (FW): while it is the average of the internal values if x < a but x + h ≥ a: The restriction of the residual of the boundary condition is trivial:

Interpolation operator
The interpolation operator acts on the error e 2h , which is continuous across the boundary, and therefore no modification from the standard linear interpolation is needed: Definition 1 We call Boundary Condition Multigrid (BCMG) method the multigrid technique proposed in this paper.

Numerical results in 1D
The efficiency of the multigrid is measured through the asymptotic convergence factor In order to avoid numerical issues due to machine precisions, the convergence factor is measured for the homogeneous problem (1) with f = g 0 = g a = 0 and with an initial guess of the multigrid cycle different from zero, say u (0) i is randomly uniformly distributed between −1 and 1 for all i. The asymptotic convergence factor ρ is then approximated by ρ = 1 5 15 k=11 ρ (m) . In Fig. 2 we plot the convergence factors versus ∆τ for different values of ϑ for the TGCS with N = 512, ν 1 = 2 and ν 2 = 1, for Dirichlet (left) and Neumann (right) boundary conditions. We report the results obtained with s = 3, but very similar results are obtained for s = 2. The dashed line represents the optimal convergence factor ρ OP T ≈ 0.08, that is the convergence factor achieved when ϑ = 0 and eliminated boundary conditions are implemented [TOS00] (i.e. the convergence factor of internal relaxations only). If the multigrid does not converge, we cut the convergence factor at ρ = 1 for graphical purposes. We observe that ρ is close to the optimal one (the convergence factor is not degraded by boundary effects of the relaxation of boundary conditions) or the multigrid does not converge (when ∆τ is too large). The plot shows that optimality is reached even for very small values of ∆τ , indicating that no relaxation is needed for the residual of the boundary condition. We point out however that this is a specific feature of the one-dimensional case, see also [BL11]. In fact in this case, the boundary residual is transferred as it is to the coarser grid and the error equation is solved exactly there. As we will see, this cannot be extended to the case in more space dimensions. The red dashed line is the optimal convergence factor without ghost nodes.

Two-dimensional case
In higher dimensions, the representation of the residual of boundary conditions on the coarser grid becomes a key point for the efficiency of the multigrid, since it is not an isolated value as in 1D, but it is instead represented by a set of ghost values just outside the domain and coarsening happens on ghost nodes too.
For all 2D problems we use a uniform Cartesian grid with h = ∆x = ∆y, albeit the method can be easily extended to anisotropic grids with ∆x = ∆y.
We consider the elliptic problem We implement an eliminated boundary conditions approach on Γ B (similarly to the 1D case (1) for x = x 0 ) and a ghost-point technique on Γ, explained in §3.2.

Set-up of 2D and 3D numerical tests and computation of the convergence factor
Unless otherwise stated, we perform in all numerical tests the TGCS with ν 1 = 2, ν 2 = 1, GS-LEX as a smoother and FW restriction operator. The convergence factor is measured for the homogeneous problem (8) with f = g B = g Γ = 0 and stencil ϑh with an initial guess of the multigrid cycle different from zero and in particular randomly uniformly distributed u ij ∼ (−1, 1) for all i and j. The asymptotic convergence factor ρ is then approximated by ρ = 1 5 15 k=11 ρ (m) .

Rectangular domains
Consider the rectangular domain Fig. 3. For the multigrid method, we use the technique and notations introduced in §2.1 for the 1D case. For the 2D case, we have to define the relaxation operator, and the transfer (restriction and interpolation) operators.
Relaxation scheme. The value on grid points of Γ B is set from the boundary condition, so u The value on the internal grid points is obtained by the standard Gauss-Seidel lexicographic (GS-LEX) iteration: , for i, j = 1, . . . , N − 1.
(9) Finally, a fictitious-time relaxation (similar to the 1D case (5)) is performed on the boundary Γ, i.e. x = x V , for which the ghost points are X N j = (x N , y j ) for j = 1, . . . , N − 1: Restriction operator: Full-Weighting (FW) With the notation T v ij , where T is a 3 × 3 matrix with elements T rs for r, s, = −1, 0, 1, we denote the value: Let u * ij be the approximation after ν 1 relaxations on the fine grid. The residual after the internal relaxations (9) is transferred to the coarser grid by the usual Full-Weighting (FW) 2D ninepoint stencil (see e.g. [TOS00] for more details): The residual after ghost value iterations (10) scales differently than the internal one and is thus transferred via the 1D operator acting only on ghost values: Note that in the present test case setup, FW on internal equations never mixes internal with ghost residuals.
Interpolation operator. Standard 2D bilinear interpolation operator is applied all over the grid (internal and ghost grid points).

Dirichlet boundary conditions: Boundary Local Fourier Analysis (BLFA)
For Dirichlet boundary conditions on x = x V , the coefficients c −2 , c −1 and c 0 are defined in (3). Assume a stencil size of s = 3. The problem has (N − 1) 2 internal points (x i , y j ) for i, j = 1, . . . , N − 1 and N − 1 ghost points (x N , y j ) for j = 1, . . . , N − 1, for a total of (N − 1)N active points. The discrete problem can be written as a linear system M h u h = f h , with: is the identity matrix, The relaxation scheme can be written in matrix form as , with the preconditioner matrix where To discriminate the contribution of the boundary relaxation from the internal equations, we analyze an alternative relaxation scheme obtained with In other words, the internal variables are solved in block and the ghost values u (m+1) N,j are updated by the relaxation (10). This alternative relaxation scheme is of course impractical and used here only for theoretical studies.
We call REL-ALL the TGCS associated with (11) and REL-BC the one associated with (12).
In Fig. 4 we compare the convergence factors obtained with REL-ALL and REL-BC for different values of ϑ ∈ [0, 1]. We see that, for different values of ϑ (different colors), the convergence factor obtained with REL-BC (solid line) performs similarly to REL-ALL (dashed line) or better. If it does not converge, the convergence factor is cut artificially at ρ = 1 for graphical purposes. Approximately, we observe that: where ρ OPT is the convergence factor obtained for regular rectangular domains Ω = [−1, 1] 2 where the ghost point technique is not employed. For ν 1 = 2, ν 2 = 1, GS-LEX as a smoother and FW restriction operator we have as predicted by the LFA (see [TOS00, Table 4.1]). Fig. 4 suggests that for each value of ϑ, there exists a range of possible values of ∆τ such that ρ REL-ALL ≈ ρ OPT . Therefore, it is not necessary to find ∆τ that minimizes ρ REL-BC , it is instead sufficient to identify a value of ∆τ in that range. However, it is clear from the figure that ∆τ ≈ 1.75 would yield ρ REL-ALL ≈ 0.3 for any ϑ (see Remark 4 for an explanation of these values), but no value of ∆τ can yield ρ REL-ALL ≈ ρ OPT for any ϑ.
We want to find ∆τ such that the smoothing property of the relaxation scheme is optimal. To this end, we want to quantify how much the relaxation scheme reduces the high-oscillatory Fourier modes along the tangential direction of the boundary. Let be a Fourier mode along the tangential (vertical) direction, with α ∈ [0, π].
Observe that we use i to denote the imaginary unit i 2 = −1, to be distinguished from i, that will be used later to denote some indices. The Fourier mode ϕ is a highly-oscillatory mode if α ∈ [π/2, π].
To perform the BLFA, we work on a semi-infinite domain: and we analyze the smoothing of REL-BC along the tangential direction. As explained above, the relaxation REL-BC is performed in two steps. The first step consists of solving the equations on the internal grid points in block when the ghost values are set to ϕ. Without loss of generality, we can work with the homogeneous problem, i.e. f i,j = 0, which means solving the problem The exact solution of this difference equation is: where The second step consists of relaxing the ghost values by (10) with g j = 0: =⇒ u (1) (20) The right-hand side of (20) is the amplification factor of the tangential Fourier mode ϕ = e iαy/h for a given ϑ ∈ [0, 1) and a given ∆τ > 0, that we denote as: where the coefficients c −k are given by (3) and are Finding the optimal ∆τ for a given ϑ requires to minimize the largest amplification factor for the high-oscillatory modes α ∈ [π/2, π] and thus Lemma 1 Let G 0 (α, ·) : (π/2, π) → R be the continuous function of α defined in (21) with coefficients (22). Then G 0 does not have stationary points in (π/2, π).
Proposition 1 The optimal ∆τ defined in (23) is obtained by the following formula: where G 0 is defined in (21) with coefficients (22).
In Fig. 6 we observe a rather satisfactory match between the optimal ∆τ obtained from numerical results and from the BLFA (25). . We observe that this choice will lead to a completely wrong estimate of ∆τ OPT . In fact, one can say that numerical and exact solution would differ by u EXA ij − u ij = O(h 2 ), but this is not true for a fixed value of α, since the boundary condition scale with h, namely (15) depends on h. This means that u EXA would only be true if the boundary condition did not scale with h, meaning that (from (15)) α/h = const., that is not the case of our interest for the BLFA. The discrepancy between the optimal ∆τ OPT obtained with (18) and the one obtained with p(α) = α is observed numerically in Fig. 6 for the 2D case and Fig. 17 for the 3D case.
In order to show that high frequency components of the error are dumped much quicker when ∆τ = ∆τ OPT , we compute the Discrete Fourier Transform of the residual along the ghost points (namely along the tangential direction of the boundary) and plot the components of the residual versus the frequency of the Fourier modes (15). In detail, let r BDY ∈ R N −1 be the residual along the ghost points defined as: (r BDY ) j = r(x N , y j ) for j = 1, . . . , N − 1.
Let y ∈ R N −1 be the y−coordinate vector with components y j = j N , for j = 1, . . . , N − 1. The Fourier modes (15) are discretized as The residual r B can be decomposed as a linear combination of Fourier modes where the components R k ∈ C are given by the Discrete Fourier Transform formula: Since r BDY is a real vector, then R k1 = R k2 if α k1 = −α k2 , where R k1 denotes the complex conjugate of R k1 . Therefore, we can restrict our study to α k ∈ (0, π). The Fourier modes ϕ (α k ) (y) are divided in low-frequency modes (α k ∈ (0, π/2)) and high-frequency modes (α k ∈ (π/2, π)). For each value of α k ∈ (0, π) we plot 2 |R k | so that we include the contribution of the symmetric mode in (−π, 0). Fig. 7 shows the amplitudes 2 |R k | versus α k ∈ (0, π) where r BDY is the residual on the ghost points after 10 relaxations on the fine grid, starting from a random initial guess whose values are randomly uniformly distributed between −1 and 1. A vertical line at α = π/2 separates the low frequency components (left) from the high frequency components (right). We repeat the test for ϑ = 0.2 (first row), ϑ = 0.5 (second row) and ϑ = 0.8 (third row). For each value of ϑ we choose three values of ∆τ : the optimal value (middle column), a lower value (left column) and a higher value (right column). We observe that the smoothing property is clearly satisfied only when ∆τ = ∆τ OPT , namely the ratio between high-and low-frequency component amplitudes is smaller. The plots are obtained for N = 128.

Neumann boundary conditions: BLFA
We use the same approach as for Dirichlet boundary conditions. We replace the relaxation (19) with the Neumann boundary relaxation: The right-hand side of (26) is the amplification factor along the tangential direction of the Fourier mode ϕ = e iαy/h for the Neumann boundary condition. Here where The definition of the optimal ∆τ is the same as (23).

Proof. We have
To complete the proof, it is sufficient to show that G 0 (·, ϑ) does not have stationary points in (0, 1). We have since p(α) > 0. Therefore, there are no stationary point of G 0 (·, ϑ) in (0, 1).

Proposition 2
The optimal ∆τ for Neumann boundary conditions is obtained by the following formula: where G 0 is defined in (27) with coefficients (28).
The proof is analogous to Proposition 1.

Non-rectangular domains
For general domains we adopt the ghost-point method and the level-set approach used in [CR13]. In that paper, the authors adopted a multigrid approach with a fictitious-time dependent relaxation scheme using a constant ∆τ along the boundary. Let us briefly recall the details. As for the rectangular case, let where N ij = {(x i±1 , y j ), (x i , y j±1 )} is the set of four neighbour grid points of (x i , y j ); I h is the set of inactive points, defined as I h = D h \(Ω h ∪ Γ h ). We order the internal and ghost points with, for example, a lexicographic order (a different order can also be chosen, as the method presented in this paper does not rely on a specific ordering). Then Ω h = {X 1 , X 2 , . . . , X N I } and Γ h = {X N I +1 , X N I +2 , . . . , X N I +N G }, where N I and N G are the total number of internal and ghost points, respectively. The discrete solution of (8) is represented by a column vector u h = [U 1 , U 2 , . . . , U N I +N G ] T . For simplicity, we refer to the discrete value of the solution in a grid point as either U k or u i,j , depending on the context. For each internal grid point X k ∈ Ω h we discretize −∆u = f using standard five-point stencil central differencing: For each ghost point G = X k ∈ Γ h we discretize the boundary condition Bu = g Γ as follows (details can be found in [CR13]). Firstly, we compute  the orthogonal projection Q of G onto the boundary Γ (see Fig. 8) by solving φ(G − η h n) = 0 for η ∈ [0, 1), where n is the outward unit normal vector that can be approximated by a finite-difference discretization of n = ∇φ/|∇φ| centered at the ghost point. Secondly, we identify the 3 × 3 stencil ST G that encloses Q and has G in one of its four corners. In detail, ST G is taken in an upwind fashion with respect to the normal direction (see Fig. 8). Such a large stencil ensures second order accuracy for the solution u and for the reconstruction of the gradient ∇u for both Dirichlet and Neumann boundary conditions. If one is interested only in the accuracy of the solution and not of the gradient, a 2 × 2 stencil will suffice to achieve second order accuracy for Dirichlet boundary conditions (see [CR13]). Finally, the discretized boundary condition is whereũ h is the biquadratic interpolant on the stencil ST G . Observe that, provided that the grid is sufficiently fine, ST G does not contain any inactive point. Equations (31) cannot be explicitly solved for the ghost values, since the stencil ST G can contain ghost points other than G and then these equations would be fully coupled each other. Equations (30) and (31) for internal and ghost values are then collected into a sparse linear system, that we denote as This linear system is solved by a multigrid approach. The main structure of the TGCS is the same as in the beginning of §2.1. Smoother and transfer operators are defined as follows.

Relaxation scheme.
As in the rectangular case, the smoother consists of a Gauss-Seidel scheme on internal points and of a relaxation of the boundary conditions on each ghost point G: u In [CR13], the value of ∆τ was chosen to satisfy a stability condition and ensure convergence of the iterative scheme to the numerical solution. However, choosing a fixed ∆τ in 2D and 3D (e.g. ∆τ = 1 and ∆τ = h for Dirichlet and Neumann boundary conditions, respectively) for all ghost iterations, although it would lead to a convergent scheme, it would degrade the performance of the multigrid strategy since boundary effects will be quickly propagated into the interior of the domain. This effect is mainly caused by the fact that adjacent ghost values are updated at different speeds depending on their distance from the boundary, generating highly-oscillatory modes that are not easily dumped from the multigrid strategy. On the other hand, for rectangular domains ghost points are located at the same distance from the boundary and a degradation effect is also observed if ∆τ is not properly chosen, as observed in §3.2. A possible solution to overcome the issue is presented in [CR13] by over-relaxing the ghost value equations, namely adding some extra relaxation sweeps on ghost points per each interior one. The computational effort introduced by over relaxations is negligible in uniform grids, but it might become more evident in adaptive grids where the majority of grid points could be concentrated close to the boundary. In this paper, we adopt a heuristic approach based on the BLFA performed in §3.2 for rectangular domains. In detail, we perform the ghost iterations (32) by choosing (see (25)): where · 2 is the Euclidean norm (distance between G and Q) and ∆τ OPT (θ) is defined in (25) and (29) for Dirichlet and Neumann boundary conditions, respectively.
Remark 2 For a curved domain, the values ofθ are fully distributed between 0 and 1, highlighting the benefit of using a variable ∆τ that depends onθ (see Fig. 10). We point out that the case ofθ close to 1 is observed only for very specific configurations. Therefore, the distribution is almost uniform forθ far from 1 and smoothly decreases whenθ approaches 1.
Remark 3 We observe that (21) and (27) can be thought as quadratic interpolations of exp(p(α)x) and p(α) exp(p(α)x), respectively, on the boundary x = x V . This suggests that we can use G 0 (α, ϑ) = exp(−p(α)ϑ) and G 0 (α, ϑ) = p(α) exp(−p(α)ϑ) for Dirichlet and Neumann boundary conditions instead of (21) and (27). In fact, (21) and (27) where computed from the exact discretizations on the rectangular domain and there is no specific benefit to use them in a non-rectangular domain, where the discretization is more articulated. The approximations proposed in this remark are therefore adopted for non-rectangular domains and we observed numerically that they actually perform slightly better than (21) and (27) (not shown).

Transfer operators.
The residual of internal iterations is transferred to the coarser grid by the FW operator (7) if all nine points of the FW stencil are all internal points, otherwise it is performed as an average of the residual only on the internal points of the nine-point FW stencil (see Fig. 9). The rationale is that the residual of internal and ghost iterations scale with a different power of h and then one should not mix them up, otherwise efficiency and sometimes even convergence cannot be guaranteed.
The residual of the ghost iterations is firstly extrapolated along the normal direction to the inactive points and then transferred to the coarser grid in the same way as the internal residual, namely by performing the FW if the nine points are all ghost or inactive points or as an average of external points only if any of the nine points is an internal one. Extrapolation to inactive points can be achieved, for example, by performing few iterations to solve ∂r ∂σ + ∂r ∂n = 0 on each inactive point (or at least in a narrow band close to the boundary), discretized by Euler explicit in σ and central difference for the normal derivative.
Interpolation from the coarse to the fine grid is performed by the standard bilinear interpolation operator on the whole computational domain without distinguishing between internal and ghost points (after extrapolating the ghost values along the normal direction to the inactive points). The rational is that the error to be interpolated is the solution of the coarse grid residual equation and then it is expected to scale with the same order for internal and ghost points.
Remark 4 If we want to use a constant ∆τ for all ghost points and aim at achieving, if possible, a convergence factor of ρ TARGET , the best that we can infer from the numerical results of the rectangular case is checking whether there exists a value of ∆τ such that ρ REL-BC (ϑ, ∆τ ) ≤ ρ TARGET for all ϑ ∈ [0, 1), where ρ REL-BC (ϑ, ∆τ ) is the convergence factor of REL-BC in the rectangular domain of §3.2. From Fig. 4, we can see that the minimum ρ TARGET for which we have a value of ∆τ that satisfies the condition above is ρ TARGET ≈ 0.3, achieved for ∆τ ≈ 1.75 for Dirichlet boundary conditions and ρ TARGET ≈ 0.15, achieved for h −1 ∆τ ≈ 1 for Neumann boundary conditions. This is confirmed by the numerical results presented in this paper on non-rectangular domains, for example Figs. 12, 14 and 15. In particular, the use of the optimal ∆τ proposed in the BCMG method is more important for Dirichlet than for Neumann boundary conditions.
In the left panel of Fig. 12 we plot the convergence factors of the TGCS (with set-up from §3.1) for different values of N . The solid line represents the test when ∆τ is constant for all ghost points, while the dashed line is the test with the optimal ∆τ (θ) defined in (33). If we choose a constant ∆τ for all ghost points, the best convergence factor that we can achieve is ρ ≈ 0.3 obtained with ∆τ ≈ 1.75 (see Remark 4 for an explanation of these values). On the other hand, the implementation of a variable ∆τ (θ) defined in (33) considerably improves the efficiency. The plot in fact suggests that the convergence factor approaches (when N goes to infinity) the optimal one predicted by the LFA for the standard case Ω = [−1, 1] 2 , which is ρ OPT ≈ 0.119 (when ν 1 = 2, ν 2 = 1 and FW restriction operator is adopted). Similar results are observe for the flower-shaped domain (right panel of Fig. 12), although the convergence factors ρ for ∆τ OPT are larger than the ones observed in a circular domain. However, we believe that ρ approaches the optimal value ρ OPT ≈ 0.119 when N goes to infinity. The discrepancy between the two tests can be explained by the fact that the boundary curvatures are not the same. In other words, the convergence factors between that two tests can be compared only when the ratio between boundary curvature and spatial step is similar. Also, it might be that the heuristic approach is a first approximation that neglects some additional features that are not accounted for in the present BCMG method, such as the relation with the local curvature of the boundary. A comprehensive theoretical analysis of the non-rectangular case is much more cumbersome and out of scope here.
In summary, we can suggest that the convergence factors of the TGCS is not much degraded by the ghost-point technique, provided that N is sufficiently large (or, in other words, the resolution is greater near the boundary as ensured, for example, in the AMR framework).

Comparison with the AGMG method
As a benchmark, we compare the performance of our BCMG method with that of an advanced algebraic multigrid method, namely, AGMG, cf. [Not,Not10,NN12,Not12]. AGMG is an aggregation-based algebraic multigrid method that performs the so-called K-cycle, i.e., two Krylov accelerated iterations at each intermediate level, and one forward and one backward Gauss-Seidel sweep for pre-and post-smoothing, respectively. AGMG is a powerful software designed to work as a black-box algorithm to solve efficiently sparse linear systems arising from discretizations of elliptic-like differential equations. Great flexibility is allowed, since it does not rely on the type of mesh (structured/unstructured, fitted/unfitted, etc.) or the particular discretization of the boundary conditions. This flexibility can come to a cost of under-performing against solvers that are designed for a specific problem. In Fig. 13 we compare AGMG with BCMG, showing that BCMG performs better than AGMG both on the circular and flower-shaped domains of Fig. 11. The asymptotic slope of the lines is an approximation of the asymptotic convergence factor, that is observed numerically to be around 0.2 for BCMG and 0.5 for AGMG for the circular domain, while it is around 0.3 for BCMG and 0.8 for AGMG for the flower-shaped domain. Moreover, the initial slope of BCMG is even steeper, suggesting that the initial convergence factor is dictated by the internal relaxations, while the asymptotic convergence factor is affected by small boundary effects (this behaviour was actually observed numerically in all tests performed in this paper). On the other hand, AGMG starts slowly, with a shallower slope compared to later iteration steps. All these aspects suggest that BCMG will fall below a desired tolerance on the residual much earlier than AGMG, as shown in Fig. 13.

Numerical results in 2D: Neumann boundary conditions
If Neumann boundary conditions are prescribed on the whole boundary, then the problem does not admit a unique solution. In particular, it can admit either infinite solutions (up to an additive constant) or no solutions (when the boundary condition and the source are not compatible). This property is also reflected at the discrete level, where the matrix of coefficients is singular. The standard technique to approach these problems is to deal with an augmented linear system, with an additional unknown (that represents a measure of the compatibility condition misfit) and an additional equation (usually the average of the solution must vanish to guarantee uniqueness), leading to a non-singular matrix of coefficients. These aspects, together with the standard multigrid techniques to face these scenarios, are detailed, for example, in [TOS00]. In this section we want to focus on the optimality of the relaxation parameters for Neumann boundary conditions, without generating unnecessary complications from the singularity of the problem that can be treated by standard techniques as anticipated above. To this purpose, we circumvent the singularity aspect in two ways: (i) by prescribing mixed boundary conditions, namely both Dirichlet and Neumann boundary conditions are applied on two different parts of the boundary, respectively (Test1N); (ii) by adding a reaction term to the Poisson equation, namely the differential equation of (8) becomes: Test1N can be achieved in the problem (8) considering a domain for which Γ B = ∅ and on which Dirichlet boundary conditions are prescribed, while we apply Neumann boundary conditions on Γ. In particular, we choose a levelset φ that represents a straight line passing through the point of coordinates (1 − 7h/10, 0) and with a negative slope of −1/ √ 3: With this choice, the distance between ghost points and the boundary is not constant; in fact it is highly varying along the tangential direction, ensuring that the main difficulties posed by non-rectangular domains are present in this numerical test. The domain is represented in Fig. 14, left panel, while in the right panel we plot the convergence factors of the TGCS (with set-up from §3.1) for different values of N . The solid line represents the test when ∆τ is constant for all ghost points, while the dashed line is with the optimal ∆τ (θ) defined in Eq. (33). With a constant ∆τ for all ghost points, the best convergence factor that we can achieve is ρ ≈ 0.25 obtained with h −1 ∆τ ≈ 1 (see Remark 4 for an explanation of these values). Efficiency is improved with a variable ∆τ (θ) as defined in (33), with a convergence factor approaching the optimal one (14) predicted by the LFA for the standard case Ω = [−1, 1] 2 , that is ρ OPT ≈ 0.119. In Test2N we consider a reaction-diffusion equation with homogeneous Neumann boundary conditions on the boundary of the circular domain (34) and the flower-shaped domain (35). Results are presented in Fig. 15, with a similar behaviour of Test1N discussed above, although with slightly larger asymptotic convergence factors, especially for the flower-shaped domain. As for the Dirichlet case, this is due to the high curvature of the boundary and the convergence factor would improve for finer grids.

Three-Dimensional case
As in the previous section, we use a uniform Cartesian grid with h = ∆x = ∆y = ∆z. The computational domain is D = [−1, 1] 3 , covered by a grid with (N + 1) 3 grid points X ijk = (x i , y j , z k ), where x i = −1 + i h, y j = −1 + j h and z k = −1 + k h for i, j, k = 0, . . . , N and h = 2/N . The elliptic problem is formally the same as (8).
As in 2D, we aim at comparing the convergence factors with the optimal one ρ OPT that is obtained for regular rectangular domains Ω = [−1, 1] 3 where the ghost point technique is not necessary. For ν 1 = 2, ν 2 = 1, GS-LEX as a smoother and FW restriction operator we have obtained numerically that We mimic the approach adopted for the analysis in two space dimension, but in 3D the boundary is a surface and then the residual of the boundary conditions is a discrete approximation of a two-dimensional function instead of a one-dimensional function as in the 2D case. This means that the BLFA will be based on a two-dimensional analysis.

Rectangular domains
The boundary to be treated with the ghost point technique, in this case is a plane x = x V . In particular we consider the rectangular domain where 1 − h < x V ≤ 1. We call ϑ = (1 − x V )/h. Then 0 < ϑ ≤ 1. The multigrid ingredients, such as relaxation scheme, restriction and interpolation operators, can be easily extended from the 2D treatment described in §3.2. Full details are provided in [CR13].
Remark 5 A similar argument of Remark 1 applies to the 3D case, where the exact solution is obtained with p(α 1 , α 2 ) = α 1 + α 2 . The discrepancy between the optimal ∆τ OPT is observed in Fig. 17.

Neumann boundary conditions: BLFA
The BLFA for the Neumann boundary condition is analogous to §4.1.1 for the 3D aspect of the analysis and to §3.2.2 for the treatment of the boundary condition. Therefore, using the same approaches the following proposition can be easily proved for the Neumann case.
Proposition 4 The optimal ∆τ for Neumann boundary conditions in 3D is obtained by the following formula: with G min 0 (ϑ) and G max 0 (ϑ) defined as in (44) and G 0 defined in (40) with coefficients (28).

Non-rectangular domains
The approach is easily extended from the 2D treatment of general domains of §3.3. Relaxation of boundary conditions is performed by (32) with (33). The optimal relaxation parameter is determined by (43) and (47) for Dirichlet and Neumann boundary conditions, respectively. Also, we use the same approach suggested in Remark 3 for the 2D case.

Numerical results in 3D: Dirichlet boundary conditions
In this section we test the multigrid efficiency on an ellipsoidal domain with Dirichlet boundary conditions. The level-set function is φ(x, y, z) = x − √ 2/20 0.686 In the left panel of Fig. 19 we plot the convergence factors of the TGCS (with set-up from §3.1) for different values of N . The solid line represents the test when ∆τ is constant for all ghost points, while the dashed line is the test with the optimal ∆τ (θ) defined in (33) with (43). If we choose a constant ∆τ for all ghost points, the best convergence factor that we can achieve is ρ ≈ 0.5 obtained with ∆τ ≈ 1.75 (a possible explanation would be based on the same argument that was discussed in Remark 4 for the 2D case). A variable ∆τ (θ) defined in (33) with (43) improves the efficiency, approaching the optimal convergence factor for the standard case Ω = [−1, 1] 3 , that is ρ OPT ≈ 0.155 (see (36)).

Numerical results in 3D: Neumann boundary conditions
Similarly to §3.3.5 (Test2N) for the 2D case, we solve a reaction-diffusion problem with Neumann boundary conditions on the whole boundary. The domain is the same ellipsoid of §4.2.1. In the right panel of Fig. 19 we plot the convergence factors of the TGCS (with set-up from §3.1). Also in this case we observe an improvement of the convergence factor when the optimal ∆τ (θ) is chosen (defined in (33) with (47)). The efficiency is considerably improved with respect to the case of constant ∆τ . The discrepancy of results between constant and optimal ∆τ strategies is even more pronounced than the 2D case. Therefore, the implication of adopting an optimal ∆τ that varies according to the distance of ghost points from the boundary remains a reasonable strategy to achieve an appropriate efficiency in view of future applications, such as AMR on larger 3D problems.

Conclusions
We have presented a Boundary Local Fourier Analysis (BLFA) theory to design an efficient relaxation scheme that is able to smooth the residual of the boundary conditions along the tangential direction, without affecting the smoothing performance of the internal equations. The smoother is embedded in a multigrid framework and the convergence factor is not degraded from boundary effects.
The model problem is addressed with unfitted boundary methods and a ghost-point technique is adopted to discretize the boundary conditions. It is known that equations for the boundary conditions pose problems to iterative solvers since their residuals behave differently than the residuals for internal points. In the context of smoothers for multigrid methods, boundary conditions can be relaxed by using a fictitious time evolution, where the time step ∆τ is a relaxation parameter.
In [CR13] the authors adopted a constant relaxation parameter for all ghost points and performed some over-relaxations on the boundary in order to mitigate the boundary effects. In this paper, we substantially improve the efficiency by adopting a variable ∆τ that depends on the distance from the ghost point to the boundary and is able to smooths the residual along the tangential direction. The optimal value is found theoretically by a BLFA and confirmed from numerical experiments. The improvement of the efficiency is justified when the majority of grid points is located in the vicinity of the boundary, such as in the AMR framework. In this paper we perform the analysis on a uniform Cartesian grid and we leave the application to AMR approaches to future efforts.
For rectangular domains where a boundary does not pass through grid points, the prediction from the BLFA is excellent, both for Dirichlet and Neumann boundary conditions. For curved boundary, there might be a small degra-dation of the convergence factor and this becomes more evident for complexshaped domains such as the flower-shaped one. We believe that this degradation will vanish when the grid is sufficiently refined. Also, the theoretical argument that is adopted for non-rectangular domain is a heuristic extension of the analysis performed for rectangular domains. Improvement of the theoretical parameters for the non-rectangular case can be performed by a more comprehensive analysis that accounts, for example, for the local curvature of the boundary.
We finally point out that the BLFA method presented here does not rely on the internal smoother, namely on the smoother adopted for the internal equation (weighted-Jacobi, Gauss-Seidel Red-Black, Gauss-Seidel LEX, etc.). In other words, the scope of this paper is to provide a general method to find optimal relaxation parameters for the boundary conditions regardless of the other ingredients of the multigrid and such that the convergence factor is not degraded from boundary effects. For a specific problem, a proper choice of the other multigrid components (internal smoother, transfer operators, etc.) must be performed carefully in order to optimize the efficiency. Also, the quantitative values shown in this paper, such as those of Remark 4, strongly depend on these choices. Finally, the BLFA will only have the task of embedding the relaxation of boundary conditions into the multigrid method while maintaining the optimal convergence factor.