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

%: ----------------------- introduction file header -----------------------
\chapter{Introduction}\label{C:introduction}

% the code below specifies where the figures are stored
\ifpdf
    \graphicspath{{1_introduction/figures/PNG/}{1_introduction/figures/PDF/}{1_introduction/figures/}}
\else
    \graphicspath{{1_introduction/figures/EPS/}{1_introduction/figures/}}
\fi

%: ----------------------- HELP: references
% References can be links to figures, tables, sections, or references.
% For figures, tables, and text you define the target of the link with \label{XYZ}. Then you call cross-link with the command \ref{XYZ}, as above
% Citations are bound in a very similar way with \cite{XYZ}. You store your references in a BibTex file with a programme like BibDesk.

Reliability of numerical methods and simulations is becoming a crucial issue as scientific computing has an increasing role in different fields of engineering, science and, more generally, society. The quality of numerical results needs to be certified not only for the mathematical investigation of the different methods, but also for assessing its impact in practice (see \cite{LibroSIAM, NASAReport}). As a matter of fact, scientific computing is extensively used beyond traditional engineering studies (civil, aerospace, automotive), in geophysics, sport, finance, medicine. A careful analysis  of the impact of errors and uncertainties is therefore not only needed for a theoretical understanding (and ranking) of the methods but also for the reliability of numerical simulations as a tool for supporting, for instance, the decision making of a medical doctor. 

More specifically, in the last 20 years numerical methods for incompressible fluid dynamics have been fostered by applications to the cardiovascular system (see e.g. \cite{cardiov,RRGT,taylor.eal:predictive}). Mathematical and numerical models have been progressively  used as a tool for supporting medical research in the cardiovascular science. In this particular field, whose social relevance is testified by the impact of cardiovascular diseases in the industrialized Countries\footnote{According to the statistics reported by the Center for Disease Control {\tt www.cdc.gov/nchs/data/nvsr/nvsr55/nvsr55\_19.pdf}.}, numerical simulations are increasingly   considered as a valuable tool for investigating physio-pathological dynamics or even designing therapies  \cite{Karn,RRGT,Complex,cardiov,Tay2,Tay1}. {\it In silico} experiments can provide remarkable insights into a physio-pathological process completing more traditional {\it in vitro} and {\it in vivo} investigations. Numerical models have been playing the role of ``individual based" simulators, able to give a dynamical representation of the biology of a specific patient as a support to the prognostic activity. 

At the same time, the need for  quantitative responses for diagnostic purposes has strongly stimulated the design of new methods and instruments for measurements and imaging. On the one hand, we can simulate large 3-dimensional portions of the cardiovascular system of a real patient properly including simplified models
for the peripheral sites (see e.g. \cite{taylor.eal:predictive,Tay2,Karn,cardiov,Complex}). On the other, images and measures nowadays provide doctors and bioengineers with  a huge amount of data. As an example, Figure \ref{F:MRIsection} shows an MRI \cite{brummer} of the ascending aorta where blood velocity measurements have been collected. Obviously these data offer new possible benchmarks for the numerical simulations (see e.g. \cite{heys2010}).  However, beyond the validation, it is possible to merge simulations and measures by means of more sophisticated numerical techniques. A strong integration of data and numerical modeling is expected to bring a great benefit to the development of mathematical and computational tools with a clinical impact. This procedure is called {\it data assimilation} (\gls{DA})  (see e.g. \cite{navon}). With this name we mean the ensemble of methods for merging observed (generally sparse and noisy)  information into a numerical model based on the approximation of physical and constitutive laws. In other words, the goal of \gls{DA} is to merge together heterogeneous (in nature, quality, and density) sources of information in order to
retrieve a consistent state for phenomena of interest. In particular, in hemodynamics, numerical simulations are improved by the merging of data that allow to include effects otherwise difficult to model (at the qualitative or quantitative level), such as the presence of tissues surrounding an artery or the motion of the heart affecting the aortic dynamics. Also, since measures are in general affected by noise, the assimilation based on physical and constitutive laws introduces a sophisticated filter, forcing the consistency with basic principles.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.4\textwidth]{field.png} &
\includegraphics[width=0.3\textwidth]{MRIsection.png}
\end{tabular}
\end{center}
\caption{View of blood measured velocities in MRIs of the ascending aorta.}
\label{F:MRIsection}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Even though the term \gls{DA} has often been used for time dependent phenomena, in this thesis (as done in \cite{heys2010}) we consider also stationary phenomena and we refer to \gls{DA} also as a technique for the inclusion of data which are sparse in space and collected at the same instant of time. In the next sections we present a brief history of \gls{DA} for fluid geophysics applications (oceanography and meteorology) mainly based on the excellent review by Blum et al. in \cite{navon} and we summarize the state of the art of \gls{DA} in hemodynamics applications.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Methods for DA in fluid geophysics} 
The mathematical modeling of geophysical flows has experienced a tremendous development during the last decades, mainly due to the growth of the available computing resources and to the improvement of networks of observations. Therefore, it is clear that modeling and simulating have to take into account observations; these are irregularly distributed in space and time, and they have different structures of random error. Errors arise from various sources such as instrumental noise, environmental noise, sampling, and interpretation of sensor measurements; according to the context, a different confidence will be attributed to the measurements. For instance, all oceanic dynamical models are imperfect, with errors due to the approximate physics (or biology or chemistry) that governs the explicit evolution of the state variables and that parameterizes the interaction of the state variables and the discretization of dynamics into a numerical model \cite{wang_data_2000}. 

