A Semi-Lagrangian Multiscale Framework for Advection-Dominant Problems

Abstract

We introduce a new parallelizable numerical multiscale method for advection-dominated problems as they often occur in engineering and geosciences. State of the art multiscale simulation methods work well in situations in which stationary and elliptic scenarios prevail but are prone to fail when the model involves dominant lower order terms which is common in applications. We suggest to overcome the associated difficulties through a reconstruction of subgrid variations into a modified basis by solving many independent (local) inverse problems that are constructed in a semi-Lagrangian step. Globally the method looks like a Eulerian method with multiscale stabilized basis. The method is extensible to other types of Galerkin methods, higher dimensions, nonlinear problems and can potentially work with real data. We provide examples inspired by tracer transport in climate systems in one and two dimensions and numerically compare our method to standard methods.

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 heterogeneous micro scale behavior. On the one hand resolving the microscopic processes would be the safest choice but such a strategy is often 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 can have edge lengths around 200 km 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 parametrizations. 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 many more undesired effects such as phase errors. This demands for mathematically rigorously justified multiscale methods.

Many multiscale methods have been proposed in the past two decades. Among them are the multiscale Eulerian Lagrangian localized adjoint methods (MsELLAM) which constitute a space-time finite element framework, see [9, 10, 20, 35]. Homogenization is an originally analytical tool to find effective models of otherwise heterogeneous models in the limit of a large scale separation [5, 8, 27]. The heterogeneous multiscale method (HMM) was pioneered by Weinan and Enquist [13, 36] and refers to a rich zoo of methods [1,2,3, 12, 18, 19]. Variational multiscale methods (VMM) have been developed in the 1990s by Hughes and collaborators, see [4, 24, 25, 33]. 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.

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]. It shares ideas with the partition of unity method [31]. The MsFEM in its current form was introduced in [15, 22, 23]. 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 we refer the reader to [14, 17] for a review.

Many of the aforementioned 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 first or zero 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 [34] 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]. None of the mentioned methods fully remedies the difficulty of the undesired loss of the localization principle (inherent to advection-dominated problems) that is necessary to obtain weakly coupled problems on the subscale. This is the starting point of our work.

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 [14, 17] with semi-Lagrangian methods [7, 11] 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 heterogeneous background velocities and diffusivities in both non-conservative and conservative form.

Conformal MsFEM techniques for advection-dominated tracer transport were already explored in our previous work [34] in one spatial dimension on a transient advection-diffusion equation. 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 boundary conditions on coarse flow characteristics and has the effect that basis functions do not develop boundary layers that are not there 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 come to close.

In order to circumvent these problems we suggest a new idea based on a semi-Lagrangian framework that locally in time constructs a multiscale basis on a fixed Eulerian grid via many local and independent inverse problems. The construction is done in a semi-Lagrangian fashion on the subgrid scale whereas the macroscopic scale is conveniently treated as completely 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.

All reconstructions of the basis cells are independent from each other and can be performed in parallel while the global time step with the modified basis is numerically cheap. Note that although we formulate the algorithm globally in an implicit form it can well be formulated explicitly which can make the global step more efficient.

Our new approach is generic, conceptionally new, and has several advantages. First, we can effectively incorporate subgrid behavior of the solution into the multiscale basis via the solution of local inverse problems (and possibly even real data). Secondly, the idea works in any dimension. Numerical tests show that it is accurate in both $L^2$ and $H^1$ since it represents subgrid features correctly and can handle problems that involve an additional reaction term. This consequently includes conservation problems.

2 Semi-Lagrangian Multiscale Reconstruction

We will briefly outline and demonstrate our ideas on an advection-diffusion equation (ADE) with periodic boundaries as a model problem. For the sake of brevity we will focus on $d=1$ dimension for the presentation of the method and mention strategies for its generalization to $d=2$ dimensions on

$$\begin{aligned} \begin{aligned} \partial _tu + c_\delta \partial _x u&= \partial _x\left( A_\varepsilon \partial _x u\right) + f \quad \text {in } [0,1]\times [0,T] \\ u(x,0)&= u_0(x) \\ \end{aligned} \end{aligned}$$

