Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling

TL;DR

This paper introduces an adaptive model reduction method for transport-dominated problems that efficiently captures propagating coherent structures by local basis updates and adaptive sampling, achieving significant runtime speedups.

Contribution

It proposes a novel adaptive sampling scheme and local basis update strategy tailored for transport-dominated phenomena, improving efficiency over traditional methods.

Findings

Achieves up to tenfold speedup compared to full models.

Effectively captures wave-type structures with fewer samples.

Demonstrates applicability on rocket engine combustor model.

Abstract

This work presents a model reduction approach for problems with coherent structures that propagate over time such as convection-dominated flows and wave-type phenomena. Traditional model reduction methods have difficulties with these transport-dominated problems because propagating coherent structures typically introduce high-dimensional features that require high-dimensional approximation spaces. The approach proposed in this work exploits the locality in space and time of propagating coherent structures to derive efficient reduced models. Full-model solutions are approximated locally in time via local reduced spaces that are adapted with basis updates during time stepping. The basis updates are derived from querying the full model at a few selected spatial coordinates. A core contribution of this work is an adaptive sampling scheme for selecting at which components to query the full model to compute basis updates. The presented analysis shows that, in probability, the more local the coherent structure is in space, the fewer full-model samples are required to adapt the reduced basis with the proposed adaptive sampling scheme. Numerical results on benchmark examples with interacting wave-type structures and time-varying transport speeds and on a model combustor of a single-element rocket engine demonstrate the wide applicability of the proposed approach and runtime speedups of up to one order of magnitude compared to full models and traditional reduced models.

Equations

$${\mathbf{q}}_k(\mathbf{\mu })=\mathbf{f}({\mathbf{q}}_{k+1}(\mathbf{\mu }),\mathbf{\mu }),k=1,\dots ,K,$$

$$\mathbf{Q}(\mathbf{\mu })=[{\mathbf{q}}_1(\mathbf{\mu }),\dots ,{\mathbf{q}}_K(\mathbf{\mu })]\in {\mathbb{R}}^{N\times K},$$

$$\mathbf{Q}=[\mathbf{Q}({\mathbf{\mu }}_1),\dots ,\mathbf{Q}({\mathbf{\mu }}_M)]\in {\mathbb{R}}^{N\times MK},$$

$$\tilde{\mathbf{f}}(\mathbf{q}(\mathbf{\mu });\mathbf{\mu })=({\mathbf{P}}^T\mathbf{U}{)}^{-1}{\mathbf{P}}^T\mathbf{f}(\mathbf{q}(\mathbf{\mu });\mathbf{\mu }),$$

$${\tilde{\mathbf{q}}}_k(\mathbf{\mu })=\tilde{\mathbf{f}}({\tilde{\mathbf{q}}}_{k+1}(\mathbf{\mu });\mathbf{\mu }),k=1,\dots ,K,$$

$${\partial }_{t}q(x,t)+\mu {\partial }_{x}q(x,t)=0,x\in \Omega ,$$

$$q(x,0)=\frac{1}{\sqrt{\pi 0.02}}\exp \left(-\frac{x^2}{0.0002}\right),$$

$${\mathbf{r}}_k={\mathbf{q}}_k-{\mathbf{U}}_k({\mathbf{P}}_k^T{\mathbf{U}}_k{)}^{-1}{\mathbf{P}}_k^T{\mathbf{q}}_k.$$

$${\mathbf{U}}_{k+1}={\mathbf{U}}_k+{\mathbf{\alpha }}_k{\mathbf{\beta }}_k^T,$$

$$\Vert {\mathbf{S}}_k^T\left(({\mathbf{U}}_k+{\mathbf{\alpha }}_k{\mathbf{\beta }}_k^T){\mathbf{C}}_k-{\mathbf{F}}_k\right){\Vert }_F^2,$$

$${\mathbf{S}}_k^T{\hat{\mathbf{q}}}_k={\mathbf{S}}_k^T\mathbf{f}({\mathbf{U}}_k{\tilde{\mathbf{q}}}_k),{\breve{\mathbf{S}}}_k^T{\hat{\mathbf{q}}}_k={\breve{\mathbf{S}}}_k^T{\mathbf{U}}_k({\mathbf{P}}_k^T{\mathbf{U}}_k{)}^{-1}{\mathbf{P}}_k^T\mathbf{f}({\mathbf{U}}_k{\tilde{\mathbf{q}}}_k),$$

$${\mathbf{F}}_k=[{\hat{\mathbf{q}}}_{k-w+1},\dots ,{\hat{\mathbf{q}}}_k].$$

$$d(\bar{\mathcal{U}},\mathcal{U})=\Vert \bar{\mathbf{U}}-\mathbf{U}{\mathbf{U}}^T\bar{\mathbf{U}}{\Vert }_F^2,$$

$${\Vert {\mathbf{S}}_k^T\left({\mathbf{U}}_{k+1}{\mathbf{C}}_k-{\mathbf{F}}_k\right)\Vert }_F^2=\Vert {\mathbf{S}}_k^T{\mathbf{R}}_k{\Vert }_F^2-\sum _{i=1}^r{\sigma }_i^2,$$

$$d({\mathcal{U}}_{k+1},{\bar{\mathcal{U}}}_{k+1})\le \frac{{\rho }_k^2}{{\sigma }_{\text{min}}^2({\mathbf{F}}_k)}$$

$${\rho }_k^2=\Vert {\breve{\mathbf{S}}}_k^T{\mathbf{R}}_k{\Vert }_F^2+\sum _{i=r+1}^{\bar{r}}{\sigma }_i^2,$$

$${\breve{\mathbf{S}}}_k^T{\mathbf{U}}_{k+1}={\breve{\mathbf{S}}}_k^T{\mathbf{U}}_k.$$

$$\Vert {\mathbf{U}}_{k+1}{\mathbf{C}}_k-{\mathbf{F}}_k{\Vert }_F^2=$$

$$=$$

$$\Vert {\mathbf{F}}_k-{\mathbf{U}}_{k+1}{\mathbf{U}}_{k+1}^T{\mathbf{F}}_k{\Vert }_F^2\le \Vert {\mathbf{U}}_{k+1}{\mathbf{C}}_k-{\mathbf{F}}_k{\Vert }_F^2.$$

$$\Vert {\mathbf{F}}_k-{\mathbf{U}}_{k+1}{\mathbf{U}}_{k+1}^T{\mathbf{F}}_k{\Vert }_F^2\le {\rho }_k^2.$$

$$\Vert {\mathbf{F}}_k-{\mathbf{U}}_{k+1}{\mathbf{U}}_{k+1}^T{\mathbf{F}}_k{\Vert }_F^2=$$

$$\ge$$

$$=$$

$$r_i=\Vert {\mathbf{R}}_k^T{\mathbf{e}}_i{\Vert }_2^2$$

$${\mathbf{S}}_k=[{\mathbf{e}}_{j_1},\dots ,{\mathbf{e}}_{j_m}].$$

$$r_{j_i}\le c_1{\mathrm{e}}^{-c_2{i}^a},$$

$$d({\bar{\mathcal{U}}}_{k+1},{\mathcal{U}}_{k+1})\le c_3{\mathrm{e}}^{-c_2{m}^a},$$

$$m\ge \min \left\{N,{\left(-\frac{1}{c_2}\log \left(\frac{ϵ}{c_3}\right)\right)}^{1/a}\right\}$$

$${\rho }_k^2=\Vert {\breve{\mathbf{S}}}_k^T{\mathbf{R}}_k{\Vert }_F^2,$$

$$\Vert {\breve{\mathbf{S}}}_k^T{\mathbf{R}}_k{\Vert }_F^2=$$

Full text

1 Introduction

Reduced models of large-scale systems of equations are typically constructed in a one-time high-cost offline phase and then used in an online phase to repeatedly compute accurate approximations of the full-model solutions with significantly lower costs. In projection-based model reduction [49, 2, 6], a low-dimensional space is constructed that approximates well the high-dimensional solution space of the full model. Then, the full-model equations are projected onto the low-dimensional space and reduced solutions are computed by solving the projected equations. However, solutions of full models that describe transport-dominated behavior, e.g., convection-dominated flows, wave-type phenomena, shock propagation, typically exhibit high-dimensional features, which means that no low-dimensional space (“linear approximation”) exists in which the full-model solutions can be approximated well; the Kolmogorov n-width is high. Thus, traditional model reduction methods fail for transport-dominated problems [36, 15, 54]. In this work, we exploit that transport-dominated problems typically have a rich structure that is local in nature, which we leverage to derive efficient reduced models.

Model reduction projects the full-model equations on a low-dimensional—reduced–space that is spanned by a set of basis vectors. There are many methods for constructing reduced spaces, including proper orthogonal decomposition (POD) [7, 52], balanced truncation [30, 31], the reduced basis method [44, 56, 50, 25, 49], and interpolatory model reduction methods [24, 2]. In the following, we will be mostly concerned with full-model equations that are nonlinear in the states, where projecting the full-model equations onto a reduced space is typically insufficient to obtain a reduced model that is faster to solve than the full model, because the nonlinear terms entail computations that scale with the dimension of the full-model solution space. One remedy is empirical interpolation, which approximates nonlinear terms by evaluating them at a few, carefully selected interpolation points and approximating all other components via interpolation in a low-dimensional space [5, 23, 17]. We will build on the discrete empirical interpolation method (DEIM) [14, 16], which is the discrete counterpart of empirical interpolation. Note that there are other techniques for nonlinear model reduction, e.g., missing point estimation [3] and the Gauss-Newton with approximated tensors (GNAT) method [13]. All these methods build on the assumption that there is a low-dimensional space in which the full-model solutions can be approximated well, which is violated in case of transport-dominated problems. However, the solutions of transport-dominated problems typically are low dimensional if considered locally in time, which we exploit by approximating the full-model solutions in local low-dimensional spaces that are adapted via low-rank basis updates over time [42, 65]. The basis updates are derived by querying the full model at selected points in the spatial domain. We derive an adaptive sampling scheme that selects where to query the full model to compute the basis updates. An analysis is developed that shows that if the reduced-model residual is local in the spatial domain—which we observe for transport-dominated problems—then only few sampling points are necessary to adapt the local spaces with the proposed adaptive sampling scheme. Even though online operations of the proposed approach incur costs that depend on the dimension of the full-model states, the dimension of the adapted spaces can be selected lower than with static reduced models, which makes the proposed approach computationally efficient and leads to speedups of up to one order of magnitude in our numerical experiments. We summarize the key contributions of this work as follows:

• •

