Semi-Lagrangian Subgrid Reconstruction for Advection-Dominant Multiscale Problems

by   Konrad Simon, et al.
University of Hamburg

We introduce a new framework of numerical multiscale methods for advection-dominated problems motivated by climate sciences. Current numerical multiscale methods (MsFEM) work well on stationary elliptic problems but have difficulties when the model involves dominant lower order terms. Our idea to overcome the assocociated difficulties is a semi-Lagrangian based reconstruction of subgrid variablity into a multiscale basis by solving many local inverse problems. Globally the method looks like a Eulerian method with multiscale stabilized basis. We show example runs in one and two dimensions and a comparison to standard methods to support our ideas and discuss possible extensions to other types of Galerkin methods, higher dimensions and nonlinear problems.



There are no comments yet.


page 19

page 20

page 21


Algebraic Multiscale Method for two–dimensional elliptic problems

We introduce an algebraic multiscale method for two–dimensional problems...

Numerical homogenization for non-linear monotone elliptic problems

In this work we introduce and analyze a new multiscale method for non-li...

Constraint Energy Minimizing Generalized Multiscale Discontinuous Galerkin Method

Numerical simulation of flow problems and wave propagation in heterogene...

A Variational View on Statistical Multiscale Estimation

We present a unifying view on various statistical estimation techniques ...

Bridging the Multiscale Hybrid-Mixed and Multiscale Hybrid High-Order methods

We establish the equivalence between the Multiscale Hybrid-Mixed (MHM) a...

Manifold Learning and Nonlinear Homogenization

We describe an efficient domain decomposition-based framework for nonlin...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

1.1 Motivation and Overview

Simulating complex physical processes at macroscopic coarse scales poses many problems to engineers and scientists. Such simulations strive to reflect the effective behavior of observables involved at large scales even if the processes are partly driven by highly heterogenous micro scale behavior. On the one hand hand resolving the microscopic processes would be the safest choice but such a strategy is prohibitive since it would be computationally expensive. On the other hand microscopic processes significantly influence the macroscopic behavior and can not be neglected.

Incorporating micro scale effects into macro simulations in a mathematically consistent way is a challenging task. There exist many scenarios in different disciplines of science that are faced with such challenges. In fully coupled paleo climate simulations, i.e., climate simulations over more than hundred thousand years a typical grid cell has edge lengths around 200 kilometres and more. Consequently, subgrid processes such as heterogeneously distributed and moving ice shields are not or just insufficiently resolved [32]

. These subgrid processes are usually taken care of by so-called parameterizations. One can imagine this as small micro scale simulations that are then coupled to the prognostic variables such as wind speed, temperature and pressure on the coarse grid (scale of the dynamical core). This coupling from fine to coarse scales is being referred to as upscaling and is unfortunately often done in rather heuristic ways. This leads to wrong macroscopic quantities like wrong pressures and eventually even to wrong wind directions and more undesired effects such as phase errors.

The increasing complexity of earth system models (ESMs), in general, demands for mathematically consistent upscaling of parametrized processes that occur or are modelled on very different scales relative to the dynamical core. A first example was already mentioned: In global climate simulations a sea ice model computes the average ice cover for each coarse cell. This information then enters heat fluxes between ocean and atmosphere since sea ice forms an interface between them and is hence modeled as an averaged diffusive process. It is known from homogenization theory that simple averaging of heterogenous diffusive quantities does not reflect the correct coarse scale diffusion [8].

Another, quite contrary, approach to perfom multiscale simulations is adaptive mesh refinement (AMR), see [7]

. The idea of AMR is to initially run a coarse scale simulation and then assess the quality of the solution locally by means of mathematically or physically based a posteriori estimators. The coarse mesh is then refined in regions where the quality of the solution is not sufficient or coarsened where the solution is smooth. This approach is different from the above upscaling since micro scales are locally resolved – it is therefore sometimes referred to as downscaling. Such an approach is attractive and is being used in practice but it is not feasible in situations in which heterogeneities in the solution are expected to occur globally.

One early attempt to use multiscale methods are Eulerian Lagrangian localized adjoint methods (ELLAM), which constitute a space-time finite element framework, see [9, 20], and form the basis of their multiscale version MsELLAM [36, 10]. Such methods rely on operator splitting and basis functions are required to satisfy an adjoint equation. For a review see [34].