(1)

and

$$\begin{aligned} \begin{aligned} \partial _tu + \partial _x \left( c_\delta u\right)&= \partial _x\left( A_\varepsilon \partial _x u\right) + f \quad \text {in } [0,1]\times [0,T] \\ u(x,0)&= u_0(x) \\ \end{aligned} \end{aligned}$$

(2)

where $c_\delta (x,t)$ is the background velocity, $A_\varepsilon (x,t)$ is a diffusivity coefficient, and f and $u_0$ are some smooth external forcing and initial condition.

The indices $\delta >0$ and $\varepsilon >0$ indicate that both quantities may have large variations on small scales that are not resolved on coarse scales $H>0$ of our multiscale method. We will also work locally on a scale $h\ll H$ that can resolve the variations in the coefficients. Furthermore, we assume that $c_\delta \gg A_\varepsilon$ (see remark below) and that $c_\delta$ is well-behaved, for example $C^{1,1}$ in space and continuous in time. This assumption is still practical since, for example, in long term climate simulations $c_\delta$ represents as a relatively smooth coarse grid function. Depending on the application variations in $c_\delta$ may be resolved on scale H or not. The latter is often the case in subsurface flow problems whereas the former is the usual case in tracer transport in climate simulations. The diffusivity which we assume to be $A_\varepsilon \in L_t^\infty L_x^\infty$ often comes from parametrized processes and is furthermore assumed to be positive definite (uniformly in $\varepsilon$ and point-wise in $\mathbf {x}$). Note, that (1) does not conserve the tracer u in contrast to Eq. (2).

Remark. Advection-dominance of a flow (what we sloppily expressed by $c_\delta \gg A_\varepsilon$) is usually expressed by a dimensionless number – the Péclet number $\mathrm {Pe}$ – 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 needto be more precise with what me mean by that. Here we assume that $\mathrm {Pe}$ is high on average, i.e., $\mathrm {Pe}=\frac{\left\| c_\delta \right\| _{L^2}L}{\left\| A_\varepsilon \right\| _{L^2}}$ where L is a characteristic length. We take L to be the length of the computational domain and construct our test examples such that the $\mathrm {Pe}$ is high on average but also such that relevant phenomena will occur far below the reference.

Fig. 1.

Global coarse mesh and local fine mesh in 1D.

We focus on outlining our method in one dimension since the idea is 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 $\mathcal {T}_H$ of width H and on each cell $K\in \mathcal {T}_H$ of the coarse mesh we have a fine mesh $\mathcal {T}_h^K$ of width $h\ll H$. On the coarse mesh we have a multiscale basis $\varphi _i^{H,\mathrm {ms}}:[0,1]\times [0,T]\rightarrow \mathbb {R}$, $i=1,\dots ,N_H$, where $N_H$ is the number of nodes of $\mathcal {T}_H$. This basis depends on space and time and will be constructed so that we will obtain a spatially $H^1$-conformal (multiscale) finite element space. Note that this is a suitable space for problems (1) and (2) since they have unique solutions $u\in L^2([0,T], H^1([0,1]))$ with $\partial _t u\in L^2([0,T], H^{-1}([0,1]))$ and hence in $u\in C([0,T], L^2([0,1]))$ [16]. The initial condition $u^0$ can therefore be assumed to be in $L^2([0,1])$. The fine mesh on each cell $K\in \mathcal {T}_H$ is used to represent the basis locally, see Fig. 1.

Standard MsFEM methods are designed so that the basis functions solve the PDE model locally with appropriate boundary conditions (often Dirichlet). Their choice is crucial. If we now replace the standard basis by these functions the true local asymptotic structure of the global solution is represented since the global solution is a linear combination of modified local basis functions that do reflect the asymptotics. This works for stationary elliptic problems and for parabolic problems as long as there is no advective term involved. The reason is that advective terms like in our model problems prevent a basis constructed by a standard MsFEM to contain the correct asymptotics since flow of information is artificially blocked at coarse cell boundaries. This forms artificial steep boundary layers in basis functions and displays local behavior that is not there in the actual global flow. For transient problems another difficulty is that the local asymptotics around a point depends on the entire domain of dependence of this point and hence there must be some “memory” on the basis.

