Weinan E. Principles of Multiscale Modeling

Подождите немного. Документ загружается.

3.3. ADAPTIVE MESH REFINEMENT 117

level is 9(N/s)s

= 9Ns since s

operations are required to calculate the interaction

between one box and its neighbors.

It should be noted that while the multi-pole expansion is oriented towards the sources,

the local expansion is oriented toward the targets. In addition, when performing the ﬁnal

summation, the focus of the O(N log N) algorithm is on the targets: For a given target,

the contributions from the sources in its interaction list are summed up from scale to

scale. In FMM, however, the focus is on the sources. At each scale, for each box, one

asks what the contributions are for the sources in that box to the boxes in its interaction

list. Previously we were examining the interaction list of the targets. Now we examine

the interaction list of the sources. This dual viewpoint is very important for generalizing

analysis-based fast summation algorithms to other settings.

3.3 Adaptive mesh reﬁnement

Adaptivity is a central issue in modeling and computation. There are diﬀerent kinds

of adaptivity: Adaptively choosing the models, discretization methods, stopping criteria,

time step size, and the mesh. Well-known examples of adaptive numerical techniques

include:

1. Adaptive numerical quadrature in which new quadrature points are added based

on some local error indicators [15].

2. ODE solvers with adaptive time step size selection [15].

3. Adaptive mesh reﬁnement for solving PDEs [1].

Here we will focus on adaptive mesh reﬁnement for solving PDEs. This has been devel-

oped along two diﬀerent lines:

1. Adaptive ﬁnite element methods pioneered by Babuska and Rheinboldt [3, 4], which

are based on rigorous a posteriori error estimates.

2. Adaptive ﬁnite diﬀerence or ﬁnite volume methods for nonlinear time-dependent

PDEs such as the ones that arise in gas dynamics, pioneered by Berger, Collela,

Oliger, et al. [8]. In this case, it is much harder to develop rigorous a posteriori

error estimates. Thus, quantities such as the gradients of the numerical solutions

are often used as the error indicators.

118 CHAPTER 3. CLASSICAL MULTISCALE ALGORITHMS

Recently, both strategies have been extended to include the possibility of adaptively

selecting the physical models [1, 23], a topic that we will return to in Chapter 7.

3.3.1 A posteriori error estimates and local error indicators

The ﬁrst thing we need for an adaptive mesh reﬁnement scheme is some local error

indicators that are easily computable. This is the purpose of developing a posteriori

error estimates. In numerical analysis, we are usually more familiar with the a priori

error estimates, which are typically in the form:

ku − u

≤ Chkuk

. (3.3.1)

Here kvk

and kvk

are the H

and H

norms of the function v respectively. (3.3.1) tells

us that the convergence rate is ﬁrst order. But as a quantitative measure of the error, it

has limited value since we do not know the norm of u th at appears at the right hand side

of (3.3.1). In addition, (3.3.1) gives an estimate of the global error, not the local error

around each element.

The purpose of a posteriori error estimates is to ﬁnd easily computable quantities

)} that gives an estimate of the error near the element K:

ku − u

1,K

∼ e

) (3.3.2)

This kind of work was pioneered by Babuska and Rheinboldt [3, 4]. A rather compre-

hensive treatment of the subject is presented in [1].

To give a typical example of such a posteriori error estimates, let us consider the

second order elliptic equation in two dimension:

Lu(x) = −∇ · (a(x)∇u(x)) = f(x), x ∈ Ω

with Dirichlet boundary condition u(x) = 0, x ∈ ∂Ω. Consider a piecewise linear ﬁnite

element method on a regular mesh and denote by u

the ﬁnite element solution. For

each element K, deﬁne the residual: r

= Lu

− f. For each edge Γ that is shared by a

pair of elements K and K

′

, deﬁne J

) to be the jump of n ·(a∇u

) across that edge,

where n is the normal to the edge (the sign of n does not matter).

Deﬁne

ε(u

) =

dx +

1/2

(3.3.3)

