Convergence properties of local defect correction algorithm for the boundary element method

: Sometimes boundary value problems have isolated regions where the solution changes rapidly. Therefore, when solving numerically, one needs a ﬁne grid to capture the high activity. The ﬁne grid can be implemented as a composite coarse-ﬁne grid or as a global ﬁne grid. One cheaper way of obtaining the composite grid solution is the use of the local defect correction technique. The technique is an algorithm that combines a global coarse grid solution and a local ﬁne grid solution in an iterative way to estimate the solution on the corresponding composite grid. The algorithm is relatively new and its convergence properties have not been studied for the boundary element method. In this paper the objective is to determine convergence properties of the algorithm for the boundary element method. First, we formulate the algorithm as a ﬁxed point iterative scheme, which has also not been done before for the boundary element method, and then study the properties of the iteration matrix. Results show that we can always expect convergence. Therefore, the algorithm opens up a real alternative for application in the boundary element method for problems with localised regions of high activity.


Introduction
Often boundary value problems have small localised regions of high activity where the solution varies rapidly compared to the rest of the domain.This behaviour may be due to boundary conditions or due to an irregular boundary.One therefore has to use relatively fine meshes to capture the high activity.Since the activity is localised, one may also choose to solve on a uniform structured grid.The size of each grid is chosen in agreement with the activity of the solution in that part of the domain.The solution is thus approximated on a composite grid, which is the union of the various uniform local grids.One way of approximating this composite grid solution that is simple and less complex than directly solving on the composite grid is by Local Defect Correction (LDC).This approach is well developed and extensively studied in methods like the finite difference methods (FDM) [Anthonissen, Bennett and Smooke (2005); Ferket and Reusken (1996); Hackbusch (1984)].In Kakuba et al. [Kakuba and Anthonissen (2014)], a local defect correction algorithm for the boundary element method (BEM) is presented.It was shown therein that the algorithm is cheaper than solving directly on the composite grid.The BEM, being a method with full matrices, would benefit much from such an algorithm when solving problems with localised regions of high activity that occur at the boundary.The LDC technique employs a uniform global coarse grid that covers the whole boundary and then a uniform local fine grid in a small part of the boundary that contains the high activity.In Ferket et al. [Ferket and Reusken (1996); Hackbusch (1984)] LDC has been shown to be a useful way of approximating the composite grid solution in which a global coarse grid solution is improved by a local fine grid solution, through a process whereby the right hand side of the global coarse grid problem system of equations is corrected by the defect of a local fine grid approximation.The properties for this method in FDM have been well studied, see for instance [Anthonissen (2001); Ferket and Reusken (1996); Hackbusch (1984); Minero, Anthonissen and Mattheij (2006)].Literature is still scanty on the LDC technique for BEM, but error studies, especially for direct multigrid applications, are increasing [Qu and Cui (2014)].A good overview is also given in Feischl et al. [Feischl, Führer, Heuer et al. (2015)].The LDC technique will be particularly useful in applications for problems that exhibit multiscales in behaviour.An example is the modelling of impressed current cathodic protection systems where the potential problem is characterised by small regions of high activity, cathodes and anodes.In this paper we study convergence properties of the LDC algorithm for BEM which, to the best of our knowledge, has not been done.In Kakuba et al. [Kakuba and Anthonissen (2014)], the algorithm is presented using integral formulations.In this paper we base on the integral formulation to derive the corresponding fixed point iterative scheme in matrix form.This formulation assumes the boundary and the discretisation are such that, on refinement, the nodes of the global grid in the refinement area also belong to the fine grid.For the sake of completeness, we first briefly outline the development of BEM in Section 2. To build the discussion for the convergence properties, the important steps of the LDC technique for BEM, as published in Kakuba et al. [Kakuba and Anthonissen (2014)] are presented in Section 3, with an example.In Section 4, the LDC algorithm is formulated as a fixed point iterative scheme on the basis of which convergence properties are tested using examples.We conclude the study with a summary of the findings, in Section 5.

