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.
.5pt \loheadK. Simon, J. Behrens \roheadMsFEM for transient ADEs \rofoot\pagemark \lofoot \cfoot
Multiscale finite elements through advection-induced coordinates for transient advection-diffusion equations
Konrad Simon and Jörn Behrens
University of Hamburg, Department of Mathematics,
Grindelberg 5, 20144 Hamburg, Germany
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 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 [Bensoussan2011, Jikov2012, Pavliotis2008] was originally developed as an analytical tool to study the behavior of solutions to equations with coefficients that rapidly oscillate on fast scales in the limit . 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 [Weinan2003, Weinan2011, Abdulle2012] 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 [Cheng2010, Wang2009] is based on ELLAM [Celia1990, Herrera1993, Russell2002] 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 [Brezzi2000, Hughes1995, Hughes1998]. 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 [Arbogast2000, Arbogast2006], see also [Graham2012]
The MsFEM is closely related to the variational multiscale method and has originally been introduced in [Babuska1983, Babuska1994] and later been substantially refined in [Hou1997, Hou1999]. It relies on the use of non-polynomial 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 [Efendiev2009] 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 [Jablonowski2004, Behrens2006]. 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 [Behrens2006, Graham2012, Efendiev2009, Weinan2003, Bensoussan2011] and to references therein.
In this work we introduce a generalization of the classical MsFEM introduced in [Hou1997] 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 one-dimensional 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.
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.
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.
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 [Efendiev2009, Graham2012]. In contrast to the HMM [Weinan2003] 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
on the unit interval in a periodic setting with rapidly oscillating coefficients , and smooth periodic . In order for (1) to be well-posed we assume that . Then, by the theory of parabolic equations [Evans10] the weak form of problem (1), i.e., find with
has a unique solution that is in and its derivative is in . By compactness we additionally have and hence the initial condition 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 and with and or , where denotes a the magnitude of a resolved numerical scale. By this we mean that oscillations in the data can not be resolved by , 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 in the background velocity the standard case is that . 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 , where is a characteristic length scale, since in our case can be is a very local quantity as well. Instead we are interested in the regime where is large in average, i.e., we suggest to look at with a large (normalized) -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 into a number coarse cells , of mesh width . With each Eulerian node , , of the coarse mesh we associate a basis function that satisfies , 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 we then seek basis functions , , that satisfy
in the weak sense where is the -th -finite element basis function on the coarse scale. Note that these problems have unique weak solutions with and hence . 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 such that
with the initial condition where is the projection onto . The space is given by
The nodal basis functions 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 has to be computed since the other (unique) solution can be computed as , , due to linearity. Consequently the basis functions fulfill , 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
Here we used the Einstein’s sum convention. With this ansatz we design our global weak form as: Find such that
with the initial condition where denotes the projection onto the subspace . Using the method of lines the global finite element space can be formulated as
Note that a basis function is glued together from two parts that are element of the space . 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 -basis functions is possible.
which can be expressed as an ODE:
where , and
and with the initial condition . Note that the second term 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 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 [Efendiev2009].
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 , 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
A little bit more general is the choice
where , respectively , and
The subscale problem for the basis functions (3) needs to be transformed as well: Find a (weak) solution , , such that
A cell then essentially moves with the flow described by the transform without being deformed.
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 . We would nevertheless like to keep possible small scale variations in to be as general as possible. Variations at a scale and larger would have the consequence that the average velocity in each cell 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 denote the characteristics emerging from . The evolution is governed by the set of ODEs that reads
This evolution induces a differentiable transform between -coordinates and -coordinates as long as characteristics do not intersect. The latter is the case if is locally Lipschitz-continuous in and continuous in . 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
Hence, we get
is the determinant of the Jacobian in (21).
Looking again at (19) we see that
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
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 (variational) solution , , such that
In each of the coarse cells we now have to compute two additional quantities, i.e., and . Computation of the former is based on linear interpolation in space between the left and right boundary nodes of which follow the characteristics starting at and . (more details in the paragraph about the implementation). For the latter, i.e., , we use the fact that at each coarse node. Inside the cells we can again interpolate linearly the values and . Note that now at all coarse nodes globally in time. This means that characteristics of the transformed system starting at the coarse nodes are separatrices of the dynamical system and hence no flow across coarse cell boundaries is possible in these coordinates, i.e., the flow is separated by the cells .
The global weak form can be formulated as: Find such that
where is given by (29). The initial condition is where denotes the projection onto the subspace . The global finite element space can be formulated as
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 , , and cells , with mesh size . Now one decides for the transform. For each node either uses (14) or (15) or one solves the ODE (19) on the time interval using a suitable 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 -finite element formulation. For this we create in each cell a uniform mesh with cells , and nodes , , with and . In cell the discretized ODE for the vector of degrees of freedom of the basis function reads
where denotes the nodal values of the standard basis. The system matrices are given by
where denote the basis of standard that consists of the common hat functions. Note again that the computation of and involves terms of the form
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 needs to be computed. By linearity the other (unique) solution is then given by , . 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 -by- (note the periodic boundary conditions):
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 each global basis function is a linear combination of fine scale (standard) basis functions:
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
where is the stored (advection) matrix of the coarse cell 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 be continuous and that 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
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 coarse elements and 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 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 with . We fix these parameters unless we explicitly say something else.
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
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 and for a number of different . 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 -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 .
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
Here the advective term is on average stronger than the diffusive term by a factor of approximately (comparing the mean values) and by a factor of in regions where he diffusivity has a minimum. Oscillations in the background velocity with can be resolved on a grid with only elements whereas oscillations of the diffusivity can not be resolved. Figure 7 shows snapshots of the solution for the case and the case . 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 , 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 we observe a similar to slightly better performance of the MF-MsFEM.