% this file is called up by thesis.tex
% content in this file will be fed into the main document

%: ----------------------- name of chapter  -------------------------
\chapter{DA for the unsteady NSE} \label{C:da_unsteady}% top level followed by section, subsection


%: ----------------------- paths to graphics ------------------------

% change according to folder and file names
\ifpdf
    \graphicspath{{5/figures/PNG/}{5/figures/PDF/}{5/figures/}}
\else
    \graphicspath{{5/figures/EPS/}{5/figures/}}
\fi

%: ----------------------- contents from here ------------------------

In previous chapters we focused our attention on the assimilation of data collected at the same time instant and we developed an assimilation procedure for steady phenomena. In the context of hemodynamics this is clearly limiting; during one heartbeat the flow dynamics might vary significantly. If we want to capture the whole behavior of blood circulation we need to account for snapshots of the blood velocity in different stages of the cardiac cycle. As an example, in the case study discussed in Chapter \ref{C:introduction} (presented in \cite{brummer} by Dr. Brummer and collaborators) using phase-velocity encoded MRIs of the aorta, 20 time points per cardiac cycle were collected on several internal slices of the vessel.

As anticipated in Chapter \ref{C:introduction} there are plenty of methods to treat time dependent data. In this chapter we consider variational approaches inspired by the steady formulation presented in Chapter \ref{C:da_nse}; in this context it is worth spending a few words on the issue of solving the adjoint problem (see Section \ref{ss:introduction_methods_control} on variational approaches). In fact, the treatment of adjoint equations is one of the major issues in using control-based methods for unsteady problems. These methods are formulated as the minimization of the misfit between data and velocity subject to the state equations. When we solve this problem finding the stationary points of the Lagrangian functional we have to deal with a coupled system involving the state and adjoint equations and the optimality condition. The most important observation is that the adjoint equations are posed backward in time, while the state ones are forward in time. This means that, when discretizing the system, it is not possible to ``march in time'' \cite{gunzburger_perspectives_2002}; thus, solutions to the optimality system at all time levels are coupled to each other. This fact makes the monolithic solution to the system, the so called {\it one-shot method} \cite{gunzburger_perspectives_2002}, in general prohibitively expensive, even for simple differential equations. 

In this chapter we first work on simple test cases where one-shot methods are still affordable and, then, we propose a variational method formulated so to circumvent the solution of the unsteady adjoint problem, we will refer to this method as {\it intermittent} \gls{DA} meaning that the optimization is intermittent in time and, hence, local. In literature several strategies are available; the simplest is to introduce iterative methods for uncoupling the state and adjoint equations. In this case, given a guess for the control variable, state and adjoint are sequentially solved and used in the optimality equation to update the control variable. This is to be done until some criterion for the convergence is satisfied; this process might be very slow or even diverge. Other methods include sensitivity-based and adjoint-based optimization methods, see \cite{gunzburger_perspectives_2002} for an extensive discussion.

This chapter is organized as follows. In Section \ref{S:unsteady_lin} we treat the unsteady linearized \gls{NSE}. For non-periodic problems we propose two  variational formulations, one-shot and intermittent; while in the periodic case we consider a one-shot method. In Section \ref{S:unsteady_lin_numres} we show numerical results stressing consistency and accuracy of the proposed schemes. 
In Section \ref{S:unsteady_nonlin} we present an intermittent algorithm for the non-periodic unsteady \gls{NSE} and we propose a generalization to the periodic case. Also, on the basis of such algorithms, we propose a statistical approach to \gls{DA} relying on schemes introduced in Chapter \ref{C:da_bayesian}. In Section \ref{S:unsteady_nonlin_numres} we compare the performance of deterministic and statistical estimators and we test consistency ad accuracy of the periodic approach.
Numerical results in this chapter are preliminary and represent a first attempt to the solution of this complex problem. The aim is to conduct a preliminary investigation which may guide the choice of effective methods.