3.3. ADAPTIVE MESH REFINEMENT 119

Figure 3.5: Reﬁning the elements

Here the ﬁrst sum is done over all elements, and h

is the diameter of the element K.

The second sum is over all interior edges an d h

is the length of the edge Γ. It can be

shown [1, 21, 38] under various conditions that

ε(u

) ≤ ku − u

≤ C

ε(u

) (3.3.4)

where C

and C

are constants. Moreover, each term inside the expression of ε(u

) is a

measure of the local error near that element. Therefore, these quantities can be used as

a criterion for mesh reﬁnement.

Other types of a posteriori error estimators are discussed in [1] (see also [5, 43]).

In some cases, our interest is not on obtaining an accurate approximation of the

whole solution, but rather some quantities associated with the solution. For example,

when modeling ﬂow past obstacles, we might be most interested in the total drag or lift

on the body. In this case, we might want to reﬁne more in places, say near the body,

120 CHAPTER 3. CLASSICAL MULTISCALE ALGORITHMS

where the accuracy of the approximate solution is more imp ortant for the quantities we

are interested in. The procedure outlined above has been extended to this case [1]. In

such a “goal-oriented adaptive mesh reﬁnement scheme”, the error indicators are usually

obtained from solving an adjoint problem.

3.3.2 The moving mesh method

Next we discuss brieﬂy an alternative procedure for adapting the mesh, the moving

mesh method. This is a way of redistributing the grid points, so that more computational

eﬀort is directed at places where it is needed, e.g. places where the solution has large

gradients. The moving mesh method has its origin in the moving ﬁnite element method

[32, 17]. Current versions of the moving mesh method has largely followed the framework

proposed by Winslow, in which the mesh is regulated using a monitor function [39], which

is a carefully designed function of the numerical solution. The quality of the numerical

mesh is improved by minimizing the monitor function.

There are two slightly diﬀerent ways of implementing the moving mesh method. One

is to introduce an auxiliary computational domain on which discretization is carried out,

usually on a regular grid. The computational domain and the physical domain are related

by a mapping, which is obtained adaptively using the monitor function. This approach

can be regarded as an adaptive (or dynamic) change of variable. One very elegant way

of implementing this is to formulate the original PDE in the computational domain, and

supplement the resulting PDE with an additional equation for the mesh [29, 33, 34].

Denote by α the variable in the computational domain and x the variable in the

physical domain. They are related by the function α = α(x). We look for the function

α such that a regular grid in the α domain is mapped to the desired grid in the physical

domain, with grid points clustered in the areas of interest. This can be done using a

variational formulation with the functional

E(α) =

w(x)

|∇

α|

dx, (3.3.5)

Here w = w(x) is the monitor function. The Euler- Lagrange equation for the minimiza-

tion of this functional is:

∇ ·



w(x)

∇



= 0 (3.3.6)

3.3. ADAPTIVE MESH REFINEMENT 121

In one dimension, this reduces to

w(x)

dα

= C. (3.3.7)

Assume that we have a uniform grid in the α domain, with α

= n∆α, n = 1, 2, ··· , N.

Let α

= α(x

). Integrating over [α

n−1

, α

] gives

n−1

w(x)dx = C∆α (3.3.8)

i.e., ∆x · w(x) ≈

n−1

w(x)dx is approximately a constant. The grid points are concen-

trated in the regions where the monitor function w is large.

The other approach is to work directly with the original physical domain, and changes

the mesh from time to time using the monitor function. At each step, one ﬁrst minimizes

the functional (3.3.5) to ﬁnd the new mesh. This is followed by an interpolation step to

provide the values of the numerical solution on the new mesh. This approach has been

explored extensively in [29] and subsequent works.

In cases when a truly high quality mesh is required, as is the case when the solution

is singular, an iterative procedure is used to improve the quality of the mesh [34].

The advantage of the ﬁrst approach is that one works with a standard regular grid

on the computational domain. There are no complications coming from, for example,

numerical interp olation. The disadvantage is that by moving away from the physical

