Weinan E. Principles of Multiscale Modeling

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

6.1. SOME CLASSICAL EXAMPLES 257

υ υ υ

∆t

∆x

Figure 6.2: Schematics for the derivation of kinetic schemes: A ﬁnite volume method is

imposed in the x−t domain, and the kinetic equation is solved (e.g. analytically) over the

shaded region to give the ﬂuxes needed in the ﬁnite volume method. The v axis indicates

the extra velocity variable in th e kinetic model, which represents the microstructure for

the present problem.

258 CHAPTER 6. CAPTURING THE MACROSCALE BEHAVIOR

6.1.4 Cloud-resolving convection parametrization

Cloud-resolving convection parametrization (CRCP), or super-parametrization, is a

technique introdu ced by Grabowski and Smolarkiewicz [42] in which the large scale trop-

ical dynamics of the atmosphere is captured by embedding cloud scale models at each

grid point of a large scale model. The macroscale model is a large eddy simulation model.

However, in contrast to standard large eddy simulation, the eﬀects of the sub-grid scales

are modeled not by empirical models, but rather through the introduction of explicit

cloud scale models. We will omit the details here and refer the interested reader to the

original article [42]. For more recent development on this subject, see [75].

6.2 Multi-grid and the equation-free approach

Given the success of these diﬀerent multiscale methods, it is natural to ask whether

we can formulate a general framework. In the case of solving diﬀerential equations,

such general strategies have proven to be very useful. A good example is the ﬁnite

element method which has provided not only a general designing principle for numerical

algorithms but also general guidelines for error analysis. In recent years, several attempts

have b een made to construct general strategies for multiscale methods. In the following,

we will discuss three examples of such general strategies, the extended multi-grid method,

the heterogeneous multiscale method, and the equation-free approach.

6.2.1 Extended multi-grid method

One of the ﬁrst attempts to construct a general framework for multiscale modeling

was the extension of the multi-grid met hod by Achi Brandt. In its original form [4, 44],

the multi-grid method is an eﬃcient way of solving the algebraic equations obtained from

discretizing PDEs (see Section 3.1). The objective was to ﬁnd accurate solutions of the

PDE. Since the 1990’s, Brandt and others have discussed the possibility of extending the

traditional multi-grid method in a number of directions [5]:

1. The models used at the diﬀerent scales can be quite general and quite diﬀerent in

nature. For example, one may use Monte Carlo methods or molecular dynamics at

the small scale, and traditional continuum models at the large scale.

6.2. MULTI-GRID AND THE EQUATION-FREE APPROACH 259

2. The traditional multi-grid method is an eﬃcient solver for the ﬁne scale details of

the problem. In the extended version, one might be only interested in capturing

the macroscale behavior using macroscale models. This is particularly useful in

cases when closed-form macroscale models are not available, and it is desirable to

have numerical algorithms that are capable of capturing the macroscale behavior

directly. In this case, one might even be able to reconstruct the eﬀective macroscale

model at the end of the computation. To quote from [5]: “At suﬃcient coarse level,

this entire algorithm eﬀectively p roduces macroscopic ‘equations’ for the simulated

system .... This can yield a macroscopic numerical description for the ﬂuid even for

those cases where the traditional derivation of closed-form diﬀerential equations is

inapplicable.”

3. One may limit the microscopic simulation at the ﬁne levels to small spatial-temporal

domains. As was remarked by Brandt: “few sweeps are enough, due to the fast

CMC (conditional Monte Carlo) equilibration. This fast equilibration also implies

that the interpolation can be done just over a restricted subdomain, serving as

window: In t he window interior ﬁne-level equilibration is reached.”

4. Renormalization group ideas are integrated into the new multi-grid method.

The general procedure is quite similar to standard multi-grid methods. At each cycle,

one repeats the following steps at the diﬀerent levels of grids [5]:

1. Interpolation: The current values of the state variables on the coarse level are used

to initialize or constraint the (micro) model at the ﬁne level.

2. Equilibration: Run the ﬁne-level model (on small windows) for a number of steps.

3. Restriction (projection): Project the variables at the ﬁne level back to the coarse

level.

This is a two-level procedure. Generalization to multiple levels is straightforward. Unlike