%\begin{figure}[t]
%\begin{center}
%\begin{tabular}{cc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
%\includegraphics[width=0.45 \textwidth]{sequential.pdf} &  
%\includegraphics[width=0.45\textwidth]{non_sequential.pdf} \\
%\end{tabular}
%\end{center}
%\caption{Four possible schemes for the assimilation of data in time.}
%\label{F:unsteady_approaches}
%\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The linearized problem}\label{S:unsteady_lin}
We consider the unsteady Stokes flow in $\dom$, a domain in $\mathbb{R}^d$ ($d = 2,\,3$) featuring an inflow section $\Gamma_{in}$, an outflow section $\Gamma_{out}$ and a wall $\Gamma_{wall}$. 
Variables of interest are velocity $\uu(\xx,t)\in L^2((t_0,\,T];{\bf H}^1_{\Gamma_{wall}}(\dom))$ and pressure 
$p(\xx,t)\in L^2((t_0,\,T];L^2(\dom))$, which satisfy the unsteady Stokes equations. We assume to have $N_s\cdot N_t$ velocity measures $\dd_n\in\mathbb{R}^{N_s}$ available at some sites $\xx^m_{i,n}\in\dom$ for $i=1,\,...,\,N_s$ at time $t_n$, $n=1,\,...,\,N_t$. We formulate the unsteady assimilation problem as a control problem of the form
\begin{equation}\label{E:min_cont_unsteady}
\begin{array}{l}
\min\limits_{\hv} \mathcal{J}(\uu,\,\hv) = dist\big(\uu_n, \, \dd_n;\, \small n=1,...,N_t \normalsize \big) + \mathcal{R}(\hv_n;\, \small n=1,...,N_t \normalsize) \\
[3mm]
\hbox{s.t.} \;\; 
\left\{
\begin{array}{ll}
\uu_t -\nu\;\gr\cdot(\gr\uu + \gr \uu^{\rm T}) + \gr p=\mathbf{s} & \textrm{in}\;\; \dom, \;\;{\rm for} \;\;t\in (t_0,\,T],\\
[3mm]
\gr\cdot\uu = 0                             								         & \textrm{in}\;\; \dom, \;\;{\rm for} \;\;t\in (t_0,\,T],\\
[3mm]
\uu = \mathbf{0} 																	 & \hbox{on}\;\;\Gamma_{wall}, \;\;{\rm for} \;\;t\in (t_0,\,T],\\
[2mm]
\displaystyle-\nu \,(\gr\uu+\gr\uu^\TT)\mathbf{n}+ p\,\mathbf{n}=\mathbf{h} & \textrm{on}\;\;\Gamma_{in}, \;\;{\rm for} \;\;t\in (t_0,\,T],\\
[3mm]
\displaystyle - \nu \,(\gr\uu+\gr\uu^\TT) \mathbf{n}+p\,\mathbf{n}=\mathbf{g}& \textrm{on}\;\;\Gamma_{out}, \;\;{\rm for} \;\;t\in (t_0,\,T],\\
[3mm]
\uu(\xx,t_0) = \uu_0(\xx)																 &	{\rm in}\;\;\dom.
\end{array}\right. 
\end{array}
\end{equation}
Here $\uu_n = \uu(\xx,t_n)$ and $\hv_n = \hv(\xx,t_n)$; in this case, the distance function $dist(\cdot,\cdot)$ depends on velocity and data at measurement instants and will be specified later on. Other functions and parameters are defined as for the steady case. For the reasons specified in Chapter \ref{C:da_nse} we use the normal stress $\hv$ as the control variable. 

One of the major issues in dealing with time dependent phenomena is the prescription of the initial condition; in fact, in practical applications the velocity at time $t_0$ is not available. Thus, we proceed following two strategies. The first consists in using a reasonable and feasible guess on the initial data; the second in assuming periodicity in time, setting $\uu_0 = \uu(\xx,\tau)$, $\tau$ being the period. The latter stems from the physics of the problem; in fact, it is known that the cardiovascular system is characterized by periodic phenomena and in such cases $\tau$ is the the length of the cardiac cycle. Another possibility, often used in meteorology, is to include the initial condition in the set of control variables. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Non-periodic formulation}
As in the steady case, we consider a \gls{DO} technique (see Chapters \ref{C:preliminary} and \ref{C:da_nse}); the space discretization is performed by means of the \gls{FE} method and the time discretization by means of the backward Euler (\gls{BE}) scheme \cite{qss}. Let ${\rm M}_{u,p}\in \mathbb{R}^{N_u+N_p,N_u+N_p}$ be the matrix that extends the velocity mass matrix, ${\rm M}$, to the entire set of \gls{DOF}s,
\begin{equation}\label{E:Mup}
{\rm M}_{u,p}=\left[\begin{array}{c c}
{\rm M}	&	{\rm O} 	\\
{\rm O}			&	{\rm O}		\\
\end{array}\right].
\end{equation}
In this context we assume the velocity mass matrix ${\rm M}$ to account for the time step, i.e. 
\begin{equation}\label{E:matrixM_velocity}
\left[{\rm M}\right]_{ij} =  \dfrac{1}{\Delta t}\int\limits_{\dom}\displaystyle\boldsymbol{\varphi}_{i}\boldsymbol{\varphi}_{j} \, {\rm d \xx} \quad i,\,j = 1,\,...,\,N_{u};
\end{equation}
where, $\boldsymbol{\varphi}$ is the Lagrangian \gls{FE} basis function for the velocity.
The space-time discretization of the constraints in \eqref{E:min_cont_unsteady}, for all $n=1,\,...,\,N_f$ reads
\begin{equation}\label{E:unsteadys_discr_compact}
{\rm S}\VV_n + {\rm M}_{u,p}(\VV_n-\VV_{n-1})= \textrm{R}^{\TT}_{in}{\rm M}_{in}\HH_{n}+ \FF_n.
\end{equation}
where, $N_f = \frac{T}{\Delta t}$, $\Delta t$ being the time discretization step. $\VV_n$ and $\HH_n$ are the discretization of velocity and pressure and of normal stress at time $t_n$. Matrices $\rm S$, $\RR_{in}$ and ${\rm M}_{in}$ have been introduced in Section \ref{S:da_nse_DO}. 
For the sake of simplicity we assume that the set of measurement instants $I_{meas}$ is contained in the set of time discretization instants $I_{discr}$, i.e. $I_{meas} = \{t_j\}_1^{N_t} \subseteq I_{discr} = \{t_n\}_1^{N_f}$; this assumption, in general, is also realistic. In the solution of the non-periodic assimilation problem we consider two possible cases: $I_{meas} = I_{discr}$ and $I_{meas} \subset I_{discr}$.
%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Formulation I} When discretization and measurement instants coincide we consider a one-shot method solved with a \gls{DO} approach; here, the constraints are the state equations discretized in space and time. We formulate the problem as follows,  
\begin{equation}\label{E:opt_unsteady_formI}
\begin{array}{l}
\min\limits_{\HH_{n=1,...,N_f}}\;\, \sum\limits_{n=1}^{N_f}\dfrac{1}{2} \|{\rm D}_n\VV_n - \dd_n\|_2^2 + \dfrac{\alpha}{2} \|{\rm L}\HH_n\|_2^2\\
\textrm{s.t.} \quad ({\rm S} + {\rm M}_{u,p})\VV_n = {\rm R}^\TT{\rm M}_{in}\HH_n + {\rm M}_{u,p}\VV_{n-1}+\FF_n,
\quad \forall\,n = 1,\,...,\,N_f.
\end{array}
\end{equation}
Here, the subscript $n$ corresponds to the time instant $t_n$ in $I_{discr}$. We rewrite \eqref{E:opt_unsteady_formI} in a more compact form
\begin{equation}\label{E:opt_unsteady_oneshot}
\begin{array}{l}
\min\limits_{\HH} \dfrac{1}{2} \|\overline{\rm D}\overline{\VV} - \overline{\dd}\|_2^2 + \dfrac{\alpha}{2} \|\overline{\rm L}\overline{\HH}\|_2^2\\
\textrm{s.t.} \quad \overline{\rm S}\overline{\VV} = \overline{\rm R}\overline{\HH} +\overline{\FF},
\end{array}
\end{equation}
where
\begin{equation}\label{E:Slarge_unsteady}
\overline{\rm S}=\left[\begin{array}{c c c c c}
{\rm S}+{\rm M}_{u,p} 			& {\rm O} 	& ... 			&    				&  	\\
-{\rm M}_{u,p} 	& {\rm S}+{\rm M}_{u,p} 	& {\rm O} 	& ... 				& 	\\
 ... 				& 				& 				& 					& 	\\
 ... 				& 				& 				& 					& 	\\ 
 					&				& ... 			& -{\rm M}_{u,p} & 	{\rm S}+{\rm M}_{u,p}
\end{array}\right];
\end{equation}
\begin{equation}\label{E:Rlarge_unsteady}
\overline{\rm R}=\left[\begin{array}{c c c c c}
{\rm R}_{in}^{\TT} {\rm M}_{in} & {\rm O} 						& ... 		 	&    	&  	\\
{\rm O}						& {\rm R}_{in}^{\TT} {\rm M}_{in} 	& {\rm O} 	& ... 	& 	\\
...								&  									&   			& 		& 	\\
...								&					  				&   			& 		& 	\\ 
 								&				 					& 				& ... 	& {\rm R}_{in}^{\TT} {\rm M}_{in}
\end{array}\right];
\end{equation}
\begin{equation}\label{E:def_vect_unsteady}
\begin{array}{l}
\overline{\rm D}=diag\left( {\rm D}_1,...,{\rm D}_{N_t} \right);\\