domain, the structure of the original PDE, such as conservation form, might be lost. The

second approach, on the other hand, is exactly complementary to the ﬁrst.

Finding the right monitor function is a crucial component of the moving mesh m ethod.

The commonly used monitor functions are in the form, for example,

w(x) =

1 + u(x)

+ |∇u(x)|

(3.3.9)

where u is the numerical solution.

Compared with the adaptive mesh reﬁnement strategies discussed earlier, the moving

mesh method has the advantage that it essentially works with a regular and smooth

mesh. Therefore it has the nice properties of a regular, smooth mesh, such as the super-

convergence properties. On th e other hand, the moving mesh method is in general less

ﬂexible, since its control of the mesh is rather indirect, through the monitor function.

122 CHAPTER 3. CLASSICAL MULTISCALE ALGORITHMS

3.4 Domain decomposition methods

The intuitive ideas behind the domain decomposition method (DDM) is quite simp le:

The computational domain is divided into smaller subdomains and a computational

strategy is established based on solving the given problem on each subdomain and making

sure that the solutions on diﬀerent subdomains match. There are several advantages for

such a strategy:

1. It is very well suited for parallel computation.

2. DDM can be used as a localization strategy, as in electronic structure analysis [40].

3. According to the nature of the solution on the diﬀerent subdomains, one may

use diﬀerent models or diﬀerent algorithms. For this reason, DDM has been used

extensively as a platform for developing algorithms with multi-physics models.

The main challenge in DDM is to make sure that the solutions on diﬀerent subdomains

match with each other. To discuss how this is done, we will distinguish two diﬀerent cases,

when the subdomains do or do not overlap. See Figure 3.6.

3.4.1 Non-overlapping domain decomposition methods

We ﬁrst discuss the case when the subdomains do not overlap. We will use the

example:







−△u(x) = f (x), x ∈ Ω

u(x) = 0, x ∈ ∂Ω

(3.4.1)

to illustrate the main issues.

After discretization on Ω, say using a ﬁnite diﬀerence method, we obtain a linear

system:

Au = b. (3.4.2)

To solve this linear problem, we divide th e domain Ω into a union of two non-overlapping

subdomains Ω

and Ω

, see Figure 3.6. We denote by u

and u

the solution in Ω

and

Ω

respectively, and by w the solution on the interface between Ω

and Ω

. The width of

3.4. DOMAIN DECOMPOSITION METHODS 123

Ω

Overlapping domain decomposition

Nonoverlapping domain decomposition

Figure 3.6: Illustration of overlapping domain decomposition (courtesy of Xiang Zhou).

the interfacial region is chosen so that the components u

and u

do not interact directly,

i.e. if we let u = (u

, u

, w)

, then (3.4.2) can be written as







0 A































. (3.4.3)

Eliminating u

and u

, we obtain

Cw = d (3.4.4)

where d = b

− A

−1

− A

−1

, C is the Schur complement of A

C = A

− A

−1

− A

−1

C is usually a full matrix and computing C explicitly is costly. Fortunately there are

many numerical algorithms for solving (3.4.4) which only require evaluating the action

of C on a given vector, not C itself. One example of such algorithms is the the conjugate

gradient method. From (3.4.1), we see that evaluating Cv for a given vector v amounts

to solving the linear systems of equations on the subdomains Ω

and Ω

Cv = A

v − A

− A

(3.4.5)

where v

and v

are solutions of

= A

v, A

= A

v (3.4.6)

124 CHAPTER 3. CLASSICAL MULTISCALE ALGORITHMS

These two equations can be solved independently, thereby decoupling the the solutions

on the two domains Ω

and Ω

. The same can be said about evaluating d.

To understand the nature of C, it helps to introduce the Dirichlet to Neumann oper-

ator (DtN) for a domain Ω and a subset of the boundary ∂Ω, denotd by Γ. Consider the

problem

− △u(x) = f(x), x ∈ Ω (3.4.7)

u(x) = 0, x ∈ ∂Ω\Γ (3.4.8)