There exist many other multiscale methods. Homogenization is an originally analytical tool to find effective models of otherwise heterogenous models in the limit of a large scale separation [5, 8, 27]. Unfortunately it is often too difficult to find such an effective equation but there exist numerical techniques that aim at designing numerical algorithms to effectively capture the behavior of the solution to the unknown homogenized problem.

The heterogenous multiscale method (HMM) was pioneered by E and Enquist [12, 37] and refers to a rich population of variants [3, 18, 17]. This method emphasizes important principles in the design of multiscale methods such as the choice of macro- and micro-models and their communication. For reviews and further references see [2, 11] and [1] for a discussion of the HMM with special focus on PDEs.

Variational multiscale methods (VMM) have been developed in the 1990s by Hughes and collaborators, see [24, 25]. The spirit of the method lies in a decomposition of the solution space and in the design of variational forms that reflect the relevant scale interactions between the spaces. Many versions of it exist and we point the reader to recent reviews [4, 33] and the vast literature therein.

The present work is inspired by multiscale finite element methods (MsFEM) which can be seen as a special variant of the VMM. The idea of this method is to introduce subgrid variations into basis functions and can be dated back to works by Babuška, Caloz and Osborn [6] and shares ideas with the partition of unity method [31]. The MsFEM in its current form was introduced in [23, 22, 14]. The essential idea of the method is to capture the local asymptotic structure of the solution through adding problem dependent bubble correctors to a standard basis and use these as ansatz and trial functions. Many variations of this method exist and refer the reader to [13, 16] for a review.

Many of the afore mentioned methods have the advantage that they work well for elliptic or parabolic problems and that they are accessible to an analysis. The difficulty in many applications on the other hand is their advection or reaction dominated character, i.e., the dynamics is often driven by low order terms. This poses major difficulties to numerical multiscale methods. Multiscale finite element methods naively applied will not converge to any reasonable solution since basis functions will exhibit artificial boundary layers that are not present in the actual physical flow. Ideas to tackle this problem are based on combining transient multiscale methods with Lagrangian frameworks [35] or with stabilization methods for stationary problems, see [28] for an overview. A HMM based idea for incompressible turbulent flows can be found in [29]. For a method based on the VMM see [30].

1.2 Contribution

Our main contribution is a framework of numerical methods for advection-dominated flows which by reconstructing subgrid variations on local basis functions aims at reflecting the local asymptotic structure of solutions correctly. The idea combines multiscale Galerkin methods with semi-Lagrangian methods by locally solving an inverse problem for the basis representation of solutions that is adapted to the actual flow scenario. We demonstrate the idea on one and two-dimensional advection-diffusion equations with heterogenous background velocities and diffusivities in both non-conservative and conservative form.

Applying standard MsFEMs directly leads to failure since boundary conditions (BCs) for the modified basis functions can not be prescribed arbitrarily. Suitable BCs on the subgrid scale are an essential ingredient for many multiscale methods. A wrong placement of Dirichlet BCs, for example, leads to boundary layers in basis functions that are not there in the real large scale flow, i.e., the way information propagates in advection-dominant flow needs to be respected and not artificially blocked. Other choices of BCs usually complicate the enforcement of conformity or obfuscate additional assumptions on the problem structure (such as local periodicity).

Conformal MsFEM techniques for advection-dominated tracer transport were already explored in our previous work in one spatial dimension on a transient advection-diffusion equation [35]. The finding is that one has to follow a Lagrangian point of view on coarse scales such that flow is “invisible”. On fine scales one can then simplify Lagrangian transforms in order to make advective effects locally milder without going to a fully Lagrangian setting. This amounts to prescribing Dirichlet BCs on coarse flow characteristics and has the effect that basis functions do not develop spurious boundary layers not present in the actual flow. While this work gave some useful insights it is unfortunately not feasible for practical applications since it suffers from several weaknesses. First, it is not directly generalizable to higher dimensions. Secondly, it needs assumptions on the background velocity that are not necessarily fulfilled in practical applications to ensure that coarse scale characteristics do not merge.

In order to circumvent these problems we suggest here a new idea based on a semi-Lagrangian framework that locally in time constructs a multiscale basis on a fixed Eulerian grid. The construction is done in a semi-Lagrangian fashion on the subgrid scale whereas the macroscopic scale is conveniently treated fully Eulerian. This is in complete contrast to our previous work but still respects that information in advection-dominant flows is “mostly” propagated along flow characteristics.