The results of [34] suggest that a coarse numerical splitting of the domain must correspond to a reasonable physical splitting of the problem. Instead of a fully Lagrangian method on coarse scales a semi-Lagrangian method is used to construct a basis to circumvent the difficulties of pure 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

Suppose we already 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 $V^H(t) \subset H^1([0,1])$ in which the solution u is sought. We denote this finite-dimensional subspace as

(3)

First, we expand the solution $u^H(x,t)$ in terms of the basis at time $t\in [0,T]$, i.e.,

$$\begin{aligned} u^H(x,t) = \sum _{j=0}^{N_H} u_j^H(t)\varphi _j^{H,\mathrm {ms}}(x, t) \,. \end{aligned}$$

(4)

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

$$\begin{aligned} \begin{aligned} \mathbf {M}(t)\frac{\,\mathrm {d}}{\,\mathrm {d}t}\mathbf {u}^{H}(t) + \mathbf {N}(t)\mathbf {u}^{H}(t)&= \mathbf {A}(t)\mathbf {u}^{H}(t) + \mathbf {f}^{H}(t) \\ \mathbf {u}^H(0)&= \mathbf {u}^{H,0} \end{aligned} \end{aligned}$$

(5)

where

$$\begin{aligned} \mathbf {A}_{ij}(t) = \!\!\int _{[0,1]} \varphi _i^{H,\mathrm {ms}}(x, t) \cdot A_\varepsilon (x, t)\partial _x \varphi _j^{H,\mathrm {ms}}(x, t) - \varphi _i^{H,\mathrm {ms}}(x, t)c_\delta (x, t)\partial _x\varphi _j^{H,\mathrm {ms}}(x, t) \,\mathrm {d}x \end{aligned}$$

(6)

for (1) and

$$\begin{aligned} \mathbf {A}_{ij}(t) = \!\!\int _{[0,1]} \varphi _i^{H,\mathrm {ms}}(x, t) \cdot A_\varepsilon (x, t)\partial _x \varphi _j^{H,\mathrm {ms}}(x, t) + \partial _x\varphi _i^{H,\mathrm {ms}}(x, t)c_\delta (x, t) \varphi _j^{H,\mathrm {ms}}(x, t) \,\mathrm {d}x \end{aligned}$$

(7)

for (2). The mass matrix is given by

$$\begin{aligned} \mathbf {M}_{ij}(t) = \int _{[0,1]}\varphi _i^{H,\mathrm {ms}}(x, t)\varphi _j^{H,\mathrm {ms}}(x, t) \,\mathrm {d}x \,, \end{aligned}$$

(8)

$\mathbf {f}^{H}(t)$ contains forcing and boundary conditions and the initial condition $\mathbf {u}^{H,0}$ is the projection of $u^0\in L^2([0,1])$ onto $V^H(0)$. Note that (5) contains a derivative of the mass matrix:

$$\begin{aligned} \mathbf {N}_{ij}(t) = \int _{[0,1]}\varphi _i^{H,\mathrm {ms}}(x, t) \partial _t\varphi _j^{H,\mathrm {ms}}(x, t) \,\mathrm {d}x \,. \end{aligned}$$

(9)

This is necessary since the basis functions depend on time and since we discretized in space first.

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

$$\begin{aligned} \mathbf {M}(t^{n})\mathbf {u}^{n+1} = \mathbf {M}(t^{n})\mathbf {u}^{n} + \delta t \left[ \mathbf {A}(t^{n+1})\mathbf {u}^{n+1} - \mathbf {N}(t^{n+1})\mathbf {u}^{n+1} + \mathbf {f}^{H}(t^{n}) \right] \,. \end{aligned}$$

(10)

Other time discretization schemes, in particular, explicit schemes are of course possible but may involve conditions on the time step size that originate from the (space-time local) transformation to Lagrangian coordinates. We pass on elaborating this here for the sake of brevity. For didactic reasons we therefore choose to present the algorithm in an implicit version. The next step is to show how to construct the multiscale basis.