Originally, the problem of \gls{DA} was to determine the initial condition of a dynamical system from observations. Since the same mathematical tools may be used, \gls{DA} also includes the estimation of model parameters, boundary conditions or the state itself. The development of \gls{DA} methodology has mainly experienced three stages: objective analysis, statistical interpolation, stochastic methods and variational analysis. Simple analysis methods were mostly used in the 50's, when computers were unavailable or at the beginning stage; these techniques were the earliest bases of \gls{DA}. In the 60's and 70's, statistical considerations were introduced into the atmospheric \gls{DA}. Based on these considerations, some forms of optimal interpolation were used to assimilate observations into forecast models. In the same decades the Kalman filter was introduced and efficiently improved (EKF, EnKF, UKF). In the 80's and 90's, atmospheric \gls{DA} switched to variational methods.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Objective analysis and optimal interpolation} 
Objective analysis, or interpolation, was introduced, in the context of meteorology, by Panofsky \cite{panofsky} in 1949; the 50's featured the advent of computers and the evolution of numerical analysis. This method interpolates the measurements from the points of observation toward the grid points. Specifically, it is based on a technique of 2-dimensional least squares fitting; it consists in expanding the state variables in a series of polynomials about the observation point, minimizing the square of their difference with the observed values. This procedure was perfected by the idea of Bergth\'orsson and D\"{o}\"{o}s \cite{doos} in 1955, who claimed that it is not possible to get a reasonable estimate only by interpolation of observations; thus, they included a mathematical model in the assimilation process. The method is based on the assumption that some observations are available prior to the assimilation instant (12 hours in advance); then, the estimate is a weighted mean of the current observations and the state forecast, i.e. the state field obtained by advancing in time from the previous observations via a mathematical (numerical) model. Even though the method is simple to implement, the retrieved fields will not be necessarily consistent from the physical viewpoint. Despite this, these methods were the most commonly used in operational meteorology until the last decade.

Optimal interpolation, also referred to as statistical interpolation in meteorology, can be seen as a simplification of the Kalman Filter (even though it was introduced before it), described in the next section. The assimilated solution is a linear combination of the deterministic solution with the misfit between data and predicted state; as opposed to Kalman approaches, the matrix weighting the misfits is empirically assigned (see \cite{kalnay} for further details).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Kalman filter and its extensions}
The Kalman filter\footnote{Here, with filter we mean a process that removes from a signal some unwanted components, e.g. the noise, as in our case.} (KF) was developed by Kalman in 1960 \cite{kalman} as a new approach to linear filtering and prediction problems. It is a recursive filter for the estimation of the state of a dynamic system from incomplete and noisy measurements. Because of its limitations in terms of models (designed for linear systems) and memory requirements for the storage of structures involved, several modifications and improvements have been developed, such as the Extended KF (EKF), the Ensemble KF (EnKF) and the Unscented KF (UKF). We report briefly the formulation of the basic KF and we touch on few details of its extensions.

The KF is a scheme designed for linear dynamic systems and it is based on Bayesian estimation; it generates the sequential, unbiased, {\it minimum error variance} estimate based on a linear combination of all past measurements and dynamics. 
%Fundamental stages of the scheme are {\bf a)} the forecast of the state and of its error covariance and {\bf b)} the combination of data and forecast and the error update; this consists in the linear combination of the dynamical forecast with the difference between data and prediction. The weighting matrix is the so called Kalman gain matrix \cite{robinson}. 
Following the formulation reported in \cite{navon}, we summarize the procedure in four simple steps. 