\overline{\rm L}=diag\left( {\rm L},...,{\rm L} \right);\\

\overline{\dd}=\left[\dd_1\;\;\dd_2\;\; ...\;\; \dd_{N_t} \right]^{\TT};\\

\overline{\VV}=\left[\VV_1 \;\; \VV_2 \;\; ... \;\; \VV_{N_t} \right]^{\TT};\\

\overline{\HH}=\left[\HH_1 \;\; \HH_2 \;\; ... \;\; \HH_{N_t} \right]^{\TT};\\

\overline{\FF}=\left[ \FF_1 + {\rm M}_{u,p}\VV_0 \;\;\;\; \FF_2 \;\; ... \;\; \FF_{N_t} \right]^{\TT}.
\end{array}
\end{equation}
Here, ${\rm D}_n$ is a selection matrix such that ${\rm D}_n\VV_n$ corresponds to the evaluation of the velocity vector at time $t_n$ on the measurement sites $\xx_{i,n}^m$, for $i= 1,\,...,\,N_s$; $\rm L$ is the discretization of the gradient.

Problem \eqref{E:opt_unsteady_oneshot} is an algebraic optimization problem with linear constraints; it has the same form of problem \eqref{E:min_discr} introduced in Chapter \ref{C:da_nse} for the steady case. For this reason, we solve it using the same approach, i.e. the reduced Hessian method. Matrix dimensions in this extended formulation might be large; nevertheless, $\overline{\rm S}$ is lower block triangular, with same diagonal and off-diagonal blocks. Thus, memory requirements and computational effort are affordable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Formulation II} When measurement instants are a subset of time discretization steps we refer to an intermittent assimilation technique consisting in two stages: assimilating and advancing in time. This formulation has affordable memory requirements and it is easy to be implemented. The following algorithm summarizes the two steps. For $n=1,\,...,\,N_f$, given a guess for velocity and pressure at time $t=t_0$, say $\VV_0$
\begin{description}
\item {\bf if} $(t_n\in I_{meas})$ solve (with the reduced Hessian method): 
\begin{equation}\label{E:assimilating_stepII}
\begin{array}{l}
\min\limits_{\HH_n} \dfrac{1}{2} \|{\rm D}_n{\VV}_n - \dd_n\|_2^2 + \dfrac{\alpha}{2} \|{\rm L}{\HH}_n\|_2^2\\
\textrm{s.t.} \quad{\rm S}\VV_n + {\rm M}_{u,p}(\VV_n - \VV_{n-1})= {\rm R}_{in}^{\TT}{\rm M}_{in}\HH_n + \FF_n;
\end{array}
\end{equation}
\item {\bf else} solve the unsteady Stokes problem \eqref{E:unsteadys_discr_compact} using suitable boundary conditions.
\end{description} 
The choice of boundary conditions for the advancing stage depends on how much we trust available measures. We consider the following options at time $t_n\in I_{discr}$, such that $t_n\in[t_{j-1}, \;t_{j}]$, where $t_j$ and $t_{j-1} \in I_{meas}$, are two subsequent measurements instants. 
\begin{enumerate}
\item {\it Dirichlet} with data $\Phi(\dd_j, \, \dd_{j-1}) = \frac{t_n- t_{j-1}}{t_j-t_{j-1}}\dd_{j-1} + \frac{t_j - t_n}{t_j-t_{j-1}}\dd_{j} $ on $\Gamma_{in}$.
\item {\it Neumann} with data $\HH_{j-1}$, the value of the control variable resulting from the previous assimilation step. 
\end{enumerate}