This is the same as (3.4.1), except that the boundary condition is speciﬁed only on ∂Ω\Γ,

therefore the problem is not fully determined. Now given a function w deﬁned on Γ. We

can obtain a unique solution of (3.4.7) with the additional condition that u = w on Γ.

Denote this solution by U

. The Dirichlet to Neumann operator acting on w is deﬁned

as the normal derivative of U

on Γ. Hence it is a map between functions deﬁned on Γ.

Consider the special case when Ω = R

= {(x

, x

), x

> 0} and Γ = ∂Ω = R

. Given

w on ∂Ω = R

, we can solve (3.4.7) with f = 0 and w as the boundary condition u = w

on ∂Ω using Fourier transform in the x

direction:

, x

) = e

−|k

ˆw(k

)

where k

is the variable for the Fourier transform in the x

direction,

f denotes the

Fourier transform of f. Hence, on the Fourier side, we have

∂

∂x

, x

= |k

|ˆw(k

)

From this we see t hat the Fourier symbol of the Dirichlet to Neumann operator is |k

This means in particular that it is eﬀectively a ﬁrst order operator.

Now we can understand the nature of C as follows. Given some candidate function

w, we can then solve for the corresponding u

and u

, from the ﬁrst two equations in

(3.4.3):

= b

− A

w (3.4.9)

= b

− A

w (3.4.10)

and u

have the same values at the interface, but their normal d erivatives do not

necessarily match. To obtain a solution of the original problem over the whole domain

3.4. DOMAIN DECOMPOSITION METHODS 125

Ω, we should require that u

and u

have the same normal derivative at the interface.

This is the equation (3.4.4) that w solves. Therefore C can be regarded as a discretization

of the diﬀerence of the Dirichlet to Neumann operator (DtN) on Ω

and Ω

. Since the

DtN operator is a ﬁrst order operator, the condition number of C should scale as O(h

−1

)

where h is the grid spacing. This is an improvement over the condition number of the

original problem (3.4.3) which is usually O(h

−2

3.4.2 Overlapping domain decomposition methods

Another option is to allow the subdomains to overlap, as indicated in Figure 3.6.

One can then set up an iterative procedure along the following lines: Starting with some

(arbitrary) boundary condition on ∂Ω

, one solves for u

(over Ω

). From u

, one can

obtain some boundary condition on ∂Ω

and use it to solve for u

(over Ω

). From u

one can obtain some bou ndary condition on ∂Ω

and use it to solve for u

. From there

one goes on to solve for u

as before. Continuing this procedure, one obt ains u

, u

, ···.

This is the well-known Schwarz iteration.

The main issue is how to get the boundary conditions. The boundary conditions can

either be Dirichlet or Neumann type, or other boundary cond itions. The original Schwarz

iteration uses a simple Dirichlet boundary condition. If the boundary conditions are not

chosen properly, such an iteration may not converge. Th e following is a simple example

that illustrates this point.

Consider the problem −u

′′

= 0 with boundary condition u(0) = 1 and u(1) = 0. The

exact solution is u(x) = 1 − x. We deﬁne Ω

= [0, b] and Ω

= [a, 1] with the overlap

region [a, b], where 0 < a < b < 1. It is easy to see that if we use the Dirichlet boundary

condition at a and b, the algorithm converges. The convergence rate depends on the size

of the overlapping region (see Section 7.4)

r =

1 − b

1 − a

But if we use instead the Neumann boundary condition at a and b, it is easy to see that

, u

} do not change and the procedure does not give the right solution. We will return

to this example in Chapter 7.

126 CHAPTER 3. CLASSICAL MULTISCALE ALGORITHMS

Figure 3.7: The overlapping domain decomposition for one dimension problem u

′′

= 0.

The dotted lines are the solutions of u

1,2

with Dirichlet condition for local problems; the

dash (initial ﬂux u

′

is zero) and dot-dash (initial ﬂux u

′

is −0.3) lines are the solutions

of local problems with Neumann condition (courtesy of Xiang Zhou).