Multiscale finite elements through
advection-induced coordinates for transient
advection-diffusion equations
Konrad Simon∗ and Jörn Behrens∗

arXiv:1802.07684v1 [math.NA] 21 Feb 2018

∗

University of Hamburg, Department of Mathematics,
Grindelberg 5, 20144 Hamburg, Germany

Abstract. Long simulation times in climate sciences typically require coarse grids due to
computational constraints. Nonetheless, unresolved subscale information significantly influences the prognostic variables and can not be neglected for reliable long term simulations. This
is typically done via parametrizations but their coupling to the coarse grid variables often involves simple heuristics. We explore a novel up-scaling approach inspired by multi-scale finite
element methods. These methods are well established in porous media applications, where
mostly stationary or quasi stationary situations prevail. In advection-dominated problems
arising in climate simulations the approach needs to be adjusted. We do so by performing coordinate transforms that make the effect of transport milder in the vicinity of coarse
element boundaries. The idea of our method is quite general and we demonstrate it as a
proof-of-concept on a one-dimensional passive advection-diffusion equation with oscillatory
background velocity and diffusion.

1 Introduction
Motivation and Overview. Geophysical processes in our atmosphere and in the oceans involve
many different spatial and temporal scales. Reliable climate simulations hence need to take into
account the interaction of many scales since these are coupled and influence each other. Resolving
all relevant scales and their interactions poses immense computational requirements. Even on
modern high-performance computers computational constraints force us to make compromises,
i.e., to use a limited grid resolution which then often neglects certain scale interactions or to
model interactions in a heuristic manner.
For example long-term climate simulations as done in the simulation of paleo climate use grid
resolutions of approximately 200km and more. This of course ignores fine scale processes that can
not be resolved by this resolution although they significantly influence prognostic variables on the
coarse grid. Examples include (but are not limited to) moving ice-shields, land-sea boundaries,
flow over rough orography, cloud physics and precipitation. None of these processes is resolved
by the grid and current climate simulations cope with this by using so-called parametrizations.
These can be seen as replacements or simplifications of subgrid processes or processes that are too
complex to be taken into account on the prognostic (coarse) scale. Coupling to the coarse grid is
often done heuristically and can even cause convergence issues when refining the computational
grid.
It is our aim to improve the process of transferring information from the subgrid scale to the
coarse grid in a mathematically consistent way. Numerical multiscale modelling offers promising
mathematical frameworks to achieve this goal. Such methods are already quite established in
other communities such as the porous media community but rarely made their way to climate
simulations. Most mentionable here are methods based on homogenization, the heterogeneous

K. Simon, J. Behrens

MsFEM for transient ADEs

multiscale method (HMM), multiscale methods based on the Eulerian-Lagrangian localized adjoint
method (ELLAM and MsELLAM), variational multiscale methods and multiscale finite element
methods (MsFEM).
Homogenization theory [7, 22, 23] was originally developed as an analytical tool to study the
behavior of solutions to equations with coefficients that rapidly oscillate on fast scales O(ε) in
the limit ε → 0. An equation with effective coefficients is developed that correctly describes
the influence of fast scales as the scale separation grows. The process of incorporating fine scale
information into coarser scales is usually referred to as upscaling. This does, however, not resolve
fine scales. Also analytical results for effective equations are very difficult to obtain and available
only in a limited number of cases.
The HMM introduces a framework for a large range of multiscale problems [12, 11, 1] rather than
a concrete method. In contrast to many other techniques it constitutes a top-down approach based
on blending very different natures of macro and micro models. One simulates a usually incomplete
macro-scale model by choosing an appropriate method (e.g., a finite element method when dealing
with a variational problem). Missing data is incorporated wherever it is needed by performing
constrained micro-scale simulations. The micro-scale data then have to be incorporated into a
properly chosen macro solver. The HMM offers a very general approach but needs a careful choice
of the macro and micro models and their interaction, i.e., the compression and reconstruction of
information.
MsELLAM [10, 25] is based on ELLAM [9, 16, 24] which uses a space time finite element method
(FEM). Basis functions that satisfy the adjoint equation are constructed locally in time, which
leads to the elimination of non-boundary terms in the weak form. Back-tracking characteristics
makes it possible to deal with advection dominated problems but leads to discontinuities in time.
Variational multiscale methods first emerged in [8, 19, 20]. The idea here is to decompose the
solution space into a coarse and a fine scale part and to then separate the usual variational form
into a part that is tested with coarse scale test functions and a part that is tested with fine scale
test functions. This approach, in contrast to the HMM and to homogenization, resolves fine scale
features of the numerical solution. Mixed variational methods have later been introduced in [2, 3],
see also [15]
The MsFEM is closely related to the variational multiscale method and has originally been
introduced in [5, 4] and later been substantially refined in [17, 18]. It relies on the use of nonpolynomial basis functions that reflect the fine scale behavior of the solution. Basis functions are
usually constructed by requiring them to satisfy the equation to leading order. Many variations
have been introduced since then, see [13] for a survey. The method is very attractive since it allows
for massive parallelization and is therefore very scalable.
All the above mentioned methods work well for elliptic problems and stand in contrast to
classical adaptive mesh refinement techniques (AMR) which essentially constitute a down-scaling
approach rather than an upscaling framework by locally resolving “interesting” regions of the
grid [21, 6]. They are proven to work well also for hyperbolic problems. Our list of methods is by
far not complete. We therefore refer in particular to more comprehensive texts on different parts
of multiscale modeling [6, 15, 13, 12, 7] and to references therein.
Contribution. In this work we introduce a generalization of the classical MsFEM introduced
in [17] to suggest an idea for dealing with problems that arise in long-term simulations of climate. The partial differential equations used for these simulations are mostly time-dependent and
are dominated by advection. We therefore demonstrate our idea in a simple setting on a onedimensional advection diffusion equation on a periodic domain that is dominated by the advective
term with fine scale diffusion. In such a scenario the classical MsFEM which is designed for elliptic
problems fails since it relies on a decomposition of the computational domain into a number of
coarse blocks. These coarse blocks are the support of the modified basis functions but blocks are
essentially decoupled from each other. Hence information travelling through the entire domain
is blocked at the coarse block boundaries. This means that boundary conditions of the basis
functions need to be chosen with care. Also the basis functions need to be transient themselves.

2

K. Simon, J. Behrens

MsFEM for transient ADEs

Choosing appropriate boundary conditions for the local problems (here for the basis functions)
is a general problem in multiscale computations. Posing Dirichlet boundary conditions as in the
classical MsFEM on the coarse block boundaries basis functions can develop boundary layers that
do not represent any features that are apparent in the actual flow. We intend to circumvent this
problem not by a different choice of boundary conditions but by posing Dirichlet conditions on
curves that account for the advective part of the equation.
Our key idea is based on the fact that for simple background velocity fields one can transform
an advection-diffusion equation into a pure diffusion equation by following the characteristics of
the advective part, i.e., we move to a Lagrangian setting where we only “see” diffusion across the
characteristics. Such a transform can be done without any danger of crossing characteristics if the
background velocity does not depend on space (it may depend on time though). If the background
velocity depends on space we either average in space and follow the mean flow characteristics or
we follow the characteristics at the coarse grid boundaries and interpolate them inside the blocks.
Each procedure introduces a transform of the equation that reduces the influence of the advective
term.
For a well-behaved background velocity, for example taken from a coarse grid in a climate
simulation, characteristics do not come too close. Inside each coarse block interpolation of the
characteristics emerging at the coarse grid boundary ensures that (interpolated) characteristics
do not cross or come too close. This amounts to a setting that is not fully Lagrangian but still
makes the advective term milder. The transformed equation then is very close to a pure parabolic
setting inside a coarse block and purely parabolic in the vicinity of the boundary of the block.
The latter ensures that no (advective) features cross the boundary and the coarse grid blocks are
not only computationally but also physically decoupled. Posing Dirichlet boundary conditions on
coarse blocks that move with the flow is therefore useful since information flow is not artificially
blocked as with fixed coarse blocks.
From the computational point of view our method is quite attractive. It is composed of two
parts: an offline phase that precomputes the basis functions and an online phase that uses these
basis functions to compute the coarse solution. The overhead of precomputing the basis functions
in each coarse block can further be reduced by parallelization as in the classical elliptic MsFEM.
The online phase is approximately as fast as a low resolution standard FEM but still reveals fine
scale features of a highly resolved solution and is therefore much more accurate than a standard
FEM.
Also note that although we suggest to use the precomputed basis functions in a finite element
framework that the idea is much more generic. Such modified basis functions can potentially be
used in a different global framework that uses a finite volume or discontinuous Galerkin formulation.
We are of course aware of the fact that methods for a simple model like the advection-diffusion
equation do not immediately generalize to more important and harder problems like the primitive
equations. Nevertheless, we strive to suggest ideas for dealing with the interaction of different
parts of such simulations, i.e., between advective and (parametrized) diffusive features that have
a multiscale character. Besides that our method can potentially be useful for simpler problems
dealing with passive tracer transport. We also point out problems with this approach and with
its generalization. At this point we are unfortunately not able to give a profound analysis of the
method but we show test cases that will demonstrate the superiority of our modified MsFEM over
standard methods in certain flow regimes. A rigourous analysis is left for future work.
Outline. This work is organized as follows. In Section 2 we introduce the modified MsFEM
method. This is done systematically starting from an easy setting to the general one. Section 3
shows tests, each of them revealing a different aspect of the modified MsFEM. We conclude with
a discussion of strengths and limitations of our method in Section 4.

3

K. Simon, J. Behrens

MsFEM for transient ADEs

2 Description of the Method
In this section we will give an overview of our method which is based on multiscale finite element
methods. We first introduce the general idea of time-dependent basis functions and formulate
the discrete equations. Since this idea does not work in general we then propose transforms
of the original equation and reformulate the discretized equation. We develop these transforms
hierarchically to provide the reader with a better understanding for the method.
Multiscale finite element methods for elliptic problems have been intensively investigated by the
porous media community [13, 15]. In contrast to the HMM [12] they are designed not only to
represent the coarse scale features of the solution to the problem correctly but also reveal the fine
structure of the solution since the fine scale behavior is resolved by the basis functions. We strive
to adopt this approach to one-dimensional transient advection-diffusion equations of the form
∂t u(x, t) + c(x, t)∂x u(x, t) = ∂x (µ(x, t)∂x u(x, t)) + g(x) ,

