key: cord-0045682-1uwvrm7y
authors: Kużelewski, Andrzej; Zieniuk, Eugeniusz; Bołtuć, Agnieszka; Szerszeń, Krzystof
title: Modified Binary Tree in the Fast PIES for 2D Problems with Complex Shapes
date: 2020-06-15
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50417-5_1
sha: 67b3816a0c739bc8abc0e3e8bb74c01459d9b565
doc_id: 45682
cord_uid: 1uwvrm7y

The paper presents a modified binary tree in the fast multipole method (FMM) included into the modified parametric integral equations system (PIES), called the fast PIES, in solving potential 2D boundary value problems with complex shapes. The modified binary tree proposed in this paper is built based on a one-dimensional reference system contrary to a quad-tree (based on a two-dimensional reference system) which is applied in the fast multipole boundary element method (FM-BEM). Application of the proposed tree allows reducing the number of numerical computations performed during its construction and fast multipole calculations in the fast PIES. The proposed modification of the tree in the fast PIES allows obtaining accurate solutions in engineering problems with complex shapes on a standard personal computer in a short time.

The fast multipole method (FMM) was initially proposed by Rokhlin [1] to accelerate the solution of 2D potential problems using boundary integral equations (BIE). Its main advantage is the reduction of the computational complexity of the matrix-vector multiplication from O(N 2 ) to O(N ) (where N is the size of the matrix) combined with an iterative solver. In papers [2, 3] , Greengard refined the algorithm by adding decomposition of the domain using the hierarchical tree structure. It significantly reduces the utilization of random access memory (RAM) in computer. The use of the FMM for solving 2D and 3D boundary value problems (BVPs) is well-documented [4] [5] [6] , also in application to the boundary element method (BEM) [7] .

Mentioned above BEM together with the finite element method (FEM) [8] are well-established (might be called classic or conventional) methods of modelling and solving boundary problems. However, new approaches are also being developed to eliminate the disadvantages of classical methods (the time-consuming discretization of a boundary or a domain resulting in a large number of finite or boundary elements). That group includes, among others, meshless methods [9] and still being developed the parametric integral equations system (PIES) [10] .

The authors of this paper still working on development and application of the PIES in modelling and solving BVPs. The PIES has been used to solve potential [11] , elasticity [12] or acoustics problems [13] . The PIES includes in its mathematical formalism the shape of the boundary of considered problem [10] , therefore it does not require discretization of the domain or the boundary, contrary to element methods. The shape of the boundary is directly included into the PIES kernels using well-known functions from computer graphics -curves for 2D and surface patches for 3D problems. The accuracy of solutions can be improved by changing the number of collocation points only, without any interference in the shape of the modelled boundary. The efficiency of modelling and high accuracy of solutions obtained using the PIES has been confirmed in previous studies (e.g. [11] [12] [13] ). The authors of this paper also proposed extensions of the PIES method for uncertainly defined [14, 15] or transient [16] problems.

Conventional PIES, similarly to the BEM, produces non-symmetrical dense matrices using O(N 2 ) operations and to solve the problem using direct solver, it needs another O(N 3 ) operations (where N is the number of system equations). Therefore, solving engineering problems with complex shapes requires a lot of RAM and time-consuming computations. Application of OpenMP [17] or acceleration of solving the PIES using CUDA [18, 19] allows for a significant reduction of time of computations. Unfortunately, the problem of limited resources of RAM in a personal computer (PC) still exists. Therefore, mentioned above parallelization techniques do not allow for efficient solving of problems on a PC.

The authors of this paper presented a way of including the FMM, based on a binary tree, into modified PIES in [20] . Application of the FMM increased the difficulty of implementation of so-called the fast PIES. However, the proposed approach gives accurate solutions in a short time for examples with a quite simple shape of the boundary. Some assumptions of the FMM may result in obtaining incorrect solutions, especially for complex shapes of a boundary. Our research has shown the need to modify the binary tree and to change the FMM algorithm to consider assumptions mentioned above for complex shapes of a boundary.

The main goal of this paper is to present the modified binary tree in the FMM used to accelerate numerical calculations in the PIES and to reduce utilization of RAM for BVPs with complex shapes of a boundary. The efficiency and accuracy of the proposed fast PIES are tested on 2D potential BVPs.