The choice of 1. corresponds to assuming that data are not significantly affected by noise. On the other hand the choice of 2. corresponds to assuming that the normal stress does not vary considerably between two subsequent measurement instants. This assumption might be reasonable especially in applications where data are collected frequently.
\begin{remark}
The optimality result in Proposition \ref{P:suff_cond_Qmatrix}, Chapter \ref{C:da_nse}, is promptly extended to this unsteady formulation at each optimization step, i.e. equation \eqref{E:assimilating_stepII}. It is easy to show that Lemma \ref{L:Svw} is applicable to this optimization problem. In fact, as stated in Remark \ref{R:lemma}, such lemma is not specific for the steady Stokes problem but for any \gls{PDE} whose discretization matrix corresponding to a homogeneous Dirichlet problem on $\Gamma_{in}$ is non-singular; which is the case for equation \eqref{E:assimilating_stepII}.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Periodic formulation}
If a guess on the initial velocity is not available, we considere a periodicity constraint, in view of our application. For the sake of simplicity we assume $I_{meas} = I_{discr}$; thus, the formulation is a modified version of problem \eqref{E:opt_unsteady_oneshot} with constraint
\begin{equation}\label{E:unsteady_state_periodic}
\overline{\rm S}_p\overline{\VV} = \overline{\rm R}\overline{\HH} +\overline{\FF},
\end{equation}
where
\begin{equation}\label{E:S_periodic}
\overline{\rm S}_p=\left[\begin{array}{c c c c c}
{\rm S}+{\rm M}_{u,p} 			& \textrm{O} 	& ... 			&    				&  	-\textrm{M}_{u,p}\\
-\textrm{M}_{u,p} 	& \textrm{S}+{\rm M}_{u,p} 	& \textrm{O} 	& ... 				& 	\\
 ... 				& 				& 				& 					& 	\\
 ... 				& 				& 				& 					& 	\\ 
 					&				& ... 			& -\textrm{M}_{u,p} & 	\textrm{S}+{\rm M}_{u,p}
\end{array}\right].
\end{equation}
As opposed to $\overline{\rm S}$, this matrix is not block triangular; the system is coupled. This fact makes the numerical solution of the problem more challenging. In the following section we provide details for the choice of suitable preconditioners, which are mandatory for the numerical solution of the system, and we discuss difficulties caused by the complex structure of the problem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{remark}\label{R:V0_VFpenalization} Another possibility for forcing periodicity is to consider formulation \eqref{E:opt_unsteady_oneshot} adding to the functional the penalizing term $\beta\|\UU_0 - \UU_{\tau}\|$, $\tau$ being the length of the period and $\beta$ a penalizing parameter chosen a priori.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical results}\label{S:unsteady_lin_numres}
We present numerical results obtained for non-periodic and periodic problems. For Formulation I, we consider different choices of $\Delta t$; for Formulation II we compare different choices of boundary conditions on $\Gamma_{in}$. For the periodic case we discuss the choice of suitable preconditioners and we present comparison results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Non-periodic case} 
We consider 2-dimensional rectangular geometries as in Figure \ref{F:rect_triang}, computational grids are generated with {\tt FreeFem++}, see Section \ref{S:preliminary_results}. We choose $\mathbf{s}=\mathbf{0}$, $\dom=[0,\;5]\times[-0.5,\;0.5]$, $\mathbf{g}=\mathbf{0}$, $\uu_0(\xx)=\mathbf{0}$ and we assume that data are available on three internal layers $\{x=1,\;x=2.5,\;x=4\}$ in correspondence of grid nodes.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.9 \textwidth]{meshFE.png} 
\end{center}
\caption{Rectangular triangulation on $\dom$.}
\label{F:rect_triang}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Space discretization is performed by means of the \gls{FE} method using compatible \gls{FE} spaces P1bubble-P1 for velocity and pressure. For the time discretization we rely on the \gls{BE} scheme. 

For data generation and error analysis we consider a reference solution, say $\UU_{\rm FE}$, as done in Section \ref{S:da_nse_real_geometries}. This is a \gls{FE} solution computed on a very fine grid (240$\times$48 structured triangulation generated with {\tt FreeFem++}) using $\Delta t = 0.02$ and \gls{FE} spaces P1bubble-P1 for velocity and pressure. On $\Gamma_{wall}$ we prescribe a homogeneous Dirichlet condition, on $\Gamma_{out}$ a homogeneous Neumann condition and on $\Gamma_{in}$ a time dependent Dirichlet condition given by $\uu|_{\Gamma_{in}}(t) = [5t(1-4y^2)\;,\; 0\;]^{\TT}|_{\Gamma_{in}}$. Data are generated adding uniform random noise to $\UU_{\rm FE}$ with \gls{SNR} = 100. The choice of such small amount of noise is due to the fact that, at this stage, our investigation is mainly focused on the consistency of the methods. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Error analysis} The index of accuracy used in these test cases is the velocity relative error at time $T=t_{N_f}$ with respect to the reference solution $\UU_{\rm FE}$; this is defined as $E_U(N_f) = \frac{\|\UU_{N_f}-\UU_{{\rm FE},N_f}\|_2}{\|\UU_{{\rm FE},N_f}\|_2}$. The assimilation problem is solved in the time interval $[t_0,\,T] = [0, \, 1]$.

Concerning the choice of the regularization parameter $\alpha$ we utilize the \gls{DP} method introduced in Chapter \ref{C:da_nse} (see Section \ref{ss:da_nse_lin_results}) and the generalized cross validation (\gls{GCV}) method. The main idea behind \gls{GCV} is that if one data point is removed from the set of observations, then a good regularized solution should predict that missing point considerably well \cite{Hansen}. 

More precisely, we let ${\rm G}(\alpha)$ be the so called {\it influence} matrix defined, in our case, as ${\rm G}(\alpha) = \overline{\rm Z}(\overline{\rm Z}^\TT \overline{\rm Z}+N_s\alpha\overline{\rm L}^\TT\overline{\rm L})^{-1}\overline{\rm Z}^\TT$, where $\overline{\rm Z}= \overline{\rm D}\overline{\rm S}^{-1}\overline{\rm R}$. The \gls{GCV} choice of $\alpha$ is the minimizer of the \gls{GCV} function 
\begin{equation}\label{E:GCVfunction}
\mathcal{G}(\alpha) = \dfrac{N_s^{-1}\|({\rm I}-{\rm G}(\alpha))\overline{\dd}\|_2^2}{N_s^{-1}tr({\rm I}-{\rm G}(\alpha))}.
\end{equation}
Here, $\rm I$ is the identity matrix and $tr(\,\cdot\,)$ is the trace operator. The function $\mathcal{G}$ can be expressed as the weighted sum of squared misfits between data and predicted velocity $\overline{\rm D}_k\overline{\VV}_k \overline{\dd}_k$, where $\overline{\VV}_k$ is the regularized solution obtained leaving out the $k$-th data point (site), for $k=1,\,...,\,N_s$. Note that while \gls{DP} requires the norm of the noise to be available, \gls{GCV} does not. 

However, values of the regularization parameter are quite low (see Table \ref{T:formulation1}) and the smallest eigenvalue of the reduced Hessian is of the order of $10^{-2}$;  even though we did not prove an optimality result for Formulation I, this fact might be due to data location, which satisfies conditions of Proposition \ref{P:suff_cond_Qmatrix}.