(x, t) ∈ I × (0, T ]

t ∈ (0, T ]

u(0, t) = u(1, t) ,

(1)

x∈I

u(x, 0) = f (x) ,

on the unit interval I = [0, 1] in a periodic setting with rapidly oscillating coefficients c, µ ∈
2
L∞ ([0, T ], L∞
per (I)), and smooth periodic f, g ∈ Lper (I). In order for (1) to be well-posed we
assume that µ > 0. Then, by the theory of parabolic equations [14] the weak form of problem (1),
i.e., find u with
d
dt

Z

Z
u(x, t)ϕ(x) dx + c(x, t)∂x u(x, t)ϕ(x) dx
I
I
Z
Z
1
= − µ(x, t)∂x u(x, t)∂x ϕ(x) dx + g(x)ϕ(x) dx ∀ϕ ∈ Hper
(I) , t ∈ (0, T ] , (2)
I

I

1
−1
has a unique solution u that is in L2 ([0, T ], Hper
(I)) and its derivative ∂t u is in L2 ([0, T ], Hper
(I)).
2
By compactness we additionally have u ∈ C([0, T ], Lper (I)) and hence the initial condition u(x, 0) =
f (x) makes sense. We will impose further assumptions on the data in the sequel.

Regime of the Data. Motivated by typical problems that arise in the parametrization of subgrid
processes in climate simulations we assume that µ = µ(x/ε, t) and c = c(x, x/δ, t) with ε  H
and δ  H or δ & H, where H denotes a the magnitude of a resolved numerical scale. By this
we mean that oscillations in the data can not be resolved by H, i.e., one can even think of , δ
as vectors that represent regimes of scales in the data. However, we suppress the , δ notation
unless necessary otherwise. The reader should note that although we keep the idea of possibly fine
oscillatory scales O(δ) in the background velocity the standard case is that δ & H. This is because
background velocity data usually comes from the coarse simulation scale whereas µ comes from
parametrizations and therefore lives on the subgrid scale. However, to be as general as possible
we would like to keep both options, i.e., δ to denote resolved or unresolved scales.
Furthermore, we are interested in a regime that is dominated by the advective term. Since
the coefficients are assumed to be oscillatory it is less enlightening to use a single Péclet number
Pe = kckL∞ L/µ, where L is a characteristic length scale, since in our case c can be is a very local
quantity as well. Instead we are interested in the regime where Pe(x) = c(x, t)L/µ(x, t) is large in
average, i.e., we suggest to look at Pe with a large (normalized) L1 -norm. This would mean that
the tracer distribution over time is on average dominated by advection, i.e., for the background
velocity we assume that it is dominated by its mean.
The Basic Idea of the Method. Multiscale finite elements consist of two main components. A
global formulation and modified basis functions, i.e., a localization. In our approach we decompose
the domain I into a number coarse cells Ki , i = 1, . . . , N − 1 of mesh width H. With each
Eulerian node xm , m = 1, . . . , N , of the coarse mesh we associate a basis function ϕms
m that

4

K. Simon, J. Behrens

MsFEM for transient ADEs

satisfies ϕms
m (xn ) = δmn , i.e., globally we use a conformal finite element formulation. Standard
finite elements interpolate the basis functions between the nodes in a prescribed functional way,
often polynomial. For the multiscale finite element basis functions used in our global formulation,
i.e., the second constituent, we interpolate in such a way that the basis contains information about
the fine scale structure of the problem at hand. In the porous media community the common way
is to require the basis to satisfy the homogeneous equation locally, i.e., in each coarse grid block,
to the highest order. We aim to adopt this idea to transient equations that contain advective
terms which is the major difficulty.
In each coarse block Ki we then seek basis functions ϕli , l = 1, 2, that satisfy
∂t ϕli (x, t) + c(x, t)∂x ϕli (x, t) = ∂x (µ(x, t)∂x ϕli (x, t)) ,
ϕli (0, t)
ϕli (1, t)
ϕli (x, 0)

=
=
=

ϕli,0 (0) ,
ϕli,0 (1) ,
ϕli,0 (x) ,

(x, t) ∈ Ki × (0, T ]
t ∈ (0, T ]
t ∈ (0, T ]

(3)

x ∈ Ki

in the weak sense where ϕli,0 is the l-th P1 -finite element basis function on the coarse scale.
Note that these problems have unique weak solutions ϕli ∈ L2 ([0, T ], ϕli,0 + H01 (I)) with ∂t ϕli ∈
L2 ([0, T ], H −1 (I)) and hence ϕli ∈ C([0, T ], L2 (I)). To numerically solve these problems we employ
the method of lines, i.e., we discretize in space and then solve the resulting ODE. Using a conformal
standard FEM in space the problem reads: Find ϕli ∈ C 1 ([0, T ], V h (Ki )) such that
Z
I

∂t ϕl,h
i (x, t)ψ(x) dx +

Z

c(x, t)∂x ϕl,h
i (x, t)ψ(x) dx
I
Z
h
= − µ(x, t)∂x ϕl,h
i (x, t)∂x ψ(x) dx ∀ψ ∈ V (Ki ) , t ∈ (0, T ] (4)
I

h
l
with the initial condition φl,h
i (x, 0) = PV h (Ki ) ϕi,0 (x) where PV h is the projection onto V (Ki ).
The space V h (Ki ) is given by

V h (Ki ) = ψ ∈ H 1 (Ki ) | ψ is piece-wise polynomial and continuous .
(5)

The nodal basis functions ϕm are then constructed by “gluing” the complementing parts of two
basis functions of two adjacent cells, i.e., the basis functions have the same boundary conditions
at the intersection node like in a standard FEM. For an illustration see Figure 1.
The basis functions need to be pre-computed in an offline-phase. This can be done in parallel for
each cell since the local problems do not depend on each other. In each of the small local problems
actually only one basis function ϕli has to be computed since the other (unique) solution can be

Figure 1: Illustration of basis functions. Standard FEM basis functions ϕi,0 (x, t) = ϕi,0 (x) (left) have constant shape in
time and serve as initial conditions for multiscale FEM basis functions ϕi (x, t) (right) which does depend on time. The
latter is glued together from parts that evolve in each adjacent cell according to the homogeneous equation. It carries the
local structure of the flow in space and time.

5

K. Simon, J. Behrens

MsFEM for transient ADEs

computed as ϕki (x, t) = 1 − ϕli (x, t), k 6= l, due to linearity. Consequently the basis functions fulfill
PN
j=1 ϕj (x, t) = 1, i.e., they constitute a decomposition of unity.
For the global formulation, in contrast to standard finite elements where basis functions do not
depend on time, we seek solutions of the form
uH (x, t) = uH
j (t)ϕj (x, t) .

(6)

Here we used the Einstein’s sum convention. With this ansatz we design our global weak form as:
H
Find uH ∈ Vper
such that
Z

∂t uH (x, t)ϕ(x, t) dx +

I

Z

c(x, t)∂x uH (x, t)ϕ(x, t) dx
I
Z
H
= − µ(x, t)∂x uH (x, t)∂x ϕ(x, t) dx ∀ϕ ∈ Vper
(t) , t ∈ (0, T ] (7)
I

H (0) f where PV H (0) denotes the projection onto the
with the initial condition uH (x, 0) = PVper
per
H
2
H
subspace Vper (0) ⊂ Lper (I). Using the method of lines the global finite element space Vper
can be
formulated as

H
1
H
Vper
= uH ∈ C 1 ([0, T ], Hper
(I)) | uH (t) ∈ Vper
(t)
(8)

with
H
Vper
(t) = spanm {ϕm (·, t) | ϕm satisfies the discretized equation (3) at time t ∈ [0, T ]}

(9)

Note that a basis function ϕm is glued together from two parts that are element of the space
V h (Ki ). This renders the method conformal in space. We use the same constructed basis as basis
of the test function space. However, a Petrov-Galerkin formulation using, e.g., standard P1 -basis
functions is possible.
Using expression (6) as ansatz for the solution in (7) we get
Z
Z

H
ϕi (x, t)∂t uj (t)ϕj (x, t) dx + ϕi (x, t)c(x, t)uH
j (t)∂x ϕj (x, t) dx (10)
I
I
Z
Z
= − ∂x ϕi (x, t)µ(x, t)uH
ϕi (x, t)g(x) dx ∀i = 1, . . . , N , t ∈ (0, T ] . (11)
j (t)∂x ϕj (x, t) dx +
I

I

which can be expressed as an ODE:
d H
u (t) + N (t)uH (t) + A(t)uH (t) = D(t)uH (t) + G(t) ,
(12)
dt
R
H
T
where uH (t) = (uH
1 (t), . . . , uN (t)) , G(t)i = I ϕi (x, t)g(x) dx and
Z
Z
M (t)ij = ϕi (x, t)ϕj (x, t) dx , N (t)ij = ϕi (x, t)∂t ϕj (x, t) dx ,
I
I
Z
Z
(13)
A(t)ij = ϕi (x, t)c(x, t)∂x ϕj (x, t) dx , and D(t)ij = − ∂x ϕi (x, t)µ(x, t)∂x ϕj (x, t) dx
M (t)

I

I

R

and with the initial condition u(0)i = I ϕ(x, 0)f (x) dx. Note that the second term N (t) in (12)
appears since the basis depends on time. Also note that this is a purely Eulerian setting.
Now, one could simply solve this ODE using a suitable time integrator. Unfortunately, this
method fails due to the dominance of the advective term. Simply taking small time steps does
not remedy the problem since the problem lies in the way we compute the solutions to the local
problems. The basis functions are supposed to carry information about the fine scale structure of
the solution but note that the local problems are posed with Dirichlet boundary conditions. Such
type of boundary condition essentially blocks advected information at coarse cell boundaries since

6

K. Simon, J. Behrens

