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

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

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

%: ----------------------- name of chapter  -------------------------
\chapter{A deterministic approach to DA}\label{C:da_nse} % top level followed by section, subsection
Preliminary results on the Stokes problem reported in Chapter \ref{C:preliminary} provide the ground for the design of a robust, noise filtering method for the estimation of blood flow from noisy velocity measures and numerical simulations. The discretize then optimize (\gls{DO}) technique is the basis for the development of a flexible and reliable assimilation tool.

Objectives are challenging: we want to construct an accurate estimator which works in case of highly noisy data and complex geometries regardless (as much as possible) of the dynamics, turbulent or laminar, of the flow. 
We want to provide conditions that guarantee good mathematical properties to the estimator, including existence and uniqueness and its stability with respect to data perturbation.

In this chapter, mainly based on \cite{md_av3} and \cite{md_av2}, we present the core of the assimilation procedure. We show how our estimator is able to accomplish such goals, fundamental in the design of a prediction tool for the medical practice. We introduce a robust and noise filtering \gls{DO} method for the linearized Navier-Stokes equations (\gls{NSE}) and we provide a possible alternative to common Tikhonov-like regularization methods. The latter are in general mandatory in the treatment of ill-posed inverse problems, and, however, might be computationally demanding. Our regularization approach is based on an optimality result for the control problem; we find conditions on the number and location of the data which ensure the existence of a unique estimator. The extension of the method to the \gls{NSE} is an iterative algorithm that still features the desired properties, also in presence of advection dominated phenomena, i.e. high Reynolds number. The application to non-trivial domains approximating cardiovascular geometries gives promising results in both the prediction of the flow and of flow related variables of high medical interest, namely the wall shear stress.

In Section \ref{S:da_nse_mathematical} we present the continuous formulation of the \gls{DA} problem for the \gls{NSE}; in Section \ref{S:da_nse_DO} we describe in detail the discrete formulation applied to the linearization of the \gls{NSE}, we prove the existence and uniqueness of the assimilated solution under a condition on the location of the measurements. Moreover, we discuss the sharpness of this condition and introduce the associated regularization technique. Numerical results for noise-free data are consistent with the theory. With noisy data results point out the noise-filtering action of the assimilation procedure. In Section \ref{S:da_nse_nonlinear} we generalize the technique to the steady \gls{NSE}; numerical results show the robustness of the proposed iterative algorithm. In Section \ref{S:da_nse_real_geometries} we use the \gls{DA} estimator on more complex test cases; simulations on 2-dimensional non-trivial domains and axisymmetric geometries show the flexibility, accuracy and reliability of the method in the prediction of velocity and velocity related variables.

%: ----------------------- contents from here ------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mathematical formulation of the DA problem}\label{S:da_nse_mathematical}
Following the results obtained in Chapter \ref{C:preliminary} we formulate the assimilation technique as a control problem where the misfit between the predicted velocity and the data is minimized.

Let us denote by $\dom$ a domain in $\mathbb{R}^{d}$ ($d=2,\,3$; in real applications $d=3$). We assume (see Figure \ref{F:MRIsection} and \ref{F:domain}) that the domain of interest $\dom$ features an inflow boundary $\Gamma_{in}$, an outflow boundary $\Gamma_{out}$ and the physical wall of the vessel $\Gamma_{wall}$. $\Gamma_{in}$ and $\Gamma_{out}$ can possibly consist of several sections. Variables of interest are velocity $\uu(\xx)\in{\bf H}_{\Gamma_{wall}}^1(\dom)$ and pressure $p(\xx)\in L^2(\dom)$ which are assumed to obey the \gls{NSE} in $\dom$. Also, we assume to have $N_s$ velocity measures $\dd$ available at some sites\footnote{Notice that we use the word ``sites'' for the location of measurements, as opposed to the word ``nodes'' for points where velocities are computed. In general sites and nodes are different, but we do not exclude that the intersection of sites set and nodes set in non-empty.} $\xx_i^{m}\in \dom$, for $i=1,\,...,\,dN_s$. We formulate the \gls{DA} problem as an inverse problem of the form
\begin{equation}\label{E:min_cont}
\begin{array}{l}
\min\limits_{\hv} \mathcal{J}(\uu,\,\hv) = dist(f(\uu), \, \dd) + \mathcal{R}(\hv) \\
[3mm]
\hbox{s.t.} \;\; 
\left\{
\begin{array}{ll}
-\nu\;\gr\cdot(\gr\uu + \gr \uu^{\rm T}) + (\uu\cdot\gr)\uu + \gr p=\mathbf{s}       & \textrm{in}\;\; \dom,\\
[3mm]
\gr\cdot\uu = 0                             								         & \textrm{in}\;\; \dom,\\
[3mm]
\uu = \mathbf{0} 																	 & \hbox{on}\;\;\Gamma_{wall},\\
[2mm]
\displaystyle-\nu \,(\gr\uu+\gr\uu^\TT)\mathbf{n}+ p\,\mathbf{n}=\mathbf{h} & \textrm{on}\;\;\Gamma_{in},\\
[3mm]
\displaystyle - \nu \,(\gr\uu+\gr\uu^\TT)\mathbf{n}+p\,\mathbf{n}=\mathbf{g}& \textrm{on}\;\;\Gamma_{out}.
\end{array}\right. 
\end{array}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=.35]{vessels_with_data.pdf}
\caption{Possible 2-dimensional and 3-dimensional domains with associated measured sparse data.}
\label{F:domain}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
Here, $dist(\cdot,\,\cdot)$ is a distance which will be specified later on; $f$ is an observation operator for the velocity field; ${\bf s}\in\mathbf{H}^{-1}(\dom)$ is the source term; ${\bf g}\in\mathbf{H}^{-1/2}(\Gamma_{out})$ is the Neumann boundary condition on the outflow section; $\mathcal{R}$ is a regularization term added to improve properties of the functional to be minimized (as we will see later). 
A Newtonian rheology is supposed to hold, since it is a common assumption in large and medium blood vessels \cite{cardiov}; $\nu$ is the kinematic viscosity. The choice of homogeneous Dirichlet boundary conditions on $\Gamma_{wall}$ is due to the fact that we consider fixed geometries. We use the inflow normal stress data ${\bf h}\in\mathbf{H}^{-1/2}(\Gamma_{out})$ as control variable; this choice is arbitrary but can be justified by the fact that an accurate knowledge of inflow-outflow boundary conditions is quite problematic in real applications, since these are ``artificial'' boundaries (see \cite{Heywood,veneziani_phd}). 
We present the method for the linearized equations at first; here, the nonlinear convection term $(\uu\cdot\gr)\uu$ is replaced with $(\boldsymbol{\beta}\cdot \gr)\uu$, where $\boldsymbol{\beta}$ is a known advection field. Then, we proceed by treating the nonlinear case in Section \ref{S:da_nse_nonlinear}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discretize then optimize method for the linear problem}\label{S:da_nse_DO}
For the numerical solution of the optimization problem we focus on the linearized problem introduced in the previous section. The choice of the optimization scheme relies upon the comparison analysis presented in Chapter \ref{C:preliminary} and in \cite{md_av1} where the \gls{DO} technique proved to be the most efficient method.  In this procedure we first discretize the continuous formulation by means of the finite element (\gls{FE}) method; this yields the following algebraic optimization problem (for $\alpha\geq 0$)
\begin{equation}\label{E:min_discr}
\begin{array}{l}
\min\limits_{\HH} J(\VV,\,\HH) = \displaystyle\frac{1}{2}\|{\rm D}\VV - \dd\|^2_2 + \displaystyle\frac{\alpha}{2}\|{\rm L} \HH\|^2_2\\
[3mm]
\hbox{s.t.} \;\; {\rm S}\VV = \RR_{in}^{\TT}{\rm M}_{in} \HH + \FF.
\end{array}
\end{equation}
For $\alpha> 0$, the term $\frac{\alpha}{2}\|{\rm L} \HH\|^2_2$ has a regularization purpose (Tikhonov regularization, see \cite{Hansen}) and it is intended to improve the convexity of the functional to be minimized. More in general, it improves the conditioning of the problem, that can be affected by the presence of noise. L is selected such that L$^{\TT}$L is positive definite. Different techniques can be used for an optimal choice of the regularization parameter $\alpha$; among the others, we mention the discrepancy principle (\gls{DP}) and the generalized cross validation (GCV) \cite{Hansen}. $\VV=[\UU\;\PP]\in\mathbb{R}^{N_u+N_p}$ is the vector of discretized velocity $\UU\in\mathbb{R}^{N_u}$ and pressure $\PP\in\mathbb{R}^{N_p}$ belonging to inf-sup compatible \gls{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_{in}$ is the number of degrees of freedom (\gls{DOF}s) of the velocity on $\Gamma_{in}$. We introduce the matrix $\rm S$, a modification of the matrix introduced in Chapter \ref{C:preliminary}.
\begin{equation}\label{E:notations}
{\rm S}=\left[\begin{array}{cc}
{\rm C}+{\rm A} & {\rm B}^\TT\\
{\rm B} & {\rm O} 
\end{array}\right].  
\end{equation}
Here, we let $\rm C$, $\rm A$ $\in \mathbb{R}^{N_u,N_u}$ and $\rm B$ $\in \mathbb{R}^{N_p,N_u}$ be the discretization of the diffusion, advection 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}_j^{\TT}) \cdot\gr\boldsymbol{\varphi}_i \, {\rm d}\xx,  &  \quad i,\,j = 1,\,...,\,N_u, \\
\left[{\rm A}\right]_{ij} & = & \int\limits_{\dom} \boldsymbol{\beta} \cdot \gr\boldsymbol{\varphi}_j \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 $\boldsymbol{\varphi}$ and $\psi$ are the Lagrangian \gls{FE} basis functions associated with velocity and pressure respectively.