Convention. Quantities marked with a tilde like $\tilde{x}$ signalize (semi-)Lagrangian quantities.

Fig. 2.

The fine mesh in each cell $K\in \mathcal {T}_H$ is traced back one time step where the known solution can be used to reconstruct a basis representation of the solution.

2.2 The Reconstruction Mesh

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 time step. This can be used to construct an Eulerian multiscale basis: we trace back an Eulerian cell $K\in \mathcal {T}_H$ at time $t^{n+1}$ on which the solution and the basis are unknown to the previous time step $t^n$. This gives a distorted cell $\tilde{K}$ over which the solution $u^n$ is known but not the multiscale basis $\tilde{\varphi }_i$, $i=1,2$.

In order to find the points where transported information originates we trace back all nodes in $\mathcal {T}_H^K$ from time $t^{n+1}$ to $t^{n}$. For this one simply needs to solve an ODE with the time-reversed velocity field that reads

$$\begin{aligned} \begin{aligned} \frac{\,\mathrm {d}}{\,\mathrm {d}t}\tilde{x}_l(t)&= -c_\delta (\tilde{x}_l(t),-t) \,, \quad t\in [-t^{n+1}, -t^n] \\ \tilde{x}_l(-t^{n+1})&= x_l \end{aligned} \end{aligned}$$

(11)

for each $x_l$ and then take $\tilde{x}_l=\tilde{x}_l(-t^n)$, see Fig. 2 for an illustration. This procedure is standard in semi-Lagrangian schemes and can be parallelized.

2.3 Basis Reconstruction

After tracing back each point $x_l$ of $K\in \mathcal {T}_H$ to its origin $\tilde{x}_l$ in a distorted coarse cell $\tilde{K}\in \mathcal {T}_H$ we need to reconstruct a local representation of the (known) solution $u^n$ on $\tilde{K}$:

$$\begin{aligned} \left. u^n(x)\right| _{\tilde{K}} = u^n(\tilde{x}_j)\tilde{\varphi }_{K,1}(x) + u^n(\tilde{x}_{j+1})\tilde{\varphi }_{K,2}(x) \end{aligned}$$

(12)

where $\tilde{x}_j$ and $\tilde{x}_{j+1}$ are the boundary points of $\tilde{K}$. 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 $t_n$ at all. We solve this problem by solving an inverse problem for the basis to modify the local basis representation. The idea is to fit a linear combination of the basis locally such that $u^n$ is optimally represented, i.e., we solve

(13)

Fig. 3.

Left: An oscillatory function (black) is being approximated by a standard linear basis (red) on an interval [a, b] compared to a modified basis (blue) that solves (13). The regularization parameters were taken as $\alpha _i=0.1$. Right: Comparison of the standard basis to the modified basis. The modified basis neither constitutes a partition of unity nor is it positive. (Color figure online)

The operators $\mathcal {R}_i : C^0(\tilde{K})\rightarrow \mathbb {R}$ denote regularizers weighted by positive numbers $\alpha _i\in \mathbb {R}$. A simple regularizer that we found useful in one spatial dimension is a penalization of the quadratic mean deviation of the modified basis function from the standard linear basis function, i.e., we use

$$\begin{aligned} \mathcal {R}_i(\tilde{\varphi }_{\tilde{K},i}) = \left\| \tilde{\varphi }_{\tilde{K},i} - \tilde{\varphi }_{\tilde{K},i}^0 \right\| _{L^2}^2 \end{aligned}$$

(14)

where $\tilde{\varphi }_{\tilde{K},i}$ denotes the t-th standard (linear) basis on $\tilde{K}$. In a spatially discrete version this system is linear and small and will be cheap to solve. A suitable choice of a regularizer depends on the problem at hand. Figure 3 illustrates the effect of a local reconstruction of a basis compared to a representation with a standard basis.

2.4 Basis Propagation