Providing an analysis that relates the locality of the residual in the spatial domain to the number of sampling points required for adapting reduced bases with the proposed approach.

• •

Deriving a sampling scheme to select which components of the full model to sample for deriving basis updates.

• •

Demonstrating the efficiency of the proposed approach on numerical examples with time-varying coefficients and nonlinear dynamics by reporting speedups of up to one order of magnitude compared to full-model solutions and traditional, static reduced models.

For reviewing the literature on model reduction for transport-dominated problems, we broadly distinguish between three lines of research. Literature that we categorize into the first line of research is primarily interested in the stability of reduced models of transport-dominated phenomena. The work [15] formulates the greedy basis construction methods, developed in the reduced basis community [44, 56], in special norms that are better suited for transport-dominated problems to obtain stable approximations. The authors of [55, 62] consider time-space discretizations to guarantee stability of the reduced models. Then, there is a large body of work on closure modeling [58, 19, 37, 38, 26, 48] where methods for stabilizing reduced models have been developed. Into the second line of research, we categorize work that transforms the full model to recover a low-rank structure that can then be exploited with a reduced space. The work on transformations in the model reduction community seems to have originated from [36], where the transport dynamics are separated from the rest of the dynamics via freezing. In [47, 46], transport maps are constructed that reverse the transport and so recover a low-rank structure. Approaches have been developed that aim to find transformations numerically via, e.g., optimization [9]. Shifted POD [45] recovers the shift due to the transport and applies POD after having reversed the shift. The shift operator introduced in [45] is time dependent and separates different transports to be applicable to, e.g., multiple waves traveling at different wave speeds. In [59, 60], snapshots are transformed to better align them before the reduced bases are constructed. The third line of research on model reduction for transport-dominated problems constructs low-dimensional spaces that explicitly target transport-dominated problems. There are approaches that exploit local structure, such as the approach introduced in [54] that builds on the spatial locality of shock fronts. The work [1] constructs bases via $L_{1}$ optimization and uses a dictionary approach. In [12], a special type of adaptation is developed that enriches the reduced space during the online phase if an error indicator signals a high error. In [27], reduced bases are adapted over time via an auxiliary equation that describes the dynamics of the bases. Under certain conditions, the approach can be seen as an approximation of Lax pairs discussed in [22]. Dynamic low-rank approximation [28, 51, 34, 32, 33] adapts bases in time similar to the proposed approach here; however, the bases in dynamic low-rank approximation approaches are constructed so that they are valid for all parameters in a given parameter domain, which is problematic if the solution manifold contains high-dimensional features because the propagating coherent structure depends on parameters, as in, e.g., [54, Example 2.5]. In contrast, the approach developed in this work is local in the parameter domain, additionally to being local in time and space. Furthermore, locality in the parameter domain enables deriving basis updates from only few samples of the full model.

This work is organized as follows. Section 2 sets up the problem and gives the problem formulation. Section 3 demonstrates the local structure in transport-dominated problems and proposes an approach to adaptively sample the full model to construct basis updates. Numerical results in Section 4 on benchmark examples and on a model combustor of a single-element rocket engine demonstrate that the proposed approach achieves significant speedups and is applicable to a wide range of problems. In particular, the numerical results indicate that the proposed approach faithfully approximates interactions between propagating coherent structures traveling at different speeds. Concluding remarks are provided in Section 5.

2 Preliminaries

We briefly discuss model reduction with empirical interpolation in Section 2.1 and then demonstrate on a toy example in Section 2.2 why these traditional model reduction methods fail for problems exhibiting transport-dominated phenomena.

2.1 Model reduction with empirical interpolation

Consider the system of discretized equations that is obtained after discretizing a partial differential equation (PDE) in space and time

[TABLE]

where $\bm{q}_{k}(\bm{\mu})\in\mathbb{R}^{N}$ is the $N$-dimensional state at time $k$ and parameter $\bm{\mu}\in\mathcal{D}$, with parameter domain $\mathcal{D}$. The number of time steps is $K\in\mathbb{N}$. The function $\bm{f}:\mathbb{R}^{N}\times\mathcal{D}\to\mathbb{R}^{N}$ describes the operators of the discretized PDE and typically is nonlinear in the state $\bm{q}_{k+1}(\bm{\mu})$. The time discretization is implicit in time, which means that at each time step $k=1,\dots,K$, a potentially nonlinear, large-scale system of equations has to be solved, e.g., with Newton’s method. The formulation of the full model (1) is different from the formulation in [14] where the linear and nonlinear terms of the equations are treated separately.

To derive a reduced model of the full model (1) with empirical interpolation [5, 14, 16], consider the trajectory at parameter $\bm{\mu}\in\mathcal{D}$

[TABLE]

which is the matrix with the states $\bm{q}_{1}(\bm{\mu}),\dots,\bm{q}_{K}(\bm{\mu})$ as columns. Let the columns of $\bm{U}=[\bm{u}_{1},\dots,\bm{u}_{n}]\in\mathbb{R}^{N\times n}$ be the POD basis of dimension $n\ll N$ obtained from the snapshot matrix

[TABLE]

with parameters $\bm{\mu}_{1},\dots,\bm{\mu}_{M}\in\mathcal{D}$. The space spanned by the columns of $\bm{U}$ is denoted as $\mathcal{U}\subset\mathbb{R}^{N}$ and is a subspace of $\mathbb{R}^{N}$. The critical assumption of traditional model reduction is that the singular values of $\bm{Q}$ decay fast so that only few basis vectors are necessary to approximate well the columns of $\bm{Q}$ in the corresponding space $\mathcal{U}$. Following QDEIM, introduced in [16], select the interpolation points $p_{1},\dots,p_{n}\in\{1,\dots,N\}$ and define the corresponding interpolation points matrix $\bm{P}=[\bm{e}_{p_{1}},\dots,\bm{e}_{p_{n}}]\in\mathbb{R}^{N\times n}$, where $\bm{e}_{p_{i}}\in\mathbb{R}^{N}$ is the $p_{i}$-th canonical unit vector with entry 1 at the $p_{i}$-th component and entry 0 at all other components. Define

[TABLE]

so that $\bm{U}\tilde{\bm{f}}(\bm{q}(\bm{\mu});\bm{\mu})$ is the DEIM approximation of $\bm{f}(\bm{q}(\bm{\mu});\bm{\mu})$ at state $\bm{q}(\bm{\mu})\in\mathbb{R}^{N}$ and parameter $\bm{\mu}\in\mathcal{D}$. Note that computing $\bm{P}^{T}\bm{f}(\bm{q}(\bm{\mu});\bm{\mu})$ typically requires evaluating $\bm{f}(\bm{q}(\bm{\mu});\bm{\mu})$ at the $n$ interpolation points $p_{1},\dots,p_{n}$ only, see [5, 14, 16]. We overload the notation $\tilde{\bm{f}}$ in the following so that if we have a reduced state $\tilde{\bm{q}}(\bm{\mu})\in\mathbb{R}^{n}$, then $\tilde{\bm{f}}(\tilde{\bm{q}}(\bm{\mu});\bm{\mu})=(\bm{P}^{T}\bm{U})^{-1}\bm{P}^{T}\bm{f}(\bm{U}\tilde{\bm{q}}(\bm{\mu});\bm{\mu})$. The reduced model corresponding to $\tilde{\bm{f}}$ is

[TABLE]

with the reduced trajectory $\tilde{\bm{Q}}(\bm{\mu})=[\tilde{\bm{q}}_{1}(\bm{\mu}),\dots,\tilde{\bm{q}}_{K}(\bm{\mu})]\in\mathbb{R}^{n\times K}$. Note that we approximate the state and the nonlinear function in the same space $\mathcal{U}$, which is in contrast to the original use of DEIM in [14] and similar to model reduction via missing point estimation [3]. Once a reduced model (2) is constructed in the offline phase, it is solved in the online phase. The one-time high costs of constructing the reduced model are compensated by approximating the full-model solutions with reduced-model solutions for a large number of parameters online, see, e.g., [49, 2, 6, 43] for details on the wide range of outer-loop and many-query applications where such an offline/online splitting is beneficial.

2.2 Problem formulation

It has been observed that states of problems with transport-dominated behavior can require DEIM (reduced) spaces $\mathcal{U}$ with high dimensions, see, e.g., [36, 15, 54]. The following toy example illustrates this by demonstrating the slow decay of the singular values of the snapshot matrix of solutions of the advection equation. Let $\Omega=(-1,1)\subset\mathbb{R}$ be the spatial domain and consider the advection equation

[TABLE]

with periodic boundary conditions $q(1,t)=q(-1,t)$ and time $t\in[0,\infty)$. Set the initial condition to

[TABLE]

which is the probability density function of a normal random variable with standard deviation 0.01 and mean 0, see Figure 1a. Discretize (3) with a second-order upwind scheme and $N=8192$ inner grid points in the spatial domain and time step size $\delta t=10^{-6}$ and end time $T=0.08$. The singular values of the trajectory $\bm{Q}(\mu)$ are plotted in Figure 1b for $\mu=10$. According to the decay of the singular values, a DEIM space of dimension $n\geq 200$ is necessary to approximate the trajectory $\bm{Q}(\mu)$ with a projection error below $10^{-15}$ in the Euclidean norm, which is a rather slow decay compared to the fast decay observed in many other problems [49, 6]. This example demonstrates that efficient reduced models of transport-dominated problems will have to exploit structure beyond the classical low-rank structure that traditional model reduction methods rely on.

3 Exploiting local structure via online adaptive basis updates

We propose AADEIM (adaptive bases and adaptive sampling) to exploit local structure to construct online adaptive reduced models of full models that exhibit transport-dominated behavior. We focus on two types of local structure: Local low-rankness that we exploit via online adaptive bases and local coherence that enables updating the bases from few samples of the full model. Section 3.1 describes local low-rankness and local coherence in more detail. Section 3.2 discusses the adaptation of the DEIM basis to exploit local low-rank structure and Section 3.3 shows that only few samples are necessary to derive basis updates if the adapted DEIM spaces are locally coherent. Algorithm 1 in Section 3.4 summarizes the proposed approach and provides practical considerations.

In most of this section, we drop the dependency on the parameter $\bm{\mu}$ of the function $\bm{f}$, the state $\bm{q}_{k}$ at time step $k$, and the trajectory $\bm{Q}$, as well as their reduced counterparts. Parametrization of the proposed reduced models based on AADEIM is discussed in Section 3.4.

3.1 Local structure in transport-dominated problems