Let us recall some notations introduced in the previous chapter. We let ${\rm Q}\in\mathbb{R}^{dN_s,N_u}$ be a matrix such that $[{\rm Q} \UU]_i$ is the numerical solution evaluated at the data sites and corresponding to $[\dd]_i$, for $i=1,\,...,\,d N_s$. When measurement sites coincide with \gls{DOF}s $[{\rm Q} \UU]_i$ is simply the corresponding entry of $\UU$; when measurement sites do not correspond to \gls{DOF}s points, $[{\rm Q} \UU]_i$ is the value of the \gls{FE} numerical solution evaluated at $\xx^m_i$ obtained by a weighted sum of the \gls{FE} Lagrangian basis functions.
As an example, let us consider the grid of Figure \ref{F:matrixQ}, with a P1bubble-P1 discretization for the velocity and pressure. Two possible sites location on $\Gamma_{in}$ are illustrated. The row of the matrix ${\rm Q}$ (relative to the first velocity component, denoted by ${\rm Q}_x$) in correspondence of the site $\xx^*$, is also reported for each case. While in the first case at row $*$ we have one non-zero entry, in the second case, we have two non-zero entries corresponding the grid nodes of the edge $\xx^*$ belongs to. 
%{\it In fact, if we consider the marked site ($\xx^*$) of Figure \ref{F:matrixQ} (right), we have 
%\begin{equation}\label{E:u_weights}
%\uu(\xx^*) = \UU_{a} \varphi_{a}(x^*,\,y^*) + \UU_{b} \varphi_{b}(x^*,\,y^*),
%\end{equation}
%where $\varphi_{a\, (b)}(x^*,\,y^*)$ is the lagrangian \gls{FE} basis function for the velocity relative to the node $\xx_a$ ($\xx_b$). Consequently, in the matrix Q at the row $*$ we have a non-zero entry at columns $a$ and $b$. Similar arguments can be used for Q$_y$ (and Q$_z$ for 3D problems). }{\tt forse questa parte la leverei.}
We extend the selection matrix to pressure \gls{DOF}s, introducing the matrix ${\rm D}=[{\rm Q}\;\;\OO]$.

We let R$_{in,u}\in\mathbb{R}^{N_{in},N_u}$ be 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 \RR_{in} \in \mathbb{R}^{N_{in},N_u+N_p}.
\end{equation}

We 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}$.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t] 
\centering
\includegraphics[scale=0.25]{sites_and_Q.png}
\caption{Possible choices of sites location and corresponding rows for the matrix Q (only Q$_x$ is reported) in correspondence of site $\xx^*$ . Left: sites coincident with grid nodes. Right: sites non-coincident with grid nodes.}
\label{F:matrixQ}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
For the solution of problem \eqref{E:min_discr}, we use the reduced Hessian method to solve the system necessary conditions for optimality induced by the Lagrangian
\begin{equation}\label{E:lagr}
L(\VV,\HH,\boldsymbol{\Lambda})=\displaystyle\frac{1}{2}\|\textrm{D}\VV-\dd\|_2^2+\frac{\alpha}{2}\;\|\textrm{L}\HH\|_2^2+\boldsymbol{\Lambda}^\TT(\textrm{S}\VV-\RR_{in}^{\TT}{\rm M}_{in}\HH - \FF),
\end{equation}
where $\boldsymbol{\Lambda}\in \mathbb{R}^{N_u+N_p}$ is the discrete Lagrange multiplier. The set of necessary conditions for optimality is given by the KKT system (see e.g. \cite{nocedal_numerical_2000}).
\begin{equation}\label{E:lagr_necessary}
\left\{\begin{array}{l}
\dfrac{\partial L}{\partial \VV} = \textrm{D}^\TT(\textrm{D}\VV-\dd)+\textrm{S}^\TT\boldsymbol{\Lambda}=\mathbf{0}\\[3mm]
\dfrac{\partial L}{\partial \HH} = \alpha\textrm{L}^\TT\textrm{L}\HH-{\rm M}_{in}^{\TT}\RR_{in}\boldsymbol{\Lambda}=\mathbf{0}\\[3mm]
\dfrac{\partial L}{\partial \boldsymbol{\Lambda}} = \textrm{S}\VV-\RR_{in}^{\TT}{\rm M}_{in}\HH - \FF=\mathbf{0}.
\end{array}\right.
\end{equation}
By defining ${\rm Z} = {\rm D}{\rm S}^{-1}\RR_{in}^{\TT}{\rm M}_{in}$ and ${\rm W} = {\rm Z}^{\TT}{\rm Z}+a\,{\rm L}^{\TT}{\rm L}$ the reduced system reads
\begin{equation}\label{E:control_J}
{\rm W}\HH=\textrm{Z}^\TT(\dd - \DD {\rm S}^{-1}\FF),
\end{equation}
where ${\rm W}$ is the reduced Hessian associated with ${\rm L}$ and ${\rm Z} = \frac{{\rm d} ({\rm D}\VV)} {{\rm d} \HH}$ is the so-called \textit{sensitivity} matrix \cite{nocedal_numerical_2000}. Spectral properties of Z determine the conditioning and well-posedness of the problem.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Non-singularity of the reduced Hessian}\label{ss:well-posedness}
Sufficient and necessary conditions for the existence of a unique minimizer of problem \eqref{E:min_discr} are based on optimization theory.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Sufficient and necessary conditions}\label{S:wellp_1}
We analyze the existence of a unique minimizer for the optimization problem \eqref{E:min_discr} under suitable assumptions on number and location of the data. It is worth mentioning that the results we present are strictly related to the specific choice of the control variable; a different choice (e.g. the pressure drop, the inflow/outflow velocity, the outflow Neumann boundary condition, etc...) would imply different optimality conditions. 

According to the formulation presented in Section \ref{S:da_nse_DO} the \gls{DO} technique yields an algebraic optimization problem; necessary (sufficient) conditions for optimality are satisfied when the reduced Hessian matrix {\rm W} is positive semi-definite (definite). In our case, it is enough to show that ${\rm W}$ is non-singular, i.e. system \eqref{E:lagr_necessary} has a unique solution. 

Let $\vv\in\mathbb{R}^{N_u+N_p}$ be the vector such that 
\begin{equation}\label{E:v}
\vv = {\rm S}^{-1}\RR_{in}^{\TT}{\rm M}_{in} \xx,
\end{equation}
where $\xx\in \mathbb{R}^{N_{in}}$ is a generic non-null vector. Then, we can write  
\begin{equation}\label{E:posdef_product_comp}
\xx^{\rm T}{\rm W}\xx  = \vv^{\rm T} {\rm D}^{\rm T}{\rm D} \vv + \alpha\, \xx^{\rm T} {\rm L}^{\rm T}{\rm L} \xx.
\end{equation}
Since $\xx\neq\mathbf{0}$, ${\rm D}^{\rm T}{\rm D} \geq 0$ and  ${\rm L}^{\rm T}{\rm L}>0$, we conclude that using a regularization term, i.e. $\alpha>0$, \textit{always} yields the well-posedness of the problem, no matter how the measures are located in the domain. Since regularization might be computationally expensive in searching the optimal regularization parameter, we prove a non-singularity result also for $\alpha=0$. 

From now on, we assume $\alpha=0$ and we focus on the first term in \eqref{E:posdef_product_comp}. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{proposition}\label{P:suff_cond_nullrange}
For $\alpha=0$, ${\rm W}$ is non-singular $\Leftrightarrow$ 
\begin{equation}\label{E:suff_cond_nullrange}  
Null({\rm D}) \cap Range({\rm S}^{-1}\RR_{in}^{\TT}{\rm M}_{in}) = \{\mathbf{0}\}.
\end{equation}
\end{proposition}

\noindent {\it Proof} $\Leftarrow$ Take $\xx\neq\mathbf{0}$; then, $\vv={\rm S}^{-1} {\rm R}_{in}^{\TT} {\rm M}_{in} \xx \neq \mathbf{0}$ since ${\rm M}_{in}$ and ${\rm S}$ are non-singular and $\RR_{in}$ is full rank. Under the condition \eqref{E:suff_cond_nullrange} $\DD \vv\neq\mathbf{0}$; hence, $\xx^{\TT}{\rm W} \xx>0$ for all $\xx\neq\mathbf{0}$. $\Rightarrow$ By contradiction, assume \eqref{E:suff_cond_nullrange} does not hold. Then, there exists $\vv\neq\mathbf{0}$, defined as in \eqref{E:v}, such that ${\rm D}\vv=\mathbf{0}$, which implies that $\exists\, \xx\neq \mathbf{0}$ such that $\xx^{\TT}{\rm W}\xx=\mathbf{0}$, so ${\rm W}$ has a non-trivial kernel.

We need to analyze how the fulfillment of \eqref{E:suff_cond_nullrange} is subordinated to the reciprocal location (and number) of sites and nodes. In view of Proposition \ref{P:suff_cond_Qmatrix}, we prove the following Lemma.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{lemma}\label{L:Svw}
For $\vv$ defined as in \eqref{E:v}, ${\rm R}_{in}\vv\neq\mathbf{0} \;\;\forall\; \xx\in\mathbb{R}^{N_{in}}$ such that $\xx\neq\mathbf{0}$. 
\end{lemma}

\noindent {\it Proof}  Let  $\ww={\rm R}^{\TT}_{in}{\rm M}_{in} \xx$; then, definition \eqref{E:v} implies ${\rm S}\vv = \ww$. Without loss of generality, assume that inflow \gls{DOF}s are numbered in the first $N_{in}$ entries such that 
\begin{equation}
\ww=\left[\begin{array}{c}
{\rm M}_{in} \xx\\
\mathbf{0}
\end{array}\right].
\end{equation}
 We rewrite matrix S and vector $\vv$ accordingly as
