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

\ifpdf
    \graphicspath{{2/figures/PNG/}{2/figures/PDF/}{2/figures/}}
\else
    \graphicspath{{2/figures/EPS/}{2/figures/}}
\fi
%: ----------------------- name of chapter  -------------------------
\chapter{Preliminary analysis}\label{C:preliminary} % top level followed by section, subsection
In Chapter \ref{C:introduction} we briefly presented the wide literature on \gls{DA} for fluid geophysics and the state of the art of emerging techniques for hemodynamics. As we approach the problem of assimilating blood velocity from noisy measures, in this chapter we perform some preliminary tests on candidate methods with the purpose of selecting the most suitable in terms of accuracy and efficiency. Although limiting, for the sake of simplicity, we focus on deterministic approaches and we refer to Chapter \ref{C:da_bayesian} for a deeper investigation of stochastic methods.

We present our preliminary analysis, see \cite{md_av1}, conducted on the Stokes problem, a simplification of the \gls{NSE}, where the nonlinearity is removed. At this stage we relax several aspects of the original problem; in fact, the analysis is limited to the 2-dimensional steady case on ``toy'' problems with synthetic data. 

In Section \ref{S:preliminary_methods} we introduce three approaches to the assimilation ranging from interpolation techniques to least squares methods. We present their formulation adapted to hemodynamics equations and we discuss potential advantages and drawbacks we might incur in. In Section \ref{S:preliminary_results} we give details regarding numerical test cases, we compare the quality of numerical results and we motivate the choice of the assimilation method that we propose in this work. In Section \ref{S:preliminary_furtheranl} we present a more detailed error analysis for the selected method, still on the Stokes problem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Candidate methods}\label{S:preliminary_methods} 
In this section we present the assimilation problem for the Stokes model and we describe a typical data configuration for the problem we are interested in. We propose three approaches and we highlight advantages and potential weaknesses.

Let us denote by $\Omega\subset\mathbb{R}^2$ the domain of interest, i.e. the vessel; $\Gamma_{in}$ and $\Gamma_{out}$ the inflow and outflow sections; $\Gamma_{wall}$ the vessel wall. Let $L^2(\Omega)$ be the Hilbert space of square summable functions, ${\bf H}^1(\Omega)$ the space of vector functions whose components are in $L^2(\Omega)$ together with their first derivatives and, ${\bf H}_{\Gamma}^1$ the corresponding space of ${\bf H}^1$ functions with null trace on the portion $\Gamma$ of $\partial \Omega$ . We denote by $\uu(\xx) \in {\bf H}_{{\Gamma}_{wall}}^1(\Omega)$ and $p(\xx) \in L^2(\Omega)$ velocity and pressure fields respectively. We assume to have $N_s$ velocity measures $\dd$ available at some sites  $\xx_i^{m}\in \dom$, for $i=1,\,...,\,2N_s$; these measurement points lie on some internal layers, which we denote $\Gamma_{d}$ (see Figure \ref{domain}), with $\Gamma_{in} \subset \Gamma_d$. In this preliminary study, we assume sites to correspond to grid nodes, in Chapters \ref{C:da_nse}, \ref{C:da_bayesian} and \ref{C:da_unsteady} this assumption is relaxed and we consider sites and nodes belonging to different sets. Velocity and pressure are assumed to fulfill the Stokes equations 
\begin{equation}\label{E:steady_stokes}
\left\{\begin{array}{ll}
-\nu\Delta\uu  + \gr p=\mathbf{0}       					 & \textrm{in}\;\; \dom,\\
[3mm]
\gr\cdot\uu = 0                             								         & \textrm{in}\;\; \dom,\\
[3mm]
\uu = \mathbf{0} 																	 & \hbox{on}\;\;\Gamma_{wall},\\
[2mm]
\displaystyle-\nu\gr\uu \, \mathbf{n}+ p\,\mathbf{n}=\mathbf{h} & \textrm{on}\;\;\Gamma_{in},\\
[3mm]
\displaystyle - \nu \gr\uu\, \mathbf{n}+p\,\mathbf{n}=\mathbf{0}& \textrm{on}\;\;\Gamma_{out}.
\end{array}\right.
\end{equation}
Here, $\mathbf{h}(\xx)\in {\bf H}^{-1/2} (\Gamma_{in})$ is the Neumann data at the inflow boundary, assumed unknown; $\nu$ is the kinematic viscosity. We aim at solving the following problem: {\it find $\uu$ and $p$ so that $\uu$ fits the available measures $\dd$ (in some sense) under the constraints (\ref{E:steady_stokes}).}

As discussed in the introduction, different methods can be devised for solving this problem, with different ways of fitting the data. Here, we address in particular three techniques, which we refer to as {\it splitting}, {\it nudging} (or dynamic relaxation) and {\it control based} techniques. The former is specifically devised for this problem and relies upon the assumption that available measures lie on internal layers of the domain of interest, i.e. points $\xx_i^m$ belong to a set of transversal sections of the vascular district at hand, as done in \cite{heys2010} (see Figure \ref{domain}). Although this is a quite restrictive assumption, it reflects the way the measurements are collected and post processed (see Figure \ref{F:MRIsection}). The key idea of this approach is the interpolation of available data. The second method is based on a more general approach, that relies upon the theory of control for distributed systems. The third one is a nudging technique reinterpreted for our specific problem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Splitting techniques}\label{ss:splitting}
We split the domain into sub domains (see Figures \ref{F:MRIsection} and \ref{domain}),  where $\Gamma_{in}$, $\Gamma_{out}$ and the measurement sections act as domain boundaries. Measures are used as boundary data for solving the global problem by sub domains. More precisely, let us assume to have $s$ sub domains $\Omega_i$ ($i=0,\ldots,s-1$), such that data lie on the interfaces. Let us denote by $\Gamma_i$, with $i=0,\,\ldots,\,s$, the layers, so that $\Gamma_0=\Gamma_{in}$ and $\Gamma_s=\Gamma_{out}$. Let $\dd_i$, $i=0,\ldots,s-1$, be the set of measures on $\Gamma_i$, or more precisely an interpolating function (e.g. piecewise linear interpolation) over $\Gamma_i$ of the data. 
%We rewrite the velocity as $\uu = \sum_{i=0}^{s-1}\uu_{0,i} + {\cal E}_i \dd$ where $\uu \in {\bf H}^1_{{\Gamma}_{wall}} (\Omega)$, $\uu_{0,i} \in {\bf H}^1_{\partial {\Omega}_i} (\Omega_i)$, and ${\cal E}_i \dd\in {\bf H}^1(\Omega_i)$ is an arbitrary extension in $\Omega_i$ of the data $\dd_i$ and $\dd_{i+1}$. We want to solve the following problem. Find $\uu$ such that for any $\vv_i \in {\bf H}^1_{\partial {\Omega}_i} (\Omega_i)$ and $q_i \in L^2(\Omega_i)$
%\begin{equation}\label{eq:split}
%\sum_{i=0}^{s-1}\left(\; \int_{{\dom}_i} \nu \gr \uu \cdot \gr \vv_i \,{\rm d} \xx  - \int_{{\dom}_i} p \gr \cdot \vv_i \,{\rm d} \xx -
%\int_{{\dom}_i} q_i \gr \cdot \uu \,{\rm d} \xx \right) = 0.
%\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\figuremacroW{domain}{Domain splitting}{Schematic representation of the domain with velocity measures (left) and its sub domain splitting (right).}{1.}
In principle, we can distinguish two approaches, according to the sequences {\it split then discretize} or {\it discretize then split}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Discretize then split: matrix updating (MU)}\label{p:matrix_updating} 
In this case, using P1bubble-P1 \gls{FE} spaces for velocity and pressure (inf-sup compatible spaces, see e.g. \cite{quarteroni_numerical_1994}, chap. 9, 10) respectively, we first discretize problem \eqref{E:steady_stokes} to the form ${\rm S} \VV = {\bf F}$, where $\VV$ is the vector of the nodal unknowns; ${\rm S}$ is the Stokes matrix and ${\bf F}$ comes from the prescription of boundary conditions. Then, the measurements are regarded as conditions in the form
\begin{equation}\label{eq:new1}
\VV(\mathbf{x}_i)= {\bf d}_i
\end{equation}
where $i$ ranges over the nodes where measures are available. Conditions (\ref{eq:new1}) are forced into the linear system similarly to Dirichlet conditions (see \cite{quarteroni_numerical_models}, chapter 11), i.e. $i$-th equations of the linear system corresponding to measurements points are replaced by (\ref{eq:new1}). In practice, we set the off-diagonal entries of the $i$-th rows to 0 and the diagonal entry to 1. The corresponding entries of the right hand side, ${\bf F}_i$ are set equal to ${\bf d}_i$. Then, ${\bf F}_j$, $j\neq i$ are updated accordingly, by setting ${\bf F}_j = {\bf F}_j - \sum_i [{\rm S}]_{ji} {\bf d}_i$, see \cite{quarteroni_numerical_models}. 
%This perturbation is expected to affect the condition number of the system. The coefficient $c$ is tuned to minimize this effect. In Table \ref{T:conditioning} we report the condition number with the choice of $c=\Delta$ on two different meshes, being $\Delta$ the mesh size. After the update the condition number increases of about one order of magnitude (setting $c=1$, $c=\Delta^2$ or $\Delta^3$ yields same results). {\it This fact has the clear drawback of requiring a higher computational effort in solving the linear systems. Even more important is the fact that the solution of the system becomes more sensitive to the perturbations on the data and in particular to the noise in the measures}.{\tt to be rephrased}. Despite of its simplicity, this approach is therefore anticipated to be extremely sensitive to the noise, see Section \ref{S:preliminary_results}.

