Approximation of skewed interfaces with tensor-based model reduction procedures: application to the reduced basis hierarchical model reduction approach

Approximation of skewed interfaces with tensor-based model reduction procedures: application to the reduced basis hierarchical model reduction approach

Mario Ohlberger111Institute for Computational and Applied Mathematics, University of Muenster, Einsteinstr. 62, 48149 Münster, Germany,,   and Kathrin Smetana

In this article we introduce a procedure, which allows to recover the potentially very good approximation properties of tensor-based model reduction procedures for the solution of partial differential equations in the presence of interfaces or strong gradients in the solution which are skewed with respect to the coordinate axes. The two key ideas are the location of the interface either by solving a lower-dimensional partial differential equation or by using data functions and the subsequent removal of the interface of the solution by choosing the determined interface as the lifting function of the Dirichlet boundary conditions. We demonstrate in numerical experiments for linear elliptic equations and the reduced basis-hierarchical model reduction approach that the proposed procedure locates the interface well and yields a significantly improved convergence behavior even in the case when we only consider an approximation of the interface.

Keywords: dimensional reduction, tensor-based model reduction, hierarchical model reduction, reduced basis methods, proper generalized decomposition, adaptive modelling

AMS Subject Classification: 65N30,35C20,35J25,35J60

1 Introduction

Fluid flow problems such as subsurface flow or blood flow problems often feature one distinguished (dominant) direction along which the essential dynamics develop. Therefore, tensor-based model reduction procedures such as the proper generalized decomposition (PGD) method, the hierarchical model reduction (HMR), and the reduced basis-hierarchical model reduction (RB-HMR) approach are well suited to compute an efficient and accurate approximation of the full(-dimensional) solution of the underlying partial differential equation (PDE). The common idea of such tensor-based model reduction procedures is to approximate the full solution by a truncated tensor product decomposition of the form , where lie in the computational domain and are associated with different coordinate axes. The resulting model reduction approaches then differ from one another in the way the tensor products , are computed.

In the PGD method, introduced in [1, 20], the tensor products are determined by iteratively solving the Euler-Lagrange equations associated with the considered problem. Alternatively, and in some cases equivalently, they may be computed as the minimizer of the variational functional corresponding to the considered PDE [19, 5]. For an overview on the PGD method we refer to [6, 7].

In contrast, the HMR approach, introduced in [31, 32, 33] and studied in a more general geometric setting in [12, 26], considers a reduced space which is a combination of the full (Finite Element) solution space along the dominant (flow) direction with a reduction space spanned by orthonormal basis functions in the so-called transverse direction. The function then solves a reduced problem obtained by a Galerkin projection onto the reduced space. While in [31, 32, 33, 12, 26] the reduction space is chosen a priori as the span of trigonometric or Legendre polynomials, a highly nonlinear approximation is employed for the construction in the RB-HMR approach [22, 23, 29, 28]. To this end, first a parametrized problem in the transverse direction is derived from the full dimensional problem, where the parameters reflect the influence from the unknown solution in the dominant direction. Then, reduced basis (RB) techniques [24, 27] are applied for the efficient construction of the reduction space from snapshots of the parametrized transverse problem, exploiting their good approximation properties [10, 17]. Thus, both in the construction of the solution manifold of the parametrized lower-dimensional problem and in the subsequent choice of the basis functions, information on the full solution is included to obtain a fast convergence of the reduced solution to the full one. In general, this yields an improved convergence rate compared to a priori chosen reduction spaces [23, 28].

In spite of their mentioned good performance for say fluid flow problems, the approximation capacity of tensor-based model reduction procedures suffers considerably if the target solution exhibits an interface, i.e. a steep gradient or even a discontinuity, which is skewed with respect to the coordinate axes. Such behavior can often be encountered in fluid flow problems and particularly in subsurface flow, where, depending on the permeability of the soil, the saturation profile may form a skewed interface along the water table. This deteriorated convergence behavior is due to the fact that for a full approximation of the skewed interface the saturation or concentration profile in each point in the dominant direction has to be included. In this article we introduce a new ansatz to tackle this problem. We propose to first approximately locate the interface by solving a lower-dimensional model or for simple model problems to infer the location of the interface from data functions. We assume that we have Dirichlet data available at the positions where the interface intersects the boundary of the considered computational domain. Thus we can then infer an approximation of the shape of the interface from the known Dirichlet boundary conditions. Otherwise an approximate shape of the interface can be computed in a preprocessing step. Finally, we prescribe the obtained saturation or concentration profile as the lifting function of the Dirichlet boundary conditions. In this way, we hope to remove the part of the full solution, which causes the bad convergence rate from the approximation process and therefore significantly improve the convergence behavior of the employed tensor-based model reduction approach. This will be demonstrated in numerical experiments.

