---
abstract: 'We consider model order reduction by proper orthogonal decomposition (POD) for parametrized partial differential equations, where the underlying snapshots are computed with adaptive finite elements. We address computational and theoretical issues arising from the fact that the snapshots are members of different finite element spaces. We propose a method to create a POD-Galerkin model without interpolating the snapshots onto their common finite element mesh. The error of the reduced-order solution is not necessarily Galerkin orthogonal to the reduced space created from space-adapted snapshot. We analyze how this influences the error assessment for POD-Galerkin models of linear elliptic boundary value problems. As a numerical example we consider a two-dimensional convection-diffusion equation with a parametrized convective direction. To illustrate the applicability of our techniques to non-linear time-dependent problems, we present a test case of a two-dimensional viscous Burgers equation with parametrized initial data.'
address:
- |
    Graduate School Computational Engineering, Technische Universität Darmstadt,\
    Dolivostr. 15, D-64293 Darmstadt, Germany
- |
    Department of Mathematics, Technische Universität Darmstadt,\
    Dolivostr. 15, D-64293 Darmstadt, Germany
- |
    Graduate School Energy Science and Engineering, Technische Universität Darmstadt,\
    Jovanka-Bontschits-Str. 2, 64287 Darmstadt
author:
- Sebastian Ullmann
- Marko Rotkvic
- Jens Lang
bibliography:
- './adaptive-pod-paper.bib'
title: |
    POD-Galerkin reduced-order modeling with\
    adaptive finite element snapshots
---

proper orthogonal decomposition, adaptive finite elements, model order reduction, reduced basis method

Introduction
============

Model order reduction is a tool to decrease the computational cost for applications where a parametrized PDE problem needs to be solved multiple times for different parameter values. Therefore model order reduction is often studied in the context of optimal control [@GhiglieriUlbrich2014; @NegriEA2013; @TroltzschVolkwein2009] or uncertainty quantification [@ChenQuarteroni2015; @HaasdonkEA2013; @UllmannLang2014]. Snapshot-based model order reduction requires a set of representative samples of the solution, which need to be computed in advance. The solution of the reduced-order model is then represented as a linear combination of these snapshots. The respective coefficients are determined by means of a Galerkin projection, based on a weak form of the governing equations. In this way, the reduced-order model inherits both the spatial structure of typical solutions as well as the underlying physics. Introductions to snapshot-based model order reduction are provided by the textbooks [@HesthavenEA2016; @QuarteroniEA2016].

Standard techniques for model order reduction assume that all snapshots use one and the same spatial mesh. We refer to this case as *static* snapshot computations. In contrast, with *adaptive* snapshot computations we mean that each snapshot may use a different mesh. Combining space-adaptive simulations with model order reduction can be an advantage from two points of view: Firstly, introducing spatial adaptivity to the set-up phase of a reduced-order model can decrease the total computation time if the solution contains local features depending on the parameter. By adapting the mesh to the features at a given parameter value, less degrees of freedom are required to obtain a certain accuracy. Secondly, introducing model order reduction to space-adaptive simulations promises computational speed-up in cases where the solution is susceptible to an approximation in a low-dimensional linear space and needs to be evaluated for multiple different parameter values.

One possible way to implement a POD-Galerkin reduced-order model from snapshots of different discretization spaces is to express the snapshots as elements of some common discretization space. Then one can use standard methods to create a POD basis and compute a respective Galerkin projection of the solution. In the case of strongly varying local refinements, however, a good common discretization space may be relatively high-dimensional, which makes it unattractive from a computational point of view. We show that by expressing the POD basis in terms of the snapshots, it is possible to avoid forming the common discretization space explicitly.

Model order reduction with spatial adaptivity has been studied in [@AliUrbanSteih20XX; @AliUrban2016] for snapshot computations with adaptive wavelets and in [@Yano2016] for snapshot computations with adaptive mixed finite elements. The main issue addressed in these publications is the assessment of the error between the reduced-order solution and the infinite-dimensional true solution. In the case of static snapshot computations, this problem can be circumvented by assuming a sufficiently fine snapshot discretization. Then the error between the reduced-order solution and a corresponding discrete solution can be estimated with the help of the discrete residual. For the case of adaptive snapshot computations, [@AliUrbanSteih20XX; @AliUrban2016] use wavelet techniques to estimate the required dual norm of the continuous residual. In contrast, [@Yano2016] derives a bound for the dual norm of the continuous residual from a special mixed finite element and reduced basis formulation.

The references [@AliUrbanSteih20XX; @AliUrban2016; @Yano2016] focus on model order reduction by the greedy reduced basis method [@PrudhommeEA2002]. In this paper, however, we consider an alternative approach, namely proper orthogonal decomposition (POD) [@HolmesEA1996; @Sirovich1987]. The major difference between both methods lies in the construction of the reduced space used as a test and trial space in a Galerkin procedure. Both methods require a fixed set of training parameters to be chosen in advance, where for time dependent problems, time is viewed as a parameter. From the training parameter set, the greedy reduced basis method selects a set of parameter values in an iterative way using an error estimator and uses the span of the corresponding snapshots as a reduced space. This can save computation time by avoiding the computation of snapshots which have not been selected. The resulting reduced space is close to the one which minimizes the maximum approximation error over the training set. In contrast, POD forms a reduced basis by linearly combining the snapshots corresponding to all training parameter values. The linear combination is done in a way which minimizes the mean square approximation error over the training set, but at the cost of computing all training snapshots. If the dimension of the reduced spaces are increased, both the greedy and the POD space eventually become equal to the span of the snapshots corresponding to all training parameter values.

A POD of snapshots resulting from a static spatial discretization can be implemented in terms of a truncated singular value decomposition [@KunischVolkwein1999]. Therefore, an error estimator is not necessary for creating a POD reduced basis. Because snapshots have to be computed for all training parameters, POD is often applied to time-dependent problems, where snapshot data arise as a by-product of the numerical time stepping scheme. We note that in this context, POD with time-adaptive snapshots has been studied in [@AllaEA2015]. POD with one-dimensional space-adaptive snapshots has already been addressed in [@Lass2014], where the POD computation relies on a polynomial approximation of the snapshots. In contrast, we focus on the two-dimensional case and present a method which does not require an intermediate approximation of the snapshots.

A major difference between greedy and POD reduced basis methods for space adaptive snapshots is caused by the relation between the snapshots and the reduced basis functions: In the greedy reduced basis method, the reduced space is formed by linear combinations of snapshots, while the snapshots are themselves elements of the reduced space. In the POD reduced basis method, the reduced space is also formed by linear combinations of snapshots, but the snapshots are not elements of the reduced space, in general. The difference between the snapshots and their closest approximation in the POD space can be measured in terms of the truncated singular values. This has consequences for the error assessment of POD-Galerkin schemes in presence of space-adapted snapshots. While the main difficulties in the greedy reduced basis method arises from the fact that the error of the reduced-order solution is not necessarily orthogonal to the reduced space anymore, the POD reduced basis method is additionally subject to a truncation error.

This paper is structured as follows: In section \[sec:Pod\], we introduce proper orthogonal decomposition for adaptive finite element snapshots. We propose methods to efficiently compute POD bases for adaptive finite element discretizations with nested refinement. A POD-Galerkin reduced-order model based on adaptive snapshots is formulated in section \[sec:elliptic\] for an elliptic boundary value problem. We prove error statements for the reduced-order solution in presence of adaptive snapshots and compare the results to snapshot computations on a static mesh. The methods and analytic results are illustrated in section \[sec:ConvectionDiffusion\] with a numerical test case involving a linear convection-diffusion equation with parametrized convective direction. The applicability to non-linear time-dependent problems is suggested by the results of section \[sec:Burgers\], which features a Burgers problem with parametrized initial condition.

