\RequirePackage[l2tabu, orthodox]{nag}  % checks for obsolete LaTeX packages and outdated commands - See more at: http://www.howtotex.com/packages/9-essential-latex-packages-everyone-should-use/#sthash.QU9O2S6g.dpuf

\documentclass[10pt]{article}

% This first part of the file is called the PREAMBLE. It includes
% customizations and command definitions. The preamble is everything
% between \documentclass and \begin{document}.

\usepackage[margin=1.2in]{geometry}  % set the margins to 1in on all sides
\usepackage{graphicx}       % to include figures
\usepackage{amsmath,amssymb}% great math stuff
\usepackage{amsfonts}       % for blackboard bold, etc
\usepackage{amsthm}         % better theorem environments
% \usepackage{times}          % to change default font to Times

\usepackage{microtype}      % must be after font: improves the spacing between words and letters

\usepackage{varwidth}       % to include a box frame outside algorithm
\usepackage{algorithm}
\usepackage{algpseudocode}  % pseudocode / algorithm
\usepackage{caption}        % provide caption for algorithm
\usepackage{url}            % cite url in reference

\usepackage{subfigure}      % provide subfigures
%\usepackage{mystyle}       % include mystyle.sty file
\usepackage{authblk}        % for many different authors and affiliations

\usepackage{epstopdf}       % import figures using pdflatex

\usepackage{color, soul}    % to highlight text
\usepackage{enumerate}      % different format for enumeration

\usepackage{hyperref}       % internal hyperlink to contents
\hypersetup{
    colorlinks=true,   % boxed links
%    colorlinks=false,   % colored links
    citecolor=blue,     % citation links (bibliography)
    filecolor=blue,     % color of file links
    linkcolor=red,    % internal links (sections, pages, etc.)
    urlcolor=blue}      % color of URL links (mail, web)

%\everymath{\displaystyle}

\usepackage{cleveref}   % must be last