\begin{equation}\label{E:rewrittenA}
{\rm S}=\left[\begin{array}{c c}
{\rm A}_{11} & {\rm A}_{12} \\ 
{\rm A}_{21} & {\rm A}_{22} \\
\end{array}\right], 
\qquad  \qquad
\vv=\left[\begin{array}{c}
\vv_1 \\
\vv_2
\end{array}\right],
\end{equation}
where A$_{11}\in\mathbb{R}^{N_{in},N_{in}}$ and $\vv_1 \in N_{in}$. Note that A$_{22}$ is the matrix of a generalized Stokes problem with homogeneous Dirichlet boundary conditions on $\Gamma_{in}\cup\Gamma_{wall}$ and Neumann conditions on $\Gamma_{out}$, so it is non-singular \cite{quarteroni_numerical_1994}. Assume by contradiction that ${\rm R}_{in}\vv=\mathbf{0}$. This implies that $\vv_1 = \mathbf{0}$. System ${\rm S}\vv=\ww$ reads:
\begin{equation}\label{Svw_rewritten}
\left[\begin{array}{c c}
{\rm A}_{11} & {\rm A}_{12} \\ 
{\rm A}_{21} & {\rm A}_{22}
\end{array}\right]\left[\begin{array}{c}
\mathbf{0}\\
\,\vv_2
\end{array}\right]=
\left[\begin{array}{c}
{\rm M}_{in} \xx\\
\mathbf{0}
\end{array}\right].
\end{equation}
We get A$_{22}\vv_2 = \mathbf{0}$ so that $\vv_2 = \mathbf{0}$, which implies that ${\rm M}_{in} \xx = \mathbf{0}$. This is a contradiction since ${\rm M}_{in}$ is non-singular and $\xx\neq \mathbf{0}$. \\

\begin{remark}\label{R:lemma}
Note that Lemma \ref{L:Svw} is not strictly specific for the Stokes problem; in fact, it typically holds also for \gls{FE} discretization of several differential problems, e.g. elliptic \gls{PDE}s. As a matter of fact, we only require that the discretization matrix corresponding to the homogeneous Dirichlet problem on $\Gamma_{in}$ is non-singular. Despite this, numerical experiments in Section \ref{ss:da_nse_lin_results} show that the sufficient condition, in the subsequent Proposition \ref{P:suff_cond_Qmatrix}, is sharp.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{proposition}\label{P:suff_cond_Qmatrix}
A sufficient condition for the existence of a unique minimizer for problem \eqref{E:min_discr} is that there exists a subset of sites such that the restriction $\QQ_r$ of $\QQ$ to such sites fulfills 
\begin{equation}\label{E:suff_cond_Qmatrix}
{\rm Q}_r (I-{\rm R}_{in,u}^{\TT}{\rm R}_{in,u})=\OO \qquad \hbox{and} \qquad rank({\rm Q}_r) = N_{in}.
\end{equation}
\end{proposition}

\noindent {\it Proof} Notice that \eqref{E:suff_cond_Qmatrix} implies that ${\rm Q}_r{\rm R}^{\TT}_{in,u}$, where ${\rm R}_{in,u}\in\mathbb{R}^{N_{in},N_u}$, has full rank; in fact 
\begin{equation}\label{E:fullrank}
\begin{array}{l}
N_{in} = rank({\rm Q}_r) = rank({\rm Q}_r{\rm R}^{\TT}_{in,u}{\rm R}_{in,u})  
\leq \min\{rank({\rm Q}_r{\rm R}^{\TT}_{in,u}), rank({\rm R}_{in,u})\}\\ 
[3mm]
\qquad = \min\{rank({\rm Q}_r{\rm R}^{\TT}_{in,u}),\,N_{in}\}, 
\end{array}
\end{equation}
hence $rank({\rm Q}_r{\rm R}^{\TT}_{in,u})$ cannot be less than $N_{in}$. Let us denote by $\DD_r$ the sub-matrix of D corresponding to $\QQ_r$, so that $\DD_r=[{\rm Q}_r \;\; \OO]$. 

We prove that, under condition \eqref{E:suff_cond_Qmatrix}, for $\vv={\rm S}^{-1}{\rm R}^{\TT}_{in}{\rm M}_{in}\xx$ and $\xx\neq\mathbf{0}$, $\DD_r\vv\neq\mathbf{0}$. 

As a consequence of Lemma \ref{L:Svw}, we have that ${\rm R}_{in}\vv \neq \mathbf{0}$, this implies ${\rm R}_{in,u}\vv_u \neq \mathbf{0}$, where $\vv_u$ is the restriction of vector $\vv$ to velocity \gls{DOF}s. Consequently, 
\begin{equation}\label{E:Dv}
{\rm D}_r \vv= {\rm Q}_r\vv_u = {\rm Q}_r (I - {\rm R}_{in,u}^{\TT}{\rm R}_{in,u} + {\rm R}_{in,u}^{\TT}{\rm R}_{in,u})\vv_u = {\rm Q}_r{\rm R}_{in,u}^{\TT}\,{\rm R}_{in,u}\vv_u.
\end{equation} 
Since ${\rm Q}_r{\rm R}_{in,u}^{\TT}$ has rank $N_{in}$, D$_r\vv\neq\mathbf{0}$ and the thesis follows.\\

Proposition \ref{P:suff_cond_Qmatrix} corresponds to design the selection matrix D such that its restriction to sites on $\Gamma_{in}$ has $N_{in}$ independent rows. As an example, in Figure \ref{F:meas_bounds}, we report two possible configurations of sites selection using P1bubble \gls{FE} for the velocity. Note that, with this choice of \gls{FE} spaces, sufficient conditions are satisfied when we have sites on each DOF on $\Gamma_{in}$. This fact opens the question of mesh refinement. For $\Delta\rightarrow 0$ we gain accuracy but, in the limit, we loose well-posedness (since $\Delta = 0$ corresponds to having an infinite number of sites on $\Gamma_{in}$). However, the presence of the noise makes the refinement, after a certain level, unnecessary since the accuracy is determined by the amount of noise only. More in general, this suggests that the selection of the parameter $\alpha$ for the regularizing term should depend on $\Delta$.
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{remark}
Conditions for the existence and uniqueness of the optimal solution might be unexpected. However, we notice that in the limit case of data given on \gls{DOF}s of $\Gamma_{in}$ only (i.e. ${\rm D} = {\rm R}_{in}$) the problem is equivalent to the solution of the linearized \gls{NSE} where Dirichlet boundary conditions on the inflow section are prescribed in a weak sense.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{remark}
Proposition \ref{P:suff_cond_Qmatrix} states that it is mandatory to have velocity measures at the inflow section; in our application this is not prohibitive. In fact, it is always possible to collect velocity data on internal layers transversal to the blood flow. Nevertheless, we cannot expect data to satisfy sufficient conditions for optimality (i.e. to be located on velocity \gls{DOF}s); in the next section we present a technique to overcome this problem.
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{remark}
Lemma \ref{L:Svw} shows that, given a non-homogeneous Neumann boundary condition on $\Gamma_{in}$, the corresponding computed velocity restricted to inflow \gls{DOF}s, $\vv_{in}$, has at least one non-zero entry. In order to relax conditions \eqref{E:suff_cond_Qmatrix} (i.e. $Rank({\rm Q}_r)=\overline{N}$, $\overline{N}<N_{in}$), an upper bound to the number of zero entries of $\vv_{in}$ would be needed. However, this is not possible since we can always find a Neumann condition on $\Gamma_{in}$ such that the corresponding $\vv_{in}$ has only one non-zero entry (consider the Dirichlet-to-Neumann map with $\left.\vv\right|_{\Gamma_{in}}=\vv_{in}$).  
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=0.30]{sufficient_m.png}
\caption{Example of sites location sufficient for theoptimality on grid nodes (left) and not on grid nodes (right) using P1bubble FE for the velocity.}
\label{F:meas_bounds}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Forcing optimality via interpolation}\label{S:wellp_2}
The existence of a unique minimizer with a regularization term $\alpha>0$ does not suffer from any constraint on sites and nodes. Despite this, there can be some drawbacks, among which we mention the fact that it alters the solution of the problem (especially when $\alpha \gg 1$). Also, Tikhonov regularization requires the choice of the optimal parameter $\alpha$, which might be computationally expensive, e.g. the \gls{DP}, see next section, requires the iterative solution of linear systems of the form ${\rm W}\HH = \textrm{Z}^\TT(\dd - \DD {\rm S}^{-1}\FF)$ \cite{Hansen}. 
There is however another possible way for forcing the optimality, exploiting the result of Proposition \ref{P:suff_cond_Qmatrix}. 

Let us assume that some data are available at the inflow, not necessarily fulfilling sufficient conditions. We can extend the given data to the entire set of \gls{DOF}s on the inflow boundary by interpolation (e.g. piece-wise linear) of available data and Dirichlet boundary conditions on $\Gamma_{wall}$. In this way, inflow data are virtually extended to inflow \gls{DOF}s and condition \eqref{E:suff_cond_Qmatrix} is automatically fulfilled. 
More specifically, if we assume data $\widetilde{\dd}$ to be given, then, after interpolation, we can write $\dd = \Pi \widetilde{\dd}$, where $\Pi\in \mathbb{R}^{N_s+N_{in},N_s}$ is the interpolation matrix associated with piece-wise linear Lagrangian polynomials. 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Example} Assume $\Gamma_{in}$ to be parallel to the $y$ axis in the reference plane $(x,\,y)$, then, $\Pi$ is defined as follows. 
\begin{equation}\label{E:interp_matrix}
\begin{array}{ll}
\displaystyle \left[\Pi\right]_{i,j}= -\frac{y^m_{j+1}- y^{DOF}_i}{\Delta_j}; & \displaystyle \left[ \Pi\right]_{i,j+1}= \frac{y^m_{j}-y^{DOF}_i}{\Delta_j} \\
[4mm]
\;\;\;\quad\quad{\rm for}\;\; y_i^{DOF}\in [y_j^m,\,y_{j+1}^m] & i = 1,\,...,\,N_{in},\; j = 1,\,...,\,N_{s,in}; \\
[5mm]
\left[\Pi\right]_{i+N_{in},i}= 1 &   i = 1,\,...,\,N_s.\\
\end{array}
\end{equation}
Here, $y^{m}_i$ is such that $\xx_i^m = (x_i^m,\,y_i^m)$ lies on $\Gamma_{in}$; $N_{s,in}$ is the number of sites on $\Gamma_{in}$ and $\Delta_j = y_{j+1}^m - y_{j}^m$. Also, $y^{DOF}_i$ such that $\xx_i^{DOF} = (x_i^{DOF},\,y_i^{DOF})$ corresponds to a DOF on $\Gamma_{in}$. \\