Consider the toy example introduced in Section 2.2. Let $w\in\mathbb{N}$ be a window size and consider the local trajectory $\bm{Q}_{k}=[\bm{q}_{k-w+1},\dots,\bm{q}_{k}]\in\mathbb{R}^{N\times w}$ at time step $k$, which consists of the $w$ states from time steps $k-w+1$ to time step $k$. Figure 2a compares the singular values of the local trajectory $\bm{Q}_{k}$ to the singular values of the global trajectory $\bm{Q}$ for $w=500$ and $k=1$. The results indicate that the singular values of the local trajectory $\bm{Q}_{k}$ decay orders of magnitude faster than the singular values of the global trajectory $\bm{Q}$. We call this behavior—that the singular values of local trajectories decay fast while the singular values of the global trajectory decay slowly—local low-rank structure. By adapting the DEIM spaces, we will exploit this local low-rank structure in Section 3.2. Similar observations about local low-rank structure are exploited in [28, 51, 34, 32, 33, 22, 39, 41], cf. Section 1. In Appendix A, we analyze analytically an example to demonstrate its low-rank structure.

Let us now consider the local DEIM space $\mathcal{U}_{k}$ of dimension $n=3$ and the corresponding local interpolation points matrix $\bm{P}_{k}$ obtained from the local trajectory $\bm{Q}_{k}$. The residual of approximating the state $\bm{q}_{k}$ at time step $k=500$ with DEIM in $\mathcal{U}_{k}$ is

[TABLE]

Figure 2b shows the decay of the sorted squared components of $\bm{r}_{k}$. The squared residual decays fast in this example, which means that the residual is localized in a few components only. Thus, the residual is low at most of the $N=8192$ components. Intuitively, such a fast decay of the residual means that the basis matrix $\bm{U}_{k}$ of the local DEIM space needs to be corrected at a few components only. We will show that such a fast decay of the residual can be the result of a local coherence structure of the local DEIM spaces, which we will make precise and will exploit with adaptive sampling in Section 3.3.

3.2 Exploiting local low-rank structure: Basis updates

To exploit local low-rank structure as described in Section 3.1, we adapt the DEIM spaces and the DEIM interpolation points during the time steps $k=1,\dots,K$ in the online phase. The adaptation is initialized with the DEIM basis matrix $\bm{U}_{1}$ and interpolation points matrix $\bm{P}_{1}$ at time step $k=1$. Then, at each time step $k=1,\dots,K$, the DEIM basis matrix $\bm{U}_{k}$ is adapted via an additive low-rank update to $\bm{U}_{k+1}$. The update is based on ADEIM [42]. The interpolation points matrix $\bm{P}_{k+1}$ is derived with QDEIM from the adapted basis matrix $\bm{U}_{k+1}$. In the following description, the DEIM interpolant is adapted at each time step $k=1,\dots,K$ for ease of exposition, even though all of the following directly applies to situations where the adaptation is performed at selected time steps only, e.g., every other time step or via a criterion that decides adaptively when to update the DEIM interpolant.

3.2.1 Adaptation with ADEIM

Let $k$ be the current time step and $\bm{U}_{k}$ and $\bm{P}_{k}$ the DEIM basis matrix and the DEIM interpolation points matrix, respectively. Consider the matrix $\bm{F}_{k}\in\mathbb{R}^{N\times w}$ and the coefficient matrix $\bm{C}_{k}=(\bm{P}_{k}^{T}\bm{U}_{k})^{-1}\bm{P}_{k}^{T}\bm{F}_{k}$. The residual of the DEIM approximation of the columns of $\bm{F}_{k}$ is $\bm{R}_{k}=\bm{U}_{k}\bm{C}_{k}-\bm{F}_{k}$. Let now $\bm{S}_{k}=[\bm{e}_{s_{1}^{(k)}},\dots,\bm{e}_{s_{m}^{(k)}}]\in\mathbb{R}^{N\times m}$ be the sampling points matrix corresponding to the $s_{1}^{(k)},\dots,s_{m}^{(k)}\in\{1,\dots,N\}$ samplings points with $m>n$. ADEIM [42] adapts the DEIM basis matrix $\bm{U}_{k}$ to $\bm{U}_{k+1}$ with a low-rank update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}\in\mathbb{R}^{N\times n}$

[TABLE]

with $\bm{\alpha}_{k}\in\mathbb{R}^{N\times r},\bm{\beta}_{k}\in\mathbb{R}^{n\times r}$ with rank $r\in\mathbb{N}$. The ADEIM update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ is the rank-$r$ matrix that minimizes the residual $\bm{U}_{k+1}\bm{C}_{k}-\bm{F}_{k}$ at the sampling points $\bm{S}_{k}$ in the Frobenius norm, which means that the ADEIM update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ minimizes

[TABLE]

see [42] for details on how to compute $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ and the computational costs of computing the update. Minimizing (4) with respect to the update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ keeps the coefficient matrix $\bm{C}_{k}$ fixed, which avoids introducing nonlinear interactions between $\bm{C}_{k+1}$ and $\bm{U}_{k+1}$ and so allows an efficient computation of the update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ as described in [42]. Note, first, that $\|\bm{S}_{k}^{T}(\bm{U}_{k+1}\bm{C}_{k+1}-\bm{F}_{k})\|_{F}^{2}$ with coefficient matrix $\bm{C}_{k+1}$ is upper bounded by (4), which means that by minimizing (4) one minimizes an upper bound of the error that takes the updated coefficient matrix $\bm{C}_{k+1}$ into account. Second, the basis update is performed many times (e.g., every other time step), and so the overall approach alternates between updating the basis by minimizing (4) and updating the coefficient matrix when making a time step, which requires computing $\bm{C}_{k+1}$ corresponding to the updated basis matrix $\bm{U}_{k+1}$; we refer to [42] for details.

Critical for the adaptation is the matrix $\bm{F}_{k}$, because ADEIM adapts the space $\mathcal{U}_{k}$ such that the residual of approximating the columns of $\bm{F}_{k}$ is minimized at the sampling points. A potentially good choice for the columns of $\bm{F}_{k}$ would be the states $\bm{q}_{k-w+1},\dots,\bm{q}_{k}$ of the full model at time steps $k-w+1,\dots,k$; however, the states of the full model are unavailable because their availability would mean the full model has been solved. Instead, we set the columns of $\bm{F}_{k}$ as follows. Let $\breve{\bm{S}}_{k}\in\mathbb{R}^{N\times(N-m)}$ be the complementary sampling points matrix derived from the points $\{1,\dots,N\}\setminus\{s_{1}^{(k)},\dots,s_{m}^{(k)}\}$ that have not been selected as sampling points. Let further $\hat{\bm{q}}_{k}\in\mathbb{R}^{N}$ be the vector with

[TABLE]

which means that the components of $\hat{\bm{q}}_{k}$ corresponding to the sampling points in $\bm{S}_{k}$ are equal to the components of $\bm{f}(\bm{U}_{k}\tilde{\bm{q}}_{k})$ and all other components are approximated via DEIM given by $\bm{U}_{k}$ and $\bm{P}_{k}$. Note that $\tilde{\bm{q}}_{k}$ is the reduced state at time $k$ and that $\bm{f}$ is the function that defines the full model (1). The matrix $\bm{F}_{k}$ that we use in the following for adaptation is

[TABLE]

The vectors $\hat{\bm{q}}_{k-w+1},\dots,\hat{\bm{q}}_{k}$ serve as surrogates of the full-model states $\bm{q}_{k-w},\dots,\bm{q}_{k-1}$ to which the DEIM space is adapted, because the full-model states $\bm{q}_{k+1}$ satisfy $\bm{q}_{k}=\bm{f}(\bm{q}_{k+1})$ for $k=1,\dots,K$, see equation (1).

3.2.2 Analysis of the error of the adapted space

We now provide an analysis of the ADEIM adaptation. Let $\mathcal{U}$ and $\bar{\mathcal{U}}$ be two $n$-dimensional subspaces of $\mathbb{R}^{N}$. We measure the distance between $\mathcal{U}$ and $\bar{\mathcal{U}}$ as

[TABLE]

where $\bm{U}$ and $\bar{\bm{U}}$ are orthonormal basis matrices of $\mathcal{U}$ and $\bar{\mathcal{U}}$, respectively. The distance $d(\bar{\mathcal{U}},\mathcal{U})$ is symmetric and invariant under orthogonal basis transformations, which is true because $d(\bar{\mathcal{U}},\mathcal{U})=n-\|\bm{U}^{T}\bar{\bm{U}}\|_{F}^{2}$ holds. The $n$-dimensional subspace of $\mathbb{R}^{N}$ to which we want to adapt at iteration $k$ is denoted as $\bar{\mathcal{U}}_{k+1}$ and the adapted space is $\mathcal{U}_{k+1}$. The following lemma quantifies the reduction of the residual after an ADEIM update and establishes Proposition 8 that bounds the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space $\mathcal{U}_{k+1}$ with respect to the space $\bar{\mathcal{U}}_{k+1}$.

Lemma 1**.**

Let $\bm{P}^{T}_{k}\bm{F}_{k}$ have rank $n$ and let $\bm{C}_{k}=(\bm{P}^{T}_{k}\bm{U}_{k})^{-1}\bm{P}_{k}^{T}\bm{F}_{k}$ be the coefficient matrix and $\bm{R}_{k}=\bm{U}_{k}\bm{C}_{k}-\bm{F}_{k}$ the residual matrix. Let further $\bm{S}_{k}$ be the sampling points matrix and $\bm{U}_{k+1}=\bm{U}_{k}+\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ the adapted basis matrix of the adapted space $\mathcal{U}_{k+1}$. Let $\bar{r}$ be the rank of $\bm{S}_{k}^{T}\bm{R}_{k}$ and let $r\in\mathbb{N}$ with $r\leq\bar{r}$ be the rank of the update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$. Then,

[TABLE]

where $\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{\bar{r}}>0$ are the singular values of $\bm{S}_{k}^{T}\bm{R}_{k}$.

Proof.

Because $\bm{P}^{T}_{k}\bm{F}_{k}$ has rank $n$, and because $\bm{P}^{T}_{k}\bm{U}_{k}$ is regular given the interpolation points are selected by standard empirical interpolation algorithms [5, 14, 16], the coefficient matrix $\bm{C}_{k}$ has full row rank. Then, the statement of this lemma follows from [42, Lemma 3.5] by transforming the generalized symmetric positive definite eigenproblem into the symmetric eigenproblem with matrix $(\bm{S}_{k}^{T}\bm{R}_{k})^{T}(\bm{S}_{k}^{T}\bm{R}_{k})$, which then is equivalent to computing the singular value decomposition of $\bm{S}_{k}^{T}\bm{R}_{k}$ so that the squared singular values of $\bm{S}_{k}^{T}\bm{R}_{k}$ are equal to the eigenvalues of the generalized eigenproblem. This shows this lemma by using the last identity of [42, proof of Lemma 3.5]. ∎