Conventional PIES for 2D potential problems is presented by the following formula [10]: where: l = 1, 2, ..., n, s l−1 ≤ s ≤ s l , s j−1 ≤ s ≤ s j , s l−1 and s j−1 correspond to the beginning of l-th and j-th segment, while s l and s j to their ends, J j (s) is the Jacobian, n is the number of parametric segments that creates boundary of domain in parametric reference system s and s (presented in Fig. 1 ). 

where: S (1) = S Boundary functions u j (s) and p j (s) in (1) are approximated by the following series:

where u The pseudospectral method was applied to solve the PIES (1). Therefore, the PIES is transformed into the system of algebraic equations Ax = b. The oldest implementation of the PIES uses Gaussian elimination with pivot to solve the system. Newer version, as well as the fast PIES, uses very fast and popular iterative solver based on the generalized minimal residual method (GMRES) [21] . Application of the FMM allows reducing computational time of matrixvector multiplication (frequently used by GMRES solver) from order O(N 2 ) to O(N ). There is also no need to store in the RAM entire matrix [6] .

The FMM in the fast PIES is composed of two main steps: the upward and the downward pass (a full description is presented in [20] ). The upward pass uses procedures of kernels expansion and moment-to-moment translation, while in the downward pass moment-to-local and local-to-local translations are applied.

The most important information about these procedures are described below.

Moments. The first procedure of the fast PIES is an expansion of kernels in Taylor series. The direct inclusion of the FMM [1] into the PIES is problematic due to the calculation of subsequent derivatives of kernels (2) for the Taylor series approximation. In the paper [20] , the authors presented a modification of the PIES kernels by complex analysis. The PIES with modified kernels have the following form [20] :

where: -is the real part of complex number, the parametric functions, which describe the shape of the boundary in the PIES, can be defined in complex notation as τ = S (c)

j (s), (c) -means complex variable and an indeterminate (the imaginary unit) i = √ −1. The complex form of kernels is as follows [20] :

where n (c) = n (1) + in (2) -the complex notation of normal vector to the curve created segment j.

We assume that the point s c (corresponding to the complex point τ c ) is close to the observation point s ob and the point s el (corresponding to the complex point τ el ) is close to the collocation point s col (presented in Fig. 2 ). If |s ob −s c | |s col − s c | and |τ ob − τ c | |τ col − τ c |, then kernels (5) can be expanded about the point τ c using the Taylor series expansion. Therefore, we obtained the multipole moments M k (τ c ) and N k (τ c ) [20] : (6) and the following form of approximated Eq. (4) [20] :

where:

Moments are calculated once only and they are independent of τ (that is also from s). Fig. 2 ) without reuse of integration. For this purpose, moment-to-moment translation was used [20] :

Moment-to-Local Translation (M2L) and Local Expansion. The next step of the fast PIES is local expansion around the point s el (corresponding to the complex point τ el ) (presented in Fig. 2 ). If |s col − s el | |s c − s el | and |τ col − τ el | |τ c − τ el |, then U k and P k in (7) can be expanded about the point τ el using the Taylor series expansion, hence [20] :

where:

This procedure allows the transformation of moments collected at a point s c (τ c ) to a local expansion point s el (τ el ) and is called moment-to-local translation (M2L) [20] .

During modelling complex shapes of the boundary, fulfilment of the condition |τ col − τ el | |τ c − τ el | might be not possible to meet for some cells. The results obtained in this case are subject to large errors. The authors of this paper proposed the way of elimination such cases -in the FMM algorithm they are considered as neighbouring cells (a more detailed description of the neighbourhood of the cell is described in Sect. 2.3).

Similarly to the M2M, moments at a point s el (τ el ) can be efficiently recalculated for a nearby point s el (corresponding to the complex point τ el ) (presented in Fig. 2 ). For this purpose, a transformation called local-to-local (L2L) translation (similar to M2M) is used [20] :