\newcommand{\green}[1]{\textcolor{green}{#1}}   % short for green text

\renewcommand{\thefootnote}{\fnsymbol{footnote}}    % footnote symbols

\newcommand{\pr} {{\rm pr}}     % roman font of the word "pr"
\newcommand{\du} {{\rm du}}     % roman font of the word "du"
\newcommand{\del} {\partial}    % short of \partial
\newcommand{\til} {\tilde}      % short of \partial
%\newcommand{\st}[2]{{\bf #1_{\it #2}^{st}}}   % short for space-time



% various theorems, numbered by section

\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{conj}[thm]{Conjecture}
\newtheorem{remark}[thm]{Remark}

\DeclareMathOperator{\id}{id}

%\newcommand{\bd}[1]{\mathbf{#1}}  % for bolding symbols
%\newcommand{\RR}{\mathbb{R}}      % for Real numbers
%\newcommand{\ZZ}{\mathbb{Z}}      % for Integers
%\newcommand{\col}[1]{\left[\begin{matrix} #1 \end{matrix} \right]}
%\newcommand{\comb}[2]{\binom{#1^2 + #2^2}{#1+#2}}



\makeatletter
\def\BState{\State\hskip-\ALG@thistlm}
\makeatother





\begin{document}

%%% to include all references in the bib file appear (without cite to them)
%\nocite{*}


\title{\Large \bf An efficient goal-oriented sampling strategy using reduced basis method for parametrized elastodynamic problems}

\author[1]{K. C. Hoang}
\author[1]{P. Kerfriden}
\author[2]{B. C. Khoo}
\author[1,3]{S. P. A. Bordas}

\affil[1]{School of Engineering, Cardiff University, The Parade, CF24 3AA Cardiff, United Kingdom}
\affil[2]{Singapore MIT-Alliance, National University of Singapore, Singapore}
\affil[3]{Facult\'{e} des Sciences, de la Technologie et de la Communication, University of Luxembourg, 6 rue Richard Coudenhove-Kalergi, L-1359 Luxembourg}
\affil[ ]{ \textit{Emails: \{hoangkhacchi, pierre.kerfriden, stephane.bordas\}@gmail.com, mpekbc@nus.edu.sg} }


%\author{K. C. Hoang, P. Kerfriden and S. P. A. Bordas\thanks{Email: \{hoangkhacchi, pierre.kerfriden, stephane.bordas\}@gmail.com} \vspace{9pt}\\
%Cardiff University, School of Engineering, \\
%The Parade, CF24 3AA Cardiff, United Kingdom }

\date{}

\maketitle

\begin{abstract}

\noindent In this paper, we study the class of linear elastodynamic problems with affine parameter dependence using a goal-oriented approach by finite element (FE) and reduced basis (RB) methods. The main contribution of this paper is the ``goal-oriented'' proper orthogonal decomposition (POD)--Greedy sampling strategy within the RB approximation context. The proposed sampling strategy looks for the parameter points such that the output error approximation will be minimized by Greedy iterations. In estimating such output error approximation, the standard POD--Greedy algorithm is invoked to provide enriched RB approximations for the FE outputs. We propose a so-called ``cross-validation'' process to choose adaptively the dimension of the enriched RB space corresponding with the dimension of the RB space under consideration. Numerical results show that the new goal-oriented POD--Greedy sampling procedure with the cross-validation process improves significantly the space-time output computations in comparison with the ones computed by the standard POD--Greedy algorithm. The method is thus ideally suited for repeated, rapid and reliable evaluations of input-output relationships in the space-time setting.


% \noindent In this paper, we study the class of linear elastodynamic problems with affine parameter dependence using a goal-oriented approach by finite element (FE) and reduced basis (RB) methods. The main contribution of this paper is the ``goal-oriented'' proper orthogonal decomposition (POD)--Greedy sampling strategy within the RB approximation context. The proposed sampling strategy looks for the parameter points such that the output error approximation will be minimized by Greedy iterations. This is different from and improves further the standard POD--Greedy algorithm already developed, where the solution error (or, solution error indicator) will be minimized. In estimating the output error approximation within this proposed algorithm, the standard POD--Greedy algorithm is used to provide adaptively good approximations for the FE output computations. Numerical results show that the new goal-oriented POD--Greedy sampling procedure improves significantly the space-time output computations in comparison with the ones computed by the standard POD--Greedy algorithm. The method is thus ideally suited for repeated, rapid and reliable evaluations of input-output relationships in the space-time setting.


\bigskip


\noindent \textbf{Keywords}: reduced basis method; goal-oriented POD--Greedy algorithm; wave equation; goal-oriented asymptotic error; cross-validation


% modify the below abstract on 15/02/2014

%\noindent In this paper, we study the class of linear elastodynamic problems with affine parameter dependence using a goal-oriented approach by finite element (FE) and reduced basis (RB) methods. The main contribution of this paper is the ``goal-oriented'' proper orthogonal decomposition (POD)--Greedy sampling strategy within the RB approximation context. First, we introduce the RB approximation method: Galerkin projection onto a space $Y_N$ spanned by solutions of the governing partial differential equation (PDE) at $N$ selected points in parameter space. This set of $N$ parameter points is constructed by the standard POD--Greedy sampling strategy already developed. Second, based on the affine parameter dependence, we deploy the offline-online computational procedures: in the offline stage, we generate the RB space; in the online stage, given a new parameter value, we calculate rapidly and accurately the space-time RB quantity of interest and its associated asymptotic error. The proposed goal-oriented POD--Greedy sampling strategy can now be implemented and will look for the parameter points such that it minimizes this (asymptotic) output error rather than the solution error (or, error indicator which is the dual norm of residual) as in the standard POD--Greedy strategy. Numerical results show that the new goal-oriented POD--Greedy sampling procedure improves significantly the accuracy of the space-time output computation in comparison with the standard POD--Greedy one. The method is thus ideally suited for repeated, rapid and reliable evaluation of input-output relationships within the space-time setting.




% the following abstract is for the paper which has both forward and inverse analysis

%   This paper improves significantly the method developed earlier which is based on the reduced basis method and the Levenberg-Marquardt-Fletcher algorithm to identify the ``unknown'' material properties: Young's modulus and stiffness-proportional Rayleigh damping coefficient of the interfacial tissue between a dental implant and the surrounding bones. In the forward problem, a finite element approximation for a three-dimensional dental implant-bone model is first built. We then propose a new ``goal-oriented'' Proper Orthogonal Decomposition (POD)--Greedy algorithm in combination with the Galerkin projection to establish the reduced basis approximations to enable extremely fast and reliable computation of displacement responses for a range of material properties. This new algorithm makes use of the information of dual or adjoint problem similarly to the dual-weighted residual approach. Numerical results shows that the new ``goal-oriented'' POD--Greedy sampling procedure improves considerably the accuracy of the output displacement responses in comparison with the standard POD--Greedy one. In the inverse analysis, we propose combining the Genetic Algorithm (GA) algorithm with the Levenberg-Marquardt-Fletcher one as an efficient inverse solver. This inverse solver will use the reduced basis approximations as a fast forward solver to enable rapid identification of the unknown material properties. Numerical results show that the approach is very efficient and robust against highly contaminated noises.
%  \bigskip
%   \noindent \textbf{Keywords}: second-order hyperbolic partial differential equations; reduced basis method; goal-oriented POD--Greedy algorithm; dual problem; inverse analysis; Levenberg--Marquardt algorithm; material characterization


\end{abstract}




\section{Introduction}\label{sec_Intro}

The design, optimization and control procedures of engineering problems often require several forms of performance measures or outputs -- such as displacements, heat fluxes or flowrates. Generally, these outputs are functions of field variables such as displacements, temperature or velocities which are usually governed by a partial differential equation (PDE). The parameter or input will frequently define a particular configuration of the model problem. Therefore, the relevant system behavior will be described by an implicit input-output relationship; where its computation requires the solution of the underlying parameter-PDE (or $\mu$PDE). We pursue model order reduction (MOR) methods (i.e., snapshots-POD \cite{Sirovich1987, Kerfriden2010, Kerfriden2011a, Kerfriden2011b, Kerfriden2012, kerfriden2014certification} and RB \cite{Rozza2008, Hoang2013}) which permits the \textit{efficient} and \textit{reliable} evaluation of this PDE-induced input-output relationship in \textit{many query} and \textit{real-time} contexts.

%the RB method \cite{Rozza2008, Hoang2013} which permits the \textit{efficient} and \textit{reliable} evaluation of this PDE-induced input-output relationship in \textit{many query} and \textit{real-time} contexts.


The RB method was first introduced in the late 1970s for nonlinear analysis of structures and has been investigated and developed more broadly \cite{Noor1980}. Recently, the RB method was well developed for various kinds and classes of parametrized PDEs such as: the eigenvalue problems \cite{Machiels2000}, the coercive/non-coercive affine/non-affine linear/nonlinear elliptic PDEs \cite{Rozza2008, Veroy2003}, the coercive/non-coercive affine/non-affine linear/nonlinear parabolic PDEs \cite{Grepl2005, Grepl2007}, the coercive affine linear hyperbolic PDEs \cite{Tan2007, Hoang2013}, and several highly nonlinear problems such as Burger's equation \cite{Nguyen2009, jung2009reduced} and Boussinesq equation \cite{Knezevic2011}. For the linear wave equation, the RB method and associated \textit{a posteriori} error estimation was developed successfully with some levels \cite{Tan2007, Huynh2011b, HoangPhDThesis2012}; however, in the RB context none of these works have focused on constructing optimally goal-oriented RB basis functions.


Goal-oriented error estimates in the context of FE analysis have been investigated deeply and widely \cite{Gratsch2005, Ladeveze2012, Ruter2006, Gonzalez2012, Gonzalez2013, Larsson2001, bordas2007derivative, duflot2008posteriori} (we only cite some typical works as there many on this topic). For the wave equation, the most well-known method is the dual-weighted residual (DWR) one which was proposed by Rannacher and co-workers \cite{Bangerth1999, Bangerth2001, Becker2001, Bangerth2010}. In those works, the authors quantified the \textit{a posteriori} error of the interest output in order to finer locally the finite element mesh in an adaptive manner. The final goal is to minimize computational efforts and maximize the accuracy of the interest output in an adaptive and controllable manner. In particular, the DWR method makes use of an auxiliary dual (or sensitivity) equation to derive an \textit{a posteriori} error expression for the interest output from the primal residual and the dual solution of that dual equation in space-time setting.


Goal-oriented error estimates in the context of MOR is currently an active research topic and has been investigated by several authors. In this regard, the construction procedure of these goal-oriented MOR basis functions is the key issue. For instance, Liu et al. \cite{Liu2008} used a Greedy algorithm to construct goal-oriented RB basis functions based on asymptotic output errors \cite{Prudhomme2002}; the surrogate RB model was then used in an inverse analysis. Chen et al. \cite{chen2013accurate} developed hybrid and goal-oriented Greedy sampling algorithms to compute failure probability for PDEs with random input data. In another work, Urban et al. \cite{urban2014greedy} developed a goal-oriented sampling strategy which consists of solving an optimization problem and a goal-oriented Greedy sampling to find the optimal parameter samples to best approximate interest outputs. We note that the representative works mentioned above are for stationary and steady problems.



For dynamic problems, goal-oriented sampling strategy for MOR were also addressed by several authors. Meyer and Matthies \cite{Meyer2003} combined the DWR with the snapshots-POD method to solve a nonlinear dynamics problem. They quantified the \textit{a posteriori} error approximation from the contributions of all POD snapshots; then the MOR basis functions are built (based on these POD snapshots) by keeping only the snapshots that caused large errors and removing all the ones which caused smaller errors. In another well-known approach by Bui \cite{Bui2007} and Willcox \cite{Willcox2005}, they solved a PDE-constrained optimization problem to find the optimally goal-oriented set of basis functions. In this way, the optimally goal-oriented basis functions are found such that they minimize the true output errors (with appropriate regularization techniques) and subject to equilibrium PDE-constraints \cite{Nocedal1999, Heinkenschloss2008}.


In general, those two above approaches are optimal. However, their computational cost are very expensive since one has to compute all the FEM solutions/outputs over the entire parameter domain (all POD snapshots -- for the former approach), and in every iteration within optimization solvers (for the latter approach); and hence, it would limit the number of input parameters in comparison with the RB approach.


%It is thus very natural to combine goal-oriented estimation with RB approximations \cite{Grepl2005, Meyer2003} to improve the accuracy of the interest output within RB context. Bringing the ``DWR idea'' to the RB context, we can also obtain an error expression which relates the output error with the primal residual of the dual solution in space-time setting. Namely, the error between the (space-time) FE and RB outputs will be exactly equal to the primal residual of the FE dual solution \cite{Bangerth2010}. However, to make this error expression truly \textit{a posteriori}, ones prefer to approximate the FE dual solution by some RB ones rather than using it directly within the expression \cite{Grepl2005, Meyer2003}. Hence, the error expression now becomes an \textit{a posteriori} error approximation since the equal sign (in the expression) now becomes the approximation one. This approach is generally good; but it suffers from two main drawbacks: first, it is only an error approximation which is not a true upper error bound and second, underlying concepts and computational effort are high/complicated since we need to solve additionally a dual problem and all related issues. We also implemented this approach recently in the work \cite{Hoang2013acme}.


%Here, we also give a brief review on goal-oriented estimations within reduced order modeling context using the snapshots-POD method \cite{Sirovich1987, Kerfriden2010, Kerfriden2011a, Kerfriden2011b, Kerfriden2012}. Goal-oriented estimations using the snapshots-POD method were implemented in \cite{Willcox2005, Bui2007}. In this approach, the authors solved an PDE-constrained optimization problem to find the optimal set of basis functions. In particular, the optimal basis functions are found such that they minimize the true output errors (with appropriate regularization techniques) and subject to equilibrium PDE-constraints \cite{Nocedal1999, Heinkenschloss2008}. This approach is optimal, however, it is very expensive since ones have to compute all the FEM solutions/outputs in every iteration within optimization solvers; and hence, it would limit the number of input parameters in comparison with the RB approach.


%After implementing the standard POD--Greedy sampling procedure, the set of basis functions (which is optimal in the sense of energy norm of the solution error) and the asymptotic output errors will be available. We propose to implement additionally the new ``goal-oriented'' POD--Greedy sampling procedure, which will pick up the parameter points such that the asymptotic output errors are minimized.
%In order to compute this output error approximation, the standard POD--Greedy algorithm will be invoked to provide a good approximation for the FE output computation.



In this work, we aim to build an optimally goal-oriented set of MOR basis functions without computing and storing all the POD snapshots as the aforementioned approaches. The best way to do that is by using the RB method with Greedy sampling strategy (see, for instance \cite{Veroy2003, Rozza2008, Hoang2013}). Thanks to the Greedy iterations, the proposed algorithm now looks for the parameter points such that the output error approximation will be minimized. For the linear wave equation, this idea is novel and further develops the idea of the standard POD--Greedy sampling procedure currently used \cite{Haasdonk2008, Knezevic2011, Hoang2013}, where the algorithm will pick up optimally all parameter points such that the error (or error indicator) of the field variable is minimized. By this way, we expect to improve significantly the accuracy of the RB output functional computations; but consequently, we might lose the rapid convergence rate of the field variable as in the standard POD--Greedy algorithm. In fact, as we can see later in the numerical results section, the convergent rate of the field variable by the two algorithms are quite similar; while the convergent rate of the output by the goal-oriented POD--Greedy algorithm is faster than that of the standard POD--Greedy one\footnote{In subsequent sections, for simplicity we shall call the ``standard algorithm'' to mention the standard POD--Greedy algorithm \cite{Haasdonk2008}, and the ``goal-oriented algorithm'' to mention the proposed goal-oriented POD--Greedy algorithm, respectively.}.


%We also note that in implementing the proposed algorithm, the standard algorithm will be required since it provides good approximations for the FE output computations.


In particular, the output error approximation used in this work is a kind of the asymptotic output error (such as in \cite{Liu2008, Prudhomme2002}) where the FE outputs will be approximated by the enriched RB outputs. The standard POD--Greedy algorithm will be invoked to compute such enriched RB outputs. Heuristically, the dimension of the enriched RB space will be usually set to two times larger than that of the RB space under consideration (see \cite{Liu2008, Prudhomme2002}). In this work, however, we devise a simple yet efficient algorithm called ``cross-validation'' process to find out adaptively the enriched RB dimension corresponding with each RB dimension under consideration. That process is performed within the offline stage of the proposed goal-oriented algorithm. We also note that this output error approximation will be used as both offline and online error indicators in the offline and online computational stages, respectively.


%These enriched RB outputs are computed by using standard POD--Greedy algorithm.
%
%As described in the works \cite{Liu2008, Prudhomme2002}, the dimension this enriched RB space will be set heuristically to two times larger than that of the RB space under consideration


%Specifically, in estimating the output error approximation, the FE outputs will be approximated by the RB outputs which are computed by using standard POD--Greedy algorithm. Heuristically, the dimension of this (enriched) RB space will be set to two times larger than that of the RB space under consideration (such as in \cite{Liu2008, Prudhomme2002} mentioned above). In this work, we devise a simple yet efficient algorithm called ``cross-validation'' process to find out adaptively this enriched RB dimension corresponding with each RB dimension under consideration.





The potential context for the proposed goal-oriented algorithm is described as follows. Suppose that one considers the linear parametrized wave equation with several different quantities of interest, and one wants to estimate the RB approximations of these quantities of interest. Note that these quantities of interest (or some of them) might be \textit{a priori} unknown, i.e., they may exist at the time of consideration, or they can appear afterwards depending on one's needs. Clearly, the standard POD--Greedy algorithm is not sufficient for this situation as it only provides the best approximations for the solution (or field variable) -- and not for these quantities of interest. Our proposal is as follows. The standard algorithm is implemented \textit{first} and only \textit{once} to create standard RB spaces\footnote{We will see later that these standard RB spaces have quite high dimensions (higher than that of the goal-oriented RB spaces) but are still very small compared to the FE space dimension.} to be used in the output error estimation of the proposed goal-oriented algorithm afterwards. Then, for each particular quantity of interest, the proposed goal-oriented algorithm will be performed once to build goal-oriented RB spaces corresponding to that quantity. By this way, the goal-oriented RB spaces are optimal for the quantity of interest under consideration (or in other words, best approximate this quantity of interest); and thus are much better than the standard RB spaces created by the standard algorithm.

%(And hence, the goal-oriented algorithm can be performed \textit{many times} depending on the particular quantities of interest in consideration.)





%we devise a simpler approach that addresses the two aforementioned drawbacks. First, we use asymptotic output errors which are also not truly upper error bounds but asymptotically converged to the true error; their performances are well comparable with the \textit{a posteriori} DWR error expression. Second, the approach completely do not use any information of the dual problem and hence simpler and easier to implement. Based on the standard POD--Greedy sampling procedure and the asymptotic output errors, we propose the new ``goal-oriented'' POD--Greedy sampling procedure, which will pick up the parameter points such that the error (or error approximation) of the output functional is minimized.





The paper is organized as follows. In Section \ref{sec_ProbStat}, we introduce necessary definitions, concepts and notations and then state the problem using a semidiscrete approach: fully discretizing in space using Galerkin FEM and marching in time using Newmark's trapezoidal rule. In Section \ref{sec_RB}, we describe various topics related to the RB methodology: approximation, the standard versus goal-oriented algorithms, error estimations and offline-online computational procedure. In Section \ref{sec_NumEx}, we verify the performance of the proposed algorithm by investigating numerically two problems: a 2D linear elastodynamic problem and a 3D dental implant simulation problem \cite{Hoang2013}. Finally, we provide some concluding remarks in Section \ref{sec_Conclusion}.






\section{Problem statement}\label{sec_ProbStat}


\subsection{Abstract formulation}\label{subsec_ProbStat_AbstForm}

We consider a spatial domain $\Omega \in \mathbb{R}^d$ ($d=1,2,3$) with Lipschitz continuous boundary $\partial \Omega$. We denote the Dirichlet portion of the boundary by $\Gamma_{D,i}, 1 \le i \le d$. We then introduce the Hilbert spaces

\begin{subequations}\label{eq:ProbStat_one}
\begin{equation}
Y^e = \{ v \equiv (v_1,\ldots,v_d) \in (H^1(\Omega))^d \; | \; v_i=0 \quad {\rm on} \quad \Gamma_{D,i}, i=1,\ldots,d \}, \label{eq:ProbStat_one-a}
\end{equation}
\begin{equation}
X^e = (L^2(\Omega))^d. \label{eq:ProbStat_one-b}
\end{equation}
\end{subequations}


Here, $H^1(\Omega) = \{ v \in L^2(\Omega) \; | \; \nabla v \in (L^2(\Omega))^d \}$ where $L^2(\Omega)$ is the space of square-integrable functions over $\Omega$. We equip our spaces with inner products and associated norms $(\cdot,\cdot)_{Y^e}$ $((\cdot,\cdot)_{X^e})$ and $\|\cdot\|_{Y^e}=\sqrt{(\cdot,\cdot)_{Y^e}}$ $(\|\cdot\|_{X^e}=\sqrt{(\cdot,\cdot)_{X^e}})$, respectively; a typical choice is

\begin{subequations}\label{eq:ProbStat_two}
\begin{equation}
(w,v)_{Y^e} = \sum_{i,j} \int_{\Omega} \frac{\partial w_i}{\partial x_j} \frac{\partial v_i}{\partial x_j} + w_i v_i, \label{eq:ProbStat_two-a}
\end{equation}
\begin{equation}
(w,v)_{X^e} = \sum_{i} \int_{\Omega} w_i v_i, \label{eq:ProbStat_two-b}
\end{equation}
\end{subequations}

where the summation over spatial dimensions $1 \le i,j \le d$ is assumed throughout this paper.




We define an input parameter set $\mathcal{D} \in \mathbb{R}^P$, a typical point in which shall be denoted $\mu \equiv (\mu_1,\ldots,\mu_P)$. We then define the parametrized bilinear forms $a$ in $Y^e$, $a:Y^e \times Y^e \times \mathcal{D} \rightarrow \mathbb{R}$; $m,c,f,\ell$ are parametrized continuous bilinear and linear forms in $X^e$, $m: X^e \times X^e \times \mathcal{D} \rightarrow \mathbb{R}$, $c:X^e \times X^e \times \mathcal{D} \rightarrow \mathbb{R}$, $f: X^e \times \mathcal{D} \rightarrow \mathbb{R}$ and $\ell: X^e \rightarrow \mathbb{R}$.



The ``exact'' continuous problem is stated as follows: given a parameter $\mu \in \mathcal{D} \subset \mathbb{R}^{P}$, the field variable $u^e(x,t;\mu) \in Y^e$ satisfies the weak form of the $\mu$-parametrized hyperbolic PDE (assume Rayleigh damping)

\begin{multline}\label{eq:ProbStat_three}
m \left( \frac{\partial^2 u^e(x,t;\mu)}{\partial t^2},v;\mu \right) + c \left( \frac{\partial u^e(x,t;\mu)}{\partial t},v;\mu \right) + a \left( u^e(x,t;\mu),v;\mu \right) = g(t) f(v;\mu), \\
\forall v \in Y^e, \; t \in [0,T], \; \mu \in \mathcal{D},
\end{multline}


\noindent with initial conditions: $u^e(x,0;\mu)=0$, $\frac{\partial u^e(x,0;\mu)}{\partial t}=0$.



In the above equation, $x$ denotes the coordinate of a point in the domain $\Omega$, $t$ is the time variable, $[0,T]$ is a finite time interval and the explicit forms of $a$, $m$, $c$ and $f$ could be defined as: $\forall w,v \in Y^e, \mu \in \mathcal{D}$,

\begin{subequations}\label{eq:ProbStat_four}
\begin{equation}
m(w,v;\mu) = \sum_{i} \int_{\Omega} \rho v_i w_i,\label{eq:ProbStat_four-a}
\end{equation}
\begin{equation}
c(w,v;\mu) = \sum_{i} \int_{\Omega} \alpha \rho v_i w_i + \sum_{i,j,k,l} \int_{\Omega} \beta \frac{\partial v_i}{\partial x_j} C_{ijkl} \frac{\partial w_k}{\partial x_l}, \label{eq:ProbStat_four-b}
\end{equation}
\begin{equation}
a(w,v;\mu) = \sum_{i,j,k,l} \int_{\Omega} \frac{\partial v_i}{\partial x_j} C_{ijkl} \frac{\partial w_k}{\partial x_l}, \label{eq:ProbStat_four-c}
\end{equation}
\begin{equation}
f(v;\mu) = \sum_{i} \int_{\Omega} b_i v_i + \sum_{i} \int_{\Gamma_N} v_i \phi_i. \label{eq:ProbStat_four-d}
\end{equation}
\end{subequations}


Where, $\rho$ is the mass density; $\alpha$ is the mass-proportional Rayleigh damping coefficient; $\beta$ is the stiffness-proportional Rayleigh damping coefficient; $C_{ijkl}$ is the material elasticity tensor; $b$ is a body force and $\phi$ is a surface traction applied to a region of the domain $\Omega$; $g(t)$ is the time history associated with the external loading $f(v;\mu)$; $\Gamma_D$ and $\Gamma_N$ are the Dirichlet and Neumann boundaries, respectively. We note that the input parameter $\mu$ could appear in (not limited to) either $\rho(\mu)$, $\alpha(\mu)$, $\beta(\mu)$, $C_{ijkl}(\mu)$ and even $b(\mu)$, $\phi(\mu)$ and $g(t;\mu)$.




We then evaluate a quantity of interest (output) from

\begin{equation}\label{eq:ProbStat_five}
s^e(\mu) = \int_{0}^{T} \int_{\Gamma_{\rm o}} u^e(x,t;\mu) \, \Sigma(x,t) \, dx \, dt, = \int_{0}^{T} \ell(u^e(x,t;\mu)) \, dt.
\end{equation}

Here, $\Gamma_{\rm o}$ are some (output) spatial regions of interest and $\Sigma(x,t)$ is an extractor which depends on the view position of an ``observer'' in the space-time domain; and $\ell(u^e(x,t;\mu)) = \int_{\Gamma_{\rm o}} u^e(x,t;\mu) \, \Sigma(x,t) dx$.



We shall assume that the bilinear forms $a(\cdot,\cdot;\mu)$ and $m(\cdot,\cdot;\mu)$ are continuous,

\begin{subequations}\label{eq:ProbStat_five_tmp1}
\begin{equation}
a(w,v;\mu) \le \gamma \|w\|_{Y^e} \|v\|_{Y^e} \le \gamma_0 \|w\|_{Y^e} \|v\|_{Y^e}, \qquad \forall w,v \in Y^e, \forall \mu \in \mathcal{D}, \label{eq:ProbStat_five_tmp1-a}
\end{equation}
\begin{equation}
m(w,v;\mu) \le \varrho \|w\|_{X^e} \|v\|_{X^e} \le \varrho_0 \|w\|_{X^e} \|v\|_{X^e}, \qquad \forall w,v \in Y^e, \forall \mu \in \mathcal{D},
\label{eq:ProbStat_five_tmp1-b}
\end{equation}
\end{subequations}


coercive,

\begin{subequations}\label{eq:ProbStat_five_tmp2}
\begin{equation}
0 \le \alpha_0 \le \alpha(\mu) \equiv \inf \limits_{v \in Y^e} \frac{a(v,v;\mu)}{\|v\|^2_{Y^e}}, \qquad \forall \mu \in \mathcal{D}, \label{eq:ProbStat_five_tmp2-a}
\end{equation}
\begin{equation}
0 \le \sigma_0 \le \sigma(\mu) \equiv \inf \limits_{v \in Y^e} \frac{m(v,v;\mu)}{\|v\|^2_{X^e}}, \qquad \forall \mu \in \mathcal{D}; \label{eq:ProbStat_five_tmp2-b}
\end{equation}
\end{subequations}


and symmetric $a(v,w;\mu)=a(w,v;\mu), \forall w,v \in Y^e, \forall \mu \in \mathcal{D}$ and $m(v,w;\mu)=m(w,v;\mu), \forall w,v \in Y^e, \forall \mu \in \mathcal{D}$. (We (plausibly) suppose that $\gamma_0$, $\varrho_0$, $\alpha_0$ and $\sigma_0$ may be chosen independent of $\mathcal{N}$.) In addition, the linear forms $f(v): Y^e \rightarrow \mathbb{R}$ and $\ell(v): Y^e \rightarrow \mathbb{R}$ are assumed to be bounded with respect to $\|\cdot\|_{Y^e}$ and $\|\cdot\|_{X^e}$, respectively. Under these conditions, there exists a unique so-called ``weak'' (or ``variational'') solution $u^e(x,t;\mu) \in Y^e$ of the equation \eqref{eq:ProbStat_three} \cite{Bangerth2010, lions1971optimal}.



We shall make an important assumption, that is, $a$, $m$, $c$ and $f$ depend affinely on the parameter $\mu$ and thus can be expressed as

\begin{subequations}\label{eq:ProbStat_six}
\begin{equation}
m(w,v;\mu) = \sum_{q=1}^{Q_m} \Theta^q_m(\mu)m^q(w,v), \quad \forall w,v \in Y^e, \mu \in \mathcal{D}, \label{eq:ProbStat_six-a}
\end{equation}
\begin{equation}
c(w,v;\mu) = \sum_{q=1}^{Q_c} \Theta^q_c(\mu)c^q(w,v), \quad \forall w,v \in Y^e, \mu \in \mathcal{D}, \label{eq:ProbStat_six-b}
\end{equation}
\begin{equation}
a(w,v;\mu) = \sum_{q=1}^{Q_a} \Theta^q_a(\mu)a^q(w,v), \quad \forall w,v \in Y^e, \mu \in \mathcal{D}, \label{eq:ProbStat_six-c}
\end{equation}
\begin{equation}
f(v;\mu) = \sum_{q=1}^{Q_f} \Theta^q_f(\mu)f^q(v), \quad \forall v \in Y^e, \mu \in \mathcal{D}, \label{eq:ProbStat_six-d}
\end{equation}
\end{subequations}


\noindent for some (preferably) small integers $Q_{m,c,a,f}$. Here, the smooth functions $\Theta^q_{m,c,a,f}(\mu): \mathcal{D} \rightarrow \mathbb{R}$ depend on $\mu$, but the bilinear and linear forms $m^q$, $c^q$, $a^q$ and $f^q$ do \textit{not} depend on $\mu$.





\subsection{Finite element discretization}\label{subsec_ProbStat_FEdiscretize}

We shall use the ``method of lines'' approach: fully discretize in space using Galerkin FE and discretize in time using Newmark's trapezoidal scheme $\left( \gamma^N=\frac{1}{2}, \beta^N=\frac{1}{4} \right)$. We introduce a reference finite element approximation space $Y \subset Y^e (\subset X^e)$ of dimension $\mathcal{N}$; we further define $X \equiv X^e$. Note that $Y$ and $X$ shall inherit the inner product and norm from $Y^e$ and $X^e$, respectively. For time integration: we divide $I = [0,T]$ into $K$ subintervals of equal length $\Delta t = \frac{T}{K}$ and define $t^k = k \Delta t$, $0 \le k \le K$. Recall that the Newmark's trapezoidal scheme is implicit and unconditionally stable. Furthermore, $\Delta t$ and the FE mesh size will be chosen such that they satisfy the solvability, stability and accuracy conditions following \cite{bathe1996finite} (Chapter 9.4.4) or \cite{Hughes1987} (Chapter 9.1). Thus, our ``true'' finite element approximation $u(x,t^k;\mu) \left(\equiv u^k(\mu)\right) \in Y$ to the ``exact'' problem is equivalent to solving $(K-1)$ following elliptic problems \cite{Hughes1987}:

\begin{equation}\label{eq:ProbStat_seven}
\mathcal{A} \left( u^{k+1}(\mu),v;\mu \right) = \mathcal{F} \left( v \right), \quad \forall v \in Y, \; \mu \in \mathcal{D}, \; 1 \le k \le K-1,
\end{equation}

\noindent where \cite{Daniel1997}

\begin{equation}\label{eq:ProbStat_eight}
\begin{cases}
\mathcal{A} \left( u^{k+1}(\mu),v;\mu \right) &= \displaystyle \frac{1}{\Delta t^2} m(u^{k+1}(\mu),v;\mu) + \frac{1}{2 \Delta t} c(u^{k+1}(\mu),v;\mu) + \frac{1}{4} a(u^{k+1}(\mu),v;\mu), \vspace{9pt} \\
\mathcal{F} \left( v \right) &= \displaystyle -\frac{1}{\Delta t^2} m(u^{k-1} (\mu),v;\mu) + \frac{1}{2 \Delta t} c(u^{k-1} (\mu),v;\mu) - \frac{1}{4} a(u^{k-1}(\mu),v;\mu) \vspace{9pt} \\
& \quad \displaystyle + \frac{2}{\Delta t^2} m(u^{k}(\mu),v;\mu) - \frac{1}{2} a(u^{k}(\mu),v;\mu) + g^{eq}(t^k) f(v;\mu), \vspace{9pt} \\
g^{eq}(t^k) &= \displaystyle \frac{1}{4} g(t^{k-1}) + \frac{1}{2} g(t^k) + \frac{1}{4} g(t^{k+1}),
\end{cases}
\end{equation}


%\begin{multline}\label{eq:ProbStat_six}
%\frac{1}{\Delta t^2} m(u(\mu,t^{k+1}),v;\mu) + \frac{1}{2 \Delta t} c(u(\mu,t^{k+1}),v;\mu) + \frac{1}{4} a(u(\mu,t^{k+1}),v;\mu) \\
%= - \frac{1}{\Delta t^2} m(u(\mu,t^{k-1}),v;\mu) + \frac{1}{2 \Delta t} c(u(\mu,t^{k-1}),v;\mu) - \frac{1}{4} a(u(\mu,t^{k-1}),v;\mu) \\
%+ \frac{2}{\Delta t^2} m(u(\mu,t^{k}),v;\mu) - \frac{1}{2} a(u(\mu,t^{k}),v;\mu) + g^{eq}(t^k) f(v;\mu), \quad \forall v \in Y, 1 \leq k \leq K-1,
%\end{multline}


\noindent with initial conditions\footnote{In order to start the procedure \eqref{eq:ProbStat_seven}, $u^1(\mu)$ is computed as on page 491 of \cite{Hughes1987}.}: $u^0(\mu)=0$, $\frac{\partial u^0(\mu)}{\partial t}=0$; we then evaluate the output of interest from (using the trapezoidal rule for integral approximation)

\begin{equation}\label{eq:ProbStat_nine}
s(\mu) = \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \ell(u(x,t;\mu)) \, dt \approx \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \ell(u^k(\mu)) + \ell(u^{k+1}(\mu)) \right).
\end{equation}


%For simplicity, we directly consider a time-discrete framework associated to the time interval $I = [0,T]$. We divide $I$ into $K$ subintervals of equal length $\Delta t = \frac{T}{K}$ and define $t^k = k \Delta t$, $0 \le k \le K$. We shall consider the Newmark's trapezoidal scheme with coefficients $\left( \gamma_s=\frac{1}{2}, \beta_s=\frac{1}{4} \right)$ for the time integration. Our results must be stable as $\Delta t \rightarrow 0, K \rightarrow \infty$.
%
%
%
%The ``exact'' semi-discrete problem is stated follows: given a parameter $\mu \in \mathcal{D} \subset \mathbb{R}^{P}$, we evaluate the (space-time) output of interest
%
%\begin{equation}\label{eq:ProbStat_three}
%s(\mu) = \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \int_{\Gamma_{\rm o}} u^e(x,t;\mu) \, \Sigma(x,t) dx dt = \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \ell(u^e(x,t;\mu)) dt,
%\end{equation}
%
%\noindent where the field variable, $u^e(\mu,t^k) \in Y^e$, $1 \le k \le K$, satisfies the weak form of the $\mu$-parametrized hyperbolic PDE \cite{Daniel1997}
%
%\begin{multline}\label{eq:ProbStat_four}
%\frac{1}{\Delta t^2} m(u^e(\mu,t^{k+1}),v;\mu) + \frac{1}{2 \Delta t} c(u^e(\mu,t^{k+1}),v;\mu) + \frac{1}{4} a(u^e(\mu,t^{k+1}),v;\mu) \\
%= - \frac{1}{\Delta t^2} m(u^e(\mu,t^{k-1}),v;\mu) + \frac{1}{2 \Delta t} c(u^e(\mu,t^{k-1}),v;\mu) - \frac{1}{4} a(u^e(\mu,t^{k-1}),v;\mu) \\
%+ \frac{2}{\Delta t^2} m(u^e(\mu,t^{k}),v;\mu) - \frac{1}{2} a(u^e(\mu,t^{k}),v;\mu) + g^{eq}(t^k) f(v;\mu), \quad \forall v \in Y^e, 1 \leq k \leq K-1,
%\end{multline}
%
%\noindent with
%
%\begin{equation}\label{eq:ProbStat_five}
%g^{eq}(t^k) = \frac{1}{4} g(t^{k-1}) + \frac{1}{2} g(t^k) + \frac{1}{4} g(t^{k+1}), \quad 1 \le k \le K-1,
%\end{equation}
%
%\noindent and initial conditions: $u^e(\mu,t^0)=0$, $\frac{\partial u^e(\mu,t^0)}{\partial t}=0$.
%
%
%Here, $\Gamma_{\rm o}$ is some (output) spatial regions of interest and $\Sigma(x,t)$ is an extractor which depends on the view position of an ``observer'' in the space-time domain; and $\ell(u^e(x,t;\mu)) = \int_{\Gamma_{\rm o}} u^e(x,t;\mu) \, \Sigma(x,t) dx$.


%In equation \eqref{eq:ProbStat_eight}, the FE mass and stiffness matrices are symmetric positive semidefinite


Clearly, with the well-conditions (i.e., symmetric positive definiteness) of the FE mass and stiffness matrices as well as of the initial values, the linear system \eqref{eq:ProbStat_eight} possesses a unique solution.


The RB approximation shall be built upon our reference finite element approximation, and the RB error will thus be evaluated with respect to $u^k(\mu) \in Y$. Clearly, our methods must remain computationally efficient and stable as $\mathcal{N} \rightarrow \infty$.



%We shall make the following assumptions. First, we assume that the bilinear forms $a(\cdot,\cdot;\mu)$ and $m(\cdot,\cdot;\mu)$ are continuous,
%
%\begin{subequations}\label{eq:ProbStat_eight}
%\begin{equation}
%a(w,v;\mu) \leq \gamma\|w\|_Y \|v\|_Y \leq \gamma_0\|w\|_Y \|v\|_Y, \quad \forall w,v \in Y, \forall \mu \in \mathcal{D},    \label{eq:ProbStat_eight-a}
%\end{equation}
%\begin{equation}
%m(w,v;\mu) \leq \rho\|w\|_X \|v\|_X \leq \rho_0\|w\|_X \|v\|_X, \quad \forall w,v \in Y, \forall \mu \in \mathcal{D},    \label{eq:ProbStat_eight-b}
%\end{equation}
%\end{subequations}
%
%\noindent coercive,
%
%\begin{subequations}\label{eq:ProbStat_nine}
%\begin{equation}
%0 \leq \alpha_0 \leq \alpha(\mu) \equiv \inf_{v \in Y} \frac{a(v,v;\mu)}{\|v\|^2_Y}, \quad \forall \mu \in \mathcal{D},    \label{eq:ProbStat_nine-a}
%\end{equation}
%\begin{equation}
%0 \leq \sigma_0 \leq \sigma(\mu) \equiv \inf_{v \in Y} \frac{m(v,v;\mu)}{\|v\|^2_X}, \quad \forall \mu \in \mathcal{D};    \label{eq:ProbStat_nine-b}
%\end{equation}
%\end{subequations}
%
%
%\noindent and symmetric $a(v,w;\mu)=a(w,v;\mu), \forall w,v \in Y, \forall \mu \in \mathcal{D}$, and $m(v,w;\mu)=m(w,v;\mu), \forall w,v \in X, \forall \mu \in \mathcal{D}$. (We (plausibly) suppose that $\gamma_0, \rho_0, \alpha_0$ and $\sigma_0$ may be chosen independent of $\mathcal{N}$ \cite{Grepl2005}). We also require that the linear forms $f(\cdot;\mu): Y \rightarrow \mathbb{R}$ and $\ell(\cdot): Y \rightarrow \mathbb{R}$ be bounded with respect to $\|\cdot\|_Y$ and $\|\cdot\|_X$, respectively.







Finally, note that our linear and bilinear forms are independent of time -- the system is thus linear time-invariant (LTI) \cite{Grepl2005}. We shall point out that one application which satisfies this property is the dental implant problem \cite{Hoang2013, Zaw2009}.





\subsection{Dealing with unknown loading}\label{subsec_ProbStat_ImpRes}

In many dynamical systems, generally, the applied force to excite the system (e.g., $g(t^k)$ in \eqref{eq:ProbStat_eight} or $g(t)$ in \eqref{eq:ProbStat_three}) is not known in advance and thus we cannot solve \eqref{eq:ProbStat_seven} for $u^{k+1}(\mu)$. In such situations, fortunately, we may appeal to the LTI property to justify an impulse approach as described now \cite{Grepl2005}. We note from the Duhamel's principle that the solution of any LTI system can be written as the convolution of the impulse response with the control input: for any control input $g_{\rm any}(t^k)$, we can obtain its corresponding solution $u^k_{\rm any}(\mu)$, $1 \leq k \le K$ from

\begin{equation}\label{eq:ProbStat_ten}
u^k_{\rm any}(\mu)=\sum_{j=1}^{k} u^{k-j+1}_{\rm unit}(\mu) \, g_{\rm any}(t^j), \quad 1 \le k \le K,
\end{equation}

\noindent where $u^k_{\rm unit}(\mu)$ is the solution of \eqref{eq:ProbStat_seven} for a unit impulse control input $g_{\rm unit}(t^k)=\delta_{1k}$, $1 \le k \le K$ ($\delta$ is the Kronecker delta symbol). Therefore, it is sufficient to build the RB basis functions for the problem based on this impulse response.






\section{Reduced basis approximation}\label{sec_RB}


Two key properties of the RB methodology will be recalled as follows. First, our attention is restricted to a \textit{smooth} and \textit{low-dimensional} manifold instead of the very high-dimensional FE space. Namely, the field variable $u^k(\mu), 1 \le k \le K$ does not belong to the very high-dimensional FE space; rather it resides, or ``evolves'' on a much lower dimensional manifold which is induced by the parametric dependence over the parameter domain \cite{Grepl2005}. Therefore, by restricting our attention to this manifold, we can adequately approximate the field variable by a space of dimension $N \ll \mathcal{N}$. Second, the parametric setting of the PDE \eqref{eq:ProbStat_three} enables to split the computational procedure into two stages: an extensive/expensive \textit{Offline} stage performed once to prepare all necessary data for numerous input-output calculations in the \textit{Online} stage afterwards. Details of these computations will be explained in subsequent sections.





%\hl{The well-posedness of RB here!}
%
% As described in many works on RB methodology, for instance \cite{Rozza2008, quarteroni2011certified}, the critical observation is that the field variable $u^k(\mu), 1 \le k \le K$ does not belong to a very high dimensional finite element space $Y$; rather it resides, or ``evolves'' on a much lower dimensional manifold which is induced by the parametric and temporal dependence \cite{Grepl2005}. Therefore, by restricting our attention to this manifold, we can adequately approximate the field variable by a space of dimension $N \ll \mathcal{N}$. This ...



\subsection{Approximation}\label{subsec_RB_Apprx}

We introduce the set of samples $S_* = \{\mu_1 \in \mathcal{D}, \mu_2 \in \mathcal{D}, \ldots, \mu_N \in \mathcal{D}\}, 1 \le N \le N_{\max} $, and associated nested Lagrangian RB spaces $Y_N = {\rm span} \{\zeta_n, 1 \le n \le N\}, 1 \le N \le N_{\max}$, where $\zeta_n \in Y_N, 1 \le n \le N_{\max}$ are mutually $(\cdot,\cdot)_Y$ -- orthogonal RB basis functions. The sets $S_*$ and $Y_N$ shall be constructed correspondingly by the standard and goal-oriented POD--Greedy algorithms described in Section \ref{subsec_RB_GOPODGreedy} afterwards.


Our reduced basis approximation $u^k_N(\mu)$ to $u^k(\mu)$ is then obtained by a standard Galerkin projection: given $\mu \in \mathcal{D}$, we now look for $u^k_N(\mu) \in Y_N$ that satisfies

\begin{equation}\label{eq:RB_one}
\mathcal{A} \left( u^{k+1}_N(\mu),v;\mu \right) = \mathcal{F} \left( v \right), \quad \forall v \in Y_N, \; \mu \in \mathcal{D}, \; 1 \le k \le K-1,
\end{equation}

\noindent where

\begin{equation}\label{eq:RB_two}
\begin{cases}
\mathcal{A} \left( u^{k+1}_N(\mu),v;\mu \right) &= \displaystyle \frac{1}{\Delta t^2} m(u^{k+1}_N(\mu),v;\mu) + \frac{1}{2 \Delta t} c(u^{k+1}_N(\mu),v;\mu) + \frac{1}{4} a(u^{k+1}_N(\mu),v;\mu), \vspace{9pt} \\
\mathcal{F} \left( v \right) &= \displaystyle -\frac{1}{\Delta t^2} m(u^{k-1}_N (\mu),v;\mu) + \frac{1}{2 \Delta t} c(u^{k-1}_N (\mu),v;\mu) - \frac{1}{4} a(u^{k-1}_N(\mu),v;\mu) \vspace{9pt} \\
& \quad \displaystyle + \frac{2}{\Delta t^2} m(u^{k}_N(\mu),v;\mu) - \frac{1}{2} a(u^{k}_N(\mu),v;\mu) + g^{eq}(t^k) f(v;\mu),
\end{cases}
\end{equation}


%\begin{multline}\label{eq:RB_one}
%\frac{1}{\Delta t^2} m(u_N(\mu,t^{k+1}),v;\mu) + \frac{1}{2 \Delta t} c(u_N(\mu,t^{k+1}),v;\mu) + \frac{1}{4} a(u_N(\mu,t^{k+1}),v;\mu) \\
%= - \frac{1}{\Delta t^2} m(u_N(\mu,t^{k-1}),v;\mu) + \frac{1}{2 \Delta t} c(u_N(\mu,t^{k-1}),v;\mu) - \frac{1}{4} a(u_N(\mu,t^{k-1}),v;\mu) \\
%+ \frac{2}{\Delta t^2} m(u_N(\mu,t^{k}),v;\mu) - \frac{1}{2} a(u_N(\mu,t^{k}),v;\mu) + g^{eq}(t^k) f(v;\mu), \quad \forall v \in Y_N, 1 \leq k \leq K-1,
%\end{multline}


\noindent with the initial conditions: $u^0_N(\mu)=0$, $\frac{\partial u^0_N(\mu)}{\partial t}=0$; we then evaluate the output estimate, $s_N(\mu)$, from

\begin{equation}\label{eq:RB_three}
s_N(\mu) = \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \ell(u_N(x,t;\mu)) \, dt \approx \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \ell(u_N^k(\mu)) + \ell(u_N^{k+1}(\mu)) \right).
\end{equation}





\subsection{Goal-oriented POD--Greedy sampling procedure}\label{subsec_RB_GOPODGreedy}

\subsubsection{The proper orthogonal decomposition}\label{subsubsec_RB_GOPODGreedy_POD}

We aim to generate an optimal (in the mean square error sense) basis set $\{ \zeta_m \}_{m=1}^{M}$ from any given set of $M_{\max} (\ge M)$ snapshots $\{\xi_k\}_{k=1}^{M_{\max}}$. To do this, let $V_M = {\rm span} \{v_1, \ldots, v_M\} \subset {\rm span} \{\xi_1, \ldots, \xi_{M_{\max}} \}$ be an arbitrary space of dimension $M$. We assume that the basis $\{v_1, \ldots, v_M\}$ is orthonormal such that $ (v_n,v_m) = \delta_{nm}, 1 \le n,m \le M$ ($(\cdot,\cdot)$ denotes an appropriate inner product and $\delta_{nm}$ is the Kronecker delta symbol). The POD space, $W_M = {\rm span} \{\zeta_1, \ldots, \zeta_M\} $ is defined as

\begin{equation}\label{eq:RB_four}
W_M = \arg \min_{V_M \subset \textrm{span} \{\xi_1,\ldots,\xi_{M_{\max}}\}} \left( \displaystyle \frac{1}{M_{\max}} \displaystyle \sum_{k=1}^{M_{\max}} \inf_{\underline{\alpha}^k \in \mathbb{R}^M} \biggl\| \xi_k - \sum_{m=1}^M \alpha_m^k v_m \biggr\|^2 \right).
\end{equation}


In essence, the POD space $W_M$ which is extracted from the given set of snapshots $\{\xi_k\}_{k=1}^{M_{\max}}$ is the space that best approximate this given set of snapshots and can be written as $W_M={\rm POD} \left( \{\xi_1,\ldots, \xi_{M_{\max}}\},M \right)$. We can construct this POD space by using the method of snapshots\footnote{Some books which investigate thoroughly this POD subject can be found in, for instance, \cite{volkwein2013proper, jolliffe2005principal}.} which is presented concisely in the Appendix of \cite{Nguyen2008}.




\subsubsection{Goal-oriented POD--Greedy algorithm}\label{subsubsec_RB_GOPODGreedy_GO}


%%%%%%%%%%%%% table 1 %%%%%%%%%%%%%
\begin{table}[h!]
\begin{center}
  {\begin{tabular}{| c | l | l |}
  \hline \rule{0pt}{3.0ex} \hspace{-2.2mm}
  (T\ref{tab1}a)   &   Set $Y^{\rm go}_N=0$   &   Set $Y^{\rm st}_N=0$ \\ [2ex]

  (T\ref{tab1}b)   &   Set $\mu^{\rm go}_{*}=\mu_{0}$   &   Set $\mu^{\rm st}_{*}=\mu_{0}$  \\ [2ex]

  (T\ref{tab1}c)   &   While $N \le N^{\rm go}_{\max}$   &   While $N \le N^{\rm st}_{\max}$  \\ [2ex]

  (T\ref{tab1}d)   &   \qquad $\mathcal{W}^{\rm go}  = \left\{e^{\rm go}_{\rm proj}(\mu^{\rm go}_{*},t^k), \, 0 \le k \le K \right\}$;   &  \qquad $\mathcal{W}^{\rm st} = \left\{e^{\rm st}_{\rm proj}(\mu^{\rm st}_{*},t^k), \, 0 \le k \le K \right\}$;   \\ [2ex]

  (T\ref{tab1}e)   &   \qquad $Y^{\rm go}_{N+M} \longleftarrow Y^{\rm go}_N \bigoplus {\rm POD}(\mathcal{W}^{\rm go},M)$;   &   \qquad $Y^{\rm st}_{N+M} \longleftarrow Y^{\rm st}_N \bigoplus {\rm POD}(\mathcal{W}^{\rm st},M)$; \\ [2ex]

  (T\ref{tab1}f)   &   \qquad $N \longleftarrow N + M$;    &   \qquad $N \longleftarrow N + M$;  \\ [2ex]

  (T\ref{tab1}g)   &   \qquad Find $\tilde{N}$ s.t. $\forall \mu \in \Xi^{\rm st}_{n} \subset \Xi^{\rm st}_{n+1} \; (\subset S^{\rm st}_*)$,  &   \qquad $\mu^{\rm st}_{*}= \arg \max\limits_{\mu \in \Xi_{\rm train}} \left\{ \frac{\boldsymbol{\Delta_u(\mu)}}{\sqrt{\sum_{k=1}^{K}\|u^{\rm st}_N(\mu,t^k)\|^2_Y }} \right\} $;   \\ [2ex]

   & \qquad \qquad $\eta_T \le \left| \frac{ \boldsymbol{\Delta_s(\mu)} }{s(\mu) - s^{\rm go}_N(\mu)} \right| \le 2-\eta_T$; & \qquad $S^{\rm st}_* \longleftarrow S^{\rm st}_* \bigcup \left\{ \mu^{\rm st}_* \right\} $;  \\ [2ex]

  (T\ref{tab1}h)   &   \qquad $\mu^{\rm go}_{*}= \arg \max\limits_{\mu \in \Xi_{\rm train}} \left\{ \left| \frac{ \boldsymbol{\Delta_s(\mu)} }{s^{\rm st}_{\tilde{N}}(\mu)} \right| \right\} $;   &  end.  \\ [2ex]

  (T\ref{tab1}i)   &   \qquad $S^{\rm go}_* \longleftarrow S^{\rm go}_* \bigcup \left\{ \mu^{\rm go}_* \right\} $;   &      \\ [2ex]

  (T\ref{tab1}j)   &       end.    &   \\
  \hline \rule{0pt}{6.0ex}

  (T\ref{tab1}k)   &   $\boldsymbol{\Delta_s(\mu)} = s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_N(\mu)$   &   $\boldsymbol{\Delta_u(\mu)} = \sqrt{\sum_{k=1}^{K} \|\mathcal{R}^{\rm st} (v;\mu,t^k)\|^2_{Y'} } $  \\ [4ex]
  \hline
\end{tabular}
\caption{(Left) The proposed goal-oriented POD--Greedy sampling algorithm and (Right) the standard POD--Greedy sampling algorithm. (The terms $\Delta_u(\mu)$ and $\Delta_s(\mu)$ are printed in boldface to highlight the difference between the two algorithms.)}
\label{tab1}}
%\tabnote{$^{\rm a}$This is not as accurate, owing to numerical
%error.} \tabnote{$^{\rm b}$An example table footnote to show the
%text turning over when a long footnote is inserted.}
\end{center}
\end{table}



%\begin{multline}\label{eq:RB_four}
%\mathcal{R}(v;\mu,t^k) = g^{eq}(t^k)f(v;\mu) \\
%- \frac{1}{\Delta t^2} \left( m(u_N(\mu,t^{k+1}),v;\mu) - 2m(u_N(\mu,t^{k}),v;\mu) + m(u_N(\mu,t^{k-1}),v;\mu) \right) \\
%- \frac{1}{\Delta t} \left( \frac{1}{2} c(u_N(\mu,t^{k+1}),v;\mu) - \frac{1}{2} c(u_N(\mu,t^{k-1}),v;\mu) \right) \\
%- \left( \frac{1}{4} a(u_N(\mu,t^{k+1}),v;\mu) + \frac{1}{2} a(u_N(\mu,t^{k}),v;\mu) + \frac{1}{4} a(u_N(\mu,t^{k-1}),v;\mu) \right);
%\end{multline}


%%%%%%%%% The ``cross-validation'' process %%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\resizebox*{5cm}{!}{\includegraphics{figs/diagram1.eps}}
%\caption{The ``cross-validation'' process to choose $\tilde{N}$.}
%\label{figtmp1}
%\end{center}
%\end{figure}



%%%%%%%%% Algorithm 1 %%%%%%%%%
\begin{center}
\noindent\fbox{%
\begin{varwidth}{\dimexpr\linewidth-2\fboxsep-2\fboxrule\relax}
\begin{algorithmic}[1]
\Statex \centerline{Algorithm 1: The ``cross-validation'' process.}
\Statex \textbf{INPUT}: $\Xi^{\rm st}_n, N$
\Statex \textbf{OUTPUT}: $\Xi^{\rm st}_n, \tilde{N}_n$
\Statex
\While {true}
    \State Given $\Xi^{\rm st}_{n}$;
    \State Compute $\tilde{N}_n$ from this $\Xi^{\rm st}_{n}$; (i.e., call the Algorithm 2 below)
    \State Create $\Xi^{\rm st}_{n+1}$;
    \State Check \eqref{eq:RB_five_tmp} with $\tilde{N}_n$ over $\Xi^{\rm st}_{n+1}$:

    \If{\eqref{eq:RB_five_tmp} holds $\forall \mu \in \Xi^{\rm st}_{n+1}$}       \Comment{if \eqref{eq:RB_five_tmp} holds}
        \State Get $\left\{ \tilde{N}_n, \Xi^{\rm st}_{n} \right\}$;
        \State Exit \textbf{while} loop;
    \Else                                                               \Comment{if \eqref{eq:RB_five_tmp} is violated}
        \State $\Xi^{\rm st}_{n} \leftarrow \Xi^{\rm st}_{n+1}$;
    \EndIf
\EndWhile
\end{algorithmic}
\end{varwidth}}
\end{center}



%%%%%%%%% Algorithm 2 %%%%%%%%%
\begin{center}
\noindent\fbox{%
\begin{varwidth}{\dimexpr\linewidth-2\fboxsep-2\fboxrule\relax}
\begin{algorithmic}[1]
\Statex \centerline{Algorithm 2: Function to compute $\tilde{N}$ based on input $\Xi^{\rm st}$ and $N$.}
\Statex \textbf{INPUT}: $\Xi^{\rm st}, N$
\Statex \textbf{OUTPUT}: $\tilde{N}$
\Statex
\State Start with $\tilde{N}=2N$;
\While{true}
\State Check \eqref{eq:RB_five_tmp} with $\tilde{N}$ over $\Xi^{\rm st}$:
\If{any $\mu \in \Xi^{\rm st}$ violates \eqref{eq:RB_five_tmp}}  \Comment{if \eqref{eq:RB_five_tmp} is violated}
    \If{$\tilde{N} < N^{\rm st}_{\max}$}
        \State $\tilde{N} \leftarrow \tilde{N} + 1$;
    \Else
        \State Run the standard POD--Greedy algorithm to increase $N^{\rm st}_{\max}$: $N^{\rm st}_{\max} \leftarrow N^{\rm st}_{\max}+1$;
    \EndIf
\EndIf

\If{\eqref{eq:RB_five_tmp} holds $\forall \mu \in \Xi^{\rm st}$}     \Comment{if \eqref{eq:RB_five_tmp} holds}
    \State Get $\tilde{N}$;
    \State Exit \textbf{while} loop;
\EndIf
\EndWhile
\end{algorithmic}
\end{varwidth}}
\end{center}




We now discuss the POD--Greedy algorithms \cite{Haasdonk2008, Hoang2013, eftang2011hp, haasdonk2013convergence} to construct the nested sets $S_*$ and $Y_N$ of interest. Let $\Xi_{\rm train}$ be a finite set of the parameters in $\mathcal{D}$ ($\Xi_{\rm train} \subset \mathcal{D}$); and $S_*$ denote the set of greedily selected parameters in $\Xi_{\rm train}$. Initialize $S_* = \{\mu_0\}$, where $\mu_0$ is an arbitrarily chosen parameter. Let $e_{\rm proj}(\mu,t^k) = u^k(\mu) - {\rm proj}_{Y_N}u^k(\mu)$, where ${\rm proj}_{Y_N}u^k(\mu)$ is the $Y_N$-orthogonal projection of $u^k(\mu)$ into the $Y_N$ space. Our proposed goal-oriented and the standard algorithms are presented simultaneously in Table \ref{tab1}. Note that the superscript ``st'' denotes the standard and ``go'' denotes the goal-oriented algorithms, respectively.

\begin{enumerate}[a)]

  \item Let us first describe the standard POD--Greedy algorithm in the right column of Table \ref{tab1}. In particular, at each Greedy iteration, one first solves \eqref{eq:ProbStat_seven} to obtain the ``true'' FE solution $\left\{ u^k(\mu^{\rm st}_*), 0 \le k \le K \right\}$; then computes the projection error to form the snapshots set $\mathcal{W}^{\rm st}$ in (T\ref{tab1}d). Next, one does snapshots-POD analysis \cite{Nguyen2008} on this set to build/extract $M$ ``new'' basis functions and adds them to the ``current'' RB basis functions set $Y^{\rm st}_N$ (T\ref{tab1}e). The RB dimension $N$ is also updated correspondingly in (T\ref{tab1}f). Now, based on these newly updated bases $Y^{\rm st}_N$, the error indicator $\frac{\Delta_u(\mu)}{\sqrt{\sum_{k=1}^{K}\|u^{\rm st}_N(\mu,t^k)\|^2_Y }}$ will be computed exhaustively for all $\mu \in \Xi_{\rm train}$ to look for the worst sampling point $\mu^{\rm st}_*$ and add it to the set $S^{\rm st}_*$ (T\ref{tab1}g). This whole procedure is iterated until it satisfies some stopping criteria. Note in the computation of $\Delta_u(\mu)$ (T\ref{tab1}k), the term $\|\mathcal{R}^{\rm st} (v;\mu,t^k)\|_{Y'}$ is the dual norm of the associated residual of equation \eqref{eq:RB_one}, namely,

      \begin{equation}\label{eq:RB_five}
      \mathcal{R}^{\rm st} (v;\mu,t^k) = \mathcal{F} \left( v \right) - \mathcal{A} \left( u^{k+1}_N(\mu),v;\mu \right), \quad \forall v \in Y_N, \; \mu \in \mathcal{D}, \; 1 \le k \le K-1.
      \end{equation}

      Detailed computation of the term $\|\mathcal{R}^{\rm st} (v;\mu,t^k)\|_{Y'}$ can be found in, for instance, \cite{Hoang2013, Rozza2008, Grepl2005}. In essence, the term $\frac{\Delta_u (\mu)}{\sqrt{\sum_{k=1}^{K}\|u^{\rm st}_N(\mu,t^k)\|^2_Y }}$ is the ratio of the dual norm of the residual to the norm of the RB solution. Thus, this term is roughly considered as an error indicator for the error in the solution (or field variable). This is one special choice of many ones and is usually used in the current standard POD--Greedy algorithm \cite{Haasdonk2008, Knezevic2011, Hoang2013}.



  \item We now consider the goal-oriented algorithm in the left column of Table \ref{tab1}. All the main steps of this algorithm are exactly the same as that of the standard algorithm, except that the error indicator is now defined as $\left| \frac{ \Delta_s(\mu) }{s^{\rm st}_{\tilde{N}}(\mu)} \right|$ in (T\ref{tab1}h). Essentially, this term is an asymptotic error approximation for the true output error (i.e., $s(\mu) - s_N(\mu) \approx s^{\rm st}_{\tilde{N}}(\mu) - s_N(\mu)$). Thus, the main idea of this goal-oriented algorithm is that the Greedy iterations now try to minimize the output error indicator rather than the solution error indicator as in the standard algorithm above. From another viewpoint, we can think of this goal-oriented algorithm as a special version of the standard algorithm using a relative output error approximation as an error indicator.



  \item In the goal-oriented algorithm, the computation of $\Delta_s(\mu)$ requires a good approximation $s^{\rm st}_{\tilde{N}}(\mu)$ for the FE output $s(\mu)$. To cope with this situation, we use the asymptotic output error approximation which makes use of an enriched RB output computation. Specifically, we propose to implement the standard algorithm \textit{in advance} to obtain enriched RB spaces $Y^{\rm st}_{N_{\max}}$ which are ready for the computation of $s^{\rm st}_{\tilde{N}}(\mu)$ ($\tilde{N} \ge N$) in the goal-oriented algorithm \textit{afterwards}. Hence, here $s^{\rm st}_{\tilde{N}}(\mu)$ implies the enriched RB output computed by standard algorithm (``st'') using $\tilde{N}$ basis functions.


      The use of asymptotic output error approximation $\Delta_s(\mu)$ is not new in the literature. For example, in the works \cite{Liu2008, Prudhomme2002} the authors used heuristically $\tilde{N}=2N$ as an adequate choice for the output error approximation. In this work, however, we propose a new algorithm called a ``cross-validation'' process to choose adaptively $\tilde{N}$ for each particular $N$ in the offline stage of the goal-oriented algorithm. Note that these found pairs $(N,\tilde{N})$ will also be used in the error approximation in the online stage later. The ``cross-validation'' process is presented in line (T\ref{tab1}g) and is detailed in Algorithm 1; however, we postpone its explanation until point f) below.



  \item The main idea of the ``cross-validation'' process is to find a sufficient $\tilde{N}$ (for an $N$ under consideration) such that the effectivity $\left| \frac{ s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_N(\mu) }{s(\mu) - s^{\rm go}_N(\mu)} \right|$ satisfies: $\forall \mu \in \Xi^{\rm try}$,

      \begin{equation*}\label{eq:RB_five_tmp}
      \eta_T \le \left| \frac{ s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_N(\mu) }{s(\mu) - s^{\rm go}_N(\mu)} \right| \le 2-\eta_T, \tag{$\clubsuit$}
      \end{equation*}

      where $\eta_T$ is a user-prescribed effectivity (say, 0.8 or 0.9), and $\Xi^{\rm try}$ is an arbitrary set of parameters ($s^{\rm go}_N(\mu)$ is a goal-oriented RB output using $N$ basis functions). Because the FE output $s(\mu)$ appears in \eqref{eq:RB_five_tmp}, we think of using $\Xi^{\rm try} \subset S^{\rm st}_*$ as all FE solutions/outputs are available for all $\mu \in S^{\rm st}_*$ (recall that the standard algorithm was \textit{already} implemented). Hence, we will use the notation $\Xi^{\rm st} (\equiv \Xi^{\rm try}) \subset S^{\rm st}_*$ to reflect this idea.



  \item We first explain the Algorithm 2 in detail. Essentially, algorithm 2 is an iteration process to find a proper $\tilde{N}$ which satisfies \eqref{eq:RB_five_tmp} for a given sample set $\Xi^{\rm st}$ and a given $N$. The iteration procedure starts with $\tilde{N}=2N$; and $\tilde{N}$ will be increased if \eqref{eq:RB_five_tmp} is violated by any $\mu \in \Xi^{\rm st}$. The procedure will stop when \eqref{eq:RB_five_tmp} holds true $\forall \mu \in \Xi^{\rm st}$. Within Algorithm 2, if $\tilde{N}$ exceeds $N^{\rm st}_{\max}$ the standard algorithm will be called and implemented to increase $N^{\rm st}_{\max}$ accordingly. Note that the standard algorithm will continue to run from the previous $N^{\rm st}_{\max}$; so in summary, it is deemed to run ``once'' but continuously in different stages when necessary.



  \item Let us now describe the ``cross-validation'' process in Algorithm 1. In fact, finding the proper size of $\Xi^{\rm st}$ is not a trivial task: fixing $\Xi^{\rm st} = S^{\rm st}_*$ is not an efficient way, and we also don't want to use one more parameter to tune this setting (the only one parameter for the GO algorithm is $\eta_T$). We thus propose an adaptive strategy to choose $\left| \Xi^{\rm st} \right|$ as follows (Algorithm 1). Suppose that at the Greedy iteration $N$ with a given set $\Xi^{\rm st}_n \subset S^{\rm st}_*$, we can find the corresponding $\tilde{N}$ thanks to Algorithm 2 above. We will request further that the currently found $\tilde{N}$ also needs to verify \eqref{eq:RB_five_tmp} over $\Xi^{\rm st}_{n+1}$, where $\left| \Xi^{\rm st}_{n+1} \right| = \left| \Xi^{\rm st}_{n} \right| + \Delta n_{\rm sample}$, and $\Delta n_{\rm sample}$ is a number of next sample points taken from the set $S^{\rm st}_*$. Otherwise, $\Xi^{\rm st}_{n+1}$ is assigned to $\Xi^{\rm st}_{n}$ (i.e., $\Xi^{\rm st}_{n}$ is enriched now) and the procedure is repeated until the found $\tilde{N}$ satisfies \eqref{eq:RB_five_tmp} over both sets $\Xi^{\rm st}_{n}$ and $\Xi^{\rm st}_{n+1}$ ($\subset S^{\rm st}_*$). By this way, we can start the sampling procedure with fairly small $\left| \Xi^{\rm st}_{n} \right|$ and let it ``evolve'' automatically without the necessity of any control or adjustment. Following the same process, $\Xi^{\rm st}_{n}$ and $\Xi^{\rm st}_{n+1}$ will also be used to find $\tilde{N}$ at the next Greedy iteration $N+1$; and they will be enriched appropriately when necessary.


      For example, the standard algorithm is implemented \textit{first} and \textit{once} for quite large $N^{\rm st}_{\max}$, say, $N^{\rm st}_{\max}=200$ and hence $\left|S^{\rm st}_*\right|=200$ (using $M=1$). Then for an arbitrary quantity of interest, the goal-oriented algorithm will be implemented accordingly. Consider the ``cross-validation'' process at the first Greedy iteration $N^{\rm go}=1$, we can choose $\Xi^{\rm st}_1 \subset \Xi^{\rm st}_2 \subset S^{\rm st}_*$ such that $\left| \Xi^{\rm st}_1 \right|=10$ and $\left| \Xi^{\rm st}_2 \right|=20$ first sample points of $S^{\rm st}_*$, respectively (hence, $\Delta n_{\rm sample}=10$). Based on this $\Xi^{\rm st}_1$ set, Algorithm 2 is invoked to find the corresponding $\tilde{N}_1$. Next, \eqref{eq:RB_five_tmp} is checked with $\tilde{N}_1$ over $\Xi^{\rm st}_2$: if it holds true $\forall \mu \in \Xi^{\rm st}_2$, then $\left( N^{\rm go}=1,\tilde{N}_1 \right)$ will be the necessarily found pair; the algorithm will quit the ``cross-validation'' process and continue with step (T\ref{tab1}h). Otherwise, the algorithm will enrich $\Xi^{\rm st}_1$ and $\Xi^{\rm st}_2$ such that $\left| \Xi^{\rm st}_1 \right|=20$ and $\left| \Xi^{\rm st}_2 \right|=30$ first sample points from $S^{\rm st}_*$, and repeat the computational procedure until the right pair $\left( N^{\rm go}=1,\tilde{N} \right)$ is found. $\Xi^{\rm st}_1$ and $\Xi^{\rm st}_2$ are also used to check \eqref{eq:RB_five_tmp} at the next Greedy iteration in a completely similar manner.


  \item Finally, we close this subsection with one remark on the possible value range of $\eta_T$. In fact, we cannot choose $\eta_T$ to be too high, i.e., too close to 1. The reason is that the convergence of the GO algorithm depends on the convergence of the standard one; and generally the standard algorithm will stall/flat with some $N \ge N^{\rm st}_{\max}$ (i.e., its error indicator and RB true error cannot decrease further as it reaches machine accuracy $\approx 10^{-8}$). If $\eta_T$ is too close to 1, say 0.95, the cross-validation process will break down (infinite loop in Algorithm 1) since it cannot find the suitable $\tilde{N}$ to satisfy ($\clubsuit$) over $\Xi^{\rm st}_{n}$ and $\Xi^{\rm st}_{n+1}$; and there is no way to handle this situation. Therefore, it is practical to choose a modest $\eta_T$, and in fact we can do that easily based on the convergence history of the standard algorithm which is implemented \textit{in advance}. Indeed, through two numerical experiments in Section \ref{sec_NumEx} later, $\eta_T=0.8$ is the maximum possible choice and it gives the best performance among all GO $\eta_T$ algorithms. On the contrary, low $\eta_T$ (i.e., close to 0) poses no difficulty for the cross-validation process since its corresponding $\tilde{N}$ generally will be smaller than that of higher $\eta_T$, and hence low $\eta_T$ is ``safer'' than high $\eta_T$ regarding $\tilde{N}$.




% The ``cross-validation'' process is repeated similarly at later Greedy iterations with $\Xi^{\rm st}_1$ and $\Xi^{\rm st}_2$.


% Note that from Algorithm 1 and 2 above, we see that there are two ways from which the goal-oriented algorithm can call the standard algorithm (i.e., to enlarge the standard RB space): either by increasing $\tilde{N}$ or by extending $\Xi^{\rm st}_n$. For simplicity, in this work we usually choose $\left| \Xi^{\rm st}_1 \right|$ and $\Delta n_{\rm sample}$ quite small compared with $\left|S^{\rm st}_*\right|$ so that the standard RB space is enriched only by $\tilde{N}$.


% verifies \eqref{eq:RB_five_tmp} not only over $\Xi^{\rm st}_n$ but also over $\Xi^{\rm st}_{n+1}$, where $\left| \Xi^{\rm st}_{n+1} \right| = \left| \Xi^{\rm st}_{n} \right| + \Delta n_{\rm sample}$. (Here, $\Delta n_{\rm sample}$ is a number of next sample points taken from the set $S^{\rm st}_*$.) Otherwise, assign $\Xi^{\rm st}_{n+1}$ to $\Xi^{\rm st}_{n}$ and repeat the procedure.


%An alternative choice is to fix $\Xi^{\rm st} = S^{\rm st}_*$, then we do not need the ``cross-validation'' process. However, we have no \textit{a priori} knowledge about how to choose $\Xi^{\rm st}$ appropriately; and we want to let $\Xi^{\rm st}$ grows adaptively and automatically whenever it needs to be enriched (i.e., violates ($\clubsuit$)) rather than fix it at constant values.


\end{enumerate}


% The offline stage of the goal-oriented algorithm will now include the offline stage of the standard algorithm (for the error estimation $\Delta_s(\mu)$); thus its offline computational cost is quite expensive. However, the online computational cost of the goal-oriented algorithm is very similar to that of the standard algorithm, i.e. $O(N^3)$ and $O(\tilde{N}^3)$ (and \textit{not} $O(\mathcal{N})$) as described in Section \ref{subsec_RB_OffOn}.


%On the contrary, the term $\left| \frac{\Delta_s(\mu)}{s^{\rm st}_{\tilde{N}}(\mu)} \right|$ is an approximation for the true output error (i.e., $s(\mu) - s_N(\mu) \approx s^{\rm st}_{\tilde{N}}(\mu) - s_N(\mu)$). Therefore, as observed from Table \ref{tab1}, the main difference between these two algorithms is that we somehow try to minimize the error indicator of the output functional ($\Delta_s (\mu)$) rather than an error indicator of the solution ($\Delta_u (\mu)$) by means of Greedy iterations. Therefore, we can think of this goal-oriented algorithm as a special version of the standard POD--Greedy algorithm using a relative output error approximation as an error indicator. By this way, we expect to improve further the accuracy of the RB output computation; but consequently, we might lose the rapid convergent rate of the solution as in the standard algorithm.


%We emphasize that in order to compute the output error approximation $\Delta_s(\mu)$, the standard POD--Greedy algorithm is invoked to provide $s^{\rm st}_{\tilde{N}}(\mu)$ $\left( \approx s(\mu) \right)$. Generally, $\tilde{N}=2N$ is used to be an adequate choice (as in \cite{Liu2008, Prudhomme2002}). However, here we propose a simple algorithm (or a ``cross-validation'' process) to choose $\tilde{N}$ adaptively as described in Algorithm 1 (note that the left column of line (T\ref{tab1}g) in Table \ref{tab1} is exactly Algorithm 1.) The idea is that we want to find $\tilde{N}$ such that the effectivity of $\Delta_s(\mu)$ lies in a range $[\eta_T, 2-\eta_T]$, $\forall \mu \in \Xi^{\rm st} \; (\subset S^{\rm st}_*)$. Here, $\eta_T$ is a prescribed effectivity (say, 0.9 or 0.8); and $\Xi^{\rm st}$ is a subset of $S^{\rm st}_*$ which is resulted from the standard algorithm (and hence, all the FE outputs of these sample points are already available for the computation ($\clubsuit$)).


%As regards the implementation of the goal-oriented algorithm, we propose to run the standard algorithm \textit{first} for quite large RB space $Y^{\rm st}_N$ (say, $N^{\rm st}_{\max}=200$ or 300), because we need the set $S^{\rm st}_*$ for the ``cross-validation'' process within the goal-oriented algorithm. Once the standard algorithm is done (i.e., $Y^{\rm st}_N$ and $S^{\rm st}_*$ are available), then the goal-oriented algorithm can be implemented. Consider the ``cross-validation'' process in Algorithm 1, we first create the sample set $\Xi^{\rm st}_{1} \subset S^{\rm st}_*$ where $\left| \Xi^{\rm st}_{1} \right|=n_{\rm sample}$ is small (say, 10 first sample points) and set $\Delta n_{\rm sample}$ (say, 10). The next sample set to validate, ($\Xi^{\rm st}_{1} \subset$) $\Xi^{\rm st}_{2}$ ($\subset S^{\rm st}_{*}$), is defined as: $\left| \Xi^{\rm st}_{2} \right| = \left| \Xi^{\rm st}_{1} \right| + \Delta n_{\rm sample}$, where $\Delta n_{\rm sample}$ is the number of next sample points taken from $S^{\rm st}_*$. The main idea of this ``cross-validation'' process is that the current found $\tilde{N}_n$ needs to validate ($\clubsuit$) not only over $\Xi^{\rm st}_{n}$ but also over $\Xi^{\rm st}_{n+1}$ (which is larger and also includes $\Xi^{\rm st}_{n}$). Otherwise, $\tilde{N}$ and $\Xi^{\rm st}$ need to be enriched accordingly\footnote{An alternative choice is to fix $\Xi^{\rm st} = S^{\rm st}_*$, then we do not need the ``cross-validation'' process. However, we have no \textit{a priori} knowledge about how to choose $\Xi^{\rm st}$ appropriately; and we want to let $\Xi^{\rm st}$ grows adaptively and automatically whenever it needs to be enriched (i.e., violates ($\clubsuit$)) rather than fix it at constant values.}.


%  \item In order to implement the goal-oriented POD--Greedy algorithm, we need to compute $s^{\rm st}_{\tilde{N}}(\mu)$ for the error estimate $\Delta_s(\mu)$. Therefore, we proposed implementing the standard POD--Greedy algorithm \textit{first} to build the basis set $Y^{\rm st}_N$; and then the goal-oriented POD--Greedy algorithm is performed \textit{subsequently}. That means the offline stage of the goal-oriented algorithm is now more expensive than that of the standard one (as it also includes the offline stage of the standard one). However, the online computational cost of the goal-oriented algorithm is completely similar to that of the standard algorithm as shown in Section \ref{subsec_RB_OffOn}.


% However, it is worth as we will see in Section \ref{sec_NumRes} later, the convergent rate of the true output error by the goal-oriented one is much faster than that of the standard one.





\subsection{Error estimations}\label{subsec_RB_EE}

\subsubsection{True errors}

We use the true errors to validate the performance of the standard and goal-oriented algorithms in the online computation stage. The relative true errors by the two algorithms for the solutions are defined as


\begin{equation}\label{eq:RB_six}
e^{\rm st}_{u}(\mu) = \sqrt{ \frac{ \sum_{k=1}^{K} \|u^{k}(\mu) - u^{{\rm st},k}_N(\mu)\|^2_Y }{ \sum_{k=1}^{K} \|u^{k}(\mu)\|^2_Y } }, \qquad
e^{\rm go}_{u}(\mu) = \sqrt{ \frac{ \sum_{k=1}^{K} \|u^{k}(\mu) - u^{{\rm go},k}_N(\mu)\|^2_Y }{ \sum_{k=1}^{K} \|u^{k}(\mu)\|^2_Y } };
\end{equation}



%\begin{subequations}\label{eq:RB_six}
%\begin{equation}
%e^{\rm st}_{u}(\mu) = \frac{ \int_{0}^{T} \|u(x,t;\mu) - u^{\rm st}_N(x,t;\mu) \|_Y \, dt}{ \int_{0}^{T} \| u(x,t;\mu) \|_Y \, dt} \approx \frac{\displaystyle \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \|u^k(\mu)-u^{{\rm st},k}_N(\mu)\|_Y + \|u^{k+1}(\mu)-u^{{\rm st},k+1}_N(\mu)\|_Y \right)}{\displaystyle \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \|u^k(\mu)\|_Y + \|u^{k+1}(\mu)\|_Y \right)}, \label{eq:RB_six-a}
%\end{equation}
%\begin{equation}
%e^{\rm go}_{u}(\mu) = \frac{ \int_{0}^{T} \|u(x,t;\mu) - u^{\rm go}_N(x,t;\mu) \|_Y \, dt}{ \int_{0}^{T} \| u(x,t;\mu) \|_Y \, dt} \approx \frac{\displaystyle \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \|u^k(\mu)-u^{{\rm go},k}_N(\mu)\|_Y + \|u^{k+1}(\mu)-u^{{\rm go},k+1}_N(\mu)\|_Y \right)}{\displaystyle \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \|u^k(\mu)\|_Y + \|u^{k+1}(\mu)\|_Y \right)}; \label{eq:RB_six-b}
%\end{equation}
%\end{subequations}


%\begin{equation}\label{eq:RB_five}
%e^{\rm st}_{u}(\mu) = \frac{ \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \|u(x,t;\mu) - u^{\rm st}_N(x,t;\mu) \|_Y \, dt}{ \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \| u(x,t;\mu) \|_Y \, dt}, \quad e^{\rm go}_{u}(\mu) = \frac{ \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \|u(x,t;\mu) - u^{\rm go}_N(x,t;\mu) \|_Y \, dt}{ \sum_{k=0}^{K-1} \int_{t^k}^{t^{k+1}} \| u(x,t;\mu) \|_Y \, dt};
%\end{equation}


\noindent and for the outputs

\begin{subequations}\label{eq:RB_seven}
\begin{equation}
e^{\rm st}_{s}(\mu) = \left| \frac{ s(\mu) - s^{\rm st}_N(\mu) }{ s(\mu) } \right| \approx \left| \frac{ \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left(\ell\left(u^k(\mu) - u^{{\rm st},k}_N(\mu)\right) + \ell\left(u^{k+1}(\mu) - u^{{\rm st},k+1}_N(\mu)\right)\right)}{ \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \ell(u^k(\mu))+\ell(u^{k+1}(\mu)) \right)} \right|, \label{eq:RB_seven-a}
\end{equation}
\begin{equation}
e^{\rm go}_{s}(\mu) = \left| \frac{ s(\mu) - s^{\rm go}_N(\mu) }{ s(\mu) } \right| \approx \left| \frac{ \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left(\ell\left(u^k(\mu) - u^{{\rm go},k}_N(\mu)\right) + \ell\left(u^{k+1}(\mu) - u^{{\rm go},k+1}_N(\mu)\right)\right)}{ \sum_{k=0}^{K-1} \frac{\Delta t}{2} \left( \ell(u^k(\mu))+\ell(u^{k+1}(\mu)) \right)} \right|, \label{eq:RB_seven-b}
\end{equation}
\end{subequations}


\noindent respectively.


In the above expressions, $u^{k}(\mu)$, $u^{{\rm st},k}_N(\mu)$ and $u^{{\rm go},k}_N(\mu)$ are the FE, standard RB and goal-oriented RB solutions; $s(\mu)$, $s^{\rm st}_N(\mu)$ and $s^{\rm go}_N(\mu)$ are the FE, standard RB and goal-oriented RB outputs, respectively.



\subsubsection{Output error approximation}

The true errors are good for validation purposes but are not of practical uses in the online stage, where one requires fast and countless online calculations. In this work, we propose to use $\Delta_s(\mu)$ as an error estimation for the output in the online computation stage (in short, $\Delta_s(\mu)$ is an error estimation in both offline and online stages). Of course this error approximation is not a rigorous upper error bound such as the \textit{a posteriori} error bounds \cite{Rozza2008, Nguyen2009, Knezevic2011}; however, there are several good reasons for using it in practice. First, the time-marching error bounds for the wave equation so far were shown to be ineffective and pessimistic due to the instability of the wave equation: exponential growing with respect to time \cite{HoangPhDThesis2012, Yano2012a, Yano2014}. (We also note that the space-time error bounds, although very promising, are still not yet derived in the literature.) Second, this error approximation converges asymptotically to the true error (thanks to $\tilde{N}$ chosen effectively by the proposed ``cross-validation'' process), and thus can approximate relatively the accuracy of the RB outputs for various choices of $\mu$. Thirdly -- most important, its computational cost is cheap: only $O(N^3) + O(\tilde{N}^3)$ as described in the next section ($N, \tilde{N} \ll \mathcal{N}$). The output error approximation $\Delta_{s}(\mu)$ in (T\ref{tab1}k) and its associated effectivity are defined as follows

\begin{equation}\label{eq:RB_eight}
\Delta^{\rm go}_{s}(\mu) = s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_{N}(\mu), \qquad
\eta^{\rm go}_{s}(\mu) = \left| \frac{ \Delta^{\rm go}_s(\mu) }{ s(\mu) - s^{\rm go}_{N}(\mu) } \right|.
\end{equation}


% \ref{subsubsec_NumEx_3dDental_NumRes}

Note that to compare the performance of output error approximation of the standard versus the goal-oriented algorithms, here we also define the output error approximation and its associated effectivity for the standard algorithm as

\begin{equation}\label{eq:RB_nine}
\Delta^{\rm st}_{s}(\mu) = s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm st}_{N}(\mu), \qquad
\eta^{\rm st}_{s}(\mu) = \left| \frac{ \Delta^{\rm st}_s(\mu) }{ s(\mu) - s^{\rm st}_{N}(\mu) } \right|,
\end{equation}

where the superscript ``st'' implies the standard algorithm. The superscript ``go'' is thus added in \eqref{eq:RB_eight} to imply the goal-oriented algorithm, respectively.



%\begin{equation}\label{eq:RB_eight}
%\eta_{s}(\mu) = \left| \frac{ s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_{N}(\mu) }{ s(\mu) - s^{\rm go}_{N}(\mu) } \right| = \frac{ \Delta_s(\mu) }{ e^{\rm go}_s (\mu) }.
%\end{equation}


%\begin{subequations}\label{eq:RB_eight}
%\begin{equation}
%\Delta_{s}(\mu) = \left| \frac{ s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_N(\mu) }{ s^{\rm st}_{\tilde{N}}(\mu) } \right|, \label{eq:RB_eight-a}
%\end{equation}
%\begin{equation}
%\eta_{s}(\mu) = \left| \frac{ s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_{N}(\mu) }{ s(\mu) - s^{\rm go}_{N}(\mu) } \right| = \frac{ \hat{e}_s(\mu) }{ e^{\rm go}_s (\mu) }, \label{eq:RB_eight-b}
%\end{equation}
%\end{subequations}
%
%\noindent respectively.



%In addition, the effectivities corresponding to these asymptotic output errors can be defined as
%
%\begin{subequations}\label{eq:RB_nine}
%\begin{equation}
%\eta^{\rm off}_{s}(\mu) = \left| \frac{ s^{\rm st}_{2N}(\mu) - s^{\rm go}_{N}(\mu) }{ s(\mu) - s^{\rm go}_{N}(\mu) } \right| = \frac{ \hat{e}^{\rm off}_s(\mu) }{ e^{\rm go}_s (\mu) },    \label{eq:RB_nine-a}
%\end{equation}
%\begin{equation}
%\eta^{\rm onl}_{s}(\mu) = \left| \frac{ s^{\rm go}_{2N}(\mu) - s^{\rm go}_{N}(\mu) }{ s(\mu) - s^{\rm go}_{N}(\mu) } \right| = \frac{ \hat{e}^{\rm onl}_s(\mu) }{ e^{\rm go}_s (\mu) }.    \label{eq:RB_nine-b}
%\end{equation}
%\end{subequations}



%We note that $\hat{e}^{\rm off}_{s}(\mu)$ is similar to $\Delta_s(\mu)$ -- the error indicator for the output in the offline stage; and $\hat{e}^{\rm onl}_{s}(\mu)$ is the useful asymptotic error practically used in the online stage of the goal-oriented POD--Greedy algorithm.




\subsection{Offline-online computational procedure}\label{subsec_RB_OffOn}

In this section, we develop offline-online computational procedures in order to fully exploit the dimension reduction of the problem \cite{Grepl2005, Hoang2013, Nguyen2009}. We note that both algorithms (standard and goal-oriented) have the same offline-online computational procedures, they are only different in the ways to build the sets $S_*$ and $Y_N$ via Greedy iterations. We first express $u^k_N(\mu)$ as:

\begin{equation}\label{eq:RB_ten}
u^k_N(\mu) = \sum_{n=1}^{N} u^k_{N \, n}(\mu) \; \zeta_n,   \qquad \forall \zeta_n \in Y_N.
\end{equation}


We then choose a test function $v = \zeta_n, 1 \le n \le N$ for the RB equation \eqref{eq:RB_one}. It then follows that $\underline{u}^k_N(\mu) = [u^k_{N \, 1}(\mu) \; u^k_{N \, 2}(\mu) \; \cdots \; u^k_{N \, N}(\mu) ]^T \in \mathbb{R}^N$ satisfies

\begin{multline}\label{eq:RB_eleven}
\left( \frac{1}{\Delta t^2} \mathbf{M}_N(\mu) + \frac{1}{2 \Delta t} \mathbf{C}_N(\mu) + \frac{1}{4} \mathbf{A}_N(\mu) \right) \underline{u}^{k+1}_N(\mu) \\
= \left( - \frac{1}{\Delta t^2} \mathbf{M}_N(\mu) + \frac{1}{2 \Delta t} \mathbf{C}_N(\mu) - \frac{1}{4} \mathbf{A}_N(\mu) \right) \underline{u}^{k-1}_N(\mu) \\
+ \left( \frac{2}{\Delta t^2} \mathbf{M}_N(\mu) - \frac{1}{2} \mathbf{A}_N(\mu) \right) \underline{u}^{k}_N(\mu) + g^{eq}(t^k) \mathbf{F}_N(\mu), \qquad 1 \leq k \leq K-1.
\end{multline}

The initial condition is treated similar to the treatment in \eqref{eq:ProbStat_seven} and \eqref{eq:RB_one}. Here, $\mathbf{C}_N(\mu)$, $\mathbf{A}_N(\mu)$, $\mathbf{M}_N(\mu) \in \mathbb{R}^{N \times N}$ are symmetric positive definite matrices\footnote{The proof of this property can be found in, for instance, Proposition 5.1, page 136 of \cite{saad2003iterative}. Thanks to this property, the stability of our proposed RB scheme will be guaranteed as a consequence.} with entries $\mathbf{C}_{N \, i,j}(\mu) = c(\zeta_i,\zeta_j;\mu)$, $\mathbf{A}_{N \, i,j}(\mu) = a(\zeta_i,\zeta_j;\mu)$, $\mathbf{M}_{N \, i,j}(\mu) = m(\zeta_i,\zeta_j;\mu)$, $1 \le i,j \le N$ and $\mathbf{F}_N \in \mathbb{R}^N$ is the RB load vector with entries $\mathbf{F}_{N \, i} = f(\zeta_i)$, $1 \le i \le N$, respectively.


The RB output is then computed from

\begin{equation}\label{eq:RB_twelve}
s_N(\mu) = \sum_{k=0}^{K-1} \frac{\Delta t}{2} \, \mathbf{L}^T_N \left(\underline{u}^k_N(\mu) + \underline{u}^{k+1}_N(\mu) \right).
\end{equation}


Invoking the affine parameter dependence \eqref{eq:ProbStat_six}, we obtain

\begin{subequations}\label{eq:RB_thirteen}
\begin{equation}
\mathbf{M}_{N \, i,j}(\mu) = m(\zeta_i,\zeta_j;\mu) = \sum_{q=1}^{Q_m} \Theta^q_m(\mu) \, m^q(\zeta_i,\zeta_j), \label{eq:RB_thirteen-a}
\end{equation}
\begin{equation}
\mathbf{C}_{N \, i,j}(\mu) = c(\zeta_i,\zeta_j;\mu) = \sum_{q=1}^{Q_c} \Theta^q_c(\mu) \, c^q(\zeta_i,\zeta_j), \label{eq:RB_thirteen-b}
\end{equation}
\begin{equation}
\mathbf{A}_{N \, i,j}(\mu) = a(\zeta_i,\zeta_j;\mu) = \sum_{q=1}^{Q_a} \Theta^q_a(\mu) \, a^q(\zeta_i,\zeta_j), \label{eq:RB_thirteen-c}
\end{equation}
\begin{equation}
\mathbf{F}_{N \, i}(\mu) = f(\zeta_i;\mu) = \sum_{q=1}^{Q_f} \Theta^q_f(\mu) \, f^q(\zeta_i), \label{eq:RB_thirteen-d}
\end{equation}
\end{subequations}


\noindent which can be written as

\begin{align}\label{eq:RB_fourteen}
\mathbf{M}_{N \, i,j}(\mu) = \sum_{q=1}^{Q_m} \Theta^q_m(\mu) \, \mathbf{M}_{N \, i,j}^q, &\qquad  \mathbf{C}_{N \, i,j}(\mu) = \sum_{q=1}^{Q_c} \Theta^q_c(\mu) \, \mathbf{C}_{N \, i,j}^q,  \nonumber \\
\mathbf{A}_{N \, i,j}(\mu) = \sum_{q=1}^{Q_a} \Theta^q_a(\mu) \, \mathbf{A}_{N \, i,j}^q, &\qquad  \mathbf{F}_{N \, i}(\mu) = \sum_{q=1}^{Q_f} \Theta^q_f(\mu) \, \mathbf{F}_{N \, i}^q,
\end{align}


\noindent where the \textit{parameter independent} quantities $\mathbf{M}_N^q$, $\mathbf{C}_N^q$, $\mathbf{A}_N^q \in \mathbb{R}^{N \times N}$, and $\mathbf{F}_N^q \in \mathbb{R}^N$ are given by


\begin{subequations}\label{eq:RB_fifteen}
\begin{equation}
\mathbf{M}_{N \, i,j}^q = m^q(\zeta_i,\zeta_j), \quad 1 \le i,j \le N_{\max}, \quad 1 \le q \le Q_m, \label{eq:RB_fifteen-a}
\end{equation}
\begin{equation}
\mathbf{C}_{N \, i,j}^q = c^q(\zeta_i,\zeta_j), \quad 1 \le i,j \le N_{\max}, \quad 1 \le q \le Q_c, \label{eq:RB_fifteen-b}
\end{equation}
\begin{equation}
\mathbf{A}_{N \, i,j}^q = a^q(\zeta_i,\zeta_j), \quad 1 \le i,j \le N_{\max}, \quad 1 \le q \le Q_a, \label{eq:RB_fifteen-c}
\end{equation}
\begin{equation}
\mathbf{F}_{N \, i}^q = f^q(\zeta_i), \quad 1 \le i \le N_{\max}, \quad 1 \le q \le Q_f, \label{eq:RB_fifteen-d}
\end{equation}
\end{subequations}

\noindent respectively.



The offline-online computational procedure is now described as follows. In the offline stage -- performed only once, we first implement the standard POD--Greedy algorithm \cite{Hoang2013}: we solve to find $Y^{\rm st}_N = \{ \zeta^{\rm st}_n, 1 \le n \le N_{\max} \}$; then compute and store the $\mu$-independent quantities in \eqref{eq:RB_fifteen} for the estimation of the RB solution and output\footnote{There are still several terms related to the computation of the dual norm of the residual $\|\mathcal{R}^{\rm st} (v;\mu,t^k)\|_{Y'}$, we do not show them here for simplicity. Detailed implementation of the standard POD--Greedy algorithm can be referred to, for instance, \cite{Hoang2013, Knezevic2011, Nguyen2009}.}. Once these RB solution and output are available, we can now implement the goal-oriented POD--Greedy algorithm. We consider each goal-oriented POD--Greedy iteration (Table \ref{tab1}) in more details. We first need to solve \eqref{eq:ProbStat_seven} for the ``true'' FE solutions; then compute the projection error in step (T\ref{tab1}d) and solve the POD/eigenvalue problem in step (T\ref{tab1}e). In addition, we have to compute $O(N^2 Q)$ $\mathcal{N}$-inner products in \eqref{eq:RB_fifteen} (we denote $Q=Q_m+Q_c+Q_a$). By approximating $s(\mu)$ via the enriched approximation $s^{\rm st}_{\tilde{N}}(\mu)$ in (T\ref{tab1}k) through the standard algorithm, we can now do an exhaustive yet cheap search over $\Xi_{\rm train}$ to look for the optimal $\mu$ in each Greedy iteration. In summary, the offline stage of the goal-oriented algorithm also includes the offline stage of the standard algorithm, and therefore it is more expensive than that of the standard algorithm. (We again emphasize that the standard algorithm is implemented only \textit{once}, and then the goal-oriented algorithm can be implemented \textit{many} times corresponding with various output functionals, respectively.)



The online stage of the goal-oriented algorithm is very similar to that of the standard algorithm \cite{Hoang2013}. In the online stage -- performed many times, for each new parameter $\mu$ -- we first assemble the RB matrices in \eqref{eq:RB_thirteen}, this requires $O(N^2 Q)$ operations. We then solve the RB governing equation \eqref{eq:RB_eleven}, the operation counts are $O(N^3+KN^2)$ as the RB matrices are generally full. Finally, we evaluate the displacement output $s_N(\mu)$ from \eqref{eq:RB_twelve} at the cost of $O(KN)$. For the error approximation (i.e., $\Delta_s(\mu) = s^{\rm st}_{\tilde{N}}(\mu) - s^{\rm go}_{N}(\mu)$), there is nothing more than computing one more output $s^{\rm st}_{\tilde{N}}(\mu)$ and then performing the associated subtraction; the cost is $O(\tilde{N}^3)$. Therefore, as required in the real-time context, the online complexity to compute the output and its associated error approximation are $O(N^3)+O(\tilde{N}^3)$ -- independent of $\mathcal{N}$; and since $N, \tilde{N} \ll \mathcal{N}$ we can expect significant computational saving in the online stage relative to the classical FE approach.





\section{Numerical examples}\label{sec_NumEx}

In this section, we will verify both POD--Greedy algorithms by investigating a simple 2D linear elastodynamic problem and a 3D dental implant model problem in the time domain. The details are described in the following.

% numerical example which is a three-dimensional dental implant model problem in the time domain. This model problem is similar to that in the work of Hoang et al. \cite{Hoang2013}. The details are described in the following.


\subsection{A 2D linear elastodynamic problem}\label{subsec_NumEx_2dProb}

\subsubsection{Finite element model and approximation}\label{subsubsec_NumEx_2dProb_ModelDescription}

%%%%%%%%%%%%%%%%% fig 2-1 %%%%%%%%%%%%%%%%%
\begin{figure}[h!]
 \centering
 \subfigure[]{ \includegraphics[height=4.9cm]{figs/2D_Model.eps} \label{fig21a}}
 \subfigure[]{ \includegraphics[height=4.8cm]{figs/2D_Mesh.eps}  \label{fig21b}}
 \caption{(a) The 2D model with its dimensionless parameters, and (b) its FE reference mesh.}
 \label{fig21}
\end{figure}


%%%%%%%%%%%%% fig 2-2 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_Load.eps}} \label{fig22a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_sFEM_Impulse.eps}} \label{fig22b} }%
\caption{(a) Time history of impulse loading, and (b) corresponding FEM output displacement in $x-$direction versus time with $\mu_{\rm test}=(1, 0.1)$. }%
\label{fig22}
\end{center}
\end{figure}



We consider a 2D plane strain model problem as in Fig.\ref{fig21a}. It is assumed that the model problem is scaled (or non-dimensionalized) from a real problem in practice and hence all the terms are dimensionless. The details of the non-dimensionalization is briefly discussed in Appendix \ref{sec_Appx}. The length and height of the model are $L=4$ and $H=1$, respectively. The model is composed of 2 subdomains $\Omega_1$ and $\Omega_2$ with two different materials: Young's moduli $E_1=1$ and $E_2 \in [0.1,10]$ and Poisson ratios $\nu_1 = \nu_2 = 0.3$, respectively. We assume Rayleigh damping for the model where the mass-proportional damping coefficients $\alpha_1 = \alpha_2 = 0$, and the stiffness-proportional damping coefficients $\beta_1 = \beta_2 \equiv \beta \in [0.05,0.5]$ such that $\textbf{C}_i = \beta_i \textbf{A}_i, \, 1 \le i \le 2$, where $\textbf{C}_i$ and $\textbf{A}_i$ are the FEM damping and stiffness matrices of each subdomain, respectively. Isotropic and homogenous material behavior is assumed for the model. We also note that the material mass densities will vanish from the weak form of the PDE due to the non-dimensionalisation (eq.\eqref{eq:Appx_five}, Appendix \ref{sec_Appx}). A surface traction which is opposite to the $x-$direction with the magnitude $\bar{\phi}=0.01$ is then applied to the right edge of the model ($\Gamma_{\rm l}$) as shown in Fig.\ref{fig21a}. As mentioned in Section \ref{subsec_ProbStat_ImpRes}, it is sufficient to perform all calculations for the impulse loading case since the solutions associated with other loading cases can be inferred by Duhamel's convolution principle. Time history of the impulse loading is also shown in Fig.\ref{fig22a}. Homogeneous Dirichlet boundary condition is enforced on the left edge ($\Gamma_D$), while zero initial conditions (i.e., zero displacement and velocity) are applied on the model. The output of interest is defined as the time integral of the average $x-$displacement along the right edge ($\Gamma_{\rm o} \equiv \Gamma_{\rm l}$) of the model.



Specifically, the explicit forms of all linear and bilinear forms associated with the 2D model problem are defined as follows

\begin{subequations}\label{eq:NumEx_one}
\begin{equation}
m(w,v) = \sum_{r=1}^{2} \sum_{i} \int_{\Omega_r} w_i v_i, \label{eq:NumEx_one-a}
\end{equation}
\begin{equation}
a(w,v;\mu) = \sum_{i,j,k,l} \int_{\Omega_1} \frac{\partial v_i}{\partial x_j} C_{ijkl}^{1} \frac{\partial w_k}{\partial x_l} + \mu_{1} \sum_{i,j,k,l} \int_{\Omega_2} \frac{\partial v_i}{\partial x_j} C_{ijkl}^2 \frac{\partial w_k}{\partial x_l}, \label{eq:NumEx_one-b}
\end{equation}
\begin{equation}
c(w,v;\mu) = \mu_2 \sum_{i,j,k,l} \int_{\Omega_1} \frac{\partial v_i}{\partial x_j} C_{ijkl}^1 \frac{\partial w_k}{\partial x_l} + \mu_{2}\mu_{1} \sum_{i,j,k,l} \int_{\Omega_2} \frac{\partial v_i}{\partial x_j} C_{ijkl}^2 \frac{\partial w_k}{\partial x_l}, \label{eq:NumEx_one-c}
\end{equation}
\begin{equation}
f(v) = \sum_{i} \int_{\Gamma_{\rm l}} v_i \, \bar{\phi}_i, \label{eq:NumEx_one-d}
\end{equation}
\begin{equation}
\ell(v) = \frac{1}{|\Gamma_{\rm o}|} \int_{\Gamma_{\rm o}} v_1, \label{eq:NumEx_one-e}
\end{equation}
\end{subequations}


for all $w, v \in Y$, $\mu \in \mathcal{D}$. Here, the parameter $\mu = (\mu_1,\mu_2) \equiv (E_2,\beta)$, where $E_2$ is the Young's modulus of the domain $\Omega_2$ and $\beta$ is the stiffness-proportional damping coefficient of both domains $\Omega_1$, $\Omega_2$. $C^r_{ijkl}$ is the constitutive elasticity tensor for isotropic materials and it is expressed in terms of the Young's modulus $E$ and Poisson's ratio $\nu$ of each region $\Omega_r, 1 \le r \le 2$, respectively. It is obvious from \eqref{eq:ProbStat_six} and \eqref{eq:NumEx_one} that the smooth functions $\Theta^1_a(\mu)=1$, $\Theta^2_a(\mu)=\mu_1$; $\Theta^1_c(\mu)=\mu_2$, $\Theta^2_c(\mu)=\mu_1 \mu_2$ depend on $\mu$ --- but the bilinear forms $a^1(w,v) = c^1(w,v) = \sum_{i,j,k,l} \int_{\Omega_1} \frac{\partial v_i}{\partial x_j} C_{ijkl}^{1} \frac{\partial w_k}{\partial x_l}$, and $a^2(w,v) = c^2(w,v) = \sum_{i,j,k,l} \int_{\Omega_2} \frac{\partial v_i}{\partial x_j} C_{ijkl}^2 \frac{\partial w_k}{\partial x_l}$ do \textit{not} depend on $\mu$.




The finite element mesh consist of 215 nodes and 376 linear triangular elements as shown in Fig.\ref{fig21b}. The FE space to approximate the 2D elastodynamic problem is of dimension $\mathcal{N}=416$. For time integration, $T=50$, $\Delta t = 0.2$, $K = \frac{T}{\Delta t}=250$. The parameter $\mu = (E_2,\beta) \in \mathcal{D}$, where the parameter domain $\mathcal{D} \equiv [0.1,10] \times [0.05,0.5] \subset \mathbb{R}^{P=2}$. The $\|\cdot\|_Y$ used in this work is defined as $\|w\|_Y^2 = a(w,w;\bar{\mu}) + m(w,w;\bar{\mu})$, where $\bar{\mu} = (1, 0.1)$; $Q_a=2$, $Q_c=2$. The entire work is performed using the software Matlab R2012b. We finally show in Fig.\ref{fig22b} the ``unit'' FEM output displacement (i.e., under the unit impulse load) in the $x-$direction versus time at $\mu_{\rm test}$, where $\mu_{\rm test} = \bar{\mu} = (1, 0.1)$.


%Fig.\ref{fig3} presents some FEM results corresponding with the test parameter $\mu_{\rm test}=(10 \times 10^6{\rm Pa}, 1 \times 10^{-5})$. In particular, Fig.\ref{fig3a} shows the ``unit'' FEM output displacement (i.e., under the unit impulse load) in the $x-$direction versus time at $\mu_{\rm test}$. In addition, Fig.\ref{fig3b} shows the FEM output displacements versus time under the ``arbitrary'' load in Fig.\ref{fig2b} by direct FEM computation and by Duhamel's convolution. It is observed that these two results match perfectly well with each other.




\subsubsection{Numerical results}\label{subsubsec_NumEx_2dProb_NumRes}


\paragraph{The impulse loading case}\label{paragraph_2dProb_ImpLoad}

For this 2D model problem, we aim to investigate the behavior of the goal-oriented algorithms with various choices of $\eta_T$ under the impulse loading regime. To start, a training sample set $\Xi_{\rm train}$ is created by a equidistant distribution over $\mathcal{D}$ with $n_{\rm train}(=30 \times 30) = 900$ samples. Note that we use $M=1$ and $N^{\rm go}_{\max}=60$ (as in Table \ref{tab1}) to terminate the iteration procedures. In the remaining sections, beside the standard and goal-oriented algorithms ($\eta_T$), we also show the results of the goal-oriented algorithm ($\tilde{N}=2N$) for comparison purpose.


\begin{enumerate}[a)]

  \item First, we implement the standard POD--Greedy algorithm (i.e., the right column of Table \ref{tab1}) \cite{Hoang2013}. The results are presented in Fig.\ref{fig23}: Fig.\ref{fig23a} shows the maximum error indicator $\Delta^{\max,{\rm rel}}_{u} = \max \limits_{\mu \in \Xi_{\rm train}} \left\{ \frac{\Delta_u(\mu)}{\sqrt{\sum_{k=1}^{K}\|u^{\rm st}_N(\mu,t^k)\|^2_Y }} \right\}$ as a function of $N$; and Fig.\ref{fig23b} shows the distribution of the sample set $S^{\rm st}_{*}$, respectively. As observe from Fig.\ref{fig23a}, the max error indicator $\Delta^{\max,{\rm rel}}_{u}$ could not decrease further when $N^{\rm st} \ge 140$, hence we set $N^{\rm st}_{\max} = 150$ for the implementation of the GO algorithms afterwards. (Of course, $N^{\rm st}_{\max}$ is still able to be enriched automatically whenever $\tilde{N} > N^{\rm st}_{\max}$ as described in Algorithm 2 of the ``cross-validation'' process, Section \ref{subsubsec_RB_GOPODGreedy_GO}.) For comparison purpose, the results associated with the first $N^{\rm st}_{\max}=60$ basis functions were plotted using different markers in Fig.\ref{fig23}, respectively.


%%%%%%%%%%%%% fig 2-3 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_StdPODGreedy_DeltaU.eps}} \label{fig23a}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_StdPODGreedy_SamplePts.eps}} \label{fig23b} }%
\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_u$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the standard POD--Greedy algorithm ($N^{\rm st}_{\max}=150$). Different markers were used for the first $N^{\rm st}_{\max}=60$ basis functions. }%
\label{fig23}
\end{center}
\end{figure}




  \item Once the sets $S^{\rm st}_{*}$ and $Y^{\rm st}_{N}$ are available, the term $s^{\rm st}_{\tilde{N}}(\mu)$ (of $\Delta_s(\mu)$ in (T\ref{tab1}k)) is now computable; and hence, it is possible to implement the goal-oriented POD--Greedy algorithms (i.e., the left column of Table \ref{tab1}). For comparison, we show the results of the goal-oriented algorithms using ($\tilde{N}=2N$), ($\tilde{N},\eta_T=0.8$), ($\tilde{N},\eta_T=0.5$) and ($\tilde{N},\eta_T=0.1$) in Fig.\ref{fig24}, respectively. Note that there is no ``cross-validation'' process for the first algorithm as we just fix $\tilde{N}=2N$ for the computation of $s^{\rm st}_{\tilde{N}}(\mu)$ in $\Delta_s(\mu)$. The latter algorithms use the proposed ``cross-validation'' process to choose adaptively $\tilde{N}$ with the settings: $\left| \Xi^{\rm st}_1 \right|=10$ and $\Delta n_{\rm sample}=10$ corresponding with $\eta_T=0.8$, 0.5 and 0.1, respectively. Fig.\ref{fig24a} presents the maximum error indicators $\Delta^{\max,{\rm rel}}_{s} = \max \limits_{\mu \in \Xi_{\rm train}} \left\{ \left| \frac{ \Delta_s(\mu) }{s^{\rm st}_{\tilde{N}}(\mu)} \right| \right\}$ as a function of $N$; and Fig.\ref{fig24b} shows the distribution of the sample sets $S^{\rm go}_{*}$. The results of the ``cross-validation'' processes ($\tilde{N},\eta_T=0.8$), ($\tilde{N},\eta_T=0.5$) and ($\tilde{N},\eta_T=0.1$) are also presented in Fig.\ref{fig25}. In particular, Fig.\ref{fig25a} presents all the $\tilde{N}$ found adaptively for each $N$; Fig.\ref{fig25b} shows the size of $\Xi^{\rm st}_n$ as a function of $N$; Fig.\ref{fig25c} and Fig.\ref{fig25d} illustrate the maximum and minimum effectivities in \eqref{eq:RB_five_tmp} over $\Xi^{\rm st}_n$, respectively. Fig.\ref{fig25} confirms the behavior of these GO algorithms: as $\eta_T$ increases close to 1, the range $[\eta_T, 2-\eta_T]$ shrinks around 1 (i.e., better test effectivities), hence the algorithm would need larger $\tilde{N}$ to hold \eqref{eq:RB_five_tmp} true.