Consider a linear model representing the observation process:
\begin{equation}\label{E:KF_obs}
y_k = H_k u_k + \varepsilon_k; 
\end{equation}
here, $k$ is a multiple of the number of time steps between consecutive observations, $u_k$ is the vector of state variables, $y_k$ is the vector of the measures and $\varepsilon_k$ is an additive random noise representing the error; this is assumed to be white with zero mean and correlation matrix $R_k$. For each time step we consider a discretized stochastic dynamic system
\begin{equation}\label{E:KF_model}
u_k = M_{k-1} u_{k-1} + \eta_{k-1}
\end{equation} 
where $M_k$ describes the model dynamics and $\eta_k$ is the random error associated with model parameters. We assume that $\eta_k$ is distributed as $\varepsilon_k$, with correlation matrix $Q_k$. 
\begin{equation}\label{E:KF_ alg}
\begin{array}{ll}
{\rm {\it Algorithm}}	&		\mathbf{1.}\;\hbox{advance in time:}\left\{
								\begin{array}{l}
								u_k^f = M_{k-1}u_{k-1}^a\\
								P^f_k = M_{k-1}P^a_{k-1}M^{\TT}_{k-1} + Q_{k-1}
								\end{array}\right.\\
								
						&		\mathbf{2.}\;\hbox{compute the Kalman gain:}\; K_k = P_k^f H^{\TT}_k(H_k P^f_k H^{\TT}_k+R_k)^{-1}\\
						
						&		\mathbf{3.}\;\hbox{state update:}\; u_k^a = u_k^f + K_k(y_k-H_ku_k^f)\\
						
						&		\mathbf{4.}\;\hbox{error correlation matrix update:}\; P_k^a = (I-K_k H_k) P^f_k\\
\end{array}
\end{equation}
Here, $P_k^f = E\{(u_k-u_k^f)(u_k-u_k^f)^{\TT}\}$ and $P_k^a = E\{(u_k-u_k^a)(u_k-u_k^a)^{\TT}\}$; the superscripts $f$ and $a$ stand for {\it forecast} and {\it assimilated} respectively and $E$ refers to the mean value for random variables.

This filter is limited to linear models, but in general many problems in fluid geophysics, including parameter identification, are nonlinear. This limitation explains the introduction of the EKF for nonlinear processes. In this case, equations \eqref{E:KF_obs} and \eqref{E:KF_model} are to be substituted with 
\begin{equation}\label{E:EKF_obs}
\begin{array}{l}
y_k = h(u_k,\varepsilon_k);\\
u_k = f(u_{k-1},\eta_{k}).
\end{array}
\end{equation} 
The KF, designed for linear problems, cannot be directly applied to equations \eqref{E:EKF_obs}. The idea of the EKF is to consider, at each iteration, the linearization of \eqref{E:EKF_obs} around the previous state. In this way, matrices $H_k$ and $M_k$ in algorithm \eqref{E:KF_ alg} are substituted with the Jacobian of $h$ and $f$ respectively, evaluated at $u_{k-1}$. This estimator is not optimal, might underestimate the covariance matrices and has high memory requirements. 

The extension to the EnKF avoids the latter using probability distributions represented by samples, i.e. covariance matrices are replaced by sample covariance matrices.

When the state transition and observation models ($f$ and $h$) are highly nonlinear the EKF might give poor performance \cite{UKF1}. This is because the covariance is propagated through linearization of the nonlinear model. The UKF represents a derivative-free alternative to the EKF and provides better performance at an equivalent computational complexity. In fact, in the UKF the state distribution is still Gaussian but it is specified using a minimal set of ``carefully chosen'' sample points or snapshots which completely capture the state mean and covariance. Furthermore, these points, when propagated through the nonlinear system, capture the mean and covariance at the subsequent step to the second order. Such sample points are chosen via {\it Unscented transform} of the state at the current time step \cite{UKF2}. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Nudging method}\label{SS:introduction_geophysics_nudging}
The nudging, or {\it dynamic relaxation}, method was introduced in 1977 by Davies and Turner in \cite{nudging} and it can be classified as a combination of Kalman filtering and a variational approach. 
%The nudging method can be classified as a variational method, continuos in time, that uses dynamic relaxation to adjust the state variable towards the observation. The assimilation is attained introducing in the governing equations a non-physical diffusive-type forcing term, with a tunable relaxation coefficient. The equations are solved in a interval of time known as the {\it pre-forecast} so that the final state is as close as possible to the measures and can be used as the initial state, automatically in dynamic balance, to produce regular forecasts \cite{haltiner}. 