in traditional multi-grid methods in which the successive coarse and ﬁne levels are simply

nested ﬁnite diﬀerence or ﬁnite element grids, here the coarse and ﬁne levels can be very

diﬀerent and very far apart.

260 CHAPTER 6. CAPTURING THE MACROSCALE BEHAVIOR

Brandt also described applications to many diﬀerent areas, including electronic struc-

ture analysis, solving integro-diﬀerential equations, modeling high frequency wave prop-

agation, Monte Carlo methods in statistical mechanics, complex ﬂuids and image pro-

cessing [5].

Brandt’s picture is very broad and very appealing. However, at an algorithmic level,

it is quite unclear to the present author how one can carry out such a program in practice,

even though a lot has been written on the subject (see, for example, [5] and the references

therein). We will return to this discussion later in the chapter.

6.2.2 The equation-free approach

The equation-free approach is another framework for addressing the same kind of

questions based on a similar philosophy. It consists of a set of techniques including

coarse bifurcation analysis, projective integrators, the gap-tooth scheme and the patch

dynamics [51]. At an abstract level, equation-free is built on a three-stage “lift-evolve-

restrict” coarse time-stepper procedure [51]. Given a macroscopic state, U(t), at some

time t, the coarse time-stepper consists of the following basic elements [51]:

1. Lift: Transform the macroscale initial data to one or more consistent microscopic

realizations.

2. Evolve: Use the microscopic simulator (the detailed time-stepper) to evolve these

realizations for the desired short macroscopic time τ , generating the value(s) u(τ, t).

3. Restrict: Obtain the restriction of u and deﬁne the coarse time-stepper solution as

U(t) = Mu(τ, t).

Here M is the restriction operator that maps micro-state variables to macro-state vari-

ables.

So far this is very similar to a two-level “interpolate-equilibrate-restrict” multi-grid

procedure. The additional component of the equation-free approach is to exploit scale

separation in the problem using interpolation (in space) and extrapolation (in time)

techniques.

Coarse projective integration

6.2. MULTI-GRID AND THE EQUATION-FREE APPROACH 261

In coarse projective integrators, one generates an ensemble of microscopic solutions

for short times, by running the microscopic solver with an ensemble of initial data that are

consistent with the current macro-state. One evaluates the average values of the coarse

variables over this ensemble. The time derivatives for the coarse variables are computed

using these averaged values and these coarse time derivatives are used to extrapolate the

coarse variables over a much larger time step.

Consider the example







= f(x

, y

, t), x

(0) = x

g(x

, y

, t), y

(0) = y

(6.2.1)

where (x, y) ∈ R

×R

. Here the macroscale variable of interest is x. For this problem,

the projective integrator proceeds in two steps:

1. Relaxation step: In the simplest case of a forward Euler scheme with micro time

step δt, we have







n,m+1

n,m

+ δtf (

n,m

, t

+ mδt)

n,m+1

n,m

δt

n,m

, t

+ mδt)

(6.2.2)

with some initial data

n,0

and

n,0

, m = 0, 1, ··· , M − 1.

2. Extrapolation step: Compute the approximate value of

and then extrapolate:

n,M

− x

n,M−1

δt

(6.2.3)

n+1

= x

+ ∆t

If one wants to construct higher order schemes, one uses higher order extrapolation:

n,M

− 2x

n,M−1

+ x

n,M−2

(δt)

(6.2.4)

n+1

= x

+ ∆t

(∆t)

and so on.

A main advantage of the projective integrators is that they are relatively simple in

the cases when they apply.

262 CHAPTER 6. CAPTURING THE MACROSCALE BEHAVIOR

Dealing with the space variable

Brandt noticed that due to the fast equilibration of the microscopic model, it might

be enough to conduct the microscopic simulations on “small windows”. The “gap-tooth

scheme” builds on a similar intuition.

The basic idea of the gap-to oth scheme is to “use the microscopic rules themselves,

in smaller parts of the domain and, through computational averaging within the sub-

domains, followed by interpolation, to (we) evaluate the coarse ﬁeld U(t, x), the time-

stepper, and the time derivative ﬁeld over the entire domain” [51].

“Given a ﬁnite dimensional representation