% $\tilde{N}$ will be increased correspondingly and hence better true error and error approximation.

% It is observed from Fig.\ref{fig24a} that the GO ($\tilde{N}, \eta_T=0.8$) algorithm is just a bit better than the GO ($\tilde{N}=2N$) one regarding maximum error indicators. However, Fig.\ref{fig24b} shows that these two algorithms are different from each other by 15/60 (=1/4) sample points. This is an important difference which will be reflected clearly in the effectivity estimation by the two algorithms afterwards.



%%%%%%%%%%%%% fig 2-4 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_GO_ErrInd_All.eps}} \label{fig24a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_GO_SamplingPts_All.eps}} \label{fig24b} }
\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_s$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the goal-oriented POD--Greedy algorithms ($\tilde{N}=2N$), ($\tilde{N},\eta_T=0.8$), ($\tilde{N},\eta_T=0.5$) and ($\tilde{N},\eta_T=0.1$).}%
\label{fig24}
\end{center}
\end{figure}


%%%%%%%%%%%%% fig 2-5 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_CVprocess_Ntilde.eps}} \label{fig25a} }
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_CVprocess_Xi_st.eps}} \label{fig25b} }
\\
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_CVprocess_maxEff.eps}} \label{fig25c} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_CVprocess_minEff.eps}} \label{fig25d} }
\caption{(a) $\tilde{N}$, (b) the size $\left| \Xi^{\rm st}_n \right|$, (c) max, and (d) min effectivities $\left| \frac{\Delta_s(\mu)}{s(\mu)-s^{\rm go}_{N}(\mu)} \right|$ in \eqref{eq:RB_five_tmp} as functions of $N$ of the GO algorithms ($\tilde{N}, \eta_T=0.8$), ($\tilde{N}, \eta_T=0.5$) and ($\tilde{N}, \eta_T=0.1$). }%
\label{fig25}
\end{center}
\end{figure}




  \item In order to evaluate the performance of the algorithms, we create a test sample set $\Xi_{\rm test} \subset \mathcal{D}$ which is a coarse subset of $\mathcal{D}$; then compute and compare mutually their RB true errors. This set has $n_{\rm test} (=10 \times 10) = 100$ sample points distributed equidistantly. We show, as functions of $N$: $e^{\max}_{u} = \max \limits_{\mu \in \Xi_{\rm test}} e_u(\mu)$ (defined in \eqref{eq:RB_six}) and $e^{\max}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} e_s(\mu)$ (defined in \eqref{eq:RB_seven}) in Fig.\ref{fig26}, respectively. Fig.\ref{fig26a} shows quite similar true solution errors by the algorithms; while from Fig.\ref{fig26b}, the goal-oriented algorithms give true output errors smaller than that computed by the standard algorithm. The differences in the RB true output errors by the standard versus GO algorithms can be up to 10 times. Therefore, we can conclude that the GO algorithms (with cross-validation processes) give RB true output errors smaller than that provided by the standard algorithm (for the same $N$ basis functions under consideration).