%The size of this coefficient is chosen empirically so to keep the nudging terms small in comparison to the dominating forcing terms (in order to avoid the rebounding effect that slows down the assimilation process), yet large enough to impact the simulation. This method is flexible and computationally cheaper than control approaches. 

According to \cite{haltiner}, we can summarize the procedure in two simple steps. Given some data at time $t=t_0$, {\bf 1)} specify an arbitrary but consistent initial condition at time  $t=t_0-T$ (where $T$ is the so called {\it pre-forecast} period); {\bf 2)} solve the governing equations during the pre-forecast, $[t_0-T, t_0]$, including a forcing term proportional to the difference between the predicted state at time $t_0$ and the observation so to \textit{dynamically relax} the solution towards the observations. 
If we assume to have a mathematical model, discretized in space (using e.g. finite differences, volumes or elements), the time dependent problem reads:
\begin{equation}
\begin{array}{ll}
\dfrac{d\UU}{dt}=F(\UU)		&	t\in(t_0-T,t_0];\\
\UU(t_0 - T)=\UU_0.
\end{array}
\end{equation}
Here, $\UU$ is the discretized state variable, and $\UU_0$ the initial condition. If $\widetilde{\UU}$ is a given measure and $\overline{\UU}$ is a projection of the observation onto the discretization space, the assimilation consists in the solution of the following initial value problem.
\begin{equation}
\begin{array}{ll}
\dfrac{d\UU}{dt}=F(\UU) + c\;(\overline{\UU}-\UU)		&	t\in(t_0-T,t_0];\\
\UU(t_0-T)=\UU_0.
\end{array}
\end{equation}
The relaxation coefficient $c$ has to be chosen so to achieve an optimal compromise between model and observations; its choice might be based on statistical considerations. If $c$ is too large, the non-physical forcing term will dominate and the predicted state would completely fit the noisy measures, amplifying the observational errors and compromising the consistency of the solution. When $c$ is too small, the observation will not significantly affect the estimate. In 1992, Zou, Navon and Le Dimet developed a parameter estimation approach in order to obtain an optimal nudging coefficient \cite{zou}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Variational approaches}\label{ss:introduction_methods_control}
Control approaches consider the minimization of the misfit between the predicted state and the observations under the constraint of the partial differential equations (\gls{PDE}s) governing the variables. Such optimization (or control) problem may be solved with a large variety of numerical methods for \gls{PDE} constrained optimization problems. In the context of fluid dynamics there is a huge literature for the solution of inverse problems, among many others we mention \cite{gunzburger_perspectives_2002, gunzburger1999sensitivities,gunzburger1991analysis-1,gunzburger1991analysis, hinze1,hinze0, ulbrich}.