Proposition 1**.**

Consider the same setup as in Lemma 1 and additionally assume that $\bm{F}_{k}$ has rank $n$ and its columns are in the $n$-dimensional space $\bar{\mathcal{U}}_{k+1}$. Let further $\mathcal{U}_{k+1}$ be the adapted space derived with the rank-$r$ ADEIM update and sampling points matrix $\bm{S}_{k}$. Then, it holds

[TABLE]

with

[TABLE]

where $\breve{\bm{S}}_{k}$ denotes the complementary sampling matrix of $\bm{S}_{k}$ and $\sigma_{\text{min}}(\bm{F}_{k})$ is the smallest non-zero singular value of $\bm{F}_{k}$.

Proof.

First, note that the matrix $\breve{\bm{S}}_{k}^{T}\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ has zero entries only, because the ADEIM update $\bm{\alpha}_{k}\bm{\beta}_{k}^{T}$ updates only rows of $\bm{U}_{k}$ corresponding to the sampling points $\bm{S}_{k}$, see [42, Lemma 3.5]. Thus, it holds

[TABLE]

Consider now the norm of the residual with respect to the adapted space $\mathcal{U}_{k+1}$

[TABLE]

Equation (10) follows because $\breve{\bm{S}}_{k}$ is the complementary sampling points matrix of $\bm{S}_{k}$. Equation (11) follows from Lemma 1. Equation (12) holds because of (9). Equation (13) uses $\bm{R}_{k}=\bm{U}_{k}\bm{C}_{k}-\bm{F}_{k}$ and that $\sum_{i=1}^{\bar{r}}\sigma_{i}^{2}=\|\bm{S}_{k}^{T}\bm{R}_{k}\|_{F}^{2}$. Consider now the projection of the columns of $\bm{F}_{k}$ onto $\mathcal{U}_{k+1}$ and observe that

[TABLE]

Note that the basis matrix $\bm{U}_{k+1}$ obtained with the ADEIM update is not necessarily orthonormal and therefore it might be necessary to consider the oblique projection $\bm{U}_{k+1}(\bm{U}_{k+1}^{T}\bm{U}_{k+1})^{-1}\bm{U}_{k+1}^{T}\bm{F}_{k}$ onto the column span $\mathcal{U}_{k+1}$ of $\bm{U}_{k+1}$ in (14); however, the projection error in the Frobenius norm is the same for both projections. With (13) and the definition of $\rho_{k}^{2}$ follows

[TABLE]

Let now $\bar{\bm{U}}_{k+1}$ be an orthonormal basis matrix of $\bar{\mathcal{U}}_{k+1}$. Since $\bar{\mathcal{U}}_{k+1}$ is spanned by the columns of $\bm{F}_{k}$, there exists a full-rank matrix $\tilde{\bm{F}}_{k}\in\mathbb{R}^{n\times w}$ such that $\bm{F}_{k}=\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$. We obtain

[TABLE]

where we used $\sigma_{\text{min}}(\tilde{\bm{F}}_{k})=\sigma_{\text{min}}(\bm{F}_{k})$, which holds because $\bar{\bm{U}}_{k+1}$ is orthonormal. The matrix $\bm{F}_{k}$ has rank $n$ and therefore there are $n$ non-zero singular values corresponding to the basis vectors in $\bar{\bm{U}}_{k+1}$. Combining (16) with (15) and the definition of $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ shows (7). Note that if $\bm{U}_{k+1}$ is not orthonormal, then the oblique projection is used in (16) to obtain $\|\bar{\bm{U}}_{k+1}-\bm{U}_{k+1}(\bm{U}_{k+1}^{T}\bm{U}_{k+1})^{-1}\bm{U}_{k+1}^{T}\bar{\bm{U}}_{k+1}\|_{F}^{2}\sigma_{\text{min}}^{2}(\bm{F}_{k})\leq\rho_{k}^{2}$, which leads to (7) as well because $\|\bar{\bm{U}}_{k+1}-\bm{U}_{k+1}(\bm{U}_{k+1}^{T}\bm{U}_{k+1})^{-1}\bm{U}_{k+1}^{T}\bar{\bm{U}}_{k+1}\|_{F}^{2}=\|\bar{\bm{U}}_{k+1}-\hat{\bm{U}}_{k+1}\hat{\bm{U}}_{k+1}^{T}\bar{\bm{U}}_{k+1}\|_{F}^{2}=d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ for an orthonormal basis matrix $\hat{\bm{U}}_{k+1}$ that spans the same column space $\mathcal{U}_{k+1}$ as $\bm{U}_{k+1}$.

∎

3.3 Exploiting local coherence: Adaptive sampling

Proposition 8 shows that the choice of the sampling points $\bm{S}_{k}$ influences the bound of the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space $\mathcal{U}_{k+1}$. In this section, we derive an adaptive sampling strategy that minimizes the upper bound derived in Proposition 8 in case of full-rank ADEIM updates. Then, we show that if the adaptive sampling strategy is used, the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ decays at least as fast with the number of sampling points $m$ as the norm of the DEIM residual, which in turn means that only few sampling points are necessary to adapt the basis if the residual decays fast.

3.3.1 Adaptive sampling strategy based on residual

Following Proposition 8 and the decay factor $\rho_{k}$ defined in (8), we select the sampling points so that $\|\breve{\bm{S}}_{k}^{T}\bm{R}_{k}\|_{F}^{2}$ is minimized. Note that $\sigma_{\text{min}}(\bm{F}_{k})$ in (7) is independent of the sampling points $\bm{S}_{k}$ and therefore it is sufficient to derive sampling points that lead to a small decay factor $\rho_{k}$.

Consider the component-wise residual

[TABLE]

and let $j_{1},\dots,j_{N}$ be an ordering of $1,\dots,N$ such that $r_{j_{1}}\geq r_{j_{2}}\geq\dots\geq r_{j_{N}}$. Then, select the first $j_{1},\dots,j_{m}$ components as sampling points and form the corresponding sampling points matrix

[TABLE]

If a full-rank ADEIM update is applied, i.e., $\bar{r}=r$ in Proposition 8, then this choice of sampling points is optimal in the sense that the bound $\rho_{k}$ is minimized.

We now show that a fast decay in the residual implies a fast decay in the error bound of $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$.

Proposition 2**.**

Consider the same setup as in Proposition 8. Let $j_{1},\dots,j_{N}$ be an ordering of $\{1,\dots,N\}$ such that the component-wise residual defined in (17) decays as

[TABLE]

with rate $a>1$ and constants $c_{1},c_{2}>0$. The error of the adapted space $\mathcal{U}_{k+1}$ is bounded as

[TABLE]

if a full-rank ADEIM update is applied and $m$ sampling points are selected via the adaptive sampling (18). The constant $c_{3}=c_{1}/((1-\mathrm{e}^{-c_{2}})\sigma_{\text{min}}^{2}(\bm{F}_{k}))$ is independent of $m$. In particular, setting the number of sampling points to

[TABLE]

guarantees $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})\leq\epsilon$ for a threshold $\epsilon>0$.

Proof.

In case of a full-rank update, the factor $\rho_{k}^{2}$ in (7) is

[TABLE]

as shown in Proposition 8. Then, we obtain with (19) and the adaptive sampling (18)

[TABLE]

where the last inequality holds because $i\geq 0,a>1$ and thus $(i+(m+1))^{a}\geq i^{a}+(m+1)^{a}$. Using that $a>1$, we obtain $i^{a}\geq i$ and $(m+1)^{a}\geq m^{a}$ and thus

[TABLE]

Set $c_{3}=c_{1}/((1-\mathrm{e}^{-c_{2}})\sigma^{2}_{\text{min}}(\bm{F}_{k}))$ to obtain (20) with Proposition 8. Setting $m$ as in (21) shows $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})\leq\epsilon$. ∎

Proposition 21 is formulated for an exponential decay of the residual (19); however, the steps in the proof of Proposition 19 seem to generalize to other decay behavior of the residual such as an algebraic decay. The key in the proof of Proposition 19 is (22), which builds an argument based on a geometric series that the sum of the norms of the residual decays exponentially in the number of sampling points $m$. In case of, e.g., an algebraic decay of the sum of norms, we expect that such an upper bound would decay algebraically. Then, the rate of the residual decay can be carried over to the decay of the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space $\mathcal{U}_{k+1}$ with mild modifications of the proof of Proposition 21.

3.3.2 Local coherence

Following, e.g., [10], define the coherence of a space $\mathcal{U}$ as

[TABLE]

where $\bm{U}$ is an orthonormal basis of $\mathcal{U}$ and $\bm{e}_{i}$ is the $i$-th canonical unit vector. Define further the local coherence as

[TABLE]

for $i=1,\dots,N$, see, e.g., [29, 63]. We refer to [11, 10, 29, 63] for details.

We now show that the rate of the decay of the local coherence of the current DEIM space $\mathcal{U}_{k}$ and of the space $\bar{\mathcal{U}}_{k+1}$ to which we want to adapt is reflected in the decay of the bound of the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space $\mathcal{U}_{k+1}$ with respect to the number of sampling points $m$. Thus, we now relate the decay of the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space to a property of the spaces $\mathcal{U}_{k}$ and $\bar{\mathcal{U}}_{k+1}$, namely their decay of the local coherence.

Lemma 2**.**

Let $j_{1},\dots,j_{N}$ be an ordering of $\{1,\dots,N\}$ such that

[TABLE]

holds for $i=1,\dots,N$, with $c_{4},\bar{c}_{4},c_{5},\bar{c}_{5}>0$ and $a,\bar{a}>1$. Let the columns of $\bm{F}_{k}$ be in $\bar{\mathcal{U}}_{k+1}$. Define $\bm{C}_{k}=(\bm{P}^{T}_{k}\bm{U}_{k})^{-1}\bm{P}_{k}^{T}\bm{F}_{k}$, then for the residual $\bm{R}_{k}=\bm{U}_{k}\bm{C}_{k}-\bm{F}_{k}$ holds

[TABLE]

with a constant $c_{7}(\Lambda_{k})>0$ that is independent of $a,\bar{a},c_{5},\bar{c}_{5}$ and that depends on $\Lambda_{k}=\|(\bm{P}^{T}_{k}\bm{U}_{k})^{-1}\|_{2}$.