The interpolated data are clearly affected by the interpolation error. More precisely, if we assume that
\begin{equation}\label{E:interp_noise_distrib_eps}
\widetilde{\dd} = \dd^{sparse}=\uu_{exact} + \boldsymbol{\varepsilon},
\end{equation}
where $\boldsymbol{\varepsilon}$ is the random noise; then, on boundary \gls{DOF}s the interpolated data can be written as
\begin{equation}\label{E:interp_noise_distrib_eta}
\dd^{interp}=\uu_{exact} + \boldsymbol{\eta}_{\boldsymbol{\varepsilon}},
\end{equation}
where $\boldsymbol{\eta}_{\boldsymbol{\varepsilon}}$ is a linear combination of noise and exact velocity on the original sites. More precisely, the statistical features of $\boldsymbol{\eta}_{\boldsymbol{\varepsilon}}$ are affected by the interpolation error and the additional error that occurs when interpolating noisy data \cite{rivlin2003}.

In this way, the functional to minimize reads
\begin{equation}\label{E:int_sparse}
 \displaystyle\frac{1}{2}\|{\rm D}^{sparse}\VV - \dd^{sparse}\|^2_2 + \frac{1}{2}\|{\rm D}^{interp}\VV - \dd^{interp}\|^2_2,
\end{equation}
where D$^{sparse}$ and D$^{interp}$ select respectively components of the velocity field on sparse sites and \gls{DOF}s on $\Gamma_{in}$. Representation \eqref{E:int_sparse} highlights the role of interpolation as a regularizing technique. 
In Section \ref{ss:da_nse_lin_results} we investigate this strategy and compare it with the classical Tikhonov regularization. A careful investigation of the effect of the interpolation error in statistical terms will be carried out in Chapter \ref{C:da_bayesian}. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical results}\label{ss:da_nse_lin_results}
In this section we test our technique on 2-dimensional cases. We consider both noise-free and noisy data, the former for a consistency analysis, the latter for a discussion on the reliability of the \gls{DA} process. At this preliminary stage data are synthetic, i.e. they do not come from measures, but from a given analytical solution suitably perturbed with additional noise. We present results for a validation of Proposition \ref{P:suff_cond_Qmatrix}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Implementation details}\label{p:implem_details}We implement the \gls{FE} method with choice of compatible \gls{FE} spaces P1bubble-P1 for velocity and pressure respectively. The mesh is generated with the software {\tt FreeFem++}. For Tikhonov regularization we use ${\rm L} = \gr_d$, the discretized gradient, and we generate the optimal parameter $\alpha$ by means of \gls{DP} \cite{Hansen}. Here, the optimal $\bar{\alpha}$ is such that $\|{\rm D} \VV_{\bar{\alpha}}-\dd\|^2_2 - \tau^2$ is minimized, being $\VV_{\bar{\alpha}}$ the optimal solution computed with $\bar{\alpha}$ and $\tau^2=\|\boldsymbol{\varepsilon}\|^2_2$, the norm of the noise. In practical applications, the \gls{DP} method is not commonly used for the search of the optimal parameter since the norm of the noise, $\tau^2$, is usually unknown. Since we have an estimate of the noise variance $\sigma^2$ (for each component of the noise vector), we use $\tau^2\cong d N_s \sigma^2$.

For the solution of the reduced system \eqref{E:control_J} and of the state and adjoint systems we use inexact Krylov subspace iterations; more specifically, we solve the reduced Hessian matrix $\rm W$ with the \gls{GMRES}R method \cite{gmresr, fgmres_gmresr}, which allows to change the preconditioner at each iteration, and the matrices ${\rm S}$ and ${\rm S}^\TT$ with the \gls{GMRES} method (inner iterations). For such nested iterations it is possible to bound the inner tolerance so that the global residual norm falls below a prescribed tolerance. This issue has been investigated in \cite{simoncini} and \cite{iteration_accuracy} for several applications in scientific computing. In the numerical simulations presented in this work outer and inner tolerances are chosen empirically; in view of applications to real data we plan to take into account an adaptive choice of the inner tolerance.

The bottleneck is the solution of linear systems associated with $\rm S$ and ${\rm S}^\TT$ where an efficient preconditioner is required; we solve these systems monolithically with the pressure-convection-diffusion (\gls{PCD}) preconditioner proposed by Kay in \cite{kay} and by Silvester et al. in \cite{silvester}. Table \ref{T:gmres_newton} reports stopping criteria used for each linear system involved. ${\rm P}({\rm S})$ stands for the preconditioned system, $\mathbf{r}_k$ is the residual at iteration $k$ and $\mathbf{rhs}$ is the right hand side of the current system. Note that these parameters refer to the solution of problems featuring non-trivial geometries or low viscosity values; in fact, such problems require an accurate solution of inner systems. 

The \gls{PCD} preconditioner exploits an approximation of the Schur complement based on considerations on the continuous operators involved. We report its formulation, details regarding motivation and performance can be found in \cite{elman_finite_2005}. We define the preconditioner P as follows: 
\begin{equation}\label{E:prec}
{\rm P}={\rm L}_p{\rm F}_p^{-1}{\rm M}_p.
\end{equation}
Here, ${\rm L}_p$ is the \gls{FE} discretization of the pressure Laplacian
\begin{equation}\label{E:Lp}
[\textrm{L}_p]_{i,j}=\int\limits_{\dom}\gr \psi_j\cdot\gr \psi_i;
 \end{equation}
being $\psi$ the pressure basis functions. ${\rm F}_p$ is the discretization of the momentum equation for the pressure variable
\begin{equation}\label{E:Cp}
[{\rm F}_p]_{i,j}=\mu\int\limits_{\dom}\gr \psi_j\cdot\gr \psi_i + \int\limits_{\dom} (\bb \cdot \gr \psi_j)\,\psi_i.
 \end{equation}
${\rm M}_p$ is the pressure mass matrix and $\bb$ is the advection field. Note that, using an iterative solver, we only specify the matrix-vector product for ${\rm P}^{-1}\mathbf{v}$; hence, we only require solvers for the pressure Laplacian and the mass matrix. As reported in \cite{elman_finite_2005} we can assemble these matrices prescribing natural boundary conditions on all boundaries but the inflow section, where a Dirichlet boundary condition is required. When solving the adjoint operator ${\rm S}^{\TT}$, we prescribe Dirichlet boundary conditions on the outflow section.

Numerical results are obtained with the C++ \gls{FE} library {\tt lifeV} and post-processed with {\tt ParaView}, see Section \ref{S:preliminary_software}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ |c | c | c | c|}
        \hline
	linear system		& solver	& tol		& stopping criteria\\
	\hline
	${\rm C}+{\rm A}$	& \gls{GMRES} 	& 1.e-9  	&  $\|\mathbf{r}_k\| / \|\mathbf{r}_0\|$\\
	\hline
	${\rm P}({\rm S})$ 	& \gls{GMRES}     & 1.e-9     &  $\|\mathbf{r}_k\| / \|\mathbf{rhs}\|$\\
    \hline
	$\rm W$           	& \gls{GMRES}R     & 1.e-6     &  $\|\mathbf{r}_k\| / \|\mathbf{rhs}\|$\\ 
	\hline
\end{tabular}
\caption{Parameters setting for the iterative solvers (implemented in the {\tt AZTEC} library).}
\label{T:gmres_newton}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Analytic test case}
Simulations presented in this section are based on the 2-dimensional Navier-Stokes flow in the domain $\dom = [-0.5,\,1.5]\times[0,\,2]$ whose analytic solution, $\uu_{anl}$, reads 
\begin{equation}\label{E:exact_ns}
\left\{\begin{array}{l}
[\uu]_1(x,y)=1-\displaystyle e^{\lambda x}\cos(2\pi y),\\
[1mm]
[\uu]_2(x,y)=\displaystyle \frac{\lambda}{2\pi}e^{\lambda x}\sin(2\pi y),\\
[1mm]
p(x,y)=\displaystyle \frac{1}{2}e^{2\lambda x}-\frac{1}{4\lambda}e^{2\lambda}.
\end{array}\right.
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.5\linewidth} 
\centering
\includegraphics[scale=0.5]{vfield035.png}
\end{minipage}
\hspace{1cm} 
\begin{minipage}[b]{0.3\linewidth}
\centering
\includegraphics[scale=0.5]{vfield005.png}
\end{minipage}
\caption{Exact velocity field in correspondence of $\nu$ = 0.035 (left) and $\nu$ = 0.005 (right).}
\label{F:velocity_nu}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Here, the kinematic viscosity of the fluid is $\nu$ = 0.035, a value featured by blood (lower parameters induce advection-dominated problems, see Figure \ref{F:velocity_nu}, right). The adimensional parameter $\lambda$ is such that $\lambda=\frac{1}{2}(\nu^{-1}-\sqrt{\nu^{-2}+16\pi^2})$. In Figure \ref{F:velocity_nu} we report the analytical velocity computed on a fine grid in correspondence of $\nu$ = 0.035 (left) and 0.005 (right). In these test cases the index of accuracy is the relative error with respect to the analytic solution, $E_U = \frac{\|\UU - \UU_{anl}\|_2}{\|\UU_{anl}\|_2}$, where $\UU_{anl}$ is the discretization of $\uu_{anl}$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Data Generation} \label{p:data_gen}We assume data to be given on the inflow boundary and on some internal points in the computational domain, not necessarily in correspondence of grid points. In the case of noise-free measures the data vector coincides with the exact solution, for the noisy case data are generated by adding Gaussian noise to the exact solution.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Noise-free data}\label{p:noise-free}We consider the linearized problem with $\boldsymbol{\beta}=\uu_{anl}$ and we analyze the error with respect to the exact solution. In this test case we set $\alpha=0$ and we consider data on \gls{DOF}s on $\Gamma_{in}$ and on three internal layers (in correspondence of grid nodes at $x$ = 0, 0.5, 1). In Figure \ref{F:lin_0noise} we report the relative error $E_U$ in logarithmic scale versus the discretization step $\Delta$. Here, $\UU_{anl}$ is the discretization of the analytic solution. Consistently with the C\'ea Lemma \cite{quarteroni_numerical_1994} we observe the expected convergence rate (for the chosen \gls{FE} spaces), $E_U= \mathcal{O}(\Delta^2)$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=0.3]{eVSh_0noise.pdf}
\caption{Relative error versus $\Delta$ for noise-free data, a reference curve, $\Delta^2$, is reported.}
\label{F:lin_0noise}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | c | c | c | c | c | c | c |  }
\hline    
        \gls{BS}		&	 \gls{IS}		& 	$K_1$		&	    $E_U$   &	\verb"it"	& 	$K_1$(reg)	&		$E_U$(reg)	&	\verb"it"(reg)	\\
       \hline
       $27$		&	0		&  	3.1e+05   	& 		0.0500	    &	19	&  	1.1e+04   	& 		0.0492	        &	16	\\

       $13$		&	0		&  	1.3e+20		&		0.0756		&	13	&  	9.7e+03  	&		0.0567			&	13	\\
	
       $9$		&	0		&  	7.9e+19		&		0.1957		&	15	&  	9.8e+03		&	 	0.1050			&	17	\\

       $27$		&	100		&  	1.3e+06		&		0.0367		&	15	&  	3.0e+04		&		0.0354			&	12\\
	
       $13$		&	114		&  	1.1e+11		&		0.0443		&	11	&  	4.6e+04 	&		0.0379			&	11	\\
        
       $9$		&	118		&  	1.5e+11		&		0.0478		&	13	&  	6.1e+04	    &		0.0348			&	12\\
       \hline