In classic FMM, the binary tree for 1D, quad-tree for 2D and octa-tree for 3D problems are applied. In this paper, 2D problems are discussed, hence in classic FMM implementation quad-tree presented in Fig. 3 [6] is used. A square surrounding the entire domain of the problem (Fig. 3a) is called level 0 cell. This cell is the parent of four cells of level 1 (so-called children cells) obtained as a result of dividing the parent cell into four identical squares. Only the cell crossing the boundary of the problem is considered to be a child cell. The parent cell of level l (l ≥ 0) is divided into children cells of a level l + 1. The division is performed until a predetermined number of elements are inside a cell (in the example in Fig. 3a -max. 2 elements whose centre is marked as a collocation point) or the assumed maximum level l has been reached. A cell without a child is called a leaf. The quad-tree presented in Fig. 3b is obtained in the described way. The presented method of the FMM tree construction is applied, among others, in the FM-BEM [6] .

In the PIES applied for the 2D problems, it is possible to implement a suitably modified and improved binary tree. The improvement is related to the way of defining problems in the PIES -the boundary of the problem is defined in the 1D parametric reference system (presented in Fig. 1 ). It is a different way of modelling the shape of a boundary than in the BEM. Also, the number of segments describing the boundary in the PIES is smaller than the number of elements in the BEM (see Fig. 4 ). It is connected with the way of defining the shape of the boundary in the PIES [10] . The proposed binary tree structure assumes that the boundary is described in the 1D reference system. A level 0 cell covers the entire boundary of the problem (presented in Fig. 4a ). This cell is the parent of two children cells of level 1 obtained as a result of dividing the parent cell into two identical segments. The parent cell of level l (l ≥ 0) is divided into children cells of a level l + 1. The division is performed until a predetermined number of segments are inside a cell (in the example in Fig. 4a -max. 2 segments) or the assumed maximum level l has been reached. In the classic FMM, the binary tree for 1D problems has the form presented in Fig. 4b . However, in the PIES 2D problem is reduced to the 1D parametric reference system, hence the first and the last boundary segment are very close to each other. Therefore, we propose to close the binary tree in the form presented in Fig. 4c . For all levels, the first and the last cells are treated as adjacent.

The algorithm of the PIES with modified kernels proceeds in several steps (flow chart presented in Fig. 5 ). The first step after initialization of the algorithm is to determine the structure of the modified binary tree. Then right-hand sided vector b is computed and GMRES is called to find solution of Ax = b. These procedures uses the modified fast multipole algorithm presented in Algorithm 1. At the end of algorithm results are presented on the screen and written to the output file.

The first run of the fast multipole procedure is used for calculating the righthand sided vector b. The fast multipole procedure is composed of two steps: the upward pass and the downward pass. At the upward pass, all moments in leaves are calculated (line 4 in Algorithm 1) for the level lev. Then, tracing the tree structure upward, moments in all parent cells are calculated up to level 2 using M2M (line 6).

In the downward pass, we trace the tree structure downward, and previously calculated moments are used in calculations. First of all, we should remind cells neighbourhood to clarify this step [20] . In modelling complex shapes of the boundary, we should modify the FMM algorithm described in [20] due to the possibility of too small distance between cells with τ c and τ el points, hence the assumption of |τ col − τ el | |τ c − τ el | used for the M2L translation is not fulfilled. The authors of this work propose to modify the algorithm by marking such cells as adjacent and not to enter them on interaction list (calculations are performed as in the case of adjacent cells). Including this modification to the fast PIES, solutions with the same accuracy as in conventional PIES are obtained.

Starting from level 2 and tracing the tree structure downward to all leaves coefficients of local expansion (the line 24 in the Algorithm 1) are computed. Coefficients at cell K − th at level i are computed as the sum of two elements:

lev -the number of tree levels min ci -the lowest number of a cell on the level i max ci -the highest number of a cell on the level i //upward pass 1: for i ← lev to 2 do 2:

for icell ← min ci to max ci do 3:

if i == lev then 4:

multipole moments(icell) 5: To solve the system of algebraic equations Ax = b, iterative GMRES solver is used. The method requires the application of multiplication of the matrix A by the vector of unknowns x, therefore the solver can be directly integrated with the fast PIES. The FMM in GMRES is performed in the same way as for vector b.

In the paper [20] , the authors presented preliminary results of the fast PIES solutions. However, solving the problem presented in this paper by the fast PIES without modification of the binary tree gives unsatisfactory results. Obtained solutions were far from expected.