Accuracy results for Formulation I are shown in Table \ref{T:formulation1}. Here, we report relative error $E_U(N_f)$ and optimal values of $\alpha$ ($\alpha_{DP}$ and $\alpha_{GCV}$ for \gls{DP} and \gls{GCV} respectively) for different choices of time steps. We consider $\Delta t$ = 0.1 and 0.2, hence, $N_t$ = $N_f$ = 5 and 10 respectively; recall that in this case $I_{meas} = I_{discr}$. In correspondence of \gls{DP} solutions are slightly more accurate, which is consistent with the fact that exact information about the noise is used. As we decrease the time step (increase the number of data in the time interval) we obtain a better accuracy; the improvement is due both to the smaller time-step and to the higher number of data.

For Formulation II we consider $\Delta t$ = 0.1 and $N_t$ = 5; this corresponds to performing one advancing step between two consecutive measurement instants. Comparison results between options 1. and 2. for the boundary condition on $\Gamma_{in}$ are reported in Table \ref{T:formulation2}. Since solutions have the same accuracy, this test case does not provide enough evidence to discard one of the two approaches. Furthermore, the accuracy is of the same order of results obtained with Formulation I with $\Delta t$ = 0.1; this fact suggests that accuracy might be mainly related to time discretization. 

Clearly, this analysis is still preliminary, but it shows that these methods are consistent in case of a low amount of noise. For drawing reasonable conclusion and choosing the most efficient technique it is fundamental to conduct a deeper investigation of the sensitivity to the noise and to the time discretization step.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | l l | c c |}
\hline 
$\Delta t$ & $\alpha_{DP}$   &  $\alpha_{GCV}$	&   $E_{U,DP}(N_f)$	&   $E_{U,GCV}(N_f)$ \\
\hline
0.2        & 3.5918e-7  	&  5.6266e-8		&	8.2962e-2       &   8.3021e-2        \\
0.1        & 1.1391e-6  	&  1e-9         	&   5.4405e-2       &   6.3660e-2        \\
\hline
\end{tabular}
\caption{Relative errors and regularization parameters for Formulation I.}
\label{T:formulation1}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| l | c |}
\hline 
Advancing 		& $E_{U}(T) $\\
\hline
{\it Dirichlet}	& 5.7602e-2\\
{\it Neumann}   & 5.6641e-2 \\
\hline
\end{tabular}
\caption{Relative errors for Formulation II in correspondence of \textit{Dirichlet} and \textit{Neumann} boundary conditions.}
\label{T:formulation2}
\end{center}
\end{table} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Time periodic solution}
In this test case we use same computational grid and software as for the non-periodic case. The reference solution is given by the following analytic solution on $\dom$: 
\begin{equation}\label{E:exact_timep}
\uu_{anl}(\xx,t)=\eta(1+\sin(100\pi t))[5t(1-4y^2)\;,\; 0\;]^{\TT}, \;\;\xx\in\dom,\;t>0.
\end{equation}
Here $\eta$ = 0.01, the period is $\tau$ = 0.02 and the source term is defined as $\mathbf{s}(\xx,\,t)=100\pi \eta \cos(100\pi t)[5t(1-4y^2)\;,\; 0\;]^{\TT}$.

We want to stress important numerical issues concerning the solution of the linear system induced by necessary conditions for optimality; we recall that matrix $\overline{\rm S}_p$ has a complex structure; hence, the choice of an effective preconditioner is a non-trivial task. We consider two options: the first consists in computing the incomplete Cholesky (\gls{IC}) factorization, it requires the explicit computation of the reduced Hessian and it is used only for comparison purposes. The other option does not require explicitly the matrix of the system and it is based on approximations of $\overline{\rm S}_p$ featuring a simpler structure. We propose two preconditioners: the first consists in using the non-periodic problem as a preconditioner for the periodic one. The second consists in building a block triangular matrix as follows. Let us rewrite matrix $\overline{\rm S}_p$ as 
\begin{equation}
\overline{\rm S}_p=\left[\begin{array}{cc}
\overline{\rm S}_{1,1} & \overline{\rm S}_{1,2}\\
\overline{\rm S}_{2,1} & \overline{\rm S}_{2,2}
\end{array}\right];
\end{equation}
where the first block row corresponds to time steps from 1 to $N_f-1$ and  the second block row accounts for velocity and pressure at the last time step. Then, the solution of system in $\overline{\rm S}_p$ is equivalent to solving a system (with right hand side modified accordingly) in $\widetilde{\rm S}_p$ defined as
\begin{equation}\label{E:blocktri}
\widetilde{\rm S}_p=\left[\begin{array}{cc}
\overline{\rm S}_{1,1} & \overline{\rm S}_{1,2}\\
\textrm{O} & \overline{\rm S}_{2,2}- \overline{\rm S}_{2,1}(\overline{\rm S}_{1,1})^{-1}\overline{\rm S}_{1,2}
\end{array}\right].
\end{equation}
Despite the block triangular structure, using \eqref{E:blocktri} is still not affordable; in order to reduce the computational effort we approximate $\overline{\rm S}_{1,1}$, in that block, with $\widetilde{\rm S}_{1,1} = diag({\rm M}_{u,p},\, ...,\, {\rm M}_{u,p})$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Error analysis and computational cost}
In this test case we use the following index of accuracy, the relative error in space and time with respect to the analytic solution
\begin{equation}
E_U = \displaystyle\sqrt{\frac{\sum_{n=1}^{N_f}\|\UU_n-\UU_{anl,n}\|^2}{\sum_{n=1}^{N_f}\|\UU_{anl,n}\|^2}};
\end{equation}
where $\UU_n$ and $\UU_{anl,n}$ are the computed velocity and discretized analytic velocity at time $t_n$. We choose SNR = 100, $N_t$ = $N_f$ = 5 and $\Delta t$ = 0.004, so to make the one-shot method affordable, as well as the Cholesky factorization. In Table \ref{T:period_lin1} we report the relative errors in correspondence of the solution of a forward problem with exact initial condition obtained with a direct solver ($E_{U,ref}$, used as a reference value of accuracy) and the one obtained with the the preconditioned conjugate gradient (\gls{PCG}) method ($E_{U,PCG}$). In solving with \gls{PCG} we use three different preconditioners: the identity, the incomplete Cholesky factor and the preconditioner induced by the approximation of $\overline{\rm S}_p$ with $\overline{\rm S}$. $E_{U,PCG}$ refers to all of them since they deliver the same solution. In Table \ref{T:period_lin2} we report CPU times for matrix assembling and solution of the linear system in correspondence of \gls{PCG} with the identity, \gls{IC} and $\overline{\rm S}$. The comparison between the second and the last column shows that the preconditioner $\overline{\rm S}$ does not bring a significant improvement. 

