A Goal-Oriented Adaptive Discrete Empirical Interpolation Method

TL;DR

This paper introduces an adaptive DEIM algorithm and a-posteriori error estimation techniques to improve the accuracy of quantities of interest in reduced order models, validated on fluid dynamics models.

Contribution

It develops a novel adaptive DEIM method combined with a-posteriori error estimation for enhanced reduced order model accuracy.

Findings

Error estimates closely match theoretical predictions

Adaptive DEIM improves quantity of interest accuracy

Validated on 1D-Burgers and Shallow Water models

Abstract

In this study we propose a-posteriori error estimation results to approximate the precision loss in quantities of interests computed using reduced order models. To generate the surrogate models we employ Proper Orthogonal Decomposition and Discrete Empirical Interpolation Method. First order expansions of the components of the quantity of interest obtained as the product between the components gradient and model residuals are summed up to generate the error estimation result. Efficient versions are derived for explicit and implicit Euler schemes and require only one reduced forward and adjoint models and high-fidelity model residuals estimation. Then we derive an adaptive DEIM algorithm to enhance the accuracy of these quantities of interests. The adaptive DEIM algorithm uses dual weighted residuals singular vectors in combination with the non-linear term basis. Both the a-posteriori error estimation results and the adaptive DEIM algorithm were assessed using the 1D-Burgers and Shallow Water Equation models and the numerical experiments shows very good agreement with the theoretical results.

Figures

Equations

$${\mathbf{x}}_{i+1}={\mathcal{M}}_{i,i+1}\left({\mathbf{x}}_i,\mu \right),i=0,..,N_t-1,\mu \in \tilde{\mathcal{P}},$$

$${\tilde{\mathbf{x}}}_{i+1}=U_{\mu }^T\cdot {\mathcal{M}}_{i,i+1}\left(U_{\mu }{\tilde{\mathbf{x}}}_i,\mu \right),i=0,..,N_t-1,\mu \in \tilde{\mathcal{P}}.$$

$$\tilde{N}(\tilde{\mathbf{x}})\approx \underset{k\times N_{\mathrm{state}}}{\underbrace{U_{\mu }^T}}\cdot \underset{N_{\mathrm{state}}\times 1}{\underbrace{\mathbf{F}(U_{\mu }{\tilde{\mathbf{x}}}_i,\mu )}},$$

$$\tilde{N}(\tilde{\mathbf{x}})\approx \underset{\mathrm{precomputed}k\times m}{\underbrace{U_{\mu }^T{V}_{\mu }(P^T{V}_{\mu }{)}^{-1}}}\cdot \underset{m\times 1}{\underbrace{P^T\mathbf{F}\left(U_{\mu }\tilde{\mathbf{x}},\mu \right)}}$$

$$\mathcal{Q}:{\mathbb{R}}^{N_{\mathrm{state}}\times \tilde{\mathcal{P}}}\to \mathbb{R},\mathcal{Q}\left({\mathbf{x}}_0,\mu \right)=\sum _{i=0}^{N_t}{r}_i\left({\mathbf{x}}_i,\mu \right),$$

$$\epsilon (\mu )=\mathcal{Q}\left(\mathbf{x},\mu \right)-\mathcal{Q}\left(\hat{\mathbf{x}},\mu \right),$$

$$\underset{{\mathbf{x}}_0}{\min }\mathcal{Q}\left({\mathbf{x}}_0,\mu \right)=\sum _{i=0}^{N_t}{r}_i\left({\mathbf{x}}_i,\mu \right),$$

$${\mathbf{x}}_{i+1}={\mathcal{M}}_{i,i+1}\left({\mathbf{x}}_i,\mu \right),i=0,..,N_t-1.$$

$$\text{Forward model: }{\mathbf{x}}_{i+1}={\mathcal{M}}_{i,i+1}\left({\mathbf{x}}_i,\mu \right),i=0,..,N_t-1,$$

$$\begin{array}{rl}\text{Adjoint model: }& {\lambda }_N=-{\left(\frac{\partial r_{N_t}}{\partial {\mathbf{x}}_{N_t}}\right)}^T({\mathbf{x}}_{N_t},\mu ),\\ & {\lambda }_i={\mathbf{M}}_{i+1,i}^{\ast }{\lambda }_{i+1}-{\left(\frac{\partial r_i}{\partial {\mathbf{x}}_i}\right)}^T({\mathbf{x}}_i,\mu ),i=N_t-1,..,0,\end{array}$$

$$\text{Reduced gradient: }{\nabla }_{{\mathbf{x}}_0^{\ast }}\mathcal{L}=-{\lambda }_0=0,$$

$$\mathcal{L}({\mathbf{x}}_{0})=\sum_{i=0}^{N_{t}}~{}r_{i}\left(\mathbf{x}_{i},\mu\right)+\sum_{i=0}^{N-1}{\mathbf{\lambda}}_{i+1}^{T}\big{(}{\mathbf{x}}_{i+1}-\mathcal{M}_{i,i+1}\left({\mathbf{x}}_{i},\mu\right)\big{)},$$

$$\begin{split}\delta\mathcal{L}({\mathbf{x}}_{0})&=-\sum_{i=0}^{N_{t}}\frac{\partial r_{i}}{\partial\mathbf{x}_{i}}(\mathbf{x}_{i},\mu)^{T}\delta\mathbf{x}_{i}+\\ &\sum_{i=0}^{N_{t}-1}\lambda_{i+1}^{T}\big{(}\delta{\mathbf{x}}_{i+1}-{\mathbf{M}}_{i,i+1}\delta{\mathbf{x}}_{i}\big{)}+\sum_{i=0}^{N_{t}-1}\delta{\mathbf{\lambda}}_{i+1}^{T}\big{(}{\mathbf{x}}_{i+1}-\mathcal{M}_{i,i+1}\left({\mathbf{x}}_{i},\mu\right)\big{)},\\ \end{split}$$

$${\mathbf{M}}_{i,i+1}=\frac{\partial {\mathcal{M}}_{i,i+1}}{\partial {\mathbf{x}}_i}({\mathbf{x}}_i,\mu )\in {\mathbb{R}}^{N_{\text{state}}\times N_{\text{state}}}.$$

$$⟨{\mathbf{M}}_{i,i+1}{z}_1,z_2{⟩}_{{\mathbb{R}}^{N_{\text{state}}}}=⟨z_1,{\mathbf{M}}_{i+1,i}^{\ast }{z}_2{⟩}_{{\mathbb{R}}^{N_{\text{state}}}},\forall z_1,z_2\in {\mathbb{R}}^{{\mathcal{N}}_{\text{state}}},$$