Alternative to our approach, in [13] an interface or shock propagating in time is included in a time-dependent basis, which is spanned by the eigenfunctions of a linear Schrödinger operator and yields a numerical approximation of a Lax pair. In [18] a reduced basis in space is constructed via a proper orthogonal decomposition of snapshots and the evolution of the coefficients in time is computed by a suitable mapping and thus in an equation-free manner. In the case of parametrized PDEs it is well-known that convection dominated evolution equations where shocks may develop are difficult to tackle with RB methods [27, 24] if linear spaces are employed. The reason for this is that similar to the setting of the skewed interface considered in this article the solution for nearly every time step has to be included in the basis, which deteriorates the approximation properties of the RB space. Therefore, in [21] a nonlinear approximation is applied by employing the method of freezing to decompose the target solution into a shape and group component. Then RB methods are applied to approximate the former while the group component say captures a drift of the interface. Also in [30] the authors propose to employ a nonlinear approximation strategy for the approximation of the solution of parametrized conservation laws in one space dimension. The approach in [30] consists of a partition of the domain induced from a suitable approximation of the shock curve such that the solution in each obtained subdomain is regular. The empirical interpolation method [2] — an interpolation strategy from the RB framework — is used to reconstruct the smooth parts of the solution in the subdomains.

The remainder of this article is organized as follows. In Section 2 we first describe our approach for the location of the interface using the example of subsurface flow and subsequently outline how the location of the interface can be inferred from data functions for linear advection-diffusion problems (Section 2.1). Afterwards, we demonstrate for linear advection-diffusion problems how the information on the location of the interface can be used to remove the interface from the model reduction procedure in Section 2.2. In Section 3 we exemplify this ansatz for the RB-HMR method and present an approach for the derivation of a lower-dimensional parametrized problem particularly suited for the presence of interfaces, which will be validated in Section 4. The capacity of the ansatz proposed in Section 2 to improve the convergence behavior is demonstrated in Section 4 for linear problems for the RB-HMR approach in several numerical experiments, including a test case, where we do not include the exact interface but only an approximation.

2 An ansatz for approximating skewed interfaces with tensor-based model reduction approaches

Let denote the computational domain with Lipschitz boundary , the Dirichlet boundary, and the Neumann boundary. We require that has positive Hausdorff measure. We assume that can be considered as a two-dimensional fiber bundle:

where and denotes the transverse fiber associated with . Note that the generalization to domains with a more complex geometry is straightforward [25]. We define for any the mapping between the fiber associated with and a reference fiber with . Furthermore, we introduce the mapping , defined as for and for . We require that is a -diffeomorphism and that the geometric transformation is differentiable with respect to .

2.1 Locating the interface

In this subsection we propose several approaches in order to (approximately) locate the (skewed) interface. First, we address our motivating example of saturated-unsaturated subsurface flow, where the interface can be located by solving a reduced model. Although the approach in this paper is mainly intended for situations such as subsurface flow, where skewed interfaces naturally occur, we can also think of some simplified settings in which more heuristic approaches can be used to locate the interface. Therefore, we address in a second step a linear advection-diffusion problem. Here, thanks to the nature of the problem, the interface is induced by the data functions.

Saturated-unsaturated subsurface flow

Let be occupied by a homogeneous soil. In the time interval we consider the Richards equation (see e.g. [3, 4]) for the water saturation and the water pressure


Here, denotes the relative permeability, the hydraulic conductivity222Following [3] the hydraulic conductivity describes the ability of the soil to conduct water through it under hydraulic gradients., and the gravity vector. We have normalized the porosity and the viscosity.