For testing the preconditioner $\widetilde{\rm S}_p$ we consider the \gls{GMRES} method. In Table \ref{T:period_lin3} we report relative errors, number of iterations and CPU times using the non-preconditioned and preconditioned \gls{GMRES}; in correspondence of different grids the methods deliver solutions with the same level of accuracy, which is mainly affected by time discretization. On the other hand, using the preconditioned \gls{GMRES} method the number of iterations and the computational time is halved compared to the non-preconditioned case. However, performance is still not satisfactory and more investigation is required. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | c c |}
\hline 
mesh 			& $E_{U,ref}$ & $E_{U,PCG}$\\
\hline
(40$\times$ 10) & 4.4e-2	& 	8.03e-2 \\
\hline
(60$\times$ 15) & 4.2e-2 		&  7.88e-2 \\
\hline
(80$\times$ 20) & 4.3e-2 		&  7.74e-2 \\
\hline
\end{tabular}
\caption{Relative errors for the reference solution and the assimilated one.}
\label{T:period_lin1}
\end{center}
\end{table}   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| l | l | l | l |}
\hline 
stage		& 	CG 		&   \gls{PCG}, \gls{IC} & \gls{PCG},  $\overline{\rm S}$ 		\\
\hline
assembling 	&	0 		& 	99 					& 0 								\\
\hline
solving 	& 	3e+3  	& 	55 					& 2e+3	 							\\
\hline
\end{tabular}
\caption{Computational times in seconds, in correspondence of a 80$\times$20 triangulation.}
\label{T:period_lin2}
\end{center}
\end{table}  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | c c c c c c |}
\hline 
mesh		   & $E_{U,G}$ 	& it$_{G}$	& CPU$_{G}$ & $E_{U,PG}$ &  it$_{PG}$ & CPU$_{PG}$\\
\hline
(40$\times$10) & 8.22e-2	& 51 		& 116		& 8.22e-2 	 & 32 		  & 86 \\
\hline
(60$\times$15) & 7.95e-2 	& 75     	& 657 		& 7.951e-2 	 & 38 		  & 394 \\
\hline
(80$\times$20) & 7.76e-2  	& 92     	& 1984 		& 7.761e-2 	 & 39 		  & 867 \\
\hline
\end{tabular}
\caption{Relative errors, number of iterations and CPU times for the non-preconditioned GMRES ($G$) and the preconditioned one ($PG$) with preconditioner $\widetilde{\rm S}_p$.}
\label{T:period_lin3}
\end{center}
\end{table}   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Formulation for the nonlinear problem}\label{S:unsteady_nonlin}
We consider the same problem setting as in the linearized case and we add the nonlinear term to the momentum equation \eqref{E:min_cont_unsteady}, which now reads
\begin{equation}\label{E:nonlinear_momentum}
\uu_t -\nu\;\gr\cdot(\gr\uu + \gr \uu^{\rm T}) + (\uu\cdot\gr)\uu + \gr p=\mathbf{s} \quad {\rm in}\; \dom,\;\;{\rm for}\;\;t\in [t_0,\,T].
\end{equation}
We design a scheme for treating both periodic and non-periodic phenomena; for the latter we assume a feasible guess on $\uu_0$ to be available. Our formulation is similar to the one introduced for the non-periodic linearized case and can be classified as intermittent \gls{DA}. We consider a \gls{DO} approach where space discretization is performed by means of the \gls{FE} method and time discretization by means of a semi-implicit method \cite{quarteroni_numerical_1994}, which reads as follows: for $n=1,\,...,\, N_f$
\begin{equation}\label{E:unsteady_semi_implicit}
{\rm S}_{n-1}\VV_n +{\rm M}_{u,p}(\VV_n-\VV_{n-1}) = {\rm R}_{in}^{\TT}{\rm M}_{in} \HH_n + \FF_n.
\end{equation}
Here, 
\begin{equation}
{\rm S}_n = \left[
\begin{array}{cc}
{\rm C}+{\rm A}_n & {\rm B}^{\TT}\\
{\rm B}			& {\rm O}
\end{array}\right];
\end{equation}
where ${\rm A}_n$ is the discrete advection operator associated with advection field $\VV_n$, defined as in equation \eqref{E:matrices_FE}. Other matrices are defined as in the linearized case.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Non-periodic case} In this case, even though not realistic, we assume that an initial velocity, $\VV_0$, is available. The following formulation is based on the assumption that data are available at each time instant of the time discretization, i.e. $I_{meas} = I_{discr}$. Let $\VV_0$ be the vector of velocity and pressure at time $t=t_0$, then for $n = 1,\, ...,\, N_f$ we solve
\begin{equation}\label{E:unsteady_algorithm}
\begin{array}{l}
\min\limits_{\HH_n} \; \dfrac{1}{2}\|{\rm D}_n\VV_n - \dd_n \|_2^2 + \dfrac{\alpha}{2} \|{\rm L}\HH_n\|_2^2\\

{\rm s.t.} \quad {\rm S}_{n-1}\VV_n +{\rm M}_{u,p}(\VV_n-\VV_{n-1}) = {\rm R}_{in}^{\TT}{\rm M}_{in} \HH_n + \FF_n.
\end{array}
\end{equation}
In our test cases we assume to have an analytic expression for velocity and pressure, say $\VV_{anl}$, and we set $\VV_0 = \VV_{anl}$. This algorithm is formally similar to \eqref{E:algorithm}, in Section \ref{S:da_nse_nonlinear}, which is designed for the steady case; hence, same numerical techniques are used for its numerical solution. Note that, with this approach we avoid the costly solution of the adjoint problem; nevertheless, the optimization in time is only local.
 