Proof.

Denote with $\bm{U}_{k}^{(i)}$, $\bar{\bm{U}}_{k+1}^{(i)}$, and $\bm{F}_{k}^{(i)}$ the $i$-th row of $\bm{U}$, $\bar{\bm{U}}_{k+1}$, and $\bm{F}_{k}$, respectively. The matrix $\bar{\bm{U}}_{k+1}$ is an orthonormal basis matrix of $\bar{\mathcal{U}}_{k+1}$. Let further $\tilde{\bm{F}}_{k}\in\mathbb{R}^{n\times w}$ be a matrix such that $\bm{F}_{k}=\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$. Note that $\|\bm{F}_{k}\|_{2}=\|\tilde{\bm{F}}_{k}\|_{2}$, because $\bar{\bm{U}}_{k+1}$ is orthonormal. Then, with $\bm{C}_{k}=(\bm{P}_{k}^{T}\bm{U}_{k})^{-1}\bm{P}_{k}^{T}\bm{F}_{k}$, follows

[TABLE]

Now make the following estimate

[TABLE]

because $\|\bm{F}_{k}\|_{2}=\|\tilde{\bm{F}}_{k}\|_{2}$. Further, we have

[TABLE]

Set now $\bar{c}_{6}=\bar{c}_{4}$ and $c_{6}=c_{4}\|(\bm{P}_{k}^{T}\bm{U}_{k})^{-1}\|_{2}^{2}$ to bound

[TABLE]

Now set $c_{7}(\Lambda_{k})=c_{6}+\bar{c}_{6}+2\sqrt{c_{6}\bar{c}_{6}}$ to obtain

[TABLE]

which shows the proposition. ∎

Since we consider subspaces of finite-dimensional spaces only, i.e., $N$ is finite, the bounds in (24) hold for any subspace by increasing the constants and choosing the rate close to 1. However, Lemma 2 still is meaningful because it shows that the constants and rates that appear in (24) are obtained in the bound of the residual in (25) as well. Thus, the local coherence structure is directly reflected in the decay of the DEIM residual. The constant $c_{7}(\Lambda_{k})$ depends on $\Lambda_{k}=\|(\bm{P}_{k}^{T}\bm{U}_{k})^{-1}\|_{2}$, which is a well-studied quantity in the context of empirical interpolation and typically found to be low in practice [5, 14, 16, 40].

The following proposition combines the local coherence structure exploited in Lemma 2 and the adaptive sampling of Proposition 21 to derive the number of sampling points $m$ that are required to achieve $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})\leq\epsilon$ in probability. We draw random samples from the space $\bar{\mathcal{U}}_{k+1}$ to avoid making assumptions on the specific right-hand side matrix $\bm{F}_{k}$, which still allows us to bound in probability the condition number of $\bm{F}_{k}$ that is required when applying Proposition 21 in the following. We draw a random sample from the space $\bar{\mathcal{U}}_{k+1}$ by taking an i.i.d. standard Gaussian matrix $\tilde{\bm{F}}_{k}$ and multiplying it with the basis matrix $\bar{\bm{U}}_{k+1}$. Other distributions over the elements of $\bar{\mathcal{U}}_{k+1}$ can be considered to derive similar statements as in the proposition below.

Proposition 3**.**

Assume that Lemma 2 applies and consider the same setting as in Proposition 8 except that $\bm{F}_{k}=\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$ with $\tilde{\bm{F}}_{k}$ being an $n\times n$ matrix with independent and identically distributed (i.i.d.) standard Gaussian entries and $\bm{P}_{k}^{T}\bar{\bm{U}}_{k+1}$ having full rank $n$. Then, $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})\leq\epsilon$ with probability at least $1-\delta$ if a full-rank update is applied and

[TABLE]

where $c_{8}(\Lambda_{k})$ is a constant independent of $m,a$, and $\bar{a}$.

Proof.

Since Lemma 2 applies, we have with (23) that

[TABLE]

With Proposition 8, which is applicable with a full-rank update because $\bm{P}^{T}_{k}\bm{F}_{k}=\bm{P}^{T}\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$ has full row rank $n$ if $\tilde{\bm{F}}_{k}$ is regular, follows

[TABLE]

where $\kappa(\bm{F}_{k})$ is the spectral condition number of $\bm{F}_{k}$. Since $\bm{F}_{k}=\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$ with $\bar{\bm{U}}_{k+1}$ orthonormal, we have $\kappa(\bm{F}_{k})=\kappa(\tilde{\bm{F}}_{k})$. We now bound the spectral condition number of $\tilde{\bm{F}}_{k}$ with high probability by exploiting that $\tilde{\bm{F}}_{k}$ has i.i.d. standard Gaussian entries. The work [18, Theorem 1.1] and [4, Theorem 1.1] show that the spectral condition number of $\tilde{\bm{F}}_{k}$ is bounded in probability as

[TABLE]

for $\delta>0$ and with a positive constant $c$ that depends $n$. Thus, we have $P[\kappa(\tilde{\bm{F}}_{k})\leq c\delta^{-1}]\geq 1-\delta$ and $P[\kappa^{2}(\tilde{\bm{F}}_{k})\leq c^{2}\delta^{-2}]\geq 1-\delta$, which leads to

[TABLE]

with probability at least $1-\delta$ and $c_{8}(\Lambda_{k})=(c_{7}(\Lambda_{k})c^{2})/(1-\mathrm{e}^{-\min\{c_{5},\bar{c}_{5}\}})$. Setting $m$ as in (26) leads to $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})\leq\epsilon$ with probability $1-\delta$. ∎

In the following, we focus on problems with moving coherent structures that are local in the spatial domain. Then, typically, the residual (17) decays fast, cf. Section 2.2, and local basis updates with AADEIM are sufficient. In contrast, if the basis updates have to take into account new features that are global in the spatial domain, then the residual (17) might decay slowly. One can then either perform a global update in the sense that the number of sampling points is $m=N$ or perform multiple local updates. Since the adaptive sampling strategy selects the points with the largest residual, multiple local updates can have a similar effect as one global update.

3.4 Practical considerations and algorithm

We now provide details on the practical implementation of the AADEIM approach and summarize AADEIM in Algorithm 1.

3.4.1 Practical considerations

The AADEIM model is initialized with a DEIM interpolant with basis matrix $\bm{U}_{1}$ and $\bm{P}_{1}$. We propose to construct $\bm{U}_{1}$ and $\bm{P}_{1}$ from a local trajectory computed with the full model. This means that the full model is solved for $w_{\text{init}}\in\mathbb{N}$ time steps, with $w_{\text{init}}\ll K$, and the initial DEIM interpolant $(\bm{U}_{1},\bm{P}_{1})$ is constructed from the corresponding full-model states. Initializing the AADEIM model with a DEIM interpolant obtained from full-model states has two benefits. First, no offline phase is necessary to initialize the AADEIM model, which means that it is unnecessary to develop (e.g., greedy) strategies to build an initial AADEIM model offline. Second, because there is no offline phase, the AADEIM model is initialized and then adapted for the parameter at hand. Thus, an explicit parametrization of the AADEIM model is unnecessary, which avoids the common challenges of parametrized model reduction [49, 6] and sampling—potentially high-dimensional—parameter spaces. However, having no offline phase also means that information available from before the online phase is ignored. For example, if only certain regions in the parameter domain lead to transport-dominated phenomena, then a reduced basis constructed in the traditional way with, e.g., greedy methods [44, 56, 49], could be used without adaptation for parameters corresponding to non-transport-dominated phenomena. Furthermore, one may argue that the initialization of the proposed approach with a few full-model states can be seen as an offline phase. The number of full-model states that are used for initialization is then an offline parameter that has to be chosen adequately. In all our experiments, the number of initial full-model states is set to a fixed value, which demonstrates that reduced models based on AADEIM seem to be robust with respect to the initial window size. A more detailed discuss is given in Section 4.2.2 and Figure 10.

The adaptive sampling strategy as described in Section 3.3.1 requires the residual $\bm{R}_{k}$ at all $N$ components. To avoid computing the residual $\bm{R}_{k}$ at all components at each adaptation iteration, we adapt the sampling points at every $z$-th iteration instead. Thus, only at every $z$-th iteration, the residual $\bm{R}_{k}$ is computed at all components to adapt the sampling points, whereas at all other iterations, the residual is computed only at the sampling points. Note that even if the sampling points are adapted at each iteration $k=1,\dots,K$, and so the residual is computed at each component, it still leads to a lower runtime to compute the ADEIM update only at $m<N$ sampling points, because the costs of computing the ADEIM update scales linearly in the number of sampling points $m$ [42].

One parameter of the proposed approach is the window size $w$. The analysis developed in Section 3.2.2 gives guidance on how to select the window size. Lemma 1 and Proposition 8 require that the columns of $\bm{F}_{k}\in\mathbb{R}^{N\times w}$ and $\bm{P}^{T}\bm{F}_{k}\in\mathbb{R}^{N\times w}$ span $n$-dimensional spaces. Thus, the condition $w\geq n$ is imposed on the window size $w$. At the same time, the discussion in Section 3.1 indicates that a smaller window size $w$ leads to faster decay of the singular values and thus to a potentially lower error of the adaptive reduced model. Given these two insights, we set $w$ at least $n$ and not much larger than $n$ in the following. In the numerical examples in Section 4.2 and Section 4.3, the window size is set to $w=n+1$, which avoids numerical problems if two columns of $\bm{F}_{k}$ are numerically linearly dependent. We refer to [42] for further discussions and more extensive numerical studies regarding the window size and its effect on ADEIM.

3.4.2 Algorithm

Algorithm 1 gives details on the AADEIM approach by summarizing time stepping of a reduced model that uses the proposed adaptive basis updates and adaptive sampling. Helper functions used in Algorithm 1 are given in Appendix B. Inputs of the algorithm are the initial condition $\bm{q}_{0}\in\mathbb{R}^{N}$, the full-model function $\bm{f}$ that describes the underlying dynamical system, and the parameter $\bm{\mu}\in\mathcal{D}$. Parameters of the approach are the dimension $n$ of the reduced space, the time step $w_{\text{init}}$ until which the full model is solved to initialize the reduced model, the window size $w$ for the adaptation, the number of sampling points $m$, the frequency $z$ of the sampling points adaptation, and the rank of the update $r$.

