Available online at www.sciencedirect.com
Procedia Computer Science 18 (2013) 1594 – 1603
2013 International Conference on Computational Science
Graph Grammar Based Direct Solver for hp-adaptive Finite Element Method with Point Singularities Arkadiusz Szymczaka a
b
, Piotr Gurgula,
a,
*
AGH University of Science and Technology, Krakow, Poland b Jagiellonian University, Krakow, Poland
Abstract In this paper we present a graph grammar based direct solver algorithm for hp-adaptive finite element method simulations with point singularities. The solver algorithm is obtained by representing computational mesh as a graph and prescribing the solver algorithm by graph grammar productions. Classical direct solvers deliver O(Np 4+N1.5) computational cost for regular 2D grids, and O(Np6+N2) for regular 3D grids, where N denotes number of degrees of freedom and p denotes the polynomial order of approximation. The solver presented in this paper delivers linear computational cost for uniform polynomial order of approximation p. For non-uniform polynomial order the computational cost is almost linear. 3 The Authors. Published by Elsevier B.V. © 2013 peer review under responsibility of the organizers of the International Conference on Computational Selection and and/or peer-review under responsibility of the organizers of 2013 the 2013 International Conference on Computational Science Science Keywords: finite element method, graph grammars, direct solvers, digital material representation.
1. Introduction The direct solver is a core part of the solution for several challenging engineering problems solved by means of the Finite Element Method (FEM) [1-3]. Exemplary applications involve generation of acoustic waves over the model of the human head [4] or borehole resistivity simulations [5]. The process of solving finite element engineering problems starts with generation of the mesh, which describes the geometry of the computational problem. Next, the physical phenomena governing the problem is described by a Partial Differential Equation (PDE) with boundary and / or initial conditions. Then, the PDE is discretized into a system of linear equations using FEM. At this point, the solver algorithm is executed in order to provide the solution to the system of
* Corresponding author. Tel.: +48 -12- 328-33-14 E-mail address:
[email protected].
1877-0509 © 2013 The Authors. Published by Elsevier B.V. Selection and peer review under responsibility of the organizers of the 2013 International Conference on Computational Science doi:10.1016/j.procs.2013.05.327
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
linear equations. Challenging engineering problems generate huge linear systems with several million unknowns, and the solver algorithm is the computationally most expensive part of the process. The multi-frontal solver is the state-of-the art algorithm for solving linear systems of equations [6, 7] using a direct method. The multi-frontal algorithm constructs an assembly tree based on the analysis of the connectivity data or the geometry of the computational mesh. Finite elements are joined into pairs and fully assembled unknowns are eliminated within frontal matrices associated to multiple branches of the tree. The process is repeated until the root of the assembly tree is reached. Finally, the common interface problem is solved and partial backward substitutions are recursively performed on the assembly tree. The computational cost of the sequential (single core) multi-frontal solver algorithm for two and three dimensional problems, requiring higher polynomial order of approximation p, solved over a regular mesh is of the order of O(Np4+N1.5) in 2D and O(Np6+N2) in 3D, where N denotes mesh size (number of unknowns, called the degrees of freedom) and the memory usage is correspondingly of the order of O(Np2+Nlog(N/p)) in 2D and O(4/3Np2) in 3D, as it has been recently shown in [8, 9]. The self-adaptive hp-FEM algorithm was developed by the group of professor Leszek Demkowicz [2-3]. It generates a sequence of meshes delivering exponential convergence of the numerical error with respect to the mesh size (the number of unknowns in the system of linear equations). This is obtained by performing a sequence of h-, p- and hp-refinements starting from a regular, initial mesh. The h-refinement breaks a finite element into smaller elements, whereas the p-refinement increases the polynomial order of approximation and the hp-refinement applies both approaches simultaneously. The exponential convergence of the hp-FEM can reduce dramatically the size of the computational mesh, since we can obtain better accuracy on small, yet highly hp-refined meshes. For example [2], the 0.01% accuracy solution (measured in H 1 norm relative error) of a two dimensional heat transfer problem with anisotropic material data requires the solution of a linear system with less than 2000 unknowns if we use the hp-FEM algorithm. On the other hand, classical FEM with regular meshes requires several million unknowns to provide the same 0.01% accuracy. This difference grows when we switch to challenging three dimensional problems [3]. In this paper we propose a graph grammar based direct solver, delivering both linear computational cost and memory usage for computational problems refined towards point singularities when the polynomial order of approximation is uniform. For the non-uniform polynomial order case, the solver cost is bounded by the linear cost for the larger grid with uniform polynomial order of approximation, being the maximum of the order used over the mesh. In other words, we present a graph grammar based solver algorithm that utilizes the special structure of the computational grid refined towards a point singularity. This paper generalizes the results for two dimensional grids with h-adaptivity only, summarized in [10]. We present also new results concerning hpadaptivity for three dimensional grids. 2. L-shape domain model problem [11,12], to test The Lthe convergence of the p- and hp-adaptive algorithms. The problem consists in solving the temperature distribution over the L-shape domain, presented in Fig. 1 with fixed zero temperature in the internal part of the boundary, and the Neumann boundary condition prescribing the heat transfer on the external boundary. There is a single singularity in the center point of the domain (the gradient of temperature goes to infinity, compare Fig. 2), so the accurate numerical solution requires a sequence of adaptations towards the central point. The problem can be summarized as follows: Find the temperature distribution u : R2
x
ux
R such that
1595
1596
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603 2 i 1
2
u xi2
0 in
(a)
(1)
(b)
Fig. 1. (a) the L-shape domain; (b) detailed view at the
u
at the central part of the mesh
zoom x 1
zoom x 10
zoom x 100
zoom x 10000
zoom x 100000
zoom x 1000000
zoom x 10000000
zoom x 100000000
solution
Fig. 2. The optimal mesh generated by the self-adaptive hp-FEM algorithm for the L-shape domain model problem. Different colors denote different polynomial orders of approximations. The final optimal mesh delivers solution with relative error below 0,001% error. The resulting solution is a temperature distribution. The bottom panels present color denoting different polynomial orders of approximation and the scale for the resulting temperature scalar field.
1597
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
with boundary conditions u 0 on
(2)
D
u g on N n with n being the unit normal outward to 2 r3
(3) vector, and
2 sin 3
(4) 2 which is defined in the radial system of coordinates with the origin point O presented in Fig 1. The formula (4) is based on the exact solution to the L-shape problem. The L-shape domain problem has been solved with both self-adaptive hp-finite element method and h-adaptive finite element method with uniform polynomial order of approximation p. For the hp-adaptive case summarized in Fig. 2, the solver delivers exponential convergence of the numerical error with respect to the mesh size. For the h-adaptive case, summarized in Fig. 3, the solver delivers algebraic convergence rate. The convergence curves are summarized in Fig. 4.
g r,
Fig. 3. Sequence of meshes generated by the h-adaptive FEM algorithm for the L shape domain model problem.
Fig. 4. Convergence history for a sequence of meshes generated by uniform h adaptation, uniform p adaptation, non-uniform h adaptation and non-uniform hp adaptation.
1598
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
3. Fichera problem The Fichera problem is the generalization of the L-shape domain problem into the three dimensions. It can be summarized in the following way: 3 x u x R over the domain presented in Fig. 3 such that Find the temperature distribution u : R 3 2u 0 in (5) 2 i 1 xi with boundary conditions u 0 on D (6)
u g on N n with n being the unit normal outward to
(7) vector, and g is the exact solution of the L shape problem.
Fig. 5. Domain for the Fichera problem
Fig. 6. The sequence of meshes generated by the self-adaptive hp-FEM for the Fichera problem. Different colors denote different polynomial orders of approximations presented in Fig. 2.
1599
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
Like the L shape model problem, the Fichera problem has also been solved with both self-adaptive hp-finite element method, and with h-adaptive finite element method with uniform polynomial order of approximation p. For the hp-adaptive case summarized in Fig. 6 the solver achieved the exponential convergence rate. For the hadaptive case, the solver delivers algebraic convergence rate only. The convergence curves are summarized in Table 1.
Fig. 7. Sequence of meshes generated by the h-adaptive FEM algorithm for the Fichera shape domain model problem. Table 1.(a) converge rate for h-adaptive FEM for uniform p=5 refined towards point singularity; (b) convergence rate for hp-adaptive FEM.
(a) hadaptive FEM
N
Error
Coarse mesh
Fine mesh
Error
1206
5.06
665
13897
9.75
8261
3.18
846
17019
6.18
13726
2.46
1093
20849
4.58
19191
2.23
1577
27775
3.55
35586
2.15
2247
36771
2.91
(b) hpadaptive FEM
4. Graph grammar based solvers
Fig. 8. Graph representation of the 1/3 of the computational mesh h-refined towards point singularity in the L-shape domain problem.
1600
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
The locally refined computational mesh can be represented as a graph, and the solver can be expressed by graph grammar productions [13,14]. The two dimensional representation shown in Fig. 8 follows the graph grammar summarized in [15], the three dimensional representation depicted in Fig. 9 follows the graph grammar summarized in [16-17].
Fig. 9. Graph representation of the 1/7 of the mesh, h-refined towards point singularity in the Fichera domain problem.
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
Fig. 10. Sequence of graph grammar productions coloring the nodes for the solver algorithm in the L shape domain problem
Fig. 11. The sequence of graph grammar productions coloring the nodes for the solver algorithm in the Fichera problem.
Some sub-graphs of the graph representation of the mesh represent particular finite elements. The solver algorithm browses the graph representation of the mesh starting from bottom elements up to the root, level by level. It generates element frontal matrices produced as the result of discretization of the computational problem, and it merges the element matrices into the single frontal matrix. First, it identifies fully assembled nodes located at each level of the graph representation of the mesh, denoted by yellow color in Fig. 10 for 2D and Fig. 11 for 3D, eliminates them, and then it repeats the process on the next level. The number of nodes that are not fully assembled, denoted by dark blue color in Fig. 10 and Fig. 11, remains constant from one level to another. This pattern for elimination ensures that the size of a single frontal matrix involved in the solver algorithm remains constant. The order of browsing of sub-graphs representing particular elements, with generation of element frontal matrices is expressed by graph grammar productions, coloring the nodes of the graph representation of the mesh, as it is presented in Fig. 10 and 11. 5. Linear computational cost and memory usage of the solver algorithm Since the cost of each step (level of the elimination tree) is constant, the total cost of the algorithm is proportional to the number of levels, which is proportional to the number of unknowns (it is implied by grid construction). As a result, for a constant polynomial order of approximation, we obtain a solver algorithm with linear computational cost with respect to the number of unknowns (see red line in Fig. 12 and Fig. 13) [10]. For non-uniform distribution of polynomial orders of approximation as it happens for the adaptive hp-refinements, the computational cost is bound by the linear cost for uniform grid with maximum order being used, as it is
1601
1602
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
presented in Fig. 12 and Fig. 13. The computational cost of the solver algorithm for hp-adaptive is bound by the linear cost for uniform p=5 until the iteration number where the higher polynomial orders of approximation p > 5 are used.
Fig. 12. Computational cost of the graph grammar based solver algorithm over the coarse mesh for the L shape domain problem for hp adaptive and h adaptive uniform p=5 case.
Fig. 13. Computational cost of the graph grammar based solver algorithm over the coarse mesh for the Fichera problem for hp adaptive and h adaptive uniform p=5 case.
6. Conclusions In this paper we presented a graph grammar model for enforcing order of elimination of degrees of freedom over graph representations of two and three dimensional grids refined towards point singularities, resulting in linear computational cost of the solver algorithm. The algorithm has been tested on the two dimensional L shape domain problem as well as on the three-dimensional Fichera model problem. The algorithm delivers linear computational cost when we restrict ourselves to h-refinements, and close to linear computational cost when the h-adaptive solver is accompanied by p-refinements. This numerical results emphasize the substantial advantage of using non-uniform locally refined grids as opposed to uniform grids, delivering computational costs of the order of O(Np4+N1.5) in 2D and O(Np6+N2) in 3D. We also plan to investigate the prospective use of an agent-based model summarized in [18, 19] for implementation of the methodology described in this paper in heterogeneous environment. Acknowledgements The work presented in this paper is supported by grant of Polish National Science Center grant no. NN 519 447 739. The work of the PG was partly supported by The European Union by means of European Social Fund, PO KL Priority IV: Higher Education and Research, Activity 4.1: Improvement and Development
Arkadiusz Szymczak et al. / Procedia Computer Science 18 (2013) 1594 – 1603
of Didactic Potential of the University and Increasing Number of Students of the Faculties Crucial for the National Economy Based on Knowledge, Subactivity 4.1.1: Improvement of the Didactic Potential of the AGH University of Science and Technology ``Human Assets'', No. UDA POKL.04.01.01-00-367/08-00. References [1] Hughes TJR, 2000. Linear static and dynamic finite element analysis. Dover Publications. [2] Demkowicz L, 2006. Computing with Hp-Adaptive Finite Elements, Vol. 1: One and Two Dimensional Elliptic and Maxwell Problems. Chapman & Hall / CRC Press. [3] Demkowicz L, Kurtz J., Pardo D, Paszynski M, Rachowicz W, Zdunek A, 2007. Computing with Hp-Adaptive Finite Elements, Vol. 2: Frontiers: Three Dimensional Elliptic and Maxwell Problems with Applications. Chapman & Hall / CRC Press. [4] Demkowicz L, Gatto P., Kurtz J M, Rachowicz W E M, Hamilton M, Champlin C, Pardo D, 2010. Modeling of bone conduction of sound in the human head using hp-finite elements: code design and verification. Computer Methods in Applied Mechanics and Engineering; 200(21-22):1757 1773. M, Schaefer R, 2009. Graph grammar driven parallel partial differential equation solver. Concurrency and Computations, Practise and Experience; 22(9):1063-1097. [6] Duff IS, Reid JK, 1984. The multifrontal solution of unsymmetric sets of linear systems. SIAM Journal on Scientific and Statistical Computing; 5:633-641. [7] Geng P, Oden TJ, van de Geijn RA, 2006. A Parallel Multifrontal Algorithm and Its Implementation, Computer Methods in Applied Mechanics and Engineering, 149:289-301. [8] Collier N, Pardo D, Dalcin L M, Calo V, 2012. The cost of continuity: A study of the performance of isogeometric finite elements using direct solvers. Computer Methods in Applied Mechanics and Engineering; 212-216:353361. [9] Pardo D M, Collier N, Alvarez J, Dalcin L, Calo VM, 2012. A survey on direct solvers for Galerkin methods. SeMA Journal, Spanish Society for Applied Mathematics; 57:107-134 [10] Gurgul P M, Szymczak A, Calo V, Pardo D, 2012. A linear complexity direct solver for h-adaptive radical grids. Submitted to Computing and Informatics. I, Guo B, 1986. The hp-version of the finite element method, Part I: The basic approximation results. Computational Mechanics; 1: 21-41. I, Guo B, 1986. The hp-version of the finite element method, Part II: General results and applications. Computational Mechanics; 1:203-220. M, Pardo D, ska A, 2010. Parallel multi-frontal solver for p adaptive finite element modeling of multi-physics computational problems; 1(1): 48-54. [14] Szymczak A, Paszynski M, Pardo D, 2010. Graph grammar based Petri net controlled direct solver algorithm; Computer Science 11(2010): 65-79. M A, 2008. Graph transformations for modeling parallel hp-adaptive finite element method, Lecture Notes in Computer Science; 4967:1313-1322. A, Grabska E, Paszynski M, 2012. A Graph Grammar Model of the hp Adaptive Three dimensional Finite Element Method, Part I. Fundamenta Informaticae;114(2):149-182. A, Grabska E, Paszynski M, 2012. A Graph Grammar Model of the hp Adaptive Three dimensional Finite Element Method, Part II. Fundamenta Informaticae; 114(2):183-201. [18] Cetnarowicz K, Gruer P, Hilaire V, 2002. A formal specification of M-agent architecture, Proc. Multi-Agent Systems CEEMAS 2001, Berlin, Heidelberg, vol. 2296, pp. 62-72. [19] Cetnarowicz K, 2009. From algorithm to agent, Computational Science ICCS 2009, LNCS 5545 Springer Verlag, pp. 825-834.
1603