$$\begin{split}&\delta\mathcal{L}({\mathbf{x}}_{0})=\delta\mathbf{x}_{N_{t}}^{T}\big{(}\frac{\partial r_{N_{t}}}{\partial\mathbf{x}_{N_{t}}}(\mathbf{x}_{N_{t}})+{\mathbf{\lambda}}_{N_{t}}\big{)}+\sum_{i=1}^{N-1}\delta{\mathbf{x}}_{i}^{T}\bigg{(}{\mathbf{\lambda}}_{i}-{\mathbf{M}}_{i+1,i}^{*}\mathbf{\lambda}_{i+1}+\frac{\partial r_{i}}{\partial\mathbf{x}_{i}}(\mathbf{x}_{i})\bigg{)}\\ &-\delta{\mathbf{x}}_{0}^{T}\bigg{(}{\mathbf{M}}_{1,0}^{*}\mathbf{\lambda}_{1}-\frac{\partial r_{0}}{\partial\mathbf{x}_{0}}(\mathbf{x}_{0})\bigg{)}+\sum_{i=0}^{N-1}\delta{\mathbf{\lambda}}_{i+1}^{T}\big{(}{\mathbf{x}}_{i+1}-\mathcal{M}_{i,i+1}\left({\mathbf{x}}_{i},\mu\right)\big{)}.\end{split}$$

$${\nabla }_{{\mathbf{x}}_0}\mathcal{Q}={\nabla }_{{\mathbf{x}}_0}\mathcal{L}=-{\lambda }_0.$$

$$\mathcal{Q}({\bar{\mathbf{x}}_{0}},\mu)-\mathcal{Q}({\mathbf{x}_{0}},\mu)\approx-{\mathbf{\lambda}}_{0}^{T}\cdot\big{(}{\bar{\mathbf{x}}_{0}}-{\mathbf{x}_{0}}\big{)}.$$

$${\mathbf{x}}_{\ell }^i={\mathcal{M}}_{i,\ell }\left({\hat{\mathbf{x}}}_i,\mu \right),\ell =i+1,\dots ,N_t,i=0,\dots ,N_t-1.$$

$$\epsilon (\mu )\approx -\sum _{i=0}^{Nt}{\hat{\lambda }}_i^T\cdot \Delta {\mathbf{x}}_i,$$

$$\Delta {\mathbf{x}}_0={\mathbf{x}}_0-{\hat{\mathbf{x}}}_0,\Delta {\mathbf{x}}_i={\mathbf{x}}_i^{i-1}-{\hat{\mathbf{x}}}_i,i=1,\dots ,N_t,$$

$${\hat{\lambda }}_{N_t}=-(\partial r_{N_t}/\partial {\mathbf{x}}_{N_t}{)}^T({\hat{\mathbf{x}}}_{Nt},\mu ).$$

$${\mathcal{Q}}^i({\hat{\mathbf{x}}}_i,\mu )=r_i({\hat{\mathbf{x}}}_i,\mu )+\sum _{j=i+1}^{N_t}{r}_j({\mathbf{x}}_j^i,\mu ),i=0,1,..,N_t-1.$$

$$\begin{array}{rl}& {\mathcal{Q}}^0({\mathbf{x}}_0,\mu )-{\mathcal{Q}}^0({\hat{\mathbf{x}}}_0,\mu )=r_0({\mathbf{x}}_0,\mu )-r_0({\hat{\mathbf{x}}}_0,\mu )+\sum _{j=1}^{Nt}{r}_j({\mathbf{x}}_j,\mu )-r_j({\mathbf{x}}_j^0,\mu )\\ & \approx -{\hat{\lambda }}_0^T\cdot \Delta {\mathbf{x}}_0,\end{array}$$

$$\begin{array}{rl}& {\mathcal{Q}}^i({\mathbf{x}}_i^{i-1},\mu )-{\mathcal{Q}}^i({\hat{\mathbf{x}}}_i,\mu )=r_i({\mathbf{x}}_i^{i-1},\mu )-r_i({\hat{\mathbf{x}}}_i,\mu )+\sum _{j=i+1}^{Nt}{r}_j({\mathbf{x}}_j^i,\mu )-r_j({\mathbf{x}}_j^i,\mu )\\ & \approx -{\hat{\lambda }}_i^T\cdot \Delta {\mathbf{x}}_i,i=1,..,N_t-1,\end{array}$$

$$r_{N_t}({\mathbf{x}}_{N_t}^{N_t-1},\mu )-r_{N_t}({\hat{\mathbf{x}}}_{N_t},\mu )\approx \frac{\partial r_{N_t}}{\partial {\mathbf{x}}_{N_t}}({\hat{\mathbf{x}}}_{Nt},\mu )\Delta {\mathbf{x}}_{Nt}.$$

$$\epsilon (\mu )\approx -\sum _{i=0}^{Nt}{\hat{\lambda }}_i^T\cdot \Delta {\mathbf{x}}_i.$$

$${\mathbf{x}}_{i+1}={\mathbf{x}}_i+h\mathbf{F}({\mathbf{x}}_i,\mu ),i=0,..,N_t-1,$$

$${\mathbf{x}}_{i+1}^i$$

$${\tilde{\mathbf{x}}}_{i+1}$$

$$\Delta {\mathbf{x}}_{i+1}$$

$${\varphi }_{i+1}(\mu )=\varphi ({\hat{\mathbf{x}}}_i,{\hat{\mathbf{x}}}_{i+1},\mu )={\hat{\mathbf{x}}}_{i+1}-{\hat{\mathbf{x}}}_i-h\mathbf{F}({\hat{\mathbf{x}}}_i,\mu ),i=0,..,N_t-1.$$

Full text

Keywords— proper orthogonal decomposition; discrete empirical interpolation method; a-posteriori error estimates; shallow water equations

1 Introduction

The major issue in large scale complex modelling is that of reducing the computational cost while preserving numerical accuracy. Among the model reduction techniques, reduced basis [7, 33, 59, 68, 24] and Proper Orthogonal Decomposition (POD) [41, 48, 39, 49] provide an efficient means of deriving the reduced basis for high-dimensional non-linear flow systems. Data analysis is conducted to extract basis functions, from experimental data or detailed simulations of high-dimensional systems, for subsequent use in Galerkin or Petrov-Galerkin projections that yield low dimensional dynamical models.

However due to the model non-linearities the computational complexity of the reduced order models still depends on the number of variables of the high-fidelity model. The current literature presents several ways to avoid this issue such as the empirical interpolation method (EIM) [7] and its discrete variant DEIM [16, 18, 17], best points interpolation method [55]. Missing point estimation [6] and Gauss-Newton with approximated tensors [14, 15] methods are relying upon the gappy POD technique [29] and were developed for the same reason. A comparative study between missing point estimation method, gappy POD and DEIM is available in [26].

A-posteriori error estimates of DEIM reduced non-linear dynamical system based on logarithmic Lipschitz constants are available in the literature [86] . Additional state space error estimates [18] are shown to be proportional to the sums of the singular values corresponding to the neglected POD basis vectors both in Galerkin projection of the reduced system and in the DEIM approximation of the non-linear term.

