A stable and accurate hybrid FE-FD method

(1)

IT 19 051

Examensarbete 30 hp June 2019

A stable and accurate hybrid FE-FD method

Tuan Anh Dao

Institutionen för informationsteknologi

(2)

(3)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:

Box 536 751 21 Uppsala Telefon:

018 – 471 30 03 Telefax:

018 – 471 30 00 Hemsida:

http://www.teknat.uu.se/student

Abstract

A stable and accurate hybrid FE-FD method

Tuan Anh Dao

We develop a hybrid method to couple finite difference methods and finite element methods in a nonconforming multiblock fashion. The aim is to optimize

computational efficiency when complex geometries present. The proposed coupling technique requires minimal changes in the existing schemes while maintaining strict stability, accuracy, and conservation. Analysis and computational results are shown for a linear problem (to the advection-diffusion equation) and a nonlinear problem (to the viscous Burger's equation) in two spatial dimensions.

Tryckt av: Reprocentralen ITC IT 19 051

Examinator: Jarmo Rantakokko Ämnesgranskare: Ken Mattsson Handledare: Murtazo Nazarov

(4)

(5)

1 Introduction

Both finite element (FE) methods and finite di↵erence (FD) methods are largely employed to obtain numerical solutions of partial di↵erential equations (PDEs), which serve many applications in engineering and especially physics (e.g., fluid dynamics, solid mechanics). The FD methods are often seen as the most efficient scheme for high-order approximations despite being restricted to structured grids. The more flexible FE methods can handle complex geometries. However, the obtained matrices can sometimes be storage demanding and the method itself is considered to be more computationally expensive. The greatest strengths of both methods can be combined by developing a hybrid scheme that allows di↵erent methods and refinement levels to be involved simultaneously. Thus, computational efficiency can be optimized by using the most suitable treatment for each portion of the domain. The major complication is known to consist in designing a method-to-method interface that results in an accurate and stable approximation. An innovative method of this type should also accommodate the demands for simplicity, freedom of node distribution, and minimal interference in the existing schemes to obtain optimal performance while preserving the native properties of the involving methods.

A well-designed numerical framework for energy-stable approximation of time-dependent problems is the combination of the summation-by-parts (SBP) operators [7] (to approximate the govern- ing equation) and the simultaneous-approximation-term (SAT) technique [3] (to impose boundary conditions). Some examples are [4, 6, 10, 13, 17, 20]. A crucial benefit of numerical schemes in this form is that proper formulation strictly prevents nonphysical energy growth, a property often referred to as “strict stability”. The vast majority of previous papers on the SBP-SAT framework were committed to FD methods due to apparent advantages of higher-order operators. Nontrivial geometries are then usually handled with curvilinear grids [17] and multiblock coupling [4]. To carry out more difficult geometries, successful attempts have been done to couple FD methods with unstructured finite volume (FV) methods (e.g., [18,19]). Stability is achieved by modifying the FV scheme at the interface. It was not until the work of [12] when nonconforming grids using mixed schemes were made possible in an energy-stable manner. Nonmatching interfaces are resolved by the newly introduced SBP-preserving interpolation operators again maintaining fundamental properties of block-to-block coupling: strict stability, accuracy and conservation. This result paves the way for many possibilities to couple di↵erent methods that can be written in SBP form. So far hybrid methods of this type are available for the coupling of FD–FD [12], FD–discontinuous Galerkin (dG) method [6] and FD–FV [10].

Attempts to combine FD methods and the FE methods have rarely been successful despite it being of great interest. Several techniques developed for the seismic wave propagation problems involving surface topography (e.g., [11,5]) use the concept of an overlapping region as the interface. With the same study objective in seismology, [9] utilizes the SBP-SAT framework for an added advantage of energy-conservation; however the coupling technique can only handle matching interface, and quadrilateral elements are required for the discretization of FE domain. Although there exists evidence that the standard Galerkin FE also has SBP properties [20], a big obstacle lies within its nondiagonal mass matrix (can be understood as the norm matrix in SBP manner). The general procedure [10] for coupling arbitrary SBP schemes is thus not directly applicable because diagonal norms are required. Some pivotal restrictions in [11, 5, 9] are also made in return for a diagonal mass matrix.

In this thesis, we propose a hybrid method to couple FD methods and FE methods in a nonconforming multiblock fashion. The approximation operators yielding from the formulation of the two methods are proved to have SBP properties. The interface continuity condition is weakly imposed by the SAT technique using the SBP-preserving interpolation operators [12] for nonconforming node distribution. The SAT treatment on the FD side is the same as for FD–FD coupling [12]

(7)

(hence FD operators of any order can be applied) while the treatment on the FE side is written in a general expression in the sense that it does not depend on the other method, the mesh type or the degree of polynomial of the solution space. We show that in both a linear case and a nonlinear case, nondiagonal norm matrices (i.e., the standard FE mass matrices) are applicable without further justification (e.g., mass lumping [20]) to the classical Galerkin formulation. We show further that strict stability, accuracy, and conservation are preserved.

An important part that could have been included in this thesis is the implementation of residual- based viscosity (RV) method [14,15,16] to handle nonsmooth solutions. RV method is a state-of- the-art stabilization technique greatly suited for the present study since an explicit expression can be derived and directly applied on both FD and FE in SBP form. However, within the limit of time, the author could not complete the correct formulation that allows the nonsmooth solution part to traverse the interface.

This thesis is organized as follows. In Section 2, we formulate a general SBP framework for the linear advection-di↵usion problem. The numerical schemes of both FE and FD are derived and proved to have SBP properties in Section 3. The proposed coupling technique is described in Section 4. We show a numerical verification of the coupling treatment and some computational results in Section 5. Section 6 demonstrates the coupling technique for a nonlinear problem (to the viscous Burger’s equation); some numerical results are also shown. Section7is the conclusion.

2 Summation-by-parts operators

2.1 The continuous problem to the linear advection-di↵usion equation

We employ the following 2D advection-di↵usion problem for the main analysis

ut+ a^Tru = "r^Tru, (x, y)2 ⌦ ⇢ R² (1) 1

2(a^Tn a^Tn )u "n^Tru = g(x, y, t), (x, y)2 @⌦, t > 0 (2) u = u0(x, y), (x, y)2 ⌦, t = 0

where a = (a1, a2)^T is a constant vector, " > 0 is a small real number, @⌦ denotes the domain boundary, and n is the normal vector pointing outward at the boundary. The notationr refers to the gradient operator r = ⇣

@

@x,_@y^@ ⌘T

. For the two-dimensional analysis of finite di↵erence schemes, it is more convenient to use the following alternative form of (1),

ut+ a1ux+ a2uy= "(uxx+ uyy), (x, y)2 ⌦ 2 R² (3) where subscript notations x, y indicate the partial derivatives with respect to the corresponding variables.

Definition 2.1. Given two real-valued vector functions u and v, we define the following continuous inner product and its corresponding norm

(u, v)⌦= Z

⌦

u^Tv dx; (u, v)@⌦= Z

@⌦

u^Tv ds;kuk²⌦= (u, u)⌦= Z

⌦

u^Tu dx.

Where no confusion may arise, subscript domain notations are skipped for simplicity (·, ·) ⌘ (·, ·)^⌦ andk·k ⌘ k·k⌦. Well-posedness can be shown for (1) with the Robin boundary condition (2) using the energy method. For continuous problems that can be written in the form ut= F (u), a basic idea of the energy method is to investigate the energy growth utilizing the equality _dt^d kuk²⌦= (u, ut)⌦+ (ut, u)⌦. A well-posed problem ensures that _dt^d kuk²⌦ is nonpositive given zero boundary data

(8)

(energy damping, if applicable, only through the boundary). Applying the following continuous integration rules

(v, a^Tru)^⌦= (v, a^Tnu)@⌦ (a^Trv, u)^⌦

(v,r^Tru)^⌦= (v, n^Tru)^@⌦ (rv, ru)^⌦, (4) we have

(u, ut)⌦= (u, a^Tnu)@⌦+ (a^Tru, u)^⌦+ "(u, n^Tru)^@⌦ "(ru, ru)^⌦ (ut, u)⌦= (a^Tru, u)^⌦+ "(n^Tru, u)^@⌦ "(ru, ru)^⌦.

Adding up the two equalities above yields d

dt

⇣kuk²⌦

⌘+ 2"kruk²⌦= (u, a^Tnu)@⌦+ 2"(u, n^Tru)^@⌦.

It is now easy to see that the boundary condition (2) imposes well-posedness to (1) since bounded energy is assured alongside. Indeed, inserting (2) with zero boundary data g = 0 gives

d dt

⇣kuk²⌦

⌘+ 2"kruk²⌦= (u, a^Tn u)@⌦. (5)

2.2 Discrete advection-di↵usion operators

In this section, we formulate a general framework that can later unify the two numerical methods.

A set of matrix operators (H, B, L) is normally required for any numerical scheme in SBP form, where

⇧ H is an integration operator over the whole domain;

⇧ B is an integration operator over the boundary;

⇧ L is a boundary selection operator.

For the current study, we restrict (H, B, L) into the assumptions: (i) H, B are symmetric, positive definite;H is of size N ⇥N where N is the number of nodes; and B is of size N^@⌦⇥N^@⌦where N@⌦

is the number of nodes on the boundary; (ii)L comprises a nonsquare zero-one matrix of which the only “one” element in each row one-to-one corresponds to a boundary node. Notice that the restriction (i) is weaker than the corresponding requirement in [10] in which, in addition, H and B need to be diagonal. This allows the use of nondiagonal norm matrices yielding from standard Galerkin FE approximations.

Definition 2.2. Given two discrete solutions u, v 2 R^N, we define the following discrete inner product and its corresponding discrete norm

(u, v)_H = u^THv; (Lu, Lv)B= (Lu)^TB(Lv); kuk²_H = (u, u)_H= u^THu.

Advantages of an SBP scheme come from the fact that its discrete operators, by design, mimic the continuous integration-by-parts. Below we define a discrete advection operator and a discrete di↵usion operator that hold properties analogous to the continuous integration rules (4).

Definition 2.3. A matrix operator H ¹Q approximating a^Tr is said to be an SBP advection operator if

Q + Q^T =L^TB Diag(a^Tn)L. (6)

Definition 2.4. A matrix operatorsH ¹R approximating r^Tr is said to be an SBP di↵usion operator if

R = L^TBS A, (7)

(9)

where A = A^T 0 and S is a consistent approximation of n^Tr at the boundary. Moreover, if there exist consistent approximationsD^1x,D^1yof _@x^@ ,_@y^@ respectively, such thatA = D1^T(I2⌦H)D¹, where I2 is the identity matrix of size 2⇥ 2, D¹ =⇥

D^T1x,D1y^T

⇤^T

, then H ¹R is called a complete SBP di↵usion operator (see [13]).

The energy method can be analogously applied to SBP schemes utilizing the properties (6), (7), and the relation _dt^d kuk²_H= u^T_tHu + u^THu^t. A strictly stable approximation ensures that _dt^d kuk²_H is nonpositive given zero boundary data.

3 The two numerical schemes in summation-by-parts form

3.1 The finite di↵erence (FD) scheme

Consider a one-dimensional domain [xl, xr] discretized by n equidistant grid points{x^l⌘ x¹, x2, . . . , xn ⌘ xr}. We first define some general one-dimensional SBP operators for FD.

Definition 3.1. [13] A finite di↵erence operator D1 = H ¹Q approximating _@x^@ using qth-order interior stencils is called a qth-order accurate first derivative SBP operator if H = H^T > 0 and Q + Q^T = Diag( 1, 0, . . . , 0, 1).

Let e1= (1, 0, . . . , 0)^T, en = (0, . . . , 0, 1)^T be column vectors of length n, and d1, dn be row vectors approximating first derivatives at x1, xn.

Definition 3.2. [13] A finite di↵erence operator D2 = H ¹( M e1d1+ endn) approximating

@²

@x² using qth-order interior stencils is called a qth-order accurate second derivative SBP operator if H = H^T > 0, M = M^T 0.

Definition 3.3. Given matrices A = (aij) of size p⇥ q and B of arbitrary size, the Kronecker product of A and B is defined as

A⌦ B = 2 66 64

a11B a12B . . . a1,qB a21B a22B . . . a2,qB

... ... . .. ... ap,1B ap,2B . . . ap,qB

3 77 75.

The following properties of the Kronecker product are frequently used.

Property 3.4. Assuming that matrices A, B, C, D are of suitable sizes, we have (i) (A⌦ B)(C ⌦ D) = (AC ⌦ BD);

(ii) (A⌦ B)^T = (A^T ⌦ B^T).

Given a column vector ¯a, for any matrices B, D of compatible dimensions, a direct consequence of Property3.4(i) is (¯a⌦ BD) = (¯a ⌦ B)D and (¯a^T ⌦ BD) = B(¯a^T ⌦ D).

FD SBP schemes for higher number of spatial dimensions are naturally extended from the one- dimensional operators. Consider a bounded rectangle domain ⌦ = [xl, xr]⇥[y^l, yr]⇢ R²illustrated in Figure1. The discrete solution is approximated on n⇥ m equidistant grid points {(xⁱ, yj), i = 1, 2, . . . , n, j = 1, 2, . . . , m}, where

xi= xl+ (i 1)hx, hx= xr xl

n 1 , yj= yl+ (j 1)hy, hy= yr yl

m 1.

(10)

East

South

North

West

n n

U U U

N

E

S W

1 2 m

Figure 1: A rectangle domain discretized by an equidistant grid

A solution vector u approximating a continuous function is then presented by u =⇥

u1,1, u1,2, . . . , u1,m, u2,1, u2,2, . . . , u2,m, . . . un,1, un,2, . . . , un,m⇤T

where ui,j approximates (xi, yj), which can be interpreted as a “vector of vectors”.

Some necessary operators for the two-dimensional approximation are listed below. Denote the identity matrix of size k⇥ k by I^k. The following finite di↵erence operators will be frequently used in the sequel,

Dx= D1⌦ I^m, Dy = In⌦ D¹, D2x= D2⌦ I^m, D2y = In⌦ D², Hx= Hx⌦ I^m, Hy = In⌦ H^y,

(8)

where Hx, Hyare the one-dimensional norm matrices by x and y direction. For ease of reading, it is worth noting that the first components in all of the above Kronecker products are of size n⇥ n, and the second components are of size m ⇥ m. The corresponding two-dimensional norm is H = HxHy = Hx⌦ H^y.

A finite di↵erence semi-discretization of (3) using diagonal norm SBP operators is given by ut=h

a1Dx a2Dy+ "D2x+ "D2y

iu + SAT. (9)

The following proposition is obvious from the way the above operators are defined, but useful for the connection with finite element schemes in the next section.

Proposition 3.1. The approximation operators used in (9) satisfy the summation-by-parts rules (6) and (7), where

H = H;

H ¹Q = a¹Dx+ a2Dy; H ¹R = D^2x+ D2y. Since this property is intended by definition. No proof is given.

The normal vectors at the north, east, south, west boundaries are nN = (0, 1)^T, nE= (1, 0)^T, nS = (0, 1)^T and nW = ( 1, 0)^T, respectively. The boundary condition (2) can be written more

(11)

precisely at the boundaries 8>

>>

><

>>

:

North : a2 |a²|

2 u "uy= gN(t), East :a1 |a¹|

2 u "ux= gE(t), South : a2 |a²|

2 u + "uy= gS(t), West : a1 |a¹|

2 u + "ux= gW(t).

We show that the following penalty terms weakly forcing (2) yield stability to (9), SAT = SATN + SATE+ SATS+ SATW,

SATN = ⌧N

✓a2 |a²|

2 uN "(Dyu)N gN

◆

⌦ (Hy¹em), SATE = ⌧E(H_x¹en)⌦

✓a1 |a¹|

2 uE "(Dxu)E gE

◆ , SATS = ⌧S

✓ a2 |a²|

2 uS+ "(Dyu)S gS

◆

⌦ (Hy¹e1), SATW = ⌧W(H_x¹e1)⌦

✓ a1 |a¹|

2 uW + "(Dxu)W gW

◆ ,

(10)

where

uN = (In⌦ e^Tm)u, uE = (e^T_n ⌦ I^m)u, uS= (In⌦ e^T1)u, uW = (e^T₁ ⌦ I^m)u, (Dyu)N = (In⌦ d^m)u, (Dxu)E = (dn⌦ I^m)u, (Dyu)S= (In⌦ d¹)u, (Dxu)W = (d1⌦ I^m)u.

(11) The following theorem completes the formulation of the finite di↵erence approximation to (1), (2).

Theorem 3.5. By choosing ⌧N = ⌧E = ⌧S = ⌧W = 1, the approximation (9) with the outer boundary SAT terms given in (10) is stable.

Proof. See AppendixA.

3.2 The finite element (FE) scheme

Weak formulation

Multiplying both sides of (1) by a test function v 2 V := {v(·, t) : kvk²L²(⌦)+krvk²L²(⌦) <1}

and integrating by parts, a weak formulation can be written as: find u2 V such that

(@tu, v) + (a^Tru, v) = ("r^Tru, v) 8v 2 V. (12) Partition the boundary @⌦ into inflow part and outflow part, @⌦ = in[ ^out, where

in={(x, y) 2 @⌦ | a^Tn < 0} and out ={(x, y) 2 @⌦ | a^Tn 0}.

For example, if a 0, it is easy to see that in = S [ ^W consists of the South and West boundaries and out = N [ ^E consists of the North and East boundaries, denoted respectively.

The boundary condition (2) becomes

(a^Tnu "n^Tru = g on in,

"n^Tru = g on out. (13)

(12)

Applying the Green’s formula on (12), we obtain (@tu, v) + (a^Tru, v) + "(ru, rv)

Z

@⌦

("n^Tru)vds = 0, 8v 2 V.

Inserting (13) gives

(@tu, v) = (a^Tru, v) "(ru, rv) + Z

in

(a^Tnu g)vds Z

out

gvds, 8v 2 V.

Galerkin finite element approximation

A finite element approximation can then be formulated as: find uh 2 V^h := {v(·, t) 2 C⁰(⌦) | v|^K 2 P¹(K),8K 2 T^h} such that

(@tuh, v) = (a^Tru^h, v) "(ru^h,rv) + Z

in

(a^Tnuh g)vds Z

out

gvds, 8v 2 V^h. (14) where Th ={Kⁱ}^Ni=1 ⇡ ⌦, N is the number of nodes the mesh and P¹(Ki) is the space of linear interpolations on Ki. We can write Vh= span{'ⁱ}^Ni=1, where 'i, i = 1, . . . , N are piecewise linear basis functions, 'i = 1 at the i^th vertex, zero elsewhere. Since uh2 V^h, there exist real numbers {⇠ⁱ}^Ni=1 such that uh =

XN j=1

⇠j(t)'j(x, y). Inserting this into (14) and letting v = 'i, i = 1, . . . , N gives

M ˙¯⇠(t) = C ¯⇠(t) "A ¯⇠(t) + Rⁱⁿ_M⇠(t)¯ r, (15) where the over dot denotes the time derivative; M, C, A, Rⁱⁿ_M, r are respectively called the mass matrix, the convection matrix, the sti↵ness matrix, the Robin boundary mass matrix, and the Robin boundary vector; to be more precise, for i, j = 1, . . . , N ,

Mij= ('j, 'i), Cij= (a^Tr'^j, 'i), Aij= (r'^j,r'ⁱ), Rⁱⁿ_{M ij}=

Z

in

a^Tn'j'ids, ri=

Z

@⌦

g'ids.

The assembly of the above operators is described in [8]. Interestingly, if we treat (15) in an SBP manner by considering u(t) ⌘ ¯⇠(t) as an approximation vector of the solution at mesh points and applying the energy method with zero boundary data (i.e., writing it in the form ut= F (u), multiplying both side by u^TM, integrating by parts and adding the transpose u^T_tMv), it deduces

that d

dtkuk²M+ 2"u^TAu = u^T(C + C^T)u + 2u^TRⁱⁿ_Mu.

To achieve the above equality, the symmetry and positive definiteness of M is necessary (M ¹)^T = M ¹. We now apply the continuous integration rule (4) on C + C^T to obtain

d

dtkuk²M+ 2"u^TAu = u^TR^@⌦_M u + 2u^TRⁱⁿ_Mu = u^T|R^M| u, where R^@⌦_{M ij} =

Z

@⌦

a^Tn'j'ids. In the interior, the sum C + C^T vanishes because the local contributions of adjoining triangles give opposite signs of the same value at both ends of their shared edge. Since the sti↵ness matrix A is positive definite, the above equality shows stability in the same way regularly seen in SBP finite di↵erence approximations (see AppendixA). It is also analogous to the continuous estimate (5). This observation reveals some SBP properties of (15).

(13)

Proposition 3.2. The approximation operators used in (15) satisfy the summation-by-parts rules (6) and (7), where

H = M;

H ¹Q = M ¹C;

H ¹R = M ¹A + M ¹R^@⌦_C , where R^@⌦_{C ij}=R

@⌦n^Tr'^j'ids, i, j = 1, 2, . . . , N .

Proof. We prove that the FE advection operator Q = C satisfies (6). From the above analysis, we already have C + C^T = R^@⌦_M. We need to show that R^@⌦_M =L^TB Diag(a^Tn)L with the FE formulation of B and L. For a general boundary selection operator L, each row of L must be a one-to-one mapping from its “one” element with a node on the boundary. The combination of a left multiplierL^T and a right multiplierL thus simply expands the local boundary operators into the full domain size (with zeros on the extended rows and columns). The relation (6) becomes obvious with the local boundary integration operator insertedB^ij =R

@⌦'j'ids, for i, j = 1, 2, . . . , N@⌦, where N@⌦is the number of boundary nodes.

For the di↵usion operatorR = A + R^WC, we need to show that A + R^@⌦_C =L^TBS A with the FE formulation ofB, L, D¹, S, and H. Firstly, we formulate S approximating n^Tru (i.e., given uh 2 V^h, S(uh) = n^Tru^h). We seek this approximation in the space Vh. Since n^Tru does not belong to Vh, the projection pS(uh) of S(uh) onto Vh must be used instead. Namely, we find pS(uh)2 V^h such that

(pS(uh) S(uh), v)@⌦= 0, 8v 2 V^h +

(pS(uh), v)@⌦= (n^Tru^h, v)@⌦,8v 2 V^h.

Since pS(uh)2 V^h, there exists some ¯⌘ = (⌘1, ⌘2, . . . , ⌘N) such that pS(uh) =PN

j=1⌘j'j. Inserting this and uh=PN

j=1⇠j'j, the above equation can be written in matrix form,

B⌘ = R˜ ^@⌦C ⇠,¯ (16)

where ˜B^ij =R

@⌦'j'ids, i, j = 1, 2, . . . , N and R^@⌦_C is defined above. Since obviouslyR

@⌦'j'ids = 0 if either i or j corresponds to a nonboundary node,

L^TBL⌘ = ˜B⌘ (17)

The use ofL is simply because B is defined in the SBP framework of size N^@⌦⇥N^@⌦. By definition, L⌘ = S ¯⇠. (S is in general not square because of (7).) From (16) and (17), we obtain

R^@⌦_C =L^TBS. (18)

What is left is to prove the semi-positive definiteness of A. Starting from a symmetric bilinear form in the weak formulation (u, v) := (ru, rv)^⌦, for u, v 2 V , applying the Poincar´e’s inequality yields

(v, v) =krvk²L²(⌦) ckvk²L²(⌦), 8v 2 V.

Since Vh ⇢ V , for all v^h 2 V^h, the property still holds (vh, vh) ckv^hk²L²(⌦). Given ¯⇠ = (⇠1, ⇠2, . . . , ⇠N)^T and vh=PN

j=1⇠j'j, we have

⇠¯^TA ¯⇠ = XN i=1

XN j=1

⇠iAij⇠j= XN i=1

XN j=1

(⇠i'i, ⇠j'j) = 0

@ XN i=1

⇠i'i, XN j=1

⇠j'j

1

A = (vh, vh) 0. (19)

From (18) and (19), the desired property of R is achieved. Note also that at some boundary points where n is not well-defined (e.g., corners of a rectangle domain), a boundary splitting [10]

is necessary for the above proof to be valid. However, this proof routine does not change.

(14)

4 Simultaneous-approximation-term technique for interface coupling

4.1 Continuous analysis

Consider the coupling of two domains ⌦L, ⌦R illustrated in Figure2, where u and v present the parts of solution on the left and right domains, respectively,

ut+ a^Tru = "r^Tru, (x, y)2 ⌦^L

vt+ a^Trv = "r^Trv, (x, y) 2 ⌦^R. (20)

Ω Ω

L R

Ω_I

Figure 2: Coupling of two adjoining blocks

The North, South, West boundaries of ⌦L and the North, East, South boundaries of ⌦R will be referred to as the “outer boundary” and the shared edge of the two adjoining blocks is called the

“interface”, denoted ⌦I. Applying the energy method on (20), a total energy estimate of the whole solution on ⌦L[ ⌦^R reads

d dt

⇣kuk²⌦L+kvk²⌦R

⌘+ 2"kruk²⌦L+ 2"krvk²⌦R = BT + IT,

where BT contains all the outer boundary terms (being bounded with the boundary condition (2)) and IT contains the interface terms given by

IT = (u, a^TnEu)⌦I+ 2"(u, n^T_Eru)^⌦^I (v, a^TnWv)⌦I + 2"(v, n^T_Wrv)^⌦^I

= a1kuk²⌦I + a1kvk²⌦I+ 2"(u, ux)⌦I 2"(v, vx)⌦I.

Therefore, the necessary continuity condition for an accurate coupling (also for conservation IT ⌘

0) is (

u = v ux= vx

at x2 ⌦^I. (21)

4.2 Necessary interface conditions for stability and conservation

An SBP-SAT semi-discretization of (20) is given by

ut+HL¹Q^Lu = "HL¹R^Lu + SATL+ SATIL

vt+HR¹Q^Rv = "HR¹R^Rv + SATR+ SATIR, (22) where SATR,L present the weak boundary treatment, and SATIR,IL present the weak interface treatment. Another advantage of this weakly forcing technique is that SATR,L,IR,IL each can be

(15)

treated separately by the energy method. Therefore, in the following subsections, we assume the stable outer boundary treatment as described in Section 3 and only analyze the estimate terms regarding the additional involvement of SATIR,IL.

I_R2L IL2R

Figure 3: Example of translating nodal values on a nonmatching interface using interpolation operators

The nonmatching nodes on ⌦I are resolved by the interpolation operators IL2R and IR2L, where IL2Rmaps the interface values from the left domain to the right domain, and IR2Ldoes the opposite direction (illustrated in Figure 3). For proven stability and energy conservation, the use of the following class of interpolation operators is necessary (first introduced in [12]).

Definition 4.1. LetH^Ly,Hy^R be some integration operators over the interface of the left and the right domains, respectively. {I^L2R, IR2L} is called a pair of SBP-preserving interpolation operators if

Hy^RIL2R= I_R2L^T H^Ly. (23)

For example, for FD,H^yis the Hy-norm that we defined ealier. For FV, this operator is a diagonal matrix containing the Euclid distances between two neighboring interface nodes [10]. For dG, this operator is a diagonal matrix presenting integration weights for interface edges [6]. However, it is quite intuitive that this matrix operator is not diagonal for FE. We will investigate this in more detail in Section4.4.

4.3 FD–FD coupling

We consider the case where FD schemes are used on both ⌦Land ⌦R. The coupling SAT treatment is given by

SATIL= ↵L(H_x^L ¹en)⌦ (u^I IR2LvI)

+ L(H_x^L ¹en)⌦ ((D^xu)I IR2L(Dxv)I) + L(H_x^L ¹d^T_n)⌦ (u^I IR2LvI), SATIR= ↵R(H_x^R ¹e1)⌦ (v^I IL2RuI)

+ R(H_x^R ¹e1)⌦ ((D^xv)I IL2R(Dxu)I) + R(H_x^R ¹d^T₁)⌦ (v^I IL2RuI),

(24)

where uI ⌘ u^E, vI ⌘ v^W, (Dxu)I ⌘ (D^xu)E, (Dxv)I ⌘ (D^xv)W are defined in (11).

A stable and accurate hybrid FE-FD method

Examensarbete 30 hp June 2019

A stable and accurate hybrid FE-FD method

Tuan Anh Dao

Institutionen för informationsteknologi

Abstract

A stable and accurate hybrid FE-FD method

Contents

1 Introduction

2 Summation-by-parts operators

3 The two numerical schemes in summation-by-parts form

4 Simultaneous-approximation-term technique for interface coupling