%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Periodic case} We assume that the solution is periodic with period $\tau$ and that a guess for $\VV_0$ is not available. We propose an algorithm based on the same iterative procedure as for the non-periodic case: set $\VV_0 = \VV_{N_f} = {\bf 0}$ and $N_f = \frac{\tau}{\Delta t}$, then
\begin{equation}\label{E:periodic_alg}
\begin{array}{l}
{\bf do}\\
\qquad \VV_0 = \VV_{N_f}; \; n=1;\\
\qquad {\bf for} \;\; n=1,\,...,\, N_f, \;\;{\rm solve} \\
\qquad \qquad \min\limits_{\HH_n} \; \dfrac{1}{2}\|{\rm D}_n\VV_n - \dd_n \|_2^2 + \dfrac{\alpha}{2} \|{\rm L}\HH_n\|_2^2\\
\qquad \qquad {\rm s.t.} \quad {\rm S}_{n-1}\VV_n +{\rm M}_{u,p}(\VV_n-\VV_{n-1}) = {\rm R}_{in}^{\TT}{\rm M}_{in} \HH_n + \FF_n\\
{\bf while}\; (\|\VV_0-\VV_{N_f}\|_2<\delta)\\
\end{array}
\end{equation}
This algorithm corresponds to solving repeatedly the non-periodic problem in $[0,\,\tau]$, each time improving the guess on $\VV_0$ using $\VV_{N_f}$, resulting from the previous simulation, until a convergence criterion is satisfied. As already stressed in Remark \ref{R:V0_VFpenalization}, another possible approach is to add the periodicity condition $\VV_0 = \VV_{N_f}$ to the set of constraints.
%%%%%%%%%%%%%%
\begin{remark}\label{R:extension_unsteady}
It is straightforward to generalize these algorithms to the case of $I_{meas} \subset I_{discr}$. One option is to perform interpolation of data in time in order to have approximated data at each step of the time discretization. This approach sounds promising on the basis of results obtained with interpolation in space, see Chapter \ref{C:da_nse}, Section \ref{S:wellp_2}. A second option is to proceed as in Formulation II: adding advancing steps of the form \eqref{E:unsteady_semi_implicit} between two measurement instants, i.e. in $[t_{j-1},\,t_j]$, $t_{j-1}$ and $t_j$ being two consecutive measurement instants. The formulation would read as follows. For $n=1,\,...,\,N_f$, given a guess for velocity and pressure at time $t=t_{n-1}$, say $\VV_{n-1}$
\begin{description}
\item {\bf if} $(t_n\in I_{meas})$ solve (with the reduced Hessian method): 
\begin{equation}\label{E:assimilating_step}
\begin{array}{l}
\min\limits_{\HH_n} \dfrac{1}{2} \|{\rm D}_n{\VV}_n - \dd_n\|_2^2 + \dfrac{\alpha}{2} \|{\rm L}{\HH}_n\|_2^2\\
\textrm{s.t.} \quad{\rm S}_n\VV_n + {\rm M}_{u,p}(\VV_n - \VV_{n-1})= {\rm R}_{in}^{\TT}{\rm M}_{in}\HH_n + \FF_n;
\end{array}
\end{equation}
\item {\bf else} solve the unsteady NSE using suitable boundary conditions.
\end{description} 
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Bayesian approach} The formulation presented in this sections can be generalized to a Bayesian approach to unsteady \gls{DA}. In fact, it is enough to compute either the \gls{MAP} or \gls{ML} estimator (see Chapter \ref{C:da_bayesian}, Section \ref{S:point_estimators}), instead of the deterministic one, at each time step in \eqref{E:unsteady_algorithm}. This corresponds to solving the following algorithm (for simplicity we report the non-periodic case only). Given a guess for the initial velocity $\VV_0$, for $n=1,\,...,\,N_f$ compute the \gls{MAP} estimator (as in equation \eqref{E:h_map_est}) solving 
\begin{equation}\label{E:unsteady_bayesian_algorithm}
\begin{array}{l}
(\Sigma_{pr}^{-1} + {\rm Z}_n^{\TT}\Sigma_{noise}^{-1}{\rm Z}_n)\HH_{n,MAP} = {\rm Z}_n^{\TT}\Sigma_{noise}^{-1}(\dd_n-{\bf b}_n).
\end{array}
\end{equation}
Here, ${\bf b}_n = ({\rm S}_{n-1}+{\rm M}_{u,p})^{-1}({\rm M}_{u,p}\VV_{n-1}+ \FF_n)$, ${\rm Z}_n = D_n({\rm S}_{n-1}+{\rm M}_{u,p})^{-1} {\rm R}_{in}^{\TT}{\rm M}_{in}$, $\Sigma_{prior}$ and $\Sigma_{noise}$ are the correlation matrices for the random variables for normal stress and measurement error. In the next section, deterministic and statistical estimators are compared.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical results}\label{S:unsteady_nonlin_numres}
In this section we consider the unsteady Navier-Stokes flow in a 3-dimensional cylindrical domain and we solve the problem using an axisymmetric formulation. The computational domain is a structured rectangular triangulation on the domain $[0,\,1.5]\times[0,\,6]$ corresponding to a section of a cylinder with radius 1.5 and height 6. Details on implementation and solvers can be found in Section \ref{ss:da_nse_lin_results}.

For data generation and accuracy analysis we consider the Womersley flow \cite{womersley} with period $\tau$ = 1 and viscosity $\nu=1$. We consider data on $\Gamma_{in}$ not satisfying sufficient conditions for optimality and we perform interpolation; 55 data are equally distributed along 5 internal slices not corresponding to grid nodes. Measurement noise is assumed to be normally distributed with correlation matrix defined as in \eqref{E:sigma_noise}, see Section \ref{S:bayes_num_res}. The time discretization step is $\Delta t = 0.1$ and the computational grid is a 60$\times$15 unstructured triangulation.