The parameter reduced order modeling poses even more difficulties since one may require to construct a global basis to enable both accurate and fast reduced order models along the entire parametric domain. Recently adaptive techniques capable to enhance the accuracy of the reduced order models or quantities of interest depending on such surrogate models have been proposed. For example, by assigning larger weights to samples that are more important or have a higher probability of occurring allow the construction of smaller global bases for both state and non-linear terms [20, 21]. In the data assimilation field, the sensitivities of a cost functional with respect to the time-varying model state is used to define appropriate weights and to implement a dual weighted POD method for order reduction [22]. Effective exploration of the parameter space by adaptive grids based on a-posteriori error estimates [36] brought improvement over the fixed and uniformly refined grid approaches. Dual techniques [80] for a-posteriori error estimates enabled the construction of a goal-oriented global parametric basis. The idea of a spanning ROM [81] suggests interpolating the associated high-fidelity solutions for a new parameter configuration and then extracting the singular vectors.

To generate reduced order models being simultaneously accurate and performant local strategies have been designed. The parametric or time domains are partitioned into several sub-regions and local bases are constructed for both the state variables and non-linear terms. Local reduced operators associated with the EIM [27] or DEIM [64] approximations and POD [65] or reduced basis models [25] are computed off-line and are properly selected during the on-line phase. Dictionary methods [42, 50] instead of building the reduced basis during the off-line phase construct a small parameter adapted basis using a Greedy procedure during the on-line stage. A different approach solves a residual minimization problem to obtain a low dimensional approximation of an incremental solution on-line [4] without the need of computing reduced operators during the off-line stage.

Another form of adaptivity involves interpolation of already existing bases in the matrix space [47], principal vector space [47] and in the space tangent to the Grassmann manifold [2] to generate a reduced order basis for a new parameter configuration. However these techniques still require computation of reduced operators thus they can be classified as off-line approaches. To mitigate this inefficiency directly interpolating the system matrices of the local reduced models [58] has been proposed in the space tangent to the Grassmann [3, 89] and Riemannian [23] manifolds.

Several approaches that incorporate new data online and rebuild the reduced system are available in the reduced order optimization field [5, 82, 78, 88] and parametric reduced order modeling [60].

In this paper we introduce a-posteriori error estimation results to approximate the errors in quantities of interests resulted from using reduced POD/DEIM reduced order models. First order expansions of the components of the quantity of interest obtained as the product between the components gradient and model residuals are summed up to generate the error estimation result. Efficient versions are derived for explicit and implicit Euler schemes and require only one reduced forward and adjoint models and high-fidelity model residuals estimation.

Later, we propose a novel strategy that incorporates on-line dual weighted residuals associated with some quantity of interest and modifies the location of the DEIM points corresponding to the DEIM non-linear reduced order terms approximations. In this way we enhance the accuracy of the quantity of interest computed via reduced order models. We show that our methodology works also for different parameter configurations as seen in the numerical experiments section using 1D-Burgers and Shallow Water Equations models. The key ingredients are: (1) an a-posteriori error estimation result based on necessary optimality conditions of a constrained optimization problem; (2) an adaptive greedy algorithm that makes use of both non-linear terms and dual weighted residuals bases, respectively.

In contrast to the work proposed in [62] where both non-linear reduced basis and DEIM interpolation points are adapted with additive low-rank updates, we are only adapting the location of the interpolation points by taking into account the dual weighted residuals. Our approach follows the research developed in [52] where dual weighted method is used to adaptively resize the number of basis vectors and the length of the time step to satisfy a given error tolerance.

This paper is organized as follows. Section 2 reviews the Proper Orthogonal Decomposition and Discrete Empirical Interpolation methods and Section 3 formulates the problem of constructing reduced order models with adaptive DEIM interpolation points. In Section 4 we introduce the a-posteriori error estimation approach for the error in a quantity of interest. Fast computational strategies based on a reduced order adjoint model are derived and enable the use of the a-posteriori error estimations for different parametric configurations. Section 5 presents a greedy adaptive algorithm that uses the dual weighted residuals to update the location of DEIM interpolation points. Section 6 presents numerical results with the 1D-Burgers and Shallow Water Equations models. In Section 8 we conclude the paper with some remarks.

2 Reduced order modeling

Reduced order modeling represents the dynamics of large-scale systems using only a smaller number of variables and reduced order basis functions. We consider here the construction of reduced order models via Proper Orthogonal Decomposition (POD) and Galerkin projection. To mitigate a known deficiency of POD Galerkin, we make use of Discrete Empirical Interpolation Method to approximate the non-linear terms.

2.1 Proper Orthogonal Decomposition

Proper Orthogonal Decomposition (Karhunen-Loève decomposition) has been proposed in [41, 48, 45, 57] for an infinite dimensional framework. Latter on, the finite dimensional case leads to the development of principal component analysis [40] in the statistical literature, inspired from the early work of Pearson [61] and Hotelling [39]. The approach is named empirical orthogonal function decomposition in oceanography and meteorology [49], and factor analysis in psychology and economics [31].

POD can be thought of as a Galerkin approximation in the spatial variable built from functions corresponding to the solution of the physical system at specified time instances, as obtained by the method of snapshots [72, 73, 74]. POD has been used successfully in numerous applications such as unsteady viscous fl‚ows [28], compressible flows [67], computational fluid dynamics [46, 66, 84, 56], turbulent flows [71, 83], and aerodynamics [11].

We consider a general discrete dynamical system

[TABLE]

where $\mathbf{x}_{i}\in\mathbb{R}^{N_{\rm{state}}}$, $i=0,\ldots,N_{t}$, denotes the state at time $t_{i}$, with $N_{\rm{state}}$ being the total number of discrete model variables per time step, and $N_{t}\in\mathbb{N},~{}N_{t}>0$ being the number of time steps. Here $\mu$ is a vector of parameters that characterizes the physical properties of the flow, e.g., the Coriolis force in a Shallow Water Equations model.

For a given parameter configuration we select an ensemble of $N_{t}$ time instances (snapshots) ${\bf x}_{t_{1}}^{\mu},...,{\bf x}_{t_{N_{t}}}^{\mu}\in\mathbb{R}^{N_{\rm{state}}}$ from the solution of (1). The POD method chooses an orthonormal basis $U_{\mu}=[{\bf u}_{i}^{\mu}]_{i=1,..,k}\in\mathbb{R}^{{N_{\rm{state}}}\times k}$ such that the mean square error between ${\bf x}(\mu,t_{i})$ and the projection ${\bf x}_{\textsc{pod}}^{\mu}(t_{i})=U_{\mu}\,{\bf\tilde{x}}(\mu,t_{i}),$ ${\bf\tilde{x}}(\mu,t_{i})\in\mathbb{R}^{k},~{}i=1,..,N_{t},$ is minimized on average. The POD space dimension $k\ll{N_{\rm{state}}}$ is appropriately chosen to capture the dynamics of the flow as described by Algorithm 1.

Another way to compute the POD basis is to make use of the eigenvalue decomposition applied to the correlation matrix [77, Alg. 1]. However the SVD-based POD basis construction is more computationally efficient.

The Galerkin approach projects the full order state $\mathbf{x}$ and the full order model equations (1) onto the space spanned by the POD basis elements to obtain the following reduced model:

[TABLE]

The efficiency of the POD-Galerkin techniques is limited to linear or bilinear terms, since the projected non-linear terms at every discrete time step still depend on their evaluation in the full model space, e.g.,