The construction of the basis in each cell at time is done by tracing each Eulerian cell back to time . This yields a distorted cell. A basis on this distorted cell is then reconstructed by solving an inverse problem for the representation of the global solution at the previous time step. Then the local representation of the solution on this cell is propagated forward in time to get a basis on the Eulerian grid at instead of propagating the solution itself.

All reconstructions of the basis in coarse cells are independent from each other and can be performed in parallel. The global time step with the modified basis is not critical since the coarse problem is small and since we can use algebraic constructs to make the assembly of the global problem efficient. Note that although we formulate the algorithm globally in an implicit form it can well be formulated explicitly which will allow for computational efficiency for the global step.

Our new approach has several advantages. First, we can effectively incorporate subgrid behavior of the solution into the multiscale basis via the solution of local inverse problems. Secondly, the method can handle advection-dominated flows in a parallel fashion and, furthermore, the idea works in any dimension. Numerical tests show that it is accurate in both and

since it represents subgrid variability correctly. The method can handle problems that involve an additional reaction term. This consequently includes conservation problems. Furthermore, the idea is generic and can be used for vector valued problems and problems that involve subgrid information coming from actual data.

2 The Semi-Lagrangian Basis reconstruction in One and Two Spatial Dimensions

In this section we will outline our ideas on an advection-diffusion equation (ADE) with periodic boundaries as a model problem. We will outline all ideas for didactical purposed first in and in dimensions on




where is the background velocity, is a matrix-valued diffusivity, and and

are some smooth external forcing and initial condition. We will use bold letters for the vectors and tensors independent of the dimension.

The indices and indicate that both quantities may have large variations on small scales that are not resolved on coarse scales of our multiscale method. We will also work locally on a scale that can resolve the variations in the coefficients. Furthermore, we assume that (see remark below) and that is well-behaved, for example in space and continuous in time. This assumption is often satisfied in practice for example in climate simulations where is given at the nodes of a coarse grid. Depending on the application variations in may be resolved on scale or not. The diffusivity tensor is assumed to be positive definite (uniformly in and point-wise in ) with , and is often derived from parametrized processes such as a varying sea ice distribution, convection, topographical features, or land use patterns. Note, that (1) does not conserve the tracer in contrast to equation (2) provided .

Remark 1

Advection-dominance of a flow (what we sloppily expressed by ) is usually expressed by a dimensionless number – the Péclet number which is essentially the ratio between advective and diffusive time scales. There exist several versions of this number [26]. Since for large variations of the coefficients on the subgrid scale advection-dominance is a very local property we need be more precise with what me mean by that. Here we assume that is high on average, i.e., where is a characteristic length. We take to be the length of the computational domain.

Although the general idea of reconstructing subgrid variability in basis functions is not substantially different for (1) and (2) there are some differences in the propagation of local boundary information described below that strongly influence the accuracy. We will explicitly describe these differences.

Figure 1: Global coarse mesh and local fine mesh in 1D (left) and 2D (right).

We start with outlining our method in one dimension since the idea is very simple and avoids complications that arise in higher dimensions. The idea is to represent non-resolved fine scale variations of the solution locally on a set of non-polynomial basis functions in each coarse cell. That is we fix two (Eulerian) meshes: A coarse mesh of width and on each cell of the coarse mesh we have a fine mesh of width (we assume that does not depend on ). On the coarse mesh we have a multiscale basis , , where is the number of nodes of . This basis depends on space and time and will be constructed so that we will obtain a spatially -conformal (multiscale) finite element space. Note that this is a suitable space for problems (1) and (2) since they have unique solutions with and hence in  [15]. The initial condition can therefore be assumed to be in but in later experiments we will choose it to be smooth. The fine mesh on each cell is used to represent the basis locally, see Figure 1.

Standard MsFEM methods used for porous media flows are designed in such a way that the coarse scale basis solves the PDE model locally with the same boundary conditions as standard FEM basis functions. Note that the global coarse scale MsFEM solution is a linear combination of modified local basis functions that do resolve the fine scale structure induced by the problem’s heterogenities and therefore do represent the asymptotics correctly. This works for stationary elliptic problems and for parabolic problems as long as there is no advective term involved or even dominant. The reason is that advective terms as they appear in our model problems prevent a basis constructed by a standard MsFEM technique to represent the correct asymptotics. This is since flow of information is artificially blocked at coarse cell boundaries. Due to the incorrect boundary conditions for the multiscale basis artificial steep boundary layers form in basis functions that do not exist in the actual global flow. For transient problems another difficulty ist that the local asymptotics around a point depends on the entire domain of dependence of this point and hence subscale information represented by a base function must contain memory of the entire history of the local domain of that base function.