Figure 1: Schematic picture of subsurface flow with a skewed water table .

To locate the water table, we first consider the groundwater flow equation for the piezometric head with a free surface. The latter is characterized by an atmospheric pressure and thus describes the location of the water table. Furthermore, denotes the density of the fluid and designates the gravity acceleration. Assuming the incompressibility of water and a flat bottom of the considered domain at level , the piezometric head or potential can be described by the following PDE [3]

Here, accounts for accretion and suitable initial conditions and additional boundary conditions are prescribed. Starting from (LABEL:groundwater) one can derive a dimensionally reduced model for the height of the water table by assuming a hydrostatic pressure distribution [11] or by employing an asymptotic expansion [8]. The PDE for the reduced model then reads


One possible scenario and the corresponding water table are depicted in Fig. 1.

Note that thanks to the structure of (3), solving (3) has the same computational complexity as approximating the solution of (1) in the tensor space

where and . It is therefore reasonable from a computational perspective to solve (3) in a preprocessing step, if the so gained information accelerates the convergence of the tensor-based model reduction procedure.

Employing the location of the water table or a similar interface described by the solution of (3), we can define a function which describes the corresponding saturation or concentration profile. Here, we use the shape of the interface at the Dirichlet boundary as a shape for the whole interface, where the location of the interface is given by and a corresponding indicator function , which is defined as for and for . If we do not have suitable boundary data available, we suggest to compute an approximation of the solution of the PDE on a mesh that is coarse along one direction and fine along the other direction to obtain an approximation of the shape of the interface.

Linear advection-diffusion equation

Similar to the situation above also for linear elliptic and parabolic problems interfaces occurring in the solution can be directly related to interfaces in the data functions. Here, we consider a linear advection-diffusion problem as parabolic problems can be reduced to the former case via discretizing in time. In detail, we consider the following (simplified) model problem for a global pressure


where with for constants and with and . Moreover, and are given Dirichlet and Neumann boundary conditions, respectively.

First, we consider the cases where either the right hand side or the diffusion coefficient exhibit a skewed interface. If we have then the location of the interface in the solution equals the one in the respective data functions. In detail the areas where the shape of changes from nearly flat to a steep slope and vice versa are the areas where exhibits the highest/smallest curvature. As a consequence if the interface is induced by , and is constant or varies only moderately, we expect to be able to identify the “boundary” of the interface by determining where has a maximal (or minimal) derivative in -direction (in -direction). Here, it can be inferred either from the shape of the interface at the Dirichlet boundary or from an approximation of the shape computed in a preprocessing step as outlined in the previous paragraph whether one has to search for the maximal or minimal values of the derivative. In case that the interface is induced by , and is constant or varies only moderately, we can use the smallest and largest derivative in -direction of to deduce (an approximation) of the location of the interface. A possible numerical procedure which determines the location of the interface from and in those cases is discussed below. Then, we can use this “boundary” of the interface either together with the prescribed Dirichlet data or an approximation of the shape to define a function that represents an approximation of the interface in the solution .

In case that the interface in is induced by the right hand side , we expect that for we can still exploit the procedure discussed above for even for rather strong advective fields for the following reason: As we consider steep interfaces we expect that the curvature of dominates the gradient of times the advective field, which is why we expect that we can still obtain a good approximation of the location of the interface in by considering only and neglecting . This will be demonstrated in numerical experiments in §4. In contrast if exhibits an interface we expect that unless the advective field is parallel to this interface the advective term will dominate the term and as a consequence we will observe that has (strong) boundary layers and only a rather moderate slope which does not require additional measures to improve convergence of tensor-based model reduction procedures.

Finally, if we have and an “inflow” boundary conditions on one part of the Dirichlet boundary then for constant advective fields it is possible to infer the location of the interface and thus in from . For one possible example in this context see [9].