MsFEM for transient ADEs

all coarse cells are decoupled from each other. Hence, for flow dominated by advection multiscale
basis functions constructed according to (3) would develop a boundary layer that is not present in
the global flow. The "correct" choice of type of boundary conditions is indeed a general problem
that one is faced with when solving such local problems [13].
A Modified Idea. We try to circumvent the problem above by moving to a different set of coordinates that make the advective term milder, e.g., we could follow a suitable set of curves emerging
at our coarse grid points as done in a full Lagrangian framework. Then Dirichlet boundary conditions could be posed on these curves since these naturally decouple flow regions. Intuitively this
would give better results since the transform will bring us “closer” to a parabolic setting. However,
a full Lagrangian framework is, in general, difficult to handle since it is not easy to implement
and since characteristics could potentially come very close (or even cross) which requires a special
treatment in such regions. We therefore propose a simpler method.
We start by explaining our idea in a simple setting and then generalize. For this let the background velocity c(x, t) = c(t), i.e., we consider a flow without shear. This special assumption
actually allows us to move to a full Lagrangian setting since we can simply follow characteristics
and know they will not intersect. The new coordinates (ξ, τ ) are then (implicitly) given by
Z τ
c(s) ds
x(ξ, τ ) = ξ +
(14)
0
t(ξ, τ ) = τ .
A little bit more general is the choice
Z τ
hci(s) ds ,
x(ξ, τ ) = ξ +
0

Z

1
where hci(t) =
|I|

c(x, t) dx
I

(15)

t(ξ, τ ) = τ
if c explicitly depends on space. In the new coordinates (15) equation (1) reduces to
∂τ û(ξ, τ ) + c̃(ξ, τ )∂ξ û(ξ, τ ) = ∂ξ (µ̂(ξ, τ )∂ξ û(ξ, τ )) + ĝ(ξ) ,

(ξ, τ ) ∈ I × (0, T ]

τ ∈ (0, T ]

û(0, τ ) = û(1, τ ) ,
û(ξ, 0) = fˆ(ξ) ,

(16)

ξ∈I

where û(ξ, τ ) = u(x(ξ, τ ), t(τ )), respectively ĉ, µ̂, fˆ, ĝ, and
c̃(ξ, τ ) = ĉ(ξ, τ ) − hci(τ ) .

(17)

Note that if c is constant or only depends on time, i.e., we employ (14), then c̃(ξ, τ ) = 0 and (16)
reduces to a simple diffusion equation. An illustration of the two transforms is done in Figure 2(a) and (b).
The subscale problem for the basis functions (3) needs to be transformed as well: Find a (weak)
solution ϕ̂li (ξ, τ ), l = 1, 2, such that
∂τ ϕ̂li (ξ, τ ) + c̃(ξ, τ )∂ξ ϕ̂li (ξ, τ ) = ∂ξ (µ̂(ξ, τ )∂ξ ϕ̂li (ξ, τ )) ,
ϕ̂li (0, τ )
ϕ̂li (1, τ )

=
=

ϕli,0 (0)
ϕli,0 (1)

(ξ, τ ) ∈ K̂i × (0, T ]

,

τ ∈ (0, T ]

,

τ ∈ (0, T ]

ϕ̂li (ξ, 0) = ϕli,0 (ξ) ,

(18)

ξ ∈ K̂i .

A cell K̂i then essentially moves with the flow described by the transform without being deformed.

7

K. Simon, J. Behrens

MsFEM for transient ADEs

A Generalized Modified Idea. In case of a large scale separation between the mean flow and the
oscillatory parts of the background velocity the method using transform (15) delivers good results.
However, as mentioned earlier, in many practically relevant scenarios in climate simulations the
background velocity does not exhibit such a scale separation, i.e., there are mostly (globally
resolved) scales δ & H. We would nevertheless like to keep possible small scale variations in c
to be as general as possible. Variations at a scale O(H) and larger would have the consequence
that the average velocity c̃ in each cell K̂i can potentially be large in contrast to the case where
one has a scale separation. Consequently basis functions following the evolution given by (18)
can still exhibit steep boundary layers that are not apparent in the actual flow. As a remedy
one can make an effort to stay “close” to the full Lagrangian setting locally, i.e., we can track the
characteristics (particle trajectories) of the flow starting at the coarse nodes. Boundary conditions
for the multiscale basis functions can then be posed on these characteristics. In scenarios in which
characteristics do not come too close or diverge too strongly this is, as we will show in numerical
examples, a reasonable strategy. Fortunately, in many situations in climate simulations there is
no situation in which true shocks occur, e.g., wind field data that is taken from a coarse grid does
not behave too badly.
To formalize this strategy let x(ξ, τ ) denote the characteristics emerging from ξ ∈ I. The
evolution is governed by the set of ODEs that reads
d
x(ξ, τ ) = c(x(ξ, τ ), τ ) , x(ξ, 0) = ξ
dτ
d
t(ξ, τ ) = 1 , t(ξ, 0) = 0 .
dτ

(19)

This evolution induces a differentiable transform between (x, t)-coordinates and (ξ, τ )-coordinates
as long as characteristics do not intersect. The latter is the case if c(x, t) is locally Lipschitzcontinuous in x and continuous in t. If we were to use transform (19) everywhere in our domain
(what we do not pursue) this would simply mean a change from Eulerian to Lagrangian variables.
The following computations are formal and only hold as long as trajectories do not cross. As a
first step we will rewrite the equations in characteristic coordinates. For the gradients we have
the relation
 ∂ξ ∂ξ 
∂t
∇(x,t) (·) = ∇(ξ,τ ) (·) ∂x
(20)
∂τ
∂τ
∂x

∂t

∂x
∂ξ
∂t
∂ξ

∂x
∂τ
∂t
∂τ

and
∇(ξ,τ ) (·) = ∇(x,t) (·)

!
(21)

interpolation

(a)

(b)

(c)

Figure 2: Illustration of different transforms. (a) If the background velocity c(x, t) is constant transform (14) gives a pure
diffusion equation. (b) A spatially dependent c is averaged over the entire domain and one follows the characteristics
of the mean flow by transform (15). In the special case c(x, t) = c(t) then the average equals c(t) and the transformed
equation is purely parabolic. (c) Transform (19) follows the true characteristics that emerge from each coarse node ξi and
interpolates them linearly inside the coarse cells.

8

K. Simon, J. Behrens

MsFEM for transient ADEs

Hence, we get
∂ξ
∂τ
+ ∂τ û(ξ, τ )
∂x
∂x
∂ξ
∂τ
∂t u(x, t) = ∂ξ û(ξ, τ )
+ ∂τ û(ξ, τ )
.
∂t
∂t

∂x u(x, t) = ∂ξ û(ξ, τ )

(22)

By means of the inverse function theorem the Jacobian in (21) is (locally) the inverse of the
Jacobian in (20). Using Kramer’s rule the inverse is given by
 ∂t

 ∂ξ ∂ξ 
1
− ∂x
∂τ
∂τ
∂x
∂t
(23)
=
∂x
∂t
∂τ
∂τ
D − ∂ξ
∂ξ
∂x
∂t
where
D=

∂x ∂t
∂x ∂t
−
∂ξ ∂τ
∂τ ∂ξ

(24)

is the determinant of the Jacobian in (21).
Now, using (19) in (22) we see that
  −1

∂x
∂x
û(ξ, τ ) .
∂t u(x, t) + c(x, t)∂x u(x, t) = ∂τ û(ξ, τ ) + ĉ(ξ, τ ) −
∂τ
∂ξ

(25)

Looking again at (19) we see that
∂x
= ĉ(ξ, τ )
(26)
∂τ
and hence the velocity term in the transformed equation would vanish. The advantage of such
a full Lagrangian setting would be that one has then a purely parabolic equation as in the case
when the background velocity is constant or solely depends on time. But transform (19) breaks
down in the presence of shocks. However, true shocks can not occur since the overall equation we
are dealing with is parabolic (although dominated by advection) and strictly speaking there are
no characteristics. We would nonetheless still like to reduce the advective effects and retain the
“nice” part of the equation. This is where our framework starts to differ from the full Lagrangian
setting.
The idea is now, as mentioned earlier, to pose Dirichlet boundary conditions on the characteristics in order to solve subscale problems similar to (18). In order to do this we only need the
coarse cell boundaries to be characteristics of the advective term. Inside each coarse cell we just
need to reduce the advective effects of the original equation. This will then prevent flow across
coarse cell boundaries since near the boundary the transformed equation will be nearly parabolic.
Further, the solution to the equation inside the coarse cell will reflect all the subscale features of
the solution of the original equation.
The diffusive term transforms to
"
 −1 #  −1
∂x
∂x
∂x (µ(x, t)∂x u(x, t)) = ∂ξ µ̂(ξ, τ )∂ξ û(ξ, τ )
.
(27)
∂ξ
∂ξ
Collecting (25) and (27) together with the boundary conditions we are now ready to introduce
the transformed problem for the multiscale basis functions as: By the method of lines find a

9

K. Simon, J. Behrens

MsFEM for transient ADEs

(variational) solution ϕ̂li ∈ C 1 ([0, T ], V h (K̂i )), l = 1, 2, such that
# 
"
 −1
−1
∂x
∂x
l
l
l
∂τ ϕ̂i (ξ, τ ) + c̃(ξ, τ )∂ξ ϕ̂i (ξ, τ ) = ∂ξ µ̂(ξ, τ )
∂ξ ϕ̂i (ξ, τ )
∂ξ
∂ξ
+ ĝ(ξ, τ )
ϕ̂li (0, τ ) = ϕli,0 (0) ,

τ ∈ (0, T ]

ϕ̂li (1, τ ) = ϕli,0 (1) ,

τ ∈ (0, T ]

ϕ̂li (ξ, 0)
where

and

, (ξ, τ ) ∈ K̂i × (0, T ]

=

ϕli,0 (ξ)

(28)

ξ ∈ K̂i

,

  −1

∂x
∂x
c̃(ξ, τ ) = ĉ(ξ, τ ) −
∂τ
∂ξ

(29)