% As expected from Fig.\ref{fig26b}, the GO ($\tilde{N}, \eta_T=0.8$) algorithm is a bit better than the GO ($\tilde{N}=2N$) algorithm because the former uses more RB basis functions ($\tilde{N}$) than the latter with the same $N$ (see Fig.\ref{fig25a}).


%%%%%%%%%%%%% fig 2-6 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_TrueErr_Sol.eps}} \label{fig26a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_TrueErr_Output.eps}} \label{fig26b} }%
\caption{Comparison of maximum (relative) RB true errors by standard and all goal-oriented POD--Greedy algorithms over $\Xi_{\rm test}$: (a) solution and (b) output.}%
\label{fig26}
\end{center}
\end{figure}




  \item Lastly, we evaluate the performance of the sampling algorithms via their (practical) output error approximations. For comparison, the output error approximation \eqref{eq:RB_nine} by the standard algorithm with ($\tilde{N},\eta_T=0.8$) is used. (Namely, the pairs $(N,\tilde{N})$ in \eqref{eq:RB_nine} is taken following and after implementing the GO ($\tilde{N},\eta_T=0.8$) algorithm\footnote{The aim of this special test is that for the same online computational effort $O(N^3)+O(\tilde{N}^3)$, we want to compare the error approximations and associated effectivities by the standard vs. the GO ($\tilde{N},\eta_T=0.8$) algorithms.}.) In particular, Fig.\ref{fig27b} shows the graphs of $\tilde{\Delta}^{\rm max,rel}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} \left| \frac{\Delta_s(\mu)}{s(\mu)} \right|$ (defined in \eqref{eq:RB_eight}, \eqref{eq:RB_nine}) as functions of $N$ by the algorithms: standard ($\tilde{N},\eta_T=0.8$), GO ($\tilde{N}=2N$), GO ($\tilde{N},\eta_T=0.8$), GO ($\tilde{N},\eta_T=0.5$) and GO ($\tilde{N},\eta_T=0.1$), respectively. Fig.\ref{fig27a} is just a repetition of Fig.\ref{fig26b} -- which shows the max relative RB true errors in order to compare with the max output error approximations in Fig.\ref{fig27b}. The associated maximum/minimum effectivities of this error approximation (i.e., $\eta^{\max}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} \eta_s(\mu)$, $\eta^{\min}_{s} = \min \limits_{\mu \in \Xi_{\rm test}} \eta_s(\mu)$, defined in \eqref{eq:RB_eight}, \eqref{eq:RB_nine}) are also shown in Fig.\ref{fig27c} and Fig.\ref{fig27d}, respectively. As observed from Fig.\ref{fig27a} and Fig.\ref{fig27b}, the GO algorithms are superior to the standard algorithm in both RB true error and output error approximation. Regarding effectivities, Fig.\ref{fig27c} and Fig.\ref{fig27d} show that the GO ($\tilde{N},\eta_T=0.8$) algorithm gives the best effectivities among all the algorithms. These results confirm the good performance of the ``cross-validation'' process: the GO ($\tilde{N},\eta_T=0.8$) algorithm not only decreases the RB true error compared with the standard algorithm, but it also improves significantly the output effectivity estimation compared with other GO algorithms (especially the GO ($\tilde{N}=2N$) one).