A first attempt to bypass these difficulties was to pose boundary conditions for the basis on suitable space-time curves, across which naturally no information is propagated on the coarse scales. Such curves are Lagrangian paths. Therefore, a Lagrangian framework on the coarse scale and an “almost Lagranian” setting on fine scales was proposed by us earlier [35]. Unfortunately, this method is not feasible since it can not be generalized to higher dimensions and needs restrictive assumptions on the flow field.

Nonetheless, the results of [35] suggest that a coarse numerical splitting of the domain should correspond to a resonable physical splitting of the problem. Instead of a fully Lagrangian method on coarse scales semi-Lagrangian techniques to build the basis can circumvent the difficulties of Lagrangian techniques. But these are only local in time and therefore they do not take into account the entire domain of dependence of a point. We show how to deal with this in the following. First, we start with the global problem.

2.1 The global time step in 1D/2D

Suppose we know a set of multiscale basis functions in a conformal finite element setting, i.e., we approximate the global solution at each time step in a spatially coarse subspace in which the solution is sought (almost everywhere). We denote this finite-dimensional subspace as


First, we expand the solution in terms of the basis at time , i.e.,


Then we test with the modified basis and integrate by parts. Therefore, the spatially discrete version of both problem (1) and (2) becomes the ODE




for (1) and


for (2). The mass matrix is given by


contains forcing and boundary conditions and the initial condition is the projection of onto . Note that (5) contains a derivative of the mass matrix:


This is necessary since the basis functions depend on time and since we discretized in space first. The reader will notice that using Rothe’s method of lines, i.e., discretizing in time first, it is not clear what basis function to use for testing and this a priori leads to a different linear system.

Figure 2: The fine mesh in each cell is traced back one time step where the known solution can be used to reconstruct a basis representation of the solution.

For the time discretization we simply use the implicit Euler method. The discrete ODE then reads


Other time discretization schemes, in particular, explicit schemes are of course possible. For didactic reasons we choose to present the algorithm in an implicit version. We will come back to that later. The next step is to show how to construct the multiscale

Figure 3: 2D illustration similar to the one shown in Figure 2. Each coarse cell is together with its fine meshbeing traced back one time step at which the known global solution can be used to reconstruct a basis.


Convention. In both 1D and 2D we use bold letters like to express potentially vector valued quantities although they would be scalar in one dimension in all formulas and figures. We do so since we would like to avoid confusing the reader with different notations that are merely due to different dimensionality although the basic ideas are the same. Difficulties that arise in more than one dimension will be pointed out explicitly. For the sake of consistency we also use instead of in one dimension. Quantities marked with a tilde like signalize (semi-)Lagrangian quantities.

2.2 The Reconstruction Mesh in 1D/2D

Our idea combines the advantage of both semi-Lagrangian and multiscale methods to account for dominant advection. The reconstruction method is based on the simple observation that local information of the entire domain of dependence is still contained in the global solution at the previous timestep. This can be used to construct an Eulerian multiscale basis: we trace back an Eulerian cell at time on which the solution and the basis are unknown to the previous time step . This gives a distorted cell over which the solution is known but not the multiscale basis , .

In order to find the points where transported information originates we trace back all nodes in from time to . For this one simply needs to solve an ODE with the time-reversed velocity field that reads


for each and then take , see Figures 2 and 3 for an illustration. This procedure is standard in semi-Lagrangian schemes and can be parallelized.

2.3 Basis Reconstruction in 1D

After tracing back each point of to its origin in a distorted coarse cell we need to reconstruct a local representation of the (known) solution on :

Figure 4: Left: Illustration of an oscillatory function (black) being approximated by a standard linear basis (red) on a single cell compared to a modified basis (blue) that solves (13). The regularization parameters were taken as . Right: Comparison of the standard basis to the modified basis. The modified basis does neither constitute a partition of unity not is it positive.