n
o
V h (K̂i ) = ψ ∈ H 1 (K̂i ) | ϕ̂ is piece-wise polynomial and continuous .

(30)

In each of the N −1 coarse cells K̂i = [ξi , ξi+1 ] we now have to compute two additional quantities,
i.e., ∂x/∂ξ and ∂x/∂τ . Computation of the former is based on linear interpolation in space between
the left and right boundary nodes of K̂i which follow the characteristics starting at ξi and ξi+1 .
(more details in the paragraph about the implementation). For the latter, i.e., ∂x/∂τ , we use the
fact that ∂x/∂τ (ξi , τ ) = ĉ(ξi , τ ) at each coarse node. Inside the cells K̂i we can again interpolate
linearly the values ĉ(ξi , τ ) and ĉ(ξi+1 , τ ). Note that now c̃(ξi , τ ) = 0 at all coarse nodes ξi globally
in time. This means that characteristics of the transformed system starting at the coarse nodes
are separatrices of the dynamical system d/dsξ(s) = c̃(ξ(s), s) and hence no flow across coarse cell
boundaries is possible in these coordinates, i.e., the flow is separated by the cells K̂i .
H
such that
The global weak form can be formulated as: Find û(ξ, τ ) ∈ V̂per
Z

∂τ ûH (ξ, τ )ϕ̂H (ξ, τ ) dξ +
I
Z "
=−
I

Z

c̃(ξ, τ )∂ξ uH (ξ, τ )ϕ̂H (ξ, τ ) dξ
 −1 # "
 −1 #
∂x
∂x
µ̂(ξ, τ )∂ξ ûH (ξ, τ )
∂ξ ϕ̂H (ξ, τ )
dξ
∂ξ
∂ξ
Z
H
+ ĝ(ξ, τ )ϕ̂H (ξ, τ ) dξ ∀ϕ̂H ∈ Vper
(τ ) , τ ∈ (0, T ] (31)
I

I

where c̃ is given by (29). The initial condition is ûH (ξ, 0) = PV̂ H (0) f where PV̂ H (0) denotes
per

per

H
H
the projection onto the subspace V̂per
(0) ⊂ L2per (I). The global finite element space V̂per
can be
formulated as
n
o
H
1
H
(32)
V̂per
= û ∈ C 1 ([0, T ], Hper
(I)) | û(τ ) ∈ V̂per
(τ )

with
H
V̂per
(τ ) = spanm {ϕ̂m | ϕ̂m satisfies the discretized equation (28) at time τ ∈ [0, T ]} .

(33)

Note that this formulation includes the ones using the transforms (14) or (15) and that it is not
fully Lagrangian if one uses (19).
The Algorithm and some Notes on the Implementation. In this paragraph we give a brief
overview of steps to be taken for the implementation of our above described method.
Step 1. The first step is to set up a uniform mesh with nodes ξm , m = 1, . . . , N , and cells K̂i ,
i = 1, . . . , N − 1 with mesh size , δ & h. Now one decides for the transform. For each node ξm
either uses (14) or (15) or one solves the ODE (19) on the time interval [0, T ] using a suitable

10

K. Simon, J. Behrens

MsFEM for transient ADEs

solver. For our examples we used an adaptive Runge-Kutta-4/5 solver.
Step 2. One now needs to solve for the basis in an offline phase. This can be done in parallel
since all local problems are independent of each other. As indicated earlier we solve equation (28)
using a conformal standard P1 -finite element formulation. For this we create in each cell K̂i a
uniform mesh with cells k̂ji , i = 1, . . . , Nf − 1 and nodes ξnm , n = 1, . . . , Nf , with ξ1m = ξm and
m
ξN
= ξm+1 . In cell K̂i the discretized ODE for the vector ϕ̂l,H
i (τ ) of degrees of freedom of the
f
basis function ϕ̂l,H
i (τ ) reads
M̂ K̂i

d l,H
l,H
l,H
K̂i l,H
ϕ̂ (τ ) + Â(τ )K̂i ϕ̂l,H
i (τ ) = D̂(τ ) ϕ̂i (τ ) , ϕ̂i (0) = ϕi,0
dτ i

(34)

where ϕl,h
i,0 denotes the nodal values of the standard basis. The system matrices are given by
K̂i
M̂kl
i
Â(τ )K̂
kl

i
D̂(τ )K̂
kl

Z
=

ψkh (ξ)ψlh (ξ) dξ ,

ZK̂i

ψkh (ξ)c̃(ξ, τ )∂ξ ψlh (ξ) dξ ,
K̂i
 −1 #
 −1 # "
Z "
∂x
∂x
h
h
ψl (ξ)
dξ
=−
µ̂(ξ, τ )∂ξ ψk (ξ)
∂ξ
∂ξ
K̂i
=

(35)

where ψkh denote the basis of standard that consists of the common hat functions. Note again that
the computation of A(τ ) and D(τ ) involves terms of the form
∂x
∂τ


and

∂x
∂ξ

−1
.

(36)

Time stepping is done using a second-order explicit scheme for the advective term and an (implicit) Crank-Nicolson scheme for the diffusion. Note that the advection matrix and the diffusion
matrix need to be assembled at each time step. All system matrices in (35) are stored for each
time step (for a reason see Step 3).
In each cell actually only one basis function ϕ̂li needs to be computed. By linearity the other
(unique) solution is then given by ϕki (x, t) = 1 − ϕli (x, t), k 6= l. Hence the constructed basis
functions globally constitute a decomposition of unity.
Step 3. The global system (31) is solved in an online phase. The discretized equations form a
system of ODEs of size (N − 1)-by-(N − 1) (note the periodic boundary conditions):
M̂ (τ )

d H
û (τ ) + N̂ (τ )ûH (τ ) + Â(τ )ûH (τ ) = D̂(τ )ûH (τ ) + Ĝ(τ )
dτ

(37)

ϕ̂H
i (ξ, τ )ĝ(ξ, τ ) dξ,
Z
Z
H
H
M̂ (τ )ij = ϕ̂H
(ξ,
τ
)
ϕ̂
(ξ,
τ
)
dξ
,
N̂
(τ
)
=
ϕ̂H
ij
i
j
i (ξ, τ )∂t ϕ̂j (ξ, τ ) dξ ,
I
I
Z
Z
H
H
H
Â(τ )ij = ϕ̂i (ξ, τ )c̃(ξ, τ )∂ξ ϕ̂j (ξ, τ ) dξ , and D̂(τ )ij = − ∂ξ ϕ̂H
i (ξ, τ )µ(ξ, τ )∂ξ ϕ̂j (ξ, τ ) dξ .

where Ĝ(τ ) =

I

R

I

I

(38)
For the time stepping we again use the same second-order explicit scheme for the advective term
and the (implicit) Crank-Nicolson scheme for the diffusion that we already used for the basis.
Since all the matrices in (38) need to be assembled at each time step solving the global system can
potentially become very expensive. By using a trick looping over all coarse cells and then over all
fine cells can be circumvented since we stored all system matrices in (35). To see this note that

11

K. Simon, J. Behrens

MsFEM for transient ADEs

each global basis function is a linear combination of fine scale (standard) basis functions:
h
ϕ̂H
i (ξ, τ ) = αij (τ )ψ̂j (ξ)

(39)

As a generic example we show this for the advective term. By virtue of (39) we have for the global
element (advection) matrix for the element K̂n
Z

H
Â(τ )nkl =
ϕH
(ξ,
τ
)c̃(ξ,
τ
)∂
ϕ
(ξ,
τ
)
dξ
ξ j
i
K̂
kl
Z n
(40)
h
h
=
αik (τ )ψ̂k (ξ)c̃(ξ, τ )∂ξ αjl (τ )ψ̂l (ξ) dξ
K̂n

= αÂK̂n (τ )αT
where ÂK̂n (τ ) is the stored (advection) matrix of the coarse cell K̂n given by (35). Hence, assembling the global system at each time step breaks down in collecting information from each coarse
block and a simple matrix multiplication. This renders the online phase very fast. Also, to ensure
the unique solvability of the characteristic equation we additionally need to assume that c(x, ·) be
continuous and that c(·, t) be Lipschitz-continuous.

3 Numerical Tests
In this section we demonstrate our method on a number of examples. For all examples presented
here we use the same initial condition. In our case this is a normalized Gaussian wave package,
see Figure 3 which is given by


(x − ν)2
1
(41)
f (x) = √ exp −
2σ 2
σ 2π
All data and all results will be set and shown from an Eulerian point of view. Since no analytical reference solutions are available we will compute a reference solution by a high-resolution
standard FEM (750 elements, first order in
space) in mean flow coordinates (15), i.e., the
advection dominance is already reduced. For
our multiscale FEM we will use N = 10 coarse
elements and Nf = 75 fine elements in each
coarse cell. Basis functions are computed with
a standard FEM. To show the advantage of
our multiscale FEM we will compare it to a
low-resolution standard FEM that uses mean
flow coordinates like the reference solution and
N = 10 elements, i.e., as many elements as the
multiscale method uses on the coarse scale. All
simulations will be carried out on the time intervall [0, T ] with T = 1. We fix these paramFigure 3: Initial condition for all tests is given by (41).
eters unless we explicitly say something else.
750 elements
10 elements

4.0

3.5

3.0

amplitude

2.5

2.0

1.5

1.0

0.5

0.0

0.0

0.2

0.4

0.6

0.8

1.0

x

Case 1. We start with a simple setting. This test is a direct generalization of existing methods
for elliptic equations. We have time dependent coefficients given by
c(t) = 5 cos(10πt) ,
µ(x, t) = 5(t + 1) (0.01 + 0.0099 cos(2πkx)) .

(42)

12

K. Simon, J. Behrens

MsFEM for transient ADEs

L2 (I)
k
k
k
k

= 15
= 30
= 45
= 60

FEM
1.036 · 10−4
1.770 · 10−4
2.054 · 10−4
2.370 · 10−4

MsFEM
1.848 · 10−5
1.437 · 10−5
1.160 · 10−5
1.143 · 10−5

1.0

0.8

0.6
Time

