A suboptimal strategy for two dimensional applications of DG methods①

HUANG Rixin,TAN Yonghua,WU Baoyuan,LI Guangxi

(1.Xi'an Aerospace Propulsion Institute, Xi'an 710100, China;2.Academy of Aerospace Propulsion Technology, Xi'an 710100, China)

Abstract：Serendipity square elements of high degree up to 10 times are constructed as a strategy for two dimensional applications of nodal DG methods. By pre-imposing the Legendre-Gauss-Lobatto (LGL) quadrature points on the borders to keep good interpolation properties and maintain boundary conforming characteristics, the novel construction method evolves a global-like solution of a constrained nonlinear optimization problem to maximize the absolute value of the Vandermonde determinant of the point set. The final point set with certain symmetry property produces low Lebesgue constants which indicate low interpolation errors, and the constants fall among the ranges of those of different point sets already known in literature. Compared with the constant metric elements equipped with nodal points of compact pattern, only two additional points are introduced for the newly presented strategy and that is also why the strategy is called suboptimal. On the other hand, the new strategy has a much smaller scale of nodal points than the traditional tensor product points, therefore it remarkably saves computing and storage resources and is more suitable for application.

Key words：Serendipity element；discontinuous Galerkin method；Interpolation；Lebesgue constant；Genetic algorithm

0 Introduction

Thanks to the available compact nodal point sets, the nodal DG (discontinuous Galerkin) methods can be expressed in a very simple form, especially for elements of constant metric during coordinate transformation, such as triangle and tetrahedron elements. But for non-constant metric elements, such as quadrilateral element, things are different; because the existing compact point sets are not conforming on boundaries, so the simplicity of the DG method is destroyed. The boundary conforming point sets, such as the tensor product type, have an unacceptably larger points scale than the compact counterpart. Therefore, how to construct the point sets of both storage saving and boundary conforming is necessary.

(1)

(2)

(3)

Good point sets have already been known for square element. The Fekete points of the cube (including square) are analytically known as tensor product of LGL points[1]. Padua points[2-4]and points optimized[5]from Padua or tensorial points are AFPs. These point sets all have low Lebesgue constants. However, the tensor product points (TPPs) have a dimension of (n+1)d, so there will be too many nodes to be affordable in practical applications. The Padua-type AFPs have the most compact node pattern, but the node distribution on edges is not conforming, which blocks a wider use of boundary conforming methods. Besides, the nodes on edges are far away from the Fekete points, so they are unfavorable for surface integrations. To remedy the shortages of the previously known point sets, a new type of serendipity square elements based on Fekete points (SFPs) is developed.

1 SFPs for square

Elements of serendipity family are commonly used for structural mechanics in the continuous finite element methods. For degrees not beyond cubic, all indispensable nodes are placed only along the edges, which is more preferable to those of Lagrange family. It is therefore easy to develop a systematic algorithm of generating shape functions and very desirable if different degrees transit between elements. But immediately, a quartic serendipity element needs an interior node added in the center to fulfill the expansion completeness. And higher degrees require more internal nodes or nodeless variables to carry appropriate information of polynomials[6]. Unfortunately, systematic studies of the serendipity family are rarely seen in literature. According to the latest progression[7], the family has had a dimension-independent definition; nevertheless it is still not evident for the node patterns of higher-degree (> 4) cases.

For square elementI2whereI=[-1,1] is the unit interval, the complete polynomial space suggested by the Pascal triangle is

Pn(I2)=span{xpyp}p,q∈γ′

γ={(p,q)|0≤p,q≤n,p+q≤n}

(4)

which is the ideal space for two-dimensional case. To keep the boundaries conforming, the polynomials produced by the boundaries in the full tensor space should be retained, consequently two polynomialsxnyandxynare introduced by the boundary points[8]and cannot be removed for requirements of completeness. Therefore, the complete polynomial space for the serendipity square are taken to be the form of

Sn(I2)=Pn(I2)?span[xny,xyn]

(5)

and the dimension is

(6)

where (·) is the binomial. For convenience, the serendipity space is decomposed into interior and boundary contribution, denoted by the superscriptsiandbrespectively

(7)

The decomposition does not require that the boundary modes have zero support in the interior region and the interior modes have zero support on the boundary.γof the interior space can be written as

γi={(p,q)|2≤p,q≤n-2,p+q≤n}

(8)

and the dimension is (n-2)(n-3)/2. The complementary space is the boundary space.

(9)

2 Global optimization

A subroutine of open source named Pakaia[10]in Fortran 77 from high altitude observatory (HAO) of national center for atmospheric research is used to implement optimizations. The original codes are extended to achieve double precision. And the random number generator is replaced by the Fortran intrinsic function. Bigger maximum mutation rate is adopted instead of the default value 0.25. Pakaia provides a mechanism of dynamically controlling the mutation rate by monitoring the difference in fitness between the best and median in the population. Making use of the roulette wheel algorithm, the selection is rank-based and stochastic. The encoding within Pikaia is based on a decimal alphabet made of the 10 simple integers (0 through 9), and in practical optimizations, the number of significant digits ranges from 5 to 7, depending on the problem scale. Concerning the inequality constraint, a penalty set in the fitness function can work well, but actually nothing to be done in this application, because all the nodes fall into the inner region of the square naturally. It is found that even for the same problem, the implementations produce different 'optimal' results among which there are small discrepancies. Therefore only the best one is accepted.

3 Initial individual with good fitness

In the original GAs, all of the individuals in the initial population are created randomly, whereas in this paper, the initial population is slightly altered by adding an excellent individual in it. With this little trick, individuals in all of the evolution will mature within fewer generations. A method of producing the initial individual with excellent genes is presented. The essence of this method is extracting the points that contribute most to the determinant of the Vandermonde matrix from the interior points of a Padua-type point set as interior points of SFPs. Letrbe the Padua-Jacobi points

(10)

whereJdenotes the set of Jacobi-Gauss-Lobatto points inI, and

(11)

(12)

Denote the diagonal submatrices byVbandVirespectively, it is obvious that

(13)

As the edge nodes are fixed as LGL points, to maximize the determinantVnis equivalent to maximize the determinant ofVi. A greedy algorithm in literature[11] is used here for selecting N-4n columns with maximal associated volume, and this algorithm can be implemented by the QR factorization with column pivoting. Before the interior points are selected,VnandViare underdetermined, so one can resort to the LAPACK routine dgeqp3 to compute the QR factorization. Throughout this paper, the Chebyshev bases are adopted to produce well-conditioning Vandermonde matrices.

4 Numerical results

The constructed serendipity square elements of degree from 5 up to 10 are listed in Fig.1. The circles in Fig.1 illustrate the patterns of the optimal SFPs. The dashed lines indicate a symmetric distribution of the interior points along the diagonal line. And an arrow-like shape divides the points into simple stratified groups whose peaks lie on the diagonal line. Inspired by this diagonally symmetric pattern, we can conjecture that the SFPs may be not unique and different arrow directions suggest other possible sets of SFPs. Two meaningful parameters, the number of arrows (also the number of the points lying on the diagonal line)Narwand the number of points on each arrowNparw, lists in Table 1 for different degrees. There are similarities between any two neighboring degrees. Either the neighbors have the same number of arrows, or have the same number of arrow points. And it is easy to conclude that

(14)

For degrees greater than 10, the diagonal symmetry property of SFPs can help us to reduce the dimension of the variable vector considerably, which will be examined in our future work.

(a) n=5 (b) n=6 (c) n=7

(d) n=8 (e) n=9 (f) n=10

Table 1 Relevant parameters for the interior SFPs of different degrees

Table 2 Quality metrics of some interpolation sets in square element

The metrics in Table 2 suggest the SFPs have a moderate quality among the given interpolation sets. For square element, the best interpolation set is the tensor product points, whose determinant absolute values are extremely large, while spectral condition numbers and Lebesgue constants are the lowest. The Lebesgue constants of the Lebesgue points are slightly greater than those of the tensorial points, but the Lebesgue points possess the most compact form. The absolute values of the determinants of Lebesgue points and the greedy points have the same order of magnitude which is two orders lower than the SFPs, however their Lebesgue constants have big discrepancies. This case reminds us that higher determinant absolute values do not guarantee lower Lebesgue constants, other factors such as the node patterns may need to be considered. Though being as compact as the Lebesgue points, the greedy points have the worst quality, so they are always improved by further iterations or as a start guess of other optimization methods.

In Fig.2, we have plotted the Lebesgue constants of the given interpolation sets. It manifests a linear increase of the Lebesgue constants with degree for the SFPs. For the tensorial points and the Lebesgue points, it has been reported as a sublinear growing,Λn=O[log2(n)][1,6]. The greedy points show a super linear growing tendency.

Fig.2 Lebesgue constants of different point sets

5 Conclusions

We have constructed the serendipity square elements of high degree from 5 up to 10 as a suboptimal strategy in DG applications. These elements have node patterns as SFPs whose interior points are optimal for maximizing the absolute value of the Vandermonde matrix in combination with the edge nodes pre-fixed as LGL points. We have used GA to implement a global optimization of the SFPs, and presented a new method to generate an initial individual with good fitness to evolve the generations faster. Though we are not sure whether the SFPs is the global optimization, by increasing the mutation rate, the optimization solution has more global meaning mathematically. Practical computations have shown that the newly presented initializing strategy can accelerate the evolution convergence.

Low condition numbers and Lebesgue constants suggest good inversion accuracy and interpolation stability of the constructed serendipity elements. Compared to other good point sets, these serendipity elements are both compact and conforming, therefore they are very promising in the applications of the finite element methods, especially for the nodal DG methods. The SFPs show a diagonal symmetric pattern, and this pattern may guide us to explore element of higher degree and higher dimension. Considering the triangle element is still the most optimal strategy, the computational domains should be triangulated as more as possible. While on the regions of boundary layer, one can use the serendipity elements as a suboptimal candidate to mimic the behaviors of the low Reynold viscous flows. Further validations of this strategy will be seen in our future work.