Proper orthogonal decomposition {#sec:Pod}
===============================

We consider snapshot-based model order reduction, where the solution of a PDE problem is represented in the space spanned by a set of reduced basis functions obtained by linearly combining a set of snapshots. Such reduced-basis functions typically have a global support and contain information about expected spatial structures of the solution.

One method to compute reduced basis functions from snapshots is the proper orthogonal decomposition [@HolmesEA1996; @Sirovich1987]. If the snapshots correspond to coefficient vectors of a finite element discretization on a fixed grid, a POD can be given in terms of the snapshot coefficient matrix [@KunischVolkwein1999]. In the following, we introduce a method that computes a POD of snapshots stemming from an adaptive finite element simulation, where a snapshot matrix can not be created in straight-forward manner.

Method of snapshots {#sec:MethodOfSnapshots}
-------------------

We consider a PDE problem defined over some bounded open spatial domain $\Omega$ and some time and/or parameter domain ${S}$. In particular, we are interested in parametrized elliptic boundary value problems, where ${S}$ is a parameter domain, and in parametrized parabolic initial boundary value problems, where ${S}$ is a tensor product of a time interval and a parameter domain.

Let ${V}$ be the infinite-dimensional Hilbert space used to characterize the solution as a function of space. A typical example is $V=H^1_0(\Omega)$, the Sobolev space of $L^2(\Omega)$ functions with weak first derivatives in $L^2(\Omega)$ and boundary values vanishing in the sense of traces. We denote the $V$-scalar product by $(\cdot,\cdot)_{{V}}$ and the $V$-norm by $\|\cdot\|_{{V}}$.

A proper orthogonal decomposition of snapshots $u_1,\dots,u_N\in V$ can be defined in terms of a system of minimization problems [@KunischVolkwein2001]: Find functions $\phi_1,\dots,\phi_N\in V$ which solve the minimization problems $$\begin{aligned}
\label{eq:PODminimization}
  &\min_{\phi_1,\dots,\phi_R\in V}\sum_{n=1}^N\Big\|u_n-\sum_{k=1}^R(u_n,\phi_k)_V\phi_k\Big\|^2_V,\quad(\phi_i,\phi_j)_V=\delta_{ij},\quad i,j=1,\dots,R\end{aligned}$$ for all $R=1,\dots,N$. The solutions can be computed by an eigenvalue decomposition of the matrix containing the mutual $V$-inner products of $u_1,\dots,u_N$. This approach is often called the *method of snapshots* [@Sirovich1987]. The eigenvalue problem can be written in components as follows: For given ${u_{1}},\dots,{u_{N}}\in{V}$, find $\lambda\in\mathbb R$ and $\vec a=({a_{1}},\dots,{a_{N}})^T\in\mathbb R^{N}$, such that $$\begin{aligned}
  \sum_{j=1}^N(u_i,u_j)_{{V}}{a_{j}}&=\lambda {a_{i}},\qquad i=1,\dots,N.\end{aligned}$$ After defining the snapshot Gramian matrix ${\mathcal G}= (g_{ij})$ with $g_{ij} = (u_i,u_j)_{{V}}$ for $i,j=1,\dots,N$, the matrix form of the set of equations is given by ${\mathcal G}\vec a=\lambda \vec a$. The eigenvalue decomposition of the snapshot Gramian results in eigenvalues $\lambda_1,\dots,\lambda_{N}\in\mathbb R$ and eigenvectors ${\vec a^{1}},\dots,{\vec a^{N}}\in\mathbb R^N$. We order the eigenvalues such that $\lambda_1\geq\dots\geq\lambda_{{D}}>0 = \lambda_{{D}+1}=\dots=\lambda_{N}$ and write the eigenvectors in components as ${\vec a^{r}}=({a^{r}_{1}},\dots,{a^{r}_{N}})^T$ for $r=1,\dots,N$. Then the first POD basis functions are given by linear combinations of snapshots, $$\begin{aligned}
\label{eq:PODbasis}
  \phi_r=\sum_{n=1}^N {u_{n}}\frac{{a^{r}_{n}}}{\sqrt{\lambda_r}},\qquad r=1,\dots,{D}.\end{aligned}$$ The space ${V^{R}}=\operatorname{span}(\phi_1,\dots,\phi_R)$ is called a POD space of dimension $R$ for any $1\leq R \leq D$.

POD reduced-order modeling tries to approximate a solution $u:{S}\rightarrow{V}$ with a function ${u^{R}}:{S}\rightarrow{V^{R}}$ defined by $$\begin{aligned}
\label{eq:PODapproximation}
  {u^{R}} = \sum_{r=1}^R\phi_r{b^{r}},\end{aligned}$$ where $\vec b=({b^{1}},\dots,{b^{R}})^T:{S}\rightarrow\mathbb R^R$ is a POD coefficient vector. Combining and gives a POD approximation in terms of the snapshots, $$\begin{aligned}
\label{eq:PODapproximation2}
  {u^{R}} = \sum_{n=1}^N {u_{n}}\sum_{r=1}^R\frac{{a^{r}_{n}}{b^{r}}}{\sqrt{\lambda_r}}.\end{aligned}$$

One particular choice of POD coefficients is implied by the POD minimization problem and given by a ${V}$-orthogonal projection of $u$ onto ${V^{R}}$: $$\begin{aligned}
  P^Ru := \sum_{r=1}^R\phi_r(\phi_r,u)_V\quad\forall u\in{V},\quad R=1,\dots,D,\end{aligned}$$ which means $b^r=(\phi_r,u)_{V}$ for $r=1,\dots,R$. The orthogonal projection is used as a reference solution later on, because it gives a POD representation with optimal coefficients: $$\begin{aligned}
  \|u-P^Ru\|_{V}&=\inf_{v\in{V^{R}}}\|u-v\|_{V}\quad\forall u\in{V}.\end{aligned}$$ The error of the POD projection of the snapshots can be computed from the POD eigenvalues, $$\begin{aligned}
\label{eq:sumEig}
  \sum_{n=1}^N \|u_n-P^Ru_n\|_V^2 = \sum_{n=R+1}^D \lambda_n,\end{aligned}$$ which implies that the POD projection error of the snapshots decreases monotonically with the POD dimension and that $u_n=P^Du_n$ for $n=1,\dots,N$. Together with this means $$\begin{aligned}
\label{eq:equalSpans}
  \operatorname{span}({u_{1}},\dots,{u_{N}})=\operatorname{span}(\phi_1,\dots,\phi_{{D}}).\end{aligned}$$

Adaptive snapshot spaces
------------------------

In order to compute a set of snapshots, we discretize our PDE problem of interest with adaptive finite elements in space. Let ${V_{1}},\dots,{V_{N}}\subset{V}$ be adapted finite element spaces, so that ${u_{1}}\in{V_{1}},\dots,{u_{N}}\in{V_{N}}$. Let $M_{1},\dots,M_{N}$ be the dimensions of the respective spaces. We focus on $h$-adaptive Lagrangian finite elements with a fixed polynomial degree, so that each snapshot finite element space is defined by a triangulation.

For discretizations on a fixed triangulation, one can represent linear combinations of snapshots by linear combinations of finite element coefficient vectors. In order to do this for adaptive spatial discretizations, however, one must first express the snapshots in terms of a suitable common finite element basis. Therefore, we introduce a space ${V_+}\subset V$ with finite dimension $M_+$, on which we impose two properties:

1.  ${V_+}$ is a finite element space of the same type as ${V_{1}},\dots,{V_{N}}$,

2.  ${V_{1}}+\dots+{V_{N}}\subset{V_+}$ in terms of a vector sum.

A consequence of the first property is that after interpolating all snapshots onto ${V_+}$, we can work with them in the same way as if they were computed on a fixed triangulation. The second property ensures that the error between any snapshot and its representation in ${V_+}$ is zero.

In general, setting ${V_+}={V_{1}}+\dots+{V_{N}}$ would be too restrictive in the sense that it does not necessarily fulfill the first property. Consider, for example, the case where ${V_{1}},\dots,{V_{N}}$ are linear Lagrangian finite element spaces defined over different triangulations of a common spatial domain, like in . While the functions in ${V_{1}}+\dots+{V_{N}}$ are still piecewise linear, they do not always correspond to a finite element discretization on a triangulation. Still, by adding degrees of freedom one can find a triangulation and a respective linear Lagrangian finite element space ${V_+}$ containing ${V_{1}}+\dots+{V_{N}}$.

![Illustration of meshes corresponding to general finite element spaces $V_1$ and $V_2$, their vector sum $V_1+V_2$ and a common finite element space $V_+$ obtained by adding nodes and edges.[]{data-label="fig:SketchFEUnion"}](meshOverlayA.pdf "fig:") ![Illustration of meshes corresponding to general finite element spaces $V_1$ and $V_2$, their vector sum $V_1+V_2$ and a common finite element space $V_+$ obtained by adding nodes and edges.[]{data-label="fig:SketchFEUnion"}](meshOverlayB.pdf "fig:") ![Illustration of meshes corresponding to general finite element spaces $V_1$ and $V_2$, their vector sum $V_1+V_2$ and a common finite element space $V_+$ obtained by adding nodes and edges.[]{data-label="fig:SketchFEUnion"}](meshOverlayAB.pdf "fig:") ![Illustration of meshes corresponding to general finite element spaces $V_1$ and $V_2$, their vector sum $V_1+V_2$ and a common finite element space $V_+$ obtained by adding nodes and edges.[]{data-label="fig:SketchFEUnion"}](meshOverlayC.pdf "fig:")

A more convenient situation is encountered if the snapshots are adapted with the newest vertex bisection algorithm starting from a common initial triangulation. It is known that the smallest common refinement of two such meshes is their overlay [@CasconEA2008; @Stevenson2007], which implies ${V_{1}}+\dots+{V_{N}}={V_+}$. A sketch is given in . Moreover, the mesh of ${V_+}$ can be found by repeated local refinements of any snapshot mesh. Consequently, to interpolate a function from $V_n$ to ${V_+}$ for any $n=1,\dots,N$, one only needs to interpolate the function between successive refinement steps. Because of this favorable property, we focus on refinement by newest vertex bisection in our numerical examples. However, the theory does not depend on this decision.

![Illustration of meshes resulting from refinement by newest vertex refinement based on a common initial triangulation corresponding to a finite element space $V_0$. Dotted lines indicate the next possible refinement. Refining the upper left triangle of the initial mesh results in $V_1$. Refining the upper right triangle of the initial mesh results in $V_2$. The common finite element space $V_+$ equals the overlay of both refined meshes, and therefore $V_+=V_1+V_2$.[]{data-label="fig:SketchFEUnionNVB"}](meshOverlayNVB_A.pdf "fig:") ![Illustration of meshes resulting from refinement by newest vertex refinement based on a common initial triangulation corresponding to a finite element space $V_0$. Dotted lines indicate the next possible refinement. Refining the upper left triangle of the initial mesh results in $V_1$. Refining the upper right triangle of the initial mesh results in $V_2$. The common finite element space $V_+$ equals the overlay of both refined meshes, and therefore $V_+=V_1+V_2$.[]{data-label="fig:SketchFEUnionNVB"}](meshOverlayNVB_B.pdf "fig:") ![Illustration of meshes resulting from refinement by newest vertex refinement based on a common initial triangulation corresponding to a finite element space $V_0$. Dotted lines indicate the next possible refinement. Refining the upper left triangle of the initial mesh results in $V_1$. Refining the upper right triangle of the initial mesh results in $V_2$. The common finite element space $V_+$ equals the overlay of both refined meshes, and therefore $V_+=V_1+V_2$.[]{data-label="fig:SketchFEUnionNVB"}](meshOverlayNVB_C.pdf "fig:") ![Illustration of meshes resulting from refinement by newest vertex refinement based on a common initial triangulation corresponding to a finite element space $V_0$. Dotted lines indicate the next possible refinement. Refining the upper left triangle of the initial mesh results in $V_1$. Refining the upper right triangle of the initial mesh results in $V_2$. The common finite element space $V_+$ equals the overlay of both refined meshes, and therefore $V_+=V_1+V_2$.[]{data-label="fig:SketchFEUnionNVB"}](meshOverlayNVB_BC.pdf "fig:")

Besides a common finite element space of all snapshots, it will be useful to have common finite element spaces of subsets of snapshots available. To this end, we extend our notation in the following way: Let $V_{n_1},\dots,V_{n_K}$ for $K>1$ define a $K$-tuple of snapshot finite element spaces with $1\leq n_k\leq N$ for all $k=1,\dots,K$. We define $V_{n_1\dots n_K}$ as a finite element space of the same type as $V_{n_1},\dots,V_{n_K}$, with $V_{n_1}+\dots+V_{n_K}\subset V_{n_1\dots n_K}$. As a special case we denote ${V_+}= V_{1\dots N}$.

Gramian of adapted snapshots {#sec:Gramian}
----------------------------

The first step in the computation of a POD with the method of snapshots is the creation of the snapshot Gramian, see . For the case of space adapted snapshots, we consider two options: The first option is to represent all snapshots as members of a common finite element space of all snapshots. The second option is to represent pairs of snapshots as members of common finite element spaces of these pairs.

At first we provide an implementation for adaptive finite element snapshots in terms of a common finite element space of *all* snapshots, where we choose ${u_{1}},\dots,{u_{N}}\in{V_+}$. We collect the finite element coefficients of the snapshots with respect to a basis of ${V_+}$ in a set of snapshot coefficient vectors ${\mathcal U_{1}},\dots,{\mathcal U_{N}}\in\mathbb R^{M_+}$ and define a snapshot matrix ${\mathcal U}=({\mathcal U_{1}},\dots,{\mathcal U_{N}})\in\mathbb R^{M_+\times N}$. Let ${\mathcal M}_+$ be the matrix associated with the ${V}$-inner product of functions in ${V_+}$, so that for $u_i,u_j\in{V_+}$ we have $(u_i,u_j)_{{V}}={\mathcal U_{i}}^T{\mathcal M}_+{\mathcal U_{j}}$. Then the snapshot Gramian matrix is given by ${\mathcal G}= {\mathcal U}^T{\mathcal M}_+{\mathcal U}$.

Now we reformulate the computation of the snapshot Gramian so that we only need to create common finite element spaces of *pairs* of snapshots. We consider the computation of a single entry of the snapshot Gramian matrix for a pair consisting of ${u_{i}}\in{V_{i}}$ and ${u_{j}}\in{V_{j}}$. Let ${\mathcal U^{ij}_{i}}$ and ${\mathcal U^{ij}_{j}}$ be the finite element coefficients of ${u_{i}}$ and ${u_{j}}$ with respect to a basis of their common finite element space ${V_{ij}}$ and let ${\mathcal M}_{ij}$ be the matrix associated with the ${V}$-inner product of functions in ${V_{ij}}$, so that $({u_{i}},{u_{j}})_{V}=({\mathcal U^{ij}_{i}})^T{\mathcal M}_{ij}{\mathcal U^{ij}_{j}}$. This means that the Gramian matrix can be filled by creating finite element spaces ${V_{ij}}$ of all pairs of snapshots.

We have presented two ways of creating the snapshot Gramian matrix for the eigenvalue decomposition associated with the POD method of snapshots. In any case, due to the properties of the common finite element spaces we obtain the exact Gramians. The advantage of the first method is that only a single common finite element space has to be created. A possible disadvantage is that the dimension of this space may be very high. In the second method a larger number of lower-dimensional finite element spaces must be created.

POD basis functions and approximation
-------------------------------------

The POD basis functions are determined as linear combinations of snapshots by . If the snapshots are represented as members of ${V_+}$, the POD basis functions are automatically members of ${V_+}$ and can be computed by linearly combining the snapshot finite element coefficient vectors corresponding to a basis of ${V_+}$. If the snapshots are represented as members of their original adapted finite element spaces ${V_{1}},\dots,{V_{N}}$, the POD basis functions can be implicitly defined as linear combinations of snapshots. In this way, forming a basis of ${V_+}$ can be avoided, but applying a linear operator to a single POD basis function means applying this operator to all snapshots. Following this idea, there are multiple ways to represent a POD approximation: in terms of ${V_+}$, or in terms of a linear combination of POD basis functions by , or in terms of snapshots by . By expressing the POD approximation and the POD basis in terms of the snapshots, one can formulate POD Galerkin models based on adaptive snapshots without the need to create the common discretization space.

The theoretical results regarding the $V$-orthogonal POD projection in the last paragraph of have been stated in the context of functions in $V$. Therefore, these results do not depend on whether the snapshots have been computed with a static or an adaptive discretization. Still, the $V$-orthogonal projection requires knowledge of the function to be projected, and is therefore only valuable as a reference.

A different scenario is the computation of POD coefficients by a reduced-order model obtained via Galerkin projection. Here, knowledge of the solution is not necessary to obtain POD coefficients. As we will see in the following section, however, the snapshot discretization influences the accuracy of the POD approximation.

POD Galerkin reduced-order modeling for an elliptic PDE {#sec:elliptic}
=======================================================

We highlight the principal differences between POD Galerkin reduced-order modeling for static and adapted snapshots with an example of a parametrized elliptic boundary value problem. For the case $R=D$ we can use results from the greedy reduced basis theory [@AliUrbanSteih20XX], because according to the POD space equals the span of the snapshots. For $R<D$ we need to take the additional POD truncation error into account. The goal of the error assessment is to understand in which way the POD truncation and the mismatch between the snapshot finite element spaces contribute to the error of the reduced solution.

Weak formulation
----------------

Let ${\mu}\in {S}$ be a parameter vector in a domain ${S}\subset\mathbb R^K$. We define a parametrized bilinear form $a(\cdot,\cdot;{\mu})\colon{V}\times{V}\rightarrow\mathbb R$ which is uniformly coercive with coercivity constant $$\alpha({\mu}) = \inf_{v\in{V}\setminus\{0\}}\frac{a(v,v;{\mu})}{\|v\|_V^2}\geq\alpha>0$$ and uniformly continuous with continuity constant $$\gamma({\mu}) = \sup_{v,w\in{V}\setminus\{0\}}\frac{a(v,w;{\mu})}{\|v\|_V\|w\|_V}\leq\gamma<\infty.$$ We also define a linear form $f(\cdot;{\mu})\colon{V}\rightarrow\mathbb R$ which is uniformly continuous with continuity constant $$\delta({\mu}) = \sup_{v\in{V}\setminus\{0\}}\frac{f(v;{\mu})}{\|v\|_V}\leq\delta<\infty.$$

The parametrized elliptic PDE problem is now formulated as follows: For ${\mu}\in{S}$, find $u({\mu})\in{V}$ such that $$\begin{aligned}
{2}\label{eq:ellipticPDE}
  a(u({\mu}),v;{\mu}) &= f(v;{\mu})&\quad&\forall v\in{V}.\end{aligned}$$ The Lax-Milgram theorem guarantees well-posedness of this problem under the given continuity and coercivity assumptions. Its solution is called *true solution* in the following.

Snapshot computation {#sec:snapshots}
--------------------

To provide snapshots for the POD computation, we introduce a collection of discretized PDE problems associated with a given discrete training set ${S}_N=\{{\mu}_1,\dots,{\mu}_N\}$ with ${\mu}_1,\dots,{\mu}_N\in{S}$. For each $\mu\in{S}_N$ we solve the respective PDE problem with an adaptive discretization scheme, which leads to the snapshot spaces ${V_{1}},\dots,{V_{N}}\subset{V}$. The snapshots are solutions of the following discretized version of : For each $n=1,\dots,N$, find ${u_{n}}\in{V_{n}}$ such that $$\begin{aligned}
{2}
  a({u_{n}},v;{\mu}_n) &= f(v;{\mu}_n)&\quad&\forall v\in{V_{n}}.\end{aligned}$$

For any $n=1,\dots,N$, the subspace property ${V_{n}}\subset{V}$ leads to Galerkin orthogonality between the error $u({\mu}_n)-{u_{n}}$ and the discrete space ${V_{n}}$, which means $$\begin{aligned}
{2}
  a(u({\mu}_n)-{u_{n}},v;{\mu}_n)=0&\quad&\forall v\in{V_{n}}.\end{aligned}$$ We imply a Céa lemma stating $$\begin{aligned}
{2}\label{eq:ceaAdaptive}
  \|u({\mu}_n)-{u_{n}}\|_V\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u({\mu}_n)-P_nu({\mu}_n)\|_V,&\quad& n=1,\dots,N,\end{aligned}$$ where $P_n$ denotes the $V$-orthogonal projection onto the snapshot discretization space ${V_{n}}$, so that $$\|u({\mu}_n)-P_nu({\mu}_n)\|_V = \inf_{v\in{V_{n}}}\|u({\mu}_n)-v\|_V.$$

Reduced-order model {#sec:rom}
-------------------

Using the methods from , we create a POD space ${V^{R}}\subset{V}$ from the snapshots. The respective POD-Galerkin reduced-order model of is formulated as follows: For ${\mu}\in{S}$, find ${u^{R}}({\mu})\in{V^{R}}$ such that $$\begin{aligned}
{2}\label{eq:rom}
  a({u^{R}}({\mu}),v;{\mu}) &= f(v;{\mu})&\quad&\forall v\in{V^{R}}.\end{aligned}$$

The considerations regarding the computation of the snapshot Gramian in give rise to two ways to implementing a reduced-order representation of . In the first approach, the snapshots and reduced basis functions are interpreted as elements of ${V_+}$. Substituting in and testing against the POD basis functions leads to the following implementation: For ${\mu}\in{S}$, find $\vec b({\mu}):{S}\rightarrow \mathbb R^{R}$ such that $$\begin{aligned}
{2}
  \sum_{i=1}^{R}a(\phi_i,\phi_r;{\mu})b^i({\mu}) &= f(\phi_r;{\mu})&\quad&r=1,\dots,R.\end{aligned}$$ This requires building ${V_+}$ and the respective finite element operators, which may be expensive in some cases.

As an alternative, one can substitute in and test against the POD basis functions represented in terms of the snapshots via , which leads to the following implementation: For ${\mu}\in{S}$, find $\vec b({\mu}):{S}\rightarrow \mathbb R^{R}$ such that $$\begin{aligned}
{2}
  \sum_{i=1}^{R}\sum_{m,n=1}^N\frac{{a^{i}_{m}}}{\sqrt{\lambda_i}}a( {u_{m}}, {u_{n}};{\mu})\frac{{a^{r}_{n}}}{\sqrt{\lambda_r}}b^i({\mu}) &= \sum_{n=1}^Nf( {u_{n}};{\mu})\frac{{a^{r}_{n}}}{\sqrt{\lambda_r}}&\quad&r=1,\dots,R.\end{aligned}$$ This requires evaluating the bilinear form for all pairs of snapshots. The cost is similar to creating the snapshot Gramian needed for the POD computation.

In any case, further information about the dependence of the linear and bilinear forms on the parameter ${\mu}$ is needed to obtain a reduced-order model which can be evaluated for any ${\mu}\in{S}$ at a cost which does not depend on the number of spatial degrees of freedom. An example is given in .

Error assessment
----------------

For the error assessment of ${u^{R}}$ from , we first recall the main results for static discretizations. Then we study the adaptive case and point out major differences compared to the static case. We restrict our attention to ${\mu}\in{S}_N$. For practical applications this means that a sufficiently rich snapshot set is assumed, so that the error from discretizing ${S}$ is negligible.

### Static discretization {#sec:errorStatic}

Assume ${V^{R}}\subset{V_{n}}$ for $n=1,\dots,N$. This assumption holds if the snapshots have been computed with static finite elements.

At first we study the error between the snapshots and the reduced-order solution evaluated at the corresponding training parameter values. From ${V^{R}}\subset{V_{n}}$, we can derive a Galerkin orthogonality between this error and the POD space, $$\begin{aligned}
{2}
  a({u_{n}}-{u^{R}}({\mu}_n),v;{\mu}_n)=0&\quad&\forall v\in{V^{R}},\quad n=1,\dots,N.\end{aligned}$$ The Céa lemma, following from coercivity and continuity, states the relation to the POD approximation error, $$\begin{aligned}
{2}
  \|{u_{n}}-{u^{R}}({\mu}_n)\|_V\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u_{n}}-P^R{u_{n}}\|_V,&\quad&n=1,\dots,N.\end{aligned}$$ Moreover, implies ${u^{D}}({\mu}_n) = {u_{n}}$ for $n=1,\dots,N$. The fact that the snapshots are recovered for large enough $R$ is called asymptotic snapshot reproducibility.

For the error with respect to the true solution we use a triangle inequality and respective Céa lemmas to obtain $$\begin{aligned}
{2}
  \|u({\mu}_n)-{u^{R}}({\mu}_n)\|_V&\leq \|u({\mu}_n)-{u_{n}}\|_{V}+\|{u_{n}}-{u^{R}}({\mu}_n)\|_V\\
                              &\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u({\mu}_n)-P_nu({\mu}_n)\|_V+\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u_{n}}-P^R{u_{n}}\|_V\end{aligned}$$ for $n=1,\dots,N$. For maximum efficiency, the errors of the finite element discretization and the POD truncation should be balanced. It should be noticed, however, that the first term stems from the off-line discretization, which is only relevant for the setup of the reduced-order model, while the second term stems from the on-line discretization, which is also relevant for the evaluation time of the reduced-order model.

### Adaptive discretization {#sec:errorAdaptive}

We derive error inequalities similar to the ones in , but for the more general case, where the snapshots are members of different finite element spaces. The assumption ${V^{R}}\subset{V_{n}}$ for $n=1,\dots,N$ is usually not satisfied in the adaptive case, we only have ${V^{R}}\subset{V_+}$ and ${V_{n}}\subset{V_+}$. As a consequence, we are not able to use a Galerkin orthogonality between the reduced-order error ${u_{n}}-{u^{R}}({\mu}_n)$ and the reduced space ${V^{R}}$ for $n=1,\dots,N$.

We start with the error between the solution of the reduced-order model and the true solution. Due to ${V^{R}}\subset{V}$, for any ${\mu}\in{S}$ one can derive a Galerkin orthogonality $$\begin{aligned}
{2}
  a(u({\mu})-{u^{R}}({\mu}),v;{\mu})=0&\quad&\forall v\in{V^{R}}\end{aligned}$$ and a corresponding Céa lemma $$\begin{aligned}
{2}
  \|u({\mu})-{u^{R}}({\mu})\|_V\leq \frac{\gamma({\mu})}{\alpha({\mu})}\|u({\mu})-P^Ru({\mu})\|_V.\end{aligned}$$ We split the right-hand side of the Céa lemma into contributions from the snapshot computation and from the POD truncation. To exclude the error associated with the discretization of the parameter domain, we consider only ${\mu}\in{S}_N$. The derivation starts with adding a zero to the right-hand side of the Céa lemma for the reduced-order model and subsequently uses triangle inequalities and the properties of orthogonal projections, $$\begin{aligned}
{1}
  \|u({\mu}_n)-{u^{R}}({\mu}_n)\|_V&\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u({\mu}_n)-u_n+u_n-P^Ru_n+P^Ru_n-P^Ru({\mu}_n)\|_V\\
  &\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u_n-P^Ru_n\|_V
  +    \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|(I-P^R)(u({\mu}_n)-u_n)\|_V\\
  &\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u_n-P^Ru_n\|_V
  +    \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u({\mu}_n)-u_n\|_V\end{aligned}$$ for $n=1,\dots,N$. Using we obtain $$\begin{aligned}
{1}
  \|u({\mu}_n)-{u^{R}}({\mu}_n)\|_V&\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u_n-P^Ru_n\|_V + \frac{\gamma({\mu}_n)^2}{\alpha({\mu}_n)^2}\|u({\mu}_n)-P_nu({\mu}_n)\|_V\end{aligned}$$ for $n=1,\dots,N$. This means that for parameter values in ${S}_N$, the error between the true and the reduced-order solution can be split into contributions from the projection of the respective snapshot onto the POD space and from the projection of the true solution onto the respective snapshot finite element space. In absence of POD truncation, i.e. for $R=D$, the POD projection error vanishes and we obtain a variant of known results from greedy reduced basis theory [@AliUrbanSteih20XX].

Because in general ${V^{R}}\not\subset{V_{n}}$ for $n=1,\dots,N$, we are not able to derive a Céa lemma for ${u_{n}}-{u^{R}}({\mu}_n)$. A straight-forward approach is using the results from above to obtain $$\begin{aligned}
  &\|{u_{n}}-{u^{R}}({\mu}_n)\|_V\\
  &\qquad\leq\|{u_{n}}-u({\mu}_n)\|_V+\|u({\mu}_n)-{u^{R}}({\mu}_n)\|_V\\  
  &\qquad\leq\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\left(1+\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\right)\|u({\mu}_n)-P_nu({\mu}_n)\|_V+\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|u_n-P^Ru_n\|_V.\end{aligned}$$ Alternatively, we can use another result adapted from the literature [@AliUrbanSteih20XX]: Due to coercivity, continuity, $u_n\in{V^{D}}$ and Galerkin orthogonality between ${V^{D}}$ and the error between the true solution and the solution of the reduced-order model we have $$\begin{aligned}
{1}
  \alpha({\mu}_n)\|{u_{n}}-{u^{D}}({\mu}_n)\|_V^2&\leq a({u_{n}}-{u^{D}}({\mu}_n),{u_{n}}-{u^{D}}({\mu}_n);{\mu}_n)\\
  &= a({u_{n}}-u({\mu}_n),{u_{n}}-{u^{D}}({\mu}_n);{\mu})\\
  &\leq \gamma({\mu}_n)\|{u_{n}}-u({\mu}_n)\|_V\|{u_{n}}-{u^{D}}({\mu}_n)\|_V,\quad n=1,\dots,N,\end{aligned}$$ so that $$\begin{aligned}
{1}
  \|{u_{n}}-{u^{D}}({\mu}_n)\|_V&\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u_{n}}-u({\mu}_n)\|_V,\quad n=1,\dots,N.\end{aligned}$$ Because ${V^{R}}\subset{V^{D}}$ for $R\leq D$ we obtain $$\begin{aligned}
  &\|{u_{n}}-{u^{R}}({\mu}_n)\|_V\\
  &\qquad\leq\|{u_{n}}-{u^{D}}({\mu}_n)\|_V+\|{u^{D}}({\mu}_n)-{u^{R}}({\mu}_n)\|_V\\  
  &\qquad\leq \frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u_{n}}-u({\mu}_n)\|_V +\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u^{D}}({\mu}_n)-P^R{u^{D}}({\mu}_n)\|_V\\
  &\qquad\leq\frac{\gamma({\mu}_n)^2}{\alpha({\mu}_n)^2}\|u({\mu}_n)-P_nu({\mu}_n)\|_V+\frac{\gamma({\mu}_n)}{\alpha({\mu}_n)}\|{u^{D}}({\mu}_n)-P^R{u^{D}}({\mu}_n)\|_V,\end{aligned}$$ with the help of the Céa lemma of the snapshot computation. In any case, the error between a snapshot and the solution of the reduced-order model at the corresponding parameter value contains a component from the finite element computation, which was not present for static snapshots.

Numerical example of a convection-diffusion equation {#sec:ConvectionDiffusion}
====================================================

We apply POD model order reduction to a two-dimensional convection-diffusion problem, where the transport direction serves as a parameter and adaptive finite element snapshots are taken over the parameter interval. The test case illustrates the computational methods introduced in and the theoretical results derived in .

Problem setting {#sec:ConvectionDiffusionProblem}
---------------

We consider a parametrized boundary value problem based on a convection-diffusion equation in two dimensions, $$\begin{aligned}
{2}
  v_x(\mu)\partial_x u + v_y(\mu)\partial_y u - \nu\partial_{xx} u - \nu\partial_{yy} u &= 1 &\quad&\forall (\vec x,{\mu})\in \Omega\times{S},\end{aligned}$$ with solution $u(\vec x,{\mu}):\Omega\times{S}\rightarrow \mathbb R$ for a spatial domain $\Omega=(0,1)\times(0,1)$, a parameter interval ${S}=[0,1]$ and a diffusivity of $\nu=0.01$. The dependence of the convective velocity components $v_x(\mu)$ and $v_y(\mu)$ on the parameter $\mu$ is given by $v_x = \cos(0.25\pi\mu)$ and $v_y = \sin(0.25\pi\mu)$, which means that only the direction of the velocity vector is varied. At the boundary $\partial\Omega$ of the spatial domain $\Omega$, we specify homogeneous Dirichlet conditions, so that $$\begin{aligned}
{2}
  u(\vec x,{\mu}) &= 0 &\quad&\forall (\vec x,{\mu})\in \partial\Omega\times{S}.\end{aligned}$$ A spatial weak form of the parametrized convection-diffusion problem is given by , where $V:=H^1_0(\Omega)$ and the bilinear and linear forms are defined by $$\begin{aligned}
  a(w,\psi;\mu)   &:= \int_\Omega v_x(\mu)\partial_x w\,\psi + v_y(\mu)\partial_y w\,\psi + \nu\partial_x w\,\partial_x \psi + \nu\partial_y w\,\partial_y \psi\; \mathrm d\vec x,\\
  f(\psi) &:= \int_\Omega \psi\; \mathrm d\vec x.\end{aligned}$$

Discretization {#sec:convectionDiffusionDiscretization}
--------------

The general procedure for the snapshot computation and the discretization of the parameter domain has been described in . Here, we focus on the weak form of the convection-diffusion equation. To provide a complete description of the snapshot computation, the choice of ${V_{1}},\dots,{V_{{N}}}$ is characterized in the following.

For the numerical discretization, we use a custom implementation [@Ullmann2016] of piecewise linear Lagrangian finite elements on triangular meshes. As an error indicator, we employ the Matlab function `pdejmps` [@pdetool], defining an error measure for a triangle $K$ by $$E(K) = \left(\frac12\sum_{\tau\in\partial K}h_\tau^2[\vec n_\tau\cdot(\nu\nabla u_h)]^2\right)^\frac12,$$ where $\tau$ denotes an edge along the triangle boundary $\partial K$, $h_\tau$ is the length of the edge, $\vec n_\tau$ is the unit normal, and $[\cdot]$ denotes a jump across the edge. The refinement loop is started from a coarse initial mesh which is identical for all snapshot simulations. In each step, the triangles with the largest error contributions are refined by newest vertex bisection. The marking of triangles and the termination of the refinement loop is done with Matlab `pdeadgsc` [@pdetool].

Numerical solutions for varying parameter values are presented in . Adapted spatial meshes are plotted in together with an overlay of all snapshot meshes constituting ${V_+}$ for further reduced-order modeling.

\

Reduced-order modeling
----------------------

In order to generate snapshots for subsequent reduced-order modeling, we solve the problem for 33 parameter values distributed equidistantly over the parameter interval and store the respective solutions. We compute a POD of these snapshots and create Galerkin reduced-order models of varying dimension, using the same techniques as described in . We can rewrite the model as an equation for the POD coefficient vector: Find $\vec b({\mu}):{S}\rightarrow \mathbb R^{R}$ such that $$\begin{aligned}
{2}
  (A_x^Rv_x(\mu) + A_y^Rv_y(\mu) + A^R_\nu) \vec b &= \vec F^R.\end{aligned}$$ Expressions for the constant model coefficient matrices $A_x^R$, $A_y^R$, $A_\nu^R$ and the right-hand side vector $\vec F^R$ follow from substitution of the POD expansion into and testing against the POD basis functions. In this example, refinement was mostly necessary near the boundary. Therefore, forming the common finite element space of all snapshots was not very costly. Consequently, the snapshot Gramian and the coefficients of the reduced-order model were created using ${V_+}$.

Results {#sec:convectionDiffusionResults}
-------

In order to test how well the reduced-order models reproduce the underlying snapshots, we solve the models for the parameter values in the training set. The error between the snapshots and the solution of the POD-Galerkin model is named $\epsilon_\text{ROM}$. It can be viewed as a benchmark for the Galerkin reduced-order model under the assumption that sufficient snapshot data has been used to create the POD. It is measured in the relative norm induced by the POD, so that $$\label{eq:e_ROM}
  \epsilon_\text{ROM} = \sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}} - {u^{{R}}}({\mu}_{n})\|_V^2}\bigg/\sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}}\|_V^2}.$$