\end{tabular}
\caption{Condition number, relative error and number of iterations for  different choices of sites with SNR = 20.}
\label{T:cond1}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Verification of the optimality condition}\label{p:wellp} With the purpose of validating and testing the sharpness of Proposition \ref{P:suff_cond_Qmatrix}, we focus on the conditioning of the problem and we inspect singular values of ${\rm W}$. We choose $\Delta$ = 0.076 and \gls{SNR} = 20. This value corresponds to adding a low amount of noise and it is chosen for testing the method. Later on in the chapter realistic values, such as \gls{SNR} = 8, are also considered. In this set of simulations $N_{in}=27$. 
In Table \ref{T:cond1} we report the condition number $K_1$ of ${\rm W}$, for $\alpha =$ 0 and $\alpha =$ 1e-2, for different values of $N_s$; on the inflow boundary sites (\gls{BS}) correspond to grid nodes, while inside the domain they do not necessarily. The value of the regularization parameter is fixed a priori on an empirical basis. In Figures \ref{F:sv_noreg_int} and \ref{F:sv_noreg_int2} we report in a semilogarithmic scale the singular values of the reduced Hessian ${\rm W}$, properly rescaled, with and without regularization. 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.40 \textwidth]{sv_27_0.pdf} &  \includegraphics[width=0.4\textwidth]{sv_13_0.pdf} \\
\includegraphics[width=0.40\textwidth]{sv_27_60.pdf} & \includegraphics[width=0.40\textwidth]{sv_13_60.pdf} \\
\end{tabular}
\end{center}
\caption{Singular values of ${\rm W}$ with $\alpha=0$ for IS =0 and 60 and BS = 27 and 13.}
\label{F:sv_noreg_int}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.40\textwidth]{sv_27_0reg.pdf} &  \includegraphics[width=0.4\textwidth]{sv_13_0reg.pdf} \\
\includegraphics[width=0.40\textwidth]{sv_27_60reg.pdf} & \includegraphics[width=0.40\textwidth]{sv_13_60reg.pdf} \\
\end{tabular}
\end{center}
\caption{Singular values of ${\rm W}$  with $\alpha$ = 1e-2 for IS = 0 and 60 and BS = 27 and 13.}
\label{F:sv_noreg_int2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us focus on the case of data {\it not} satisfying optimality conditions. When the number of internal sites (\gls{IS}) = 0 and $\alpha = 0$ the distribution of the singular values features a gap and we can classify the problem as \emph{rank-deficient} (see \cite{Hansen}). In these cases it is not worth using Tikhonov regularization; in fact, in a certain sense, regularization induces a shift of the eigenvalues (see Figure \ref{F:sv_noreg_int2}, top right). It is rather effective to identify such gap using methods such as truncated singular value decomposition (TSVD) or modified TSVD, see \cite{Hansen}. In our case the system is solved by means of the \gls{GMRES} method setting empirically a tolerance so to stop iterations before the gap.
When conditions are not satisfied, \gls{IS} $\neq 0$ and $\alpha=0$ there is no gap in the singular values distribution and the problem can be classified as \emph{discrete ill-posed} \cite{Hansen}. Also, the addition of internal data does not improve significantly the condition number and the relative error.

We also investigate the behavior of the procedure for different values of the kinematic viscosity $\nu$, i.e. for different Reynolds number. The choice of $\nu_1=0.1$ and $\nu_2=0.01$ is driven by the fact that $\nu_2$ induces ill-conditioning of the forward problem associated with the matrix $\rm S$. 

From results reported in Table \ref{T:cond2} we can infer two considerations. With high Reynolds number the accuracy of the recovered solution is worse, this is due to numerical instability which arises with problems dominated by advection phenomena. Despite this, the condition number of the control problem is lower than in the stable case. This is consistent with the theory presented by Lions in \cite{LionsPNAS} and by Gunzburger in \cite{gunzburger_perspectives_2002}, Chapter 1. In these works the key idea is that the more a system is unstable (and more sensitive to perturbations), the simplest, or the cheapest, it is to achieve exact or approximate controllability.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| c | l | c | c | c |}
\hline    
        \gls{BS}		& 	\gls{IS}             	& $K_1$				&       $E_U$ 		&	\verb"it"	\\
        \hline
        27	&   0					&  	 1.3e+06   		& 		0.0380	    &	17	\\
        
        13	&	0					&  	 2.0e+21		&		0.0802		&	12	\\
	
        9	&	0					&  	 5.7e+21		&		0.2566		&	15	\\

        27	&	27{\small $_{in}$}	&  	 2.5e+06		&		0.0303		&	15	\\
        	
        27	&	27{\small$_{out}$}	&  	 1.6e+06		&		0.0348		&	14	\\

        27	&	27{\small$_{mid}$}	&  	 1.9e+06		&		0.0331		&	14	\\    
        \hline
\end{tabular}
\hspace{0.cm}
\begin{tabular}{| c | l | c | c | c |}
	\hline    
    \gls{BS}		& 	\gls{IS}          		& $K_1$		&       $E_U$  		&	\verb"it"	\\
	\hline
	27		&	0					& 1.8e+05	&		0.0616		&	21	\\
	
    13		&	0					& 4.0e+19	&		0.0877		&	17	\\
        
    9		&	0					& 7.1e+19	&		0.1372		&	14	\\
	
    27		&	27{\small $_{in}$}	& 3.7e+05	&		0.0574		&	19	\\
	
    27		&	27{\small$_{out}$}	& 3.7e+05	&		0.0498		&	19	\\

	27		&	27{\small$_{mid}$}	& 3.7e+05	&		0.0461		&	18	\\    	
	\hline
\end{tabular}
\caption{Condition number, relative error and number of iterations with SNR = 20 for different choices of sites and viscosity: $\nu$ = 0.1 (left) and $\nu$ = 0.01 (right). Subscripts $in$, $out$ and $mid$ refer respectively to IS close to the inflow, outflow and middle of the domain.}
\label{T:cond2}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%
Also, results of Table \ref{T:cond2} are consistent with the result reported in \cite{Veneziani20071685} which states that the noise affecting the data decays exponentially with respect to the distance from the inflow boundary. With a low Reynolds number the introduction of sites located close to $\Gamma_{in}$ improves the accuracy; for high Reynold number the solution of the forward problem has a slower decay as we move far form $\Gamma_{in}$; this fact implies that internal sites significantly determine the accuracy of the solution. 
\begin{remark}
When the value of the kinematic viscosity induces advection-dominated phenomena we can stabilize the formulation adding artificial terms to the weak form of the problem, see \cite{quarteroni_numerical_1994} chap. 6. These terms modify the matrix $\SS$ and in some cases they might also circumvent the inf-sup condition so to allow the choice of incompatible \gls{FE} spaces. Using a simple and straightforward method we stabilize the problem and check possible effects on the \gls{DA} process. We consider the {\it Streamline Diffusion} method; we add artificial viscosity in the predominant direction of the advection field. Formally we add the term
\begin{equation}\label{E:streamline}
s(\uu,\vv)=\xi(\Delta,\bb)\displaystyle\int\limits_{\dom}\left(\bb\cdot\gr\uu\right)\left(\bb\cdot\gr\vv\right)\;\textrm{d}\xx;
\end{equation}
where the stabilization coefficient $\xi$ is defined as $\xi(\Delta,\bb)=\frac{\Delta}{2\|\bb\|_2}$. According to the Strang Lemma for generalized Galerkin methods (see \cite{quarteroni_numerical_1994}), the expected convergence rate for the velocity discretization error is only linear. Figure \ref{F:streamline} displays results obtained with $\nu$ = 0.005, \gls{SNR} = $\infty$ (left) and \gls{SNR} = 8 (right) and data on $\Gamma_{in}$ satisfying well-posedness conditions. In case of noise-free data, the discretization error $E_U$ versus $\Delta$, as expected, has linear convergence rate. In the noisy case, convergence rate is of the order of $\mathcal{O}(\sqrt{\Delta}) \sim \mathcal{O}(\sqrt{N_s^{-1}})$, as observed for the Stokes problem in Chapter \ref{C:preliminary}. Again, we note that, for a given set of $N$ independent identically distributed random variables, the standard deviation is proportional to $N^{-0.5}$; this explains our results. We conclude that including artificial stabilizing term does not compromise the assimilation process, indeed it improves the conditioning of the forward problem.  
\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.45\linewidth} 
\centering
\includegraphics[scale=0.27]{eVSh_0noiseSTAB.pdf}
\end{minipage}
\hspace{1cm} 
\begin{minipage}[b]{0.3\linewidth}
\centering
\includegraphics[scale=0.26]{eVSh_12noiseSTAB.pdf}
\end{minipage}
\caption{Relative errors versus $\Delta$ for noise-free (left) and noisy data (right) problem. Reference curves, $\Delta^2$ and $\Delta$ (left) and $\sqrt{\Delta}$ (right), are reported for a comparison.}
\label{F:streamline}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The interpolation strategy}\label{p:interpolation} We report numerical results on the interpolation technique addressed in Section \ref{S:wellp_2}; our main aim is to investigate its reliability and competitiveness with respect to the usual Tikhonov regularization. Also, we point out some facts related to the dependence of the conditioning of the problem on sites location as a further confirmation of results presented in Section \ref{S:wellp_1}.