In general, the main components in the formulation of a variational problem are a state variable $u$, a deterministic model, a control variable $h$, a set of observations and a cost function $\mathcal{J}(u,h;t)$. The model describing the evolution of the variables involved in the process is, in general, a system of nonlinear differential equations
\begin{equation}
\left\{\begin{array}{l}
\mathcal{A}(u,h;t)=0,\\
u(t_0) = u_0.
\end{array}\right.
\end{equation}
Here, $\mathcal{A}$ is a differential operator and $u_0$ is the initial condition. Most of the time the control variable $h$ is the initial condition (e.g. in meteorology), but it might also consist in model parameters or boundary conditions. In general, we assume that for a fixed value of $h$ the model is well-posed. Additional constraints on the control variable are mandatory in relation to the physical quantities described. The measures, say $\widetilde{\UU}$, are considered discrete in space and time and they do not belong to the same space of the state variable. The cost function $\mathcal{J}$ measures the discrepancy between the state variable and the measurements
\begin{equation}
\mathcal{J}(u,h;t)= \int\limits_0^T dist(u;\,\widetilde{\UU})\; dt + \mathcal{R}(h;t).
\end{equation}
The function $dist$ will be specified later and it might include information on stochastic properties of the field. The additional term $\mathcal{R}$ is usually introduced with regularization purposes to prevent potential ill-posedness of the formulation and/or to impose smoothness of the control variable; in fact, when the objective functions do not depend explicitly on the control variable, the solution of the problem could lead to unbounded control variables \cite{gunzburger_perspectives_2002}. 
Optimization methods can be grouped in two main categories: the {\it optimize then discretize} (OD) and the {\it discretize then optimize} (DO) approach. 
%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The OD approach} We introduce the functional 
\begin{equation}
\mathcal{J}(u,h;t)= \int\limits_0^T \|u(h)-\tilde{u}\|\; dt + \mathcal{R}(h;t).
\end{equation}
where $\tilde{u}$ is a projection of $\widetilde{\UU}$ into the state functional space $\mathcal{H}$. In order to solve the optimization problem we introduce the \textit{Lagrangian} function $\mathcal{L}$
\begin{equation}\label{E:lagr}
\mathcal{L}(u,h,\lambda)=\mathcal{J}(u,h;t)-\lambda^*\mathcal{A}(u,h;t).
\end{equation}
Here $^*$ represents the inner product in $\mathcal{H}$ and $\lambda$ is the Lagrange multiplier. The optimality is reached when we find $u$, $h$ and $\lambda$ such that $\mathcal{L}$ is stationary, i.e. when they solve the system of first order necessary conditions \cite{gunzburger_perspectives_2002}:
\begin{equation}\label{E:nec_cond}
\left\{\begin{array}{ll}
\displaystyle\frac{\partial \mathcal{L}}{\partial \lambda}= \mathcal{A}(u,h;t)=0 & \hbox{\textit{state equation}},\\
[3mm]
\displaystyle\frac{\partial \mathcal{L}}{\partial u} = \left(\displaystyle\frac{\partial \mathcal{A}}{\partial u}\right)^*\lambda - \left(\displaystyle\frac{\partial \mathcal{J}}{\partial u}\right)^*=0 & \hbox{\textit{adjoint equation}}, \\
[3mm]
\displaystyle\frac{\partial \mathcal{L}}{\partial h}=\left(\displaystyle\frac{\partial \mathcal{A}}{\partial h}\right)^*\lambda=\left(\displaystyle\frac{\partial \mathcal{J}}{\partial h}\right)^*=0 & \hbox{\textit{optimality condition}} ;
\end{array}\right.
\end{equation} 
where we make the assumption that each argument is independent from the other. This last system is typically huge and difficult to be numerically solved monolithically (the so called ``one-shot'' approach); iterative procedures are usually preferred \cite{gunzburger_perspectives_2002,haber_optimization_2000}. 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.
%%%%%%%%%%%%%%%%%%%%%
\paragraph{The DO approach}  In this case we first discretize the functional and the state equations, getting
\begin{equation}
\left\{\begin{array}{l}
J(\UU,\,\HH)=\|Q(\UU)-\widetilde{\UU}\| + R(\HH);\\
[2mm]
A(\UU,\, \HH) = 0.
\end{array}\right.
\end{equation}
Here, $\UU$ and $\HH$ are the discretized state and control variables, $Q$ is the operator that maps the discretized state into the data space, $R$ is a discrete regularizing operator and $A$ is the discretized differential operator. To solve this algebraic optimization problem we introduce the discrete Lagrangian function
\begin{equation}
L(\UU,\,\HH,\,\boldsymbol{\Lambda}) = J(\UU,\,\HH) - \boldsymbol{\Lambda}^{\TT} A(\UU,\,\HH),
\end{equation}
where $\boldsymbol{\Lambda}$ is the Lagrange multiplier. Optimality is reached when we find $\UU$, $\HH$ and $\boldsymbol{\Lambda}$ such that 
\begin{equation}\label{E:nec_condDO}
\left\{\begin{array}{ll}
\displaystyle\frac{\partial L}{\partial \boldsymbol{\Lambda}}= A(\UU, \, \HH) =0 & \hbox{\textit{state equation}},\\
[3mm]
\displaystyle\frac{\partial L}{\partial \UU} = \left(\displaystyle\frac{\partial A}{\partial \UU}\right)^{\TT}\boldsymbol{\Lambda} - \left(\displaystyle\frac{\partial J}{\partial \UU}\right)^{\TT}=0 & \hbox{\textit{adjoint equation}}, \\
[3mm]
\displaystyle\frac{\partial L}{\partial \HH}=\left(\displaystyle\frac{\partial A}{\partial \HH}\right)^{\TT}\boldsymbol{\Lambda}-\left(\displaystyle\frac{\partial J}{\partial \HH}\right)^{\TT}=0 & \hbox{\textit{optimality condition}}.
\end{array}\right.
\end{equation} 
As for system \eqref{E:nec_cond}, there are plenty of methods for the solution of \eqref{E:nec_condDO}, see e.g. \cite{gunzburger_perspectives_2002, haber_optimization_2000}.
Note that, the DO approach allows to obtain the optimality system in a more straightforward way than working with unknowns belonging to infinite-dimensional functional spaces (for an extensive discussion see \cite{gunzburger_perspectives_2002}). The following example describes how to apply both schemes to a simple \gls{PDE}. 
%%%%%%%%%%%%%%%%%%%%%%%% 
\paragraph{Application to a simple differential problem} 
We consider the following steady state equation:
\begin{equation}\label{E:contr_ex}
\left\{\begin{array}{l l}
-\mu u''=hf & \hbox{in } \dom\\
u=0 & \hbox{on } \partial\dom;
\end{array}\right.
\end{equation}
where $\mu$ and $f$ are a given functions defined on $\overline{\dom}$, and the following objective function:
\begin{equation}\label{E:contr_ex_obj}
\mathcal{J}(u,h)=\displaystyle\frac{1}{2}\int\limits_{\dom}(u-\tilde{u})\hbox{d}\xx + \displaystyle\frac{\sigma}{2}h^2.
\end{equation}
Here $\displaystyle\frac{\sigma}{2}h^2$ is a penalization term introduced in order to limit the growth of the control variable and $\tilde{u}$ represents observed data assumed to belong to the same space of the state variable $u$. The OD formulation reads 
\begin{center} 
\textit{find $u$ such that $\mathcal{J}$ is minimized subject to \eqref{E:contr_ex}}.
\end{center} 
The Lagrangian of the problem reads:
\begin{equation}\label{E:contr_ex_lagrOD}
\mathcal{L}(u,h,\lambda)=\mathcal{J}(u,h)+\int\limits_{\dom}\lambda\left(\mu u''-h\,f\right)\hbox{d}\xx.
\end{equation} 
In order to obtain the optimality solution, the system of necessary conditions \eqref{E:nec_cond} is derived and then discretized (by means of any discretization technique).

The first step of the DO method is the discretization of both the objective function, say $J(\UU,\HH)$, and of the constraints, say $\textrm{A}\UU=\hbox{F}\HH$. Here, $\UU$ and $\HH$ are the discretized state and control variable, ${\rm A}$ and ${\rm F}$ are the discretized differential and source operators. The second step is the solution of the problem 
\begin{center} 
\textit{find $\UU$ such that} $J(\UU,\HH)$ \textit{is minimized subject to} $\textrm{A}\UU=\hbox{F}\HH$;
\end{center}
This is an algebraic optimization problem which can be solved by finding the stationary points of the Lagrangian function, which reads:
\begin{equation}
L(\UU,\HH,\boldsymbol{\Lambda})=J(\UU,\HH)-\boldsymbol{\Lambda}^{\TT}\left(\hbox{A}\UU-\hbox{F}\HH\right).
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Recent advances in DA for the cardiovascular system} 
The tremendous increase of data gathering and the massive change in computational capabilities in the last few years make the \gls{DA} problem a very challenging task and stimulate the design of new, efficient and competitive computational tools for the prediction of the state and the evolution of complex systems, such as the cardiovascular one. The enormous amount of information available to us requires the combination of advanced statistical techniques with precise and sophisticated mathematical models, perfected and largely validated in recent years \cite{Tay2,Tay1,cardiov,RRGT,taylor.eal:predictive}. These new approaches are developments of well-established \gls{DA} techniques and, furthermore, they feature an overlap between estimation theory, control theory and stochastic approaches. Moreover, they benefit from the continuous improvement of numerical methods for order reduction \cite{ROUKF, bertoglio} and advanced discretization techniques \cite{dwight2010,heys2010}. Among the others, we mention and briefly describe the most relevant works.

%\paragraph{Bayesian inference} In 2009 Cotter et al. \cite{cotter} present a Bayesian approach to inverse problems applied to fluid dynamics. They develop a \gls{DA} framework at a continuos level showing the well-posedness of the problem and the existence of a unique point estimator, the {\it maximum a posteriori estimator}, \gls{MAP}. They also show how that abstract theory may be applied to concrete fluid dynamics problems. In Chapter \ref{C:bayes} we present a detailed description of the Bayesian approach to inverse problems applied to hemodynamics. A very good introduction to statistical inverse problems can be found in the book by Kaipio and Somersalo, chpt 4 \cite{somersalo}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Least-squares finite elements} In 2010 Heys et al. propose a weighted least-squares finite element (WLSFE) method  for the inclusion of velocity data in fluid dynamics equations, mainly inspired by hemodynamics applications \cite{heys2010}. They consider particle imaging velocimetry that allows to collect noisy velocity data on 2-dimensional internal layers in a 3-dimensional domain. On the basis of the LSFE method \cite{bochev1,bochev2, gunzburgerLSFEM}, which features flexibility in the enforcement of various boundary conditions, they developed a weighted formulation based on the level of accuracy with which the boundary values are known. Given a system of discretized \gls{PDE}s ${\rm A}\UU = \FF$ the WLSFE solution is the minimum of the functional
\begin{equation}
J = \|{\rm A}\UU-\FF\|_{\dom} + \|\UU-\HH\|_{\partial\dom} + \|\UU-\widetilde{\UU}_1\|_{\Gamma_1} + \|\UU-\widetilde{\UU}_2\|_{\Gamma_2} + ...
\end{equation}
where ${\rm A}$ is a discrete differential operator (corresponding generally to the \gls{NSE}), $\FF$ a forcing term, $\HH$ the boundary data and $\widetilde{\UU}_i$ available noisy data on the layer $\Gamma_i$. Using additional terms which penalize the difference between the numerical solution and the experimental data, the assimilation is weakly enforced. Thus \gls{DA} costs little more than a standard LSFEM calculation. The method can be seen as nudging (for the modification of the model residual), but also as a control approach (for the minimization of a functional). Note that, LSFEM, originally designed for finding the solution of \gls{PDE}s, has been applied to inverse problems since the 90's \cite{bedivan} for the numerical solution of \gls{PDE} constrained control problems; main developers have been Bochev and Gunzburger \cite{bochev1,bochev2, gunzburgerLSFEM}. 
In the application to the \gls{NSE}, for recasting the system into a linearized first order differential system, the authors consider a non-primitive variable set. This approach, as opposed to primitive variables formulations, might be less conducive to the straightforward inclusion of available ``boundary'' conditions. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Bayesian interpretation of the LSFEM} In the same year Dwight \cite{dwight2010} presents a reinterpretation of the WLSFEM in terms of assimilation of data; in fact, Heys and collaborators do not present their work in an inverse problem framework. Dwight shows that the WLSFE solution can be interpreted as the {\it maximum a posteriori} (\gls{MAP}) estimator in a variational Bayesian approach to \gls{DA} (an extensive description of the Bayesian method is given in Chapter \ref{C:da_bayesian}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Reduced-order unscented Kalman filter} In 2010, based on the UKF, presented in the previous section, Moireau and Chapelle \cite{ROUKF} propose a reduced-order (RO) UKF applied to parameter identification in the context of cardiac biomechanics. Their filtering strategy, usable for any choice of sampling points distribution, provides a tractable filtering algorithm that can be used with large-dimensional systems when the uncertainty space is of reduced size. Such algorithm invokes the original dynamical and observation operators, i.e. it does not require the computation of the tangent operator. Specifically, the covariance matrices are factorized in a form such that the costly computations are performed on a reduced-order matrix of the size of the uncertainty space.\\

In recent years many works investigated efficient and robust coupling algorithms for fluid-structure interaction (FSI) in large vessels. To impact the clinical practice such techniques have to be patient-specific; for this reason, model parameters have to be adapted to available clinical measurements. Two methods for FSI parameter estimation have been proposed in recent publications using Kalman filtering \cite{bertoglio} and control theory \cite{MPAVCV}. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{ROUKF applied to fluid-structure interaction problems} In a very recent work by Bertoglio et al. \cite{bertoglio} the ROUKF is used for the estimation of constitutive parameters in FSI systems. In particular, from wall displacement measurements extracted from medical images, they predict the artery wall stiffness. Using this sequential algorithm, the model prediction is improved at each time step analyzing the difference of the measures and model output. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Control theory approach to the estimation of FSI parameters} The same problem has been addressed by Perego et al. in \cite{MPAVCV}. Here, the estimation of the Young modulus is pursued using an OD technique in control theory. The problem of assimilating vessel displacements from registered medical images is formulated as the minimization of the difference between observation and computed variables tuning the variable representing the stiffness itself.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Thesis outline}
In this work we propose a \gls{DA} approach for including noisy and sparse measurements of the velocity into the simulation of the \gls{NSE}, driven by hemodynamics applications. This inverse, potentially ill-posed, problem is addressed using two approaches oriented to finding point and spread estimates for the velocity field. The purpose is {\bf 1)} recovering an accurate and noise filtered approximation of the blood flow and {\bf 2)} quantifying the reliability of numerical results. The first goal is achieved by solving an optimization problem where we use a DO technique to minimize the misfit between the recovered velocity field and the data, subject to the incompressible \gls{NSE}. The control variable is the normal stress at the inflow boundary of the domain of interest, the vessel inlet section; the reason for this choice resides in the biomedical application and will be specified in Chapter \ref{C:da_nse}. The \gls{DA} procedure is a combination of the Newton method for handling the nonlinear term of the \gls{NSE} and the \gls{DA} procedure designed and tested for the linearized problem. Then, the uncertainty quantification is addressed using a Bayesian approach to inverse problems recast in a variational formulation. The prediction of probability density functions for velocity and normal stress allows us to estimate how likely are those variables to take values in specific intervals. This approach, that accounts for the nature of the measurement noise, is tested on 2-dimensional and axisymmetric 3-dimensional geometries and compared with deterministic approaches.