After having reconstructed a suitable basis on each coarse cell $\tilde{K}$ we have an $H^1$-conformal basis. This basis, however is a basis at time step $t^n$ and does not live on the coarse Eulerian grid $\mathcal {T}_H$ that we initially fixed. The step to take now is to construct a basis at $t^{n+1}$ on $\mathcal {T}_H$. This is done by evolving the basis according to the model at hand with a vanishing external forcing. Note, however, that we compute the basis at $t^{n+1}$ along Lagrangian trajectories starting from $t^n$, i.e., we need to transform the original model. Eq. (1) becomes

$$\begin{aligned} \begin{aligned} \partial _t\varphi _{K,i}&= \tilde{\partial }_x\left( \tilde{A}_\varepsilon \tilde{\partial }_x \varphi _{K,i} \right) \quad \text {in } \tilde{K}\times [ t^n,t^{n+1} ] \\ \varphi _{K,i}(\tilde{x}_j,t)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}_j) \,, \quad t\in [ t^n,t^{n+1} ] \\ \varphi _{K,i}(\tilde{x}_{j+1},t)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}_{j+1}) \,, \quad t\in [ t^n,t^{n+1} ] \\ \varphi _{K,i}(\tilde{x},t^n)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}) \end{aligned} \end{aligned}$$

(15)

and Eq. (2) transforms into

$$\begin{aligned} \begin{aligned} \partial _t\varphi _{K,i} + \left( \tilde{\partial }_x \tilde{c}_\delta \right) \varphi _{K,i}&= \tilde{\partial }_x\left( \tilde{A}_\varepsilon \tilde{\partial }_x \varphi _{K,i}\right) \quad \text {in } \tilde{K}\times [t^n,t^{n+1}] \\ \varphi _{K,i}(\tilde{x}_j,t)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}_j) \,, \quad t\in [ t^n,t^{n+1} ] \\ \varphi _{K,i}(\tilde{x}_{j+1},t)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}_{j+1}) \,, \quad t\in [ t^n,t^{n+1} ] \\ \varphi _{K,i}(\tilde{x},t^n)&= \tilde{\varphi }_{\tilde{K},i}(\tilde{x}) \,. \end{aligned} \end{aligned}$$

(16)

These evolution equations are solved on $\tilde{K}$, i.e., on the element $K\in \mathcal {T}_H$ traced back in time. Advection is “invisible” in these coordinates. The end state $\varphi _{K,i}(\tilde{x},t^{n+1})$ on $\tilde{K}$ can then be transformed onto the Eulerian element $K\in \mathcal {T}_H$ to obtain the desired basis function $\varphi _{K,i}(x,t^{n+1})\sim \varphi _{K,i}^{n+1}(x)$ at thenext time step. Corresponding basis functions in neighboring cells can then be glued together to obtain a modified global basis $\varphi _i^{H,\mathrm {ms}}$, $i=1,\dots ,N_H$. This way we get a basis of a subspace of $H^1$ 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 Fig. 4.

Fig. 4.

The basis reconstructed according to (13) at time $t^n$ is propagated forward to time $t^{n+1}$ according to (15) or (16).

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 [34] where the global step is Lagrangian and the local step is “almost”-Lagrangian. Note that the steps to reconstruct the multiscale basis are embarrassingly parallel and all constitute of small problems.

2.5 Basis Reconstruction and Propagation in 2D

The ideas of the above method can be transferred to two dimensions. The reconstruction in two dimensions though is different. We intend to briefly give the reader an idea of the differences without describing the details in order to point out difficulties in the generalization to higher dimensions.

Suppose that we want to reconstruct a basis on cell $K\in \mathcal {T}_H$ at time $t^{n+1}$. We trace back the cell as described in (11) and denote the distorted cell at time $t^n$ by $\tilde{K}$ and its edges by $\tilde{\varGamma }$. As in one spatial dimension, to construct a $H^1$-conformal basis we need to ensure continuity of the basis across coarse cell boundaries. This can be achieved by first reconstructing the solution at the previous time step $t^n$ with a basis representation on each edge, i.e., we solve first an inverse problem on each deformed edge $\tilde{\varGamma }$ similar to (13).