Line 2 solves the full model until time step $w_{\text{init}}$ to compute the corresponding trajectory $\bm{Q}_{w_{\text{init}}}=[\bm{q}_{1},\dots,\bm{q}_{w_{\text{init}}}]\in\mathbb{R}^{N\times w_{\text{init}}}$. Lines 3–6 initialize the reduced model by constructing the basis matrix $\bm{U}_{k}$ from $\bm{Q}_{w_{\text{init}}}$ and the vector of interpolation points $\bm{p}_{k}$. The loop in line 7 iterates over the time steps $k=w_{\text{init}}+1,\dots,K$ at which the reduced model is solved instead of the full model. Line 8 finds the reduced state $\tilde{\bm{q}}_{k}$ that satisfies the reduced model with respect to $\tilde{\bm{f}}_{k}$ corresponding to the current DEIM interpolant with basis matrix $\bm{U}_{k}$ and interpolation points $\bm{p}_{k}$. The if branch on line 10 decides if either the sampling points are adapted or the sampling points from the previous iteration are reused. If the sampling points are adapted, then the full-model function $\bm{f}$ is evaluated at all components and $\bm{S}_{k}$ is derived via the adaptive sampling strategy. Note that $\bm{R}_{k}.\widehat{\leavevmode\nobreak\ \leavevmode\nobreak\ }2$ is Matlab notation and means that the entries of $\bm{R}_{k}$ are squared. If the sampling points are not adapted, then the full-model $\bm{f}$ is evaluated at $\bm{s}_{k-1}$ (see line 17) and all other components are approximated (see line 18). The basis and the interpolation points are adapted with ADEIM at line 21. After the time stepping, the trajectory $\bm{Q}\in\mathbb{R}^{N\times K}$ is returned. The first $w_{\text{init}}$ columns in $\bm{Q}$ are computed with the full model and the subsequent $K-w_{\text{init}}+1$ columns are computed with the adaptive reduced model.

3.4.3 Costs and online efficiency of approach

The proposed approach as summarized in Algorithm 1 has online computational costs that scale with the dimension $N$ of the full model. Thus, the proposed approach is not online efficient in the sense of traditional model reduction; however, the proposed approach still has the potential to achieve speedups compared to static reduced models and full models because operations with costs that scale with $N$ are used economically and the online adaptation allows operating in lower dimensional reduced spaces than with static reduced models: First, the full model is solved for $w_{\text{init}}$ time steps at line 2 in Algorithm 1. The number of initial time steps $w_{\text{init}}$ is typically chosen orders of magnitude smaller than $K$, and thus only few time steps are computed with the full model. Furthermore, since $w_{\text{init}}\ll K$, constructing the POD basis at line 4 is typically cheap as well, even though the costs scale with $N$. Second, adapting the sampling points at lines 11–14 requires evaluating the full-model function $\bm{f}$ at all $N$ components, which incurs costs that scale with $N$. The sampling points are adapted in $K/z$ iterations, which means that $K/z$ evaluations of the full-model function $\bm{f}$ at all $N$ components are necessary. This shows that the speedup of the proposed approach is limited by the number of times the sampling points are adapted; however, evaluating $\bm{f}$ is significantly cheaper than performing a time step with the full model, because the latter requires solving a nonlinear system of $N$ equations, whereas the former typically requires a function evaluation only. Third, the online adaptation typically enables operating in lower dimensional reduced spaces than with static reduced models if the full model exhibits transport-dominated phenomena. Thus, the computational costs of the nonlinear solve for computing the reduced state at the next time step (line 8 in Algorithm 1) is kept low, which has the potential to lead to speedups compared to static reduced models and full models as demonstrated by the numerical results in Section 4.

4 Numerical results

This section demonstrates model reduction based on AADEIM on three numerical examples. First, the toy example based on the advection equation introduced in Section 2.2 is revisited in Section 4.1. Second, AADEIM is demonstrated on the Burgers’ equation with a setup that leads to two interacting waves and time-varying viscosity and time-varying transport-direction coefficients. Third, we consider a model of a rocket combustor and demonstrate that the AADEIM approach achieves significant speedups in contrast to static reduced models that take even longer to run than the full model. All runtime results are computed with a Matlab 2017b implementation.

4.1 Advection equation

Consider the same setup as in Section 2.2. Set $k=25$ and $w=25$ and let $\mathcal{U}_{k}$ be the $n=3$ dimensional DEIM space derived from $[\bm{q}_{k-w+1},\dots,\bm{q}_{k}]\in\mathbb{R}^{N\times w}$ and let $\bm{U}_{k}$ and $\bm{P}_{k}$ be the corresponding basis matrix and interpolation points matrix, respectively. Let now $\bar{\mathcal{U}}_{k+1}$ be the $n=3$ dimensional space derived from $[\bm{q}_{76},\dots,\bm{q}_{100}]$. Figure 3 shows the local coherence $\gamma_{i}(\mathcal{U}_{k})$ and $\gamma_{i}(\bar{\mathcal{U}}_{k+1})$ for $i=1,\dots,1500$; the local coherence is sorted. The corresponding dashed curve are the bounds as in Lemma 2. Let now $\tilde{\bm{F}}_{k}\in\mathbb{R}^{n\times n}$ have i.i.d. standard Gaussian entries and consider $\bm{F}_{k}=\bar{\bm{U}}_{k+1}\tilde{\bm{F}}_{k}$ and the corresponding residual $\bm{R}_{k}=\bm{F}_{k}-\bm{U}_{k}(\bm{P}_{k}^{T}\bm{U}_{k})^{-1}\bm{P}_{k}^{T}\bm{F}_{k}$. The row-wise squared norm $\|\bm{R}^{T}_{k}\bm{e}_{i}\|_{2}^{2}$ of the residual $\bm{R}_{k}$ is plotted in Figure 3, together with the bound given by Lemma 2. The results provide evidence that the decay of the residual is inherited from the decay of the local coherence as shown in Lemma 2.

Consider now Figure 4 that shows the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ of the adapted space $\mathcal{U}_{k+1}$, for ADEIM updates with rank $r=1,2,3$, against the number of sampling points $m$. First, observe that if the rank is $r<3$, then the bound of $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ levels off because the singular values $\sigma_{r+1},\dots,\sigma_{\bar{r}}$ dominate the decay factor $\rho_{k}$ (8). Second, the error $d(\bar{\mathcal{U}}_{k+1},\mathcal{U}_{k+1})$ is bounded by $\rho^{2}_{k}/\sigma_{\text{min}}^{2}(\bm{F}_{k})$ as proved in Proposition 8. Finally, the results show that the fast decay of the error of the adapted space with respect to the number of sampling points $m$ is inherited from the fast decay of the residual (shown as the dashed curve in Figure 4), which demonstrates Proposition 21.

4.2 Burgers’ equation with time-varying viscosity

We now apply AADEIM to the Burgers’ equation with time-varying viscosity and a transport direction that changes with time.

4.2.1 Problem setup

Let $\Omega=(-1,1)\subset\mathbb{R}$ be the spatial domain, set $\bar{\Omega}=[-1,1]$, and let $T=1.5$ be end time. Consider now the Burgers’ equation

[TABLE]

with time $t\in[0,T]$, the solution function $q:\bar{\Omega}\times[0,T]\to\mathbb{R}$, the time-varying viscosity $\nu:[0,T]\to\mathbb{R}$, and the transport direction $\eta:[0,T]\to\{-1,1\}$. The viscosity is

[TABLE]

where $\mu\in\mathcal{D}\subset\mathbb{R}$ is the nominal viscosity parameter. The transport direction is

[TABLE]

The viscosity and transport direction change over time, as shown in Figure 5 for $\mu=3\times 10^{-3}$. The PDE (27) is closed with Dirichlet boundary conditions and the initial condition given by

[TABLE]

The full model is obtained by discretizing the spatial domain $\Omega$ of the PDE (27) on an equidistant grid with $N=1024$ inner grid points. A second-order finite-difference scheme is used. The time domain is discretized with the implicit Euler method and time step size $\delta t=5\times 10^{-5}$. In each time step, Newton’s method is used to solve the corresponding system of nonlinear equations of the form (1). At each time step, 15 Newton iterations are performed with a step size of 1. Jacobians are derived analytically and returned in assembled form to the Newton routine. The linear solves in each Newton step are performed with Matlab’s backslash operator. No parallelization is used besides Matlab’s built-in parallelization for the backslash operator. The same linear-algebra routines, Newton scheme, and parallelization are used for time-stepping in all models. The solution of the full model for $\mu=3\times 10^{-3}$ is shown in Figure 6.

Figure 7a reports the decay of the singular values corresponding to local trajectories of length $w=50$ at time $t=0.77,1.12,1.50$, respectively, and parameter $\mu=3\times 10^{-3}$. The results show an orders of magnitude faster decay of the singular values of the local trajectories than of the global trajectory, which indicates that this problem exhibits a local low-rank structure. Figure 7b shows the decay of the squared residual of approximating the solution at $t=0.77,1.12,1.50$ at the corresponding local spaces of dimension $n=8$ of the trajectories for which the singular values are plotted in Figure 7a. A fast decay of the norm of the rows of the residual is observed.

4.2.2 Performance of AADEIM

We now compare the runtime and accuracy of static reduced models, AADEIM models, and the full model. The static reduced models are derived from the trajectories corresponding to the parameters $\mu\in\{5\times 10^{-3},10^{-3},5\times 10^{-4}\}$, following the procedure outline in Section 2.1. The dimension of the DEIM space is $n$ and the interpolation points are selected with QDEIM [16]. The AADEIM models are initialized with the first $w_{\text{init}}=100$ states computed with the full model. The dimension of the DEIM spaces of the AADEIM models is $n=8$. The interpolation points are selected with QDEIM [16]. The DEIM space $\bm{U}_{k}$ and the DEIM interpolation points $\bm{P}_{k}$ are adapted every other time step. The interpolation points $\bm{P}_{k}$ are adapted by applying QDEIM to the adapted basis. Only updates with rank $r=1$ are applied. The window size is $w=n+1=9$. The sampling points are obtained with the adaptive sampling strategy described in Section 3.3. The sampling points are adapted every $z=5$ iteration, except otherwise noted. The error of the static and the AADEIM models are measured in the Frobenius norm

[TABLE]

where $\bm{Q}(\mu)$ is the trajectory obtained with the full model at parameter $\mu$ and $\tilde{\bm{Q}}(\mu)$ is the trajectory obtained either with a static reduced model or an AADEIM model. For computing $\operatorname{err}(\tilde{\bm{Q}}(\mu))$, only states at time steps $k=1,\dots,1000$ and then for $k>1000$ at every 50-th time step are taken into account.