In Table \ref{T:interp_reg} in correspondence of different choices of location and number of sites we report the relative error and the number of \gls{GMRES} iterations for two different grids. Here, in the first column, ``grid'' stands for ``points located on grid nodes'', i.e. on \gls{DOF}s; note that $N_{in}$ = 12 (right) and 27 (left). The case where neither Tikhonov regularization or interpolation are performed is reported as a reference test. $\alpha_{DP}$ is the optimal regularization parameter obtained using \gls{DP}. The activation of the interpolation process is determined by the flag \verb"i": when 0 no interpolation is performed and $\dd$ corresponds to given data; when 1, $\dd$ corresponds to both the original data and the interpolated values on grid nodes. From these results, we can infer the following considerations.
%According to the results presented in Section \ref{S:wellp_1}, adding a regularization term yields the well-posedness and a better accuracy.

The interpolation process has the effect of yielding the well-posedness of the problem, 
%as we can notice by comparing the case \verb"i"=0 with \verb"i"=1; 
moreover, interpolation can improve both the conditioning and the accuracy of the computed solution. We note that interpolation is a competitive technique with respect to common regularization methods in terms of accuracy. Tikhonov regularization is not a consistent technique, even in absence of noise; in fact, it modifies the nature of the problem adding an artificial term to the formulation. As opposed to this, using interpolation with noise free data, the assimilated solution is as accurate as a \gls{FE} solution of the forward problem with exact boundary conditions. Furthermore, it requires much less computational time. In fact, we recall that the \gls{DP} requires the iterative solution of reduced systems for retrieving the optimal parameter $\alpha_{DP}$.

In order to get the picture of what has been implemented, in Figure \ref{F:8data_solved} we report the sparse and interpolated data (left) and the corresponding assimilated solution (right). 
%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{remark}
%Here we look for the optimal $\bar{\alpha}$ such that $\|{\rm D} \VV_{\bar{\alpha}}-\dd\|^2_2 - \tau^2$ is minimized, being $\VV_{\bar{\alpha}}$ the solution computed with $\bar{\alpha}$ and $\tau^2=\|\boldsymbol{\varepsilon}\|^2_2$, the norm of the noise, usually unknown. 
%\end{remark}
%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{| l | c | c | c | c |  }
    	\hline
        \gls{BS}	& 	\verb"i"	& 	 $\alpha_{DP}$	& $E_U$		&	\verb"it"	\\
        \hline
        11 (grid)		&	0			&	   -		&  	0.108 	& 	14	\\        
        5 $\;\;$(grid)	&	0			&	   -		&  	0.831	&	9	\\	
        6 				&	0			&	   -		&  	 0.182	&	10	\\	
        6 				&	0			&	   0.047	&  	 0.150	&	11	\\    
    	6 		 		&	1			&	   -		&  	 0.167	&	14	\\        
 %      6    			&	1			&	0.047$^*$   &  	 0.161	&	14	\\	
        3 $\;\;$(grid)	&	0			&	   -		&  	 0.888	&	6	\\	
        4  				&	0			&	   -		&  	 0.312	&	8	\\    
        4  		 		&	0			&	   0.056	&  	 0.214	&	14	\\    
	    4 		  		&	1			&	   -		&  	 0.262	&	13	\\         
%		4 		  		&	1			&	0.056$^*$	&  	 0.223	&	13	\\
    	\hline     
         \end{tabular}
\hspace{0.cm}
\begin{tabular}{| l | c | c | c | c |  }
    	\hline
        \gls{BS}	& 	\verb"i"	& 	 $\alpha_{DP}$	& $E_U$	&	\verb"it"	\\
        \hline
        27 (grid)		&	0			&   -			& 0.050 & 	19	\\
        13 (grid)		&	0			&   -			& 0.076	&	13	\\
        14 				&	0			&   -			& 0.068 &	14	\\
        14				& 	0			&   0.021		& 0.061 &	15	\\
	    14 	 			&	1			&   -			& 0.059	&	18	\\
        %14         	&	1			&   0.021$^*$	& 0.056	&	16	\\
        9 $\;\;$(grid)	&	0			&   -			& 0.196	&	15	\\
        8  				&	0			&   -			& 0.199	&	11	\\
        8  				&	0			&	0.038		& 0.137	&	18	\\
        8	  			&	1			&	-			& 0.139	&	17	\\
	 	%8		  		&	1			&	0.038$^*$	& 0.129	&	17	\\
         \hline
\end{tabular}
\caption{Relative error and number of GMRES iterations (GMRES tolerance tol=1.e-6) in correspondence of $\Delta = 1/6$ (left) and $1/14$ (right). IS = 0 and SNR = 20.}
\label{T:interp_reg}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.3\linewidth} 
\centering
\includegraphics[scale=0.32]{8interp8.png}
\end{minipage}
\hspace{0.01cm} 
\begin{minipage}[t]{0.7\linewidth}
\centering
\includegraphics[scale=0.32]{pvd5_8interp8_reg.png}
\end{minipage}
\caption{On the left, actual data used in the DA process (original and interpolated). On the right, pressure and velocity regularized solution and data.}
\label{F:8data_solved}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Iterative procedure for the nonlinear problem}\label{S:da_nse_nonlinear}
We consider the nonlinear advection term $(\uu\cdot\gr)\uu$ and we solve a nonlinear \gls{PDE} constrained optimization problem. In literature several methods are available for the solution of the discretized problem; among those we mention the Newton method \cite{nocedal_numerical_2000} and its variants (common optimization methods will be addressed in future works with comparison purposes). Here, we address this issue by combining the \gls{DA} procedure for the linear case presented in the previous section with classical fixed point linearization schemes. In particular, we refer to Picard and Newton methods \cite{quarteroni_numerical_1994}. The \gls{DA} assimilation problem is solved iteratively as follows. Let $\UU_k$ be a given guess for the velocity field at iteration $k+1$. The discrete iterative procedure reads
\begin{equation}\label{E:algorithm}
\begin{array}{l}
\min\limits_{\HH_k}  \displaystyle\frac{1}{2}\|{\rm D}\VV_{k+1} - \dd\|^2_2 + \displaystyle\frac{\alpha}{2}\|{\rm L} \HH_{k+1}\|^2_2 \\
[3mm]
\hbox{s.t.} \quad {\rm S}_k\VV_{k+1} = \RR_{in}^{\TT}{\rm M}_{in} \HH_{k+1} + \FF_k
\end{array}
\end{equation}
up to fulfillment of a convergence criterion (typically $\|\VV_k - \VV_{k+1}\|\leq \delta$, being $\delta$ a user defined tolerance chosen empirically in our case). We use the following notation
\begin{equation}\label{E:discr_lin_ns}
{\rm S}_k=\left[\begin{array}{cc}
{\rm C} + {\rm A}_k & {\rm B}^\TT\\
{\rm B} & {\rm O} 
\end{array}\right], \qquad \hbox{and} \qquad \FF_k = \FF + w \mathbf{Y}_k.
\end{equation}
Here, A$_k$ comes from the discretization of $(\overline{\uu}_k\cdot \gr)\uu_{k+1} + w (\uu_{k+1}\cdot \gr)\overline{\uu}_k$ ($w=0$ for Picard method, $1$ for Newton); $\mathbf{Y}_k$ is the discretization of $(\overline{\uu}_k\cdot \gr)\overline{\uu}_k$. Here $\overline{\uu}_k$ is defined as $\vartheta\uu_{k-1} + (1-\vartheta) \uu_k$, being $\vartheta\in[0,\,1]$, $w$ is a relaxation parameter, also chosen empirically. 
A necessary condition for the convergence of algorithm \eqref{E:algorithm} is the existence of a unique solution, at each iteration of the optimization problem, which has been proved in Section \ref{ss:well-posedness}. In Figure \ref{F:ns_convergence_delta_conv} (left) we show numerical evidence of the convergence: in a logarithmic scale we report, at each iteration of algorithm \eqref{E:algorithm} with $w$ = 0.8, the relative error, with respect to a reference solution, defined as $E_{U} = \frac{\|\UU_k - \UU_{anl}\|_2}{\|\UU_{anl}\|_2}$, $\UU_{anl}$ being as in \eqref{E:exact_ns}. For increasing number of iterations the graph shows a convergence to a limit value in correspondence of different values of \gls{SNR}. Also, numerical tests show that the number of iterations is almost $\Delta$-independent, where $\Delta$ is the discretization step used for the \gls{FE} approximation. Figure \ref{F:ns_convergence_delta_conv} (right) reports relative errors $E_{U}$ and number of iterations \verb"it" for different values of $\Delta$. In the table \gls{SNR} and $N_s$ are fixed.