[TABLE]

where $\mathbf{F}:\mathbb{R}^{N_{\rm{state}}}\to\mathbb{R}^{N_{\rm{state}}}$ is the non-linear component of the model (1). In case of polynomial non-linearities, the tensorial POD technique [77] can be employed to remove the dependence on the dimension of the full order system by manipulating the order of computing.

To mitigate this inefficiency we make use of the Discrete Empirical Interpolation Method [17, 75] that can handle efficiently any type of non-linearity and provides a natural framework for adaptive reduced order modeling.

2.2 Discrete Empirical Interpolation Method

The empirical interpolation method [7, 8] and its discrete version DEIM [17] were designed to estimate non-linear terms allowing for an effectively affine offline-online computational decomposition in the context of reduced order models. Both interpolation methods provide an efficient way to approximate non-linear functions in the continuous and discrete frameworks. They were successfully used in the POD framework with finite difference, finite element, and finite volume discretization methods. A description of EIM in connection with the reduced basis framework and a-posteriori error bounds can be found in [51, 34].

The DEIM implementation is based on a POD approach combined with a greedy algorithm. For $m\ll N_{\rm{state}}$, the POD/DEIM approximation of (3) is

[TABLE]

where $V_{\mu}\in\mathbb{R}^{N_{\rm{state}}\times m}$ collects the first $m$ POD basis modes of non-linear function ${\bf F}$, while $P\in\mathbb{R}^{N_{\rm{state}}\times m}$ is the DEIM interpolation selection matrix. The core of the DEIM procedure is an iterative greedy procedure given in Algorithm 2 that inductively constructs $P$ from the linearly independent set $V_{\mu}$ [17, section 3]. Specifically, Algorithm 2 determines the indices $\rho_{1},\rho_{2},...,\rho_{m}$, and constructs $P=[e_{\rho_{1}},..,e_{\rho_{m}}]\in\mathbb{R}^{n\times m}$, where $e_{\rho_{i}}\in\mathbb{R}^{n}$ is the $\rho_{i}$-th column of the identity matrix.

The algorithm first searches for the largest absolute value entry of the first POD basis $|{\bf v}_{1}|$, and the corresponding index represents the first DEIM interpolation index $\rho_{1}\in\{1,2,..,n\}$. The remaining interpolation indices $\rho_{\ell},~{}\ell=2,3..,m$, are selected so that each of them corresponds to the entry with the largest absolute value in $|{\bf r}|$. The vector ${\bf r}$ can be viewed as the residual or the error between the next basis vector ${\bf v}_{\ell},~{}l=2,3..,m$ and its approximation $V_{\mu}{\bf c}$ from interpolating the previous basis vectors $\{{\bf v}_{1},{\bf v}_{2},..,{\bf v}_{\ell-1}\}$ at the previously determined indices ${\rho_{1}},{\rho_{2}},..,{\rho_{\ell-1}}$. The linear independence of the vectors $\{{\bf v}_{\ell}\}_{\ell=1}^{m}$ guarantees that ${\bf r}$ is a nonzero vector at each iteration, and that the output indices $\{\rho_{\ell}\}_{\ell=1}^{m}$ are not repeating.

The spectrum analysis of the non-linear snapshots correlation matrix offers guidance to choosing the number of DEIM interpolation points. We are particularly interested in the eigenvalues rate of descent. Lemma $(3.2)$ in [17, section 3.2] provides an error bound for the DEIM approximation (4) that can be estimated by the largest POD eigenvalue of the snapshots correlation matrix not taken into account by POD basis $V_{\mu}$.

Recent developments of DEIM include rigorous state space error bounds [18], a-posteriori error estimation [85], and applications to 1D FitzHugh-Nagumo model [17], 1D simulating neurons model [43], 1D non-linear thermal model [38], 1D Burgers equation [1, 16], 2D non-linear miscible viscous fingering in porous medium [19], oil reservoirs models [79], and 2D Swallow Water Equations model (SWE) [75]. We emphasize that only few POD/DEIM studies with finite elements or finite volume methods have been performed, e.g., for electrical networks [37] and for a 2D ignition and detonation problem [10]. Flow past a cylinder simulations using a hybrid reduced approach combining the quadratic expansion method and DEIM are performed in [87].

3 Problem formulation

We are interested in a particular aspect of the solution of (1) defined by the smooth scalar function $\mathcal{Q}(\mathbf{x},\mu)$ with

[TABLE]

with $r_{i}:\mathbb{R}^{N_{\rm{state}}\times\mathcal{\tilde{P}}}\to\mathbb{R},~{}i=0,\ldots,N_{t}$. We call (5) the quantity of interest (QoI).

A POD/DEIM reduced order model described in (2) and (4) is capable to decrease the computational complexity of evaluating $\mathcal{Q}$ in (5) by orders of magnitudes. However the reduced order approximation leads to an error in the computed QoI denoted by

[TABLE]

where $\widehat{\mathbf{x}}=\{\widehat{\mathbf{x}}_{0},~{}\widehat{\mathbf{x}}_{1},..,\widehat{\mathbf{x}}_{N_{t}}\}$ is the reduced order solution projected onto the full space, i.e. $\widehat{\mathbf{x}}_{i}=U_{\mu}\,\mathbf{\tilde{x}}_{i}$, $i=0,..,N_{t}$.

We seek to construct adaptive reduced order models that allow to accurately estimate the QoI (5). Specifically, we propose a procedure to adaptively modify the locations of DEIM interpolation points associated with the model non-linear term (4) such as to decrease the error in the QoI (6). Our approach uses of an efficient a-posteriori error estimate based on the solution of a reduced order adjoint model, and the residual of the high-fidelity model computed with the projected reduced order solution. The dual weighted residuals (the Hadamard products of the projected reduced order adjoint solution and the model residuals) are used to design the new locations for the DEIM points to increase the accuracy of $\mathcal{Q}\left(\widehat{\mathbf{x}},\mu\right)$.

4 A-posteriori error estimates

4.1 Gradient of the quantity of interest

Our a-posteriori error estimate framework requires the first order necessary optimality conditions of the following constrained optimization problem

[TABLE]

Theorem 1** (Minimization of the QoI)**

Assume that the model operators $\mathcal{M}_{i,i+1}:\mathbb{R}^{N_{\rm{state}}\times\mathcal{\tilde{P}}}\to\mathbb{R}^{N_{\rm{state}}},~{}i=0,..,N_{t}-1,$ are of class $C^{1}$, and the scalar functions $r_{i}:\mathbb{R}^{N_{\rm{state}}\times\mathcal{\tilde{P}}}\to\mathbb{R},~{}i=0,..,N_{t},$ belong to the class $C^{1}$ in the first variable. Then the first order necessary optimality conditions of the problem (7) are given by:

[TABLE]

where ${\mathbf{M}}_{i+1,i}^{*}$ is the Jacobian matrix of $\mathcal{M}_{i,i+1}$.

Proof 1

Using the Lagrange multiplier technique the constrained optimization problem (7) is replaced with the unconstrained optimization of the following Lagrangian function, $\mathcal{L}:\mathbb{R}^{N_{\rm{state}}}\to\mathbb{R}$