In Chapter \ref{C:preliminary}, based on \cite{md_av1}, we present preliminary results aimed at finding the most accurate and efficient approach to the solution of the \gls{DA} problem. Different methods presented in previous sections are considered and compared in terms of accuracy with respect to reference solutions and computational time. The investigation is conducted on the linearization of the \gls{NSE}, namely the Stokes problem. This preliminary analysis results in the choice of the DO method.

In Chapter \ref{C:da_nse}, based on \cite{md_av3} and \cite{md_av2}, we present the core of our methodology. We design a \gls{DA} procedure where a DO approach is used to minimize the difference between measures and assimilated solutions for the generalized Stokes problem. The combination of such procedure with the Newton method for the solution of the forward \gls{NSE} yields an algorithm of assimilation for velocity data. We discuss conditions on the location of measurements that guarantee the existence and uniqueness of an optimal solution for the linearized problem and we give numerical proof of convergence of the algorithm for the nonlinear case. Numerical results, with both noise-free and noisy data, certify the theoretical analysis. Also, we investigate the sensitivity to the noise of non-primitive variables of medical interest.

In Chapter \ref{C:da_bayesian} we introduce statistical inverse problems and the Bayesian variational formulation for the discretized variables. This approach yields a probability density function (\gls{PDF}) for velocity and normal stress which can be used to estimate confidence intervals for such variables. Moreover, statistical parameters of the \gls{PDF}s are used as point estimators for the control variable (the normal stress on the inflow boundary), specifically the \gls{MAP} and the {\it maximum likelihood} (\gls{ML}) estimators. We also provide details regarding a regularization technique based on the optimality result of Chapter \ref{C:da_nse} and we discuss how it affects statistical properties of the measurement noise. Then, we show numerical results in 2-dimensional and 3-dimensional geometries for data affected by Gaussian noise. 

In Chapter \ref{C:da_unsteady} we propose two formulations for the unsteady (linearized and nonlinear) problem. Due to the hemodynamics application we consider time periodic phenomena and we investigate possible techniques for addressing this issue. We show numerical results on 2-dimensional and 3-dimensional geometries where we test the accuracy of the assimilation with respect to reference solutions.

Concluding remarks, future work guidelines and insights are reported in Chapter \ref{C:conclusion}. 

%\figuremacroW{largepotato}{A common glucose polymers}{The figure shows starch granules in potato cells, taken from \href{http://molecularexpressions.com/micro/gallery/burgersnfries/burgersnfries4.html}{Molecular Expressions}.}{.8}

%: ----------------------- HELP: adding figures with macros
% \figuremacro{1}{2}{3}{4} calls up a series of commands formating your image.
% 1 = name of the file without extension; PNG, JPEG is ok; GIF doesn't work
% 2 = title of the figure AND the name of the label for cross-linking
% 3 = caption text for the figure
% 4 = size relative to text width which is 1; use this to reduce figures