We close this subsection with the proposal of a numerical procedure that determines the location of the interface in if this interface is either induced by or . To this end we introduce a fine partition of with elements and a coarse partition of with elements. Here, it is recommended that the partition of has approximately the same number of elements as the mesh that is employed for the FE computations in the model reduction procedure, while, due to computational feasibility, the partition should be significantly coarser as the mesh of the FE approximation used in the model reduction procedure. Then we evaluate or on the partition , determine for each grid point in -direction the elements in -direction with the highest and/or smallest derivative in -direction (or vice versa), and interpolate between the midpoints of those elements to obtain an approximative location of the interface in the solution . Note that by using the coarse partition and an associated Finite Element (FE) space one can compute a (rough) approximation of the solution of the PDE which can then be used to infer an approximation of the shape of the interface.

2.2 Removing the interface from the model reduction procedure

For the sake of clarity we restrict ourselves for the rest of this article to the model problem (4). The ideas in the nonlinear and time-dependent setting are essentially the same. To remove the skewed interface from the model reduction procedure we propose to prescribe the saturation or concentration profile of the preceding subsection as the lifting function of the Dirichlet boundary conditions. In detail we define the solution space such that and consider the following full problem:



The full solution is given as .

Next, we recall that the spaces and introduced in §2.1 satisfy and , supposing compatibility with the boundary conditions prescribed on . We approximate by a linear combination of tensor products


where and , , and define the reduced solution . Depending on the employed tensor-based model reduction method the function solves a reduced problem obtained by a Galerkin projection as for instance in the HMR [26] or the RB-HMR approach (cf. [23] and Section 3). Alternatively in the PGD method (see for instance [19]) the lower dimensional functions and , either minimize a variational functional333 Using the PGD method (cf. [19]), we may obtain by computing sequentially the pair of functions , , as the solution of the minimization problem where shall be the variational functional whose minimizer is the unique solution of (5) for a symmetric bilinear form . or solve the associated Euler-Lagrange equations (see for instance [19, 1, 6]).

Let us assume for a moment that describes the exact location of the interface. Let us furthermore assume that the solution of (5) for , meaning in the case that no interface is present can be approximated exponentially fast by a tensor-based approximation as in (6). Thanks to (5) we hope that we may then also find for any function an approximation (6) which converges with the same or a slightly deteriorated exponential rate in to the solution of (5). In that sense we hope to be able to recover a possibly exponential convergence rate of a tensor-based model reduction procedure in the case of a skewed interface, which will be verified in the numerical experiments in §4. Note that it depends on the applied tensor-based model reduction approach and the underlying problem whether an exponential rate can be realized or not.

For future reference we close this section by introducing some notations. We denote by the symmetric part of and define a -inner product and the induced -norm as and Finally, we define the coercivity and the continuity constants of the bilinear form with respect to the -norm as and

3 Exemplification for the RB-HMR approach

In this section we exemplify the ansatz for the treatment of skewed interfaces proposed in the previous section 2 for the RB-HMR approach introduced in [23, 22]. To obtain a good approximation of solutions exhibiting a skewed interface, we have to eliminate the interface in the solutions of the lower-dimensional problem in the transverse direction. It is therefore crucial to reproduce the balance of the relative terms in the equation of the full problem (5). This is difficult to realize using the approach introduced in [23] as choosing the evaluation of the unknown part of the solution in -direction as a parameter allows too much variation in the scaling of the respective terms to counterbalance them. Thus, we present in §3.2 a new approach for the derivation of a lower dimensional problem based on a FE discretization of the full problem and exemplify it for an advection-diffusion equation in §3.3. We also briefly describe the algorithms for the construction of the reduction space introduced in [23] and comment on necessary adaptations due to the exchange of the parametrized 1D problem. We begin this section by formulating the reduced problem of the RB-HMR approach for the full problem (5).

3.1 Formulation of the reduced problem

We assume orthonormality of the set of functions with respect to the -inner product on and define the reduced space


By using the Galerkin projection we obtain the reduced problem:


which can be rewritten as: Find , such that

To compute an approximation of the coefficient functions , , we introduce a subdivision of with elements of width and maximal step size . We also introduce a corresponding conforming FE space with and basis , . Combining with the reduction space , we define the discrete reduced space

and obtain the discrete reduced problem: Find , , such that


where the discrete reduced solution is defined as for .

3.2 Derivation of a parametrized 1D problem in transverse direction

First, we introduce a subdivision of with elements of width and maximal step size . Furthermore, we introduce an associated conforming FE space with , and basis . Using the FE spaces and , where, denotes the set of polynomials of order in one variable, we may consider the following reference FE approximation of the full problem (5): Find , , such that