of the coarse solution (e.g. nodal

values, cell averages, spectral co eﬃcients, coeﬃcients for ﬁnite elements or empirical basis

functions) the steps of the gap-tooth scheme are the following.

1. Boundary conditions. Construct boundary conditions for each small box b ased on

the coarse representation





2. Lift. Use lifting to map the coarse representation





to initial data for each

small box.

3. Evolve. Solve the detailed equation (the microscale model) for time t ∈ [0, τ] in

each small box y ∈ [0, h] ≡ [x

− h/2, x

+ h/2] with the boundary conditions and

initial data given by steps (1) and (2).

4. Restrict. Deﬁne the representation of th e coarse solution at th e next time level by

restricting the solutions of the detailed equation in the boxes at t = τ” [51].

One interesting application of the gap-tooth scheme is presented in [38] in which

a particle model was used to capture the macroscale dynamics of the viscous Burgers

equation. The microscale mo del is a one-dimensional Brownian dynamics model:

˙x

(t) = u

(t) + ˙w

(t), (6.2.5)

where x

(t) is the position of the j-th particle at time t, the { ˙w

}’s are ind ependent

white n oises, u

(t) is the drift velocity, which is deﬁned to be the density of particles

over an interval of length δ centered at x

(t). To apply the gap-tooth scheme, microscale

simulation boxes of size δ are put uniformly on the real axis. Particle simulations of

(6.2.5) are carried out inside the boxes. At each macroscopic time step n, particles are

initialized according to the known local (macroscopic) density at that time step. The

6.2. MULTI-GRID AND THE EQUATION-FREE APPROACH 263

boundary condition for the particle simulation is imposed as follows. For the k-th box,

there is an outgoing ﬂux (number of particles that exit the box per unit time) and

there is an incoming ﬂux at the left and right side of the box. These are denoted by

ℓ,+

, j

ℓ,−

, j

r,+

, j

r,−

respectively. The outgoing ﬂuxes j

r,+

, and j

ℓ,−

can be measured from

the simulation. The incoming ﬂuxes j

ℓ,+

, j

r,−

will have to be prescribed. One possibility

is to set

ℓ,+

= αj

r,+

k−1

+ (1 − α)j

r,+

(6.2.6)

r,−

= αj

ℓ,−

k+1

+ (1 − α)j

ℓ,−

(6.2.7)

where α = δ/∆x, δ is the size of the box, ∆x is the distance between the centers of

adjacent boxes. To implement this boundary condition, particles that come out of a box

are either “tele-transported” to the next b ox, or put back into the same box from the

other side, according to a probability distribution speciﬁed by the boundary condition

shown above. The above corresponds to a linear interpolation. Quadratic interpolations

are also discussed in [38].

To update the particle density at the next time step t

n+1

= (n + 1)∆t, the particle

densities at that instant of time inside the boxes are computed, and interpolated in space

to give the macroscopic density proﬁle at the new time step.

Patch dynamics

Patch dynamics is a combination of the projective integrators and the gap-tooth

scheme. It consists of the following [51]:

1. Short time steps. Repeat the gap-tooth time stepper a few times, i.e. compute patch

boundary conditions, lift to patches, evolve microscopically and then restrict.

2. Extrapolate. Advance coarse ﬁelds for a long time step into the future through

projective integration. This ﬁrst involves estimation of the time-derivatives for the

coarse ﬁeld variables, using the successively reported coarse ﬁelds followed by a

large projective step.

An example of patch dynamics will be presented in Section 6.7.

264 CHAPTER 6. CAPTURING THE MACROSCALE BEHAVIOR

6.3 The heterogeneous multiscale method

6.3.1 The main components of HMM

We now turn to the framework of the heterogeneous multiscale method (HMM). The

general setup is as follows. At the macroscopic level, we have an incomplete macroscale

model:

∂

U = L(U; D) (6.3.1)

where D denotes the data needed in order for the macroscale model to be complete. For

example, for complex ﬂuids, U might be the macroscopic velocity ﬁeld and D might be

the stress tensor. In addition, we also have a microscopic model:

∂

u = L(u; U). (6.3.2)

Here the macroscale variable U may enter the system as constraint. We may also write

the macro and micro models abstractly as:

F (U, D) = 0 (6.3.3)

f(u, d) = 0, d = d(U)

d is the data needed in order to set up th e microscale model. For example, if the mi-

croscale mo d el is the NV T ensemble of molecular dynamics, d might be the temperature.

(6.3.1) or (6.3.3) represent the knowledge we have about the possible form of the

eﬀective macroscale model.

The general philosophy of HMM is to solve the incomplete macroscale model by

extracting the needed data from the microscale model. The coupling between the macro

and micro models are done in such a way that the macro-state provides th e constraints

for setting up the micro model and the micro model provides the needed constitutive

data D for the macro model. The two main components of HMM are:

1. A macroscopic solver. Based on whatever knowledge that is available on the

macroscale b ehavior of the system, we make an assumption about the form of

the macroscale model, from which we select a suitable macroscale solver. For ex-

ample, if we are dealing with a variational problem, we may use a ﬁnite element

method as the macroscale solver.

6.3. THE HETEROGENEOUS MULTISCALE METHOD 265

U F(U, D) = 0

f(u, d) = 0

compression

reconstruction

constraints

data estimation

Figure 6.3: Schematics of the HMM framework

2. A procedure for estimating the missing macroscale data D using the microscale

model. This is typically done in two steps:

(a) Constrained mi croscale simulation: At each point where some macroscale data

is needed, perform a series of microscopic simulations which are constrained

so that they are consistent with the local value of the macro variable.

(b) Data processing: Use the results from the microscopic simulations to extract

the macroscale data needed in the macroscale solver.

Consider for example the case when the microscale model is an elliptic equation with

multiscale coeﬃcients:

−∇ · (a

(x)∇)u

(x) = f (x) (6.3.4)

Assume that the macroscale model is of the form

−∇ · (A(x)∇)U (x) = f (x) (6.3.5)

Naturally, as the macroscale solver, we choose standard ﬁnite element methods, e.g. the

piecewise linear ﬁnite element method, over a coarse mesh. The data that need to be

estimated is the stiﬀness matrix for the ﬁnite element method. If A = A(x) were known,

we would simply follow standard practice and use numerical quadrature to compute the

elements in the stiﬀness matrix. Since A is not known, we set up a microscale simulation

around each quadrature point in order to estimate the needed function value at that

quadrature point. The details of this procedure will be discussed in Chapter 8.

266 CHAPTER 6. CAPTURING THE MACROSCALE BEHAVIOR

As a second example, consider incompressible polymeric ﬂuid ﬂows for which the

macroscale model is a continuum model for the macroscale velocity ﬁeld U in the form:

(∂

U + (U · ∇)U) = ∇ · σ

∇ · U = 0

These are simply statements of the conservation of momentum and mass, for a ﬂuid of

constant density ρ

. The u nknown data is the stress σ: D = σ.

Let us assume that the micro model is a molecular dynamics model for the particles

that make up the ﬂuid:

= f

, j = 1, 2, ··· , N (6.3.6)

Here m

, y

are respectively the mass and position of the j-th particle, f

is the force

acting on the j-th particle.

Given that the macroscale model is in the form of an incompressible ﬂow equation, it is

natural to select the projection method as the macro-solver [11]. In the implementation of

the projection method, we will need the values of σ at the appropriate grid p oints. These

are the data th at need to be estimated. At this point, we have to make an assumption

about what σ depends on, this enters in the constraints that we put on the microscale

model. Let us assume that

σ = σ(∇U) (6.3.7)

We will constrain the molecular dynamics in such a way that the average strain rate

is given by the value of ∇U at the relevant grid point. In general, implementing such

constraints is the most diﬃcult step in HMM. For the present example, one possible

strategy is discussed in [63].

From the results of molecular dynamics, we need to extract the values of the needed

components of the stress. For t his purpose, we need a formula that expresses stress in

terms of the output of the molecular dynamics. This can be obtained by modifying the

Irving-Kirkwood formula [49]. These details will be explained in Section 6.6.

In summary, three main ingredients are required in HMM:

1. a macro-solver, here the projection method;

2. a micro-solver, here the constrained molecular dynamics;

3. a data estimator, here the modiﬁed Irving-Kirkwood formula.