The error between the snapshots and their orthogonal projection on the POD is named $\epsilon_\text{POD}$. It can be viewed as a benchmark for the POD in the sense that it characterizes the ability of the snapshot data to be represented in a low-dimensional space. It is measured in the relative norm induced by the POD, which enables an alternative expression in terms of the POD eigenvalues via : $$\label{eq:e_POD}
  \epsilon_\text{POD} = \sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}} - P^Ru_{{n}}\|_V^2}\bigg/\sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}}\|_V^2} = \sqrt{\sum_{n=R+1}^D\lambda_n}\bigg/\sqrt{\sum_{n=1}^D\lambda_n}.$$

To compare the approximation errors resulting from the model order reduction with the finite element discretization error, we computed reference snapshots $u_1^\text{ref},\dots,u_N^\text{ref}$ with a stricter spatial tolerance. The resulting estimated finite element discretization error of the original snapshots is measured in the same norm as above, so that $$\label{eq:e_FEM}
  \epsilon_\text{FEM} = \sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}}^\text{ref} - u_{{n}}\|_V^2}\bigg/\sqrt{\sum_{{n}=1}^{{N}}\|u_{{n}}^\text{ref}\|_V^2}.$$

The errors defined in – are plotted in (left). We observe that the POD projection error $\epsilon_\text{POD}$ decreases monotonically, which is implied by its representation in terms of the POD eigenvalues. The convergence of the error of the reduced-order model $\epsilon_\text{ROM}$, however, stagnates at some level.