where and are the boundary points of . In one dimension one can of course choose a representation using the standard basis of hat functions but this would not incorporate subgrid information at step at all. Even in the forward advection step (explained below) that is being carried out in the next step subgrid information on the basis will not contain the correct subgrid variability because the information contained in the basis does not take into account the entire domain of dependence of . We choose to bypass this problem by solving an inverse problem for the basis to adjusts the representation. The idea is to fit a linear combination of the basis locally such that is optimally represented, i.e., we solve


The functions are regularizers weighted by positive numbers . A simple regularizer that we found useful in one spatial dimension is a penalization of the deviation of the modified basis function from the standard linear basis function with the same boundary values in the quadratic mean, i.e., we use


where denotes the -th standard (linear) basis on . This amounts to a system of linear equations that can be computed explicitly. In a spatially discrete version this system will be small and cheap to solve. A suitable choice of a regularizer

Figure 5: The basis reconstructed according to (13) at time is propagated forward to time according to (15) or  (16).

depends on the problem at hand. For problems in two dimensions we choose a different one (see below). Figure 4 illustrates the effect of a local reconstruction of a basis compared to a representation with a standard basis.

2.4 Basis Propagation in 1D

After having reconstructed a suitable basis on each coarse cell we have an -conformal basis. This basis, however is a basis at time step and does not live on the coarse Eulerian grid that we initially fixed. The step to take now is to construct a basis at on . This is done similarly to [35], i.e., we evolve the basis according to the model at hand with a vanishing external forcing. Note, however, that we compute the basis at along Lagrangian trajectories starting from , i.e., we need to transform the original model. Equation (1) becomes


transforms into


Note that these evolution equations are solved on , i.e., on the element traced back in time. Advection is “invisible” in these coordinates. The end state on can then be transformed onto the Eulerian element to obtain the desired basis function at the next time step. Corresponding basis functions in neighboring cells can then be glued together to obtain a modified global basis , . This way we get a basis of a subspace of that is neither a partition of unity nor is it necessarily positive. Nonetheless, it is adjusted to the problem and the data at hand. The propagation step is illustrated in Figure 5.

1 input : Problem parameters for (1) or (2), , ,
output : Weights of multiscale solution and set of multiscale basis functions
1 begin
2       Initialize a coarse mesh ;
3       On each cell initialize a fine mesh ;
       /* We need to reconstruct a basis at time . Note that no trace back and no propagation are needed. */
4       for  do in parallel
5             Reconstruct the optimal basis representation of according to (13);
7       end
      /* Now the time stepping starts. We compute the solution at . */
8       for  to  do
9             for  do in parallel
10                   Trace back each node in one time step from to according to (11);
11                   Reconstruct the optimal basis representation of according to (13);
12                   Propagate the optimal basis forward onto according to (15) or (16);
14             end
15            Assemble the global (coarse) system matrices for (10) ;
16             Make a global backward Euler time step using (10);
18       end for
19      Postprocess the solution;
20       Return;
22 end
Algorithm 1 Sketch of reconstruction algorithm for a 1D problem. Implicit version using backward Euler. Note that all loops over cells can be parallelized.

Using our method we reconstruct and advect the representation of the global solution first and then the solution itself using the modified representation. The global step is completely Eulerian while the local reconstruction step is semi-Lagrangian in contrast to [35] where the global step is Lagrangian and and the local step is “almost”-Lagrangian. The essential steps described above are summerized in Algorithm 1. Note that the steps to reconstruct the multiscale basis are embarrassingly parallel and consist of small local problems.

2.5 Basis Reconstruction in 2D

The basis reconstruction step in two dimensions is slightly different. We separate the description of the procedure from the reconstruction in one dimension to make the reader aware of difficulties that arise when increasing the space dimension. Keeping the difficulties in mind one can then easily generalize the procedure even to 3D. We will comment on that again later. It is worth mentioning that these difficulties do not need to arise in non-conformal settings but since we agreed on a conformal method we will be consistent here.

As the reader may have guessed to construct a -conformal basis we need to ensure that the reconstructed global basis is continuous across coarse cell boundaries. This can be achieved by first looping over all coarse edges of the traced back mesh and reconstructing the solution at the previous time step with a basis representation on each edge, i.e., we solve first an inverse problem of the kind (13).


Note that the regularizer (14) needs to be replaced since the edge can be curved (it is unclear what in this case is). We use a harmonic prior for the reconstructed basis