The boundary element method
Consider a closed domain Ω ⊂ R 2 , with boundary Γ.Denote by n the outward unit normal at Γ.We consider on Ω the two dimensional potential problem for which Dirichlet, Neumann or mixed boundary conditions may be defined as illustrated in Fig. 1.
where the coefficient c(s) is given by, where r is the variable field point, s is a fixed point, χ a coordinate along the boundary, Ω c is the complement of Ω in R 2 and α(s) is the internal angle at s.The derivation of this boundary integral equation (BIE) is readily available [Katsikadelis (2002); Paris and Canas (1997); Pozrikidis (2002)].The function v is the fundamental solution of the Laplace equation in R 2 given by: According to Eq. ( 2), the solution u can be computed at any point in the domain if we know u and its normal derivative everywhere at the boundary.The goal in BEM therefore is to look for the missing data at the boundary, which is the function u at the Neumann boundary or the normal derivative ∂u/∂n at the Dirichlet boundary.At the boundary, the discretised BIE leads to the linear system of equations where i, j = 1, 2, . . ., N , Γ j , j = 1, 2, . . ., N is a discretisation of Γ and N is the number of grid elements Γ j used in the discretisation.We have also introduced the vectors where q is defined as Using boundary conditions in Eq. ( 5) leads to the square system The vector x contains the unknown values of either u or q in the grid elements nodes at the boundary.The solution of the system (9) gives a BEM approximation of the unknowns in x in the grid nodes at the boundary.We denote by x L a BEM approximation on a grid of size L. Thus, u L j (or q L j ) is a BEM approximation of u j (or q j ), using a grid of size L. Solving Eq. ( 9) gives the unknown boundary quantities of u and q.Therefore, all the boundary quantities are available and the solution u i at any point r i ∈ Ω can then be computed using the identity (10)

local defect correction
In this section we briefly present the local defect correction algorithm that was introduced in Kakuba et al. [Kakuba and Anthonissen (2014)].The presentation focusses on the important steps of the algorithm that are necessary in the development of the same, as a fixed point iterative scheme that is discussed in this paper.Consider the potential problem: where The exact solution of this problem, shown in Fig. 2, has a small area close to the boundary where it changes rapidly.As a result, the solution u(r) at the boundary has a region of high activity in a small part of the boundary.The LDC algorithm distinguishes a small where |Γ L j | = L for all j = 1, 2, . . ., N .The number of elements of the local problem is denoted N local , and the individual elements of the local grid are denoted Γ l local .Thus, the local problem is solved on Ω local , using a uniform local fine grid Γ l local , of N local elements each of size l covering Γ local .That is, active [Kakuba and Anthonissen (2014)] constant elements that we use here, the collocation nodes are the midpoints of the elements.Let denote the set of coarse grid nodes.Let be the set of the local fine grid nodes.Then we denote the nodes of the local problem that belong to the active part of the boundary as r l active .We assume that all the grid nodes of r L ∩ r l active belong to r l active , Fig. 4. The composite grid nodes r l,L are the union r L ∪ r l active of the global coarse grid nodes r L and the active local fine grid nodes r l active .First, we have the following discretised integral equation on which leads to the initial global coarse grid system of equations Note that on the right hand side of Eq. ( 17), the function q(r) is given and, therefore, to make sure we are only measuring the errors due to the discretisation of u, we minimise the interpolation error [Kita and Kamiya (2001)] by evaluating that integral without discretising q, but by using an appropriate integration rule to evaluate that integral of the product of q and v. Once Eq. ( 18) is solved, the solution is used to complete the formulation of the local problem by computing u on Γ inside to be used as the Dirichlet boundary condition.The local problem on Ω local satisfies the same operator as in the global problem.Since Γ active ⊂ Γ, the boundary conditions on Γ active are the same as those in the global problem.The local problem is solved on a finer grid and the solution is used to estimate the defect as by the following explanation.
Consider the coarse grid discretisation Eq. ( 17).If we knew the exact solution u j := u(r j ) in the nodes, then using it in Eq. ( 17) would give: where d L i is the local defect for the i-th equation.We also have the exact BIE obtained by using the undiscretised exact function u on the elements, as Subtracting Eq. ( 20) from Eq. ( 19) gives where d L i is the defect of the i−th equation.If we had this defect, we would add it to the system Eq.( 18) and solve to obtain the exact solution at the nodes.The contribution to the defect d L i , from element j, is given by: so that the defect for the i−th equation is given by: The contributions to this defect are significant in the active region, and are assumed negligible elsewhere.However, the exact function u is not known and we cannot therefore use Eq. ( 21) to compute the defect d L i .What we can instead do is to formulate and solve a local problem using a fine grid in the active region.The solution to this problem is better than the coarse grid solution because of the fine grid used.At this point, the best solution available is and the best approximations to the integrals in Eq. ( 22) are In Eq. ( 25), we assume that in the local fine grid Γ l active , a global coarse grid element  Therefore, using the initial fine grid solution, the initial best approximation of the defect per element is for elements in the active region of the boundary, whereas it is assumed zero for elements outside the active region.We can then compute the defect The formulation described above can be summarised in matrix form as follows.On the local domain Ω local we have the local fine grid problem where and the vectors on Γ l local have been partitioned as The solution u l 0active , on the fine grid, which is expected to be more accurate than the coarse grid solution, is then used to approximate the defect for the global coarse grid problem So we solve the updated system Solving the system in Eq. ( 31) gives the updated coarse grid solution u L 1 .At this stage we use the fine grid solution on Γ l active and the global coarse grid solution to form a composite grid solution u l,L as The composite grid solution in Eq. ( 32) can now be used to compute better boundary conditions b l 1local on Γ inside , and then form and solve the updated fine grid problem Then we obtain the updated composite grid solution given by: This step marks the end of one complete cycle of the LDC algorithm.The iteration process is summarised in Algorithm 1, whose more detailed presentation is discussed in Kakuba et al. [Kakuba and Anthonissen (2014)], and results on an example are shown in Fig. 6 and Fig. 7, where functions are plotted only for the side of the square boundary with the active region.The authors in Kita et al. [Kita and Kamiya (2001)] present a broader review of adaptive error refinement techniques.In this paper, we compare the numerical solution with the exact solution to make conclusions on the error, a kind of error measurement similar to the one attributed to Mullen et al. [Mullen and Rencis (1985)].