% quite similar to the standard algorithm and that both of them completely beat the GO ($\tilde{N}=2N$) one. These results confirm the good performance of the ``cross-validation'' process: the GO ($\tilde{N},\eta_T=0.8$) algorithm decreases the RB true error compared with the standard algorithm, and it also improves significantly the output effectivity estimation compared with the GO ($\tilde{N}=2N$) algorithm. It also seems that (without proof) $\Delta^{\rm go}_s(\mu)$ by the GO ($\tilde{N}=2N$) algorithm does not grow exponentially with time.



%  As observed, we see that the proposed error approximation $\Delta_s(\mu)$ is very close to the true output error (Fig.\ref{fig13a}) and that their effectivities are really good (Fig.\ref{fig13b}). In fact, the Fig.\ref{fig12b} and Fig.\ref{fig13b} confirm the good performance of the ``cross-validation'' process: the GO ($\tilde{N},\eta_T=0.8$) algorithm improves significantly the output effectivity estimation compared with the GO ($\tilde{N}=2N$) one. It also seems that (without proof) $\Delta_s(\mu)$ does not grow exponentially with time.

%      that its maximum effectivity is about $O(10)$ which is very good for a 3D time-dependent dynamic problem (Fig.\ref{fig13b}).

%%%%%%%%%%%%% fig 2-7 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_Impulse__TrueErrs.eps}} \label{fig27a} }
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_Impulse__ErrApprx.eps}} \label{fig27b} } \\
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_Impulse__MaxEff.eps}} \label{fig27c}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/2D_Impulse__MinEff.eps}} \label{fig27d}}
\caption{ ``Impulse load'' case: comparison of (a) max relative RB true errors, (b) max relative error approximations, (c) corresponding max effectivities and (d) corresponding min effectivities of all sampling algorithms over $\Xi_{\rm test}$. }
\label{fig27}
\end{center}
\end{figure}