The problem of temperature distribution (modelled by Laplace's equation) in heat-sink is considered. The shape of the boundary and boundary conditions are shown in Fig. 6 . The boundary is composed of 716 linear segments. The same number of collocation points (from 4 to 8) is defined on each segment, and finally, we should solve the system of 2864 to 5728 algebraic equations. The problem is solved by conventional and fast PIES. PC based on Intel Core i5-4590S with 8 GB RAM and g++ 7.4.0 compiler with -O2 optimization on 64-bit Linux operation system (Ubuntu, kernel 5.0.0) and LAPACK 3.10.3 library [22] is used during tests. First of all, a comparison of CPU time and RAM utilization between a different number of the modified tree levels in the fast PIES is performed to find the optimal value of tree levels. Taylor series composed of 25 terms approximated the fast PIES kernels. The value of tolerance (convergence criterion) of the GMRES was equal to 10 −8 .

As can be seen from Fig. 7 , both CPU time and memory utilization decrease with the growing number of tree levels reaching the minimum value for a tree with 5-6 levels regardless of the number of collocation points. In our studies, the tree with 6 levels is adopted.

The research involved a comparison of the speed and RAM utilization between conventional, the fast PIES and the fast multipole BEM (application from [6] ). In the fast multipole BEM two meshes are used: 2864 and 5728 elements. The value of tolerance (convergence criterion) of the GMRES and the number of terms in Taylor series are the same as previous, i.e. 10 −8 and 25 respectively. Solutions are presented for two versions of conventional PIES As can be seen from Table 1 , the fast PIES is significantly faster than both versions of conventional PIES and about 2 times faster than the fast multipole BEM. The fast PIES also needs up to 4 times less RAM than the PIES i and about 8 times less than the PIES d . Obtained solutions are practically the same as in conventional versions of the PIES -MSE between the PIES d and the fast PIES does not exceed 2.31 · 10 −10 . The MSE between the fast PIES and the fast multipole BEM has higher value. However, our previous studies proved that the PIES is more accurate than the BEM.

Graphical comparison of CPU time and RAM utilization between all the PIES methods is presented in Fig. 8. 

The paper presents the fast PIES based on the FMM with the modified binary tree in solving potential BVPs with complex shapes. The proposed modified binary tree is built based on the 1D parametric reference system. Such modification of the tree in the fast PIES allows obtaining very accurate solutions for engineering problems with complex shapes on standard PC in a short computational time. It also allows solving problems in which the condition of the appropriate distance between the centres of the leaves is not fulfilled.

The numerical test shows a reduction of the computation time and RAM utilization of the fast PIES compared to conventional ones as well as the fast multipole BEM. The speed-up of computations between the fast and conventional PIES increases with the size of solving problem, while the accuracy of solutions is almost the same.

Rapid solution of integral equations of classical potential theory

A fast algorithm for particle simulations

The Rapid Evaluation of Potential Fields in Particle Systems

FMM-accelerated source-model technique for manyscatterer problems

Fast multipole accelerated boundary integral equation methods

The fast multipole boundary element method for potential problems: a tutorial

Boundary Element Techniques, Theory and Applications in Engineering

The Finite Element Method

Adaptive methodology for meshless finite point method

Hermite curves in the modification of integral equations for potential boundary-value problems

Modeling the shape of boundary using NURBS curves directly in modified boundary integral equations for Laplace's equation

Algorithms of identification of multiconnected boundary geometry and material parameters in problems described by Navier-Lamé equation using the PIES

Triangular Bézier patches in modelling smooth boundary surface in exterior Helmholtz problems solved by PIES

The influence of interval arithmetic on the shape of uncertainly defined domains modelled by closed curves

Modification of interval arithmetic for modelling and solving uncertainly defined problems by interval parametric integral equations system

Parametric integral equations systems in 2D transient heat conduction analysis

OpenMP for 3D potential boundary value problems solved by PIES

Application of CUDA for acceleration of calculations in boundary value problems solving using PIES

Acceleration of integration in parametric integral equations system using CUDA

The fast parametric integral equations system in an acceleration of solving polygonal potential boundary value problems

GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems

Acknowledgments. The scientific work is founded by resources of Ministry of Science and Higher Education for research, granted to the Institute of Informatics, University of Bialystok.