with a low weight in (13). The operator denotes the Laplace-Beltrami operator induced by the standard Laplace operator with the trace topology of the respective (traced back) edge , i.e., is the metric tensor. We pass on providing details here.

The edge reconstructed basis then serves as boundary value for the cell basis reconstruction. The optimization problem to solve on the traced back cell then reads


However, the essential task is to ensure conformity of the global basis by first reconstructing representations on all edges and then inside the cells .

In three dimensions it would be necessary to first reconstruct edge then face and only then the interior representations. This sounds potentially expensive but it is embarrassingly parallel since all reconstructions are independent.

2.6 Basis Propagation in 2D

Again, for didactical reasons we make the reader aware of differences to the propagation of the basis in contrast to one spatial dimension. As in Section 2.5 we need to preserve conformity of the global basis. This can easily be done by fixing the boundary values in the propagation step of the basis functions similar to the propagation in 1D, see Section 2.4.

However, this strategy can result in numerical errors from two sources. Recall that the goal of our multiscale method is to represent the local asymptotic structure of the solution correctly by putting subgrid variations into the basis. This can be interpreted as adding a non-polynomial corrector function in each coarse cell to a standard FEM basis which changes the boundary conditions. These boundary conditions are themselves subject to evolution and keeping them fixed in time will not reflect any system intrinsic dynamics.

The first source of a numerical error is a resonance error. This error usually becomes dominant if the scale of oscillations in the diffusion term is close to the scales resolved by the coarse grid. This is well-documented in the literature on multiscale FEMs for stationary elliptic problems that are usually not dominated by lower order terms. Note, that in practical problems it is often unclear if a scale separation exists or it can not be quantified which makes it difficult to identify a resonance regime [19, 23].

Another source for not capturing the asymptotic structure with fixed boundary values even if there was a scale separation between coarse mesh and diffusion oscillations (which can not be guaranteed) appears in the conservative case (1). This is due to the reactive term


that appears also in (16). This term compensates divergence, i.e., it triggers a reaction on the boundary wherever the divergence of the background velocity does not vanish. Essentially that means that this term is very local but can be quite strong, i.e., it is of the order of if is of order 1 but oscillates on a scale of order . This term must be taken into account when evolving a reconstructed basis.

Fortunately, both errors can be dealt with by employing two strategies. The first one is oversampling, i.e., the local reconstruction and basis propagation on a coarse cell is performed on a slightly larger domain than . After that the basis is being truncated to the relevant domain . Using this strategy is generally possible but results in a non-conformal technique and hence we will not use it here [19]. Instead we will solve a reduced evolution problem on the edges that predicts the boundary values and then evolve the reconstructed interior basis for which we then know the boundary values in time. The edge evolution strategy is in general only necessary in the case of a missing scale separation for the diffusion as explained above or in the conservative case. Note that both errors do not appear in one dimension. We summarize this discussion below.

Strategy without Edge Evolution

In this case we keep the boundary values of the basis that we reconstructed at the previous time step fixed. Similarly to (15) this amounts to solving the evolution problem


for the -th basis function.

1 input : Problem parameters for (1) or (2), , ,
output : Weights of multiscale solution and set of multiscale basis functions
1 begin
2       Initialize a coarse mesh ;
3       On each cell initialize a fine mesh ;
       /* We need to reconstruct a basis at time . Note that no trace back and no propagation are needed. */
4       for  do in parallel
5             Reconstruct the optimal basis representation of according to (13);
7       end
      /* Now the time stepping starts. We compute the solution at . */
8       for  to  do
9             for  do in parallel
10                   Trace back each node in one time step from to according to (11);
11                   Reconstruct the boundary condition of optimal basis representation of according to (17), see Section 2.5;
12                   Reconstruct the optimal basis representation of according to (19), see Section 2.5;
13                   if Problem is not in conservation form then
14                         Propagate the optimal basis forward onto according to (21);
16                   else
17                         Propagate the boundary conditions of the optimal basis forward onto according to (22);
18                         Propagate the optimal basis forward onto according to (23);
20                   end if
22             end
23            Assemble the global (coarse) system matrices for (10) ;
24             Make a global backward Euler time step using (10);
26       end for
27      Postprocess the solution;
28       Return;
30 end
Algorithm 2 Sketch of modified reconstruction algorithm for a 2D problem. Implicit version using backward Euler.