Note that the background velocity only depends on time.
Therefore, characteristics can never cross and hence
transforms (14), (15) and (19) are identical. Snapshots
of the solution and the corresponding errors to the reference solution are shown at four different time stamps
in Figure 5 and error graphs are shown in Figure 6. The
characteristics induced by the background velocity are
shown in Figure 4 and Table 1 summarizes the errors
in L2 (I) and L∞ (I) for a number of different k. Since
the multiscale solution resolves the local structure of the
solution within coarse elements we also show the error
of the derivative, i.e., the H 1 (I)-error. This can also
indirectly be seen in Figure 5 since the error of the multiscale solution is smoother than the one of the standard
solution. The time step was taken to be δt = 10−3 .

0.4

0.2

0.0

0.0

0.2

0.4

0.6

0.8

1.0

X

Figure 4: Characteristics emerging at each coarse
node induced by the velocity coefficient in (42)
with Characteristics originating at ξ = 0 and
ξ = 1 are identical due to periodic boundary conditions.

L∞ (I)
FEM
MsFEM
1.141 · 10−4 2.258 · 10−5
1.809 · 10−4 1.785 · 10−5
2.031 · 10−4 1.460 · 10−5
2.293 · 10−4 1.405 · 10−5

H 1 (I)
FEM
5.622 · 10−4
6.329 · 10−4
6.299 · 10−4
6.523 · 10−4

MsFEM
2.160 · 10−4
1.706 · 10−4
1.556 · 10−4
1.887 · 10−4

Table 1: Relative error of case 1 at T = 1. Coefficients are given by (42).

Case 2. This test is to compare the two MsFEM versions using either transform (15), i.e., the
mean flow transform (MF-MsFEM), or transform (19) which we will refer to as the characteristic
MsFEM (Char-MsFEM). We test the influence of oscillations in the background velocity that
are either in the resolved part of the spectrum or in the unresolved part. By this we mean that
we compare scenarios in which oscillations can or can not be resolved by the coarse grid. The
coefficients are given by
c(x, t) = 10 + cos(2kπx) ,
µ(x, t) = 5(t + 1) (0.01 + 0.0099 cos(60πx)) .

(43)

Here the advective term is on average stronger than the diffusive term by a factor of approximately
5 ∗ 102 (comparing the mean values) and by a factor of 105 in regions where he diffusivity has
a minimum. Oscillations in the background velocity with k < 5 can be resolved on a grid with
only N = 10 elements whereas oscillations of the diffusivity can not be resolved. Figure 7 shows
snapshots of the solution for the case k = 3 and the case k = 60. Errors are shown in Figure 8. In
the first case the Char-MsFEM shows superior performance over both the standard method and
the MF-MsFEM. In the case of a larger frequency separation between the mean velocity and the
oscillatory part one can see that the MF-MsFEM performs slightly better than the Char-MsFEM.
Errors are summarized in Table 2. Note that the interesting scenario for us is the one where k < 5,
i.e., the background velocity can be resolved by the coarse grid since typically in climate simulations
velocities are taken from a coarse mesh and the diffusion parameter includes parametrizations of
processes that are not resolved. In such a scenario our Char-MsFEM outperforms the MF-MsFEM
whereas if k > 5 we observe a similar to slightly better performance of the MF-MsFEM.

13

K. Simon, J. Behrens

MsFEM for transient ADEs

k=3
k = 15
k = 30
k = 60

FEM
4.052 · 10−5
4.080 · 10−5
1.171 · 10−4
4.765 · 10−5

k=3
k = 15
k = 30
k = 60

FEM
5.258 · 10−5
5.071 · 10−5
1.385 · 10−4
5.969 · 10−5

k=3
k = 15
k = 30
k = 60

FEM
5.201 · 10−4
5.185 · 10−4
7.323 · 10−4
5.558 · 10−4

L2 (I)
MF-MsFEM
4.284 · 10−5
3.205 · 10−6
1.634 · 10−5
7.046 · 10−6
L∞ (I)
MF-MsFEM
5.404 · 10−5
6.818 · 10−6
2.332 · 10−5
9.408 · 10−6
H 1 (I)
MF-MsFEM
4.108 · 10−4
2.539 · 10−4
3.467 · 10−4
2.459 · 10−4

Char-MsFEM
1.038 · 10−5
6.787 · 10−6
7.543 · 10−6
2.868 · 10−5
Char-MsFEM
1.743 · 10−5
1.542 · 10−5
1.087 · 10−5
3.863 · 10−5
Char-MsFEM
2.987 · 10−4
2.609 · 10−4
2.665 · 10−4
5.210 · 10−4

Table 2: Relative error of case 2 at T = 1. Coefficients are given by (43).

v=4
v=8
v = 12
v = 16

L2 (I)
FEM
Char-MsFEM
2.688 · 10−4
5.635 · 10−5
−4
1.722 · 10
4.731 · 10−5
−4
1.286 · 10
2.213 · 10−5
−4
9.850 · 10
2.736 · 10−5

L∞ (I)
FEM
Char-MsFEM
2.692 · 10−4
6.779 · 10−5
−4
1.684 · 10
6.077 · 10−5
−4
1.281 · 10
2.848 · 10−5
−4
9.924 · 10
3.874 · 10−5

Table 3: Relative error of case 3 at T = 1. Coefficients are given by (44).

Case 3. Here we demonstrate the advantage of the Char-MsFEM in regimes dominated by advection where we have frequencies in the background velocity that are resolved and that are unresolved
by the coarse grid. We do not show a comparison here to the MF-MsFEM since this method fails
if resolved frequencies are apparent. Furthermore, we use an oscillatory diffusion coefficient as we
did in previous other tests. Coefficients for this test are given by
c(x, t) = v + 1.5 cos(2πx) + 0.5 cos(60πx) ,
µ(x, t) = 5(t + 1) (0.01 + 0.0099 cos(50πx)) .

(44)

Snapshots of the solution and errors for v = 4 are shown in Figure 9 as well as a summary of
errors in Table 3. From the snapshots it is clearly visible that the standard FEM does not capture
the correct large scale behavior, i.e., it is too fast.

14

4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

4

2
1
0
0.2

0.4

0.6

0.8

4

1.0

1

0.0

0.2

0.4

3
2
1
0
0.2

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.4

0.6

0.8

1.0

3
2
1

0.2

0.4

0.6

0.8

absolute error
1.0

0.2

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow

0.0

0.2

0.4

time passed: 0.25 / 1.00

0.6

0.8

1.0

time passed: 0.50 / 1.00

4

4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

1.0
reference
MsFEM, mean flow

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow

0.0

2
1
0

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

0.4

0.6

0.8

4

reference
MsFEM, mean flow

3

multiscale

multiscale

0.8

0

0.0

2
1
0

1.0
reference
MsFEM, mean flow

3
2
1
0

0.0

0.2

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.4

0.6

0.8

1.0

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow

0.0

0.2

0.4

0.6

0.8

absolute error

absolute error

0.6

4

reference
MsFEM, mean flow

multiscale

multiscale

2

0

0.0

absolute error

reference
Char-FEM, Elements: 10, order: 1

3

1.0

0.2

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow

0.0

0.2

0.4

time passed: 0.75 / 1.00

0.6

0.8

1.0

time passed: 1.00 / 1.00

Figure 5: Snapshots of the solution of case 1 with data given by (42) at time stamps t = 0.25, t = 0.5, = 0.75 and T = 1
for k = 30. The upper graphs show a low-resolution standard FEM (red) compared to a high-resolution standard FEM
reference solution (black). The second row of each snapshot shows our MsFEM (blue) compared to the reference solution.
The third row of each snapshot shows the absolute difference of the low-resolution FEM solution (red) and of the MsFEM
solution (blue) to the reference solution. The absolute error of the MsFEM is clearly lower and smoother that the error of
the standard FEM.

Reference solution: Char-FEM, Elements: 750, order: 1

3

Reference solution: Char-FEM, Elements: 750, order: 1

10

10

10

4

5

0.0

0.2

0.4

0.6
time

0.8

1.0

Reference solution: Char-FEM, Elements: 750, order: 1

3

standard FEM (low res)
MsFEM, mean flow

standard FEM (low res)
MsFEM, mean flow

error

error

standard FEM (low res)
MsFEM, mean flow

10

10

10

3

10

4

error

10

4

5

0.0

0.2

0.4

0.6
time

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

time

Figure 6: Errors of the solution of case 1 with data given by (42) for k = 30. Errors of a low-resolution standard FEM
(red) are compared to errors of our MsFEM (blue) with respect to a high-resolution reference solution. The L2 -error and
the L∞ -error of the MsFEM are approximately one order of magnitude better than the errors of the standard FEM. The
error of the derivative of the MsFEM is approximately half an order of magnitude better than the error of the standard
FEM.

4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

4

2
1
0
0.4

0.6

0.8

multiscale

4

1.0

0.0

2
1
0
0.4

0.6

0.8

1
0
0.6

0.8

1.0

0.2

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.2

0
0.4

0.6

0.8

2
1
0.6

0.8

1.0

0.2

0.4

0.6

0.8

0.2

1
0.2

0.4

0.6

1

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

4

2

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.2

0.4

0.6

0.8

1.0

0.0

multiscale

3

0.2

0.4

0.6

0.8

4

reference
MsFEM, mean flow

2
1

1.0
reference
MsFEM, mean flow

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

multiscale

3

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

2
1
0

1.0
reference
MsFEM, characteristics

3
2
1
0

0.2

0.4

0.6

0.8

1.0

0.0

0.8
absolute error

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.6
0.4
0.2
0.2

0.4

0.6

0.8

4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.2
0.0

0.2

0.4

0.6

0.8

1.0

time passed: 0.50 / 1.00

4
standard

reference
Char-FEM, Elements: 10, order: 1

3

0.4

0.4

0.0

1.0

0.2

0.6

time passed: 0.25 / 1.00

standard

1.0

0.2

0

2
1
0

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

3

0.2

0.4

0.6

0.8

4

reference
MsFEM, mean flow

multiscale

multiscale

0.8
reference
MsFEM, characteristics

(a)

4

2
1
0

1.0
reference
MsFEM, mean flow

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

3

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