[TABLE]

where ${\mathbf{\lambda}}_{i}\in\mathbb{R}^{N_{\textnormal{state}}}$ is the Lagrange multipliers vector at observation time $t_{i}$.

An infinitesimal change in $\mathcal{L}$ due to an infinitesimal change $\delta{\mathbf{x}}_{0}$ in ${\mathbf{x}}_{0}$ is

[TABLE]

where $\delta{\mathbf{x}}_{i}=\frac{\partial{\mathbf{x}}_{i}}{\partial{\mathbf{x}}_{0}}\delta{\mathbf{x}}_{0}$ and ${\mathbf{M}}_{i,i+1}$ is the Jacobian matrix of $\mathcal{M}_{i,i+1}$ with respect to $\mathbf{x}_{i}$ for all time instances $t_{i}$, $i=0,\dots,N_{t}-1$, at ${\bf x}_{i}$,

[TABLE]

The corresponding adjoint operator is ${\bf M}_{i+1,i}^{*}\in\mathbb{R}^{N_{\textnormal{state}}\times N_{\textnormal{state}}}$ and satisfies

[TABLE]

where $\langle\cdot,\cdot\rangle_{\mathbb{R}^{N_{\textnormal{state}}}},$ is the corresponding Euclidian product. This means that the adjoint operator of ${\mathbf{M}}_{i,i+1}$ is nothing but its transpose. To underline that the corresponding adjoint model is running backward in time we denote it by ${\bf M}_{i+1,i}^{*}$.

After rearranging the above equation and using again the definition of the adjoint operators we obtain

[TABLE]

By setting to zero the perturbations with respect to $\delta{\mathbf{\lambda}}_{i}$ and $\delta\mathbf{x}_{i}$, $i=0,\ldots,N_{t}$, the first order necessary optimality conditions (8) are obtained.

While the gradient is equal to zero at the optimum, the formula (8c) can still be used to evaluate the gradient of $\mathcal{Q}$ away from the optimum

[TABLE]

Then for ${\bar{\mathbf{x}}_{0}}\in\mathbb{R}^{N_{\rm{state}}}$ in a vicinity of any $\mathbf{x}_{0}\in\mathbb{R}^{N_{\rm{state}}}$, the following first order approximation holds

[TABLE]

4.2 Error estimates

We now introduce the first a-posteriori error estimation result.

Theorem 2** (A-posteriori error estimation)**

Let $\mathbf{x}=\{\mathbf{x}_{0},~{}\mathbf{x}_{1},\ldots,~{}\mathbf{x}_{N_{t}}\}$ be the solution of the high-resolution (1) model, and $\widehat{\mathbf{x}}=\{\widehat{\mathbf{x}}_{0},\widehat{\mathbf{x}}_{1},\ldots,\widehat{\mathbf{x}}_{N_{t}}\}$ the projection of reduced order model solution(2) onto the full space. Moreover, let $\mathbf{x}^{i}=\{{\widehat{\mathbf{x}}_{i}},~{}\mathbf{x}^{i}_{i+1},\ldots,\mathbf{x}^{i}_{N_{t}}\}$ be the partial trajectories obtained via the full model (1) using as initial conditions the solution of reduced order model (2) at time $t_{i}$ projected onto the full space, i.e. ${\widehat{\mathbf{x}}_{i}}=U\,{\mathbf{\tilde{x}}_{i}}$ and

[TABLE]

The partial trajectory $\mathbf{x}^{i}$ contains only $N_{t}-i+1$ time steps.

Assume that the reduced order model solution $\widehat{\mathbf{x}}_{i}$ is in a neighborhood of the high-resolution model solution $\mathbf{x}_{i},~{}i=0,\ldots,N_{t}$ and if the assumptions of Theorem 1 hold, then the first order estimation of the error between the quantity of interest computed with the high-fidelity and reduced order models (6) is given by

[TABLE]

where the model residuals are

[TABLE]

and ${\widehat{\mathbf{\lambda}}}_{i},~{}i=0,\ldots,N_{t}-1$, is the solution of the full adjoint model (8b) linearized about the trajectory $\mathbf{x}^{i},~{}i=0,\ldots,N_{t}-1.$ At the final time step

[TABLE]

Proof 2

For each partial trajectory $\mathbf{x}^{i}$ we introduce the associated quantity of interest

[TABLE]

Since the trajectories are computed using the high-fidelity model, the gradient of each quantity of interest with respect to $\widehat{\mathbf{x}}_{i}$ is given by a relation analogous to (12). Since the reduced order model solution $\widehat{\mathbf{x}}_{i}$ lies in a neighborhood of the high-resolution model solution $\mathbf{x}_{i},~{}i=0,\ldots,N_{t}$, then from equation (13), we obtain the following first order approximations

[TABLE]

[TABLE]

where $\widehat{\lambda}_{i},~{}i=0,..,N_{t}-1$ is the solution of the adjoint model (8b) linearized at the trajectory $\mathbf{x}^{i},~{}i=0,\ldots,N_{t}-1.$

Since $r_{N_{t}}$ is continuously differentiable the following first order estimation holds

[TABLE]

According to (8b) we have $\widehat{\lambda}_{N_{t}}=-(\partial r_{N_{t}}/\partial\mathbf{x}_{N_{t}})^{T}(\widehat{\mathbf{x}}_{Nt},\mu)$ and adding all three equations (18) - (20), we obtain

[TABLE]

Figure 1 explains the procedure at an intuitive level. The discrepancy in the quantity of interest (5) computed with the high-resolution model solution (bottom trajectory) and projected reduced order model solution (top trajectory) can be estimated by adding all the intermediate dot products.

4.3 Efficient computation of the error estimates

The a-posteriori error estimation result developed in Theorem 2 requires one full, and several partial high-fidelity model runs along with their associated high-fidelity adjoint model runs. Here we propose efficient a-posteriori estimation results for a general explicit and implicit time integration schemes.

Explicit Euler scheme

The general model introduced in (1) can be described by

[TABLE]

where $h$ is the selected discrete time step.

Using the reduced order solution $\mathbf{\tilde{x}}_{i}$ of model (2) we perform one time step integration with both high-fidelity and reduced order models

[TABLE]

By multiplying (24) with $U_{\mu}$ from the left and subtracting the result from (23) we obtain

[TABLE]

where $\phi:\mathbb{R}^{N}_{\rm{state}}\times\mathbb{R}^{N}_{\rm{state}}\times\mathcal{\tilde{P}}\to\mathbb{R}^{N}_{\rm{state}}$ is the residual associated with the explicit full model

[TABLE]

Remark 1