\end{enumerate}






\subsection{A 3D dental implant model problem}\label{subsec_NumEx_3dDental}

\subsubsection{Finite element model and approximation}\label{subsubsec_NumEx_3dDental_ModelDescription}

%%%%%%%%%%%%%%%%% fig 1 %%%%%%%%%%%%%%%%%
\begin{figure}[h!]
 \centering
 \subfigure[]{ \includegraphics[height=5.5cm]{figs/3D_FEMmodel.eps} \label{fig1a}}
 \subfigure[]{ \includegraphics[height=6.3cm]{figs/3D_MeshModel.eps} \label{fig1b}}
 \caption{(a) The 3D simplified FEM model with sectional view, and (b) meshing in ABAQUS.}
 \label{fig1}
\end{figure}



%%%%%%%%%%%%%%%%% fig 2 %%%%%%%%%%%%%%%%%
\begin{figure}[h!]
 \centering
 \subfigure[]{ \includegraphics[height=6cm]{figs/3D_BoundaryCondition.eps} \label{fig2a}}
 \subfigure[]{ \includegraphics[height=6cm]{figs/3D_Load.eps} \label{fig2b}}
 \caption{(a) Output area, applied load {\bf F} and boundary condition, and (b) time history of an impulse and an ``arbitrary'' load.}
 \label{fig2}
\end{figure}


We now consider a simplified 3D dental implant-bone model in Fig.\ref{fig1a}. The geometry of the simplified dental implant-bone model is constructed by using the software SolidWorks 2010. The physical domain $\Omega$ consists of five regions: the outermost cortical bone $\Omega_1$, the cancellous bone $\Omega_2$, the interfacial tissue $\Omega_3$, the dental implant $\Omega_4$ and the stainless steel screw $\Omega_5$. The 3D simplified model is then meshed using the software ABAQUS/CAE version 6.10-1 (Fig.\ref{fig1b}). A dynamic force opposite to the $x-$direction is then applied to a prescribed area on the body of the screw as shown in Fig.\ref{fig2a}. As mentioned in Sections \ref{subsec_ProbStat_ImpRes} and \ref{subsubsec_NumEx_2dProb_ModelDescription}, all computations and simulations will be performed for the impulse loading case thanks to the Duhamel's convolution principle \eqref{eq:ProbStat_ten}. Fig.\ref{fig2b} shows the time history of both loading cases used in this work: the impulse load and an ``arbitrary'' load. The Dirichlet boundary condition $(\partial \Omega^D)$ is specified in the bottom-half of the simplified model as illustrated in Fig.\ref{fig2a}; while zero initial conditions (i.e., zero displacement and velocity) are applied on the model. The output of interest is defined as the average $x-$displacement responses of a prescribed area on the head of the screw (Fig.\ref{fig2a}). The finite element mesh consists of 9479 nodes and 50388 four-node tetrahedral solid elements. The coinciding nodes of the contact surfaces between different regions (the regions $\Omega_1$, $\Omega_2$, $\Omega_3$, $\Omega_4$, $\Omega_5$) are assumed to be rigidly fixed, i.e. the displacements in the $x-$, $y-$ and $z-$directions are all set to be the same for the same coinciding nodes.