Due to the local convergence of the Newton method for the forward problem, a common procedure is to perform a few Picard iterations and use the resulting velocity as initial guess (we refer to this method as Picard-Newton). The activation of the Newton method is driven by the parameter $w\in[0,\,1]$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{tabular}{ccc}
\begin{minipage}{.4\textwidth} 
%\centering
\includegraphics[scale = .4]{newton_eConv_legend.png} 
\end{minipage}
& \hspace{1.7cm} &
\begin{minipage}{.05\textwidth}
	\begin{center}
	\begin{tabular}{ |c | c | c |}
	  		\hline
			$\Delta$ 	& $E_{U}$ & \verb"it" \\
			\hline
			1.7e-1		& 0.0534 	& 9\\
			7.1e-2 		& 0.0288 	& 9\\
			5.0e-2 		& 0.0241 	& 9\\
			\hline
	\end{tabular}
	\end{center}
\end{minipage}
\end{tabular}
\caption{Left: relative error $E_{\UU}$ versus number of iterations in correspondence of different values of SNR. Right: Relative error $E_{\UU}$ and number of iterations of algorithm \eqref{E:algorithm} for SNR = 20 and $N_s$ = 160.}
\label{F:ns_convergence_delta_conv}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical results}
In this Section we test the \gls{DA} procedure for the \gls{NSE} on 2-dimensional and axisymmetric cases. With noise-free and noisy data, we conduct a consistency analysis, and we discuss the filtering properties of the method on an analytic test case. Synthetic data are generated either from a given analytical solution or a numerical solution of the forward problem. We also show results obtained on 2-dimensional non-rectangular domains, retrieving from the assimilated velocities a flow-related  variable of medical relevance, namely the wall shear stress (\gls{WSS}). Finally, we present results obtained with an axisymmetric formulation of the problem in a cylindrical domain.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Noise-free data}\label{p:noise-free} We analyze the error with respect to the analytic solution. We set $\alpha=0$ and consider data on \gls{DOF}s on $\Gamma_{in}$ and on three internal layers ($x= 0,\,0.5,\,1$).

In Table \ref{T:ns_newton} (left) we report for the Picard-Newton method the error $E_U$ and the number of Picard-Newton iterations (the number of Newton iterations to get convergence is reported in brackets) in correspondence of different mesh sizes $\Delta$. Same results are reported in Figure \ref{F:pvd_ns} (left) in a logarithmic scale with a reference quadratic curve, the expected \gls{FE} convergence rate is recovered. Moreover, we observe that the number of iterations needed for the convergence is almost mesh independent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ | c | c | c | }
    \hline
	$\Delta$ & $E_{U}$ 		& \verb"it"\\ \hline
	0.160    &   4.97e-2	&  5 (3)\\
	0.072 	 &   9.49e-3	&  4 (2)\\
	0.050	 &   4.73e-3	&  4 (2)\\
	\hline
\end{tabular}
\hspace{2.cm}
\begin{tabular}{| c | c | c |}
     \hline
	$N_s$	 	&  $E_{U}$ 		& \verb"it"\\ \hline
	 22 (11+11) &  9.01e-2		&  7 (3)\\
	 44 (22+22)	&  6.01e-2     	&  9 (5)\\
	 88 (44+44)	&  4.52e-2		&  7 (3)\\ 
	\hline
\end{tabular}
\caption{On the left results for noise-free simulations; on the right results for SNR =20 and $\Delta$ = 0.160, in the first column $N_s$ = (BS $+$ IS). In the third column, in brackets, the number of Newton iterations.}
\label{T:ns_newton}
\end{center}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=0.42\textwidth]{e_vs_dx_picnew.png} & 
\includegraphics[width=0.34\textwidth]{noisyData_assVelocity.png}
\end{tabular}
\end{center}
\caption{Left: relative errors versus $\Delta$, Picard-Newton iterations with noise free data. Right: assimilated velocity with noisy data (black) for SNR =20 in a region with horizontal analytical solution.}
\label{F:pvd_ns}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Noisy data}\label{p:noisy_data} We investigate the dependence of the relative error $E_U$ on the number of available measures and on the amount of noise. Regularization parameter and data are set as in the previous paragraph.
In Table \ref{T:ns_newton} (right) we report the discretization error and the number of iterations for the Picard-Newton method in correspondence of \gls{SNR} = 20 and $\Delta$ = 0.160 for different number of sites $N_s$ = \gls{BS} + \gls{IS}, where \gls{BS} and \gls{IS} are the number of boundary and internal sites. In this case $w$ = 0.8. In Figure \ref{F:ns_convergence_noise} (left) we report the average discretization error over 100 noise realizations, ${\rm E}_{100}$, and a reference curve $\mathcal{O}(N_s^{-0.5})$. Same statistical considerations reported in Section \ref{S:preliminary_furtherresults} and \ref{ss:da_nse_lin_results} may be inferred.
In Figure \ref{F:ns_convergence_noise} (right) we report in a logarithmic scale the average error versus the inverse of \gls{SNR} for \gls{SNR} = 33.3, 20, 10, 8.3 and $\Delta$ = 0.072. We can notice a linear behavior by a comparison with the dashed reference curve. In Figure \ref{F:convergence_realiz} we report the relative error $E_{U}$ of the sample mean of the velocity over $N_r$ noise realizations: $\overline{\UU} = \frac{1}{N_r}\sum_{i=1}^{N_r} \UU_i$, where $\UU_i$ is the velocity corresponding to the $i$-th noise realization. As expected, the error converges with rate -0.5 to the discretization error for noise-free data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.4\linewidth} 
\centering
\includegraphics[scale=0.38]{e_vs_Ns_12.png}
\end{minipage}
\hspace{0.1cm} 
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.38]{e_vs_SNR_ns.png}
\end{minipage}
\caption{Relative errors $E_U$ versus $\Delta$ with noisy data (left); relative errors $E_U$ versus SNR$^{-1}$ (right) for Picard-Newton iterations.}
\label{F:ns_convergence_noise}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[scale=0.35]{e_vs_Nr_28.pdf}
\caption{Relative error of $\overline{\UU}$ versus number of noise realizations $N_r$.}
\label{F:convergence_realiz}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In Figure \ref{F:pvd_ns} (right) we report a zoom of the computed field and the noisy data in a region with horizontal analytical solution; the presence of the noise is evident in the vertical components and in low magnitude velocities; this happens since the amount of noise does not depend on the local magnitude of the velocity vector field. In correspondence of these low magnitude values the recovered assimilated field differs significantly from the data and it is closer to the exact solution. This pinpoints the role of the assimilation procedure as a sophisticated way for filtering noise from the data, based on the mathematical model of the process at hand.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Towards real geometries}\label{S:da_nse_real_geometries}
In this section we consider more complex geometries, which are intended to resemble vessels of medical interest. In Test case I and II we reproduce in 2-dimensional curved computational domains a section of the carotid and of the aorta respectively. In test case III we consider a 3-dimensional cylindrical domain treated with an axisymmetric formulation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Test case I} We assume that some measurements of the blood velocity are available on $\Gamma_{in}$ and on sparse internal sites, see Figure \ref{F:carotid} (left) for domain and data. In real applications the \gls{SNR} is strongly determined by the biomedical tools used to observe the data; when using data processed from a scan of the aorta, it can be as low as 10 for flow measures (personal communication of Dr. M. Brummer, Emory CHOA). In Figure \ref{F:carotid} (left) a comparison between reference velocity and noisy data is reported. Data on $\Gamma_{in}$ do not satisfy conditions for the well-posedness, consequently interpolation is performed. In this test case, since we are not provided with real data and an analytic solution is not available, we consider synthetic data generated from a reference solution, $\UU_{\rm FE}$. This is a \gls{FE} solution computed on very fine grid (300\% of nodes used for the test case) where the discretization error is considered fairly small. For the computation of $\UU_{\rm FE}$ we prescribe the following boundary conditions: 
\begin{equation}\label{E:bc_aorta}
\uu=\mathbf{0} \;\;\hbox{on}\,\Gamma_{wall},\;
-\nu \,(\gr\uu+\gr\uu^\TT) \mathbf{n}+ p\,\mathbf{n}=\mathbf{0} \;\; \hbox{on}\,\Gamma_{out},\;
\uu=[0\;\;g(x)]^{\TT} \;\; \hbox{on}\,\Gamma_{in},
\end{equation}
where $g(x)=k\,x\,(2r-x)$, $\nu=3.5\,\frac{mm^2}{s}$, $2r=7\,mm$ ($\Gamma_{in}$ length) and $k=12\,ms^{-1}$.

In the generation process, we introduce a projection matrix, ${\rm P}$ from the fine to the coarse grid. Data are generated adding random noise to P$\UU_{\rm FE}$ on selected sites. The reference solution is affected both by the \gls{FE} approximation error and the projection one; however, this fact does not penalize the reliability of results since approximation and projection errors can be regarded as part of the noise. As index of accuracy we consider the ratio { \small $E_{\UU}=\frac{\|\UU-{\rm P}\UU_{\rm FE}\|_2}{\|{\rm P}\UU_{\rm FE}\|_2}$}.

In Figure \ref{F:carotid} (center) we report the magnitude of velocity field assimilated from noisy data with \gls{SNR} = 20. In this case, with Tikhonov regularization ($\alpha$ = 0.1), $E_{\UU}$ = 8.074e-2 after 60 Picard-Newton iterations. We note that the noise mainly affects the components of the velocity which are transversal to the flow; despite this, the recovered field is close to the reference solution. This fact confirms the noise-filtering property of the proposed method. It is interesting to observe the streamlines in the bulb area (see Figure \ref{F:carotid}, right) where we expect to have recirculation and vortex formation.