Clearly, in each sub domain, the incompressibility condition might be violated since data at the inflow and outflow boundary might not fulfill the compatibility condition
\begin{equation}\label{E:compatibility}
\int_{\Gamma_i} \dd_i \cdot {\bf n}\, {\rm d}\gamma =  \int_{\Gamma_{i-1}} \dd_{i-1} \cdot {\bf n}\, {\rm d}\gamma.
\end{equation}
A possible way to avoid this problem, which has not been investigated in this work, is to ``adjust'' properly the data so to obtain \eqref{E:compatibility}, as a pre-processing step.
Nevertheless, the updated discretization matrix is not singular (though its conditioning is higher than the original one). This can be explained by the fact that the incompressibility condition is forced in a weak sense in the entire domain, therefore on each sub domain it can be violated. 
%Note that a different choice of \gls{FE} spaces might lead to non-singular matrices\footnote{This is the case, as an example, of Raviart-Thomas \gls{FE} spaces, which guarantee the mass conservation in each element \cite{FEbook}.}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{table}[t]
%\begin{center}
%\begin{tabular}{| c | c | c|}
%        \hline
%        mesh & no update & $c=\Delta$\\
%        \hline
%        40$\times$ 10   &   2.58e3 &  1.25e4  \\
%        \hline
%        80$\times$ 20   &   1.05e4 &  8.30e4   \\
%        \hline
%\end{tabular}
%\caption{1-norm condition number estimate of S.}
%\label{T:conditioning}
%\end{center}
%\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Split then discretize: domain splitting}\label{p:domain_splitting} In this case we first split the differential problem in $s$ sub domains $\dom_i$, $i = 0,\,...,\,s-1$ and we consider each of them separately. 
%Formally, the subproblem on $\Gamma_{i}$ corresponds to setting in \eqref{eq:split} $\vv_j = {\bf 0}$ and $q_j = 0$, for all $j = 0,\,...,\,s-1$, $j\neq i$. 
Then, we perform the discretization of each problem; in the first $s-1$ sub domains we have Dirichlet conditions over the entire $\partial \Omega_i$, in $\dom_{s-1}$ we have a homogeneous Neumann condition on $\Gamma_s = \Gamma_{out}$. Note that, while in $\dom_{s-1}$ the pressure is uniquely determined, in $\dom_i$, $i = 0,\,...,\,s-2$, it is determined up to a constant (since we have a Dirichlet problem). We set such constant prescribing continuity for the pressure in $\dom_i$ with respect to the one in $\dom_{i+1}$ for $i = 0,\,...,\,s-2$.
%so that the pressure is determined up to a constant. We select a specific pressure field by taking a null average pressure or fixing the pressure in one node. In principle, this method is fast and prone to parallelism, since each subproblem is solved independently of the others. 

