A spectral element method for meshes with skinny elements
AA SPECTRAL ELEMENT METHOD FOR MESHES WITH SKINNYELEMENTS
AARON YEISERMIT, CAMBRIDGE, MA 02139, US ( [email protected] )ADVISOR: ALEX TOWNSENDDEPARTMENT OF MATHEMATICS, CORNELL UNIVERSITY, ITHACA, NY 14853.
Abstract.
When numerically solving partial differential equations (PDEs), the first step isoften to discretize the geometry using a mesh and to solve a corresponding discretization of thePDE. Standard finite and spectral element methods require that the underlying mesh has no skinnyelements for numerical stability. Here, we develop a novel spectral element method that is numericallystable on meshes that contain skinny elements, while also allowing for high degree polynomials oneach element. Our method is particularly useful for PDEs for which anisotropic mesh elements arebeneficial and we demonstrate it with a Navier–Stokes simulation. Code for our method can be foundat https://github.com/ay2718/spectral-pde-solver.
Key words. spectral element method, skinny elements, ultraspherical spectral method
AMS subject classifications.
1. Introduction.
A central idea for numerically solving partial differential equa-tions (PDEs) on complex geometries is to represent the domain by a mesh, i.e., acollection of nonoverlapping triangles, quadrilaterals, or polygons. The mesh is thenused by a spectral element method to represent the solution on each element by ahigh degree polynomial [10, 15]. Engineering experience suggests that meshes con-taining skinny or skewed elements are undesirable and should be avoided, if possi-ble [13]. However, from a mathematical point-of-view, long and thin elements canbetter resolve anisotropic solutions [13] that appear in computational fluid dynamicsand problems involving boundary layers.In this paper, we develop a spectral element method that is numerically sta-ble on meshes containing elements with large aspect ratios and also allows for highdegree polynomials on each element. In particular, the numerical stability of ourmethod is observed to be independent of the properties of the mesh and polyno-mial degrees used. This has several computational advantages: (1) A practical hp -refinement strategy does not need to take into consideration the quality of the finalmesh, (2) Anisotropic solutions can be represented with fewer degrees of freedom, and(3) Low-quality meshes can be used without concern for numerical stability.Over the last two decades, there has been considerable research on mesh genera-tion algorithms that avoid unnecessarily skinny elements [11, 14]. While these meshgeneration algorithms are robust and effective, they can be the dominating computa-tional cost of numerically solving a PDE. By allowing low-quality meshes, our methodalleviates this current computational burden on mesh generation algorithms. This isespecially beneficial for time-dependent PDE in which an adaptive mesh is desired,for example, to resolve the propagation of a moving wave front.There are a handful of mesh generation algorithms that have been recently de-signed for anisotropic solutions that generate meshes with skinny elements [12]. It istherefore timely to have a general-purpose element method that allows for an accu-rate solution on such meshes. For time-dependent PDEs, such as the Navier–Stokesequations, skinny elements can relax the Courant–Friedrichs–Lewy time-step restric-tions, leading to significant computational savings in long-time incompressible fluidsimulations. a r X i v : . [ m a t h . A P ] M a r standard element method constructs a global discretization matrix with a po-tentially large condition number when given a mesh containing skinny elements. Thisis because the global discretization matrix is assembled element-by-element from a lo-cal discretization of the PDE on each element. Since the elements may have disparateareas, the resulting discretization matrix can be poorly scaled. One can improve thecondition number of the global discretization matrix by scaling the local discretiza-tion matrices by the element’s area. Since we are using a spectral element method, wecan go a little further and scale the local discretization matrices by the Jacobian of abilinear transformation that maps the element to a reference domain (see Section 3).The solutions of uniformly elliptic PDEs with Dirichlet boundary conditions de-fined on a convex skinny element are not sensitive to perturbations in the Dirichletdata or forcing term. Since we know it is theoretically possible to solve elliptic PDEson meshes with skinny elements, it is particularly disappointing to employ a numeri-cally unstable method (see Section 2).The structure of the paper is as follows: In Section 2, we define some useful termi-nology and provide the theory to show that certain elliptic PDEs on skinny elementsare well-conditioned. Section 3 describes the process of using our spectral elementmethod on a single square element, and Section 4 extends that process to convexquadrilateral domains. In Section 5, we define a method for joining together twoquadrilateral domains along an edge, and Section 6 extends that method to joiningan entire quadrilateral mesh. Key optimizations to our method are described in Sec-tion 7. Finally we present a Navier–Stokes solver using our spectral element methodin Section 8.
2. Uniformly elliptic PDEs on skinny elements.
It is helpful to reviewsome background material on skinniness of elements and the condition number ofPDEs on skinny elements. Given an element T of any convex 2D shape, let r in denotethe inradius , i.e., the radius of largest circle that can be inscribed in T , and r out denote the min-containment radius , i.e., the radius of smallest circle that contains T .We define the skinniness of T to be s ( T ) = r in r out . For finite element methods, it is widely considered that a skinny triangle with a largeangle, i.e., θ ≈ π , is more detrimental to numerical stability than those triangleswith only acute angles. This is based on the theoretical results on interpolationand discretization errors in [2]. The spectral element method that we develop isnumerically stable on meshes containing any type of skinny quadrilateral. Therefore,it is also numerically stable on meshes containing any kind of skinny triangle, becausea triangle can be constructed out of three non-degenerate quadrilaterals by dissectingalong the medians, as in Figure 3. Consider the second-order elliptic differential operatoroperating on domain Ω. L u = d (cid:88) i,j =1 a ij ( x ) ∂ u∂x i ∂x j + d (cid:88) i =1 b i ( x ) ∂u∂x i + c ( x ) u, here the coefficients a ij , b i , and c are continuous functions on Ω. We say that L isuniformly elliptic if for some constant θ > d (cid:88) i,j =1 a ij ( x ) ξ i ξ j ≥ θ d (cid:88) i =1 ξ i , x ∈ Ω , ξ ∈ R d . Lemma
Consider the uniformly elliptic partial differential operator L u = n (cid:88) i,j =1 a ij u x i x j + n (cid:88) i =1 b i u x i , with coefficients a ij and b i continuous. Let L u = f on domain Ω with boundary ∂ Ω .If f ≥ on Ω , then the maximum value of u occurs on ∂ Ω . If f ≤ on Ω , then theminimum value of u occurs on ∂ Ω .Proof. For a proof, see Theorem 1 in Section 6.5.1 of [3].
Lemma
Consider the uniformly elliptic partial differential operator L u = n (cid:88) i,j =1 a ij u x i x j + n (cid:88) i =1 b i u x i + cu, with coefficients a ij , b i , and c continuous. Let L u = f on domain Ω with boundary ∂ Ω , and let c < on Ω . If f ≥ on Ω , then the maximum positive value of u occurson ∂ Ω . If f ≤ on Ω , then the minimum negative value of u occurs on ∂ Ω .Proof. For a proof, see Theorem 2 in Section 6.5.1 of [3].Theorem 2.3 uses Lemmas 2.1 and 2.2 to show that perturbations to the boundaryconditions of a uniformly elliptic PDE with c ≤ Theorem
Let L u = (cid:80) ni,j =1 a ij u x i x j + (cid:80) ni =1 b i u x i + cu , with c ≤ . Let u bethe solution to L u = f on domain Ω with boundary ∂ Ω , and u ( ∂ Ω) = g . If v is thesolution to L u = f with boundary conditions u ( ∂ Ω) = g + (cid:15) , then max | v − u | ≤ max | (cid:15) | .Proof. Let us write v = u + δ . Since L is a linear operator, we know that L v = L u + L uδ = f . Therefore, L δ = 0, with δ ( ∂ Ω) = (cid:15) . By Lemma 2.2, we knowthat if δ achieves a nonnegative maximum on the interior of Ω, then δ is constant.Similarly, if δ achieves a nonpositive minimum on the interior of Ω, then δ is constant.If max | δ | > max | (cid:15) | , then δ will achieve a nonpositive minimum or a nonnegativemaximum on the interior of Ω. However, Lemma 2.2 states that δ must then beconstant, so therefore max | δ | = max | (cid:15) | , hence a contradiction. Therefore, max | δ | ≤ max | (cid:15) | .Theorem 2.4 provides a weak bound for the effect of perturbations to the righthand side on uniformly elliptic PDEs with c ≤ Theorem
Let Ω ∈ R d be a bounded domain and L a second-order uniformlyelliptic differential operator with c ≤ . Let u be a weak solution of (cid:40) L u = f on Ω ,u = g on ∂ Ω , nd ˆ u be a weak solution of (cid:40) L ˆ u = f + ∆ f on Ω , ˆ u = g + ∆ g on ∂ Ω , where u, ˆ u ∈ C (Ω) ∩ C ( ∂ Ω) and ∂ Ω is the closure of Ω . Then, sup x ∈ Ω | u ( x ) − ˆ u ( x ) | ≤ sup x ∈ ∂ Ω | ∆ g ( x ) | + (cid:16) e ( β +1)diag(Ω) − (cid:17) sup x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12) ∆ f ( x ) θ ( x ) (cid:12)(cid:12)(cid:12)(cid:12) , where β = max di =1 sup x ∈ Ω | b i ( x ) /θ ( x ) | .Proof. Since u satisfies L u = f on Ω and u = g on ∂ Ω, ∆ u = ˆ u − u satisfies L (∆ u ) = ∆ f on Ω and ∆ u = ∆ g on ∂ Ω. One can apply the bound in Theorem 3.7of [4].Theorem 2.3 shows that solutions to uniformly elliptic PDEs with c ≤ Theorem
Let L u = u xx + u yy − k u , and let u be the solution to L u = f ondomain Ω with boundary ∂ Ω , and u | ∂ Ω = g . Let r be the radius of the smallest circlecompletely containing Ω . If s is the solution to L s = f + (cid:15) with boundary conditions s | ∂ Ω = g , then max | s − u | ≤ max | (cid:15) | r .Proof. First, consider the smallest disk ω with radius r that completely containsregion Ω. Let this disk have center ( a, b ). For ( x, y ) ∈ ω and c <
0, the solution t = c (cid:0) ( x − a ) + ( y − b ) − r (cid:1) is nonnegative, max t = − r c , and t xx + t yy = 4 c .Next, consider solution w such that w xx + w yy = 4 c for c < w = 0 on ∂ Ω.The solution w − t satisfies ( w − t ) xx + ( w − t ) yy = 0 and w − t ≤ ∂ Ω, since t ≥ x, y ) ∈ Ω. By Lemma 2.1, ( w − t ) ≤ x, y ) ∈ Ω. Therefore,max w ≤ − r c .Finally, we will consider the solution v such that v xx + v yy = (cid:15) with v = 0 on ∂ Ωand max | (cid:15) | = c . By Lemma 2.1, v is bounded above by w + satisfying ∇ w + = − c and w + = 0 on ∂ Ω, and v is bounded below by w − satisfying ∇ w − = c and w − = 0on ∂ Ω. We can write (cid:15) = (cid:15) + + (cid:15) − , where (cid:15) + = max( (cid:15), (cid:15) − = min( (cid:15), δ + be the solution to L δ + = (cid:15) − with zero Dirichlet boundary conditionsand let δ − be the solution to L δ − = (cid:15) + with the same boundary conditions. ByLemma 2.2, δ + ≥ δ − ≤ ∇ δ + = k δ + + (cid:15) − . Since δ + ≥
0, min( k δ + + (cid:15) − ) ≥ min( (cid:15) − ),so δ + is bounded above by the solution of ∇ v = (cid:15) − with zero Dirichlet boundaryconditions by Lemma 2.1. Similarly, δ − is bounded below by the solution of ∇ v = (cid:15) + with zero Dirichlet boundary conditions. Therefore, s = u + δ + + δ − . Thus,max | s − u | = max(max δ + , max − δ − ) ≤ max | (cid:15) | r .Theorem 2.3 proves that any change to the boundary conditions of the screenedPoisson equation will not be amplified in the solution on any domain. Theorem2.4 gives a weak bound for solving all second–order uniformly elliptic PDEs with c ≤
0, and Theorem 2.5 shows that the amplitude of a perturbation in the righthand side of the screened Poisson equation is bounded by a function of the radiusof the domain. Therefore, solving the screened Poisson equation on any domain is awell-posed problem. − − − − − − κ ( × ) 0.1597 0.3161 0.9632 1.0763 1.0832 1.0832 1.0832 Table 1
The condition number of our matrix for solving Poisson’s equation on a skinny quadrilateralwith vertices (0 , , (0 . , . . (cid:15) ) , (1 , , (1 , − . (cid:15) ) . As (cid:15) approaches 0, the condition number isbounded from above.
010 5 × -101
01 43.532.52-5 1.52 -210 5 × -101 × -15 Fig. 1 . The left graph shows a solution to Poisson’s equation on quadrilateral times aslong as it is wide. The right graph shows the error between the computed solution and the exactsolution. Even on skinny quadrilaterals, our method is still numerically stable. Our only limitationon the skinniness of mesh elements is roundoff error in actually defining the element. It seems unfortunate to have an unstable numerical algorithm for a well-conditionedproblem. Therefore, we ask ourselves the question: Can we create a numerically stablealgorithm for solving PDEs on meshes with skinny elements?
Our method represents any linear differentialoperator as a matrix operator. One method of assessing the numerical stability of themethod is to examine the condition number of the matrix. In order to demonstratethat the condition number is bounded, let us define a skinny triangle with coordinates(0 , , (1 , (cid:15) ) , (2 , − (cid:15) ). We can now define a skinny quadrilateral with vertices(0 , , (0 . , . . (cid:15) ) , (1 , , (1 , − . (cid:15) ). When we construct a differential operatormatrix on that skinny quadrilateral, we still need a preconditioner. In fact, we foundthat row scaling so that each row has a supremum norm of 1 is fairly close to optimal.Figure 1 shows the condition number κ of the normalized matrix for solving Poisson’sequation as (cid:15) approaches zero. Even as the quadrilateral becomes very skinny, thecondition number of the matrix is bounded from above. The process of creating thematrix is shown in sections 3 and 4.Figures 1 and 3 show that our method gives extremely accurate results even onmeshes with skinny triangles. Figure 2 shows that our method performs equally wellon meshes with or without skinny elements.
3. Spectral methods.
In order to develop a useful spectral method, one mustchoose the proper set of basis functions. Although wavelets and complex exponentialsare tempting possibilities, certain sets of polynomials promise simplicity and easycomputability. Although polynomial interpolation is numerically unstable using eq-uispaced points, it is numerically stable using Chebyshev points, with the n th set ofChebyshev points x j = cos (cid:16) n − − j ( n − π (cid:17) [18]. solution to Poisson’s equation on two el-ements with low aspect ratio Error between the exact solution and thecomputed solution ( < − ) A solution to Poisson’s equation on an ele-ment with high aspect ratio attached to anelement with low aspect ratio Error between the exact solution and thecomputed solution ( < − ) Fig. 2 . We can solve Poisson’s equation equally accurately on a mesh with and without skinnyelements. Each mesh is the union of the square [ − , × [ − , and the quadrilateral with vertices (1 , , (1 , − , ( − . , (cid:15) ) , (0 . , (cid:15) ) . The top graph shows a solution with (cid:15) = 0 . and the bottom showsa solution with (cid:15) = 10 − . In both cases, the error is on the order of − . Fig. 3 . Using the Schur complement method, quadrilateral domains can be stitched together tosolve differential equations on meshes. Our method is numerically stable even on a mesh containingboth skinny and fat triangles.
Chebyshev polynomials are an efficient and numerically stable basis for poly-nomial approximation [18]. We define the n th degree Chebyshev polynomial as T n ( x ) = cos( n arccos( x )), for x ∈ [ − ,
1] [8, (18.3.1)].
Spectral methods are especiallyuseful for solving linear differential equations. The goal of using spectral methods s to approximate an infinitely-dimensional linear differential operator with a finite-dimensional sparse matrix operator. We can represent a function as an infinite vectorof Chebyshev coefficients, and we can truncate that vector at an appropriate size toachieve a sufficiently accurate approximation of a smooth function. If we attemptedto create a differential operator D λ that mapped a function in the Chebyshev ba-sis to its λ th derivative in the Chebyshev basis, then that operator would be dense,as most of its matrix elements would be nonzero. Current spectral methods mapChebyshev polynomials to Chebyshev polynomials, giving dense, ill-conditioned dif-ferentiation matrices. In contrast, our method constructs D λ to map a function inChebyshev coefficients to its λ th derivative in the ultraspherical basis of parameter λ .For every positive real λ , ultraspherical polynomials of parameter λ are orthogonalon [ − ,
1] with respect to the weight function (1 − x ) λ − [8, Eq. 18.3.1]. Here, wewill ultraspherical polynomials with positive integer parameter.Using ultraspherical coefficients of parameter λ , the matrix λ th differentiationmatrix D λ is(3.1) D λ = 2 λ − ( λ − λ times (cid:122) (cid:125)(cid:124) (cid:123) · · · λ λ + 1 λ + 2 . . . . These operators are diagonal, but combining different-order derivative coefficientsleads to a problem: different-order D λ transform functions in the Chebyshev basis toultraspherical bases of different parameters. The solution is to create basis conver-sion operators, denoted as S λ , which convert between ultraspherical bases of differentorders. Specifically, S converts functions from the Chebyshev basis to the ultras-pherical basis of parameter 1, and S λ converts functions from the ultraspherical basisof parameter λ to the ultraspherical basis of parameter λ + 1. For example, we canrepresent the differential operator ∂u∂x + 2 x as the matrix operator ( D + 2 S ) u . Theultraspherical conversion operators take the form(3.2) S = 12 −
11 0 −
11 0 . . . . . .. . . , S λ = − λλ +2 λλ +1 − λλ +3 λλ +2 . . . λλ +3 . . .. . . , λ > . Since D λ and S λ are sparse for all valid λ , the final differential operator will be sparseand banded. Let L be some n -order differential operator such that L = c + c D + c D + . . . + c n D n where c , . . . , c n are constants and D n is equivalent to d n dx n . If we wishto solve the differential equation Lu = f , we create the matrix operator L . We alsocreate vectors (cid:126)u and (cid:126)f as vectors of Chebyshev coefficients approximating u ( x ) and f ( x ). Now, we can write L = c S n − · · · S + c S n − · · · S D + · · · + c n − D n − + c n D n ,and we can write the matrix equation L (cid:126)u = S n − · · · S (cid:126)f .This method is not limited to equations with constant coefficients. It is possibleto define multiplication matrices of the form M λ [ f ], where M λ [ f ] will multiply avector in the ultraspherical basis of parameter λ by function f [16]. The bandwidth f multiplication matrices increases with the degree of f , so high speed computationrequires f with low polynomial degree.In order to find a particular solution to the differential equation Lu = f , weneed boundary conditions. An n th order differential equation requires n distinctboundary conditions. Since the highest ordered coefficients of f are often very closeto zero, we can simply replace the last n rows of L and the last n elements of (cid:126)f withboundary condition rows. To set u ( x ) = a , we create the boundary condition rowof [ T ( x ) , T ( x ) , T ( x ) , . . . , T n − ( x )], and the corresponding value of (cid:126)f is set to a . Ifwe use the boundary condition u (cid:48) ( x ) = a , we create the boundary condition row of[ T (cid:48) ( x ) , T (cid:48) ( x ) , T (cid:48) ( x ) , . . . , T (cid:48) n − ( x )], and the corresponding value of (cid:126)f is set to a . Wenow have a sparse, well-conditioned differential operator, assuming that the problemis well-posed [17]. In two dimensions,we can write any polynomial p as its bivariate Chebyshev expansion p ( x, y ) = n (cid:88) i,j =0 a ij T i ( y ) T j ( x ) , x, y ∈ [ − , . To form a matrix operator, we can stack the columns of the coefficient matrix to formthe vector(3.3) u = [ a , a , a , . . . , a n − , , a , a , . . . , a n − , , . . . , a n − ,n − ] T . We use Kronecker products (denoted as ‘ ⊗ ’) in order to construct differential operatorsin two dimensions. For example ∂u∂x is represented as D ⊗ I and ∂u∂y is represented as I ⊗ D . To put together differential operators, we can use the identity ( A ⊗ B )( C ⊗ D ) = AC ⊗ BD , so to represent an operator like ∂ ∂x∂y , we use D ⊗ D . We can also addKronecker products to construct more complicated operators; for instance,(3.4) ∂ ∂x − ∂ ∂x∂y + 2 ∂∂x ⇐⇒ D ⊗ S S − S D ⊗ S D + 2 S D ⊗ S S . A derivation of this two-dimensional sparse operator construction can be found in [17].Boundary conditions in two dimensions are slightly more complicated than in onedimension. Consider solving partial differential equations on a square with resolution n × n . We now remove rows corresponding to a right-hand side term of degree n − n −
1. These 4 n − n − a, b ), we create the boundaryrow(3.5) [ T ( a ) , T ( a ) , . . . , T n − ( a )] ⊗ [ T ( b ) , T ( b ) , . . . , T n − ( b )] . In order to apply Neumann conditions, S − D u will give the first derivative ofChebyshev vector u in Chebyshev coefficients. To find the value the derivative of u atpoint a , we use [ T ( a ) , T ( a ) , . . . , T n − ( a )] S − D u . Thus, a boundary condition rowrepresenting ∂u∂x at point ( a, b ) would be(3.6) [ T ( a ) , T ( a ) , . . . , T n − ( a )] S − D ⊗ [ T ( b ) , T ( b ) , . . . , T n − ( b )]and a boundary condition row representing ∂u∂y at point ( a, b ) would be(3.7) [ T ( a ) , T ( a ) , . . . , T n − ( a )] ⊗ [ T ( b ) , T ( b ) , . . . , T n − ( b )] S − D . . Spectral methods on convex quadrilaterals. The method of constructingmatrices described above works only on square domains. We now adapt existingspectral methods to solve differential equations on quadrilaterals, and eventually,quadrilateral meshes.We use the following bilinear map from the square [ − , to the quadrilateralwith vertices ( x , y ) . . . ( x , y ) in counterclockwise order.(4.1) x = a + b r + c s + d rs, y = a + b r + c s + d rs, where ( x, y ) refers to coordinates on the quadrilateral and ( r, s ) refers to coordinateson the square.(4.2) a = 14 ( x + x + x + x ) , b = 14 ( x − x − x + x ) ,c = 14 ( x + x − x − x ) , d = 14 ( x − x + x − x ) , and a , b , c , and d are similarly defined with y , y , y , and y . The transformationfrom the quadrilateral to the square is more complicated, but it will not be needed.When the quadrilateral domain is mapped onto the square, the differential equa-tions on the quadrilateral are also distorted. Through the use of the chain rule, weobtain(4.3) u x = u r r x + u s s x ,u y = u r r y + u s s y ,u xx = u rr ( r x ) + 2 u rs ( r x s x ) + u ss ( s x ) + u r r xx + u s s xx ,u xy = u rr ( r x r y ) + u rs ( r x s y + s x r y ) + u ss ( s x s y ) + u r r xy + u s s xy ,u yy = u rr ( r y ) + 2 u rs ( r y s y ) + u ss ( s y ) + u r r yy + u s s yy . We still need to evaluate derivatives r x , r y , s x , s y , r xx , r xy , r yy , s xx , s xy , and s yy .Using the Inverse Function Theorem, we can invert the Jacobian of the transformationto obtain(4.4) (cid:18) r x r y s x s y (cid:19) = (cid:18) x r x s y r y s (cid:19) − = 1det( x, y ) (cid:18) y s − x s − y r x r (cid:19) , det( x, y ) = x r y s − x s y r . Using these equations, we can also find second derivatives r xx , r xy , r yy , s xx , s xy , and s yy . These derivatives can be represented as a bivariate cubic polynomial divided bydet( r, s ) .In order to keep our differential operators sparse, we simply multiply the equa-tions and right-hand side by det( r, s ) . That way, we only multiply the operatorby multiplication matrices with polynomial degree <
3, rather than by a dense ap-proximation of a rational function, and our differential operators remain sparse andbanded.
5. Interface conditions.
In order to actually do something useful with quadri-lateral domains, it is necessary to somehow stitch them together into a larger mesh.Here, we will first consider the case of stitching two quadrilaterals together along anedge. We found that the most efficient way to stitch together multiple domains was touse the Schur Complement Method [6]. Essentially, we first solve for the values alongthe interface between the two domains, and then we use that interface to find the ig. 4 . Let the left quadrilateral be u , the right quadrilateral be u , and the middle linebe u Γ . Only n − points on the rightmost edge of u are connected to u Γ , since the top rightcorner is excluded. Similarly, the bottom left corner is excluded from u While this arrangementmay seem unnecessarily complicated on the surface, it facilitates stitching together large numbers ofquadrilaterals. bulk solution inside of each domain. The advantage to using the Schur ComplementMethod is that once we have solved for the interface solution, we can compute thesolution on each element in parallel.In order to use the Schur Complement method on two quadrilaterals, we musthave a linear system of the form(5.1) A A A A A Γ1 A Γ2 A ΓΓ u u u Γ = f f f Γ . Here, u is the solution on the first quadrilateral, u is the solution on the secondquadrilateral, and u Γ is the solution on the interface between u and u . We willrepresent u Γ as a vector of values at Chebyshev points.We also have A and A for solving partial differential equations on quadrilateralelements. Finally, A and A constrain u Γ to u and u , and A Γ1 and A Γ2 ensurethat the solution across the interface is differentiable. A ΓΓ is zero, as it is not necessaryhere.First, we need to set A and A . As shown in Figure 4, if an edge of u isconnected to u Γ , then we connect all of the nodes on that edge of u to u Γ except forthe node on the counterclockwise corner of that edge. This arrangement ensures thatall nodes on the exterior boundary of the domain are assigned the proper boundaryconditions. It also avoids redundant boundary conditions when multiple quadrilateralsmeet at a point on the interior of the mesh.Now, we need to ensure that the solution is differentiable across the interface with A Γ1 and A Γ2 . In order to do so, we need to find the derivative perpendicular to theboundary. Let us say that the boundary has endpoints ( x , y ) and ( x , y ). Let usdefine(5.2) α = ∆ x (cid:112) ∆ x + ∆ y , β = ∆ y (cid:112) ∆ x + ∆ y . et z = βx − αy , so the derivative normal to the boundary is u z = βu x − αu y .We want this derivative in terms of r and s . Using (4.3), we have(5.3) u x = u r r x + u s s x , u y = u r r y + u s s y . From (4.4), we have(5.4) u x = 1det( r, s ) ( u r y s − u s y r ) , u y = 1det( r, s ) ( − u r x s + u s x r ) . By differentiating (4.1), we have(5.5) u x = 1det( r, s ) ( u r ( c + d r ) − u s ( b + d s )) u y = 1det( r, s ) ( − u r ( c + d r ) + u s ( b + d s )) . Finally, we can generate boundary rows for u r and u s with (3.6) and (3.7). Eachrow of A Γ1 represents the derivative of u perpendicular to the interface at one ofthe n Chebyshev points on the interface, and each row of A Γ2 represents the negativederivative of u at each Chebyshev point on the boundary.Now, we want to ensure continuity across the boundary at the two endpoints. Wecan simply replace the differentiability conditions corresponding to the endpoints ofthe interface with continuity conditions from (3.5).
6. Constructing differential operators on meshes.
The key to designing amesh solver is strict bookkeeping. Every quadrilateral in the mesh is given a uniquenumber, as is every vertex and every edge. We can define our mesh with a list ofcoordinates of vertices, and a list of the vertices contained in each triangle. It isimportant that all of the vertices are numbered in counterclockwise order. Next, eachedge must be given a unique global number. Edges also have local definitions—anedge can be expressed as two vertices, or as a quadrilateral number and a local edgenumber on the quadrilateral. We generate arrays to convert between these global andlocal definitions. Finally, we must construct an array that keeps track of whethervertices are on the inside or boundary of the domain. Each vertex is also assigneda single edge adjacent to that vertex. This “vertex list” is essential for constructingnon-singular boundary conditions.Next, we define u Γ such that elements n ( k −
1) + 1 : nk of u Γ correspond to the k th interior edge. Thus, columns n ( k −
1) + 1 : nk of A i Γ and rows n ( k −
1) + 1 : nk of A Γ i correspond to interface conditions across the k th interior edge.We now generate a matrix A i for each quadrilateral. If quadrilateral i and j share an edge with interior edge number k , then columns n ( k −
1) + 1 : nk of A i Γ and columns n ( k −
1) + 1 : nk of A j Γ are set like in the previous section. Next,rows n ( k −
1) + 1 : nk of A Γ i and rows n ( k −
1) + 1 : nk of A Γ j are set to ensuredifferentiability, with the endpoints set to continuity conditions. Finally, we check thetwo endpoints of the edge on the vertex list. If an endpoint is listed as an interiorvertex and the vertex list marks edge k , then the continuity conditions are removedfor that endpoint and replaced with differentiability conditions across the edge at thatendpoint. That way, we make sure that the solution is fully continuous and almostfully differentiable without singularities in the linear system. . Optimization. When solving systems of equations, banded matrices are typi-cally efficient to solve. Our matrices for finding solutions on quadrilaterals are “almostbanded,” meaning that most of the elements lie close to the main diagonal, but a fewrows extend the whole length of the matrix. These dense rows correspond to boundaryconditions, as defined in Section 3. Therefore, sparse LU decomposition and Gaussianelimination have large backfill—they are forced to set many zero matrix elements tononzero values, resulting in a much more computationally intensive solve. We foundthat we can use the Woodbury matrix identity to efficiently replace the rows lyingoutside of the bandwidth, creating a banded matrix that can easily be stored as anLU decomposition. The Woodbury matrix identity can compute the inverse a matrixgiven the inverse of another matrix and a rank- k correction. It satisfies(7.1) ( A + U CV ) − = A − − A − U ( C − + V A − U ) − V A − , where A is n -by- n , U is n -by- k , V is k -by- n , and C is k -by- k [5].In our case, we regard our matrix as A + U CV , where A is a matrix that isfast to solve and U CV is a low–rank correction. We generate V by stacking everyrow of A + U CV that corresponds to a boundary condition, since these rows are theonly dense rows. We define U as a sparse binary matrix with exactly one ‘1’ in eachcolumn and no more than one ‘1’ in each row, and C is simply the identity matrix.The replacement boundary condition rows are 1’s in columns corresponding to thelowest-order coefficients of the solution. Our matrix A is now sparse and banded,so fast sparse LU decomposition is possible. Therefore, once the LU decompositionis completed, back-substitution can be evaluated for each individual right-hand sideextremely rapidly, allowing a system of variables with over 60,000 degrees of freedomto be solved in under a second.Another place for optimization is the matrixΣ = A ΓΓ − (cid:88) j A Γ j A − jj A j Γ . The matrix Σ is typically quite large, so it is important to be efficiently solvable.Consider A Γ j A − jj A j Γ . This matrix has nonzero columns only in columns of A j Γ thatcontain nonzero values, and it has nonzero rows only in rows of A Γ j that containnonzero values. Fortunately, most rows of A Γ j and columns of A j Γ are zero. In orderto optimize the bandwidth of Σ, we must optimize the placement of the nonzeroelements of A Γ j and A j Γ .In general, if there are m different boundaries between quadrilaterals, and n is thediscretization size, then u Γ is composed of m blocks of n elements each, with each blockcorresponding to a boundary between two quadrilaterals. Each boundary is assigneda number, such that elements n ( a −
1) + 1 : na of u Γ correspond to boundary a . Ifboundary a is contained by quadrilateral j , then A Γ j has nonzero rows n ( a − na and A j Γ has nonzero columns n ( a −
1) + 1 : na . Thus, if two quadrilaterals containboundaries a and b , then Σ has a bandwidth of at least ( | a − b | + 1) n . Therefore, inorder to minimize the bandwidth of Σ, we minimize the maximum | a − b | for any twoboundaries contained by the same quadrilateral. The theoretical minimum is O ( √ N ),where N is the number of elements, and standard techniques can efficiently get prettyclose to this minimum. Therefore, computing the LU decomposition of Σ has a timecomplexity of O ( N n ) and a space requirement of O ( N . n ). The time complexityof back-substitution is O ( N . n ).Using the Woodbury identity for all of the matrices A jj and computing A − jj A j Γ has a time complexity of O ( N n ). Back substitution for all A jj has a time complexity rame 500 ( t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Fig. 5 . A vortex street develops around a tilted elliptical object in fluid flow. Red regions havepositive vorticity (counterclockwise) and blue regions have negative vorticity. The simulation isequivalent to air flow through a 1 cm × . m/s. Each time step is . × − s. of O ( N n ), and the storage space for the matrices is also O ( N n ). In order to solvepartial differential equations on a mesh, there is an O ( N n + N n ) precomputationtime, and for each right hand side, the solve time is O ( N . n + N n ). It may bepossible to use an iterative method to multiply a vector by Σ − , reducing the factorsof O ( N ) and O ( N . ) to O ( N ).
8. The Navier–Stokes equations.
We decided to demonstrate our methodwith the direct numerical simulation of the Navier–Stokes equations in two dimensionsat moderately high Reynolds numbers. For incompressible flows, the Navier–Stokesequation are(8.1) ∂ u ∂t + ( u · ∇ ) u + ∇ p = ∇ u, ∇ · u = 0 , where u is the velocity vector field and p is the internal pressure field. Note that theseequations set density and viscosity to 1, so the Reynolds number can be varied solelyby changing the velocity of the fluid or the dimensions of the domain. Since theseequations are coupled and nonlinear, we use the multistep method described in [7]. Wecan use a first-order projection method to linearize and decouple the Navier–Stokes quations. With no-slip boundary conditions, we have(8.2) ∇ u n +1 / − u n +1 / (cid:52) t = ( u n · ∇ ) u n − u n (cid:52) t , u n +1 / = 0 , on ∂ Ω ∇ p n +1 = u n +1 / (cid:52) t , ∂p n +1 ∂ n = 0 , on ∂ Ω u n +1 = u n +1 / − (cid:52) t ∇ p n +1 , where n is the normal vector to the boundary [7].Navier–Stokes simulations are extremely useful in the form of wind tunnel simu-lations, where a test object is placed in a steady flow and analyzed. For an incom-pressible flow, different boundary conditions are applied at the inlet, outlet, walls, andtest object. Typically, the test object will have no-slip boundary conditions, and thewalls will have free-slip boundary conditions. The inlet has a constant fluid velocityand zero pressure gradient, and the outlet has a constant pressure and zero velocitygradient.Figure 5 shows a wind tunnel simulation of an object in a moving fluid. A tiltedellipse was chosen to generate asymmetries in the flow, leading to vortex sheddingearlier in the simulation. Acknowledgments.
We thank Pavel Etingof, Slava Gerovitch, and Tanya Kho-vanova for their work on the MIT PRIMES program, which allowed us to collaboratetogether. In March 2017, the first author was awarded the second prize at the Re-generon Science Talent search worth $175,000 based on this work that will cover thetuition fees at MIT starting in August 2017. We are sincerely humbled by Regen-eron’s very generous support. We thank Grady Wright for giving us advice on theimplementation of the Navier–Stokes simulation. We also thank Dan Fortunato andHeather Wilber. This work was supported by National Science Foundation grantNo. 1645445.
REFERENCES[1]
T. Apel , Anisotropic finite elements: local estimates and applications , vol. 3, Stuttgart, Teub-ner, 1999.[2]
I. Babuˇska and A. K. Aziz , On the angle condition in the finite element method , SIAM J.Numer. Anal., 13 (1976), pp. 214–226.[3]
L. C. Evans , Partial differential equations , American Mathematical Society, 2010.[4]
D. Gilbarg and N. S. Trudinger , Elliptic partial differential equations of second order ,Springer-Verlag Berlin Heidelberg, reprint of 1998 edition, 2015.[5]
N. Higham , Accuracy and stability of numerical algorithms , Second edition, SIAM, 2002.[6]
T. Mathew , Domain decomposition methods for the numerical solution of partial differentialequations , Springer Science & Business Media, 2008.[7]
J. G. Liu , Projection methods I: Convergence and numerical boundary layers , SIAM J. Numer.Anal., 32 (1995), pp. 1017–1057.[8]
F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark , NIST Digital Library ofMathematical Functions , http://dlmf.nist.gov[9]
S. Olver and A. Townsend , A fast and well-conditioned spectral method , SIAM Review, 55(2013), pp. 462–489.[10]
A. T. Patera , A spectral element method for fluid dynamics: laminar flow in a channel ex-pansion , J. Comput. Phys., 54 (1984), pp. 468–488.[11]
J. Ruppert , A Delaunay refinement algorithm for quality 2-dimensional mesh generation. , J.Algor., 18 (1995), pp. 548–585.[12]
J. Schoen , Robust, guaranteed-quality anisotropic mesh generation , PhD Dissertation, Depart-ment of Electrical Engineering and Computer Sciences, University of California at Berkeley,2008. 1413]
J. R. Shewchuk , What is a good linear finite element? interpolation, conditioning, anisotropy,and quality measures , University of California at Berkeley, 73, 2002.[14]
J. R. Shewchuk , Delaunay refinement algorithms for triangular mesh generation. , Computa-tional Geometry, 22 (2002), pp. 21–74.[15]
P. ˇSol´ın, K. Segeth, and I. Doleˇzel , Higher-order finite element methods , Chapman &Hall/CRC Press, 2003.[16]
A. Townsend , Computing with functions in two dimensions ,PhD Dissertation Department ofMathematics, Oxford University, 2014.[17]
A. Townsend and S. Olver , The automatic solution of partial differential equations using aglobal spectral method , J. Comput. Phys., 299 (2015), pp. 106–123.[18]
L. N. Trefethen , Approximation theory and approximation practice , SIAM, 2013.
Appendix A. Determinants.
Let ( x , y ), ( x , y ), and ( x , y ) be the verticesof an arbitrary triangle in counterclockwise order. By the shoelace formula, the areaof the triangle is(A.1) A = 12 (cid:104)(cid:0) x y + x y + x y (cid:1) − (cid:0) x y + x y + x y (cid:1)(cid:105) . We can construct three quadrilaterals from the triangle by partitioning the trianglealong the line segments connecting the centroid of the triangle to the midpoints of thesides. Without loss of generality, let our quadrilateral contain vertex 1 of the triangle,so it has the following vertices in counterclockwise order:(A.2) (cid:0) x , y (cid:1) , (cid:18) x + x , y + y (cid:19) , (cid:18) x + x + x , y + y + y (cid:19) , (cid:18) x + x , y + y (cid:19) . Using equations (4.2) and (A.2), we have(A.3) x = 124 (cid:104)(cid:0) x +5 x +5 x (cid:1) + (cid:0) x − x + x (cid:1) r + (cid:0) x + x − x (cid:1) s + (cid:0) x − x − x (cid:1) rs (cid:105) , and similarly for y . We can use equations (4.1) and (4.4) to find that(A.4) det( r, s ) = (cid:0) b c − b c (cid:1) + (cid:0) b d − b d (cid:1) r + (cid:0) c d − c d (cid:1) s. Plugging in values for b , c , d , b , c , and d , we havedet( r, s ) = 1576 (cid:104)(cid:0) (4 x − x + x )(4 y + y − y ) − (4 y − y + y )(4 x + x − x ) (cid:1) + (cid:0) (4 x − x + x )(2 y − y − y ) − (4 y − y + y )(2 x − x − x ) (cid:1) r + (cid:0) (4 y + y − y )(2 x − x − x ) − (4 x + x − x )(2 y − y − y ) (cid:1) s (cid:105) . (A.5)We can note that(A.6) (cid:32) n (cid:88) i =1 v i x i (cid:33)(cid:32) n (cid:88) i =1 w i y i (cid:33) − (cid:32) n (cid:88) i =1 w i x i (cid:33)(cid:32) n (cid:88) i =1 v i y i (cid:33) = (cid:88) ≤ i