To test the competitiveness of the \gls{DA} procedure we compare the relative error of the assimilated velocity with the one of the velocity obtained from a forward simulation (here, Dirichlet boundary conditions are prescribed using the interpolated data retrieved from sparse sites on $\Gamma_{in}$) in this case $E_{\UU}$ = 15.957e-2. This comparison confirms the reliability of the \gls{DA} procedure and the gain, in accuracy, that we obtain with respect to solving a forward problem, even if computationally expensive.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cccc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.32\textwidth]{dcarotid_comp_noise.png} &  
\includegraphics[width=0.35\textwidth]{uMagn_carotid_12.png} &
\includegraphics[width=0.27\textwidth]{stream_lines_carotid_12.png} 
\end{tabular}
\end{center}
\caption{Left: actual data used in the DA process (original, in black, and interpolated). Center: magnitude of the velocity vector field. Right: detail of the stream lines.}
\label{F:carotid}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{The wall shear stress} We consider the computation of the \gls{WSS} $\boldsymbol{\tau}$ which is a quantity of medical interest. An accurate approximation of the \gls{WSS} is fundamental in the investigation of cardiovascular pathologies since it is an index of the possibility of rupture of the vessel wall and formation of stenosis \cite{cardiov}. It is defined as the tangential component of the stress exerted by the fluid (blood in this case) on the wall, formally 
\begin{equation}\label{E:wss}
\boldsymbol{\tau}= \displaystyle\nu(\gr \uu+\gr \uu^\TT) \mathbf{n} - \nu\left( \left(  (\gr \uu +\gr \uu^\TT)\mathbf{n}\right)\cdot \mathbf{n} \right)\mathbf{n},
\end{equation}
where $\mathbf{n}$ is the normal vector in correspondence of the wall. Approximations of the \gls{WSS} retrieved from indirect measurements are in general unreliable because of the post-processing numerical errors and the noise affecting the measures. Including measurements in simulations is a way for improving the reliability of computed solutions. In fact, the introduction of the numerical blood flow model results in noise filtering: thanks to the assimilation process, the recovered \gls{WSS} catches accurately the behavior of the reference one. These preliminary results pinpoint again the role of \gls{DA} as a way for filtering and eventually computing hemodynamics indexes with good accuracy.
In order to have a quantitative comparison with the reference solution we introduce an index of accuracy of the magnitude of the \gls{WSS} vector field on $\Gamma_{wall}$, { \small $E_{\rm WSS}=\frac{\|{\rm WSS}-{\rm WSS}_{\rm FE}\|_2}{\|{\rm WSS}_{\rm FE}\|_2}$}; where ${\rm WSS}_{\rm FE}$ is the value retrieved from the reference solution ${\rm P}\UU_{\rm FE}$. 

Table \ref{T:wss_axisym} (left) reports, in correspondence of decreasing \gls{SNR}, errors obtained with the assimilated velocity field, $E_{\rm WSS,DA}$, compared with those obtained from a forward simulation on the same grid with Dirichlet boundary conditions (interpolating function retrieved from sparse data) on $\Gamma_{in}$, $E_{\rm WSS,FW}$. With high \gls{SNR} the gain obtained with \gls{DA} process is not significant, as we decrease \gls{SNR} we can obtain up to the 50\% of gain with respect to the forward simulation. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[t]
\begin{center}
\begin{tabular}{ |c | c | c | }
    \hline
	\gls{SNR} &  $E_{\rm WSS,DA}$ & $E_{\rm WSS,FW}$\\
	\hline
	100 & 0.2536& 0.2667\\
	\hline
	20 & 0.2591 & 0.3030\\
	\hline
	10  & 0.2738  & 0.3861\\
	\hline
	5   &  0.3149 & 0.6114\\
	\hline
\end{tabular}
\hspace{2.cm}
\begin{tabular}{ | l | l | l | }
    \hline
	slices				&	$S$	&	$E_U$\\
	\hline
	none				&	0	&	0.01600\\
	$s_1$				&	5	&	0.00985\\
	$s_1$, $s_2$		&	5	&	0.00576\\
	$s_1$, $s_2$, $s_3$	&	5	&	0.00992\\
	$s_1$				&	10	&	0.00554\\
	$s_1$, $s_2$		&	10	&	0.00535\\	
	$s_1$, $s_2$, $s_3$	&	10	&	0.00647\\
	$s_4$, $s_5$		&	10	&	0.00724\\		
	\hline
\end{tabular}
\end{center}
\caption{Left: comparison of relative errors for the WSS computed with DA and forward solution for test case I. Right: relative errors for different $s$ and $S$ for test case II.}
\label{T:wss_axisym}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The accuracy in the computation of the \gls{WSS} can be improved; nevertheless, the comparison with results from a forward simulation highlights the strength of the assimilation procedure and the relevance of internal data. Also, numerical tests proved that, as we increase the number of internal data, we obtain more accurate solutions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Test case II}
We consider a 2-dimensional computational domain approximating a section of the aortic arc. We assume to have velocity measurements on $\Gamma_{in}$ and on three internal layers $\Gamma_{d}$, as shown in Figure \ref{F:aorta_noise} (left), with the purpose of mimicking Figure \ref{F:MRIsection}. Data are synthetic and, as in test case I, they are generated from a reference solution, $\UU_{{\rm FE}}$, computed on a very fine grid (300\% of nodes used for the test case). For its computation we prescribe the following boundary conditions
\begin{equation}\label{E:bc_aorta}
\uu=\mathbf{0} \;\hbox{on}\,\Gamma_{wall};\quad 
-\nu(\gr\uu+\gr\uu^{\TT})\mathbf{n} + p\,{\bf n}= 0\; \hbox{on}\,\Gamma_{out};\quad
\uu=[\,0\;\;(x-r_1)(r_2-x)\,]^{T} \; \hbox{on}\,\Gamma_{in},
\end{equation}
where $r1$ = 2 and $r_2$ = 6. In Figure \ref{F:aorta_noise} (right) we report the computed pressure, the velocity vector field and noisy data (\gls{SNR} = 10). The noise mainly affects the components of the velocity which are transversal to the flow; in the recovered field such components are filtered. Here, $E_{\UU}^*= 5.12$e-$2$ after $40$ Picard-Newton iterations. In Figure \ref{F:wss} we report a comparison between the \gls{WSS} computed using the reference solution P$\UU_{\rm FE}$ and the one computed upon assimilated velocities; the assimilated \gls{WSS} approximate the reference \gls{WSS} with $E_{\rm WSS,DA}$ = 0.2401.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.5\linewidth} 
\centering
\includegraphics[scale=0.31]{data_aorta_n5.png}
\end{minipage}
\hspace{0.25cm} 
\begin{minipage}[b]{0.4\linewidth}
\centering
\includegraphics[scale=0.31]{pvd_aorta_n5.png}
\end{minipage}
\caption{On the left, computational mesh reproducing the aortic arc and noise free data. 
On the right, pressure and velocity fields recovered from noise free data.}
\label{F:aorta_noise}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{minipage}[b]{0.5\linewidth} 
\centering
\includegraphics[scale=0.31]{WSSexact.png}
\end{minipage}
\hspace{0.25cm} 
\begin{minipage}[b]{0.4\linewidth}
\centering
\includegraphics[scale=0.37]{WSS.png}
\end{minipage}
\caption{Comparison between the map of the stress exerted by the blood in the arterial wall recovered from the reference solution, on the left, and from noisy velocity data of Figure \ref{F:aorta_noise} (left), on the right.}
\label{F:wss}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Test case III} We consider a 3-dimensional cylindrical domain with radius $r$ and height $h$. Exploiting symmetry, we formulate the problem in an axisymmetric formulation and we solve it problem in a rectangular domain, the shaded area in Figure \ref{F:axisym_domain_velocity} (left). We assume to have data on the inflow boundary not satisfying optimality conditions and we perform interpolation. We also assume that additional internal data are available. More precisely, we consider $s$ internal slices (sections) with $S$ data per slice, with the purpose of reproducing the setting of Figure \ref{F:MRIsection}. Data are generated adding to the Poiseuille flow random noise with \gls{SNR} = 20. In Figure \ref{F:axisym_domain_velocity} (right) we report, for $r$ = 1.5 $cm$, $h$ = 6 $cm$, $s$ = 5 and $S$ = 10, the data and the recovered vector field. The relative error with respect to the Poiseuille reference solution is $E_U$ = 0.007793. In Figure \ref{F:axisym_3d} we report, for $s$ = 5 and $S$ = 17, a 3-dimensional visualization the data and the assimilated velocity: the magnitude of the velocity is displayed in three sections in the lower part of the cylinder, the vector field in the upper ones.

We investigate the accuracy of the assimilated solution using fixed \gls{SNR} and $\Delta$, results are reported in Table \ref{T:wss_axisym} (right). In this test case we consider $s$ slices such that each slice $s_k, \; k=1,\,...,\,s$ is placed at distance $k\,cm$ from $\Gamma_{in}$. The non-monotonic behavior of the error as a function of $s$ might be surprising. The explanation resides in the presence of the noise, which has two different effects on the accuracy of the solution. On one hand additional data help the well-conditioning and accuracy (in absence of noise) of the problem, no matter their location. On the other hand, coherently with the result in \cite{Veneziani20071685}, the incidence of noise in distal sections is supposed to be more relevant than in proximal sections. As a confirmation of this fact, we observe that using $s_1$ and $s_2$ gives a more accurate solution than $s_4$ and $s_5$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}%{c{0.45\textwidth}c{0.45\textwidth}}\\
\includegraphics[width=0.35\textwidth]{cylinder.pdf} &  
\includegraphics[width=0.67\textwidth]{axisymm_vel_data.png}  
\end{tabular}
\end{center}
\caption{Left: in the axisymmetric formulation a rectangle is selected to be the computational domain. Right: data and assimilated vector field for $s$ = 5 and $S$ = 10.}
\label{F:axisym_domain_velocity}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.95\textwidth]{axi_3d.png} 
\end{center}
\caption{3-dimensional visualization the data and the assimilated velocity.}
\label{F:axisym_3d}
\end{figure}



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