However, the presence of the noise has three important drawbacks: 
\begin{itemize}
\item measurements errors affect the solution inside each domain, without filtering;
\item the stress is not continuous across the sections;
\item the divergence compatibility condition over the whole domain is violated.
\end{itemize}

Concerning the first point, it is known that perturbations on the boundary affect the solution of a Stokes problem with a decay featuring an exponential dependence on the distance from the boundary (see \cite{VV}). Therefore, the presence of noise on the interfaces brings errors inside the domain, yielding inaccurate results. For this reason we do not dwell any longer with this approach. As a matter of fact, preliminary results show that the numerical error with respect to reference solutions is one order of magnitude larger than for other methods, see \cite{md_av1}.
%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Control based methods}\label{ss:control}
In formulating control approaches we introduce the weak form of problem \eqref{E:steady_stokes}: for all $\vv \in {\bf H}_{{\Gamma}_w }^1(\Omega), q \in L^2(\Omega)$ velocity and pressure satisfy
\begin{equation}\label{E:steady_stokes2}
\int_{\dom} \nu \gr \uu \cdot \gr \vv\,{\rm d} \xx  - \int_{\dom} p \gr \cdot \vv\,{\rm d} \xx -
\int_{\dom} q \gr \cdot \uu\,{\rm d} \xx +
\int_{{\Gamma}_{in}} {\bf h}\,\vv\,{\rm} d \gamma = 0.
\end{equation}
We aim at finding the optimal boundary condition $\hv$, treated as control variable, such that some distance (to be defined) between velocity and data is minimal. The choice of $\hv$ can be motivated by the medical application and will be discussed in Chapter \ref{C:da_nse}. According to the theory presented in Chapter \ref{C:introduction}, we consider the \textit{optimize then discretize} and the \textit{discretize then optimize} approaches.
%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Optimize then discretize (OD)} \label{p:od} We define the distance between solution and data as
\begin{equation}\label{E:functional_od}
dist(\uu, \, \dd)=\frac{1}{2}\int_{\Gamma_{d}}\big(\uu-\dd\big)^2\,\vv\,\textrm{d}\gamma,
\end{equation}
for all $\vv\in {\bf H}^1_{\Gamma_{wall}}(\dom)$. Here, with an abuse of notation, $\dd$ is an extension (obtained, for instance, by interpolation) of the vector data to a continuous functional space. We solve the problem of minimizing the functional
\begin{equation}
\mathcal{J}(\uu, \,\hv) = dist(\uu, \, \dd);
\end{equation}
under the constraint of \eqref{E:steady_stokes}. The Lagrangian of the problem reads
\begin{equation}\label{E:lagr_od}
\begin{array}{l l}
\mathcal{L}(\uu,p,\boldsymbol{\lambda}_{\uu},\lambda_{p},\hv)= & \mathcal{J}(\uu,\, \hv)+ \int_{\dom}\nu\gr\uu\cdot\gr\boldsymbol{\lambda}_{\uu}\textrm{d}\xx\,+\\
& -\int_{\dom}p\gr\cdot\boldsymbol{\lambda}_{\uu}\textrm{d}\xx
-\int_{\Gamma_{in}}\hv\boldsymbol{\lambda}_{\uu}\textrm{d}\gamma
+\int_{\dom}\lambda_{p}\gr\cdot\uu\,\textrm{d}\xx;
\end{array}
\end{equation}
where $\boldsymbol{\lambda}_{\uu}$ and ${\lambda}_{p}$ are the Lagrange multipliers associated with velocity and pressure. According to \cite{gunzburger_perspectives_2002}, we derive the system of necessary conditions for optimality
\begin{equation}\label{E:od_necessary_cond}
\left\{\begin{array}{l}
\dfrac{\partial \mathcal{L}}{\partial[\boldsymbol{\lambda}_{\uu},\lambda_p]} = 
\int_{\dom}\nu\gr\uu\cdot\gr\vv\textrm{d}\xx
-\int_{\dom}p\;\gr\cdot\vv\textrm{d}\xx
-\int_{\Gamma_{in}}\hv\vv\textrm{d}\gamma
+\int{\dom}q\;\gr\cdot\uu\;\textrm{d}\xx=0, \\
[3mm]
\dfrac{\partial \mathcal{L}}{\partial[\uu,p]} = 
\int_{\Gamma_{d}}\left(\uu-\dd\right)\ww\textrm{d}\gamma+
\int_{\dom}\nu\gr\cdot\boldsymbol{\lambda}_{\uu}\gr\ww\textrm{d}\xx
+\int_{\dom}\lambda_{p}\gr\cdot\ww\textrm{d}\xx
-\int_{\dom}r\;\gr\cdot\boldsymbol{\lambda}_{u}\textrm{d}\xx=0, \\
[3mm]
\dfrac{\partial \mathcal{L}}{\partial \hv} =
\int_{\Gamma_{in}}\boldsymbol{\lambda}_{\uu}{\bf s}\;\textrm{d}\gamma=0.
\end{array}
\right.
\end{equation}
These conditions are intended in the weak sense for all $\vv\in{\bf H}_{\Gamma_{w}}^{1},\,q\in L^2, \,\ww\in{\bf H}_{\Gamma_w}^1,\,r\in L^2$ and ${\bf s}\in{\bf H}^{1/2}(\Gamma_{in})$ respectively. In the discretization step we rely on the finite element (FE) method. The discretization of \eqref{E:od_necessary_cond} yields the linear system

\begin{equation}\label{E:linsys_od}
\left[\begin{array}{l l l}
\textrm{S} & \textrm{O} & \textrm{N} \\
\textrm{J} & \textrm{S} & \textrm{O} \\
\textrm{O}  & \textrm{N}^{\TT} & \textrm{O} \\
\end{array}\right]\;
\left[\begin{array}{l}
\mathbf{V}\\
\mathbf{\Lambda}\\
\mathbf{H}\\
\end{array}\right]\;=\mathbf{F} 
\end{equation}
where
\begin{equation}\label{E:notation_od}
\textrm{S}=\left[\begin{array}{ll}
\textrm{C} & \textrm{B}^{\TT}\\
\textrm{B} & \textrm{O} 
\end{array}\right], 
\mathbf{V}=
\left[\begin{array}{l}
\mathbf{U}\\
\mathbf{P}\\
\end{array}\right],
\mathbf{\Lambda} =\left[\begin{array}{l}
\mathbf{U}_{\lambda}\\
\mathbf{P}_{\lambda}
\end{array}\right].
\end{equation}
Here, $\UU\in\mathbb{R}^{N_u}$ and $\PP\in\mathbb{R}^{N_p}$  are the discretization of velocity and pressure belonging to inf-sup compatible FE spaces (see e.g. \cite{quarteroni_numerical_1994}, chap. 9, 10), $\HH\in\mathbb{R}^{N_{in}}$ is the discretization of the control variable $\hh$, $N_{u}$ and $N_p$ are the number of degrees of freedom (\gls{DOF}s) of velocity and pressure and $N_{in}$ is the number of \gls{DOF}s of the velocity on $\Gamma_{in}$. 
Let ${\rm N} = {\rm R}_{in,u}^{\TT}{\rm M}_{in}$, where R$_{in,u}\in\mathbb{R}^{N_{in},N_u}$ is a \textit{restriction matrix} which selects the \gls{DOF}s of the velocity on $\Gamma_{in}$. The generalization of R$_{in,u}$ to velocity and pressure \gls{DOF}s reads
\begin{equation}\label{E:matrixR}
\RR_{in}=\left[\RR_{in,u}\;\;\OO\right], \quad R_{in} \in \mathbb{R}^{N_{in},N_u+N_p}.
\end{equation}

Let ${\rm M}_{in} \in \mathbb{R}^{N_{in},N_{in}}$ be the discretization of the mass operator restricted to inlet boundary nodes, specifically
\begin{equation}\label{E:matrixM_FE}
\begin{array}{llll}
\left[{\rm M}_{in}\right]_{kl}& = &\int\limits_{\Gamma_{in}}\displaystyle\left.\boldsymbol{\varphi}_{j_l}\right|_{\Gamma_{in}} \left.\boldsymbol{\varphi}_{i_k}\right|_{\Gamma_{in}} \, {\rm d \sigma}      & \quad k,\,l = 1,\,...,\,N_{in},
\end{array}
\end{equation}
where $\boldsymbol{\varphi}$ are the Lagrangian basis functions associated with the velocity and $\{i_k\}$ and $\{j_l\}$ are subsets of the velocity indexes, selecting \gls{DOF}s on $\Gamma_{in}$.

Finally, let ${\rm C}\in \mathbb{R}^{N_u,N_u}$, and ${\rm B}\in \mathbb{R}^{N_p,N_u}$ be the discretization of the diffusion and divergence operators. 
\begin{equation}\label{E:matrices_FE}
\begin{array}{llll}
\left[{\rm C}\right]_{ij} & = & \int\limits_{\dom}\nu \gr\boldsymbol{\varphi}_j \gr\boldsymbol{\varphi}_i \, {\rm d}\xx,  &  \quad i,\,j = 1,\,...,\,N_u, \\
\left[{\rm B}\right]_{ij} & =  & \int\limits_{\dom}\gr\cdot\boldsymbol{\varphi}_j \,\psi_i \, {\rm d}\xx,   & \quad  i = 1,\,...,\,N_p, \;\;  j = 1,\,...,\,N_u, \\
\end{array}
\end{equation}
where $\psi$ are the Lagrangian FE basis functions associated with the pressure. 

$\UU_{\lambda}$ and $\PP_{\lambda}$ represent the discrete Lagrangian multipliers; J and $\mathbf{F}$ account for the discretization of the term $\int_{\Gamma_{d}}\left(\uu-\dd\right)\ww\textrm{d}\gamma$.
%%%%%%%%%%%%%%%%%%%%%
\begin{remark}\label{L:regularization}
As described in Section \ref{ss:introduction_methods_control}, when using control approaches a regularization term may be added to the functional to prevent ill-posedness and improve the convexity of the problem. For reasons discussed in Chapter \ref{C:da_nse}, since the data vector is projected into a continuous functional space, the optimization problem has a unique solution; hence, here we drop the regularization term. For the same reasons we do not consider any regularization also in the discretize then optimize case.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%
\paragraph{Discretize then optimize (DO)}\label{p:do}
In this case, using the FE method we first discretize the functional and the constraints; the result is an algebraic optimization problem with linear constraints:
\begin{equation}\label{E:min_do}
\begin{array}{l}
\min\limits_{\HH} J(\VV,\,\HH) = \displaystyle\frac{1}{2}\|{\rm D}\VV - \dd\|^2_2\\
[3mm]
\hbox{s.t.} \;\; {\rm S}\VV = {\rm N} \HH;
\end{array}
\end{equation}
where the notation is the same of the previous paragraph. We let ${\rm Q} \in\mathbb{R}^{2N_s,N_u}$ be a matrix designed such that $[{\rm Q \UU}]_i$ is the numerical solution evaluated at the data sites and corresponding to $[\dd]_i$, for $i=1,\,...,\,2\, N_s$. We extend ${\rm Q}$ to pressure \gls{DOF}s introducing the matrix ${\rm D} = [{\rm Q}\;\;\OO]$. For the solution of this minimization problem, we write the Lagrangian
\begin{equation}\label{E:do_lagr}
L(\VV,\HH,\boldsymbol{\Lambda})=\displaystyle\frac{1}{2}\|{\rm D}\VV-\dd\|_2^2 + \boldsymbol{\Lambda}^{\TT}({\rm S}\VV-{\rm N}\HH),
\end{equation}
and we derive the induced system of optimality conditions, also known as KKT system \cite{nocedal_numerical_2000}
\begin{equation}\label{E:do_lagr_necessary}
\left\{\begin{array}{l}
\dfrac{\partial L}{\partial\VV} = {\rm D}^\TT({\rm D}\VV-\dd)+{\rm S}^\TT\boldsymbol{\Lambda}=\mathbf{0}\\
[3mm]
\dfrac{\partial L}{\partial\boldsymbol{\Lambda}} = {\rm S}\VV-{\rm N}\HH =\mathbf{0}\\
[3mm]
\dfrac{\partial L}{\partial\HH} = -{\rm N}^{\TT}\boldsymbol{\Lambda}=\mathbf{0}.
\end{array}\right.
\end{equation}
Upon block elimination, we get the reduced system
\begin{equation}\label{E:do_control}
{\rm Z}^{\TT}{\rm Z} = {\rm Z}^{\TT}\dd;
\end{equation}
where ${\rm Z}={\rm D}{\rm S}^{-1}{\rm N}$ is the so called \textit{sensitivity} matrix. It is defined as ${\rm Z} = \frac{{\rm d} (\textrm{D}\VV)} {{\rm d} \HH}$ and it is a measure of how much perturbations of the control variable may affect the output ${\rm D}\VV$. For this reason, spectral properties of ${\rm Z}$ determine the conditioning or well-posedness of the problem. For insights on efficient techniques for the solution of problem \eqref{E:do_lagr_necessary} see \cite{haber_optimization_2000}. 

\subsection{Dynamic relaxation}\label{SS:preliminary_methods_nudging}
We propose a reinterpretation of the nudging method, or dynamic relaxation (\gls{DR}), in a weak sense and we discuss a possible choice of the relaxation parameter. In fluid geophysics the initial condition of the forecast is the final solution of the pre-forecast period (Section \ref{SS:introduction_geophysics_nudging}). Here, we consider the solution of the steady problem, which we are interested in, as the stationary state of an unsteady problem with artificial terms enforcing data fitting. The initial condition is chosen arbitrarily, but in such a way that it is physically consistent. The formulation of the method in a weak form makes the problem more amenable to a FE discretization. 

We solve the problem of finding, for each test function $\vv(\xx,t)$ and $q(\xx,t)$ in the velocity and pressure spaces, $\uu(\xx,t)$ and $p(\xx,t)$ such that
\begin{equation}\label{E:stokes_dr}
\left\{\begin{array}{l}
\int_{\dom}\uu_t \;\vv\textrm{d}\xx
+\int_{\dom}\nu\gr\uu\gr\vv\textrm{d}\xx 
-\int_{\dom}p\;\gr\vv\textrm{d}\xx
+\int_{\dom}q\;\gr\cdot\uu\textrm{d}\xx\\
-\int_{\Gamma_{in}}\left(\nu \,\gr\uu \, \mathbf{n}+ p\,\mathbf{n}\right)\vv\textrm{d}\gamma +c\;\int_{\Gamma_{in}\cup\Gamma_{d}}(\uu-\dd)\vv\textrm{d}\gamma=0\\
[3mm]
\uu(\xx,t=0) = \uu_0.
\end{array}\right.
\end{equation}
Here, the forcing term is a boundary inner product, between the test function and the difference of data and predicted velocity, evaluated on measurement layers. We note that this term embeds a weak prescription of boundary conditions on $\Gamma_{in}$; for this reason the term $\int_{\Gamma_{in}}\left(\nu \gr\uu \,\mathbf{n}+ p\,\mathbf{n}\right)\vv\textrm{d}\gamma $ is added to the formulation. This fact suggests a way for the selection of the nudging coefficient, which can be thought of as a penalty parameter. We choose $c$ according to the theory of artificial boundary conditions presented in \cite{burman_edge_2006}. In the assumption of FE discretization, $c$ is proportional to $\mathcal{D}_{max}$, the largest triangle diameter, $c = k \nu \mathcal{D}_{max}$, where $k$ is a positive parameter and $\nu$ the viscosity. 

This approach is clearly related to the more sophisticated weighted least squares FE (WLSFE) method for \gls{DA}, introduced in \cite{heys2010} by Heys et al., already mentioned in Chapter \ref{C:introduction}. Such method, designed for steady problems, consists in the minimization of a functional given by the sum of the weak formulation of the state equations, weak prescription of boundary conditions and weighted norm of the misfit between data and velocity. The weights are chosen according to a measure of trust in the data. The substantial difference, with respect to nudging, is in setting the method as an optimization problem. This fact allows to account for the noise without altering the nature of the problem since no artificial forcing terms are added to the state equations. Details on the WLSFE can be found in Chapter \ref{C:introduction} and in \cite{heys2010}.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical results} \label{S:preliminary_results}
We report numerical results obtained testing the proposed approaches on a toy problem with synthetic data. We consider the Poiseuille flow, a steady flow governed by the Stokes equations in the rectangular domain $\dom\subset\mathbb{R}^2$, $\dom=[0,L]\times[-R,R]=[0,5]\times[-0.5,0.5]$. The exact solution reads
\begin{equation}\label{E:poiseuille_preliminary}
\left\{\begin{array}{l}
\uu_{P}(x,y)=\; [1-4y^2\;,\; 0\;]^{\TT},\\
p_P(x,y)=-8\nu (x-L);
\end{array}\right.
\end{equation} 
Figure \ref{vel1and5} reports the horizontal velocity field computed by the FE solver {\tt FreeFem++}\footnote{An open source C++ software for the numerical solution of two dimensional partial differential equations with the FE method. http://www.freefem.org.} on a fine unstructured grid.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\figuremacroW{vel1and5}{Poiseuille flow}{Horizontal velocity profile, vertical lines in correspondence of $x\in\{0,\;1,\;2.5\;4\}$ show data location.}{.85}
%\paragraph{Data generation} 
According to the problem presented in the introduction, we assume data to be given on three internal layers and on the inflow boundary in correspondence of discretization nodes, i.e. the space discretization step $\Delta$ satisfies $\Delta \propto N_s^{-1}$. Precisely, in this configuration $\Gamma_d = \left\{(x,y)\;|\;x\in\{0,\;1,\;2.5,\;4\},\;\;y\in[-R,R]\;\right\}$. Data are generated adding to the exact solution random noise, normally distributed with zero mean and variance determined by the prescription of the signal-to-noise ratio (\gls{SNR}) defined as the ratio between the peak velocity and the level of noise.

Computational grids are unstructured triangulations generated with {\tt FreeFem++}, see Figure \ref{meshFE}.  We implement the FE method with choice of compatible FE spaces P1bubble-P1 for the velocity and pressure fields respectively (see e.g. \cite{quarteroni_numerical_1994}, chap. 9, 10).
%\paragraph{Software details} This preliminary comparison is performed combining two softwares. For the matrix assembling we rely on {\tt FreeFem++}, this software provides several built-in functions for the evaluation of the variational forms induced by the FE formulation. Weak points are the storage of sparse matrices, the treatment of boundary conditions and the solution of linear systems; for these reasons, we use {\tt Matlab} for matrix operations and system solvers. For \gls{DR} simulations we rely on {\tt FreeFem++} only.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\figuremacroW{meshFE}{Computational grid}{Triangulation (80$\times$ 20) generated with {\tt FreeFem++}.}{.85}
%\paragraph{Computational details} In the case of the DO method we choose the gradient as a smoothing regularizing operator and we generate the optimal parameter by means of \gls{DP} and \gls{GCV} methods \cite{Hansen}. For reasons that will be clear in Chapter \ref{C:da_nse}, the regularization is not needed; in fact, data located on grid nodes on $\Gamma_{in}$ ensure the well-posedness of the inverse problem. As a matter of fact, with both \gls{DP} and \gls{GCV} the regularization parameter is of the order of $10^{-9}$.
At this preliminary stage it is important to answer the question {\it ``is \gls{DA} worthwhile?''}. In fact, in principle, internal measures could be disregarded and data on $\Gamma_{in}$, properly extended to grid nodes, could be used as Dirichlet condition and this would lead to a (cheap) well-posed problem. Mathematically and physically \gls{DA} is always preferable since it takes into account the mathematical model and produces feasible solutions. On the other hand, at a practical level, the issue is to find a good compromise between accuracy and efficiency. For these reasons we compare the performance of \gls{MU}, \gls{OD}, \gls{DO} and \gls{DR} with a direct approach, say Dirichlet (D), where the Stokes equations are solved prescribing the noisy velocity data on the inflow as boundary condition. Note that the \gls{DR} method deserves a separate treatment; in fact, tuning the time interval, the constant $k$ and the discretization time-step different levels of accuracy can be achieved.

We test the methods over a set of noise realizations and we consider the mean relative error as accuracy index
\begin{equation}
\overline{E}_U=\sum\limits_{i=1}^{N_r}\displaystyle\frac{\|\UU_i-\UU_{P}\|_2}{\|\UU_{P}\|_2}, 
\end{equation}
being $\UU_{P}$ the discretization of $\uu_P$ and $\UU_i$ the solution computed in correspondence of the $i$-th noise generation, out of $N_r$ = 16 with \gls{SNR} = 20 and 8; while the first value is chosen for testing the methods, the second is a reasonable \gls{SNR} for practical applications, e.g. for measures acquired via phase velocity encoded MRIs \cite{brummer}. 

For \gls{DR} simulations we set $c=k\,\nu\, \mathcal{D}_{max}$ with $k=10^4$, $\nu = 0.035$, a value featured by blood, $\Delta t=10^{-2}$, $\uu_0(\xx)=10^{-5}\uu_{P}(\xx)$. Here, $\Delta t$ is the time-step used in the time discretization scheme; in this case we use the Backward-Euler (\gls{BE}) method \cite{qss}. 
Accuracy results are reported in Figure \ref{F:comp}; for each method we report the relative error $\overline{E}_U$. The diagram on the right also displays the error associated with the \gls{DR} method after 150 time steps. Even if much cheaper in terms of computational effort, the \gls{MU} method is not accurate enough; this is due to the sensitivity to the noise, which makes this approach unreliable. The \gls{DR} method has a poor accuracy and, furthermore, requires the solution of an unsteady problem, which makes it computationally more expensive.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.48\textwidth]{comp20.pdf} &  
\includegraphics[width=0.48\textwidth]{comp8.pdf} \\
\end{tabular}
\end{center}
\caption{Comparison results for the OD, DO, D, MU and DR approaches in terms of relative error $\overline{E}_U$ with SNR = 20 (left) and 8 (right).}
\label{F:comp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Table \ref{T:DRresults} reports results obtained for \gls{SNR} = 8 with different number of iterations. In Figure \ref{F:DRgraph_solutions} we report the horizontal velocity at different iterations of the process: during the transient we notice that we do not have a physical solution due to the prescription of the data at the internal layers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | c |}
    \hline
	\verb"iter"			& $E_{U}$\\
	\hline
	50				& 1.416e-1\\
	\hline
	100				& 1.249e-1\\
    \hline
	150	 			& 1.225e-1\\  
	\hline
\end{tabular}
\caption{Relative errors in correspondence of SNR = 8 for the DR technique with different number of iterations.}
\label{T:DRresults}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{1\linewidth} 
\centering
\includegraphics[scale=1]{vel1e-2.png}
\end{minipage}
\hspace{0cm} 
\begin{minipage}[b]{1\linewidth}
\centering
\includegraphics[scale=1]{vel2e-1.png}
\end{minipage}
\hspace{0cm} 
\begin{minipage}[b]{1\linewidth}
\centering
\includegraphics[scale=1]{vel1and5.png}
\end{minipage}
\caption{Numerical solutions for the horizontal velocity field obtained at different iterations, specifically, $0$, $20$ and $150$ of the DR process.}
\label{F:DRgraph_solutions}
\end{figure} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
These results show that even with 150 steps the \gls{DR} approach is much less accurate than other methods. While \gls{DR} does not require an adjoint operator computation, the high number of iterations and the poor accuracy make this method not competitive. We note that the choice of $c$ has not been investigated in this work; we do not exclude that an optimal choice of this parameter might yield more accurate results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Selection of the DA method}
The four candidate methods, \gls{MU}, \gls{OD}, \gls{DO} and \gls{DR}, present significant differences in terms of accuracy and time saving. 

The \gls{MU} method is computationally cheap, costs a forward solve, but its accuracy is poor. In fact, the reference method, Dirichlet approach, can achieve more accurate results. This is somehow expected since data are incorporated without taking care of the presence of the noise and no filtering mechanism is active. With the same purpose of using measures as boundary data, a valid alternative would be a weighted weak prescription of boundary conditions where weights depend on how much we trust the measures. This approach results in a modified weak formulation, with additional terms of the form
\begin{equation}\label{E:weak_prescription}
\int_{\Gamma_d}{\bf w}\,(\uu-\dd)\,\vv \, {\rm d}\gamma
\end{equation}
where ${\bf w}$ is the weight function. Methods of this sort can be considered a combination of splitting and nudging techniques. A sophisticated formulation of this approach is the WLSFE approach \cite{heys2010}, already mentioned in Section \ref{SS:preliminary_methods_nudging}. Here, boundary penalty terms as in \eqref{E:weak_prescription} are added to the weak formulation; the minimizer of the resulting functional is the assimilated solution. Modification, improvement, and comparison of splitting techniques is part of future work. Nevertheless, it is important to recall that these approaches have a limiting feature; in fact, they are suitable only for problems where data are located on internal slices and do not allow for generalization to sparse data.

On the other hand, \gls{DR} is not reliable in terms of accuracy and requires the solution of a time-dependent problem, which makes the method computationally expensive. Nevertheless, these facts qualify \gls{DR} to be considered as a candidate method for approaching the assimilation problem for time dependent data.

Control approaches are the most accurate; however they require higher computational effort. In average, the computational time required by \gls{OD} is 1.4 times larger than \gls{DO}'s. With the purpose of improving computational aspects we consider the \gls{DO} approach as the most promising; in fact, we believe that accuracy has a primary importance in biomedical applications. Besides that, current techniques for optimization problems allow fast and efficient solution of inverse problems \cite{haber_optimization_2000}; improving the computational efficiency of the \gls{DO} approach is part of the current work and it is fundamental in treating real data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Further analysis of the DO method}\label{S:preliminary_furtheranl}
We conduct a deeper investigation of the \gls{DO} approach applied to the Stokes equations. In particular, we analyze the behavior of the error with respect to the analytic solution (still the Poiseuille flow) versus the discretization error, the number of data and the \gls{SNR}. We discuss implementation issues and we report numerical results which highlight the promising filtering properties of our approach. At this stage we utilize \href{www.lifev.org}{{\tt lifeV}}\footnote{\tt www.lifev.org.}, an efficient and flexible open source C++ FE library. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The software}\label{S:preliminary_software}
The FE library {\tt lifeV} is a software for scientific computing, whose development started in 1999 as a joint work project involving Ecole Polytechnique F\'ed\'erale de Lausanne (CMCS) in Switzerland, Politecnico di Milano (MOX) in Italy, INRIA (REO)
in France. Since 2008 Emory University (Department of Mathematics and Computer Science) in Georgia, US, has become part of the developers community. {\tt lifeV} is a free software, subject to the LGPL license, and it is compilable on Unix-like systems. This software implements in C++ language algorithms and data structures to approximate numerically systems of \gls{PDE}s, that arise from mathematical modeling applied to fluid-dynamics, fluid-structure interaction, flow in porous media and electrocardiology. {\tt lifeV} has been employed in joint collaborations with medical and industrial partners (for instance Emory University, School of Medicine (GA
- USA); Policlinico di Milano (Italy). The library is composed of a core part which provides an abstract framework for the implementation of Galerkin FE methods, and of a set of solvers for specific applications. Although the library was originated as a serial software, recently the developers community implemented a parallel version of the library (in this work we use the serial version). Moreover, the software interfaces with external libraries like \href{http://www.boost.org/}{Boost} for containers, basic algorithms and smart pointers, \href{http://www.netlib.org/}{BLAS and LAPACK} for standard linear algebra tools, \href{http://www.cs.sandia.gov/CRF/aztec1.html}{Aztec} or \href{http://trilinos.sandia.gov/}{Trilinos} for matrix manipulation and linear systems solution. Results obtained by the solvers can be post processed using {\tt lifefilters} methods which take care of writing results on files, in a format compatible with common visualization software like \href{http://www.paraview.org/}{Paraview}\footnote{http://www.paraview.org}, \href{http://www.ensight.com/}{Ensight}\footnote{http://www.ensight.com/} or \href{http://www.ann.jussieu.fr/˜frey/publications/RT-0253.pdf}{Medit}\footnote{http://www.ann.jussieu.fr/˜frey/publications/RT-0253.pdf}.
% -------------------
\subsection{Implementation issues}\label{S:preliminary_DOimplementation}
The core of the data assimilation solver is the solution of the reduced system \eqref{E:do_control}. The bottleneck of the computations is the Stokes system, corresponding to either ${\rm S}^{-1}$ or ${\rm S}^{-\TT}$; for its solution we do not consider a monolithic solver, instead, we perform an exact algebraic factorization and we split the original matrix into two block-triangular matrices. We rewrite the matrix S as 
\begin{equation}\label{E:exact_facto}
{\rm S}=
\left[\begin{array}{c c}
{\rm C} & {\rm B}^\TT\\
{\rm B} & {\rm O}
\end{array}\right]=
\left[\begin{array}{c c}
{\rm C} & {\rm O}\\
{\rm B} & -{\rm BC}^{-1}{\rm B}^\TT
\end{array}\right]\,
\left[\begin{array}{c c}
{\rm I} & {\rm C}^{-1}{\rm B}^\TT\\
{\rm O} & {\rm I}
\end{array}\right].
\end{equation}  
 In this way we can solve separately the two following block triangular systems:
\begin{align}	
	& \left\{
	\begin{array}{l}
	\rm C \widetilde{\UU} = \mathbf{b}_1\\
	\rm B \widetilde{\UU}- {\rm BC}^{-1}{\rm B}^\TT\widetilde{\PP} = \mathbf{b}_2;
	\end{array}
	\right.   \label{E:block_tria1}  
	\\
	& \left\{
	\begin{array}{l}
	\rm C\UU = \mathbf{b}_1 - {\rm B}^\TT \mathbf{P}\\
	 \PP=\widetilde{\PP};
	\end{array}\right. \label{E:block_tria2}
\end{align}
where $\widetilde{\PP}$ and $\widetilde{\UU}$ are two auxiliary variables; these systems are solved with iterative approaches, namely \gls{GMRES}-type solvers. 
%Computational difficulties may arise in the solution of the block $\textrm{BC}^{-1}\rm B^T$, the Schur complement. A common technique for preconditioning such block for the Stokes equations is to consider the pressure mass matrix; this is known as \textit{Uzawa} preconditioner. It is based on the spectral equivalence to the inverse of the Schur complement. As opposed to the linearized \gls{NSE}, where we have an additional advection term, and the \gls{NSE}, with the nonlinear term, the Stokes model presents less troublesome numerical issues (e.g. ill-conditioning). For this reason we focus on solver parameters and preconditioners in Chapter \ref{C:da_nse}.
% ---------------------
\subsection{Numerical results}\label{S:preliminary_furtherresults}
Using {\tt FreeFem++} we generate unstructured triangulations of dimension $n$ (number of edges on the longest side of $\dom$), we set $n$ = 40, 60, 80, 120, 160. We do not add any regularization term to the formulation; also in this case, we assume to have measures on each grid node on $\Gamma_{in}$, this assumption ensures the existence and uniqueness of a minimizer for the problem as we will see Chapter \ref{C:da_nse}.

In Figure \ref{F:stokes1} we report, for $\Delta=1/40$, the assimilated pressure and velocity fields with noise-free data; in correspondence of the internal layers the velocity matches the data (the exact solution). Figure \ref{F:stokes2} reports computed velocity and noisy data with \gls{SNR} = 10. The noise mainly affects the vertical components of the field (zero in the analytic solution); in the assimilated velocity such vertical the noise is filtered; see the zoomed area on the right.

We consider the behavior of the discretization error versus the discretization step $\Delta$ using noise-free data. This is a consistency check; numerical results, displayed in Figures \ref{F:e_vs_h_perc} (left) are consistent with the FE theory (for the specific choice of FE spaces), i.e. $E_U = \frac{\|\UU-\UU_{P}\|_2}{\|\UU_{P}\|_2} = \mathcal{O}(\Delta^2)$ and $E_P = \frac{\|\PP-\PP_{P}\|_2}{\|\PP_{P}\|_2} = \mathcal{O}(\Delta)$. 

In Figure \ref{F:e_vs_h_noise} we report in logarithmic scale $E_U$ versus $\Delta$ for \gls{SNR} =20 (left) and 8 (right). We observe a behavior proportional to $\Delta^{0.5}$ (see Table \ref{T:life_p_eh} for the complete report of the test case). To understand this result we note that, at this initial stage, the number of data on $\Gamma_{in}$ is mesh dependent, more precisely $N_s \propto \Delta^{-1}$; hence, $E_U\propto {N_s^{-1}}$. This is coherent with statistical theory: the standard deviation of the sample mean of $N$ random variables is proportional to $N^{-0.5}$. In Chapter \ref{C:da_nse} further results confirm this interpretation.

Figure \ref{F:e_vs_h_perc} (right) reports in a logarithmic scale $E_U$ versus \gls{SNR}$^{-1}$; here, the error increases linearly with the amount of noise, which is reasonable in the case of linear problems.   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | c | c | c | c | c |}
    \hline
	\gls{SNR}			&	$n$ = 40	&		$n$ = 60	&	$n$ = 80	&	$n$ = 120	&	$n$ = 160\\
	\hline
	$\infty$	&	1.1874e-02 	&   	5.1610e-03 	&   3.0476e-03 	&   1.3949e-03 	&   8.1229e-04\\
	\hline
	20			& 	3.5414e-02 	&   	2.9724e-02  &  	2.4532e-02 	&  	2.0646e-02  &  	1.7957e-02\\
    \hline
	8			& 	7.4238e-02 	&  		6.7801e-02 	&  	5.8876e-02  & 	4.6559e-02 	&  	4.1652e-02\\
	\hline
\end{tabular}
\caption{Velocity relative errors in correspondence of different values of the discretization step for different choices of SNR.}
\label{T:life_p_eh}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=0.46]{pvd_stokes.png}
\caption{Computed pressure on the background, computed velocity and the noise-free data (black vector field).}
\label{F:stokes1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=0.46]{pvd_stokes5n.png}
\caption{Computed pressure on the background, computed velocity and the noisy data with SNR = 10 (black vector field).}
\label{F:stokes2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.28]{e_vs_h_0noise.pdf}
\end{minipage}
\hspace{0.25cm} 
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.28]{e_vs_perc_randFIX.pdf}
\end{minipage}
\caption{Left: relative errors for velocity and pressure versus $\Delta$ with noise-free data. Right: mean relative errors, over 20 realizations, versus SNR$^{-1}$ = 0.03, 0.05, 0.07, 0.10, 0.12, 0.15, 0.17, $n$ = 60.}
\label{F:e_vs_h_perc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.5\linewidth} % A minipage that covers half the page
\centering
\includegraphics[scale=0.28]{e_vs_h_5noise.pdf}
\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.28]{e_vs_h_12noise.pdf}
\end{minipage}
\caption{Mean relative errors, over 20 realizations, for the velocity field versus $\Delta = 1/8,\,1/12,\,1/16,\,1/24,\,1/32$, with SNR = 20 (left) and SNR =8 (right).}
\label{F:e_vs_h_noise}
\end{figure}