Let us first consider parameter $\mu=3\times 10^{-3}$. Figure 8a shows the error $\operatorname{err}(\tilde{\bm{Q}}(\mu))$ of the static reduced models with $n\in\{125,175,225\}$ and of the AADEIM models with $m\in\{96,160,224,352,480,608\}$ and dimension $n=8$. The runtime results are obtained on compute nodes with Intel Xeon E5-1660v4 with 64GB RAM. The sampling points in the AADEIM models are adapted very $z$-th time step with $z=5$. The dashed line in Figure 8a marks the runtime of the full model. The plot shows that the static reduced model achieves a speedup compared to the full model in this example; however, more than 100 dimensions are required for the DEIM space to achieve an error below $10^{-1}$. In contrast, the AADEIM model achieves errors below $10^{-1}$ with $n=8$ dimensions and $m=96$ sampling points, which leads to about an order of magnitude speedup compared to the full model. The plot in Figure 8a further indicates that increasing the number of sampling points from $m=96$ to $m=608$ reduces the error from $10^{-1}$ to $10^{-2}$ without a significant increase of the runtime, which provides evidence that evaluating the full-model function $\bm{f}$ is significantly cheaper than solving the full model in this example. Figure 8b shows similar behavior of AADEIM for parameter $\mu=8\times 10^{-4}$. Figure 9a summarizes the runtime of the static reduced model with dimension $n=225$ and the AADEIM model with $m=608$ sampling points, which are required to achieve an error (31) below $10^{-2}$.

We now investigate the effect of adapting the sampling points. The sampling points need to be adapted during the time stepping to ensure that components of the basis matrix corresponding to a high residual are updated. The frequency of adapting the sampling points is controlled by the parameter $z$ and has to be chosen depending on the problem at hand. Just as the time step size of the full model has to be chosen such that the coherent structure moves only so far in the spatial domain per time step that the full model remains stable, the frequency of adapting the sampling points $z$ has to be set such that sampling points can move along with the coherent structure to indicate high-residual components. Thus, if too few adaptations of the sampling points are performed, then high-residual components might be missed. At the same time, too many adaptations can become costly because each adaptation of the sampling points requires an evaluation of the full-model right-hand side function at all components. In Figure 9b, we study the effect of the parameter $z$ on the error (31) of the AADEIM model for parameter $\mu=3\times 10^{-3}$. The sampling points are adapted at every time step ($z=1$), at every $3$-rd time step ($z=3$), and at every $5$-th time step ($z=5$). The number of sampling points $m$ is varied $m=96,160,224,352,480,608$. First, note that, for $z=1$, i.e., the sampling points are adapted at each time step, the runtime is almost constant for an increasing number of sampling points. Second, there is almost no difference in the error for adapting the sampling points at every $3$-rd iteration to adapting at every $5$-th iteration, which means that the sampling points are valid over several time steps in this example.

Figure 10 reports error (31) of the AADEIM model for different initial window sizes $w_{\text{init}}$. The setup of the AADEIM model is the same as for the runtime results reported in Figure 9a with $\mu=3\times 10^{-3}$, i.e., dimension $n=8$, number of sampling points $m=608$, and sampling points are adapted every $5$-th time step. The initial window size $w_{\text{init}}$ is varied in $\{25,50,75,125,150,200,225,250,275,300,500,750,1000\}$. If the window size is set too large, then the local low-rank structure of the snapshots is lost and the correspondingly larger error in the first few time steps accumulates quickly, which can be observed in Figure 10 by noting that the curve increases starting around $w_{\text{init}}=750$. If the window size is chosen too small, then the initial basis might so poorly approximate the dynamics in the first few time steps that the basis adaptation cannot be initialized. Thus, the initial window size $w_{\text{init}}$ is a tuning parameter of the proposed approach; however, as the results in Figure 10 show, the error (31) of the AADEIM model remains near $10^{-2}$ for a wide range of initial window sizes and therefore, in this example, the proposed approach is robust with respect to choices of the initial window sizes.

4.2.3 Performance of adaptive sampling strategy

We now compare adaptive sampling to random uniform sampling, which is used in, e.g., [42]. Figure 11 shows the error (31) of the AADEIM model if the sampling points are selected with the proposed adaptive sampling strategy and random uniform sampling of $\{1,\dots,N\}$. The random uniform sampling is without replacement. The curves in Figure 11 correspond to $m\in\{32,96,160,224,352,480,608,736,N\}$ sampling points. Uniform sampling with $m=32$ leads to an unstable model in case of $\mu=8\times 10^{-4}$ and therefore its error is not plotted. The results in Figure 11 indicate that adaptive sampling reduces the error (31) of the AADEIM model compared to uniform sampling. Improvements of up to two orders of magnitude can be observed. Note that the two sampling schemes coincide if all points $m=N$ are selected. Figure 12 reports the error (31) of the AADEIM model for dimensions $n\in\{4,6,8,10,12\}$ with $m\in\{96,224,480\}$ sampling points (adaptive sampling strategy). The results indicate that a larger dimension of the DEIM space leads to a lower error only if sufficiently many sampling points are selected for the adaptation. For example, the AADEIM model with dimension $n=12$ and $m=480$ sampling points achieves an about one order of magnitude lower error than the AADEIM model with the same dimension and $m=96$ sampling points.

4.2.4 Performance of AADEIM with respect to full-model dimension

We consider the same setup as described in Section 4.2.1, except that we now base the comparison on four full models corresponding to $N\in\{256,512,1024,2048\}$. For all four full models, the time-step size is set to $10^{-5}$ to avoid instabilities of the full model with dimension $N=2048$, in contrast to Section 4.2.2 where the time-step size is $5\times 10^{-5}$. We derive static reduced models and AADEIM models for each of the four full models with the same procedure as described in Section 4.2.1 and Section 4.2.2, respectively. The dimension of the static reduced models is 125 for $N=256$, 225 for $N=512$, 175 for $N=1024$, and 200 for $N=1024$ so that an error (31) of less than $10^{-2}$ is achieved for parameters $\mu=3\times 10^{-3}$ and $\mu=8\times 10^{-4}$. The dimension of the AADEIM models is 8 and the numbers of sampling points are 32, 96, 224, and 736 for full-model dimension 256, 512, 1024, and 2048, respectively, which leads to error (31) of less than $10^{-2}$ for parameter $\mu=3\times 10^{-3}$ and parameter $\mu=8\times 10^{-4}$. As in Section 4.2.2, the adaptive sampling strategy is used and the sampling points are adapted every $5$-th time step. Runtimes were measured on compute nodes with Intel Xeon E5-2690v4 2.6GHz and 32GB RAM. The results in Figure 13a show that speedups of about one order of magnitude are achieved by the AADEIM models in all four cases of $N\in\{256,512,1024,2048\}$ for parameter $\mu=3\times 10^{-3}$. Figure 13b reports the error (31) versus the ratio $m/N$. The reported results demonstrate that the error (31) is similar for all $N\in\{256,512,1024,2048\}$ with respect to the ratio $m/N$. Note that the number of grid points within a neighborhood of the coherent structure grows with the full-model dimension $N$, and thus the number of sampling points $m$ has to grow with $N$ as well. However, the results in Figure 13b further indicate that the number of sampling points $m$ has to grow linearly with $N$ only. Furthermore, note that even though the number of sampling points $m$ grows with $N$, for all $N\in\{256,512,1024,2048\}$, about the same speedup is achieved by AADEIM as shown in Figure 13a. Similar conclusions are drawn from the results corresponding to parameter $\mu=8\times 10^{-4}$ in Figure 13c and Figure 13d.

Figure 14 shows the runtime of the full, static, and AADEIM models for $N=1024$. The runtime of the full model and the static reduced model is partitioned into 4 parts that correspond to evaluating the right-hand side function (“RHS”), computing the Jacobian (“Jacobian”), and solving the corresponding linear system (“solve”) in each Newton iteration. The rest of the runtime is combined into overhead (“overhead”). The runtime of the AADEIM models is split into 7 parts that include the ones mentioned previously and additionally the runtime of sampling the right-hand side function and adapting the sampling points (“sample”), adapting the basis (“adapt $U$”), and computing the QDEIM interpolation points (“adapt $P$”). The overhead includes the runtime of initializing the AADEIM models with $w_{\text{init}}$ many time steps with the full model. The results in Figure 14 show that the full model spends most of the runtime in solving the linear system in each Newton iteration (“solve”). The AADEIM models reduce the runtime of solving the linear system drastically and so achieve about one order of magnitude speedup. When comparing the runtime of the AADEIM models with $m=96,224,736$, one sees that the share of the runtime spent on sampling (“sample”) increases. The runtime due to the adaption of the basis (“adapt $U$”) and computing the QDEIM interpolation points (“adapt $P$”) is low compared to sampling. In particular, this means that using other methods for adapting reduced bases than ADEIM will have a minor effect on the overall runtime. To summarize, the results in Figure 14 indicate that AADEIM achieves runtime speedups in this example because the time for the linear solves in the Newton iterations is reduced compared to full and static reduced models.

4.2.5 Performance of AADEIM with respect to updates with SVD

Consider a procedure Full+SVD that follows AADEIM except that the full-model right-hand side function $\bm{f}$ is evaluated at all components at all time steps and the adapted basis is computed with a thin SVD. This means that the procedure Full+SVD follows Algorithm 1 except that lines 10-19 are replaced with $\bm{F}[:,k]=\bm{f}(\bm{Q}[:,k];\bm{\mu})$ and the new basis is computed with a thin SVD from $\bm{F}_{k}$ instead of ADEIM at line 21. The interpolation points $\bm{P}_{k}$ are computed with QDEIM, cf. Algorithm 2. Thus, the procedure Full+SVD evaluates the right-hand side function $\bm{f}$ at all components rather than using local sampling to update the DEIM basis as AADEIM. The new basis is computed with Matlab’s implementation of the thin SVD.

Consider now the same setup as in Section 4.2.4 to compare Full+SVD with AADEIM. In case of AADEIM, we adapt the basis every other time step as in Section 4.2.4. In case of the procedure Full+SVD, the basis is adapted every third time step so that a comparable error (31) is achieved as with AADEIM; see Figure 15a and Figure 15c for $\mu=3\times 10^{-3}$ and $\mu=8\times 10^{-4}$, respectively. We compare the runtime $T_{A}$ of AADEIM to the runtime $T_{S}$ of Full+SVD. The runtime reduction $T_{R}$ achieved by AADEIM is

[TABLE]