Strategy with Edge Evolution

This method requires a prediction of the boundary values. The prediction is done by solving a reduced evolution problem


on each edge of . The reduced quantities and only contain the tangential parts of their counterparts and . Once this evolution problem is solved for each edge the boundary values for the full evolution problem on are known and one can solve


to get the -th basis.

Figure 6: Illustration of a reconstructed multiscale basis in 2D.

In both cases note again that the evolution equations (21) and (23) are solved on . The end state on can again be mapped onto the Eulerian element to obtain the desired basis function at the next time step.

The described modifications for 2D problems are summarized in Algorithm 2.

3 Numerical Examples in 1D

We will show several 1D examples in a non-conservative and conservative setting according to (1) and (2), respectively. For all 1D tests we use a Gaussian


with variance

centered in the middle of the domain , i.e., . The end time is set to with a time step . We show our semi-Lagrangian multiscale reconstruction method (SLMsR) with a coarse resolution in comparison to a standard FEM with the same resolution and high order quadrature. As a reference method we choose a high-resolution standard FEM with . For the multiscale method we choose a fine mesh with in each coarse cell .

Test 1

We start with four tests showing the capability of the SLMsR to capture subgrid variations correctly. Note that the coarse standard FEM has as many cells the the SLMsR has coarse cells and that it does not capture subgrid variations in the following tests. The resolution for the reference solution resolves all subgrid variations but the reader should keep in mind that practical applications do not allow the application of high-resolution methods. The coefficients




are chosen for the non-conservative equation (1) and the coefficients



Figure 7: Snapshots of the solution at , and . The colored dashed lines show the solution of the standard FEM, the colored line shows the SLMsR. The reference solution is shown in black. (a) Non-conservative equation (1) Coefficients (25). (b) Non-conservative equation (1) Coefficients (26). (c) Conservative equation (2) Coefficients (27). (d) Conservative equation (2) Coefficients (28).

conservative equation (2). The latter one is numerically more difficult when it comes to capturing fine-scale variations, i.e., if then . So if and are of the same order like in one term of a Fourier expansion of then one can expect very steep slopes in the solution. The results of the tests are shown in Figure 7 and the corresponding errors in Figure 8.

Test 2a

It is an important question how the SLMsR behaves in different regimes of data. For the non-conservative case we show two tests.

The first test demonstrates how the SLMsR behaves when all coefficients of the equation are resolved by the SLMsR, i.e., but not by the low resolution FEM that has the same resolution as the SLMsR’s coarse resolution, i.e., . For the coefficients we chose


The different spatial resolutions were taken as while the number of fine cells in each coarse cell was fixed to . The reference solution was computed with and the time step for all methods was taken as .

Error plots of the results are shown in the first row of Figure 9. The plots indeed show that the SLMsR captures the reference solution much more accurately than the standard FEM in both and even for small . Only when decreasing to a resolution that resolves all coefficients the FEM is able to represent the solution correctly and starts converging. The reader should keep in mind that in real world applications the standard FEM will be too expensive to compute while the SLMsR allows for massive parallelization. This is the regime that is of practical interest. Increasing the fine resolution in each coarse cell on the other hand does not improve the the accuracy of the SLMsR. It increases the size of the subgrid problems for the basis functions though.

Test 2b

The second test is to show that the SLMsR essentially starts behaving like a standard FEM if . To demonstrate that we took low frequency coefficients


and refined with the sequence while we we fixed the number of fine cells in each coarse cell to . Note that we start with a coarse resolution regime that almost resolves the data. The error plots for the SLMsR and the FEM indeed indicate that as decreases the SLMsR does not dramatically increase its accuracy while the FEM with linear basis functions increases its accuracy according to the standard estimates. This can be observed in the second row of Figure 9. The relative errors at and for both test problems (29) and (30) are summarized in Tables 2 and 3 in the appendix together with the estimated order of convergence (EOC) in space. We recall that the EOC can be computed as follows: Let be an (error) function that maps a high dimensional vector to a positive real for each . Then the EOC w.r.t. can then be computed as