multiscale

multiscale

0.6

0.4

0

2
1
0

1.0
reference
MsFEM, characteristics

3
2
1
0

0.2

0.4

0.6

0.8

1.0

0.0

0.8

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.6

absolute error

absolute error

0.4

2

0.0

reference
Char-FEM, Elements: 10, order: 1

3

0.4
0.2
0.0

1.0
reference
MsFEM, mean flow

0.6

standard

standard

4

0.0

0.8

3

0.0

1.0

time passed: 0.75 / 1.00

0.0

0.6

1

0.0

0.2

0.8

0.4

2

0.8

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

absolute error

absolute error

0.4

0.4

multiscale

0.2

3

0
0.2

0.6

multiscale

1

0.0

0

absolute error

2

4
multiscale

multiscale

1.0
reference
MsFEM, characteristics

3

0.0

1.0

0
0.2

4

0.0

0.8

reference
Char-FEM, Elements: 10, order: 1

3

0.0

1

0.0

0.6

4
multiscale

multiscale

1.0

2

0.0

0.4

4

reference
MsFEM, mean flow

3

0.8

1.0

0
0.2

4

0.0

0.8
Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.2

0

0.0

0.6

time passed: 0.50 / 1.00

1

0.0

0.4

0.4

0.0

2

0.0

1.0

1

0.0

1.0

reference
Char-FEM, Elements: 10, order: 1

3

0.0

0.8

2

0.6

standard

standard

4

0.8

0.6

reference
MsFEM, characteristics

time passed: 0.25 / 1.00

0.0

0.4

3

0.0

0.2

0.0

1.0
reference
MsFEM, mean flow

1

0.8

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

absolute error

absolute error

0.4

0.4

0.0

0.8

0
0.2

0.6

0.0

0.6

2

0.0

2

0.0

0.4

4
multiscale

multiscale

1.0
reference
MsFEM, characteristics

3

0.8

0.2

3

0
0.2

4

0.0

1

4

reference
MsFEM, mean flow

3

0.0

2

0
0.2

multiscale

0.0

reference
Char-FEM, Elements: 10, order: 1

3

0.2

0.4

0.6

time passed: 0.75 / 1.00

0.8

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, mean flow
MsFEM, characteristics

0.2
0.0

(b)

0.4

0.4

0.0

1.0

0.2

0.6

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

Figure 7: Comparison of snapshots of the solution of case 2 with data given by (43) at time stamps t = 0.25, t = 0.5,
= 0.75 and T = 1 for (a) k = 3 and (b) k = 60. In each image the first row represents a comparison of either a standard
FEM (red), the MF-MsFEM (blue) or a Char-MsFEM (green) to a reference solution (black). The respective absolute
errors are shown in the last row of each image.

Reference solution: Char-FEM, Elements: 750, order: 1

Reference solution: Char-FEM, Elements: 750, order: 1

standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

5

10

4

10

5

error
10

standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

error

10

4

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

time

0.6

0.8

Reference solution: Char-FEM, Elements: 750, order: 1

Reference solution: Char-FEM, Elements: 750, order: 1

standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

5

10

4

10

5

error
10

standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

error

10

4

1.0

time

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

time

0.6

0.8

1.0

time

Reference solution: Char-FEM, Elements: 750, order: 1

Reference solution: Char-FEM, Elements: 750, order: 1
standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

standard FEM (low res)
MsFEM, mean flow
MsFEM, characteristics

error

10 -3

error

10 -3

10 -4
0.0

10 -4
0.2

0.4

0.6
time

(a)

0.8

1.0

0.0

0.2

0.4

0.6

0.8

time

(b)

Figure 8: Errors of the solution of case 2 with data given by (43) for (a) k = 3 and (b) k = 60.

1.0

4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

4

2
1
0
0.0

0.4

0.6

0.8

2
1

0.0

0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1

0.0

0.2

1

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

3

multiscale

multiscale

0.6

4

reference
Char-FEM, Elements: 10, order: 1

0

2
1
0

1.0
reference
MsFEM, characteristics

3
2
1
0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

0.0

absolute error

absolute error

0.4

time passed: 0.50 / 1.00

2

0.0

1.0

2

0.0

1.0

3

0.0

0.8

3

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

standard

standard

4

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.6

reference
MsFEM, characteristics

time passed: 0.25 / 1.00

0.0

0.4

0
0.2

absolute error

absolute error

0.0

1

4

reference
MsFEM, characteristics

multiscale

multiscale

1.0

3

0.0

2

0
0.2

4

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

reference
Char-FEM, Elements: 10, order: 1

3

1.0

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

time passed: 0.75 / 1.00

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

(a)
4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

4

2
1
0
0.0

0.4

0.6

0.8

2
1

0.0

0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1

0.0

0.2

1

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

3

multiscale

multiscale

0.6

4

reference
Char-FEM, Elements: 10, order: 1

0

2
1
0

1.0
reference
MsFEM, characteristics

3
2
1
0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

0.0

absolute error

absolute error

0.4

time passed: 0.50 / 1.00

2

0.0

1.0

2

0.0

1.0

3

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.8

3

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

standard

standard

4

0.0

0.6

reference
MsFEM, characteristics

time passed: 0.25 / 1.00

0.0

0.4

0
0.2

absolute error

absolute error

0.0

1

4

reference
MsFEM, characteristics

multiscale

multiscale

1.0

3

0.0

2

0
0.2

4

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

reference
Char-FEM, Elements: 10, order: 1

3

1.0

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.0

time passed: 0.75 / 1.00

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

(b)
Figure 9: Comparison of snapshots of the solution of case 3 with data given by (44) at time stamps t = 0.25, t = 0.5,
= 0.75 and T = 1 for (a) v = 4 and (b) v = 16. Each image contains a comparison of a standard FEM (red) and the
Char-MsFEM (green) to a reference solution (black). The absolute errors with respective colors are shown in the last row
of each image.

K. Simon, J. Behrens

MsFEM for transient ADEs

Case 4. This test is to demonstrate consistency of the Char-MsFEM in the limit of a large scale
separation. This means that we will demonstrate that our method yields a solution that correctly
follows the reference solution as the diffusion coefficient becomes more and more oscillatory. It is
not our aim to find the limiting equation, i.e., a homogenized equation, and to demonstrate that
our numerical solution enjoys the correct behavior with respect to the solution of the homogenized
equation. Finding the correct homogenized equation is often a difficult task, even in simple
settings such as purely diffusive ones. Additionally to the oscillating diffusive term we will have a
background velocity that is resolved by the coarse grid since our motivation comes from climate
simulations. There the velocity data is usually taken from a coarse grid while the fine scale part
is hidden in the parametrizations. The coefficients that we use for this case are
c(x, t) = (2t + 0.5)(3 + cos(2πx) + cos(60πx)) ,
µ(x, t) = 0.01 + 0.0099 cos(2kπx) .

(45)

Note that the standard solution also runs faster than the Char-MsFEM. This effect is much more
obvious than the wrong representation of the fine scale diffusion on the coarse grid, see Figure 10.
Similarly to case 1 but less obvious the standard solution does not capture the oscillations in the
diffusive term well. This effect is less severe since the velocity coefficient “wipes” these oscillations
out. Table 4 shows the relative deviation of the maximum of the coarse standard solution and of
the Char-MsFEM with respect to the reference solution. Errors in L2 and L∞ are summarized in
Table 5.

k = 15
k = 30
k = 60
k = 90
k = 140
k = 200

FEM
3.66 · 10−2
5.67 · 10−2
7.80 · 10−2
8.52 · 10−2
1.09 · 10−1
1.47 · 10−1

Char-MsFEM
8.94 · 10−3
9.63 · 10−3
2.43 · 10−2
3.08 · 10−2
1.50 · 10−2
1.74 · 10−2

Table 4: Relative deviation of the maximum of case 4 at T = 1. Coefficients are given by (45).

k = 15
k = 30
k = 60
k = 90
k = 140
k = 200

L2 (I)
FEM
Char-MsFEM
2.137 · 10−4
5.418 · 10−5
−4
2.173 · 10
2.233 · 10−5
2.220 · 10−4
4.898 · 10−5
−4
2.240 · 10
4.585 · 10−5
−4
2.308 · 10
4.481 · 10−5
−4
2.420 · 10
4.933 · 10−5

L∞ (I)
FEM
Char-MsFEM
2.204 · 10−4
5.168 · 10−5
−4
2.196 · 10
5.142 · 10−5
2.189 · 10−4
4.731 · 10−5
−4
2.203 · 10
4.755 · 10−5
−4
2.211 · 10
5.028 · 10−5
−4
2.215 · 10
6.083 · 10−5

Table 5: Relative error of case 4 at T = 1. Coefficients are given by (45).

19

4

4

reference
Char-FEM, Elements: 10, order: 1

2
1

0

0.0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

multiscale

1

0.2

0.4

0.6

0.8

1.0

1

0.2

0.5
0.4
0.3
0.2
0.1
0.0

0.2

0.4

0.6

4

0.6

0.4
0.3
0.2

0.8

0.0

1.0

0.0

0.2

0.6

4

reference
Char-FEM, Elements: 10, order: 1

0.8

1.0

reference
Char-FEM, Elements: 10, order: 1

3
standard

standard

0.4

time passed: 0.50 / 1.00

2
1

2
1

0

0

0.0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

1.0

reference
MsFEM, characteristics

3
multiscale

3
multiscale

1.0

0.5

3

2
1
0

2
1
0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7
0.6

0.2

0.4

0.5
0.4
0.3
0.2
0.1

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7
absolute error

absolute error

0.8

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

time passed: 0.25 / 1.00

0.6
0.5
0.4
0.3
0.2
0.1

0.0

0.0

0.2

0.4

0.6

0.8

4

0.0

1.0

time passed: 0.75 / 1.00

0.0

(a)
standard

0.4

0.6

0.8

1.0

reference
Char-FEM, Elements: 10, order: 1

2

0.0

0.2

multiscale

1
0
0.4

0.6

0.8

1.0

1.0

2
1

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.5
0.4
0.3
0.2
0.1

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7
absolute error