For testing the accuracy of the solution we consider the following indexes
\begin{equation}
\begin{array}{l}
E_{abs,n} = \|\UU_n-\UU_{anl,n}\|_2;\\
[4mm]
E_U = \sqrt{\dfrac{\sum_{n=1}^{N_f}\|\UU_n-\UU_{anl,n}\|_2}{\sum_{n=1}^{N_f}\|\UU_{anl,n}\|_2}};
\end{array}
\end{equation}
where $\UU_n$ and $\UU_{anl,n}$ are the computed velocity and discrete analytic velocity at time $t=t_n$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Non-periodic solution} In this test case we set $\VV_0 = \VV_{anl,0}$, the analytic solution at time $t$ = 0, and we solve the problem in $[0,\,\tau]$ with $\tau$ = 1. As an example, in Figure \ref{F:snapshots} we report data and assimilated velocity field for $t$ = 0, 0.3, 0.7.

In Table \ref{T:wom_non_per} we report $E_U$ for the deterministic solution and for the \gls{ML} estimator in correspondence of \gls{SNR} = 20 and 10; in the third column $\gamma$ is defined as $1-\frac{E_{U,{\rm ML}}}{E_{U,{\rm det}}}$. The choice of the \gls{ML} estimator is due to the fact that, since we perform interpolation, the problem is well-posed and no regularization (or smoothness priors, see \ref{S:point_estimators}) is required. In Figure \ref{F:wom_non_per} we report $E_{abs,n}$ at each time step for both deterministic and statistical solutions. From these results we can infer that the additional information brought by the knowledge of the noise yields more accurate results. Furthermore, we observe that the statistical solution features a smoother error behavior.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{c}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.8 \textwidth]{snapshot0.png} \\  
\includegraphics[width=0.8\textwidth]{snapshot03.png} \\
\includegraphics[width=0.8\textwidth]{snapshot07.png} 
\end{tabular}
\end{center}
\caption{Assimilated velocity field and noisy data (SNR = 10) at time $t$ = 0, 0.3, 0.7 of a Womersley simulation.}
\label{F:snapshots}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ | c | c | c | c |}
	\hline
	\gls{SNR}	&	$E_{U,det}$		&	$E_{U,{\rm ML}}$	&	$\gamma$	\\
	\hline
	20	&	0.0935				&	0.0633				&	32\%	\\
	10	&	0.1608				&	0.105				&	35\%	\\
	\hline
\end{tabular}
\caption{For different values fo SNR, we report comparison results for statistical and deterministic estimators.}
\label{T:wom_non_per}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.48\textwidth]{SNR20.pdf} &  
\includegraphics[width=0.48\textwidth]{SNR10.pdf} \\
\end{tabular}
\end{center}
\caption{$E_{abs,n}$ for $n = 1,\, ..., \, N_f$ in $[0,\,2\tau]$ for the deterministic and statistical solution for SNR = 20 and 10.}
\label{F:wom_non_per}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Periodic solution} In this case we want to compare the solution obtained using $\VV_0 = \VV_{anl,0}$ and $\VV_0 = {\bf 0}$. In Table \ref{T:wom_per} we report $E_{\tau}$ obtained after 3 iterations of algorithm \eqref{E:periodic_alg}; recall that we set the initial velocity equal to the final velocity of the previous period. For the same simulations we report in Figure \ref{F:wom_per} for 3 subsequent iterations values of the absolute errors $E_{abs,0}$ and $E_{abs,anl}$ corresponding, respectively, to the solution obtained with $\VV_0={\bf 0}$ and $\VV_0=\VV_{anl,0}$. This means that the assimilation is able to capture the solution in a very small amount of time (at most two iterations of algorithm \eqref{E:periodic_alg}). 

These results, although satisfactory and promising, are based on the assumption that $I_{meas}=I_{discr}$; this is definitely not realistic. Future work includes the implementation of the techniques proposed in Remark \ref{R:extension_unsteady}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ | c | c | c | c |}
	\hline
	\gls{SNR}	&	it	&	$E_U$, $\VV_0={\bf 0}$		&	$E_U$,$\VV_0=\VV_{anl,0}$	\\	\hline
	20	&	1	&	0.1318							&	0.09985							\\	
		&	2	&	0.09456							&	0.09456							\\	
		&	3	&	0.09314							&	0.09314							\\	\hline
	10	&	1	&	0.1928							&	0.1711							\\	
		&	2	&	0.1533							&	0.1533							\\	
		&	3	&	0.1512							&	0.1512							\\	\hline
\end{tabular}
\caption{For different values of SNR we report accuracy results obtained comparing approach \eqref{E:periodic_alg} with \eqref{E:unsteady_algorithm} using $\VV_0 = \VV_{anl,0}$.}
\label{T:wom_per}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.48\textwidth]{SNR20_per.pdf} &  
\includegraphics[width=0.48\textwidth]{SNR10_per.pdf} \\
\end{tabular}
\end{center}
\caption{$E_{abs,n}$ for $n = 1,\, ..., \, N_f$ in $[0,\,\tau]$ for $\VV_0 = {\bf 0}$ and $\VV_0 = \VV_{anl,n}$ for SNR = 20 and 10.}
\label{F:wom_per}
\end{figure}
%________________________________________________________________________________________
%\begin{figure}[t]
%\begin{minipage}[b]{0.5\linewidth} % A minipage that covers half the page
%\centering
%\includegraphics[scale=0.32]{60_12_2e3.png}
%\end{minipage}
%\hspace{0.25cm} % To get a little bit of space between the figures
%\begin{minipage}[b]{0.5\linewidth}
%\centering
%\includegraphics[scale=0.32]{60_12_1e3.png}
%\end{minipage}
%\caption{Singular Values of the sensitivity matrix J$_t$ on a $60\times12$ mesh with $\Delta t=2\;10^{-3}$ (left) and $\Delta t=1\;10^{-3}$ (right). }
%\label{F:sval_unst}
%\end{figure}

% ---------------------------------------------------------------------------
%: ----------------------- end of thesis sub-document ------------------------
% ---------------------------------------------------------------------------