Figure 8: Relative errors over time in (a) and (b) to the reference solution of the tests shown in Figure 7. The dashed lines show the error of the standard FEM, the full line shows the error of the SLMsR. Color codes are the same as in Figure 7.
Figure 9: Relative errors over time for the unresolved regime of test problem (29) in are shown in (a1) and errors in are shown in (a2). The erros of the numerically resolved regime of problem (30) in are shown in (b1) and errors in are shown in (b2). The respective reference solutions were computed with and for all experiments we used . The dashed lines show the error of the standard FEM, the full line shows the error of the SLMsR. The coarse resolutions were chosen as (yellow), (red), (green), (blue), and (cyan, only second row). All errors are shown in a logarithmic scale.

Similar results in terms of accuracy can be obtained in case of a conservative equation. We pass on showing this here.

Figure 10: Comparison of SLMsR and standard FEM for randomly generated (but fixed) coefficient functions. (a) Snapshots of the solution at and . Solid black lines show the reference solution, dashed red lines show the standard solution and solid blue lines show the SLMsR. (b) Relative error plots for and . (c)

Coefficient functions: the diffusion coefficient (upper plot) was generated from a uniform distribution and scale to have minumum

and maximum . The velocity coeffcient (lower plot) we chose to be a smooth function disturbed by Gaussian noise with mean zero and variance .

Test 3

This test shows an example where both diffusion and background velocity are generated randomly. We intend to show an example of the SLMsR behaves when data is involved that does not exhibit a clear scale separation. For this we initially generate fixed mesh based functions with random nodal coefficients. In each mesh cell the functions are interpolated linearly. Note that this is not to simulate a sampled stochastic process. We simply intend not to create any scale or symmetry bias when constructing coefficient functions. The results look appealing and show a clear advantage of the SLMsR, see Figure 


4 Numerical Examples in 2D

This section is to experimentally demonstrate that the SLMsR can handle conservative and non-conservative advection-diffusion equations with dominant advection term in higher dimensions. All tests are being carried out on the torus (periodic unit square) in the time interval and use the same initial condition. As initial value we chose a normalized super-position of two non-isotropic Gaussians




All tests of the SLMsR are are perfomed on a coarse unstructured uniform triangular Delaunay mesh with coarse cells, i.e., for our triangulation (mean diameter of circumcircle of a cell). We compare the SLMsR to a standard low resolution FEM with the same resolution and to a standard high resolution FEM with approximately cells. To get a fine mesh on each coarse cell of the SLMsR we created a triangulation such that the sum of all fine cells over all coarse cells is

Figure 11: Background velocity for Test 1 and Test 3. Four vortices moving through the domain from left to right and come back to their starting points at .

approximately to get a fair comparison of the SLMsR to the low resolution standard FEM with respect to the reference solution that resolves all coefficients involved.

Test 1

Here we test our multiscale reconstruction method with a solenoidal field described by the stream function


so that


This background velocity desribes four vortices moving in time through the (periodic) domain from left to right and get back to their starting point at . Note that this velocity field involves scales that are resolved by the coarse mesh as well as scales that are not resolved, see Figure 11. Note that since equation (1) and (2) are (analytically) identical and hence we only solve (1).

Figure 12: Snapshots of the solution for Test 1 at time for (a) the low resolution standard FEM, (b) the SLMsR and (c) the reference solution.

The diffusion tensor is chosen to be


Note that in this case advection dominance is a local property and Péclet numbers are ranging from to . Snapshots of the solutions are shown in Figure 12.

Figure 13: Background velocities for Test 2a and Test 2b given by (37) and (38), respectively. Both fields are divergent.

It can be observed that the low resolution FEM does not capture the effective solution well since it diffuses too strongly while the SLMsR reasonably caputeres the effective behavior of the solution and even the fine scale structure. This can be seen also in the error plots, Figure 17.

Test 2a

For this test we use two divergent background velocities and hence we split it into two parts to distiguish between the diffusive transport problem (1) and the conservation law (2). We start with showing how the SLMsR behaves on equation (1). The velocity field is a modification of (35) and is given by


while the diffusion tensor is also given by (36). Here the standard solution does not converge to any reasonable approximation of the effective solution and no valuable understanding of the dynamics can be drawn from it. The SLMsR shows a suprisingly good approximation of the reference solution, see Figure 14.

Test 2b

Here we solve equation (2), i.e., a conservation problem. The divergent velocity field for this test describes regions of fast and slow flow with two separatrices across which there is no flow moving in time from left to right once through the domain during the time interval. The field is given by


Both velocity fields for Test 2a and Test 2b are shown in Figure 13. The diffusion tensor is again given by (