If the projected reduced order model solutions $\widehat{\mathbf{x}}$ is accurate with respect to the high-fidelity model solution $\mathbf{x}$ then the partial trajectories $\mathbf{x}^{i}=\{{\widehat{\mathbf{x}}_{i}},~{}\mathbf{x}^{i}_{i+1},..$, $\mathbf{x}^{i}_{N_{t}}\}$, $i=0,..,N_{t}-1$ can be approximated by truncated trajectories obtained using one single high-fidelity model run $\{\mathbf{x}_{0},~{}\mathbf{x}_{1},..,~{}\mathbf{x}_{N_{t}}\}$. Then for estimating ${\widehat{\mathbf{\lambda}}}_{i}^{T},~{}i=0,\ldots,N_{t}$ in (15) we need only a single high-fidelity adjoint model run instead of several partial high-fidelity adjoint model trajectories required by Theorem 2.

Remark 2

Unlike the Galerkin POD residual, the DEIM based residual (26) is not orthogonal to the reduced manifold $U_{\mu}$. Thus we can compute the mismatch in (6) by employing a reduced order adjoint model solution

[TABLE]

By incorporating the high-fidelity adjoint snapshots for the construction of the reduced order basis [78] we can design an accurate reduced order adjoint model. The error estimate proposed in (27) requires only one reduced forward and one adjoint model runs as well as evaluating the residuals in (26). Here $\mathbf{\tilde{\lambda}}_{i}\in\mathbb{R}^{k}$ denotes the solution of the reduced adjoint model at time $t_{i}$

[TABLE]

Techniques for fast evaluations of the reduced Jacobians exist in the literature and a comparison between them is available in [76].

Remark 3

The use of a reduced adjoint model (28) allows the a-posteriori error estimation result in Theorem 2 to be exploited for different parametric configuration. The more accurate the global reduced forward and adjoint models are, the more precisely the a-posteriori estimate in Theorem 2 is.

Implicit Euler scheme

The model (1) has the following form

[TABLE]

where each component of $\mathbf{F}:\mathbb{R}^{N_{\rm{state}}\times\mathcal{\tilde{P}}}\to\mathbb{R}^{N_{\rm{state}}}$ is assumed to be of class $C^{1}$ in the first variable. We take one implicit Euler step with the full and reduced order models initialized at the reduced order solution:

[TABLE]

By subtracting (31) from (30) after projecting equation (31) to the full space we obtain

[TABLE]

The error (32) is approximated to first order by

[TABLE]

and then

[TABLE]

The high-fidelity adjoint model solution $\mathbf{\lambda}_{i}$ of the implicit Euler model associated with the quantity of interest defined in (5) is

[TABLE]

As stated in Remarks 1 and 2, we can assume accurate forward and reduced order model solutions;i.e., $\lambda_{i+1}\approx{\widehat{\mathbf{\lambda}}}_{i+1}$. Left-multiplying equation (34) by ${\widehat{\mathbf{\lambda}}}_{i+1}$ leads to

[TABLE]

where the second equality follows from the adjoint model (LABEL:eqn::adjoint_implicit_Euler).

The error in the quantity of interest (15) due to the usage of the reduced order model can be also estimated by

[TABLE]

where $\phi:\mathbb{R}^{N}_{\rm{state}}\times\mathbb{R}^{N}_{\rm{state}}\times\mathcal{\tilde{P}}\to\mathbb{R}^{N}_{\rm{state}}$ is now the residual associated with the implicit full model.

[TABLE]

As in the case of the explicit model, we can compute the estimated error in (37) using only one reduced forward and adjoint model runs and the evaluation of the high-fidelity model residuals (38). The statement in Remark 3 is valid here too.

5 Adaptive location of DEIM interpolation points

Here we describe a novel algorithm for selection of the DEIM points such that the accuracy of the quantity of interest $\mathcal{Q}$ (5) evaluated using a reduced order model is increased. In [63] the adaptivity mechanism changes the non-linear term reduced basis via rank-one updates. Here we do not change the basis, but rather we adaptively relocate the DEIM interpolation points.

Our adaptive strategy uses the a-posteriori error estimation result (15) and the fast computation techniques presented in Section 4.3. The individual contribution at each spatial location and time step to the error in the quantity of interest can be calculated by using the Hadamard product $\odot$ instead of the scalar products in (27) and (37). The Hadamard products are the dual weighted residuals.

For the explicit case the dual weighted residuals are defined as

[TABLE]

while for the implicit case they are defined as

[TABLE]

The high-fidelity explicit and implicit models residuals $\phi_{i}(\mu)$ are defined in (26) and (38).

Next the dual weighted residuals are collected into a matrix $Z$ and a singular vector decomposition is applied to extract the left singular vectors denoted by $W=\{w_{0},~{}w_{1},\ldots,w_{m}\}$. Now we are ready to introduce our adaptive strategy. In addition to the non-linear basis $V_{\mu}$, the dual residual basis $W$ is used as input.

The residual $r^{w}$ estimates the error of the dual weighted residuals representations in the non-linear term subspace $V_{\mu}$. In contrast with the traditional DEIM approach, we place the points where the highest value of the combined residual $\alpha|r^{v}|+(1-\alpha)|r^{w}|$ is found. This allows to generate a satisfactory globally accurate non-linear reduced order term, while enhancing the accuracy in the spatial locations with the higher contributions to (6). Here the parameter $\alpha$ is chosen heuristically, and finding an automated selection procedure is an open problem.

6 Numerical experiments

We evaluate the a-posteriori error estimation results and adaptive DEIM algorithm for two non-linear test problems, the one-dimensional Burgers and the two-dimensional Shallow Water Equations. Proper orthogonal decomposition and discrete empirical interpolation methods are applied along with the Galerkin projection to generate the associated. We validate the a-posteriori error formulas (15) and (37) and show that the adaptive DEIM strategy proposed in Section 5 is successfully in reducing the error (6). Our experiments are performed for the same parametric configurations for the Swallow Water Equations model. Additional efforts for the 1D-Burgers model confirm the Remark 3 statement, thus enabling the use of the a-posteriori error results for parametric configurations others than the one used to construct the reduced order bases for both state variables and non-linear terms.

6.1 One-dimensional Burgers model

Burgers’ equation [12, 13] provides a simplification of the equations of fluid dynamics by omitting the pressure terms. For a given viscosity coefficient $\mu$, the evolution of the fluid velocity $u$ is given by

[TABLE]

We assume Dirichlet homogeneous boundary conditions $u(0,t)=u(L,t)=0,~{}t\in(0,t_{\textnormal{f}}]$ and as initial conditions we use a seventh degree polynomial.

The discretization uses a uniform spatial mesh with $n$ grid points and $\Delta x=L/(n-1)$, and uniform temporal mesh with $N_{t}$ grid points and $\Delta t=t_{\textnormal{f}}/(N_{t}-1)$. The discrete solution vector is denoted by ${\boldsymbol{u}}(t_{N})\approx[u(x_{i},t_{N})]\in\mathbb{R}^{N_{\rm{state}}},~{}N=1,2,\ldots,N_{t}$ (the spatial dimension is $N_{\rm{state}}=n-2$ after eliminating the known boundaries). The semi-discrete version of the 1D Burgers model (41) reads

[TABLE]

where ${\bf u}^{\prime}$ denotes the time derivative of ${\bf u}$. $A_{x}$, $A_{xx}\in\mathbb{R}^{N_{\rm{state}}\times N_{\rm{state}}}$ are the central difference first-order and second-order space derivatives operators which also account for the boundary conditions, respectively.