(i) Initialisation
-Solve the global coarse grid system -Solve for the updated coarse grid solution -Form the updated composite grid solution u l,L i,i and compute the new defect 4 Properties of the LDC for BEM algorithm as a fixed point iterative scheme We consider the following Neumann problem To formulate the LDC algorithm as a fixed point iteration, we need a vector formulation for the steps in Algorithm 3.1.The first step of the algorithm is to solve a global coarse grid problem for an initial solution x L 0 .The solution x L 0 is a vector of u's except in one node, the last node, where a Diriclet boundary condition is prescribed to ensure uniqueness of the  2) for interior points, we compute the potential u(r) on Γ inside .That is, Introducing a vector g and a matrix HL such that then we can write Eq. ( 37) as Using Eq. ( 40), we obtain Dirichlet boundary conditions on Γ inside .The boundary conditions on Γ active are the same as the given boundary conditions in the global problem since, Γ active ⊂ Γ.Using (5), we can then write the equations on Γ l local in vector form as where H l local and G l local are the BEM H and G matrices on the fine grid of the local problem boundary Γ l local , u l active and q l active are vectors on the active part Γ l active of the local problem boundary, and u l 0inside and q l 0inside are vectors on Γ l inside , the part of the local problem boundary that is inside the global domain.The vector u l 0inside is known through Eq. (40) and the vector q active is known through the boundary conditions.So, we rearrange Eq. ( 41) as The matrix H l active is a block of H l for which the column index corresponds to nodes in Γ l active .Similarly H l inside is a block of H l for which the column index corresponds to nodes in Γ l inside .The blocks G l active and G l inside are defined analogously.The quantities on the right hand side of Eq. ( 42) are all known.Let x l 0local := Then we have (47) The solution of the system in Eq. ( 47) gives us another solution u l 0active in Γ active , which should be a better approximation of u(r) than u L 0 in Γ active , because of the fine grid used.We, then, use this solution to compute the defect and update the global coarse grid solution.The defect on an element Γ L j when the collocation node is i is given by: where ∪ k Γ l active,jk = Γ L j as illustrated in Fig. ( 5).The integration in Eq. ( 48) is computed at all the global elements that lie in Γ active , so that the total defect for the i-th collocation node is Since each node communicates with the local active region through integration, the defect d 0i is computed for all the nodes.Let us introduce a matrix H, defined as Let P L,l be a restriction from the fine grid Γ l active , to the coarse grid Γ L active in Γ active .Then we can write the defect d 0 as: where the matrix Ĥ is as in Eq. ( 7) with the superscript L indicating a coarse grid of size L. Now we have the defect for all the coarse grid nodes.We update the coarse grid system in Eq. ( 36) to obtain the updated coarse grid solution x L 1 , that is, (51) At this stage we can assemble a composite grid solution on Γ l,L that consists of the initial fine grid solution and the updated coarse grid solution.So, where u L 1 c is the updated coarse grid solution on Γ c outside the active region Γ active .To complete the updated composite grid solution, we need to solve a new local problem.To this end, we use the solution in Eq. ( 52) to compute another approximation of u(r) on Γ inside .Thus, we have where the matrix Hl,L is as defined in Eq. ( 40) but on the composite grid, that is, and r i ∈ Γ inside .Then we formulate an updated system for the local problem where Solving the system in Eq. ( 55) gives an updated solution u l 1active of u(r) on Γ active .At this stage we have a completely updated composite grid solution given by: This completes the first iteration that gives us the first updated composite grid solution.The process can be repeated until there is no more change in the solution.In what follows, we formulate the above process as a fixed point iterative process.Let I l active be an identity of size N active , the number of local elements in Γ active .Then, the part of the local solution in Γ active is given by, for the i-th iteration.Consider the updated composite grid solution in Eq. ( 56) for iteration i + 1.Using Eq. ( 55) and Eq. ( 57), we have From the second block row of Eq. ( 58) we have . This is the global coarse grid solution outside the active region.For a Neumann problem, we prescribe Dirichlet boundary conditions in the last node, in order to obtain a unique we can then write Eq. ( 65) as In Eq. ( 66) the updated solution on the active grid Γ l active is expressed in terms of the previous solution and the updated solution outside the active grid.To have an expression for the iteration that takes place on the active region alone, we use Eq. ( 60) to replace u L i+1,c .Thus, where b c is the vector of boundary conditions outside the active region.Since the last entry of x L i+1 c is a q-value and that of b c is a u-value, then matrices D 1 and D 2 are the projections and we have Introducing the notation we can write Eq. ( 69) as Using Eq. ( 71) in Eq. ( 66), we obtain ) which can be written as 72) can be written as The vector v remains fixed throughout the iteration, since q l active , g, b c , remain fixed, and remains the same throughout the iteration.Eq. ( 73) expresses the iteration that takes place on the fine grid Γ l active as a fixed point iteration with iteration matrix Q defined as Thus, we have This iteration will converge if the spectral radius of the iteration matrix Q is less than unity.In Tab. 1 and Tab. 2, we have the spectral radii of Q for different combinations of L and l.
All the values are less than unity implying convergence of the fixed point algorithm.
Consider the problem in Eq. ( 11).We identify Ω local as Ω local := [0.2, 0.8] × [0, 0.4], as illustrated in Fig. 2. The LDC is then used with various sizes of coarse grid size L and fine grid size l.We expect the ratios to be less than unity as well.In Tab. 3, we have computed these ratios for five iterations and different combinations of grid sizes.The results fit our expectations, to further illustrate guaranteed convergence of the algorithm.