0.6

0.6
0.5
0.4
0.3
0.2
0.1

0.2

0.4

0.6

0.8

0.0

1.0

0.0

0.2

time passed: 0.25 / 1.00

4

0.4

0.6

0.8

1.0

time passed: 0.50 / 1.00

4

reference
Char-FEM, Elements: 10, order: 1

reference
Char-FEM, Elements: 10, order: 1

3
standard

3
2
1

2
1

0

0
0.2

0.4

0.6

4

0.8

1.0

0.0

0.2

0.4

0.6

4

reference
MsFEM, characteristics

0.8

1.0

reference
MsFEM, characteristics

3
multiscale

3
2
1
0

2
1
0

0.2

0.4

0.6

0.8

1.0

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7

0.2

0.4

0.5
0.4
0.3
0.2
0.1

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7
absolute error

0.6

0.0

0.8

reference
MsFEM, characteristics

0
0.2

0.7

0.0

0.6

3

2

0.0

0.4

4

reference
MsFEM, characteristics

3

0.0

1.0

0
0.2

4

0.0

0.8

1

0

0.0

0.6

3

1

0.0

0.4

4

reference
Char-FEM, Elements: 10, order: 1

2

0.0

0.2

time passed: 1.00 / 1.00

3
standard

0.6

0.1

0.0

multiscale

0.4

0.7
absolute error

0.6

absolute error

1.0

2

0.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.7

standard

0.8

reference
MsFEM, characteristics

0

0.0

multiscale

0.6

3

2

0

absolute error

0.4

4

reference
MsFEM, characteristics

3
multiscale

2
1

0

absolute error

reference
Char-FEM, Elements: 10, order: 1

3
standard

standard

3

0.6
0.5
0.4
0.3
0.2
0.1

0.2

0.4

0.6

time passed: 0.75 / 1.00

0.8

0.0

1.0

0.0

(b)

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

Figure 10: Comparison of snapshots of the solution of case 4 with data given by (45) at time stamps t = 0.25, t = 0.5,
= 0.75 and T = 1 for (a) k = 15 and (b) k = 140. Each image contains a comparison of a standard FEM (red) and the
Char-MsFEM (green) to a reference solution (black). The absolute errors with respective colors are shown in the last row
of each image.

4

reference
Char-FEM, Elements: 10, order: 1

3

standard

standard

4

2
1

K. Simon, J. Behrens
0

0.0

0.2

0.6

0.8

1
0
0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1

0.0

0.2

1

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

reference
Char-FEM, Elements: 10, order: 1

3
2
1
0

0.2

0.4

0.6

0.8

4

1.0

0.0

0.2

0.4

0.6

0.8

4

reference
MsFEM, characteristics

3

multiscale

multiscale

0.6

4

reference
Char-FEM, Elements: 10, order: 1

0

2
1
0

1.0
reference
MsFEM, characteristics

3
2
1
0

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

0.0

absolute error

absolute error

0.4

time passed: 0.50 / 1.00

2

0.0

1.0

2

0.0

1.0

3

0.0

0.8

3

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

standard

standard

4

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.6

reference
MsFEM, characteristics

time passed: 0.25 / 1.00

0.0

0.4

0
0.2

absolute error

absolute error

MsFEM for transient ADEs
0.2

4

reference
MsFEM, characteristics

2

0.0

1

0.0

multiscale

multiscale

1.0

3

0.0

2

0
0.4

4

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

reference
Char-FEM, Elements: 10, order: 1

3

1.0

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.0

time passed: 0.75 / 1.00

0.2

0.4

0.6

0.8

1.0

Char-FEM, Elements: 750, order: 1
MsFEM, characteristics

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

Figure 11: Comparison of snapshots of the solution of case 5 with data given by (46) at time stamps t = 0.25, t = 0.5,
= 0.75 and T = 1. Each image contains a comparison of a standard FEM (red) and the Char-MsFEM (green) to a reference
solution (black). The absolute errors with respective colors are shown in the last row of each image.

Case 5. We use a simple example to demonstrate that the Char-MsFEM also works with nonvanishing external forcing. The coefficients we choose are given by
g(x) =0.015 sin(8π)
c(x, t) = (2t + 0.5)(3 + cos(2πx) + cos(60πx)) ,

(46)