(middle) shows that the stagnation level is related to the finite element error of the snapshots, as suggested by the results of section \[sec:errorAdaptive\]. In the plot, the number of POD basis functions is kept fixed, but the number of finite element functions is varied by changing the tolerance of the refinement algorithm. One can observe that the error of the reduced-order model and the POD projection error roughly follow the finite element discretization error up to some constant.

The error depending on to the number of basis functions is shown in (right) for simulations on a fixed grid given by the overlay of all adapted snapshots. Here both the error of the POD projection and the solution of the reduced-order model converge. This can be explained by the Galerkin orthogonality property detailed in section \[sec:errorStatic\]. The convergence rate is higher in the case of a static finite element discretization. Our interpretation of this behavior is that in the adaptive case the POD basis functions of higher index start approximating spatial artifacts resulting from the different discretizations.

To underline our interpretation of the differences in the POD approximation error $\epsilon_\text{POD}$ between adaptive and static snapshots simulations, we plot a selection of POD basis functions resulting from the adapted snapshots in (top). We observe that the POD basis functions corresponding to the adaptive simulation start exhibiting local variations in the size of typical mesh cells when the index is increased. For comparison, we present a selection of POD basis functions resulting from the static snapshots in (bottom). It is reasonable that the oscillations visible in these plots are necessary to approximate the parameter-dependent physical structures of the solution with increasing accuracy, see .