The implicit Euler method is employed for time discretization and it is implemented in Matlab. The non-linear algebraic systems are solved using Newton-Raphson method and the maximum number of Newton iterations allowed each time step is set to $50$. The solution is considered accurate enough when the euclidian norm of the residual is less than $10^{-10}$.

6.2 Numerical experiments with the one-dimensional Burgers model

We propose the following configuration $L~{}=~{}1$, $~{}t_{f}~{}=~{}1$, $n~{}=~{}201$, $N_{\rm{state}}~{}=~{}199,$ $N_{t}~{}=~{}201$. The viscosity parameter is set initially to $\mu=0.1$, and the initial conditions are depicted in Figure 2a. The quantity of interest depends only on some particular components of the solution at the final time step (colored in green in Figure 2b) and it is defined below

[TABLE]

A number of $401$ snapshots are used to construct the reduced order basis for the state variable and include the solutions of both high-fidelity forward and adjoint models. For the advection non-linear term we applied the singular value decomposition and generate the reduced order basis required by the DEIM approximation using $201$ snapshots. The spectra of the snapshots matrices are illustrated in Figure 3, and for POD basis dimensions larger than $25$, the a-priori estimate suggests reduced order solutions errors smaller than $1e-5$ for the same parametric configuration. However this estimate does not include integration error and the overall error is usually underestimated.

Since the implicit Euler method was used to discretize the 1D-Burgers model, we will verify the theoretical a-posteriori error estimate result described in (37). First we compute the error induced by estimating the quantity of interest (43) using the POD/DEIM reduced order model and compare it against the dual weighted residuals sum. The formula in (37) requires only reduced order model runs and the results are depicted in Figure 4. By increasing the dimension of the reduced manifold, the error estimation gets closer to the true error as seen in Figure 4(a) and for POD basis larger than $12$ the discrepancies are lower than $1e-03$. The more accurate the underlying reduced order model the more precise the proposed a-posteriori error estimation result is. The sensitivity of the estimated error (37) with respect to the number of DEIM interpolation points is shown in 4(b) where dimension of POD basis is set to $15$. For more than $14$ DEIM points the mismatches between the true and estimated errors are smaller then $1e-4$.

Next we test the a-posteriori estimation result (37) for a different parametric configuration. As such, using bases computed for $\mu=0.1$, we aim to estimate the quantity of interest error for $\mu=0.07$. The approximation will be accurate as long the reduced order forward and adjoint models will be accurate with respect to the high fidelity models. First we set the number of DEIM points to $40$ and compute the estimates for various POD basis dimensions. A comparison against the true error is depicted in Figure 5(a). We notice larger discrepancies than in the case of using similar parametric configurations. For POD basis dimension larger than $17$, our estimation is very accurate.

By fixing the POD basis dimension to $15$ and varying the number of DEIM points, we set up another experiment and the results are shown in Figure 5(b). The a-posteriori estimation is accurate though is less precise than in the case of similar parametric configuration.

Next we will prove the usefulness of the adaptive DEIM algorithm 3 by enhancing the accuracy of the quantity of interest computed using the updated reduced order model. First we compute the dual weighted residuals (40) and perform a singular value decomposition applied to generate left singular vectors. The first $15$ of the singular vectors are stored in matrix $W$ and then used, together with the non-linear term basis $V_{\mu}$, to initiate Algorithm 3. The parameter $\alpha$ is set to $0.5$

The locations of the DEIM points computed using Algorithms 2 and 3 are shown in Figure 6. The adaptive algorithm places more points in the spatial domain $[0.05,0.1]$ of the quantity of interest marked with red in Figure 6. The locations of the adaptive DEIM points change the matrix $P$ in DEIM approximation of the non-linear term (4).

Figure 7(a) illustrates the comparison between the 1D-Burgers advection term errors computed at time $t_{2}$ using the standard vs. adaptive DEIM methods. The number of points is varied and the results show that standard DEIM approximation is more accurate over the global spatial domain $[0,1]$. This is expected since the adaptive DEIM method is tailored to improve the accuracy of the quantity of interest and not the global accuracy. The condition number of the matrix $P^{T}V_{\mu}$ (4) is increased in the case of adaptive DEIM approximaton, contributing to the global loss of accuracy. This result is shown in Figure 7(b).

The adaptive DEIM approximation of the non-linear term leads to a different reduced order solution. The quantity of interest is computed using POD/DEIM standard and adaptive models and its errors are presented in Figure 8. The adaptive model leads to more accurate quantities of interest confirming the usefulness of our a-posteriori error estimation results and adaptive DEIM algorithm.

6.3 Shallow Water Equations (SWE) model

The SWE is a popular simple model for meteorological and oceanographic problems. SWE can be used to model Rossby and Kelvin waves in the atmosphere, rivers, lakes and oceans as well as gravity waves in a smaller domain. The alternating direction fully implicit (ADI) scheme [35] considered in this paper is first order in both time and space and it is stable for large CFL condition numbers. It has been proven that the method is unconditionally stable for the linearized version of the SWE model. Other research work on this topic include efforts of Fairweather and Navon [30] and Navon and Villiers [54].

We solve the SWE model using the $\beta$-plane approximation on a rectangular domain [35]

[TABLE]

where $w=(u,v,\phi)^{T}$ is a vector function, and $u,v$ are the velocity components in the $x$ and $y$ directions, respectively. The geopotential $\phi=2\sqrt{gh}$ is computed by multiplying the depth of the fluid $h$ and the acceleration due to gravity $g$.

The matrices $A$, $B$ and $C$ are

[TABLE]

where $f$ is the Coriolis term

[TABLE]

with $\hat{f}$ and $\beta$ constants.

We assume periodic solutions in the $x$ direction for all three state variables. In the $y$ direction

[TABLE]

and Neumann boundary conditions are used for $u$ and $\phi$. The initial condition is $w(x,y,0)=\psi(x,y),~{}\psi:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R},~{}(x,y)\in[0,L]\times[0,D]$.

The space discretization is performed on a uniform mesh with $n=N_{x}\cdot N_{y}$ equidistant points on $[0,L]\times[0,D]$, with $\Delta x=L/(N_{x}-1),~{}\Delta y=D/(N_{y}-1)$. We also discretize the time interval $[0,t_{\rm f}]$ using $N_{t}$ equally distributed points and $\Delta t=t_{\rm f}/(N_{t}-1)$. The discrete solution vector is:

[TABLE]

After spatial discretization of (44) the semi-discrete SWE equations are:

[TABLE]

where $\odot$ is the component-wise multiplication operator, ${\bf u}^{\prime}$, ${\bf v}^{\prime}$, ${\bm{{\phi}}}^{\prime}$ denote semi-discrete time derivatives, and ${\bf F}={\bf f}\cdot\mathbf{1}_{N_{x}}^{T}$ stores Coriolis components ${\bf f}=[f(y_{j})]_{j=1,2,..,N_{y}}$. The non-linear terms involving derivatives in $x$ and $y$ directions, respectively, are defined as follows:

[TABLE]

Here $A_{x},A_{y}\in\mathbb{R}^{n\times n}$ are constant coefficient matrices for discrete first-order and second-order differential operators which take into account the boundary conditions.

The numerical scheme is implemented in Fortran and uses a sparse matrix environment. For operations with sparse matrices we employ the SPARSEKIT library [69], and the sparse linear systems obtained during the quasi-Newton iterations are solved using MGMRES library [9, 44, 70]. Here we do not decouple the model equations like in Stefanescu and Navon [75] where the Jacobian is either block cyclic tridiagonal or block tridiagonal.

6.4 Numerical experiments with the shallow water model

We computed the initial conditions from the initial height condition No. 1 of Grammeltvedt [32] i.e.

[TABLE]

The initial velocity fields are derived from the initial height field using the geostrophic relationship

[TABLE]

We use the following constants $L=6000km,~{}D=4400km,~{}\hat{f}=10^{-4}s^{-1}$, $\beta=1.5\cdot 10^{-11}s^{-1}m^{-1},~{}g=10ms^{-2},~{}H_{0}=2000m,~{}H_{1}=220m,~{}H_{2}=133m.$

The domain is discretized using a mesh of $31\times 17=527$ points, with $\Delta x=200$km and $\Delta y=275$km. We select the integration time window to be $24$h and we use $181$ time steps corresponding to $\Delta t=480$s.

The considered quantity of interest depends on some particular components of the geopotential $\phi$ at the final time step and it is defined below

[TABLE]

To generate the reduced order forward and adjoint SWE models, we collect $180$ snapshots from each high-fidelity model representing the state solutions $u$, $v$ and $\phi$ and their corresponding adjoint variables at every time step. For each state variable, a snapshot matrix with $360$ columns containing both forward and adjoint solutions is formed. Three POD bases are then derived following singular value decompositions. Moreover, snapshots of the six forward non-linear terms were also collected and $6$ POD bases required by the DEIM approximation were constructed. Finally, we apply a Galerkin projection to the discrete SWE model and the forward POD/DEIM reduced order model is derived. More details regarding the construction of the SWE reduced order model can be found in [77]. Applying chain derivatives, the reduced adjoint model is obtained following the ARRA method [78].

Next we asses the a-posteriori error estimation result (15) using numerical experiments. The quantity of interest (46) is computed using the forward high-fidelity and reduced order models and the true error is obtained. Using the projected solution of the reduced of adjoint model the estimated error of the quantity of interest is computed. Comparative results are shown in Figure 9. In panel (a), we set the number of DEIM points and test various POD bases dimensions. The estimated errors are very accurate even for smaller bases dimensions. Panel (b) describes a different experiment where the dimensions of POD bases are fixed to $30$ and the number of DEIM points is varied. Again, the estimated errors shows good agreement with the true.

In the sequel, we make use of the dual weighted residuals employed by the a-posteriori error estimation result (15) to generate adaptive DEIM points and decrease the quantity of interest error computed using the reduced order model. The dual weighted residuals are first computed replacing the dot product with the component-wise multiplication in the right-hand side part of the expression (15); i.e., ${\widehat{\mathbf{\lambda}}}_{i}^{T}\cdot\Delta\mathbf{x}_{i},~{}i=0,\ldots,N_{t}$. These residuals are then used together with the non-linear bases to generate the adaptive DEIM points. For example, the continuity equation dual weighted residuals are employed together with the non-linear basis of the non-linear term $F_{31}$ and the resulted adaptive DEIM points are depicted in Figure 10. We notice that more points are placed in the spatial domain of the quantity of interest delimited by the red rectangle in comparison with the standard DEIM points. Parameter $\alpha$ required by Algorithm 3 was set to $\frac{1}{2}$.

The adaptive DEIM algorithm changes the DEIM approximation of the non-linear terms since matrix $P$ (4) is modified. This leads to modified accuracy of the non-linear terms approximations over all the spatial points. The errors associated with various DEIM approximations of the non-linear terms $F_{31}$ at time step $t_{2}$ are shown in Figure 11 (a). The adaptive DEIM approximation is less accurate than the standard DEIM over the entire spatial domain. This is explained by larger condition numbers of $P^{T}V_{\mu}$ (4) generated by the adaptive DEIM algorithm as noticed in Figure 11 (b).

Next we compare the quantity of interest errors obtained using POD/DEIM and POD/adaptive DEIM reduced order models. Among $11$ values of

[TABLE]

we show only the resulted most accurate estimates in Figure 12 for two POD bases dimensions. The labels along the adaptive DEIM trajectories describe the values of $\alpha$ that generated the most accurate quantities of interest. In general, it is more difficult to enhance the quantity of interest accuracy once the reduced order state solution is very precise. This is obvious by comparing the results presented in Figure 12. Larger POD basis dimensions usually leads to more accurate ROM solutions.

7 Conclusions

In this paper we proposed a-posteriori error estimation results to estimate the errors in quantities of interests calculated using POD/DEIM reduced order models. Later, using the dual weighted residuals, we introduced an adaptive DEIM algorithm to improve the accuracy of these quantities of interests.

The first a-posteriori error estimation result sums up the first-order expansions of the quantity of interest’s components computed as the product between the gradient of each component and high-fidelity model residuals. The gradient of each component is obtained as the solution of a high-fidelity adjoint model. Efficient versions of this a-posteriori error estimation result are obtained for the explicit and implicit Euler schemes by assuming a good level of accuracy for the surrogate models solutions. Only one reduced forward and adjoint model runs and the evaluation of the high-fidelity model residuals are required for estimating the quantity of interest errors for both schemes. Reduced order models generated using a single parametric configuration were employed to estimate the quantity of interest for a different parametric configuration in the case of 1D-Burgers model. In the case of the SWE model, we tested only the same parametric configuration case. Several POD basis dimensions and numbers of DEIM points were tested and the estimated error accuracy increased with more accurate reduced order model solutions.

Next we designed an adaptive DEIM algorithm to increase the precision of the quantity of interest computed using the forward reduced order model. The singular vectors of the dual weighted residuals together with the non-linear term basis are used to generate the new DEIM points. This changes the DEIM approximation of the non-linear terms and subsequently the associated reduced order model solutions. Numerical results using the new reduced order models revealed that the accuracy of the quantities of interest was enhanced for both 1D-Burgers and SWE models.

In the future we plan to apply regression machine learning models to predict the coefficient $\alpha$ used to combine the dual weighted residuals singular vectors and non-linear basis within the adaptive DEIM algorithm. In combination with MP-LROM models [53], the a-posteriori error estimation result and the adaptive DEIM algorithm will provide accurate outcomes along the entire parametric space.

Acknowledgment

The work of Dr. Răzvan Stefanescu and Prof. Adrian Sandu was supported by the NSF CCF–1218454, AFOSR FA9550–12–1–0293–DEF, AFOSR 12-2640-06, and by the Computational Science Laboratory at Virginia Tech. Răzvan Ştefănescu would like to thank Vishwas Rao for useful conversations about error estimation results.