We assume that the regions $\Omega_i, 1 \leq i \le 5$, of the simplified model are homogeneous and isotropic. The material properties: the Young's moduli, Poisson's ratios and densities of these regions are presented in Table \ref{tab2} \cite{Wang2010}. As similar to \cite{Hoang2013}, we still use Rayleigh damping with stiffness-proportional damping coefficient $\beta_i$, $1 \le i \le 5$ (Table \ref{tab2}) such that $\textbf{C}_i = \beta_i \textbf{A}_i, \, 1 \le i \le 5$, where $\textbf{C}_i$ and $\textbf{A}_i$ are the FEM damping and stiffness matrices of each region, respectively. We also note in Table \ref{tab2} that ($E$,$\beta$) -- the Young's modulus and Rayleigh damping coefficient associated with the interfacial tissue are our sole parameters.


%%%%%%%%%%%%% table 2 %%%%%%%%%%%%%
\begin{table}[h!]
\begin{center}
  {\begin{tabular}{|c|l|r|l|r|r|}
  \hline \rule{0pt}{2.5ex}
 Domain & Layers & E (Pa) & $\nu$ & $\rho (\rm g/mm^3)$ & $\beta$ \\ [0.5ex]
  \hline \rule{0pt}{2.5ex}
 $\Omega_1$ & Cortical bone         & $2.3162 \times 10^{10}$ & 0.371  & $1.8601 \times 10^{-3}$ & $3.38 \times 10^{-6}$ \\
 $\Omega_2$ & Cancellous bone       & $8.2345 \times 10^{8}$  & 0.3136 & $7.1195 \times 10^{-4}$ & $6.76 \times 10^{-6}$ \\
 $\Omega_3$ &  Tissue                &      E                  & 0.3155 & $1.055  \times 10^{-3}$ &         $\beta$ \\
 $\Omega_4$ &  Titan implant         & $1.05  \times 10^{11}$  & 0.32   & $4.52   \times 10^{-3}$ & $5.1791 \times 10^{-10}$ \\
 $\Omega_5$ &  Stainless steel screw & $1.93  \times 10^{11}$  & 0.305  & $8.027  \times 10^{-3}$ & $2.5685 \times 10^{-8}$ \\ [0.3ex]
   \hline
\end{tabular}
\caption{Material properties of the dental implant-bone structure.}
\label{tab2}}
%\tabnote{$^{\rm a}$This is not as accurate, owing to numerical
%error.} \tabnote{$^{\rm b}$An example table footnote to show the
%text turning over when a long footnote is inserted.}
\end{center}
\end{table}



With respect to our particular dental implant problem, the actual integral forms of the linear and bilinear forms are defined as:


\begin{subequations}\label{eq:NumEx_two}
\begin{equation}
m(w,v) = \sum_{r=1}^{5} \sum_{i} \int_{\Omega_r} \rho_r w_i v_i, \label{eq:NumEx_two-a}
\end{equation}
\begin{equation}
a(w,v;\mu) = \sum_{r=1, r \neq 3}^{5} \sum_{i,j,k,l} \int_{\Omega_r} \frac{\partial v_i}{\partial x_j} C_{ijkl}^{r} \frac{\partial w_k}{\partial x_l} + \mu_{1} \sum_{i,j,k,l} \int_{\Omega_3} \frac{\partial v_i}{\partial x_j} C_{ijkl}^3 \frac{\partial w_k}{\partial x_l}, \label{eq:NumEx_two-b}
\end{equation}
\begin{equation}
c(w,v;\mu) = \sum_{r=1, r \neq 3}^{5} \beta_r \sum_{i,j,k,l} \int_{\Omega_r} \frac{\partial v_i}{\partial x_j} C_{ijkl}^r \frac{\partial w_k}{\partial x_l} + \mu_{2}\mu_{1} \sum_{i,j,k,l} \int_{\Omega_3} \frac{\partial v_i}{\partial x_j} C_{ijkl}^3 \frac{\partial w_k}{\partial x_l}, \label{eq:NumEx_two-c}
\end{equation}
\begin{equation}
f(v) = \sum_{i} \int_{\Gamma_{\rm l}} v_i \, \bar{\phi}_i, \label{eq:NumEx_two-d}
\end{equation}
\begin{equation}
\ell(v) = \frac{1}{|\Gamma_{\rm o}|} \int_{\Gamma_{\rm o}} v_1, \label{eq:NumEx_two-e}
\end{equation}
\end{subequations}


\noindent for all $w,v \in Y$, $\mu \in \mathcal{D}$. Here, the parameter $\mu=(\mu_1, \mu_2) \equiv (E,\beta)$ belongs to the region $\Omega_3$. $C_{ijkl}^r$ is the constitutive elasticity tensor for isotropic materials and it is expressed in terms of the Young's modulus $E$ and Poisson's ratio $\nu$ of each region $\Omega_r, 1 \leq r \leq 5$, respectively. $\Gamma_{\rm l}$ is the prescribed loading area (surface traction) and $\Gamma_{\rm o}$ is the prescribed output area as shown in Fig.\ref{fig2a}, respectively. It is obvious from \eqref{eq:ProbStat_six} and \eqref{eq:NumEx_two} that the smooth functions $\Theta^1_a(\mu) = 1$, $\Theta^2_a(\mu) = \mu_1$; $\Theta^1_c(\mu) = 1$, $\Theta^2_c(\mu) = \mu_1 \mu_2$ depend on $\mu$ --- but the bilinear forms $a^1(w,v) = \sum_{r=1, r \neq 3}^{5} \sum_{i,j,k,l} \int_{\Omega_r} \frac{\partial v_i}{\partial x_j} C_{ijkl}^{r} \frac{\partial w_k}{\partial x_l}$, $a^2(w,v) = c^2(w,v) = \sum_{i,j,k,l} \int_{\Omega_3} \frac{\partial v_i}{\partial x_j} C_{ijkl}^3 \frac{\partial w_k}{\partial x_l}$ and $c^1(w,v) = \sum_{r=1, r \neq 3}^{5} \beta_r \sum_{i,j,k,l} \int_{\Omega_r} \frac{\partial v_i}{\partial x_j} C_{ijkl}^{r} \frac{\partial w_k}{\partial x_l}$ do \textit{not} depend on $\mu$.




%%%%%%%%%%%%% fig 3 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_sFEM_Impulse.eps}} \label{fig3a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_sFEM_Duhamel.eps}} \label{fig3b} }%
\caption{(a) ``Unit'' FEM output displacement in $x-$direction versus time, and (b) comparison of ``arbitrary'' FEM output displacements (under the applied load in Fig.\ref{fig2b}) computed by Duhamel's convolution and direct computation, with $\mu_{\rm test}=(10 \times 10^6{\rm Pa}, 1 \times 10^{-5})$. }%
\label{fig3}
\end{center}
\end{figure}


%%%%%%%%%%%%%% fig 4 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig4a.eps}} \label{fig4a} }\\%
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig4b.eps}} \label{fig4b} }%
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig4c.eps}} \label{fig4c} }%
%\caption{Comparison of ``arbitrary'' FEM output displacements (with the applied load in Fig.\ref{fig2b}) computed by Duhamel's convolution and direct computation with respect time in the $x-$ (a), $y-$ (b), and $z-$direction (c) with $\mu_{\rm test}=(10 \times 10^6{\rm Pa}, 1 \times 10^{-5})$.}%
%\label{fig4}
%\end{center}
%\end{figure}


%%%%%%%%%%%%%% fig 5 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\resizebox*{7cm}{!}{\includegraphics{figs/fig5.eps}}
%\caption{The training set $\Xi_{\rm train}$ with $n_{\rm train}=900$ samples.}
%\label{fig5}
%\end{center}
%\end{figure}


The FE space to approximate the 3D dental implant-bone problem is of dimension $\mathcal{N}=26343$. For time integration, $T=1 \times 10^{-3}$s, $\Delta t = 2 \times 10^{-6}$s, $K=\frac{T}{\Delta t}=500$. The input parameter $\mu \equiv (E,\beta) \in \mathcal{D}$, where the parameter domain $\mathcal{D} \equiv [1 \times 10^6, 25 \times 10^6]{\rm Pa} \times [5 \times 10^{-6}, 5 \times 10^{-5}] \subset \mathbb{R}^{P=2}$. (Note that the range of $E$ of this parameter domain is nearly two times larger than that of \cite{Hoang2013}.) The $\|\cdot\|_{Y}$ norm used in this work is defined as $\|w\|_Y^2 = a(w,w;\bar{\mu}) + m(w,w;\bar{\mu})$, where $\bar{\mu} = (13 \times 10^6 {\rm Pa}, 2.75 \times 10^{-5})$ is the arithmetic average of $\mu$ in $\mathcal{D}$; $Q_a=2$, $Q_c=2$. The entire work is performed using the software Matlab R2012b. Fig.\ref{fig3} presents some FEM results corresponding with the test parameter $\mu_{\rm test}=(10 \times 10^6{\rm Pa}, 1 \times 10^{-5})$. In particular, Fig.\ref{fig3a} shows the ``unit'' FEM output displacement (i.e., under the unit impulse load) in the $x-$direction versus time at $\mu_{\rm test}$. In addition, Fig.\ref{fig3b} shows the FEM output displacements versus time under the ``arbitrary'' load (in Fig.\ref{fig2b}) by direct FEM computation and by Duhamel's convolution. It is observed that these two results match perfectly well with each other.


% The entire work is performed using the software Matlab R2012b.

% To verify our computational code (performed in Matlab R2012b), we first compare the FEM outputs computed by ABAQUS and by our code with the test parameter $\mu_{\rm test}=(10 \times 10^6{\rm Pa}, 1 \times 10^{-5})$. Fig.\ref{fig3} shows the ``unit'' output displacements (i.e., under the unit impulse load) in the $x-$, $y-$ and $z-$directions versus time at $\mu_{\rm test}$ via ABAQUS and our code, respectively. Fig.\ref{fig3} demonstrates that the FEM results by our code match very well with the results computed by ABAQUS. Next, we show in Fig.\ref{fig4} the FEM output displacements versus time under the ``arbitrary'' load in Fig.\ref{fig2b} by direct FEM computation and by Duhamel's convolution. It is observed that these two results match perfectly well with each other.





\subsubsection{Numerical results}\label{subsubsec_NumEx_3dDental_NumRes}

%\subsubsection{The impulse loading case}\label{subsubsec_NumEx_3dDental_ImpLoad}

\paragraph{The impulse loading case}\label{paragraph_3dDental_ImpLoad}


We now discuss the POD--Greedy algorithms of interest. A training sample set $\Xi_{\rm train}$ is created by a equidistant distribution over $\mathcal{D}$ with $n_{\rm train}(=30 \times 30) = 900$ samples. We still use $M=1$ and $N^{\rm go}_{\max}=60$ (Table \ref{tab1}) to terminate the iteration procedures. In the remaining sections, we only show the results of 3 typical algorithms which are the standard, the GO ($\tilde{N}, \eta_T=0.8$) and the GO ($\tilde{N}=2N$) algorithms for comparison.


\begin{enumerate}[a)]

  \item As similar to Section \ref{subsubsec_NumEx_2dProb_NumRes} above, the standard POD--Greedy algorithm is implemented first. The results are presented in Fig.\ref{fig6}: Fig.\ref{fig6a} shows the maximum error indicator $\Delta^{\max,{\rm rel}}_{u} = \max \limits_{\mu \in \Xi_{\rm train}} \left\{ \frac{\Delta_u(\mu)}{\sqrt{\sum_{k=1}^{K}\|u^{\rm st}_N(\mu,t^k)\|^2_Y }} \right\}$ as a function of $N$; and Fig.\ref{fig6b} shows the distribution of the sample set $S^{\rm st}_{*}$, respectively. $N^{\rm st}_{\max}$ is set to be 200 since the goal-oriented algorithms may need a quite large $N^{\rm st}_{\max}$ afterwards. For comparison purpose, the results associated with the first $N^{\rm st}_{\max}=60$ basis functions were also plotted using different markers as in Fig.\ref{fig6}, respectively.

% (Of course, $N^{\rm st}_{\max}$ is still able to be enriched automatically whenever $\tilde{N} > N^{\rm st}_{\max}$ as described in Algorithm 2 of the ``cross-validation'' process, Section \ref{subsubsec_RB_GOPODGreedy_GO}.)

%%%%%%%%%%%%% fig 6 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_StdPODGreedy_DeltaU.eps}} \label{fig6a}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_StdPODGreedy_SamplePts.eps}} \label{fig6b} }%
\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_u$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the standard POD--Greedy algorithm ($N^{\rm st}_{\max}=200$). Different markers were used for the first $N^{\rm st}_{\max}=60$ basis functions. }%
\label{fig6}
\end{center}
\end{figure}



%%%%%%%%%%%%%% fig 7 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig7a.eps}} \label{fig7a} }%
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig7b.eps}} \label{fig7b} }%
%\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_u$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the standard POD--Greedy algorithm ($N^{\rm st}_{\max}=60$).}%
%\label{fig7}
%\end{center}
%\end{figure}




  \item Now, the goal-oriented POD--Greedy algorithms can be implemented since the term $s^{\rm st}_{\tilde{N}}(\mu)$ is computable from the availability of the sets $S^{\rm st}_{*}$ and $Y^{\rm st}_{N}$. For comparison, we show the results of the goal-oriented algorithms using ($\tilde{N}=2N$) versus ($\tilde{N},\eta_T=0.8$) in Fig.\ref{fig8}. The settings $\left| \Xi^{\rm st}_1 \right|=10$ and $\Delta n_{\rm sample}=10$ are used for this GO ($\tilde{N},\eta_T=0.8$) algorithm. Fig.\ref{fig8a} presents the maximum error indicators $\Delta^{\max,{\rm rel}}_{s} = \max \limits_{\mu \in \Xi_{\rm train}} \left\{ \left| \frac{ \Delta_s(\mu) }{s^{\rm st}_{\tilde{N}}(\mu)} \right| \right\}$ as a function of $N$; and Fig.\ref{fig8b} shows the distribution of the sample sets $S^{\rm go}_{*}$, respectively. It is observed from Fig.\ref{fig8a} that the GO ($\tilde{N}, \eta_T=0.8$) algorithm is just a bit better than the GO ($\tilde{N}=2N$) one regarding maximum error indicators. However, Fig.\ref{fig8b} shows that these two algorithms are different from each other by 15/60 (=1/4) sample points. This is an important difference which will be reflected clearly in the effectivity estimation by the two algorithms afterwards. The results of the ``cross-validation'' process ($\tilde{N},\eta_T=0.8$) are also presented in Fig.\ref{fig10}. In particular, Fig.\ref{fig10a} presents all the $\tilde{N}$ found adaptively for each $N$; Fig.\ref{fig10b} illustrates the maximum and minimum effectivities in \eqref{eq:RB_five_tmp} over $\Xi^{\rm st}_n$; and Fig.\ref{fig10c} shows the size of $\Xi^{\rm st}_n$ as a function of $N$, respectively.

% Note that there is no ``cross-validation'' process for the former as we just fix $\tilde{N}=2N$ for the computation of $s^{\rm st}_{\tilde{N}}(\mu)$ in $\Delta_s(\mu)$. The latter uses the proposed ``cross-validation'' process to choose adaptively $\tilde{N}$ with the settings: $\eta_T=0.8$, $\left| \Xi^{\rm st}_1 \right|=10$ and $\Delta n_{\rm sample}=10$.


% \hl{Include a discussion on} $\left| \Xi^{\rm st}_{1} \right|$ and $\Delta n_{\rm sample}$ here!



%      For comparison, we first show the results of the goal-oriented algorithm with $\tilde{N}=2N$ in Fig.\ref{fig8}. (Note that there is no ``cross-validation'' process with this algorithm; we just fix $\tilde{N}=2N$ for the computation of $s^{\rm st}_{\tilde{N}}(\mu)$ in $\Delta_s(\mu)$.)

%        Then we show the results using the goal-oriented algorithm with our proposed ``cross-validation'' process to choose $\tilde{N}$ with $\eta_T=0.8$ in Fig.\ref{fig9}, respectively.

%      Consider the term $s^{\rm st}_{2N}(\mu)$, the subscript ``$2N$'' means that we approximate the FE output $s(\mu)$ by the (standard algorithm) enriched RB approximation $s^{\rm st}_{2N}(\mu)$ with $2N$ basis functions corresponding to the current $N^{\rm th}$-Greedy iteration of the goal-oriented algorithm: $N^{\rm st} = 2 N^{\rm go}$. Since $N^{\rm st}_{\max} = 120$, then we set
%
%      \begin{equation}\label{eq:NumEx_two}
%      N^{\rm st} =
%          \begin{cases}
%            2N^{\rm go} & \text{if } \quad 1 \le N^{\rm go} \le 60, \\
%            120         & \text{if } \quad 60 < N^{\rm go} \le 120,
%          \end{cases}
%      \end{equation}
%
%      to compute the term $s^{\rm st}_{2N}(\mu)$. Implementing the goal-oriented algorithm, the results are presented in Fig.\ref{fig7}:


%%%%%%%%%%%%% fig 8 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_ErrApprx_GO__2N_Ntilde.eps}} \label{fig8a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_SamplingPts__2N_Ntilde.eps}} \label{fig8b} }
\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_s$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the goal-oriented POD--Greedy algorithms ($\tilde{N}=2N$) v.s. ($\tilde{N},\eta_T=0.8$).}%
\label{fig8}
\end{center}
\end{figure}


%%%%%%%%%%%%%% fig 9 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig9a.eps}} \label{fig9a} }%
%\subfigure[]{
%\resizebox*{7cm}{!}{\includegraphics{figs/fig9b.eps}} \label{fig9b} }%
%\caption{(a) Maximum of error indicator $\Delta^{\max,{\rm rel}}_s$ over $\Xi_{\rm train}$ as a function of $N$, and (b) distribution of sampling points by the goal-oriented POD--Greedy algorithm ($\tilde{N}, \eta_T=0.8$).}%
%\label{fig9}
%\end{center}
%\end{figure}


%%%%%%%%%%%%% fig 10 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_CVprocess_Ntilde.eps}} \label{fig10a} } \\
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_CVprocess_maxminEff.eps}} \label{fig10b} }
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_CVprocess_size_Xi_st.eps}} \label{fig10c} }%
\caption{(a) $\tilde{N}$, (b) max/min effectivities $\left| \frac{\Delta_s(\mu)}{s(\mu)-s^{\rm go}_{N}(\mu)} \right|$ in \eqref{eq:RB_five_tmp}, $\forall \mu \in \Xi^{\rm st}_n$, and (c) the size $\left| \Xi^{\rm st}_n \right|$ as functions of $N$ of the algorithm GO ($\tilde{N}, \eta_T=0.8$).}%
\label{fig10}
\end{center}
\end{figure}



%The set is defined as follows: $\Xi_{\rm test} \equiv [1 \times 10^6, 25 \times 10^6]{\rm Pa} \times [5 \times 10^{-6}, 5 \times 10^{-5}]$ where $n_{\rm test}=100$ sample points distributed uniformly.


  \item In order to evaluate the performance of the algorithms, we create a test sample set $\Xi_{\rm test} \subset \mathcal{D}$ which is a coarse subset of $\mathcal{D}$; then compute and compare mutually their RB true errors. This set has $n_{\rm test} (= 10 \times 10) = 100$ sample points distributed equidistantly. We show, as functions of $N$: $e^{\max}_{u} = \max \limits_{\mu \in \Xi_{\rm test}} e_u(\mu)$ (defined in \eqref{eq:RB_six}) and $e^{\max}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} e_s(\mu)$ (defined in \eqref{eq:RB_seven}) in Fig.\ref{fig11}, respectively. We observe again from Fig.\ref{fig11b} that the goal-oriented algorithms completely prevail over the the standard algorithm as regards true output errors. The differences in the RB true output errors by the standard versus GO algorithms can be up to 10 times. As expected from Fig.\ref{fig11b}, the GO ($\tilde{N}, \eta_T=0.8$) algorithm is a bit better than the GO ($\tilde{N}=2N$) algorithm because the former uses more RB basis functions ($\tilde{N}$) than the latter with the same $N$ (see Fig.\ref{fig10a}).