The runtime reduction $T_{R}$ is shown in Figure 15b and Figure 15d. AADEIM achieves a runtime reduction of up to 30% compared to Full+SVD because AADEIM requires computing the full-model right-hand side function $\bm{f}$ at only a few components whereas Full+SVD evaluates the right-hand side function at all components. The costs of adapting the bases by either AADEIM or computing a thin SVD is negligible compared to the costs of to evaluating the right-hand side function, cf. Figure 14.

Figure 16 reports results for Full+SVD if the full-model right-hand side function $\bm{f}$ is evaluated at all components at every third time step only, rather than at all time steps as in the previous numerical experiment. Because of the fewer evaluations of $\bm{f}$, Full+SVD is 2–8% faster than AADEIM; see Figure 16a. However, the runtime decrease of Full+SVD, because of fewer evaluations of $\bm{f}$, leads to poorer updates, which is shown in Figure 16b where now Full+SVD is up to one order of magnitude worse than AADEIM with respect to error (31). Note that AADEIM evaluates $\bm{f}$ at all components at every 5-th time step to update the sampling points and takes sparse evaluations of $\bm{f}$ at all other time steps, which is the same setup as in the previous experiment.

4.3 Combustion model

In this section, we apply AADEIM to a quasi-1D version of a single-element model rocket combustor. The model we use has been developed in the works [53, 20, 21]. The goal is to approximate the growth of the amplitude of pressure oscillations at a monitoring point, which provides critical insights for designing engines that avoid combustion instabilities and unbounded growth of the amplitude of the pressure oscillations. The runtime results are obtained on compute nodes with Intel Xeon E5-1660v4 with 64GB RAM.

4.3.1 Problem setup and full model

The problem setup and full model follows [61, 57]. The problem setup consists of three parts, namely the oxidizer post, the combustion chamber, and the exit nozzle, see Figure 17. The oxidizer is induced and meets the fuel at the back-step, where it reacts instantaneously. The combustion products exit the chamber through the nozzle. The combustion follows a one-step reaction model

[TABLE]

where the fuel is gaseous methane and the oxidizer is a mixture of oxygen and water. Details on the operating conditions are given in [61, Table 2]. The governing equations are in conservative form

[TABLE]

where

[TABLE]

where $\rho$ is the density, $v$ is the velocity, $E$ is the total internal energy, $Y_{\text{ox}}$ is the oxidizer mass fraction, $A$ is the cross sectional area of the fuel duct, and $p$ is the pressure. The source terms $s_{A}$ and $s_{g}$ are given in [61, equation (2)]. The source term $s_{q}(\mu)=[0,0,q^{\prime}(\mu),0]^{T}$ models the heat release, where the parameter $\mu\in\mathcal{D}=[2,4.2]\subset\mathbb{R}$ in $q^{\prime}(\mu)$ controls the amplification of the heat release, see [61, equation (5)]. Following [61], a steady-state solution is first obtained by ignoring the source term $s_{q}(\mu)$. The steady-state solution is then used as initial condition for computing the time-dependent solution that takes the source term $s_{q}(\mu)$ into account. The initial condition is perturbed to trigger an instability, which depends on the heat-release parameter $\mu$. The spatial domain is discretized on 300 equidistant grid points. There are four degrees of freedom at each grid point (density, velocity, energy, mass fraction), and so the full model has a total of $N=1200$ degrees of freedom. The spatial discretized is a first-order finite-difference scheme. Time is discretized with a fourth-order implicit scheme and time step size $\delta t=10^{-7}$ with end time $T=10^{-1}$. Newton’s method is used to solve the corresponding system of nonlinear equations at each time step. Jacobians are derived analytically and passed in assembled form to the Newton solver. Matlab’s backslash operator is used to solve the linear system in each Newton step. The Newton scheme takes 10 iterations and uses a line search based on the Armijo condition [35, p. 33] with control parameter $10^{-4}$ as recommended in [35, p. 33]. The rest of the setup of the Newton scheme is the same as in Section 4.2. The same Newton scheme is used in all models. The pressure is monitored at spatial coordinate $x=0.3683$, see Figure 17. Figure 18 shows the pressure at the monitoring point for parameters $\mu\in\{2.4,3.0,3.8,4\}$. The solutions converge to a steady state for $\mu=2.4$. A limit cycle oscillation is entered for $\mu=3.8$. A combustion instability is observed for $\mu=4.0$.

4.3.2 Performance of AADEIM

We compare static reduced models, AADEIM models, and the full model. The static reduced model is derived from the trajectories corresponding to the parameters $\mu\in\{2,2.44,2.88,3.32,3.76,4.2\}$, which are the six equidistant parameters in the parameter domain $\mathcal{D}=[2,4.2]$. A separate DEIM basis of dimension $n$ is computed for each degree of freedom, see, e.g., [64, Section 2.2]. The DEIM interpolation points are derived with QDEIM, where $n$ points are derived for each of the four DEIM basis and then the union of all four sets of points is used as the set of DEIM interpolation points. The dimension of the AADEIM model is $n=8$ and the window size is $w=n+1$ in the following. The AADEIM model is initialized with the full-model states obtained until $t=1.6\times 10^{-4}$, where only every $50$-th state is used so that $w_{\text{init}}=4\times n=32$. The sampling points are derived for each degree of freedom separately, then they are ranked by how often each sampling point has been selected, and then the same $m$ sampling points are used for all four bases that have been selected the most. This means that a total of $4\times m$ components of the full-model residual are computed in each iteration. The basis is adapted at every time step and the sampling points are adapted very other time step ($z=2$). We measure the error of the pressure at the monitoring point. Let $\bm{y}(\mu)=[y_{1}(\mu),\dots,y_{K}(\mu)]^{T}\in\mathbb{R}^{K}$ be the trajectory of the pressure at the monitoring point computed with the full model. Then, we measure the error

[TABLE]

where $\tilde{\bm{y}}(\mu)$ is the trajectory computed with either a static or an AADEIM model. The rest of the setup is the same as in Section 4.2.

Figure 19 reports the error (33) of the reduced models and their runtime. The error and runtime of the static reduced model is plotted for $n\in\{60,70,80,90,100\}$. The results for the AADEIM model are reported for $m\in\{20,30,40,50\}$ and $n=8$. The parameter is set to $\mu=3.8$. The AADEIM model achieves a speedup of about 6 compared to the full model in this example. The static reduced model is slower than the full model. Figure 20a shows the speedup of the AADEIM model with $m=50$ and for a parameter sweep over $\mu\in\{2.4,3.0,3.8,4.0\}$. The dimension of the static reduced model is $n=225$. The static and the AADEIM model achieve about the same error for all four parameters. Note that the four parameters lead to significantly different behaviors in the solutions, see Figure 18. The AADEIM model achieves a significant speedup compared to the full model, whereas the static reduced model is slower than the full model. Figure 20b demonstrates that the proposed adaptive sampling scheme is orders of magnitude more efficient than uniform sampling without replacement. Uniform sampling requires at least $m=175$ sampling points per degree of freedom to prevent the Newton method from diverging and to achieve an error (33) of about $10^{-2}$. With the proposed adaptive sampling scheme, the AADEIM approach achieves an error (33) of about $10^{-4}$ with $m=50$ sampling points per degree of freedom.

5 Conclusions

The proposed AADEIM approach demonstrates that transport-dominated problems have a rich local structure that can be exploited to construct efficient reduced models. We exploit locality in time via adaptive basis updates and locality in space via adaptive sampling. An analysis establishes a connection between the local coherence properties of reduced spaces and the number of samples that are required to adapt the basis. The faster the local coherence decays, the fewer samples are required to adapt the reduced spaces. Reduced models built with AADEIM are implicitly parametrized, which means that there is no offline phase to construct reduced models, rather the basis is adapted online to changes in the parameter. Numerical results demonstrated that AADEIM is applicable to a wide range of problems and faithfully approximates behavior that changes significantly with parameters. At the same time, AADEIM achieves significant runtime speedups compared to full and traditional, static reduced models.

There are many directions for future work. First, the dimension of the reduced space, as well as the window size, can be adapted over time. For example, in problems with strongly time-varying coefficients or transport behavior that changes over time, time-varying windows and time-varying dimensions of the reduced spaces could help to further reduce the runtime. Second, the proposed approach relates the number of sampling points for the basis updates to the locality of the moving coherent structure. It would be interesting to compare to other adaptive basis schemes, e.g., incremental SVDs [8], and to establish how many sampling points are required per adaptation step.

Acknowledgments

The author would like to thank Nina Beranek (University of Ulm, Germany) for carefully reading an earlier version of this manuscript and for reporting typos and helpful comments.

Appendix A Analytic example with local low-rank structure

We analyze analytically an example to demonstrate the local low-rank structure of transport-dominated problems. The example follows [54, Example 2.5]. Consider the function

[TABLE]

in the spatial domain $x\in[-5,5]$ and $t\in[0,T]$ with $T=1$. There is no $n$-dimensional Lagrangian space $\mathcal{U}$ for which the error $\sup_{t\in[0,T]}\inf_{q^{*}\in\mathcal{U}}\|q(\cdot,t)-q^{*}\|_{L^{2}(-5,5)}$ decays faster than linearly in $1/\sqrt{n}$, which means there is no space spanned by $q(\cdot,t_{1}),\dots,q(\cdot,t_{n})$ with $n$ pairwise distinct $t_{i}\in[0,T]$ with $i=1,\dots,n$ that achieves a faster error decay than $1/\sqrt{n}$. Instead of considering the whole time domain $[0,T]$, let us now consider $[0,T/\zeta]$ with $\zeta>0$. Let $\mathcal{U}_{\zeta}$ be the space spanned by $q(\cdot,t_{i})$ with $t_{i}=ih$ for $i=1,\dots,n$ and $h=(T/\zeta)/n$. For $t\in[0,T/\zeta]$, the space $\mathcal{U}_{\zeta}$ achieves

[TABLE]

and thus $\sup_{t\in[0,T/\zeta]}\inf_{q^{*}\in\mathcal{U}_{\zeta}}\|q(\cdot,t)-q^{*}\|_{L^{2}(-5,5)}\leq\sqrt{h}=\sqrt{T/(\zeta n)}$. The rate of the error decay can be increased by letting $\zeta$ depend on $n$. For example, setting $\zeta=\mathrm{e}^{n}$ gives $\sup_{t\in[0,T/\zeta]}\inf_{q^{*}\in\mathcal{U}_{\zeta}}\|q(\cdot,t)-q^{*}\|_{L^{2}(-5,5)}\in\mathcal{O}(\mathrm{e}^{-n})$, which shows that a local low-rank structure can be recovered if (34) is approximated locally in time.

Appendix B Helper functions for Algorithm 1