and we define .

Next, we introduce for an arbitrary integrand of an integral the quadrature formula


where , are the weights, and , are the quadrature points. Replacing by in the bilinear form and the linear form , we obtain the approximations and . The discrete problem with quadrature then reads: Find , , such that

Figure 2: Plot of the modified basis functions with , , for a polynomial order . The ’o’ marks the removed nodes of the original triangulation . The original basis functions are plotted with a dashed line.

Next, we parametrize (11) by introducing a parameter vector with entries , , in order to find the optimal locations of the quadrature points by applying RB methods and thus to find the optimal points in for solving the lower-dimensional problem in transverse direction. The parameter domain is defined as . As solving (11) is for reasons of efficiency only feasible for small values of , the dimension of is limited to . Therefore, we have in general for functions .

To introduce a coupling between the respective functions we first replace the functions , , by basis functions associated with a new subdivision of . The latter is obtained by deleting all nodes of in the open intervals , , as depicted in Fig. 2. Here and henceforth we assume that the quadrature points are sorted in ascending order and that , where has been defined as the left interval boundary of . denotes the ceil and the floor function. Moreover, we enhance the set of quadrature points by the points , , if . Possible weights of the quadrature , with are defined as

This closes the description of the quadrature rule


Using the quadrature formula (12) instead of (10) we obtain the following coupled system of parametrized 1D partial differential equations in the transverse direction: Given any , find , , such that


for , where and , are the quadrature points of the quadrature formula . Note that (13) is a coupled system of size , where has been defined in the beginning of this subsection as the polynomial order of the FE space . We emphasize that in contrast to [23] we are solving in (13) for the unknown parts of the solution in the dominant direction via the coefficient functions and do not consider them as part of the parameter.

Note that solving (13) without an artificial coupling is equivalent to solving coupled systems of size . This may lead to rather limited variations in the solutions and hence the solution manifold, in which case the solution manifold does not contain all essential information of the full solution in transverse direction. We would hence expect a poor convergence behavior also for which is confirmed by the numerical experiments. In contrast, for and the artificial coupling suggested above we include global information from the dominant direction and hence expect a very much improved convergence behavior, which is again confirmed by the numerical experiments. This stresses the importance of introducing an artificial coupling.

Other choices of an artificial coupling are of course also possible. For instance one could only delete the nodes of in the open intervals , and thus keep the two nodes marked with a circle in Fig. 2. Then one could preserve the FE basis functions associated with the nodes and , and add an additional FE basis function(s), which couples the former. However, due to this additional FE basis function the size of system (13) would increase significantly — at most by . As it is in general only computationally feasible to solve system (13) if it is of small size, this (further) limits the number of quadrature points that can be chosen via RB methods. Since choosing as many quadrature points as possible via RB methods seems preferable as this yields nearly optimal chosen quadrature points, we suggest using the artificial coupling described above as in this case the size of system (13) only depends on . How adding more quadrature points manually without adding more basis functions improves the approximation behavior of the proposed method is subject of future research. Note however that as soon as the quadrature rule is exact no further improvement can be realized without adding new basis functions.

3.3 Example: An advection-diffusion problem

We exemplify the derivation of the coupled system of parametrized 1D partial differential equations for the model problem (5) with non-homogeneous Dirichlet boundary conditions on . For the sake of clarity we restrict our exposition to a rectangular domain , implicating and . The full space thus coincides with and the spaces and coincide with and , respectively.

By applying the quadrature formula defined in (10), we obtain the discrete problem with quadrature: Find , , such that

where the coefficients and are given by

Here we have omitted the on the integrands (cf. (10)) to simplify notations. Using the artificial coupling introduced in the previous subsection and the associated quadrature formula (12) we obtain the parametrized coupled 1D PDE in transverse direction: Given any , find , such that

where for all the coefficients and are given by

3.4 Reduced basis generation — the Adaptive-RB-HMR algorithm

In this subsection we briefly summarize the Adaptive-RB-HMR algorithm introduced in [23] which constructs the reduction space using RB sampling techniques and comment on necessary modifications due to the different parametrized 1D problem.