Note that the regularizer (14) needs to be replaced since the edge $\tilde{\varGamma }$ is usually curved. We use a harmonic prior

$$\begin{aligned} \mathcal {R}_i(\tilde{\varphi }_{\tilde{K},i}) = \left\| \varDelta _{g(\tilde{\varGamma })}\tilde{\varphi }_{\tilde{\varGamma },i} \right\| _{L^2(\tilde{\varGamma })}^2 \end{aligned}$$

(17)

with a low weight $\alpha _i$ as in (13). The operator $\varDelta _{g(\tilde{\varGamma )}}$ denotes the Laplace-Beltrami operator induced by the standard Laplace operator with the trace topology of the respective edge $\tilde{\varGamma }$, i.e., $g(\tilde{\varGamma })$ is the metric tensor. We pass on providing details here.

The edge reconstruction provides boundary values for the cell basis reconstruction. The optimization problem to solve on $\tilde{K}$ can then again be designed similar to (13) just constrained by the previously reconstructed boundary values. The essential task is to ensure conformity of the global basis by first reconstructing representations on all edges $\tilde{\varGamma }$ of $\tilde{K}$ and then inside the cells $\tilde{K}$.

In three dimensions it would be necessary to first reconstruct edges then faces and only then the interior representations. This might seem expensive but as in one dimension it is embarrassingly parallel since all reconstructions are independent and the local problems are small.

The next step is the basis propagation. The observation here is that one needs to distinguish between the conservative form (2) and the non-conservative form (1) since in the conservative form an additional local reaction term is responsible for strong local variations. Hence, reconstructed edge boundary values can vary quite strongly in a propagation step similarly to the one described in Sect. 2.4. Consequently, edge boundary values need to be adjusted in the propagation. This is done by first propagating the reconstructed edge boundary values and then using these as (time-dependent) boundary values for propagation problems similar to (15) or (16). We pass on providing technical details.

3 Numerical Examples

For all 1D tests we use a Gaussian

$$\begin{aligned} u_0(x) = \frac{1}{\sigma \sqrt{2\pi }}\exp {-\frac{(x-\mu )^2}{2\sigma ^2}} \,. \end{aligned}$$

(18)

with variance $\sigma =0.1$ centered in the middle of the domain [0, 1], i.e., $\mu =0.5$. The end time is set to $T=1$ with a time step $\delta t = 1/300$. We show our semi-Lagrangian multiscale reconstruction method (SLMsR) with a coarse resolution $H=10^{-1}$ in comparison to a standard FEM with the same resolution and high order quadrature. As a reference we choose a high-resolution standard FEM with $h_{\mathrm {ref}}=10^{-3}$. For the multiscale method we choose a fine mesh $\mathcal {T}_h^K$ with $h=10^{-2}$ in each coarse cell $K\in \mathcal {T}_H$. We would like to point out that there are no standardized test cases for our type of model. Therefore, we designed our tests in such a way that small scale effects occur in the solution below the coarse resolution and compare our multiscale methods to the performance of standard methods with the same (coarse) resolution.

Fig. 5.

Snapshots of the solution at $t=1/3$, $t=2/3$ and $T=1$. The colored dashed lines show the solution of the standard FEM (10 elements), the colored line shows the SLMsR (10 coarse elements). The reference solution is shown in black (1K elements). (a) non-conservative equation (1) coefficients (19). (b) conservative equation (2) coefficients (20). (Color figure online)

Test 1. We will show two examples in a non-conservative and conservative setting according to (1) and (2), respectively, to show the capability of the SLMsR to capture subgrid variations correctly. Note that the coarse standard FEM has as many cells as the SLMsR has coarse cells and that the standard FEM does not capture subgrid variations in the following tests which can result in aliasing and phase errors. 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

$$\begin{aligned} \begin{aligned} c_\delta (x,t)&= \frac{1}{2}\cos (2\pi t) + \frac{1}{4}\cos (6\pi t)\cos (8\pi x) + \frac{1}{8}\cos (4\pi t)\cos (62\pi x) + \frac{1}{8}\cos (150\pi x) \\ A_\varepsilon (x,t)&= 10^{-3} + 9\cdot 10^{-4}\cos (10\pi t) \cos (86\pi x) \\ \end{aligned} \end{aligned}$$