µ(x, t) = 0.01 + 0.0099 cos(2kπx) .
Snapshots are shown in Figure 11. Remarkably, the Char-MsFEM also yields good results if we
choose the above forcing to act on a scale much finer than resolved by the coarse grid. The
resulting multiscale solution then of course does not capture the fine scale forcing since the basis
functions do not (although they capture the fine scale diffusion) but we observed in tests that it
still runs correctly with the reference solution in contrast to the solution of the standard FEM.
However, we do not strive to investigate external forcing on subgrid scales and therefore pass on
showing this.
A note on convergence. In Case 4 we demonstrated that the Char-MsFEM correctly captures
fine scale features of the solution as the scale separation in the diffusion coefficient grows. Here
we show indications for convergence of the Char-MsFEM as the coarse resolution H → 0. The
classical multiscale finite element method for elliptic problems as described in [17, 18] is shown
to behave like a standard FEM as H → 0. Here we see a similar behavior but the Char-MsFEM
convergences slowly. The data for the tests is given by
α
cos(2πx)) ,
2
µ(x, t) = 0.01 + 0.0099 cos(2kπx) .
c(x, t) = (2t + 0.5)(1.5 +

(47)

21

K. Simon, J. Behrens

MsFEM for transient ADEs

The right-hand side g is taken to be zero. Snapshots for the cases α = 0 and α = 1 for k = 10
are shown for different coarse resolutions in Figure 13 as well as the respective mean square and
maximum errors in Figure 12. The influence of the smoothness of the diffusion coefficient on
the error taken with respect to a high-resolution reference solution is summarized in Table 6 for
spatially constant background velocity as well as for spatially varying velocity. The tests indicate
that the error is smaller the smoother the diffusion coefficient is relative to the coarse grid, i.e.,
the method converges faster for smooth diffusion. Clearly visible from Figure 13 is that once the
coefficients are resolved by the coarse grid the Char-MsFEM solution is slower than the reference
solution but this effect vanishes as H decreases. The time step was taken as δt = 1/1000 and the
number of fine elements in each coarse block was taken to be constant for all different H. However,
note that this regime is practically not relevant since here the standard FEM is sufficient if the
data is resolved by the coarse grid. The advantage of the Char-MsFEM is at scales when data is
not resolved by the coarse grid.

4 Discussion
We have suggested a flexible numerical multiscale method based on a finite element formulation
that uses a set of non-polynomial basis functions to solve an advection-diffusion equation with
a smooth, well-behaved and dominant background velocity and a heterogeneous diffusive term.
The model problem is motivated by problems that arise in long-term climate simulations where
prognostic variables like the background velocity live on a coarse grid while parametrized subgrid
processes are modelled in the diffusive part.
The numerical method is motivated by a classical MsFEM suggested by the porous media
community. We substantially modified the ideas to account for scenarios relevant to climate
simulations. This modification is based on transforms suggested by the background velocity. In
particular these transforms aim to make the dominant advective term milder. Our suggested
multiscale methods are not overly complicated and consist of an offline phase in which many small
local problems for the basis functions (with Dirichlet boundary conditions on curves suggested by
the transform) are solved and an online phase that uses the basis in a finite element like framework
in the global formulation. The offline phase is suitable for massive parallelization where as the
online phase is fast. Also the global form is not restricted to use a finite element framework. It
can potentially be used with finite volumes or discontinuous Galerkin formulations.
In tests we showed an improvement of numerical solutions using our new technique compared to
standard solutions on the same coarse grid. Since no analytical solutions are available we compare
our solutions to high-resolution standard FEM solutions. Tests indicate that our multiscale solutions enjoy the correct upscaling behavior in the limit of fast oscillations in the diffusion, i.e., the
faster the diffusion oscillates relative to the coarse grid the more advantageous are our multiscale
methods over standard methods.
Our multiscale techniques show good results but unfortunately their scope is limited to velocity
fields that have a dominant mean value for all times. This has to do with the dynamics of the
different characteristics. The less dominant the mean of the background velocity the closer characteristics emerging at coarse grid boundaries can come in short time. In particular time-independent
zeros in the background velocity represent stable and unstable attractors of the advective part.
Such features can make coarse blocks collapse and make the method unstable.
Furthermore, a direct generalization to higher dimensions is quite difficult. This is due to
richer dynamics in the characteristics even if the velocity field is solenoidal. Also, our methods
fail if the advection-diffusion equation is in conservation form. The reason is that a conservative
transport term of the form (cu)x decomposes into an advective part cux and a reactive part cx u
that compensates for non-divergence of the background velocity c. Such reaction terms we found
in tests enjoy a different upscaling mechanism and can not be treated with our method.
However, although our method does not work in all cases we showed that it works in relevant
scenarios and identified difficulties in their generalization. Being aware of such difficulties can
facilitate the choice of modified ingredients to account for generalizations and for the development

22

10

4

10

4

10

5

10

6

MSE
10

6

10

7

0.0

0.2

0.4

0.6

0.8

1.0

0.0

time

10

0.2

0.4

0.6

0.8

1.0

time

n=24
n=48
n=96
n=192
n=384

4

10

4

10

5

10

6

n=24
n=48
n=96
n=192
n=384

5

max error

max error

10

n=24
n=48
n=96
n=192
n=384

5

MSE

10

n=24
n=48
n=96
n=192
n=384

10

6

10

7

0.0

0.2

0.4

0.6

0.8

time

1.0

0.0

0.2

0.4

0.6

0.8

1.0

time

(a)

(b)

Figure 12: Relative errors of the solution with data given by (47) for k = 10 over time for (a) α = 0 and (b) α = 1.

k
10

N
24
48
96
192
384
15 24
48
96
192
384
20 24
48
96
192
384

kekL2
2.732 · 10−4
1.898 · 10−4
3.963 · 10−5
4.688 · 10−6
3.939 · 10−7
9.603 · 10−5
3.400 · 10−4
1.141 · 10−4
1.731 · 10−5
1.455 · 10−6
2.629 · 10−5
3.819 · 10−4
1.986 · 10−4
4.110 · 10−5
4.104 · 10−6
(a) α = 0

kekL∞
2.435 · 10−4
1.692 · 10−4
3.569 · 10−5
4.474 · 10−6
4.204 · 10−7
9.129 · 10−5
2.911 · 10−4
9.920 · 10−5
1.566 · 10−5
1.368 · 10−6
2.730 · 10−5
3.250 · 10−4
1.686 · 10−4
3.616 · 10−5
3.702 · 10−6

k
10

N
24
48
96
192
384
15 24
48
96
192
384
20 24
48
96
192
384

kekL2
2.409 · 10−4
3.140 · 10−4
1.783 · 10−4
4.447 · 10−5
1.828 · 10−5
1.647 · 10−4
5.162 · 10−4
4.159 · 10−4
9.032 · 10−5
2.001 · 10−5
7.963 · 10−5
6.185 · 10−4
6.595 · 10−4
1.817 · 10−4
2.930 · 10−5

kekL∞
2.105 · 10−4
3.096 · 10−4
1.888 · 10−4
4.640 · 10−5
1.465 · 10−5
1.335 · 10−4
4.487 · 10−4
4.128 · 10−4
9.610 · 10−5
1.965 · 10−5
6.090 · 10−5
5.368 · 10−4
6.122 · 10−4
1.847 · 10−4
3.160 · 10−5

(b) α = 1

Table 6: Relative errors of the numerical solution for data given by (47) for (a) α = 0 and (b) α = 1 as H becomes small,
i.e., N becomes large. The table shows the influence of increasing frequency of the diffusion coefficient.

reference
MsFEM, mean flow

4

multiscale

multiscale

3

2

1

2

1

0
0.0

reference
MsFEM, characteristics

4

3

0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

MsFEM, mean flow

0.7
0.6
absolute error

absolute error

0.6
0.5
0.4
0.3
0.2

0.0

0.5
0.4
0.3
0.2

0.1
0.0

0.1
0.2

0.4

0.6

0.8

0.0

1.0

0.0

time passed: 1.00 / 1.00

reference
MsFEM, mean flow

0.4

0.6

0.8

multiscale

3

2

1

2

1

0

0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

MsFEM, mean flow

0.7
0.6
absolute error

absolute error

0.6
0.5
0.4
0.3
0.2

0.5
0.4
0.3
0.2

0.1
0.0

0.1
0.2

0.4

0.6

0.8

0.0

1.0

0.0

time passed: 1.00 / 1.00

reference
MsFEM, mean flow

0.4

0.6

0.8

multiscale

3

2

1

2

1

0

0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

MsFEM, mean flow

0.7
0.6
absolute error

absolute error

0.6
0.5
0.4
0.3
0.2

0.5
0.4
0.3
0.2

0.1
0.0

0.1
0.2

0.4

0.6

0.8

0.0

1.0

0.0

time passed: 1.00 / 1.00

reference
MsFEM, mean flow

0.4

0.6

0.8

multiscale

3

2

1

2

1

0

0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

MsFEM, mean flow

0.7
0.6
absolute error

absolute error

0.6
0.5
0.4
0.3
0.2

0.5
0.4
0.3
0.2

0.1
0.0

0.1
0.2

0.4

0.6

0.8

0.0

1.0

0.0

time passed: 1.00 / 1.00

reference
MsFEM, mean flow

0.4

0.6

0.8

multiscale

3

2

1

2

1

0

0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

MsFEM, mean flow

0.7
0.6
absolute error

absolute error

0.6
0.5
0.4
0.3
0.2

0.5
0.4
0.3
0.2

0.1
0.0

1.0

MsFEM, characteristics

0.7

0.0

1.0

reference
MsFEM, characteristics

4

3
multiscale

0.2

time passed: 1.00 / 1.00

4

0.0

1.0

MsFEM, characteristics

0.7

0.0

1.0

reference
MsFEM, characteristics

4

3
multiscale

0.2

time passed: 1.00 / 1.00

4

0.0

1.0

MsFEM, characteristics

0.7

0.0

1.0

reference
MsFEM, characteristics

4

3
multiscale

0.2

time passed: 1.00 / 1.00

4

0.0

1.0

MsFEM, characteristics

0.7

0.0

1.0

reference
MsFEM, characteristics

4

3
multiscale

0.2

time passed: 1.00 / 1.00

4

0.0

1.0

MsFEM, characteristics

0.7

0.1
0.2

0.4

0.6

time passed: 1.00 / 1.00

(a)

0.8

1.0

0.0

0.0

0.2

0.4

0.6

0.8

1.0

time passed: 1.00 / 1.00

(b)

Figure 13: Comparison of snapshots of the solution with data given by (47) at time T = 1 and k = 10 for (a) α = 0 and
(b) α = 1. The multiscale solution is compared to a reference solution (black). The absolute errors are shown in the
second row of each image.

K. Simon, J. Behrens

MsFEM for transient ADEs

of suitable multiscale methods. For example, the reader may think of a combination of semiLagrangian methods to deal with the problem of converging or crossing characteristics and our
multiscale methods. This is left for future research. With this work we hope to spark a debate
on Galerkin based multiscale methods in the climate simulation community and to contribute to
their success.
Acknowledgements. This work was supported by German Federal Ministry of Education and
Research (BMBF) as Research for Sustainability initiative (FONA); www.fona.de through Palmod
project (FKZ: 01LP1513A).

References
[1] Assyr Abdulle, Weinan E, Björn Engquist, and Eric Vanden-Eijnden. The heterogeneous
multiscale method. Acta Numerica, 21:1–87, 2012.
[2] Todd Arbogast. Numerical subgrid upscaling of two-phase flow in porous media. LECTURE
NOTES IN PHYSICS-NEW YORK THEN BERLIN-, pages 35–49, 2000.
[3] Todd Arbogast and Kirsten J Boyd. Subgrid upscaling and mixed multiscale finite elements.
SIAM Journal on Numerical Analysis, 44(3):1150–1171, 2006.
[4] Ivo Babuška, Gabriel Caloz, and John E Osborn. Special finite element methods for a class of
second order elliptic problems with rough coefficients. SIAM Journal on Numerical Analysis,
31(4):945–981, 1994.
[5] Ivo Babuška and John E Osborn. Generalized finite element methods: their performance and
their relation to mixed methods. SIAM Journal on Numerical Analysis, 20(3):510–536, 1983.
[6] Jörn Behrens. Adaptive Atmospheric Modeling. Springer, 2006.
[7] Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. Asymptotic analysis for
periodic structures, volume 374. American Mathematical Soc., 2011.
[8] Franco Brezzi. Interacting with the subgrid world. CHAPMAN AND HALL CRC RESEARCH NOTES IN MATHEMATICS, pages 69–82, 2000.
[9] Michael A Celia, Thomas F Russell, Ismael Herrera, and Richard E Ewing. An eulerianlagrangian localized adjoint method for the advection-diffusion equation. Advances in Water
Resources, 13(4):187–206, 1990.
[10] Aijie Cheng, Kaixin Wang, and Hong Wang. A preliminary study on multiscale ellam schemes
for transient advection-diffusion equations. Numerical Methods for Partial Differential Equations, 26(6):1405–1419, 2010.
[11] Weinan E. Principles of multiscale modeling. Cambridge University Press, 2011.
[12] Weinan E, Bjorn Engquist, et al. The heterognous multiscale methods. Communications in
Mathematical Sciences, 1(1):87–132, 2003.
[13] Yalchin Efendiev and Thomas Y Hou. Multiscale finite element methods: theory and applications, volume 4. Springer Science & Business Media, 2009.
[14] Lawrence C. Evans. Partial differential equations. American Mathematical Society, 2010.
[15] Ivan G Graham, Thomas Y Hou, Omar Lakkis, and Robert Scheichl. Numerical analysis of
multiscale problems, volume 83. Springer Science & Business Media, 2012.

25

K. Simon, J. Behrens

MsFEM for transient ADEs

[16] Ismael Herrera, Richard E Ewing, Michael A Celia, and Thomas F Russell. Eulerianlagrangian localized adjoint method: The theoretical framework. Numerical Methods for
Partial Differential Equations, 9(4):431–457, 1993.
[17] Thomas Y Hou and Xiao-Hui Wu. A multiscale finite element method for elliptic problems
in composite materials and porous media. Journal of computational physics, 134(1):169–189,
1997.
[18] Thomas Y Hou, Xiao-Hui Wu, and Zhiqiang Cai. Convergence of a multiscale finite element
method for elliptic problems with rapidly oscillating coefficients. Mathematics of Computation
of the American Mathematical Society, 68(227):913–943, 1999.
[19] Thomas JR Hughes. Multiscale phenomena: Green’s functions, the dirichlet-to-neumann
formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer
methods in applied mechanics and engineering, 127(1-4):387–401, 1995.
[20] Thomas JR Hughes, Gonzalo R Feijóo, Luca Mazzei, and Jean-Baptiste Quincy. The variational multiscale method – a paradigm for computational mechanics. Computer methods in
applied mechanics and engineering, 166(1-2):3–24, 1998.
[21] Christiane Jablonowski. Adaptive grids in weather and climate modeling. PhD thesis, University of Michigan, 2004.
[22] Vasilii Vasil’evich Jikov, Sergei M Kozlov, and Olga Arsen’evna Oleinik. Homogenization of
differential operators and integral functionals. Springer Science & Business Media, 2012.
[23] Grigoris Pavliotis and Andrew Stuart. Multiscale methods: averaging and homogenization.
Springer Science & Business Media, 2008.
[24] Thomas F Russell and Michael A Celia. An overview of research on eulerian–lagrangian
localized adjoint methods (ellam). Advances in Water Resources, 25(8):1215–1231, 2002.
[25] Hong Wang, Yabin Ding, Kaixin Wang, Richard E Ewing, and Yalchin R Efendiev. A multiscale eulerian–lagrangian localized adjoint method for transient advection–diffusion equations
with oscillatory coefficients. Computing and Visualization in Science, 12(2):63–70, 2009.

26