%%%%%%%%%%%%% fig 11 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_TrueErr_Sol.eps}} \label{fig11a} }%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_TrueErr_Output.eps}} \label{fig11b} }%
\caption{Comparison of maximum (relative) RB true errors by standard and all goal-oriented POD--Greedy algorithms over $\Xi_{\rm test}$: (a) solution and (b) output.}%
\label{fig11}
\end{center}
\end{figure}




  \item Finally, we evaluate the performance of the sampling algorithms via their (practical) output error approximations. Again, the output error approximation \eqref{eq:RB_nine} by the standard algorithm with ($\tilde{N},\eta_T=0.8$) is used. (Namely, the pairs $(N,\tilde{N})$ in \eqref{eq:RB_nine} is taken following and after implementing the GO ($\tilde{N},\eta_T=0.8$) algorithm.) Specifically, Fig.\ref{fig12b} shows the graphs of $\tilde{\Delta}^{\rm max,rel}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} \left| \frac{\Delta_s(\mu)}{s(\mu)} \right|$ (defined in \eqref{eq:RB_eight}, \eqref{eq:RB_nine}) as functions of $N$ by 3 algorithms: standard ($\tilde{N},\eta_T=0.8$), GO ($\tilde{N},\eta_T=0.8$) and GO ($\tilde{N}=2N$), respectively. Fig.\ref{fig12a} is again a repetition of Fig.\ref{fig11b} -- which shows the max relative RB true errors to compare with the max output error approximations in Fig.\ref{fig12b}. The associated maximum/minimum effectivities of this error approximation (i.e., $\eta^{\max}_{s} = \max \limits_{\mu \in \Xi_{\rm test}} \eta_s(\mu)$, $\eta^{\min}_{s} = \min \limits_{\mu \in \Xi_{\rm test}} \eta_s(\mu)$, defined in \eqref{eq:RB_eight}, \eqref{eq:RB_nine}) are also shown in Fig.\ref{fig13a} and Fig.\ref{fig13b}, respectively. Fig.\ref{fig12a} and Fig.\ref{fig12b} again show that the GO algorithms give smaller RB true output errors than that provided by the standard algorithm (for the same $N$ basis functions). Fig.\ref{fig13a} and Fig.\ref{fig13b} also show much better effectivities of the standard and GO ($\tilde{N},\eta_T=0.8$) algorithms compared with that of the GO ($\tilde{N}=2N$) one. These results again confirm the good performance of the ``cross-validation'' process: the GO ($\tilde{N},\eta_T=0.8$) algorithm decreases the RB true error compared with the standard algorithm, and it also improves significantly the output effectivity estimation compared with the GO ($\tilde{N}=2N$) algorithm.

% It also seems that (without proof) $\Delta^{\rm go}_s(\mu)$ by the GO ($\tilde{N}=2N$) algorithm does not grow exponentially with time.



%  As observed, we see that the proposed error approximation $\Delta_s(\mu)$ is very close to the true output error (Fig.\ref{fig13a}) and that their effectivities are really good (Fig.\ref{fig13b}). In fact, the Fig.\ref{fig12b} and Fig.\ref{fig13b} confirm the good performance of the ``cross-validation'' process: the GO ($\tilde{N},\eta_T=0.8$) algorithm improves significantly the output effectivity estimation compared with the GO ($\tilde{N}=2N$) one. It also seems that (without proof) $\Delta_s(\mu)$ does not grow exponentially with time.

%      that its maximum effectivity is about $O(10)$ which is very good for a 3D time-dependent dynamic problem (Fig.\ref{fig13b}).

%%%%%%%%%%%%% fig 12 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_Impulse__TrueErrs.eps}} \label{fig12a} }
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_Impulse__ErrApprx.eps}} \label{fig12b} } \\
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_Impulse__MaxEff.eps}} \label{fig13a}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_Impulse__MinEff.eps}} \label{fig13b}}
\caption{ ``Impulse load'' case: comparison of (a) max relative RB true errors, (b) max relative error approximations, (c) corresponding max effectivities and (d) corresponding min effectivities of all sampling algorithms over $\Xi_{\rm test}$. }
\label{fig12}
\end{center}
\end{figure}


%%%%%%%%%%%%%% fig 13 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/fig13a.eps}} \label{fig13a} }%
%\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/fig13b.eps}} \label{fig13b} }%
%\caption{ (a) Comparison of max (relative) RB true output error versus max output error approximation, and (b) their corresponding max effectivity over $\Xi_{\rm test}$ using GO ($\tilde{N}, \eta_T=0.8$) algorithm. }%
%\label{fig13}
%\end{center}
%\end{figure}

\end{enumerate}







\paragraph{An arbitrary loading case}\label{paragraph_3dDental_ArbLoad}


All the results shown so far are for the unit impulse load, which is of very limited use in practice. Fortunately, as mentioned in Section \ref{subsec_ProbStat_ImpRes}, we can use the Duhamel's principle to convolute all the ``unit'' RB solutions with an ``arbitrary'' load to obtain correspondingly ``arbitrary'' RB solutions/outputs. Here, we also test and present some results of our error approximation $\Delta_s(\mu)$ with the ``arbitrary'' load defined in Fig.\ref{fig2b}. We note that the error approximation (T\ref{tab1}k) is good only for the impulse load since we built the goal-oriented RB basis functions for this load. For more complicated loadings, the error approximation might be much worse since there might be error in the convolution. Similar to paragraph d) of Section \ref{paragraph_3dDental_ImpLoad} above, Fig.\ref{fig14a} and Fig.\ref{fig14b} show the graphs of max relative RB true error and max relative error approximation; while Fig.\ref{fig15a} and Fig.\ref{fig15b} show the corresponding max/min effectivities of the 3 algorithms: standard, GO ($\tilde{N},\eta_T=0.8$) and GO ($\tilde{N}=2N$) over $\Xi_{\rm test}$, respectively. Again, we obtain quite similar phenomena as in Fig.\ref{fig12}: the GO algorithms beat the standard one in RB true error and error approximation; and the GO ($\tilde{N},\eta_T=0.8$) algorithm provides the best effectivities among all algorithms. The results again prove the preeminence of the proposed ``cross-validation'' process to choose adaptively and sufficiently $\tilde{N}$ for a particular $N$ of interest.


%Fig.\ref{fig14a} and Fig.\ref{fig14b} show the results using the GO ($\tilde{N}=2N$) algorithm; while Fig.\ref{fig15a} and Fig.\ref{fig15b} show that using the GO ($\tilde{N},\eta_T=0.8$) algorithm over $\Xi_{\rm test}$, respectively. Generally, the results are good for this ``arbitrary'' load case. Again, we observe from Fig.\ref{fig14b} and \ref{fig15b} that the proposed algorithm helps to improve significantly the effectivity measure of the problem and hence outperforms the GO ($\tilde{N}=2N$) algorithm. The ``cross-validation'' process did choose adaptively and sufficiently $\tilde{N}$ for a particular $N$ of interest (see Fig.\ref{fig10a}).


%Fig.\ref{fig14} shows the results using the GO ($\tilde{N}=2N$) algorithm, and Fig.\ref{fig15} show that using the GO ($\tilde{N},\eta_T=0.8$) algorithm, respectively. Generally, the results are quite good for this ``arbitrary'' load: the curves behave quite similarly to that of the impulse load. In particular, Fig.\ref{fig15} shows much better results than Fig.\ref{fig14} meaning that the proposed GO ($\tilde{N},\eta_T=0.8$) algorithm performs better than the GO ($\tilde{N}=2N$) algorithm (again, this is straightforward from Fig.\ref{fig10a}).


%%%%%%%%%%%%% fig 14 %%%%%%%%%%%%%
\begin{figure}[h!]
\begin{center}
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_1Gauss__TrueErrs.eps}} \label{fig14a}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_1Gauss__ErrApprx.eps}}  \label{fig14b}} \\
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_1Gauss__MaxEff.eps}}  \label{fig15a}}%
\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/3D_1Gauss__MinEff.eps}}  \label{fig15b}}%
\caption{ ``Arbitrary load'' case: comparison of (a) max relative RB true errors, (b) max relative error approximations, (c) corresponding max effectivities and (d) corresponding min effectivities of all sampling algorithms over $\Xi_{\rm test}$. (Shorter legends were used in Fig.\ref{fig15a} and Fig.\ref{fig15b} to avoid affecting the graphs.)}
\label{fig14}
\end{center}
\end{figure}


%%%%%%%%%%%%%% fig 15 %%%%%%%%%%%%%
%\begin{figure}[h!]
%\begin{center}
%\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/fig15a.eps}} \label{fig15a} }%
%\subfigure[]{\resizebox*{7cm}{!}{\includegraphics{figs/fig15b.eps}} \label{fig15b} }%
%\caption{``Arbitrary load'': (a) Comparison of max (relative) RB true output error versus max output error approximation, and (b) their corresponding max effectivity over $\Xi_{\rm test}$ using GO ($\tilde{N}, \eta_T=0.8$) algorithm. }%
%\label{fig15}
%\end{center}
%\end{figure}





Finally, regarding computational time, all computations were performed on a desktop Intel(R) Core(TM) i7-3930K CPU @3.20GHz 3.20GHz, RAM 32GB, 64-bit Operating System. The computational time for the RB solver ($t_{\rm RB(online)}$), the CPU-time for the FEM solver by our code ($t_{\rm FEM}$) and the CPU-time saving factor $\kappa=t_{\rm FEM}/t_{\rm RB(online)}$ are listed on Table \ref{tab3}, respectively. We also provide roughly the computational time for the error approximations ($\Delta_s(\mu)$) using both GO algorithms in that Table. We see that the RB solver is approximately $O(1000)$ times and the error approximation is $O(100)$ faster than the FEM solver; and thus it is clear that the RB is very efficient and reliable for solving time-dependent dynamic problems \cite{rabczuk2007three, song2008comparative, rabczuk2010three, nguyen2012extended, nguyen2010alternative}.

% \footnote{The work focuses on real-time context with many online computations, the offline stage is done once and expensive: the total computational time is approximately 2 weeks (including all FEM solutions/outputs and RB true errors of the standard and goal-oriented algorithms) on this computer.}



%%%%%%%%%%%%% table 3 %%%%%%%%%%%%%
\begin{table}[h!]
\center
\caption{Comparison of the CPU-time for a FEM and RB analysis.}
{\begin{tabular}{|c|c|c|c|c|c|}
\hline \rule{0pt}{2.5ex}
    N & $t_{\rm RB(online)}$ ({\rm sec}) & $t_{\rm FEM}$ ({\rm sec}) & $\kappa=t_{\rm FEM}/t_{\rm RB(online)}$  & $t_{\Delta_s(\mu)}$ ({\rm sec}) & $t_{\Delta_s(\mu)}$ ({\rm sec}) \\
      &     &   &   & ($\tilde{N}=2N$) &  ($\tilde{N},\eta_T=0.8$)  \\[0.5ex]
\hline \rule{0pt}{2.5ex} \hspace{-2.5mm}
   10 & 0.006757 & 29  & 4291   & 0.008329 & 0.038035  \\
   20 & 0.008391 & 29  & 3456   & 0.014242 & 0.099046  \\
   30 & 0.011123 & 29  & 2607   & 0.031723 & 0.049015  \\
   40 & 0.014531 & 29  & 1996   & 0.042440 & 0.130987  \\
   50 & 0.024381 & 29  & 1189   & 0.061122 & 0.192829  \\
   60 & 0.031196 & 29  & 930    & 0.075049 & 0.328012  \\
\hline
\end{tabular}}
%\tabnote{$^{\rm a}$This is not as accurate, owing to numerical
%error.} \tabnote{$^{\rm b}$An example table footnote to show the
%text turning over when a long footnote is inserted.}
\label{tab3}
\end{table}






\section{Conclusion}\label{sec_Conclusion}

A new goal-oriented POD--Greedy sampling algorithm was proposed. The proposed algorithm cooperates and improves further the standard POD--Greedy algorithm by using the output error approximation rather than the dual norm of residual as error indicator within the Greedy iterations. A cross-validation process is proposed to choose adaptively the dimension of the larger-dimensional auxiliary RB space in the output error approximation. The proposed strategy is verified by investigating a 2D linear plane strain problem and a 3D dental implant problem in the time domain. It is demonstrated that this type of error indicator will guide the Greedy iterations to select the parameter samples to optimize the true output error. In comparison with the standard algorithm, we conclude that our proposed algorithm performs better -- in terms of output's accuracy, and quite similar -- in terms of solution's accuracy. The proposed algorithm is applicable to various (regular) output functionals and is thus very suitable within the goal-oriented RB approximation context.



%The derived \textit{a posteriori} error representation is simple, effective, easy to implement and totally discrete although the problem is stated in space-time setting. The derivation technique is also applicable to any other finite difference schemes in time such as Crank-Nicolson, central finite difference or Backward Euler schemes, etc.





\section*{Acknowledgements}

We are sincerely grateful for the financial support of the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578) entitled ``Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery''. We also thank the anonymous reviewers for many contributive comments which improve the quality of this paper.




\bibliographystyle{ieeetr}        % plain, unsrt, ieeetr, abbrv, siam

\bibliography{myrefs_s3}




\appendix

\section{Non-dimensionalization of damped wave equation (weak form)}\label{sec_Appx}

In the following, the physical and non-dimensional terms will be denoted as $\bar{x}$ and $x$, respectively. With all the terms defined as in \eqref{eq:ProbStat_three}--\eqref{eq:ProbStat_four} and assuming zero mass-proportional Rayleigh damping coefficients, the weak form of the dimensional damped wave equation considered in this work has the following form: $\forall \bar{t} \in [0,\bar{T}]$ and a test function $\bar{v}$

\begin{equation}\label{eq:Appx_one}
\int_{\bar{\Omega}} \bar{\rho} \bar{v}_i \frac{\partial^2 \bar{u}_i}{\partial \bar{t}^2} + \bar{\beta} \frac{\partial}{\partial \bar{t}} \int_{\bar{\Omega}} \frac{\partial \bar{v}_i}{\partial \bar{x}_j} \bar{C}_{ijkl} \frac{\partial \bar{u}_k}{\partial \bar{x}_l} + \int_{\bar{\Omega}} \frac{\partial \bar{v}_i}{\partial \bar{x}_j} \bar{C}_{ijkl} \frac{\partial \bar{u}_k}{\partial \bar{x}_l} = \int_{\bar{\Gamma}_{\rm l}} \bar{v}_i \bar{\phi}_i \qquad {\rm on} \; \bar{\Omega},
\end{equation}


which subjects to the boundary conditions

%\begin{align}\label{eq:Appx_two}
%\bar{u}_i &= \bar{u}_i^{\rm d} \qquad {\rm on} \; \bar{\Gamma}_D, \\
%\bar{\sigma}_{ij} n_j &= \bar{C}_{ijkl} \left( \frac{\partial \bar{u}_k}{\partial \bar{x}_l} + \bar{\beta} \frac{\partial}{\partial \bar{t}} \frac{\partial \bar{u}_k}{\partial \bar{x}_l} \right) n_j = \bar{\phi}_i \qquad {\rm on} \; \bar{\Gamma}_{\rm l},
%\end{align}

\begin{subequations}\label{eq:Appx_two}
\begin{equation}
\bar{u}_i = \bar{u}_i^{\rm d} \qquad {\rm on} \; \bar{\Gamma}_D, \label{eq:Appx_two-a}
\end{equation}
\begin{equation}
\bar{\sigma}_{ij} n_j = \bar{C}_{ijkl} \left( \frac{\partial \bar{u}_k}{\partial \bar{x}_l} + \bar{\beta} \frac{\partial}{\partial \bar{t}} \frac{\partial \bar{u}_k}{\partial \bar{x}_l} \right) n_j = \bar{\phi}_i \qquad {\rm on} \; \bar{\Gamma}_{\rm l}, \label{eq:Appx_two-b}
\end{equation}
\end{subequations}


and initial conditions


%\begin{align}\label{eq:Appx_three}
%\bar{u}_i(\bar{t}=0) &= \bar{u}^0_i \qquad {\rm on} \; \bar{\Omega},  \\
%\frac{\partial \bar{u}_i(\bar{t}=0)}{\partial \bar{t}} &= \bar{v}^0_i \qquad {\rm on} \; \bar{\Omega}.
%\end{align}

\begin{subequations}\label{eq:Appx_three}
\begin{equation}
\bar{u}_i(\bar{t}=0) = \bar{u}^{{\rm d},0}_i \qquad {\rm on} \; \bar{\Omega}, \label{eq:Appx_three-a}
\end{equation}
\begin{equation}
\frac{\partial \bar{u}_i(\bar{t}=0)}{\partial \bar{t}} = \bar{v}^{{\rm d},0}_i \qquad {\rm on} \; \bar{\Omega}. \label{eq:Appx_three-b}
\end{equation}
\end{subequations}




Let $L^*$, $\rho^*$ and $E^*$ are correspondingly the characteristic length, mass density and stress of the model and define the non-dimensional variables as

%\begin{align}\label{eq:Appx_four}
%x_i &= \frac{\bar{x}_i}{L^*},  &  u_i &= \frac{\bar{u}_i}{L^*},  &  \phi_i &= \frac{\bar{\phi}_i}{E^*}, \\
%C_{ijkl} &= \frac{1}{E^*} \bar{C}_{ijkl},  &  t &= \sqrt{\frac{E^*}{\rho^* L^{*^2}}} \; \bar{t},  &  \beta &= \sqrt{\frac{E^*}{\rho^* L^{*^2}}} \; \bar{\beta}.
%\end{align}

\begin{subequations}\label{eq:Appx_four}
\begin{equation}
x_i = \frac{\bar{x}_i}{L^*},  \qquad  u_i = \frac{\bar{u}_i}{L^*},  \qquad  \phi_i = \frac{\bar{\phi}_i}{E^*}, \label{eq:Appx_four-a}
\end{equation}
\begin{equation}
C_{ijkl} = \frac{1}{E^*} \bar{C}_{ijkl},  \qquad  t = \sqrt{\frac{E^*}{\rho^* L^{*^2}}} \; \bar{t},  \qquad  \beta = \sqrt{\frac{E^*}{\rho^* L^{*^2}}} \; \bar{\beta}. \label{eq:Appx_four-b}
\end{equation}
\end{subequations}



Under this transformation, the non-dimensional weak form of \eqref{eq:Appx_one} is defined as: $\forall t \in [0,T]$, $T = \sqrt{\frac{E^*}{\rho^* L^{*^2}}} \; \bar{T}$ and a test function $v$

\begin{equation}\label{eq:Appx_five}
\int_{\Omega} v_i \frac{\partial^2 u_i}{\partial t^2} + \beta \frac{\partial}{\partial t} \int_{\Omega} \frac{\partial v_i}{\partial x_j} C_{ijkl} \frac{\partial u_k}{\partial x_l} + \int_{\Omega} \frac{\partial v_i}{\partial x_j} C_{ijkl} \frac{\partial u_k}{\partial x_l} = \int_{\Gamma_{\rm l}} v_i \phi_i \qquad {\rm on} \; \Omega,
\end{equation}


with boundary conditions

%\begin{align}\label{eq:Appx_six}
%u_i &= \frac{\bar{u}_i^{\rm d}}{L^*} \qquad {\rm on} \; \Gamma_D, \\
%\sigma_{ij} n_j &= C_{ijkl} \left( \frac{\partial u_k}{\partial x_l} + \beta \frac{\partial}{\partial t} \frac{\partial u_k}{\partial x_l} \right) n_j = \frac{\bar{\phi}_i}{E^*} \qquad {\rm on} \; \Gamma_{\rm l},
%\end{align}

\begin{subequations}\label{eq:Appx_six}
\begin{equation}
u_i = \frac{\bar{u}_i^{\rm d}}{L^*} \qquad {\rm on} \; \Gamma_D, \label{eq:Appx_six-a}
\end{equation}
\begin{equation}
\sigma_{ij} n_j = C_{ijkl} \left( \frac{\partial u_k}{\partial x_l} + \beta \frac{\partial}{\partial t} \frac{\partial u_k}{\partial x_l} \right) n_j = \frac{\bar{\phi}_i}{E^*} \qquad {\rm on} \; \Gamma_{\rm l}, \label{eq:Appx_six-b}
\end{equation}
\end{subequations}

and initial conditions

%\begin{align}\label{eq:Appx_seven}
%u_i(t=0) &= \frac{\bar{u}^0_i}{L^*} \qquad {\rm on} \; \Omega,  \\
%\frac{\partial u_i(t=0)}{\partial t} &= \sqrt{\frac{\rho^*}{E^*}} \; \bar{v}^0_i \qquad {\rm on} \; \Omega.
%\end{align}

\begin{subequations}\label{eq:Appx_seven}
\begin{equation}
u_i(t=0) = \frac{\bar{u}^{{\rm d},0}_i}{L^*} \qquad {\rm on} \; \Omega, \label{eq:Appx_seven-a}
\end{equation}
\begin{equation}
\frac{\partial u_i(t=0)}{\partial t} = \sqrt{\frac{\rho^*}{E^*}} \; \bar{v}^{{\rm d},0}_i \qquad {\rm on} \; \Omega. \label{eq:Appx_seven-b}
\end{equation}
\end{subequations}



\end{document}