(19)

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

$$\begin{aligned} \begin{aligned} c_\delta (x,t)&= \frac{1}{2} + \frac{1}{8}\cos (8\pi x) + \frac{1}{8}\cos (62\pi x) + \frac{1}{8}\cos (150\pi x) \\ A_\varepsilon (x,t)&= 10^{-2} + 9\cdot 10^{-3}\cos (10\pi t) \cos (86\pi x) \\ \end{aligned} \end{aligned}$$

(20)

for the conservative equation (2). The latter one is numerically more difficult when it comes to capturing fine-scale variations, i.e., if $c_\delta (x)\sim f(x/\delta )$ then $\frac{\,\mathrm {d}}{\,\mathrm {d}x}c_\delta (x)\sim \delta ^{-1}f'(x/\delta )$ and one can expect very steep slopes in the solution. The results of the tests are shown in Fig. 5. The corresponding errors in Table 1 show clearly the superior performance of the SLMsR in regimes of low coarse resolution while it performs similarly to a standard FEM as subgrid variations are resolved by $\mathcal {T}_H$. The reader may observe that the $L^2$-error of the multiscale method is superior to the standard method in the pre-asymptotic regime as well as the $H^1$-error. The latter may increase slightly with growing coarse resolution (but stays in the same order) due to a well-known resonance effect that can occur when physical scales and the coarse numerical scale are of the same order. This can be taken care of by other methods [14] such as oversampling.

Table 1. Relative errors of standard FEM and SLMsR for the non-conservative test problem (19) (left table) and for the conservative problem (20) (right table) at final time $T=1$.

Test 2. This test shows an example where both diffusion and background velocity are generated randomly. We intend to show an example of how the SLMsR behaves when data is involved that does not exhibit a clear scale separation which is a common situation in practice. 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 Fig. 6.

Fig. 6.

Comparison of SLMsR and standard FEM for randomly generated (but fixed) coefficient functions. (a) Snapshots of the solution at $t=1/2$ and $t=1$. Solid black lines show the reference solution, dashed red lines show the standard solution and solid blue solid lines show the SLMsR. (b) The velocity coeffcient we chose to be a smooth function disturbed by Gaussian noise with mean zero and variance 0.1. (c) The diffusion coefficient was generated from a uniform distribution and scaled to have minumum $10^{-5}$ and maximum $10^{-2}$. (Color figure online)

Test 3. Here we show one preliminary example of our SLMsR equation (1) with a dominant advection term in two dimensions. The test was carried out on the torus $\mathbb {T}^2$ (periodic unit square) in the time interval $t\in [0,1]$. As initial value we chose a normalized super-position of two isotropic Gaussians

$$\begin{aligned} u_0(\mathbf {x}) = \frac{1}{2\sqrt{(2\pi )^2\det (\mathbf {M})}} \sum _{i=1}^2 \exp \left\{ -\frac{1}{2}(\mathbf {x}-\varvec{\mu }_i)^T\mathbf {M}^{-1}(\mathbf {x}-\varvec{\mu }_i) \right\} \end{aligned}$$

(21)

where

$$\begin{aligned} \mathbf {M} = \begin{bmatrix} \frac{3}{100}&0 \\ 0&\frac{3}{100} \\ \end{bmatrix} \quad \text {and} \quad \varvec{\mu }_i = \begin{bmatrix} \frac{i}{3} \,,&\frac{1}{2} \end{bmatrix}^T \,. \end{aligned}$$

(22)