Despite this qualitative differences in the appearance of the POD basis functions of higher index, we stress that in this example including more than 4 basis functions does not decrease the total error of the reduced-order solution. This is because the POD errors are dominated by the finite element approximation error for $R\geq 4$ in both the adaptive and the static case, see (left).

\

Numerical example of a Burgers equation {#sec:Burgers}
=======================================

To illustrate the potential of POD-Galerkin modeling with adaptive snapshots for non-linear time-dependent problems, we apply our techniques to a Burgers equation. Adaptive finite element snapshots are taken over the time and parameter domain. Note that by choosing the parameter domain equal to a single point, a time-dependent non-parametrized Burgers problem is obtained as a special case. The test case suggests that the theoretical results derived in for a linear elliptic setting can be transferred to a parabolic setting.

Problem setting {#problem-setting}
---------------

We consider a parametrized initial-boundary value problem based on a scalar viscous Burgers equation in two dimensions, $$\begin{aligned}
{2}
  \partial_t u + u\partial_x u - \nu\partial_{xx} u - \nu\partial_{yy} u &= 0 &\quad&\forall (t,\vec x,{\mu})\in I\times \Omega\times{S},\end{aligned}$$ with solution $u(t,\vec x,{\mu}):I\times\Omega\times{S}\rightarrow \mathbb R$ for a time interval $I = (0,1.2]$, a spatial domain $\Omega=(0,1)\times(0,0.5)$, and a parameter interval ${S}= [0,1]$. The problem is a modified version of a numerical example in [@DrohmannEA2012]. A parametrized initial condition is given by $$\begin{aligned}
  u(0,\vec x,{\mu}) - u_0(\vec x){\mu}= 0\quad\forall (\vec x,{\mu})\in\Omega\times{S}\end{aligned}$$ with spatial data $u_0(x):\Omega\rightarrow \mathbb R$ given by a two-dimensional sinusoidal profile $$u_0 = 0.5+0.5\sin\big((x-y-0.75)\pi\big)\sin\big((x+y+0.25)\pi\big).$$ At the boundary we choose a combination of homogeneous Neumann and periodic conditions, so that $$\begin{aligned}
{2}
  u(t,(0,y)^T,{\mu}) - u(t,(1,y)^T,{\mu})&= 0 &\quad&\forall (t,y,{\mu})\in I\times (0,0.5)\times{S},\\
  \partial_y u(t,(x,0)^T,{\mu}) = \partial_y u(t,(x,0.5)^T,{\mu}) &= 0&&\forall(t,x,{\mu})\in I\times (0,1)\times{S}.\end{aligned}$$ We choose a viscosity constant of $\nu=0.001$. Note that for the inviscid Burgers equation with $\nu=0$ we would have $u(t,\vec x,\mu) = \mu u(t\mu,\vec x,1)$.

In the context of , we define $V\subset H^1(\Omega)$ as the Sobolev space of $L^2(\Omega)$ functions which have weak first derivatives in $L^2(\Omega)$ and which satisfy the periodic boundary conditions in a weak sense. A spatial weak form of the parametrized Burgers problem is then given as follows: For $u_0\in L^2(\Omega)$, find $u(t,{\mu}):I\times{S}\rightarrow V$, such that $$\begin{aligned}
{3}
  (u(0,{\mu})-u_0{\mu},\phi)&= 0 &&\forall\phi\in V,&&{\mu}\in{S},\\
  (\partial_t u,\phi) + a(u,\phi) + b(u,u,\phi) &= 0 &\quad&\forall\phi\in V,&\quad&(t,{\mu})\in I\times{S},\end{aligned}$$ where $(\cdot,\cdot)$ denotes the $L^2(\Omega)$-inner product and the bilinear and trilinear forms are defined by $$\begin{aligned}
  a(v,\psi)   &:= \nu\int_\Omega \partial_x v\,\partial_x \psi + \partial_y v\,\partial_y \psi\; \mathrm d\vec x,&
  b(v,w,\psi) &:= \int_\Omega v\partial_x w\,\psi\; \mathrm d\vec x.\end{aligned}$$

Discretization {#discretization}
--------------

We apply Rothe’s method to the problem, i.e. we employ a time discretization followed by space discretizations of the sequence of resulting boundary value problems. We use a constant time step size $\tau = T/({K}-1)$ with according discrete time instances $t_{k}= ({k}-1)\tau$ and semi-discretized solutions $\tilde u_{{k}}\approx u(t_{k},\cdot)$ for ${k}=1,\dots,{K}$. An implicit Euler time discretization is given as follows: For $u_0\in L^2(\Omega)$, find $\tilde u_{1}({\mu}),\dots,\tilde u_{{K}}({\mu}):{S}\rightarrow V$ such that $$\begin{aligned}
{3}
  (\tilde u_{{k}}-u_0{\mu},\phi)&= 0 &\quad&\forall\phi\in V,&\quad&{k}=1,\label{eq:semiTime1}\\
  (\tilde u_{{k}}-\tilde u_{{k}-1},\phi) + \tau a(\tilde u_{{k}},\phi) + \tau b(\tilde u_{{k}},\tilde u_{{k}},\phi) &= 0 &\quad&\forall\phi\in V,&&{k}=2,\dots,{K}.\label{eq:semiTime2}\end{aligned}$$ This semi-discretization is the starting point for further adaptive space discretization and reduced-order modeling.

In order to obtain snapshots for a reduced-order model, we discretize the parameter domain into a training set ${S}_{P}:=\{{\mu}_1,\dots,{\mu}_{{P}}\}$ with ${\mu}_1,\dots,{\mu}_{{P}}\in{S}$. We introduce a discretization in space with adaptive finite elements, so that each snapshot belongs to an individual finite element space chosen with respect to some global error tolerance. We denote the respective fully discrete solutions and finite element spaces by $u_{{k},{p}}\in V_{{k},{p}}$ for ${k}=1,\dots,{K}$ and ${p}=1,\dots,{P}$, so that $u_{{k},{p}}\approx u(t_{k},{\mu}_{p})$. A discretization of the problem with an implicit Euler scheme in time and adaptive finite elements in space for the training parameter values is given as follows: For $u_0\in L^2(\Omega)$ and ${p}=1,\dots,{P}$, find $u_{1,{p}}\in V_{1,{p}},\dots,u_{{K},{p}}\in V_{{K},{p}}$ such that $$\begin{aligned}
{3}
  (u_{k,{p}}-u_0{\mu}_{p},\phi)&= 0 &\;&\forall\phi\in V_{k,{p}},&\;&{k}=1,\\
  (u_{{k},{p}}-u_{{k}-1,{p}},\phi) + \tau a(u_{{k},{p}},\phi) + \tau b(u_{{k},{p}},u_{{k},{p}},\phi) &= 0 &\quad&\forall\phi\in V_{{k},{p}},&&{k}=2,\dots,{K}.\end{aligned}$$ The finite element discretization leads to a set of non-linear algebraic equations which can be solved in parallel for ${p}=1,\dots,{P}$ and in sequence for ${k}=1,\dots,{K}$. The non-linear equations are solved with a standard Newton method using a direct solver for the resulting sequence of linear algebraic equations in each time step.

In the presented example, the time discretization uses a step size of $\tau = 0.005$. The spatial discretization in each time step employs the same error estimation, triangle marking and refinement strategies as in . The adapted meshes and corresponding numerical solutions at various time instances are presented in for a fixed parameter choice of ${\mu}= 1$. To illustrate the parameter dependency of the solution, shows the respective numerical solutions at the initial and final times for ${\mu}= 0.5$.

\
\
\
\

Reduced-order modeling
----------------------

[For consistency, we ]{}denote the snapshots by $u_1\in V_1,\dots,u_{N}\in V_{N}$ with ${N}={K}{P}$. We set $u_{n}= u_{{k},{p}}$ and $V_{N}= V_{{k},{p}}$ with ${N}={k}+({p}-1){K}$ for ${k}=1,\dots,{K}$ and ${p}=1,\dots,{P}$. Using the methods from , we create a POD space $V^{R}\subset V$ with basis functions $\phi_1,\dots,\phi_{R}$. Substituting $V$ by $V^{R}$ in the semi-discretization – gives rise to the following discretized reduced-order model: Find $u^{R}_1({\mu}),\dots,u^{R}_{K}({\mu}):{S}\rightarrow V^{R}$ such that $$\begin{aligned}
{3}
  (u^{R}_{k}-u_0{\mu},\phi)&= 0 &\quad&\forall\phi\in V^{R},&\quad&{k}= 1,\label{eq:rom1}\\
  (u^{R}_{k}-u^{R}_{{k}-1},\phi) + \tau a(u^{R}_{{k}},\phi) + \tau b(u^{R}_{{k}},u^{R}_{{k}},\phi) &= 0 &\quad&\forall\phi\in V^{R},&\quad&{k}=2,\dots,{K}.\label{eq:rom2}\end{aligned}$$ We can rewrite the model as a set of equations for the POD coefficient vectors: Find $\vec b_1({\mu}),\dots,\vec b_{K}({\mu}):{S}\rightarrow \mathbb R^{R}$ such that $$\begin{aligned}
{2}
  M^R\vec b_{k}- \vec b_0{\mu}&= \vec 0, &\quad&{k}=1,\\
  M^R(\vec b_{k}-\vec b_{{k}-1}) + \tau A^R \vec b_{{k}} + \tau B^R(\vec b_{{k}})\vec b_{{k}} &= \vec 0, &\quad&k=2,\dots,{K}.\end{aligned}$$ Expressions for the constant model coefficient matrices $M^R$, $A^R$, $B^R$ and the initial data $\vec b_0$ follow directly from substitution of the POD expansion into – and testing against the POD basis functions. An evaluation of the term involving $B^R$ amounts to a multiplication with a tensor containing ${R}^3$ constant model coefficients.

Results {#results}
-------

We use the same error measures as in to test how well the reduced-order models are able to reproduce the snapshots. The results are shown in (left). We observe that the number of POD basis functions needed for a given relative error is increased in comparison with the convection-diffusion problem presented in the previous section. Otherwise the results are similar. In particular, the POD projection error with respect to the snapshots monotonically decreases, but the convergence of the solution of the reduced-order model stagnates at some point.

(middle) shows the errors depending on the number of finite element degrees of freedom. It can be observed that the level of stagnation is coupled to the finite element discretization error.

The results for a fixed grid, shown in (right), are qualitatively similar to the corresponding results of the convection-diffusion problem, so the same conclusions regarding the use of adaptive snapshots can be drawn in the case of a Burgers problem.

Conclusions
===========

We have extended the framework of POD-Galerkin reduced-order modeling to snapshot data obtained with adaptive finite elements, where each snapshot may be represented in a different finite element space. The considered reduced-order models rely on a representation in a common finite element space of all snapshots. Because creating such a space may be computationally demanding, we have proposed a method to create the reduced-order model without actually building the common finite element space.

The POD projection error of our method converges when the POD dimension is increased. The error between the POD Galerkin solution and the snapshots, however, contains a contribution from the spatial finite element discretization, which does not vanish when the POD dimension is increased. We have shown this effect for the case of linear elliptic boundary value problems.

Our findings are underlined with a numerical example of a convection-diffusion problem. Here we could observe that the error caused by the spatial adaptivity is dominated by the finite element error of the snapshots. Computational results for a Burgers problem suggest that the principal statements can be carried over to a broader class of problems, including non-linear parametrized parabolic PDEs.

A detailed analysis of the influence of the discretization errors onto the reduced solution for non-linear and/or time-dependent problems with space adapted snapshots is still to be done. The ultimate goal of combining adaptivity with model order reduction is to automate the creation of the reduced-order model by adapting the dimension of the reduced basis and the snapshot discretization according to some global error criterion. To reach this goal, it is still necessary to derive adequate error bounds and to combine them in a global adaptation loop.

Acknowledgments
===============

This work was supported by the Excellence Initiative of the German Federal and State Governments via the Graduate School of Excellence Computational Engineering at Technische Universität Darmstadt and the Darmstadt Graduate School of Excellence Energy Science and Engineering. The third author was partly supported by the DFG within the collaborative research center TRR154 “Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks”.