Conclusions
The boundary element method is a relatively new method, whose development started in the late 1970's although the underlying theory of integral equations could be traced to earlier decades.Its biggest computational disadvantage is that the resulting matrices are full matrices, making the method expensive.In the present study, we have successfully opened a new direction in the implementation of the method, by formulating a fixed point iterative scheme for the LDC technique and showing numerically that the algorithm converges.This approach will be very useful in the implementation of the method in problems with The solution u(r) in part of ∂Ω that borders the high activity

Figure 3 :
Figure 3: An example of a multiscaled solution with localised high activity, in 3(a) and, in 3(b), an illustration of a local problem domain.The boundary of Ω local is Γ local := Γ active ∪ Γ inside[Kakuba and Anthonissen (2014)]

Figure 4 :
Figure 4: Global coarse and local fine grids.The dots are the nodes r l local of the local fine grid Γ l local and the big open circles are the nodes r L of the global coarse grid Γ L .Node 2 belongs to r L ∩ r lactive[Kakuba and Anthonissen (2014)]

Fig
A refinement of the coarse element Γ L j in (a) above

Figure 5 :
Figure 5: A coarse element that is refined into three elements in the local fine grid, such that Γ L j = 3 ∪ k=1 Γ l active,j k [Kakuba and Anthonissen (2014)]

Figure 6 :Figure 7 :
Figure 6: Results of a typical LDC process for a Neumann problem in one iteration (The solid line is the exact solution [Kakuba and Anthonissen (2014)]) where on Γ 1 we have Neumann boundary conditions, and on Γ 2 we have Dirichlet conditions [Kakuba and Anthonissen (2014)]If Γ 1 ≡ Γ, we have a Neumann problem, and if Γ 2 ≡ Γ, we have a Dirichlet problem otherwise we have a Mixed problem.In the solution of problem (1) using BEM, the problem is formulated as an integral equation

Table 1 :
Spectral radius of the iteration matrix Q for a Neumann problem for different combinations of fine and coarse grid sizes l and L respectively when local problem domain

Table 3 :
The ratios γ i defined in Eq. (76) when we use LDC to solve the problem in Eq. (11) high activity in the boundary that demand localised regions of high resolution grids.We have shown, by numerical experiments, that the resulting algorithm converges and thus provides a real alternative to composite grid solutions.There is still need to theoretically establish the convergence results, but given the complexity of the matrices involved, this paper presents a good opening to that line of research work.