The test of the SLMsR was performed on a coarse unstructured uniform triangular Delaunay mesh with $n_c=62$ coarse cells, i.e., for our triangulation $H\sim 0.3$ (maximum 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 $n_f=63K$ 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 approximately $n_f$ 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.

Fig. 7.

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 $T=1$.

We test our multiscale reconstruction method with a solenoidal field $\mathbf {c}_\delta$ described by the stream function

$$\begin{aligned} \psi (\mathbf {x}, t) = \sin (2\pi (x_1-t))\sin (2\pi x_2) \end{aligned}$$

(23)

so that $\mathbf {c}_{\delta }(\mathbf {x}, t) = \nabla ^T\psi$.

This background velocity describes four vortices moving in time through the (periodic) domain from left to right and get back to their starting point at $T=1$. Note that this velocity field involves both scales that are resolved by the coarse mesh and scales that are not resolved, see Fig. 7. Also note that since $\nabla \cdot \mathbf {c}_\delta =0$ Eqs. (1) and (2) are (analytically) identical and hence we only solve (1).

The diffusion tensor is chosen to be

$$\begin{aligned} \mathbf {A}_\varepsilon (\mathbf {x}, t) = \frac{1}{100} \begin{bmatrix} 1-0.9999\sin (60\pi x_1)&0 \\ 0&1-0.9999\sin (60\pi x_2) \\ \end{bmatrix} \,. \end{aligned}$$

(24)

In this case advection dominance is a local property and Péclet numbers are ranging from $\mathrm {Pe}= 0$ to $\mathrm {Pe}\sim 6\cdot 10^6$. Snapshots of the solutions are shown in Fig. 8. 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 captures the effective behavior of the solution and even the fine scale structure.

4 Summary and Discussion

In this work we introduced a new idea for a framework of multiscale methods demonstrated on advection-diffusion equations that are dominated by the advective term. Such a methods are of importance, for example, in reservoir modeling and tracer transport in earth system models. The main obstacles in these applications are, first, the advection-dominance and, secondly, the multiscale characterof the background velocity and the diffusion tensor. The latter makes it impossible to simulate with standard methods due to computational constraints while simulating using standard methods with lower resolution that does not resolve variations in the coefficients leads to incorrect solutions on coarse meshes.

Fig. 8.

Snapshots of the solution (surface plot) for Test 3 at time $T=1$ for the low resolution standard FEM (62 elements), the SLMsR (62 coarse elements) and the reference solution (63 K elements).

Our idea to cope with these difficulties is inspired by ideas for semi-Lagrangian methods, ideas based on “convected fluid microstucture” as described in [21], inverse problems and multiscale finite elements [14]. At each time step we reconstruct fine scale information from the solution at the previous time step. This fine scale information enters the local representation of the solution in each coarse cell, i.e., it is added as a corrector to the local basis such that the basis representation is optimal in some sense. The reconstruction is done by solving an inverse problem with a suitable regularizer and constructs a basis that does not constitute a partition of unity (PoU) and that is made for the concrete problem at hand. The idea of adding prior knowledge about the solution to a local representation in PoU methods, however, is similar, see for example [18, 31]. After reconstructing the basis at the previous time step the basis is evolved with suitable boundary conditions to the time step the basis is sought for, i.e., we evolve the local representation of the solution rather than the solution itself. Note that the global framework of the SLMsR is completely Eulerian while only the local reconstruction step in each coarse cell is semi-Lagrangian.

One of the main features of the SLMsR is its scalability: Although it sounds expensive to trace back each coarse cell, then solve an inverse problem and then solve a PDE at each time step (the so-called offline phase) we would like to point out that these local problems are independent and usually small and therefore the offline phase is embarrassingly parallel, although we did not take advantage of that in our implementation. The global time step (online phase) also consists of a small problem and matrix assembly procedures can be made very efficient by using algebraic tricks, see [14, 23].

We would like to further emphasize the flexibility of the SLMsR. Here we presented an implicit version but explicit time stepping is possible. The method can be transferred to higher dimensions as well as it can be extended to deal with advection-diffusion-reaction problems. Furthermore, the use of inverse problems in the local steps to adjust the basis makes it generally possible to incorporate knowledge coming from measurement data. For this a thorough understanding of the data is necessary (as for any other assimilation method). The SLMsR is promising for practical applications but lots of work needs to be done on the path towards applicability. This includes a numerical analysis which we do not aim at in this work. We would like to explore that opportunity in the future.