Solving optimal control problems governed by random Navier-Stokes equations using low-rank methods
Peter Benner\footnotemark[1]
Sergey Dolgov\footnotemark[2]
Akwum Onwunta\footnotemark[1]
Martin Stoll\footnotemark[1]
Many problems in computational science and engineering are simultaneously characterized by the following challenging issues: uncertainty, nonlinearity, nonstationarity and high dimensionality. Existing numerical techniques for such models would typically require considerable computational and storage resources. This is the case, for instance, for an optimization problem governed by time-dependent Navier-Stokes equations with uncertain inputs. In particular, the stochastic Galerkin finite element method often leads to a prohibitively high dimensional saddle-point system with tensor product structure. In this paper, we approximate the solution by the low-rank Tensor Train decomposition, and present a numerically efficient algorithm to solve the optimality equations directly in the low-rank representation. We show that the solution of the vorticity minimization problem with a distributed control admits a representation with ranks that depend modestly on model and discretization parameters even for high Reynolds numbers. For lower Reynolds numbers this is also the case for a boundary control. This opens the way for a reduced-order modeling of the stochastic optimal flow control with a moderate cost at all stages.
paper_arxiv.tex
\documentclass[12pt]{article}

\usepackage{amsmath,amssymb,amsthm}
\usepackage{amsfonts}
\usepackage{color}
\usepackage{bm}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{graphicx}
\usepackage[left=2.6cm,right=2.2cm,top=2.5cm,bottom=2.5cm]{geometry}

\newcommand{\be}{\begin{eqnarray}}
\newcommand{\e}{\end{eqnarray}}
\newcommand{\bes}{\begin{eqnarray*}}
\newcommand{\es}{\end{eqnarray*}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\nb}{\nonumber}
%\newcommand{\R}{\ensuremath{\mathbb{R}}}
\newcommand{\trace}[1]{\ensuremath{\mathop{\mathrm{trace}}\left( #1 \right)}}
%\renewcommand{\vec}[1]{\ensuremath{\mathop{\mathrm{vec}}\left( #1 \right)}}

\newcommand{\R}{\ensuremath{\mathbb{R}}}
\newcommand{\K}{\mathcal{K}}
\newcommand{\Am}{\mathcal{A}}
\newcommand{\A}{\mathcal{A}}
\newcommand{\vvo}{{ \mathbf{y}_{d}}}
% \newcommand{\mmu}{\boldsymbol{\mu}}
\newcommand{\La}{L}
\newcommand{\Laa}{\mathcal{L}}
% \newcommand{\M}{\mathcal{M}}
\newcommand{\tM}{\tilde{\mathcal{M}}}
\newcommand{\N}{\mathbf{N}}
\def\L{\mathbf{L}}
\newcommand{\J}{\mathcal{J}}
\newcommand{\h}{\mathcal{H}}
\newcommand{\G}{\mathcal{G}}
\newcommand{\I}{\mathcal{I}}
\newcommand{\p}{\mathcal{P}}

\newcommand{\D}{\mathbf{D}}
\newcommand{\M}{\mathbf{M}}
\newcommand{\Ll}{\mathbf{L}}
\newcommand{\C}{\mathbf{C}}

\newcommand{\Lls}{\mathcal{L}}
\newcommand{\Ms}{\mathcal{M}}
\newcommand{\Ns}{\mathcal{N}_s}
% \newcommand{\Ll}{{L}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\Ab}{\mathbf{A}}

\newcommand{\hA}{\widehat{\mathcal{A}}}
\newcommand{\ltwo}{\ensuremath{L^{2}(\Omega)}}
\newcommand{\ua}{\ensuremath{\underline{u}_{a}}}
\newcommand{\ub}{\ensuremath{\overline{u}_{b}}}
\newcommand{\uopt}{\ensuremath{{u}^{*}}}
\newcommand{\yopt}{\ensuremath{{y}^{*}}}
\newcommand{\lopt}{\ensuremath{{\lambda}^{*}}}
\newcommand{\uh}{\mathbf{\ensuremath{{u}}}}%_{h}}}
\newcommand{\wh}{\mathbf{\ensuremath{{w}}}}%_{h}}}
\newcommand{\W}{\mathbf{\ensuremath{{W}}}}%_{h}}}
\newcommand{\yh}{\mathbf{\ensuremath{{y}}}}
\newcommand{\fh}{\ensuremath{{f}}}
\newcommand{\ybh}{\mathbf{{\bar{y}}}}
\newcommand{\lh}{\boldsymbol{\ensuremath{{p}}}}
\newcommand{\muh}{\boldsymbol{\ensuremath{{\mu}}}}
\newcommand{\etah}{\boldsymbol{\ensuremath{{\eta}}}}
\newcommand{\phk}{\ensuremath{{p}}}%^{(k-1)}}}
\newcommand{\gf}{\ensuremath{\nabla F(\uh^{(k-1)})}}
\newcommand{\inprod}[1]{\langle #1 \rangle}
\newcommand{\ww}{z}
\newcommand{\w}{\vartheta}
\newcommand{\eps}{\varepsilon}
% \newcommand{\m}{\mathbf{m}}
\newcommand{\pvint}{-\hspace{-1em}\int}%ordinary pv int after symbol
\newcommand{\bpvint}{-\hspace{-1em}\int}%pvint after bracket
% \newcommand{\ya}{\ensuremath{\underline{y}_{a}}}
% \newcommand{\yb}{\ensuremath{\overline{y}_{b}}}
\newcommand{\ya}{\ensuremath{\underline{y}}}
\newcommand{\yb}{\ensuremath{\overline{y}}}


\newcommand{\pp}{{\boldsymbol\lambda}}
\newcommand{\xx}{{\bf x}}
\newcommand{\uu}{{\bf u}}
\newcommand{\vv}{{\bf v}}
\newcommand{\mmu}{\boldsymbol{\mu}}

\newcommand{\bg}{{\bf g}}
\newcommand{\bu}{{\bf u}}
\newcommand{\bv}{{\bf v}}
\newcommand{\bz}{{\bf z}}
\newcommand{\bd}{{\bf d}}
\newcommand{\rr}{{\bf r}}
\newcommand{\bb}{{\bf b}}
\newcommand{\bp}{{\bf p}}
\newcommand{\bx}{{\bf x}}
\newcommand{\by}{{\bf y}}

\def\cg{{{\sc Cg}}}
\def\minres{{\scshape Minres}}
\def\symmlq{{ {\sc Symmlq}}}
\def\gmres{{ {\sc Gmres}}}
\def\gmresr{{ {\sc Gmresr}}}
\def\bicgstab{{ {\sc Bicgstab}}}
\def\qmr{{ {\sc Qmr}}}
\def\sqmr{{ {\sc sQmr}}}
\def\bicg{{ {\sc Bicg}}}
\def\ifiss{{ {\sc Ifiss}}}
\def\matlab{{ {\sc Matlab}}}
\def\itfqmr{{ {\sc Itfqmr}}}



\def\abs#1{\left|#1\right|} % absolute value, any size
\def\norm#1{\left\|#1\right\|} % norm, any size




\usepackage{tikz} \usetikzlibrary{positioning,shapes,calc,backgrounds,fit,plotmarks}
\usepackage{pgfplots}
\pgfplotsset{compat=newest,
every axis/.append style={axis x line=bottom,
			  scale only axis,
			  },
}
\tikzset{dashdot/.style={dash pattern=on .4pt off 3pt on 4pt off 3pt}}
\pgfplotscreateplotcyclelist{stokescontrol}{
         black, mark=o, draw, line width=1pt\\%
         blue, dotted, line width=1.6pt, mark=x\\%
         red, dashed, line width=2.4pt, mark=+\\%
         green!50!black, opacity=0.4, line width=3.0pt, mark=triangle*\\%
         orange, dashdot, line width=3.0pt, mark=*\\%
}

% % Draft
% \usepackage{showlabels,rotating}
% \renewcommand{\showlabelsetlabel}[1]
% {\begin{turn}{30}\showlabelfont \footnotesize \textcolor{blue!30!black}{#1}\end{turn}}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}

\date{March 17, 2017}

\begin{document}
% \title{A low-rank tensor approach to random Navier-Stokes optimal control problems}
\title{Solving optimal control problems governed by random Navier-Stokes equations using low-rank methods} % from Peter
% LOW-RANK SOLUTION TO AN OPTIMIZATION PROBLEM CONSTRAINED BY THE NAVIER-STOKES EQUATIONS % determ. paper
\author{Peter Benner\footnotemark[1], Sergey Dolgov\footnotemark[2], Akwum Onwunta\footnotemark[1], Martin Stoll\footnotemark[1]}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\footnotetext[2]{University of Bath,
Department of Mathematical Sciences,
Claverton Down,
BA2 7AY Bath, United Kingdom (\texttt{S.Dolgov@bath.ac.uk}). S. Dolgov gratefully acknowledges funding by the EPSRC Fellowship EP/M019004/1.}
\footnotetext[1]{Computational Methods in Systems and Control Theory, Max Planck Institute
for Dynamics of Complex Technical Systems,
Sandtorstr. 1,
39106 Magdeburg,
Germany. \textit{Email addresses:} \texttt{\{benner,onwunta,stollm\}@mpi-magdeburg.mpg.de}}

\renewcommand{\thefootnote}{\arabic{footnote}}

\maketitle

% REQUIRED
\begin{abstract}
Many  problems in computational science and engineering are simultaneously characterized by
the following  challenging issues:
uncertainty, nonlinearity, nonstationarity and high dimensionality.
Existing numerical  techniques for  such models
would typically require considerable computational and storage resources.
This is the case, for instance,
for an optimization  problem governed by  time-dependent Navier-Stokes equations 
with uncertain inputs. In particular, the stochastic Galerkin finite element
method often leads to a prohibitively high dimensional saddle-point
system with tensor product structure.  In this paper, we approximate the solution by the low-rank Tensor Train decomposition,
and present a numerically efficient algorithm to solve the optimality equations directly in the low-rank representation.
We show that the solution of the vorticity minimization problem with a distributed control admits a representation with ranks 
that depend modestly on model and discretization parameters even for high Reynolds numbers.
For lower Reynolds numbers this is also the case for a boundary control.
This opens the way for a reduced-order modeling of the stochastic optimal flow control with a moderate cost at all stages.
\end{abstract}


% REQUIRED
{\textit{Keywords: }
Stochastic Galerkin system, iterative methods, PDE-constrained optimization, saddle-point system, tensor train format, low-rank solution,  preconditioning, Schur complement. \\
}

% % REQUIRED
% \begin{AMS}
%   68Q25, 68R10, 68U05
% \end{AMS}
\section{Introduction}
We consider the numerical simulation of
optimization problems constrained by  partial differential equations (PDEs). 
This class of problems can be computationally challenging.
This is particularly so if the constraints are time-dependent PDEs since time-stepping methods quickly reach their
limitations due to the enormous demand for storage \cite{StoB13}.
The computational complexity associated with these problems is further increased when the constraints are nonlinear  \cite{ChoLu2015} 
and contain (up to countably many) parametric or uncertain inputs \cite{BSOS215,BOS15a}.
This is the core problem of this paper.
A viable solution approach to optimization problems with stochastic constraints  employs the spectral stochastic Galerkin 
finite element method (SGFEM). However, this intrusive approach\footnote{Generally speaking, SGFEM techniques are intrusive in the
sense that the codes for the associated deterministic problems
cannot be directly reused. These methods are mainly non-ensemble-based methods and require the solution of discrete systems that couple all spatial and 
probabilistic degrees of freedom.} leads  to the so-called \textit{curse of dimensionality}, in the sense
that it results in prohibitively high dimensional linear systems with tensor product structure \cite{BSOS215,BOS14,EW12}.

It is worth pursuing  computationally efficient ways to simulate optimization problems with stochastic constraints 
using SGFEMs since the Galerkin approximation yields the best approximation
with respect to the energy norm, as well as a favorable framework for error estimation \cite{BPS13}.
In order to cope with the curse of dimensionality we exploit  the underlying mathematical structure of the discretized optimality system.
We develop a low-rank technique based on recent advances in numerical
tensor methods  \cite{hackbusch-2012,bokh-surv-2014} for efficient solution of an optimization  problem governed by nonlinear PDEs
with random coefficients. More specifically, we numerically  simulate an  \underline{O}ptimization  \underline{P}roblem constrained
by  time-dependent \underline{N}avier-\underline{S}tokes equations
with random coefficients\footnote{We remark here that our approach can be easily generalized 
to other nonlinear optimization problems with stochastic constraints as well.} (\underline{OPNS}).
Our aim in this paper is to lift the curse of dimensionality inherent in the
OPNS and allow for efficient simulations of the  model on not much more than an average desktop computer.
Such simulations would enhance
the understanding of the underlying physical model as the computed data can then be used for the quantification of the statistics
of the system response. % -- a task that could  in turn result in additional enormous computational expense.


Alternative  approaches to tackle optimization problems with stochastic constraints  include  stochastic collocation schemes \cite{BW09, KHRW13, TKXP12}, 
as well as Monte Carlo methods \cite{AHU16}. These methods are essentially sampling-based and non-intrusive.  
However, for optimization problems, the SGFEM exhibits superior performance compared to the stochastic collocation method \cite{EW12}.
This is due to the fact that, unlike SGFEM,  the non-intrusivity property of the stochastic collocation method
is lost when moments of the state variable appear in the cost functional, or when the control is a
deterministic function. On the other hand, Monte Carlo methods are relatively straightforward to implement.
However, they generally converge rather slowly and do not exploit the  regularity with respect to the parameters that the
solution might have \cite{TKXP12}.


The rest of the paper is organised as  follows. First, we present in Section \ref{probst} the PDE-constrained
optimization that we would like to solve, as well as the necessary mathematical concepts and notation
on which we shall rely in the rest of our discussion. Next, we proceed to Section \ref{discrete_prob} to discuss
the SGFEM discretization of the problem. Section \ref{tt_decom} presents our low-rank iterative
solver and preconditioner which we use to tackle the high-dimensional saddle point systems arising from
the SGFEM discretization of the optimization problem. Finally, in Section \ref{experiments}, we present
numerical experiments to illustrate the performance of the low-rank approach.


\section{Problem statement and mathematical description}
\label{probst}
Let $(\Omega,\mathcal{F},{\mathbb{P}})$ be a complete probability space, where
 $\Omega$ is the set of outcomes, $\mathcal{F}\subseteq 2^{\Omega}$ is the $\sigma$-algebra of events, and $\mathbb{P}:\mathcal{F}\rightarrow [0,1]$ is an appropriate probability measure.
Let $\mathcal{D}\subset \mathbb{R}^d$ with $d\in\{1,2,3\},$ be a bounded spatial domain with a 
 piecewise Lipschitz boundary $\Upsilon.$
Moreover, we consider a finite time interval $[0,T]$.
A random field $z: \Omega \times \mathcal{D} \rightarrow \mathbb{R},$
 means that $z(\cdot,{\bf x})$ is a random variable defined on  $(\Omega,\mathcal{F},{\mathbb{P}})$ for each  ${\bf x}\in \mathcal{D}.$
We assume that $z$ belongs to the  tensor product Hilbert space $L^2(\Omega) \otimes L^2(\mathcal{D})$ which is endowed with the norm
\[
||z||_{L^2(\Omega) \otimes L^2(\mathcal{D})}:=\left(\int_\Omega ||z(\omega,\cdot)||^2_{L^2(\mathcal{D})} \;d\mathbb{P}(\omega)\right)^{\frac{1}{2}} <\infty,
\]
where ${L}^2(\Omega):={L}^2(\Omega,\mathcal{F},{\mathbb{P}}).$
 For any random variable $g$ defined on  $(\Omega,\mathcal{F},{\mathbb{P}}),$
 the mean $\mathbb{E}(g)$ of $g$ is given by
\be
\label{mean}
 \left<g\right>:=\mathbb{E}(g)=\int_\Omega g\;d\mathbb{P}(\omega)<\infty.
\e
For a Hilbert space $H$ of functions on $\mathcal{D}$ and a time interval $[0,T],$ we write  $L^2(0,T;H)$ for  the tensor product space $L^2([0,T])\otimes H.$
However, in particular, we  write $L^2(0,T;\mathcal{D})$ for $L^2(0,T;L^2(\mathcal{D}));$ we will be using these last two notations interchangeably.

% where  $\mathbb{E}(g)$ is the mean of $g$ defined by
% \be
% \label{mean}
% \left<g\right>:=\mathbb{E}(g)=\int_\Omega g\;d\mathbb{P}(\omega)<\infty.
% \e


We consider an optimization problem involving an incompressible Newtonian flow  with an uncertain inflow in a
backward step domain.
% While our approach is general,  we here focus on the numerical simulation of optimization problems
% constrained by Navier-Stokes equations with uncertain coefficients. 
More precisely, the OPNS we want to solve is given by the minimization of the total vorticity: % and the misfit with some desired state:
\begin{align}
\label{J}
% \mathcal{J} & = \frac{1-\alpha}{2}\|\nabla\times\mathbf{v}\|^2
% + \frac{\alpha}{2}\|{\bf v}- \hat{{\bf v}}\|^2 %+ \frac{\alpha}{2}||\mbox{std}({\bf v})||^2_{L^2(0,T;\mathcal{D})}
% +\frac{\beta}{2}\|\mathbf{u}\|^2,
\mathcal{J} & = \frac{1}{2}\|\nabla\times\mathbf{v}\|^2
+\frac{\beta}{2}\|\mathbf{u}\|^2,
\end{align}
where
% $\alpha\in\{0,1\}$ switches between the functionals,
${\bf v},\; {\bf u}: [0,T] \times \Omega \times \mathcal{D}\rightarrow \mathbb{R}^d$ are random fields,
representing the state (velocity) and  the control functions.
%and all norms $\|\cdot\|$ are taken in $L^2([0,T])\otimes  L^2(\Omega) \otimes L^2(\mathcal{D})$.
The regularization constant $\beta$ in \eqref{J} balances between minimization of the vorticity and penalization of the control magnitude.
The objective function $\mathcal{J}({\bf v},{\bf u}) $ is a deterministic quantity with uncertain components ${\bf v}$  and ${\bf u}$.

\begin{figure}[h!]
  \begin{center}
\caption{An uncertain inflow in a backward step domain}
\label{NS_flow}
\resizebox{\linewidth}{!}{
% \begin{tikzpicture}
%  \draw[ultra thick,brown] (0,-1) -- (0,1) node[left,midway] {$\Upsilon_{in}$}; % left
%  \draw[ultra thick] (0,1) -- (8,1) node[above,midway] {$\Upsilon_{wall}$}; % top
%  \draw[ultra thick] (0,-1) -- (8,-1) node[below,midway] {$\Upsilon_{wall}$}; % bottom
%  \draw[ultra thick,dashed] (8,-1) -- (8,1) node[right,midway] {$\Upsilon_{out}$}; % right
%
%  \draw[fill=blue] (1.75,-0.25) rectangle (2.25,0.25);
%  \draw[ultra thick] (1.75,-0.25) -- (1.75,0.25); % node[left,midway] {\small $\partial\mathcal{D}_{wall}$}; % left obs
%  \draw[ultra thick,red] (1.75,0.25) --  (2.25,0.25) node[above,midway] {\small $\Upsilon_{c}$}; % top obs
%  \draw[ultra thick] (2.25,0.25) --  (2.25,-0.25) node[right,midway] {\small $\Upsilon_{wall}$}; % right obs
%  \draw[ultra thick,red] (2.25,-0.25) -- (1.75,-0.25) node[below,midway] {\small $\Upsilon_{c}$}; % bottom obs
%
%  % Inflow arrows g=(1-y)(1+y)
%  \draw[->,brown] (0,-0.75) -- (4.3750e-01,-0.75);
%  \draw[->,brown] (0,-0.50) -- (7.5000e-01,-0.5);
%  \draw[->,brown] (0,-0.25) -- (9.3750e-01,-0.25);
%  \draw[->,brown] (0,0.00) --  (1.0000e+00,0.00);
%  \draw[->,brown] (0,0.25) --  (9.3750e-01,0.25);
%  \draw[->,brown] (0,0.50) --  (7.5000e-01,0.50);
%  \draw[->,brown] (0,0.75) --  (4.3750e-01,0.75);
%
%  % Corners
% \node[anchor=north] at (0,-1.02) {$(0,-1)$};
% \node[anchor=south] at (8,1.02) {$(8,1)$};
% \end{tikzpicture}
%
\begin{tikzpicture}{\small
 \draw[ultra thick,brown] (-1,0) -- (-1,1) node[left,midway] {$\Upsilon_{in}$}; % left
 \draw[ultra thick] (-1,1) -- (5,1) node[above,midway] {$\Upsilon_{wall}$}; % top
 \draw[ultra thick] (-1,0) -- (0,0) node[below,near start] {$\Upsilon_{wall}$}; % top of the step
 \draw[ultra thick,red,dashed] (0,0) -- (0,-1) node[right,midway] {$\Upsilon_{cnt}$};  % left of the step
 \draw[ultra thick] (0,-1) -- (5,-1) node[below,midway] {$\Upsilon_{wall}$}; % bottom
 \draw[ultra thick,dashed] (5,-1) -- (5,1) node[right,midway] {$\Upsilon_{out}$}; % right

 % Inflow arrows g=4 y (1-y)
 \draw[->,brown] (-1,0.125) -- (-0.5625,0.125);
 \draw[->,brown] (-1,0.250) -- (-0.2500,0.250);
 \draw[->,brown] (-1,0.375) -- (-0.0625,0.375);
 \draw[->,brown] (-1,0.500) -- ( 0.0000,0.500);
 \draw[->,brown] (-1,0.625) -- (-0.0625,0.625);
 \draw[->,brown] (-1,0.750) -- (-0.2500,0.750);
 \draw[->,brown] (-1,0.875) -- (-0.5625,0.875);

 % Corners
\node[anchor=south] at (-1,1.02) {$(-1,1)$};
\node[anchor=south] at (5,1.02) {$(5,1)$};
\node[anchor=north] at (0,-1.02) {$(0,-1)$};
\node[anchor=west] at (-0.05,0.05) {$(0,0)$};
}
\end{tikzpicture}
}
\end{center}
\end{figure}

The minimization of $\mathcal{J}$ is considered subject, $\mathbb{P}$-almost surely, 
to the Navier-Stokes equations\footnote{In this paper, we do not consider the
case of state- or control- or mixed control-state-constrained problems \cite{HIK02,SPW11,PST15}. 
These problems can be tackled via, for instance, semi-smooth Newton algorithms \cite{HS10, IK03, Kanzow04}.}
\be%
 \label{modelprob}
\left\{
\begin{aligned}
    \frac{\partial {\bf v}}{\partial  t}  - \nu \Delta {\bf v} + ({\bf v}\cdot \nabla) {\bf v} + \nabla p &= F_d {\bf u},
       & \mbox{in} \quad  (0,T]\times \Omega  \times  \mathcal{D},  \\
    -\nabla \cdot {\bf v} &= 0, & \mbox{in} \quad (0,T]\times \Omega  \times\mathcal{D} , \\
     {\bf v} &= {\boldsymbol \theta}, & \mbox{on} \quad  (0,T]\times  \Omega \times \Upsilon_{in},\\
     {\bf v} &= 0, & \mbox{on} \quad  (0,T]\times \Omega \times \Upsilon_{wall} ,\\
     \frac{\partial \vv}{\partial  {\bf n}}  &= 0,  & \mbox{on} \quad (0,T]\times \Omega  \times  \Upsilon_{out},\\
     \frac{\partial \vv}{\partial  {\bf n}}  &= F_b {\bf u},  & \mbox{on} \quad (0,T]\times \Omega  \times  \Upsilon_{cnt},\\
     {\bf v}(0,\cdot,\cdot) &= 0, & \mbox{in} \quad \Omega   \times \mathcal{D}.
\end{aligned}
\right.
\e
Here,
${\boldsymbol \theta}:  [0,T] \times \Omega \times \Upsilon_{in} \rightarrow \mathbb{R}^d$ is the inflow boundary condition,
$p: [0,T] \times \Omega \times \mathcal{D} \rightarrow \mathbb{R}$ is the pressure, and
${\bf n}$ is the outward-pointing  normal to the boundary.
The linear  operators $F_d$ and $F_b$ apply the control to the appropriate part of the system.
We consider two cases.
\begin{enumerate}
 \item \emph{Distributed} (or full) control: $F_d = I$, $F_b = 0$, and
 \item \emph{Boundary} control: $F_d=0$, and $F_b = I_{\Upsilon_{cnt}}$, the identity on $\Upsilon_{cnt}$ and zero otherwise.
\end{enumerate}
%Note that the condition  $ \frac{\partial \vv}{\partial  {\bf n}}  = 0$  in  (\ref{modelprob}) is particularly important when the $\Upsilon_{out}$  is very far from the location of the obstacle in the channel.
%  As shown in Figure \ref{NS_flow}, we want to manipulate the inflow at the obstacle's boundary $\Upsilon_{{c}}$
%  in such a way that the flow vorticity is minimized.
% Applying a control ${\bf u}$ to a part of the boundary means imposing Neumann boundary conditions, and this gives us a Neumann control problem in a natural way.

The  parameter  $\nu$ is the  kinematic  viscosity.
There exist variations of this problem such as a deterministic control \cite{EW12} or an uncertain domain \cite{BH14}.
% and is often modelled deterministically.
% However, its value is not always known precisely in practice.
% Rather than guessing a value it is more sensible to  model  $\nu$ as a random variable defined on $(\Omega,\mathcal{F},{\mathbb{P}}).$ More specifically, we model the kinematic viscosity $\nu$ in the state equations  as
% \be
% \label{visc}
% \nu(\omega) = \nu_0 + \nu_1\xi(\omega),  \quad\;  \nu_0, \nu_1\in\mathbb{R}^+, \quad \nu_1<\nu_0,
% \e
% where  $\xi $ is a uniformly distributed random variable with  $\xi \sim \mathcal{U}(-1,1)$.
% As a result, we have for some $0<\nu_{\min}< \nu_{\max}<+\infty$ that
% \be
% \label{probcond}
% \mathbb{P}\left( \omega\in \Omega:\nu(\omega)\in [\nu_{\min}, \nu_{\max}]  \right)=1.
% \e
% This could be interpreted as a scenario where the volume of fluid moving into a channel is uncertain due to measurement error in $\nu$ or probably some other factors  \cite{PS12}.
 
In our experiments, we use a zero initial condition, applying instead an exponentially growing stochastic inflow
\[
\boldsymbol\theta = \begin{bmatrix}\kappa({\bf x},\omega)(1-\mathrm{e}^{-t}) \\ 0\end{bmatrix}
\]
 (see Fig. \ref{NS_flow}), where $\kappa({\bf x},\omega)$ is a random field on $\Upsilon_{in}$.
We assume that $\kappa$ admits the  Karhunen-Loeve expansion (KLE):
\begin{equation}
\kappa({\bf x},\omega) = 4 x_2 (1-x_2) + \sum_{k=1}^{\infty} \frac{1}{2} \gamma^{-k} \cdot \sin(2 \pi k x_2) \cdot \xi_k(\omega).
\label{eq:kle_in}
\end{equation}
Here, the first term is the mean parabolic inflow, whereas $\xi_k$, $k=1,2,\ldots$ are independent uniformly distributed random quantities, $\xi_k \sim \mathcal{U}(-1,1)$.
The parameter $\gamma$ governs the KLE decay rate, and hence the smoothness of the field.
In computational practice, we truncate (\ref{eq:kle_in}) after $m\in \mathbb{N}$ terms such that the error is
sufficiently small:
\begin{equation}
\kappa({\bf x},\omega) \approx \kappa_m({\bf x},\xi(\omega))=  4 x_2 (1-x_2) + \sum_{k=1}^{m} \frac{1}{2} \gamma^{-k} \cdot \sin(2 \pi k x_2) \cdot \xi_k(\omega),
\label{eq:kle_in_new}
\end{equation}
where $\xi:=\{\xi_1,\xi_2,\ldots,\xi_m\}.$
This coefficient mimics the (truncated) Fourier expansion of the Matern covariance function for $\kappa.$
The assumption that  $\kappa(x,\omega)$ admits a KLE representation allows to transform
the stochastic OPNS into a parametric deterministic problem, depending on $\xi:=\{\xi_1,\xi_2,\ldots,\xi_m\}.$
For brevity, since we consider $m$ as an a priory model parameter, we always denote $\kappa_m$ by $\kappa$ in the rest of the paper.

We also assume that each random variable is  characterized by a probability density function $\rho_i:\Gamma_i\rightarrow[0,1].$
If the distribution measure of the random vector $\xi(\omega)$ is absolutely continuous with respect
to the Lebesgue measure, then there exists a joint probability density function
$\rho:\Gamma\rightarrow \mathbb{R}^+,$ where $\Gamma:=\bigotimes^m_{k=1}\Gamma_k\subset \mathbb{R}^m,  \;\;  \rho(\xi)=\prod^m_{k=1}\rho_k(\xi_k), $  and $\rho\in L^{\infty}(\Gamma).$
In particular, given the parametric representation (\ref{eq:kle_in_new}) of  $\kappa({\bf x},\omega),$
 the Doob-Dynkin Lemma \cite{BTZ04}, guarantees that the  solution, for example, $\vv$,
corresponding to the the OPNS  as given by (\ref{J})  and (\ref{modelprob})  admits  exactly the same parametrization;
that is,
\be
\label{solnv}
 \vv(t,\omega,{\bf x})=\vv(t,\xi_1(\omega),\xi_2(\omega),\ldots, \xi_m(\omega),{\bf x}).
\e
Furthermore, we can now replace the probability space $(\Omega,\mathcal{F},{\mathbb{P}})$ with $(\Omega,\mathbb{B}(\Gamma),\rho(\xi)d\xi),$ where $\mathbb{B}(\Gamma)$ denotes the
Borel $\sigma$-algebra on $\Gamma$ and $\rho(\xi)d\xi$ is the finite measure of the vector $\xi.$
Besides, denoting  the space of square-integrable random variables with respect to the density $\rho$ by ${L}_\rho^2(\Gamma),$
 we introduce the space  $L_\rho^2(\Gamma) \otimes L^2(\mathcal{D})$ equipped with the norm
\be
\label{hilrho}
||z||_{L_\rho^2(\Gamma) \otimes L^2(\mathcal{D})}:=\left(\int_\Gamma ||z(\xi,\cdot)||^2_{L^2(\mathcal{D})}\rho(\xi) \;d\xi\right)^{\frac{1}{2}} <\infty.
\e
Similarly, using equation (\ref{mean}), we have
\be
\label{exp2}
 \left<g\right>:= \mathbb{E}[g]= \int_\Gamma g(\xi) \rho(\xi) \;d\xi <\infty.
\e
% More precisely, the parameterized OPNS now reads:
% \begin{align}
% \label{J_new}
% \mathcal{J} & = \frac{1}{2}\|\nabla\times\mathbf{v}\|^2_{L^2(0,T;\mathcal{D})\otimes L^2(\Omega)}
% %\mathcal{J} &= \frac{1}{2}||{\bf v}- \bar{{\bf v}}||^2_{L^2(0,T;\mathcal{D})\otimes L^2(\Omega)} %+ \frac{\alpha}{2}||\mbox{std}({\bf v})||^2_{L^2(0,T;\mathcal{D})}
% +\frac{\beta}{2}\|\mathbf{u}\|^2_{L^2(0,T;\mathcal{D})\otimes L^2(\Omega)},
% \end{align}
% subject, almost surely to
% \be%
%  \label{modelprob_new}
% \left\{
% \begin{aligned}
%     \frac{\partial {\bf v}}{\partial  t}  - \nu \Delta {\bf v} + ({\bf v}\cdot \nabla) {\bf v} + \nabla p &= u,
%        \quad \mbox{in} \quad   (0,T]\times  \Omega \times \mathcal{D},  \\
%    -\nabla \cdot {\bf v} &= 0, \;\quad \mbox{in} \quad (0,T]\times \Omega \times\mathcal{D}, \\
%      {\bf v} &= {\boldsymbol \theta}, \quad \mbox{on} \quad  (0,T]\times  \Omega \times\Gamma_{in},\\
%      {\bf v} &= 0, \quad \mbox{on} \quad  (0,T]\times  \Omega \times\Gamma_{wall},\\
%      \frac{\partial \vv}{\partial  {\bf n}}  &= 0, \quad \mbox{on} \quad (0,T]\times  \Omega \times \Gamma_{out},\\
%      {\bf v}(\cdot,0,\cdot) &= {\bf v}_0, \;\; \mbox{in} \quad   \Omega \times\mathcal{D},
% \end{aligned}
% \right.
% \e
% 

For the  numerical simulation of the  OPNS given by \eqref{J} and \eqref{modelprob}, 
we will adopt the so-called $optimize$-$then$-$discretize$ (OTD) strategy
in which case  we first build an infinite dimensional Lagrangian and then 
consider its variation with respect to state, pressure, control, and 
two Lagrange multipliers that can be identified as the adjoint state
${{\boldsymbol{\lambda}}}$ and adjoint pressure $\mu$  \cite{book::FT2010,Hinze2000}. 
To this end, observe now that the optimization problem under consideration is nonlinear due to the nonlinearity
of the convective term $({\bf v}\cdot \nabla) {\bf v}$.
%The use of mixed finite element approximations in space leads to a nonlinear algebraic systems of equations.
Both Newton and Picard  iterations have shown
to be good iterative solvers  to tackle these nonlinear equations \cite{ESW14}. 
Since the Picard iteration has a larger radius of convergence compared to the Newton iteration,
our  choice in this contribution is the Picard iteration. We apply the
  so-called Karush-Kuhn-Tucker procedure \cite[Chapter 8.2]{ESW14} 
 to  obtain the optimality system \cite{SSNS17}:
\begin{alignat}{3}
\label{kkt_cond1}
\partial_t\vv-\nu \Delta \vv+ \left(\bar\vv \cdot \nabla\right)\vv+\bigtriangledown p &= F_d{\bf u}  &&\nonumber\\
\nabla\cdot \vv &= 0  && \nonumber\\
\mbox{on} \quad \Upsilon_{wall}: \qquad \vv &= 0  \nonumber\\\
\mbox{on} \quad \Upsilon_{in}:   \qquad {\bf v} & = {\boldsymbol \theta} \nonumber\\
\mbox{on} \quad \Upsilon_{out}: \qquad \frac{\partial \vv}{\partial  {\bf n}}  &= 0, \\
\mbox{on} \quad \Upsilon_{cnt}: \qquad  \frac{\partial \vv}{\partial  {\bf n}}  &= F_b {\bf u}, \nonumber\\
     {\bf v}(0,\cdot,\cdot) &= 0,
\nonumber&\\
\nonumber&\\
\label{kkt_cond2}
-\partial_t\pp- \nu \Delta \pp- \left(\bar\vv \cdot \nabla\right)\pp+\left(\nabla\bar\vv\right)^{T}\pp+\bigtriangledown \mu&=
-\nabla \times (\nabla \times {\bf v}) &&\nonumber\\
% -(1-\alpha)\nabla \times (\nabla \times {\bf v}) - \alpha (\vv - \hat\vv)  &&\nonumber\\
\nabla\cdot \pp&=0  &&\nonumber\\
\mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{in}: \qquad \pp&=0  \\
\mbox{on} \quad \Upsilon_{out} \cup \Upsilon_{cnt}: \qquad \frac{\partial\pp}{\partial\mathbf{n}} & =0 \nonumber\\
\pp(T,\cdot,\cdot) & = 0, \nonumber\\
\nonumber&\\
\label{kkt_cond3}
\beta \uu+(F_d+F_b)\pp &=0,  &&
\end{alignat}
where 
%for brevity we omitted initial/final and boundary conditions. 
%Here, 
${\bar \vv}$ denotes the velocity from the previous Picard iteration.
Having solved this system, we update ${\bar\vv} = {\vv}$ and so on until convergence.

\section{The discrete problem}
\label{discrete_prob}
The stochasticity of the inflow $\kappa({\bf x},\omega)$ is inherited by all solution components.
For example, the velocity becomes a function of time, $d$ spatial variables and $m$ random parameters, $\vv(t,\xi,{\bf x})$.
This makes the problem far more challenging than its deterministic counterpart.
In this paper, we employ simultaneous independent discretization of all variables, and the stochastic Galerkin method (SGFEM) w.r.t. the parameters $\xi$.

% In this paper, we assume that the state, the control, and the Lagrange multiplier are analytic functions depending on
% the uncertain parameters.
% Now, the solution (that is, velocity and pressure) of
% automatically inherits that stochasticity
% in the coefficient $\kappa({\bf x},\omega)$  and, as a consequence, the numerical solution of the OPNS is
% far more challenging than
% its deterministic counterpart.
% %Each variable $\xi_k$ is treated as an independent coordinate.
% Based on the fact that we can express the solution  in the form (\ref{solnv}),
% %  \[
% %   {\bf y}(t, \xi, x) := y(t,\xi_1,\ldots,\xi_M,x),
% %  \]
%  it is natural to
% treat $  {\bf y}(t, \xi, {\bf x}),$ a function of time, $d$ spatial variables and $ M$ random parameters,
% as a function of $1+ d + M$ variables.
% As the result, the whole problem has $1+ d + M$  dimensions, which need to be discretized simultaneously.
% In order to study efficient solvers of our OPNS model problem, we consider the stochastic, temporal, and spatial discretization via the stochastic Galerkin method, which we discuss next.
%
The  SGFEM is an intrusive approach
% \footnote{Generally,
% SGFEM techniques are intrusive in the sense that they are non-ensemble-based
% methods; that is, they require the solution of discrete systems that couple
% all spatial and parametric degrees of freedom.}
 in which  one seeks the solution in a finite-dimensional subspace
 $   \mathcal{Y}_{n_\xi} \otimes  \mathcal{X}_h \subset L_\rho^2(\Gamma) \otimes L^2(0,T;L^2(\mathcal{D}))$
consisting of tensor products of deterministic functions defined on the spatial domain
and stochastic functions defined on the probability space
\cite{GWZ14,PE09}. 
Different classes of SGFEMs are distinguished by their choices for $\mathcal{Y}_{n_\xi}.$
When all random variables $\xi_k$ are independent
and identically distributed Gaussian, the basis of multidimensional
Hermite polynomials of the total degree $n_\xi$ is called the {\it polynomial chaos}, a terminology originally introduced by Norbert
Wiener \cite{NW38} in the context of turbulence modeling.
The use of Hermite polynomials ensures that the corresponding basis functions are orthogonal with
respect to the Gaussian probability measure. This leads to sparse linear systems, a crucial property that
must be exploited for fast solution schemes \cite{PE09}.
Relying on the  fact that there exists a one-to-one correspondence between the probability density functions
of alternative distributions and the weight functions of certain orthogonal polynomials,
 the concept of Hermite polynomials chaos has  been extended to {\it generalized} polynomial chaos \cite{XK03}.
For instance, if uniform
random variables (having support on a bounded interval)
are chosen, then Legendre polynomials are the correct choice. Similarly, Jacobi polynomials go with beta-distributed random variables.
When random variables with bounded images are used, the convergence and approximation properties
of the resulting SGFEM are discussed in \cite{BTZ04}.
Some SGFEM classes use tensor products of piecewise
polynomials on the subdomains $\Gamma_k\subset \Gamma$ \cite{BTZ04, EEOS05, NT09}, where
the polynomial degree is fixed and approximation is improved by refining the partition of $\Gamma.$


Throughout this paper we use the classical, so-called {\it spectral} SGFEM (see e.g. \cite{EF07, GS91,  LKDNG03, PE09}), which employs global Legendre polynomials for uniformly distributed $\xi_k$.
We discretize each $\xi_k$ independently via polynomials of the maximal degree $n_\xi-1$.
That is, we use the space
\begin{equation*}
\mathcal{Y}_{n_\xi} = \left\{\mathrm{span} \left(\psi_j(\xi) = \prod_{k=1}^{m} \psi_{j_k} (\xi_k)\right), \quad \forall j_k = 0,\ldots,n_\xi-1, \quad j=(j_1,\ldots,j_m) \right\}.
\end{equation*}
This implies that the global discretization is done via multilinear polynomials of total degree $m (n_\xi-1)$.
The total number of basis functions in $\mathcal{Y}_{n_\xi}$ is $N_\xi=n_\xi^{m}$, which can be prohibitively large for a straightforward storage.
We circumvent this problem by using low-rank product decompositions on the discrete level in Section \ref{tt_decom}.

For the spatial discretization we use  mixed finite element spaces
with stable elements, i.e. elements that satisfy the {\it inf-sup} condition, e.g. the $Q_2$-$Q_1$ pair \cite{JS14}.
This gives us different spaces $V_h=\mbox{span}\{\phi_1,\ldots, \phi_{N_v}\} \subset L^2(0,T;H^1_0(\mathcal{D}))$ 
for the velocity and $\tilde V_h=\mbox{span}\{{\tilde\varphi}_1,\ldots, {\tilde\varphi}_{N_p}\}\subset L^2(0,T;L^2(\mathcal{D}))$ for the pressure.


The time is discretized on a uniform grid, $t_n = \tau n$, $n=1,\ldots,N_t$, such that $t_{N_t}=T$.
The fully discrete coefficients of the velocity will be denoted as $\vv_h(n,j,k)$,
and similarly the coefficients of $p,{\bf u},\pp$, and $\mu$ are denoted by $p_h, {\bf u}_h, \pp_h$ and $\mu_h$, respectively.
The solution functions are approximated as
\begin{equation}
\label{sbdiscrete}
\begin{aligned}
  \vv(t_n,\xi,{\bf x}) & \approx \sum_{j=0}^{N_\xi-1} \sum_{k=1}^{N_v}\vv_h(n,j,k) \psi_j(\xi) \phi_k({\bf x}) , \\
  p(t_n,\xi,{\bf x}) & \approx \sum_{j=0}^{N_\xi-1} \sum_{k=1}^{N_p}p_h(n,j,k) \psi_j(\xi) {\tilde\varphi}_k({\bf x}).
\end{aligned}
\end{equation}
Analogous expressions hold for ${\bf u}, \pp$ and $\mu$.


We plug these expressions into \eqref{kkt_cond1}--\eqref{kkt_cond3}, and project the optimality system onto $\mathrm{span}(\psi_j \phi_k)$ for the velocity equations and $\mathrm{span}(\psi_j \tilde\varphi_k)$ for the pressure equations.
We introduce the following matrices, composed from the corresponding bilinear forms:
\begin{itemize}
 \item $L_0(k,k') = \int\limits_{\mathcal{D}} \nabla \phi_k({\bf x}) \cdot \nabla \phi_{k'}({\bf x}) d{\bf x}$ is the Laplace matrix for a single velocity component, and $L = \mathrm{blkdiag}(L_0, L_0)$ is the Laplace matrix for both velocity components.
 \item $M_0(k,k') = \int\limits_{\mathcal{D}}  \phi_k({\bf x}) \phi_{k'}({\bf x}) d{\bf x}$ and $M = \mathrm{blkdiag}(M_0,M_0)$, the mass matrix.
 \item $B(k,k') = \int\limits_{\mathcal{D}}  \tilde\varphi_k({\bf x}) \nabla \phi_{k'}({\bf x}) d{\bf x}$, the mixed gradient matrix.
 \item $N_{n} = \mathrm{blkdiag}(N_{0,n},N_{0,n})$, the convection matrix, where
 \begin{eqnarray}
 N_{0,n}(jk,j'k') & = & \int\limits_{\Gamma} \psi_j \sum_{j''=0}^{N_\xi-1} N_{s}[\bar\vv_h(n,j'',:)](k,k') \psi_{j''} \psi_{j'} \rho d\xi, \\
 N_{s}[\bar\vv_h](k,k') & = & \int\limits_{\mathcal{D}} \phi_k({\bf x}) \sum_{k''=1}^{N_v} \phi_{k''}({\bf x}) \bar\vv_h(k'') \cdot \nabla \phi_{k'}({\bf x}) d{\bf x}.
 \label{eq:N_detailed}
 \end{eqnarray}
 \item $W_{n}(jk,j'k') = \int\limits_{\Gamma} \psi_j(\xi) \sum_{j''=0}^{N_\xi-1} W_{s}[\bar\vv_h(n,j'',:)](k,k') \psi_{j''}(\xi) \psi_{j'}(\xi) \rho(\xi)d\xi$, the adjoint convection matrix, where
\begin{equation}
 W_{s}[\bar\vv_h](k,k') = \int\limits_{\mathcal{D}} \phi_k({\bf x}) \sum_{k''=1}^{N_v}\nabla\bar\vv_h(k'') \phi_{k''}({\bf x}) \phi_{k'}({\bf x}) d{\bf x}.
\label{eq:Ws}
\end{equation}
\end{itemize}
Since the Legendre polynomials are orthogonal, a proper normalization can turn the stochastic mass matrix into an identity.


We use the implicit Euler approximation for the time derivative, although it can be replaced by higher-order schemes.
Overall, this gives the following discrete optimality equations.
\allowdisplaybreaks
\begin{alignat}{3}
\label{kkt_discr1}
\frac{\vv_h(n,j,k)-\vv_h(n-1,j,k)}{\tau} - \sum_{k'=1}^{N_v} \nu L(k,k')\vv_h(n,j,k') \nonumber\\
 + \sum_{j'=0}^{N_\xi-1}\sum_{k'=1}^{N_v} N_{n}(jk,j'k') \vv_h(n,j',k') + \sum_{k'=1}^{N_p} B(k',k)^T p_h(n,j,k') &= \sum_{k'=1}^{N_v} M_{cnt}(k,k') {\bf u}_h(n,j,k')  &&\nonumber\\
\sum_{k'=1}^{N_v} B(k,k') \vv_h(n,j,k') &= 0  && \\
\mbox{for} \quad k \in \Upsilon_{wall}: \qquad \vv_h(n,j,k) &= 0  \nonumber\\\
\mbox{for} \quad k \in \Upsilon_{in}:   \qquad \vv_h(n,j,k) & = {\boldsymbol \theta}_h(n,j,k) \nonumber\\
     \vv_h(0,j,k) &= 0,
\nonumber&\\
\nonumber&\\
\label{kkt_discr2}
\frac{\pp_h(n,j,k)-\pp_h(n+1,j,k)}{\tau} - \sum_{k'=1}^{N_v} \nu L(k,k')\pp_h(n,j,k')   \nonumber\\
 + \sum_{j'=0}^{N_\xi-1}\sum_{k'=1}^{N_v} (-N_{n}(jk,j'k') + W_{n}(jk,j'k')) \pp_h(n,j',k') \nonumber\\
 + \sum_{k'=1}^{N_p} B(k',k)^T \mu_h(n,j,k') &= \sum_{k'=1}^{N_v}   L(k,k') \vv_h(n,j,k') &\nonumber\\
\sum_{k'=1}^{N_v} B(k,k') \pp_h(n,j,k') &= 0  && \\
\mbox{for} \quad k \in \Upsilon_{wall}\cup \Upsilon_{in}: \qquad \pp_h(n,j,k)&=0  \nonumber\\
\pp_h(N_t,j,k) &= 0, \nonumber\\
\nonumber&\\
\label{kkt_discr3}
\beta \sum_{k'=1}^{N_v} M(k,k') \uu_h(n,j,k')+ \sum_{k'=1}^{N_v} M_{cnt}(k,k') \pp_h(n,j,k') &=0,  &&
\end{alignat}
where $M_{cnt} = F_d + F_b$.
%\textcolor{red}{These $\uparrow\uparrow\uparrow$ $\tau$ should be removed I guess, since either there should be $\tau$ everywhere from the integral of the Lagrange multiplier with the constraints, or we divide all equations by $\tau$.}
% Next, we plug the expressions (\ref{sbdiscrete}) into  (\ref{kkt_cond1}), (\ref{kkt_cond2}) and (\ref{kkt_cond3}) to obtain
% \be
% \label{kkt_cond_new1}
% %\left\{
% \begin{aligned}
% \sum_{i=0}^{N_{\xi}-1}\partial_t v_i\psi_i - \nu \sum_{i=0}^{N_{\xi}-1} \Delta v_i\psi_i
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1}\left(\bar v_i \cdot \nabla\right)v_j\psi_i\psi_j
% + \sum_{j=0}^{N_{\xi}-1}\bigtriangledown p_i\psi_i  &=\sum_{i=0}^{N_{\xi}-1}u_i\psi_i  &&\\
% \sum_{i=0}^{N_{\xi}-1}\nabla\cdot v_i \psi_i&=0  && \nonumber\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{out}: \qquad \sum_{i=0}^{N_{\xi}-1} v_i \psi_i &=0  \\
%  \mbox{on} \quad \Upsilon_{in}: \qquad \sum_{i=0}^{N_{\xi}-1} v_i \psi_i &= {\boldsymbol \theta}_N \\
%      \sum_{i=0}^{N_{\xi}-1} \frac{\partial v_i \psi_i }{\partial  {\bf n}}  &= 0, \\
%     \sum_{i=0}^{N_{\xi}-1} v_i(\cdot,0,\cdot)\psi_i &= 0,
% %    \nonumber&\\
% \end{aligned}
% %\right.
% \e
% %\vskip{1mm}
% \be
% \label{kkt_cond_new2}
% %\left\{
% \begin{aligned}
% -\sum_{i=0}^{N_{\xi}-1}\partial_t \lambda_i \psi_i - \nu \sum_{i=0}^{N_{\xi}-1} \Delta \lambda_i \psi_i
% - \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\bar v_i  \cdot \nabla\right)\lambda_j\psi_i\psi_j &  &&\\
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\nabla\bar v_i \right)^{T}\lambda_j\psi_i\psi_j
% +  \sum_{i=0}^{N_{\xi}-1}\bigtriangledown \mu_i\psi_i &=
% - \sum_{i=0}^{N_{\xi}-1}\nabla \times (\nabla \times v_i\psi_i)  &&\\
% \sum_{i=0}^{N_{\xi}-1}\nabla\cdot \lambda_i \psi_i&=0  &&\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{in}: \qquad \sum_{j=0}^{N_{\xi}-1}\lambda_i \psi_i &=0  \\
% \mbox{on} \quad \Upsilon_{out}: \qquad  \sum_{i=0}^{N_{\xi}-1} \frac{\partial \lambda_i \psi_i }{\partial  {\bf n}} & =0 \\
% \sum_{i=0}^{N_{\xi}-1} \lambda_i(\cdot,T,\cdot)\psi_i & = 0,
% %    \nonumber&\\
%  \end{aligned}
% %\right.
% \e
% and
% \begin{eqnarray}
% \label{kkt_cond_new3}
% \sum_{i=0}^{N_{\xi}-1}\lambda_i\psi_i &=& - \beta \sum_{i=0}^{N_{\xi}-1}u_i\psi_i.
% %\nonumber&\\
% \end{eqnarray}
% Note that in (\ref{kkt_cond_new1}),
% \begin{eqnarray}
% \label{theta_n}
% {\boldsymbol \theta}_M &=& \begin{bmatrix} \kappa_M(x,\xi)(1-\mathrm{e}^{-t}) \\ 0\end{bmatrix}=
%   \begin{bmatrix}\sum_{j=0}^{M}\kappa_j(1-\mathrm{e}^{-t}) \\ 0\end{bmatrix}
% \end{eqnarray}
% Here, $\kappa_0= 4 x_2 (1-x_2)$ and $\kappa_j= \frac{1}{2} \gamma^{-j} \cdot \sin(2 \pi j x_2) \cdot \xi_j,
% \quad j=1,\ldots,M.$
%
%  Now, eliminating $u_i$ from (\ref{kkt_cond_new3}) and performing a Galerkin projection of
%  (\ref{kkt_cond_new1}) and (\ref{kkt_cond_new2}) onto  the space $\mathcal{Y}_n\subset L_\rho^2(\Gamma)$ spanned by the $M$-dimensional
%  Legendre polyomials $\psi_k, \; k=0,1,\ldots N_{\xi}-1,$ we obtain, respectively, coupled linear primal and dual systems
% \textcolor{red}{Why we need to eliminate $u$? It's not how we proceed in practice, moreover, for a partial control the equations below are not true. We can introduce the time discretization immediately, and leave only one final instance of this large formula, \eqref{kktime2}}
%  \be
% \label{kkt_cond_neweq1}
% %\left\{
% \begin{aligned}
% c_{kk}\partial_t v_k - \nu  c_{kk}\Delta  v_k
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1}\left(\bar v_i \cdot \nabla\right) v_j c_{ijk}
% + \bigtriangledown c_{kk} p_k  &= -\frac{1}{\beta}\lambda_k  c_{kk} &&\\
% \nabla\cdot v_k &=0  && \nonumber\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{out}: \qquad  v_k &=0  \\
%  \mbox{on} \quad \Upsilon_{in}: \qquad c_{kk}v_k  &= \tilde{\boldsymbol \theta}_N \\
%      \frac{\partial v_k  }{\partial  {\bf n}}  &= 0, \\
%    v_k(\cdot,0,\cdot) &= 0,
% %    \nonumber&\\
% \end{aligned}
% %\right.
% \e
%  and
% \be
% \label{kkt_cond_neweq2}
% %\left\{
% \begin{aligned}
% -c_{kk}\partial_t \lambda_k - \nu c_{kk}\Delta \lambda_k
% - \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\bar v_i  \cdot \nabla\right)\lambda_j c_{ijk} &  &&\\
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\nabla\bar v_i \right)^{T}\lambda_j c_{ijk}
% +  \bigtriangledown c_{kk}\mu_k &=
% c_{kk}\Delta  v_k &&\\
% \nabla\cdot \lambda_k &=0  &&\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{in}: \qquad \lambda_k  &=0  \\
% \mbox{on} \quad \Upsilon_{out}: \qquad   \frac{\partial \lambda_k  }{\partial  {\bf n}} & =0 \\
%  \lambda_k(\cdot,T,\cdot) & = 0,
% %    \nonumber&\\
%  \end{aligned}
% %\right.
% \e
% where we have used the fact that
% $\nabla \times (\nabla \times v_i) = \nabla (\nabla\cdot v_i) - \nabla^2 v_i = - \Delta v_i,$
% since $v_i$ is divergence-free. Here,
% $c_{kk}:=\left<\psi^2_k(\xi)\right>$  and $c_{ijk}:=\left<\psi_i(\xi) \psi_j(\xi)\psi_k(\xi)\right>.$
% Moreover, from (\ref{theta_n}) we have in (\ref{kkt_cond_neweq1})  that
% \[
% \tilde{\boldsymbol \theta}_M =
%   \begin{bmatrix}\sum_{j=0}^{M} \left<\kappa_j\psi_k(\xi)\right>  (1-\mathrm{e}^{-t}) \\ 0\end{bmatrix}.
% \]
%
% For temporal discretization, we use implicit time stepping techniques to
% semi-discretize the optimality systems given by
% (\ref{kkt_cond_neweq1}) and (\ref{kkt_cond_neweq2})
% to avoid stability issues.
% This allows us to be able to choose the timestep size only depending on the accuracy
% demands, independent of any stability constraints. Here, we  specifically
% use the backward Euler scheme as a representative of implicit
% schemes. Needless to say, one could as well use higher-order schemes like Crank-Nicolson, but
% such schemes lead to
% a more complicated matrix structure (with a couple of nonlinear terms also on the
% offdiagonals).
% To this end, we set $t_n=n\tau,\; n=0,1,\dots, N_{t},$
% with $\tau= T/N_{t}$  so that (\ref{kkt_cond_neweq1}) and (\ref{kkt_cond_neweq2}) yield
% \be
% \label{kktime1}
% %\left\{
% \begin{aligned}
% c_{kk} \frac{v^{n}_k - v^{n-1}_k }{\tau} - \nu  c_{kk}\Delta  v^n_k
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1}\left(\bar v^n_i \cdot \nabla\right) v^n_j c_{ijk}
% + \bigtriangledown c_{kk} p^n_k  &= -\frac{1}{\beta}\lambda^n_k  c_{kk} &&\\
% \nabla\cdot v^n_k &=0  && \nonumber\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{out}: \qquad  v^n_k &=0  \\
%  \mbox{on} \quad \Upsilon_{in}: \qquad c_{kk}v^n_k  &= \tilde{\boldsymbol \theta}^n_M \\
%      \frac{\partial v^n_k  }{\partial  {\bf n}}  &= 0, \\
%    v^0_k &= 0,
% %    \nonumber&\\
% \end{aligned}
% %\right.
% \e
% and
% \be
% \label{kktime2}
% %\left\{
% \begin{aligned}
% c_{kk} \frac{\lambda^n_k - \lambda^{n+1}_k}{\tau}  - \nu c_{kk}\Delta \lambda^n_k
% - \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\bar v^n_i  \cdot \nabla\right)\lambda^n_j c_{ijk} &  &&\\
% + \sum_{i=0}^{N_{\xi}-1}\sum_{j=0}^{N_{\xi}-1} \left(\nabla\bar v^n_i \right)^{T}\lambda^n_j c_{ijk}
% +  \bigtriangledown c_{kk}\mu^n_k &=
% c_{kk}\Delta  v^n_k &&\\
% \nabla\cdot \lambda^n_k &=0  &&\\
% \mbox{on} \quad \Upsilon_{wall}\cup \Upsilon_{in}: \qquad \lambda^n_k  &=0  \\
% \mbox{on} \quad \Upsilon_{out}: \qquad   \frac{\partial \lambda^n_k  }{\partial  {\bf n}} & =0 \\
%  \lambda^{N_t}_k & = 0.
% %    \nonumber&\\
%  \end{aligned}
% %\right.
% \e
%
% Now, we discretize  (\ref{kktime1}) and (\ref{kktime2}) with $Q_2$-$Q_1$ (Taylor-Hood) elements
% to  obtain the following  saddle point system (or the Oseen equation)  at every nonlinear Picard step \cite{SSNS17}
This system can be written in a compact matrix form
\begin{align}
\label{stokespa}
 \underbrace{\left[\begin{array}{ccc}
\bm{M}_1 & 0 & -\bm{K}^{*} \\
0 & \beta\bm{M}_2 & \bm{M_3}^T \\
-\bm{K} & \bm{M_3} & 0 \\
\end{array}\right]}_{:=\mathfrak{A}} 
\underbrace{\left[\begin{array}{c}
{\bf v}_h \\
{ p}_h \\
{\bf u}_h \\
{{\boldsymbol{\lambda}}}_h \\
{\mu}_h \\
\end{array}\right]}_{:=\bf y}=\underbrace{\left[\begin{array}{c}
% \tau\alpha\bm{M}_0\hat\vv \\ %_{N_t} \\
{\bf 0} \\
{\bf 0} \\
{\bf g}_h \\
\end{array}\right]}_{:=\bf b}.
\end{align}
Since we use tensor product basis functions,
each of the block matrices in $\mathfrak{A}$ can be represented via Kronecker products.
In particular, the forward and adjoint operators for the time-stochastic-space Navier-Stokes equations write as
\begin{eqnarray}
\label{eq:K3}
% \bm{K}= I_{N_t} \otimes \left[ \nu_0G_0 \otimes \mathcal{L} +  \nu_1G_1 \otimes \mathcal{L}_0\right]  +  \bm{N}[\hat \vv] + C\tau^{-1} \otimes G_0 \otimes  \mathcal{M},
\bm{K} & = & I_{N_t} \otimes I_{n_\xi}^{\otimes m} \otimes \mathcal{L}  + C\tau^{-1} \otimes I_{n_\xi}^{\otimes m} \otimes  \mathcal{M}  +  \bm{N}[\bar \vv_h], \\
\bm{K}^* & = & I_{N_t} \otimes I_{n_\xi}^{\otimes m} \otimes \mathcal{L}  + C^\top\tau^{-1} \otimes I_{n_\xi}^{\otimes m} \otimes  \mathcal{M}  +  \bm{W}[\bar \vv_h] - \bm{N}[\bar \vv_h],
\end{eqnarray}
where \cite{BSOS215}
\begin{itemize}
\item $I_{N_t}$ is the identity  matrix and $C=\mathrm{tridiag}(-1,1,0)$
comes from the implicit Euler discretization;
\item $I_{n_\xi}$ is the identity mass matrix for each variable $\xi$.
\item $\mathcal{L} = \begin{bmatrix}\nu L & B^\top \\ B & 0 \end{bmatrix}$ is the discretization of the stationary Stokes operator,
and $\mathcal{M}=\mathrm{blkdiag}(M, 0)$ for the velocity mass matrix, extended by zeros to the pressure space;
\item $\bm{N}$ and $\bm{W}$ correspond to the convection terms. They depend on the {\it low-rank} structure of the solution $\bar\vv_h$,
which will be introduced in the next section.
\end{itemize}
Furthermore, we introduce an extended velocity Laplace matrix $\mathcal{L}_0 = \mathrm{blkdiag}(L, 0)$.
Then we have
\begin{align}
\label{stoms}
%  \bm{M}_0 & = I_{N_t}\otimes I_{n_\xi}^{\otimes m}\otimes \mathcal{M}, &\quad
%  \bm{M}_1 &=   I_{N_t}\otimes I_{n_\xi}^{\otimes m}\otimes \left((1-\alpha)\mathcal{L}_0 + \alpha \mathcal{M}\right), \nonumber\\
 \bm{M}_1 &=   I_{N_t}\otimes I_{n_\xi}^{\otimes m}\otimes \mathcal{L}_0, &
 \bm{M}_2 &=   I_{N_t}\otimes I_{n_\xi}^{\otimes m} \otimes M, & \bm{M}_3 & = I_{N_t} \otimes I_{n_\xi}^{\otimes m}\otimes \begin{bmatrix}M_{cnt} \\ 0\end{bmatrix},
\end{align}
where
\begin{enumerate}
 \item $M_{cnt}=M$ for the distributed control, and
 \item $M_{cnt}(k,k') = \left\{\begin{array}{ll}M(k,k'), & \mbox{if }k,k' \in \Upsilon_{cnt}, \\ 0, & \mbox{otherwise}, \end{array}\right.$ for the boundary control.
\end{enumerate}

% and, in the  notation of \cite{PS12},
% \begin{equation}
% \label{eq:stokes}
% \bar{\mathcal{L}} =    \left[ \begin{array}{cc}
% \mathcal{A}& \mathcal{B}^T \\
% \mathcal{B} & 0
%  \end{array} \right]
% \end{equation}
% represents an instance of the time-dependent Stokes problem with
% \begin{eqnarray}
% \label{stokA}
% \mathcal{A} = G_0 \otimes A +  G_1 \otimes \nu_1 K, \quad A = \tau^{-1}  M  + \nu_0 K   + M_{\varrho}, \qquad \mathcal{B}=G_0 \otimes B,
% \end{eqnarray}
% and $G_i(j,k) = \left< \psi_i\psi_j(\xi)\psi_{k}(\xi)\right>.$
% % Also, observe  that
% % \begin{equation}
% % \label{eq:K3new}
% % \bm{K} = I_{n_t} \otimes G_0 \otimes \begin{bmatrix}A & B^T \\ B & 0\end{bmatrix} + I_{n_t} \otimes G_1 \otimes \begin{bmatrix}\nu_1 K & 0\\ 0 & 0\end{bmatrix} + C \otimes G_0 \otimes \begin{bmatrix}\tau^{-1} M & 0 \\ 0 & 0\end{bmatrix},
% % \end{equation}
% 
% \be
% \label{largeN}
% \bm{N} = I_{N_t} \otimes G_0 \otimes N.
% \e

% \begin{remark}
%  Note that, with the exception of block (2,2), the
% third matrices  in the Kronecker representations of  all the blocks in (\ref{stokespa}) are themselves also block matrices.
% \end{remark}
The third vector ${\bf g}_h$ in the right-hand side ${\bf b}$ in \eqref{stokespa} depends on the boundary conditions.
The KLE yields the following tensor form for the inflow function,
\begin{equation}
\boldsymbol\theta_h = \mathbf{g}_t \otimes \sum_{k=0}^{m} \left[ \left(\bigotimes_{\ell=1}^{k-1} \mathbf{e}_1\right) \otimes \mathbf{e}_2 \otimes \left(\bigotimes_{\ell=k+1}^{m} \mathbf{e}_1 \right) \otimes \boldsymbol\theta_{kh}\right],
\label{eq:inflow-discr}
\end{equation}
where $\mathbf{g}_t$ is the vector of the $1-\exp(-t)$ factor sampled at all time points, $\mathbf{e}_1$ and $\mathbf{e}_2$ are the first and second unit vectors, respectively, and $\boldsymbol\theta_{kh}$ is the discretization of the $k$-th spatial factor in \eqref{eq:kle_in}.
Now, distinguishing inner and boundary degrees of freedom, we partition the PDE matrix into the corresponding blocks,
$$
\bm{K} = \begin{bmatrix} \bm{K}_{II} & \bm{K}_{IB} \\ \bm{K}_{BI} & \bm{K}_{BB} \end{bmatrix}.
$$
The actual separation of inner (I) and boundary (B) elements is done in the spatial factors only.
Then the usual finite element approach is employed: we eliminate $\bm{K}_{IB}$ and $\bm{K}_{BI}$ in the left hand side of \eqref{stokespa},
replacing $\bm{K}$ by $\bm{K}_{II}$,
and construct the right-hand side from \eqref{eq:inflow-discr} as $\mathbf{g}_h = \bm{K}_{IB} \boldsymbol\theta_h.$


%
% \section{A low-rank solution method}
% %\subsection*{Curse of dimensionality}
% It is important to note here  that, due to the tensor product structure inherent
% in the blocks of the system \eqref{stokespa}, solving it poses quite a difficult challenge,   especially when the discretization of the OPNS is performed on very fine grids.
% Notice that the coefficients of the solution components in \eqref{stokespa} can be seen as tensors, depending on three indices, $y = \left[y(i,j,k)\right]$, where $i$ runs through different time steps, $j$ enumerates different stochastic polynomials and $k$ corresponds to the spatial degrees of freedom (DOFs).
% For instance, if the stochastic, temporal, and spatial degrees of freedom are given, respectively, by $N_{\xi}=2^8,\; N_t= 2^{10}$ and $J=10^6,$ then for each Picard iteration, one would need to solve a system with $ 3JN_{\xi}N_t \approx 7.8 \times 10^{11}$ DOFs!
% As a consequence, one expects overwhelming memory and computational time requirements.
%
% Considerable data compression can, fortunately,  be achieved using {\it low-rank} decompositions of high-dimensional tensor products.
% One of these decompositions is the robust {\it Tensor Train}
% (TT) format \cite{OSE11}. In this work, we will rely on  TT algorithms to solve our linear systems. Before introducing the TT solvers, we first
%  assume that all quantities are in low-rank format:
% \begin{alignat}{3}
% \label{eq:lr-state}
% \vv_h &=\sum_{i=1}^{r_{\vv_h}}t_{i,\vv_h}\otimes w_{i,\vv_h}\otimes z_{i,\vv_h}; \quad &p_h&=\sum_{i=1}^{r_{\pp_h}}t_{i,p_h}\otimes w_{i,p_h}\otimes z_{i,p_h}; \\
% \label{eq:lr-adj}
% \pp_h &=\sum_{i=1}^{r_{\pp_h}}t_{i,\pp_h}\otimes w_{i,\pp_h}\otimes z_{i,\pp_h}; \quad &\mu_h&=\sum_{i=1}^{r_{\mu_h}}t_{i,\mu_h}\otimes w_{i,\mu_h}\otimes z_{i,\mu_h}; \\
% \label{eq:lr-con}
% \uu_h &=\sum_{i=1}^{r_{\uu_h}}t_{i,\uu_h}\otimes w_{i,\uu_h}\otimes z_{i,\uu_h}. &&
% \end{alignat}
%
% % We now get the following system at every nonlinear step
% % \begin{equation}
% %  \left[
% % \begin{array}{ccc}
% % \Theta\otimes\left(\alpha_1 \mathcal{M}+\alpha_2\mathcal{L}_0\right) &0&\mathcal{K^*}\\
% % 0&\Theta\otimes \beta \M & -\mathcal{M}_3^\top\\
% % \mathcal{K}& -\mathcal{M}_3 &0\\
% % \end{array}
% % \right]
% % \begin{bmatrix}
% % \vv_{h}\\
% % p_h\\
% % \uu_{h}\\
% % \pp_{h}\\
% % \xi_h\\
% % \end{bmatrix} =
% % \begin{bmatrix}
% % 0 \\
% % 0 \\
% % 0 \\
% % \bar{d} \\
% % 0
% % \end{bmatrix} =
% % \begin{bmatrix}
% % b_1 \\
% % 0 \\
% % b_2 \\
% % b_3 \\
% % 0
% % \end{bmatrix},
% % \label{eq:kkt-matrix}
% % \end{equation}
% % where $\mathcal{L}_0 = \begin{bmatrix}\L & 0 \\ 0 & 0\end{bmatrix}$ represents the discretized version of $\mathbf{curl}(\mathbf{curl}(\vv))$, which is just the Laplacian operator, since the velocity is divergence-free.
% % Moreover, $b$ represents the right-hand side and $\mathcal{K}$ describes the forward operator for the space-time Navier-Stokes equations
% % \begin{equation}
% % \mathcal{K}=C\otimes \mathcal{M}+I_{n}\otimes\mathcal{L}+\sum_{i=1}^{r_{\bar\vv_h}}D_i\otimes\mathcal{N}_{i}.
% % \label{eq:K_ns}
% % \end{equation}
% % The adjoint PDE represented by $\mathcal{K^*}$ contains more terms than the forward equation due to the terms $\left(\vv\cdot \nabla\right)\pp+\left(\nabla\vv\right)^{T}\pp.$
% % As in the forward problem, we assume \eqref{eq:yprev_lr},
% % so that  $\left(\bar \vv \cdot \nabla\right)\pp$ is discretized as $\sum_{i=1}^{r_{\bar\vv_h}}D_i\otimes\mathcal{N}_{i}$.
% % The term $\left(\nabla \bar y\right)^{T}\pp$ now becomes $\sum_{i=1}^{r_{\bar y}}D_i\otimes\mathbf{H}_{i}(\bar{\vv}_h)$ with
% % $\mathbf{H}_{i}(\bar{\vv}_h)$ being a matrix of entries $\int {\boldsymbol\phi}_j\cdot\nabla \left(\sum_{k=1}^m w_{i,\bar\vv_h,k} {\boldsymbol\phi}_k\right)\cdot {\boldsymbol\phi}_{j'}$ for $j,j'=1,\ldots,m$.
% % The adjoint matrix is then given by
% % $$
% % \mathcal{K^*}=C^{\top}\otimes \mathcal{M}+I_{n}\otimes\mathcal{L}-\sum_{i=1}^{r_{\bar\vv_h}}D_i\otimes\mathcal{{N}}_{i}+\sum_{i=1}^{r_{\bar\vv_h}}D_i\otimes\mathcal{{H}}_{i},
% % $$
% % where $\mathcal{{H}}_{i} = \mathrm{blkdiag}(\mathbf{H}_{i}, 0)$.
%

 
\section{Tensor Train decomposition approach}
\label{tt_decom} 

The solution $y(t,\xi_1,\ldots,\xi_m,{\bf x})$ (where $y$ stands for $\vv,p,{\bf u},\pp$ or $\mu$) is a multivariate function.
After discretization, independently in each variable $t,\xi_1,\ldots,\xi_m,{\bf x}$,
the discrete values of $y$ can be enumerated by $D=m+2$ independent 
indices\footnote{We don't separate different \emph{components} of ${\bf x}$,
which are actually \emph{dependent}, due to the domain geometry. 
Therefore, all spatial degrees of freedom are treated as one variable.}, i.e. they form a \emph{tensor}.
Storing such a tensor directly might be prohibitively expensive.
Therefore, we \emph{decompose} tensors using the separation of variables.

Given a $D$-index tensor $y(i_1,\ldots,i_D)$, its Tensor Train decomposition \cite{osel-tt-2011} is written as follows:
\begin{equation}
y(\mathbf{i}) = \sum\limits_{s_1=1}^{r_1} \cdots \sum\limits_{s_{D-1}=1}^{r_{D-1}} y^{(1)}_{s_1}(i_1) y^{(2)}_{s_1,s_2}(i_2) \cdots y^{(D)}_{s_{D-1}}(i_D),
\label{eq:tt}
\end{equation}
where $\mathbf{i}$ denotes the multi-index, $\mathbf{i} = (i_1,\ldots,i_D)$.
Each element of $y$ is represented (or approximated) by a sum of products of elements of smaller tensors $y^{(\ell)}$, $\ell=1,\ldots,D$, called TT blocks.
The auxiliary summation indices $s_1,\ldots,s_{D-1}$ are called rank indices, and their ranges $r_1,\ldots,r_{D-1}$ are called TT ranks.
The TT decomposition is also known as the Matrix Product States \cite{schollwock-2005,PerezGarcia-mps-2007}, since, omitting $s_1,\ldots,s_{D-1}$, we can say that an element $y(\mathbf{i})$ is equal to a product of $i_\ell$-dependent matrices.
The TT ranks depend on the particular tensor and accuracy, if \eqref{eq:tt} is satisfied approximately.
Denoting upper bounds $r_\ell \lesssim r$, $i_\ell\lesssim N$, $\ell=1,\ldots,D-1$, we estimate the storage cost of the TT blocks to be in $\mathcal{O}(DNr^2)$, which can be much less than the full amount $N^D$.
Usually, we refer to $r$ as the maximal TT rank.

The TT format allows to compress a given tensor up to a threshold $\varepsilon$, 
also if $y$ is already given in the TT format, but possibly with overestimated TT ranks.
This happens if we sum two arrays in the TT format, or perform the matrix multiplication.
A matrix $A$, acting on tensors of type $y$, can be seen as a $2D$-dimensional tensor and represented in a slightly different TT format,
\begin{equation}
A(\mathbf{i},\mathbf{j}) = \sum\limits_{s_1,\ldots,s_{D-1}=1}^{R_1,\ldots,R_{D-1}} A^{(1)}_{s_1}(i_1,j_1) \cdots  A^{(D)}_{s_{D-1}}(i_D,j_D),
\label{eq:ttm}
\end{equation}
where $\mathbf{i},\mathbf{j}$ are multi-indices consisting of $i_\ell$ and $j_\ell$, respectively.
This is consistent with the Kronecker product ($\otimes$) if $D=2$ and $R=1$.
The matrix-vector product $\sum_{\mathbf{j}} A(\mathbf{i},\mathbf{j}) y(\mathbf{j})$ can then be implemented in the TT format block by block \cite{osel-tt-2011,schollwock-2005}.

\subsection*{TT format of the convection matrix}
We assume now that the velocity, in particular the previous Picard iterate, is represented by a TT decomposition \eqref{eq:tt},
\begin{equation}
\bar\vv_h(n, j_1,\ldots,j_m, k) = \bar\vv^{(1)}(n) \bar\vv^{(2)}(j_1) \cdots \bar\vv^{(m+1)}(j_m) \bar\vv^{(m+2)}(k),
\label{eq:hat_v}
\end{equation}
with the ranks $r_1,\ldots,r_{m+1}$.
We are going to derive a matrix TT format \eqref{eq:ttm} for $\bm{N}[\bar\vv_h]$ from \eqref{eq:K3}.
First, we notice from \eqref{eq:N_detailed} that $\bm{N}$ depends linearly on $\bar\vv_h$.
Therefore, we can plug \eqref{eq:hat_v} into \eqref{eq:N_detailed} and distribute the summations.
This gives
\begin{eqnarray}
N_{0,n}(jk,j'k') &  = & \bar\vv^{(1)}(n) \\
            \label{eq:N_tt_xi}
               &\cdot & \left[\sum_{j''=0}^{N_\xi-1} \int_{\Gamma} \psi_j(\xi)\psi_{j'}(\xi)\psi_{j''}(\xi) \rho(\xi)d\xi \cdot  \bar\vv^{(2)}(j''_1) \cdots \bar\vv^{(m+1)}(j''_m) \right] \\
               &\cdot & \left[\sum_{k''=1}^{N_v} \int_{\mathcal{D}} \phi_k({\bf x}) \nabla \phi_{k'}({\bf x}) \phi_{k''}({\bf x}) dx \cdot \bar\vv^{(m+2)}(k'') \right].
\end{eqnarray}
The last term in this expression is a deterministic convection matrix $N_s\left[\bar\vv^{(m+2)}\right]$ \eqref{eq:N_detailed},
which can be assembled with the cost $\mathcal{O}(r_{m+1} N_v)$.
The middle term can be further factorised due to the tensor product PCE space $\mathcal{Y}_{n_\xi}$.
Indeed, introducing the matrices of triple products for all $k=1,\ldots,m,$
\begin{equation}
H_{j_k''}(j_k,j_k') = \int\limits_{\Gamma_k} \psi_{j_k}(\xi_k)\psi_{j_k'}(\xi_k)\psi_{j_k''}(\xi_k) \rho_k(\xi_k)d\xi_k, \quad H_{j_k''} \in \mathbb{R}^{n_\xi \times n_\xi},
\label{eq:H3}
\end{equation}
we see that the multi-dimensional triple product in \eqref{eq:N_tt_xi} can be written as a Kronecker product of $m$ individual matrices \eqref{eq:H3}.
Finally, we notice that $N_{0,n}$ acts independently on each time step $n$, i.e. it is diagonal w.r.t. $n$, and so is $\bm{N}$.
Overall, given the TT format of $\bar\vv$ \eqref{eq:hat_v}, there exists a TT format for the convection matrix with the same TT ranks,
\begin{equation}
\begin{split}
\bm{N}[\bar\vv_h] &= \sum_{s_1,\ldots,s_{m+1} =1}^{r_1,\ldots,r_{m+1}}\mathrm{diag}(\bar\vv^{(1)}_{s_1}) \\
       & \otimes  \left[\sum_{j_1=0}^{n_\xi-1}H_{j_1}\bar\vv^{(2)}_{s_1,s_2}(j_1)\right] \otimes \cdots \otimes \left[\sum_{j_{m}=0}^{n_\xi-1}H_{j_{m}}\bar\vv^{(m+1)}_{s_{m},s_{m+1}}(j_{m})\right] \\
       & \otimes  \mathrm{blkdiag}\left(N_s\left[\bar\vv^{(m+2)}_{s_{m+1}}\right], N_s\left[\bar\vv^{(m+2)}_{s_{m+1}}\right]\right).
\end{split}
\label{eq:conv_tt}
\end{equation}
Notice that the right-hand side $\mathbf{g}_h = \bm{K}_{IB} \boldsymbol\theta_h$ in \eqref{stokespa} must be recomputed in every Picard iteration, since $\bm{K}_{IB}$ carries the corresponding (new) part $\bm{N}_{IB}$.
The adjoint convection matrix can be constructed similarly, with $W_s$ defined in \eqref{eq:Ws},
\begin{equation*}
\begin{split}
\bm{W}[\bar\vv_h] &= \sum_{s_1,\ldots,s_{m+1} =1}^{r_1,\ldots,r_{m+1}}\mathrm{diag}(\bar\vv^{(1)}_{s_1}) \\
       & \otimes  \left[\sum_{j_1=0}^{n_\xi-1}H_{j_1}\bar\vv^{(2)}_{s_1,s_2}(j_1)\right] \otimes \cdots \otimes \left[\sum_{j_{m}=0}^{n_\xi-1}H_{j_{m}}\bar\vv^{(m+1)}_{s_{m},s_{m+1}}(j_{m})\right]
        \otimes  W_s\left[\bar\vv^{(m+2)}_{s_{m+1}}\right].
\end{split}
\end{equation*}


\subsection*{Alternating linear solver}
The Oseen equation  given by \eqref{stokespa} is a large linear system for ${\bf y}$, and needs to be solved keeping all the components in the TT format in order to keep the storage requirements low.
A state of the art approach to this problem is alternating tensor product algorithms \cite{white-dmrg-1993,schollwock-2005,holtz-ALS-DMRG-2012}.
Given the system $\mathfrak{A}{\bf y}={\bf b}$, we
iterate over $\ell=1,\ldots,D$, and seek only the elements of $y^{(\ell)}$ in each step, while the other TT blocks are fixed.
Notice that the TT format \eqref{eq:tt} is \emph{linear} with respect to the elements of each $y^{(\ell)}$, i.e., there exists a matrix $Y_\ell$ \eqref{eq:frame} such that ${\bf y} = Y_\ell y^{(\ell)}$.
% subsequently fix all TT blocks of ${\bf y}$ and update only one or two.
This renders $\mathfrak{A}{\bf y}={\bf b}$ an overdetermined system $\mathfrak{A} Y_\ell y^{(\ell)}={\bf b}$ w.r.t. the elements of $y^{(\ell)}$.
This system is resolved via a projection onto $Y_\ell$, such that $y^{(\ell)}$ is computed from a smaller system $\left(Y_\ell^\top \mathfrak{A} Y_\ell\right) y^{(\ell)} = Y_\ell^\top {\bf b}$.
Particularly efficient realizations are the Density Matrix Renormalization Group (DMRG) methods from quantum physics \cite{white-dmrg-1993,jeckelmann-dmrgsolve-2002} and the Alternating Minimal Energy (AMEn) algorithm \cite{ds-amen-2014} from the mathematical community.
The DMRG approach computes two neighboring TT blocks in each step, say, $\ell$ and $\ell+1$, which allows to adapt the TT rank $r_\ell$ to the desired accuracy.
While the  DMRG method was found to be extremely effective for spin Schroedinger eigenvalue problems, it may deliver insufficient accuracy for linear systems, especially with non-symmetric matrices.
The AMEn method is usually faster and more robust, since it seeks only one TT block in each step, but then performs an explicit augmentation of the computed TT block of the solution by a TT block of the current residual.
This allows to change TT ranks and facilitate convergence.

However, the standard AMEn method is not applicable to $\mathfrak{A}$ directly: due to indefiniteness of $\mathfrak{A}$, its Galerkin projection may be degenerate.
For example, consider
$$
\mathfrak{A} = \begin{bmatrix}1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 0\end{bmatrix} \quad \mbox{and} \quad Y_\ell = \begin{bmatrix}0 \\ 0 \\ 1\end{bmatrix}.
$$
One can readily verify that $Y_\ell^\top \mathfrak{A} Y_\ell=0$.
To cope with this issue, we employ the so-called block TT format \cite{dkos-eigb-2014} and project each \emph{submatrix} of $\mathfrak{A}$ separately.

Let us enumerate the components of \eqref{stokespa} as $y_{\iota}$, where
$y_1 = \vv_h$, $y_2=p_h$, $y_3={\bf u}_h$, $y_4=\pp_h$ and $y_5=\mu_h$.
We approximate all components simultaneously by a TT format with the same blocks except the $\ell$-th one for some $\ell=1,\ldots,D$, where the enumerator $\iota=1,\ldots,5$ appears,
\begin{equation}
y_{\iota}(\mathbf{i}) = \sum_{s_1,\ldots,s_{D-1}} y^{(1)}_{s_1}(i_1) \cdots y^{(\ell)}_{s_{\ell-1,s_\ell}}(i_\ell,\iota) \cdots y^{(D)}_{s_{D-1}}(i_D).
\label{eq:btt}
\end{equation}
Using an SVD of the core $y^{(\ell)}$, we can turn \eqref{eq:btt} into the form
$$
y_{\iota}(\mathbf{i}) = \sum_{s_1,\ldots,s_{D-1}} y^{(1)}_{s_1}(i_1) \cdots y^{(\ell)}_{s_{\ell-1},s_\ell}(i_\ell) \cdot y^{(\ell+1)}_{s_{\ell},s_{\ell+1}}(i_{\ell+1},\iota) \cdots y^{(D)}_{s_{D-1}}(i_D)
$$
or vice versa.
That is, we can \emph{move} $\iota$ to any particular TT block in the course of the alternating iteration \cite{dkos-eigb-2014}.
When the $\ell$-th block carries $\iota=1,\ldots,5$, the Oseen system \eqref{stokespa} is projected onto the other TT blocks component by component, giving a system
\begin{equation}
\label{stokespa-red}
\begin{bmatrix}
\hat M_1 & 0 & -\hat K^{*} \\
0 & \beta \hat M_2 & \hat M_3^T \\
-\hat K & \hat M_3 & 0 \\
\end{bmatrix}
y^{(\ell)}=
\begin{bmatrix}
% \alpha\tau Y_\ell^\top(\bm{M}_0 \hat \vv) \\
0 \\
0 \\
\hat g \\
\end{bmatrix}
\end{equation}
on the components of the current TT block.
Here, $\hat A = Y_\ell^T {\bm A} Y_\ell$ for $\bm A \in \{\bm{K},\bm{M}_1,\bm{M}_2,\bm{M}_3\}$ are the submatrices projected onto the Galerkin basis $Y_\ell$ composed from the frozen TT blocks,
\begin{equation}
\begin{split}
Y_\ell(\mathbf{i},s_{\ell-1} i'_\ell s_\ell) & = y^{(1)}(i_1) \cdots y^{(\ell-1)}_{:,s_{\ell-1}}(i_{\ell-1}) \cdot I_{n_\ell}(i_\ell,i_\ell') \cdot y^{(\ell+1)}_{s_\ell,:}(i_{\ell+1}) \cdots y^{(D)}(i_D),
\end{split}
\label{eq:frame}
\end{equation}
where ``$:$'' stands for the full index range.


Under certain assumptions on the original system, one can prove that the block-reduced system \eqref{stokespa-red} is not degenerate.
\begin{theorem}
Suppose that the symmetric parts of $\bm{K}$ and $\bm{K}^*$ are positive, $\bm{K}+\bm{K}^\top>0$, $\bm{K}^*+\bm{K}^{*\top}>0$.
Then the reduced matrix in \eqref{stokespa-red} is invertible.
\end{theorem}
\begin{proof}
Due to the Poincar\'e theorem, the eigenvalues of an orthogonal projection of a symmetric matrix interlace with the eigenvalues of the original matrix.
In particular,
$$
\lambda_{\min}(\hat K + \hat K^\top) = \lambda_{\min}\left(Y_k^\top (\bm{K}+\bm{K}^\top) Y_k \right) \ge \lambda_{\min}(\bm{K}+\bm{K}^\top)>0.
$$
We also use that $Y_k$ is orthogonal in the AMEn method.
So, the symmetric part of $\hat K$ (as well as of $\hat K^*$) is positive.
Moreover, by the same interlace theorem we have that $\hat M_1 \ge 0$ and $\hat M_2>0$.
Now employ \cite[Theorem 3.2]{BenGolLie05},
which says that sufficient conditions for the KKT matrix to be invertible are that
the matrix of constraints $\begin{bmatrix}-\hat K & \hat M_3\end{bmatrix}$ is full rank, and
$$
\mathrm{ker}\begin{bmatrix}\hat M_1 & 0 \\ 0 & \beta\hat M_2\end{bmatrix} \cap \mathrm{ker}\begin{bmatrix}-\hat K & \hat M_3\end{bmatrix} = \{0\}.
$$
The first condition is fulfilled since $\hat K$ is invertible.
To verify the second criterion, consider a vector in the kernel of the constraints, which has the form
$$
w = \begin{bmatrix}y \\ u\end{bmatrix} = \begin{bmatrix}\hat K^{-1} \hat M_3 u \\ u\end{bmatrix}, \qquad \forall u \neq 0.
$$
Now check if it belongs to the kernel of the other matrix:
$$
w^\top \begin{bmatrix}\hat M_1 & 0 \\ 0 & \beta\hat M_2\end{bmatrix} w =
u^\top \hat M_3^\top \hat K^{-\top} \hat M_1 \hat K^{-1} \hat M_3 u + \beta u^\top \hat M_2 u > 0,
$$
since $u^\top \hat M_2 u>0$ while the first term is non-negative.
Together with positive semi-definiteness of $\mathrm{blkdiag}(\hat M_1,~ \beta\hat M_2),$ this yields that $w$ is not in its kernel.

\hfill\end{proof}


The second modification to the AMEn algorithm is the unification of sizes and elimination of pressures.
Note that the last ($D$-th) TT block carries the spatial degrees of freedom,
and different components of $y_{\iota}$ have different numbers of spatial basis functions: $N_v$ for $\vv_h$, ${\bf u}_h$ and $\pp_h$, and $N_p$ for $p_h$ and $\mu_h$.
The Galerkin matrix $Y_\ell$ has always the same number of rows; we make it $N_t n_\xi^m N_v$, i.e. consistent with $\vv_h$, ${\bf u}_h$ and $\pp_h$.

The state and adjoint pressures are on the one hand more difficult, since they are of different sizes, and moreover make $\bm{K}$ itself saddle-point, but on the other hand they are only connected with the other components via a Kronecker-rank-1 matrix $\bm{B} = I_{N_t} \otimes I_{n_\xi}^{\otimes m} \otimes B$.
Therefore, we can compute the pressure components only within the $D$-th, the spatial, block.
When we proceed to $\ell<D$ (time and stochastic variables), we eliminate the pressures in a Gauss-Seidel fashion: having originally
$
\bm{K} = \begin{bmatrix}\bm{A} & \bm{B}^T \\ \bm{B} & 0\end{bmatrix},
$
we project only the velocity part $\hat K = Y_\ell^T \bm{A} Y_\ell$ in the left hand side,
and cast $Y_\ell^T \bm{B}^T p_h$ and $Y_\ell^T \bm{B}^T \mu_h$ to the right hand side.
The consequence is two-fold: the velocity part $\bm{A}$ is positive definite, hence so is $\hat K$, and all sizes are now consistent.
% The algorithm is thus as follows \cite{BSOS215}:
%
% 1. Project submatrices of \eqref{stokespa} onto $Y_\ell$.
%
% 2. If $\ell=D$: solve \eqref{stokespa-red} with pressure components.
%
% 3. If $\ell<D$: solve \eqref{stokespa-red} without pressure components but adding $-Y_\ell^T \bm{B}^T p$ and $-Y_\ell^T \bm{B}^T \mu$ to the right-hand side.
%
% 4. Compute SVD of $y^{(\ell)}$, switch to the next TT block ($\ell+1$ or $\ell-1$) and go to Step 1.


\subsection*{Preconditioning}
Even the reduced system \eqref{stokespa-red} can still be rather large: for $\ell=D$ for example, it is of size $(3 N_v + 2N_p) r_{m+1}$.
Therefore, we use GMRES with the matching preconditioner \cite{SW13,BSOS215}, based on approximating the Schur complement.
The reduced matrix admits a (straightforwardly verifiable) decomposition
\be
\label{shurfac}
\begin{bmatrix}
\hat M_1 & 0 & -\hat K^{*} \\
0 & \beta \hat M_2 & \hat M_3^T \\
-\hat K & \hat M_3 & 0 \\
\end{bmatrix} =
\begin{bmatrix}
I & * & * \\
 & I  \\
 &  & I \\
\end{bmatrix}
\begin{bmatrix}
 & & -\hat S  \\
 & \beta \hat M_2 & \hat M_3^T  \\
-\hat K & \hat M_3 & 0 \\
\end{bmatrix},
\e
where $\hat S = \hat K^* + \hat M_1 \hat K^{-1} \hat M_3 (\beta \hat M_2)^{-1} \hat M_3^T$.
Next, we use the second (anti-triangular) factor in  (\ref{shurfac}) as a preconditioner.
However, $\hat K$ and $\hat S$ must be approximated to make solving linear systems with them feasible.

The first matrix has the Kronecker form
$$
\hat K = \hat I \otimes \mathcal{L} + \hat C \otimes \mathcal{M} + \sum_{s=1}^{r_{m+1}} \hat D_s \otimes \mathcal{N}_s,
$$
where $\hat I$, $\hat C$ and $\hat D_s$ are $I_{N_t} \otimes I_{n_\xi}^{\otimes m}$, $C\tau^{-1} \otimes I_{n_\xi}^{\otimes m}$ and the parts of the convection matrix \eqref{eq:conv_tt}, respectively, projected via $Y_{m+2}$.
We approximate $\hat K$ by the following Sylvester operator:
$$
\tilde K = \hat I \otimes \left(\mathcal{L} + \sum_{s=1}^{r_{m+1}} \mathbb{E}\lambda(\hat D_s) \cdot \mathcal{N}_s \right) + \hat C \otimes \mathcal{M},
$$
where $\mathbb{E}\lambda(\hat D_s)$ is an average of the eigenvalues of $\hat D_s$ (it is symmetric, since so is the time-stochastic part of $\bm{N}$).
Now $\tilde K$ can be inverted by the Bartels-Stewart method, since $\hat C$ is small and can be easily Schur-factorized.

For the Schur complement $\hat S$ we use the matching factorisation,
\begin{equation}
\tilde S = \left(\hat K^* + \frac{1}{c} \hat M_1\right) \hat K^{-1} \left(\hat K + c \cdot \hat M_3 (\beta \hat M_2)^{-1} \hat M_3^T\right),
\label{eq:sc2}
\end{equation}
where $c>0$ is a normalization constant.
Traditionally it is proposed \cite{SW13,BSOS215} to take $c=\sqrt{\beta}$, and this is indeed an optimal choice for the distributed control.
For the boundary control, however, $\hat M_3 (\beta \hat M_2)^{-1} \hat M_3^T$ is significantly rank-deficient, and an optimal $c$ may differ.
We look for $c$ that minimizes the residual of the system $(\tilde S^{-1} \hat S) u = \tilde S^{-1} f$ with a random right-hand side $f$ after $5$ GMRES iterations.
Since it is difficult to differentiate the GMRES residual w.r.t. $c$, we employ the zero-order \emph{golden section} optimization algorithm \cite{Cheney}, initialized with an interval $\log_{10} c \in [\log_{10}\sqrt{\beta}-6, \log_{10}\sqrt{\beta}+6]$.
Fortunately, it is sufficient to perform this procedure in the first Picard iteration only, as the optimal $c$ does not seem to change in the latter.
To solve systems with $\tilde S$ efficiently, we also approximate the factors in brackets in \eqref{eq:sc2} by Sylvester matrices, similarly to $\tilde K$.

The full Algorithm \ref{alg} combines the standard Picard iteration for the Navier-Stokes equation \cite{ESW14} and the AMEn iteration in the block TT format \cite{BSOS215}.

\begin{algorithm}[t]
\caption{Block AMEn-Picard iteration for solving the stochastic inverse Navier-Stokes equations in the TT format}
\label{alg}
\begin{algorithmic}[1]
\State Initialize $\vv_h={\bf u}_h=\pp_h=0$, $p_h =\mu_h=0$.
\State Initialize the block TT format \eqref{eq:btt} with $\iota$ placed in $y^{(m+2)}$.
\For{iter=1,2,\ldots}
  \State Copy $\bar\vv_h = \vv_h, \bar p_h = p_h, \bar{\bf u}_h = {\bf u}_h, \bar\pp_h = \pp_h, \bar\mu_h = \mu_h$.
  \State Construct the convection matrices $\bm{N}[\bar\vv_h]$ and $\bm{W}[\bar\vv_h]$ \eqref{eq:conv_tt}.
%   \State Orthogonalise $y^{(1)},\ldots,y^{(m+1)}$ and construct $Y_{m+2}$.
  \State Construct and solve the projected system \eqref{stokespa-red} for $\ell=m+2$.
  \State Extract individual components $\vv_h=y_1,p_h=y_2,{\bf u}_h=y_3,\pp_h=y_4,\mu_h=y_5$.
  \If{$\|\vv_h-\bar\vv_h\|\le \eps \|\vv_h\|$, $\|p_h-\bar p_h\|\le \eps\|p_h\|$ and $\|{\bf u}_h-\bar {\bf u}_h\|\le \eps\|{\bf u}_h\|$}
    \State Stop.
  \EndIf
  \State Assemble $y_{\iota}$ from $y_1=\vv_h,y_2={\bf u}_h,y_3=\pp_h$ only.
  \State Apply the SVD to $y^{(m+2)}$ and move $\iota$ to $y^{(m+1)}$.
  \For{$\ell=m+1,m,\ldots,1,2,\ldots,m+1$}
    \State Construct and solve an analog of \eqref{stokespa-red} with fixed pressures,
      \begin{equation*}
      \begin{bmatrix}
      \hat M_1 & 0 & -\hat A^{*} \\
      0 & \beta \hat M_2 & \hat M_3^T \\
      -\hat A & \hat M_3 & 0 \\
      \end{bmatrix}
      y^{(\ell)}=
      \begin{bmatrix}
      -Y_\ell^T \bm{B}^T \mu_h \\
      0 \\
      \hat g - Y_\ell^T \bm{B}^T p_h \\
      \end{bmatrix}
      \end{equation*}
     \State Apply the SVD to $y^{(\ell)}$.
     \If{$\ell$ is increasing}
       \State Move $\iota$ to $y^{(\ell+1)}$.
     \Else
       \State Move $\iota$ to $y^{(\ell-1)}$.
     \EndIf
  \EndFor
\EndFor
\end{algorithmic}
\end{algorithm}

%%%%%%%%%%%%%
\section{Numerical results}
\label{experiments}
We implemented the computational codes on the basis of the Matlab TT-Toolbox \cite{tt-toolbox} and
the IFISS 3.3 toolbox \cite{ifiss3.4} and run on one core of the \texttt{otto} cluster at MPI Magdeburg,
an Intel Xeon X5650 @ 2.67GHz.
The main focus of this paper is the numerical scheme, hence we vary the model and discretization parameters one by one.
The default parameters, unless otherwise stated, are presented in Table \ref{tab:def}.

\begin{table}[h!]
\centering
\caption{Default model and discretization parameters}
\label{tab:def}
\begin{tabular}{ccccccccc}
$\eps$    & $\nu$ & $\beta$  &  $N_t$    & $T$  & $h$      & $n_\xi$ & $\gamma$ & $m$  \\
\hline
$10^{-4}$ & $0.1$ & $10^{-2}$&  $2^{10}$ & $30$ & $2^{-3}$ & $8$       & $4$      & $4$ \\
\end{tabular}
\end{table}

The natural indicators of computational complexity are the CPU time and the number of iterations, while the storage complexity of the TT representation is governed by the TT ranks.
We show the maximal TT rank and the memory saving ratio of the TT format, compared to the full tensor representation, i.e.
$$
\mbox{MemR} = \frac{N_t r_1 + \sum_{\ell=1}^m r_{\ell}r_{\ell+1} n_\xi + r_{m+1} N_v}{N_t n_\xi^m N_v}.
$$
The smaller the TT ranks, the smaller is this ratio, which indicates the storage savings.

As a quantity of interest, we compute the average squared norm of vorticity 
$\|\nabla \times \vv\|^2$, where the norm is taken in $L^2([0,T]) \times L^2(\Omega) \times L^2(\mathcal{D})$, i.e. the doubled first term of the functional \eqref{J}.

The largest problem was solved with $m=8,\; N_t=1024$ and $h=2^{-4};$ this value of $h$ gives $N_v=23042$ and $N_p=2945.$
Thus, the number of unknowns in the full KKT problem would be $(3N_v+2N_p)\cdot {n_\xi}^m\cdot N_t = 1.29\times 10^{15}.$
Note that the number of TT elements is $(3N_v+2N_p)r +  m{n_\xi}r^2 + N_tr;$
but since the first term is dominating,  it suffices to estimate the whole cost with the first term.
The maximum ranks are $70$ for the distributed control  and $120$  for the boundary control cases (see Figures \ref{fig:m} and \ref{fig:gamma_m_bnd} right, respectively).
This gives roughly $5$ and $9$ millions of unknowns.
%respectively in the spatial TT block.


% \pgfplotsset{%
% every tick label/.append style={font=\LARGE},%
%   every axis label = {font=\Huge},%
%   label style = {font=\Huge}%
% }
%
% \begin{figure}[h!]
% \centering
% \caption{Performance vs. truncation tolerance $\varepsilon$}
% \label{fig:eps}
% %\includegraphics[width=0.95\textwidth,height=0.5\textwidth]
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={-2,-3,-4,-5},
%   ymode=normal,
%   xlabel={$log_{10}\varepsilon$},
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol
%   ]
%
%   \addplot+[] coordinates {(-2,134.81) (-3,606.89) (-4,1553.59) (-5,4952.14)};
%   \node at (axis cs:-4.9,4000) [rotate=-55,anchor=west,color=black] {\LARGE CPU time};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=20,ymax=82,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(-2,20) (-3,35) (-4,65) (-5,78)};
%   \node at (axis cs:-2.1,25) [rotate=-32,anchor=east,color=blue] {\LARGE TT rank};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=0.4,ymax=1.6,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(-2,0.36) (-3,0.65) (-4,1.10) (-5,1.58)};
%   \node at (axis cs:-3.5,0.8) [rotate=-35,anchor=east,color=red] {\LARGE \% Mem};
%   \end{axis}
% \end{tikzpicture}}
% \resizebox{0.393\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={-2,-3,-4,-5},
%   ymode=normal,
%   xlabel=$log_{10}\varepsilon$,
%   ylabel=$\log_{10}\mbox{error}$,
%   y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
%   legend style={at={(0.99,0.01)},anchor=south east},
%   cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
%   ]
%   \addplot+[] coordinates {(-2,6e-3) (-3,3e-4) (-4,8e-5) (-5,6e-6)};\addlegendentry{\LARGE $\vv$};
%   \addplot+[] coordinates {(-2,24e-3) (-3,22e-4) (-4,23e-5) (-5,26e-6)};\addlegendentry{\LARGE p};
%   \addplot+[] coordinates {(-2,17e-3) (-3,12e-4) (-4,19e-5) (-5,15e-6)};\addlegendentry{\LARGE ${\bf u}$};
%   \end{axis}
% \end{tikzpicture}}
% \end{figure}
%
%
% \begin{figure}[h!]
% \centering
% \caption{Performance for different viscosity parameters. Left: $\nu_0=1/10$, $\nu_1$ is varied. Right: $\nu_1$ and $\nu_0$ are varied together as $\nu_1=0.25\nu_0$}
% \label{fig:nu}
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={2,4,6,8,10},
%   ymode=normal,
%   xlabel=$\nu_0/\nu_1$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol
%   ]
%   \addplot+[] coordinates {(2,1117.37) (2.5,726.56) (5,472.37) (10,421.40)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={2,4,6,8,10},
%   ymode=normal,ymin=30,ymax=58,
%   xlabel=$\nu_0/\nu_1$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(2,55) (2.5,49) (5,39) (10,35)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={2,4,6,8,10},
%   ymode=normal,ymin=0.5,ymax=1.2,
%   xlabel=$\nu_0/\nu_1$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(2,1.07) (2.5,0.94) (5,0.73) (10,0.65)};
%   \end{axis}
% \end{tikzpicture}}
% %
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={10,20,30,40,50},
%   ymode=normal,
%   xlabel=$1/\nu_0$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol
%   ]
%   \addplot+[] coordinates {(10,552.16) (20,1210.17) (30,2468.01) (40,3354.38) (50,5055.16)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={10,20,30,40,50},
%   ymode=normal,ymin=40,ymax=92,
%   xlabel=$1/\nu_0$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(10,43) (20,57) (30,70) (40,82) (50,87)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={10,20,30,40,50},
%   ymode=normal,ymin=0.7,ymax=1.9,
%   xlabel=$1/\nu_0$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(10,0.81) (20,1.11) (30,1.40) (40,1.68) (50,1.80)};
%   \end{axis}
% \end{tikzpicture}}
% \end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Distributed control}
\subsubsection{Discretization parameters}
In the first series of tests we consider the distributed control, applied in the right hand side of \eqref{kkt_cond1}.
We start with verifying correctness of the discretization schemes.
We show the results for different spatial mesh sizes in Fig. \ref{fig:nc}, and for different time step sizes in Fig. \ref{fig:nt}.

The number of spatial degrees of freedom grows quadratically with the number of steps in each dimension.
Since we use a direct solver for the spatial Stokes-like problems in \eqref{eq:sc2}, the computational complexity is higher;
in particular, it grows proportionally to $h^{-2.5}$ in Fig. \ref{fig:nc}.
To estimate the discretization error, we take the vorticity from the finest grid $h=2^{-4}$ as the reference value $\|\nabla \times \vv_*\|$,
and compare it with the values on coarser grids.
We observe the convergence order $h^{1.5}$,
which is reasonable for the backstep domain, 
containing a corner with an angle 
$3\pi/2;$  see e.g., \cite{babuska-regularity-1988,babuska-hp-1988}.

In order to keep the same number of stochastic variables 
throughout the spatial mesh test, we restrict it to $m=2$.
This is necessary since the higher KLE components 
in \eqref{eq:kle_in} are not resolved on the coarsest grid in 
this test,
and the problem would be in a severely pre-asymptotic regime 
otherwise.

The TT ranks grow very moderately with the grid sizes, with the asymptotic dependence close to logarithmic.
Since the full tensor storage, proportional to both $N_v$ and $N_t$ grows much faster, this leads actually to a decrease of the memory ratio with both spatial (Fig. \ref{fig:nc}) and time (Fig. \ref{fig:nt}) grid refinement.

The CPU time versus the time grid size demonstrates two regions.
For smaller numbers of time steps, solution of the spatial problems is dominating in the total complexity, while the TT ranks grow logarithmically.
This leads to a logarithmic growth of the CPU time as well.
When the number of time steps is large, the temporal problems become the most time consuming, and the CPU time starts
to grow linearly with $N_t$ (the actual estimated order is $1.2$ due to the additional logarithmic growth of the ranks).

The implicit Euler time discretization manifests an expected convergence of the first order.
Here the reference vorticity $\|\nabla \times \vv_\star\|$ is computed at the $N_t=2^{15}$ grid.


% \subsection{Varying spatial parameter}
% Three levels from ifiss starting from smallest
%
% Here we use $d=2$ to make it the same throughout all tests.

% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for distributed control with $\alpha=0.$}
% \begin{tabular}{ccccc}
% \hline
% $nc \;(J)$ &Time& Pic iter&Max Rank (r)& mem(Y)/prod(Y.n)\\
% \hline
% $2 ()$    & 64.39    &  3  & 11  & 9.95e-06  \\
% $3 ()$    & 82.16    &  3  & 22  & 5.72e-06  \\
% $4 ()$    & 316.08   &  5  & 33  & 3.42e-06  \\
% $5 ()$    & 1521.12  &  4  & 38  & 3.08e-06  \\
% $6 ()$    & 9803.10  &  4  & 40  & 2.91e-06  \\
% \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $nc \;(J)$ &Time& Pic iter  &Max Rank (r)& mem(Y)/prod(Y.n) & vort & J & mass \\
% \hline
% $3 ()$     & 59.7220    &  3  & 21  & 1.9847e-03 & 11.882014339  &   11.883688371  &  1.949985838 \\
% $4 ()$     & 171.1626   &  3  & 31  & 1.8313e-04 & 12.002581938  &   12.004670607  &  1.956032852 \\
% $5 ()$     & 881.3946   &  3  & 33  & 1.4929e-04 & 12.071493602  &   12.073847446  &  1.959280558 \\
% $6 ()$     & 5.8359e+03 &  3  & 35  & 1.4039e-04 & 12.097289795  &   12.099765617  &  1.960513784 \\
% \hline
% \end{tabular}

\begin{figure}[h!]
\caption{CPU time and convergence w.r.t. grid refinement (left) and TT rank and memory reduction ratio (right) for different spatial mesh sizes. Change of defaults: $m=2$.}
\label{fig:nc}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  xlabel={$-\log_2 h$},
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=$\log_2\mbox{CPU Time}$, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)/log10(2.0)}\pgfmathresult},
  ]

  \addplot+[] coordinates {
   (1, 59.7220   )
   (2, 171.1626  )
   (3, 881.3946  )
   (4, 5.8359e+03)}; \addlegendentry{$\log_2\mbox{CPU Time}$}
    % Just add some legend
   \addplot+[] coordinates {(0,1)}; \addlegendentry{$\log_2\left(\|\nabla \times \vv_{\star}\|^2 - \|\nabla \times \vv\|^2\right)$};

   \addplot+[line width=1pt,color=orange,solid,no marks,domain=1.5:4] {2^(2.5*x+2.4)};
   \node at (axis cs:1.7,6.5) [anchor=west,color=orange] {$2.5x+2.4$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_2\left(\|\nabla \times \vv_{\star}\|^2 - \|\nabla \times \vv\|^2\right)$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(12.097289795 - \pgfmathresult)/log10(2.0)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (1, 11.882014339)
   (2, 12.002581938)
   (3, 12.071493602)};

   \addplot+[line width=1pt,color=magenta,solid,no marks,domain=1:3.5] {12.097289795 - 2^(-1.5*x-0.5)};
   \node at (axis cs:3.2,-5.2) [anchor=west,color=magenta] {$-1.5x-0.5$};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  legend style={at={(0.99,0.50)},anchor=east},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$-\log_2 h$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (1, 21)
   (2, 31)
   (3, 33)
   (4, 35)}; \addlegendentry{max. TT rank};
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$}; % from the 2nd axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (1, 1.9847e-03)
   (2, 1.8313e-04)
   (3, 1.4929e-04)
   (4, 1.4039e-04)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}


% \subsection{Varying $N_t$}
% Three levels as well.
~

% \begin{table}[h!]
% \centering
% \caption{The temporal parameter $N_t$ for distributed control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $N_t $      & Time  & Pic iter & Max Rank (r)& norm(vort) &mem(Y)/prod(Y.n)\\
% \hline
% $2^{8}$     &  1182.1  & 2    & 41 & 12.0908  & 1.18e-05\\
% $2^{9}$     &  1216.2  & 2    & 44 & 12.0784  & 6.00e-06\\
% $2^{10}$    &  1139.5  & 2    & 47 & 12.0722  & 3.08e-06\\
% $2^{11}$    &  1152.3  & 2    & 50 & 12.0691  & 1.79e-06\\
% $2^{12}$    &  1594.8  & 2    & 54 & 12.0676  & 1.06e-06\\
% $2^{13}$    &  2213.5  & 2    & 57 & 12.0668  & 6.34e-07\\
% $2^{14}$    &  4597.8  & 2    & 61 & 12.0664  & 4.20e-07\\
% \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $N_t $      & Time       & Pic iter & Max Rank (r)& mem(Y)/prod(Y.n)  & vort & J & mass \\
% \hline
% $2^{7}$     & 978.8806   & 3        & 41          & 2.3510e-05        & 12.115638555  &   12.117999645  &  1.966325251 \\
% $2^{8}$     & 996.5336   & 3        & 43          & 1.1838e-05        & 12.090841716  &   12.093199204  &  1.962301335 \\
% $2^{9}$     & 1.1932e+03 & 3        & 46          & 6.0018e-06        & 12.078438998  &   12.080794523  &  1.960288675 \\
% $2^{10}$    & 1.2139e+03 & 3        & 50          & 3.0838e-06        & 12.072237403  &   12.074591819  &  1.959282297 \\
% $2^{11}$    & 1.2351e+03 & 3        & 55          & 1.7945e-06        & 12.069136609  &   12.071490372  &  1.958779099 \\
% $2^{12}$    & 1.5197e+03 & 3        & 60          & 1.0611e-06        & 12.067586210  &   12.069939596  &  1.958527496 \\
% $2^{13}$    & 2.2695e+03 & 3        & 65          & 6.3420e-07        & 12.066811001  &   12.069164180  &  1.958401693 \\
% $2^{14}$    & 5.2919e+03 & 3        & 71          & 4.2072e-07        & 12.066423392  &   12.068776462  &  1.958338790 \\
% \hline
% \end{tabular}

\begin{figure}[h!]
\caption{CPU time and convergence w.r.t. grid refinement (left) and TT rank and memory reduction ratio (right) for different numbers of time steps.}
\label{fig:nt}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=7,xmax=14,
  xlabel={$\log_2 N_t$},
  legend style={at={(0.60,0.99)},anchor=north},
  ylabel=$\log_2\mbox{CPU Time}$, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)/log10(2)}\pgfmathresult},
  ]

  \addplot+[] coordinates {
   (7 , 978.8806  )
   (8 , 996.5336  )
   (9 , 1.1932e+03)
   (10, 1.2139e+03)
   (11, 1.2351e+03)
   (12, 1.5197e+03)
   (13, 2.2695e+03)
   (14, 5.2919e+03)}; \addlegendentry{$\log_2\mbox{CPU Time}$};
   \addplot+[] coordinates{(0,1)}; \addlegendentry{$\log_2\left(\|\nabla \times\vv\|^2 - \|\nabla \times\vv_{\star}\|^2\right)$}; % from 2 axis

   \addplot+[line width=1pt,color=orange,solid,no marks,domain=11.5:14] {2^(1.2*x-4.4)};
   \node at (axis cs:11.8,9.8) [anchor=east,color=orange] {$1.2x-4.4$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=7,xmax=14,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_2\left(\|\nabla \times\vv\|^2 - \|\nabla \times\vv_{\star}\|^2\right)$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult - 12.066229586)/log10(2.0)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (7 ,    12.115638555)
   (8 ,    12.090841716)
   (9 ,    12.078438998)
   (10,    12.072237403)
   (11,    12.069136609)
   (12,    12.067586210)
   (13,    12.066811001)};

   \addplot+[line width=1pt,color=magenta,solid,no marks,domain=7:13.5] {12.066229586 + 2^(-1*x+2.6)};
   \node at (axis cs:7.5,-4.7) [anchor=west,color=magenta] {$-x+2.6$};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=7,xmax=14,
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$\log_{2}N_t$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (7 , 41)
   (8 , 43)
   (9 , 46)
   (10, 50)
   (11, 55)
   (12, 60)
   (13, 65)
   (14, 71)}; \addlegendentry{max. TT rank};
   \pgfplotsset{cycle list shift=1};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (7 , 2.3510e-05)
   (8 , 1.1838e-05)
   (9 , 6.0018e-06)
   (10, 3.0838e-06)
   (11, 1.7945e-06)
   (12, 1.0611e-06)
   (13, 6.3420e-07)
   (14, 4.2072e-07)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}




\subsubsection{Robustness with respect to model parameters: viscosity, KLE decay rate and dimension}
Having checked the discretization schemes, we can vary the other model parameters and see how the performance and the quantity of interest behave.
In Fig. \ref{fig:nu} we vary the viscosity from $10^{-1}$ to $10^{-4}$ and track the computational time and the number of Picard iterations (left), as well as the storage indicators and the quantity of interest (right).
We see that the CPU time grows only logarithmically with the Reynolds number, and $3$ Picard iterations are enough for the nonlinear system to converge.
This is also the case for other experiments with a \emph{sufficient} control, i.e. with the distributed control and small enough $\beta$, that was observed previously \cite{SSNS17} for a deterministic problem as well.
The TT rank (and hence memory ratio) and the vorticity remain at almost the same level (compared to variation of other parameters and other tests), which indicates that the minimal-vorticity solution is similar to the Stokes flow, fully defined by the domain and boundary conditions, rather than the viscosity.


% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for distributed control with $\alpha=0.$}
% \begin{tabular}{ccccc}
% \hline
% $\nu $&Time& Pic iter&Max Rank (r)& mem(Y)/prod(Y.n)\\
% \hline
% 1/10    &  2581.26 &  6  & 40 & 3.08e-06 \\
% 1/100    & 2081.20 &  4  & 37 & 3.14e-06 \\
% 1/500    & 4101.20 &  6  & 36 & 3.14e-06 \\
% 1/1000   & 4880.73 &  6  & 36 & 3.14e-06 \\
% 1/10000  & 3712.00 &  4  & 36 & 3.14e-06 \\
% \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $\nu $&Time& Pic iter&Max Rank (r)& mem(Y)/prod(Y.n) & vort & J & mass \\
% \hline
% 1/10     & 1.1962e+03 &  3  & 50 & 3.0838e-06 & 12.072237403  &   12.074591819  &  1.959282297 \\
% 1/50     & 1.4576e+03 &  3  & 49 & 3.1408e-06 & 12.072165240  &   12.074672861  &  1.959169096 \\
% 1/100    & 1.6360e+03 &  3  & 49 & 3.1408e-06 & 12.072162448  &   12.074675678  &  1.959154368 \\
% 1/500    & 2.0184e+03 &  3  & 49 & 3.1408e-06 & 12.072161554  &   12.074676571  &  1.959142540 \\
% 1/1000   & 2.1673e+03 &  3  & 48 & 3.1408e-06 & 12.072161527  &   12.074676597  &  1.959141060 \\
% 1/5000   & 2.4554e+03 &  3  & 48 & 3.1408e-06 & 12.072161519  &   12.074676604  &  1.959139875 \\
% 1/10000  & 2.5997e+03 &  3  & 47 & 3.1408e-06 & 12.072161519  &   12.074676604  &  1.959139727 \\
% \hline
% \end{tabular}

\begin{figure}[h!]
\caption{CPU time and the number of Picard iterations (left) and TT 
rank, memory reduction ratio and the vorticity (right) 
for different kinematic viscosities.}
\label{fig:nu}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  legend style={at={(0.99,0.01)},anchor=south east},
  xlabel={$-\log_{10}\nu$},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (1,              1.1962e+03 )
   (log10(50.0),    1.4576e+03 )
   (2,              1.6360e+03 )
   (log10(500.0),   2.0184e+03 )
   (3,              2.1673e+03 )
   (log10(5000.0),  2.4554e+03 )
   (4,              2.5997e+03 )}; \addlegendentry{CPU Time};
    % Just add some legend
   \addplot+[] coordinates {(0,0)}; \addlegendentry{Picard iterations};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  ymin=2,ymax=4,
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (1,              3 )
   (log10(50.0),    3 )
   (2,              3 )
   (log10(500.0),   3 )
   (3,              3 )
   (log10(5000.0),  3 )
   (4,              3 )};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  ymin=40,ymax=50,
  legend style={at={(0.01,0.70)},anchor=west},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$-\log_{10}\nu$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (1,              50 )
   (log10(50.0),    49 )
   (2,              49 )
   (log10(500.0),   49 )
   (3,              48 )
   (log10(5000.0),  48 )
   (4,              47 )}; \addlegendentry{max. TT rank};
   % Just for legend
   \pgfplotsset{cycle list shift=1};
   \addplot+[] coordinates {(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$};
   \addplot+[] coordinates {(0,0)}; \addlegendentry{$\|\nabla \times \vv\|^2$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  ymin=-6,ymax=-5,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (1,              3.0838e-06 )
   (log10(50.0),    3.1408e-06 )
   (2,              3.1408e-06 )
   (log10(500.0),   3.1408e-06 )
   (3,              3.1408e-06 )
   (log10(5000.0),  3.1408e-06 )
   (4,              3.1408e-06 )};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  ymin=12,ymax=12.3,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\|\nabla \times \vv\|^2$, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=green!50!black},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=3};
  \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm,color=green!50!black} }
  \addplot+[] coordinates {
   (1,              12.072237403 )
   (log10(50.0),    12.072165240 )
   (2,              12.072162448 )
   (log10(500.0),   12.072161554 )
   (3,              12.072161527 )
   (log10(5000.0),  12.072161519 )
   (4,              12.072161519 )};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

There are two parameters defining the ``randomness'' of the problem: the decay rate of the KLE coefficients $\gamma$, and the total number of KLE components after truncation $m$, see \eqref{eq:kle_in_new}.
The smaller is $\gamma$, the slower is the decay, and hence the larger is the number of ``effective'' dimensions.
This makes the problem more difficult to solve, as we observe in Fig. \ref{fig:gamma}.

% \subsection{Varying stochastic parameter $\gamma$ }

% ~
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for distributed control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
% $1 $     & 9158.50   & 12    &  61     & 13.64    & 6.08e-06  \\
% $2 $     & 1997.50   & 4     &  52     & 12.33    & 4.06e-06  \\
% $3 $     & 1571.00   & 4     &  43     & 12.11    & 3.32e-06  \\
% $4 $     & 1500.70   & 4     &  40     & 12.07    & 3.08e-06  \\
% \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &   mem(Y)/prod(Y.n)  & vort & J & mass \\
% \hline
% $1 $     & 2.5431e+03   & 4    &  91    & 6.0802e-06          & 13.634238558  &   13.637977290  &  1.965385081 \\
% $2 $     & 1.4419e+03   & 3    &  71    & 4.0608e-06          & 12.332533631  &   12.335123842  &  1.960519058 \\
% $3 $     & 1.3330e+03   & 3    &  58    & 3.3245e-06          & 12.116461259  &   12.118856602  &  1.959513328 \\
% $4 $     & 1.1806e+03   & 3    &  50    & 3.0838e-06          & 12.072237403  &   12.074591819  &  1.959282297 \\
% \hline
% \end{tabular}

\begin{figure}[h!]
\caption{CPU time and the vorticity (left) and TT rank and
the memory reduction ratio (right) for different KLE decay rates $\gamma.$}
\label{fig:gamma}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
%   ytick={4,5,6,7,8,9,10},
%   yticklabels={4000,5000,6000,7000,8000,9000,10000},
  xlabel={$\gamma$},
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  ]
  \addplot+[] coordinates {
   (1, 2.5431e+03)
   (2, 1.4419e+03)
   (3, 1.3330e+03)
   (4, 1.1806e+03)}; \addlegendentry{CPU time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\|\nabla \times\vv\|^2$}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\|\nabla \times\vv\|^2$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{log10(12.098133772 - \pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (1, 13.634238558)
   (2, 12.332533631)
   (3, 12.116461259)
   (4, 12.072237403)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$\gamma$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (1, 91)
   (2, 71)
   (3, 58)
   (4, 50)}; \addlegendentry{max. TT rank};
   \pgfplotsset{cycle list shift=1};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (1, 6.0802e-06)
   (2, 4.0608e-06)
   (3, 3.3245e-06)
   (4, 3.0838e-06)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

% \subsection{Varying stochastic parameter $m$ }

Increasing the number of stochastic variables faces the same issue as decreasing the number of spatial grid points: a coarse spatial grid may not resolve the latter KLE terms.
In order to vary $m$ up to $8$, we employ the spatial grid with $h=2^{-4}$.
The results are shown in Fig. \ref{fig:m}.
The CPU times grows milder than linearly with $m$ since (a) the most time-consuming stage is still the solution of the spatial problems, and (b) the model saturates with larger $m$ due to the decaying KLE series.
The latter phenomenon is also illustrated by convergence of the vorticity norm with $m$.

This example illustrates how the TT decomposition gets rid of the curse of dimensionality.
The total number of degrees of freedom in a full tensor would grow exponentially in $m$, while the TT ranks (and hence the TT storage) grow milder than linearly.
This leads to an exponential decay of the memory reduction ratio, which reaches 9 orders of magnitude for the largest dimension.


% ~
%
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $m$ for distributed control with $nc=6,\; \alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $m$ &Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
%  $2$     & 7083.10      &  3    & 31    & 12.10    & 1.41e-04   \\
%  $4$     & 10301.00     &  4    & 42    & 12.10    & 2.91e-06 \\
%  $6$     & 10621.00     &  4    & 51    & 12.10    & 5.34e-08   \\
%  $8$     & 11002.00     &  4    & 49    & 12.12    & 1.08e-09 \\
%  \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $m$ &Time        & Pic iter&Max Rank (Yad) &  mem(Y)/prod(Y.n) & vort & J & mass \\
% \hline
% $1$ & 4.3112e+03 & 3       & 24            & 8.7738e-04        & 12.084563366  &   12.087026021  &  1.960440640 \\
% $2$ & 6.5982e+03 & 3       & 36            & 1.4039e-04        & 12.097289795  &   12.099765617  &  1.960513784 \\
% $3$ & 7.4130e+03 & 3       & 46            & 2.1397e-05        & 12.098013242  &   12.100489680  &  1.960515592 \\
% $4$ & 7.6319e+03 & 3       & 54            & 2.9180e-06        & 12.098107570  &   12.100584075  &  1.960515721 \\
% $5$ & 8.7223e+03 & 3       & 59            & 3.9534e-07        & 12.098126714  &   12.100603229  &  1.960515738 \\
% $6$ & 8.7907e+03 & 3       & 64            & 5.3436e-08        & 12.098131755  &   12.100608273  &  1.960515741 \\
% $7$ & 9.4490e+03 & 3       & 67            & 7.6472e-09        & 12.098133309  &   12.100609827  &  1.960515741 \\
% $8$ & 9.8337e+03 & 3       & 71            & 1.0758e-09        & 12.098133772  &   12.100610291  &  1.960515742 \\
% \hline
% \end{tabular}

\begin{figure}[h!]
\caption{CPU time and the vorticity (left) and TT rank and the memory reduction ratio (right) for different numbers of stochastic variables $m$. Change of defaults: $h=2^{-4}$.}
\label{fig:m}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=8,
  ytick={4,5,6,7,8,9,10},
  yticklabels={4000,5000,6000,7000,8000,9000,10000},
  xlabel={$m$},
  legend style={at={(0.40,0.01)},anchor=south},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  ]
  \addplot+[] coordinates {
   (1, 4.3112e+03)
   (2, 6.5982e+03)
   (3, 7.4130e+03)
   (4, 7.6319e+03)
   (5, 8.7223e+03)
   (6, 8.7907e+03)
   (7, 9.4490e+03)
   (8, 9.8337e+03)}; \addlegendentry{CPU Time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\left(\|\nabla \times\vv_{\star}\|^2 - \|\nabla \times\vv\|^2\right)$}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=8,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\left(\|\nabla \times\vv_{\star}\|^2 - \|\nabla \times\vv\|^2\right)$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (1,    1.3570e-02)
   (2,    8.4398e-04)
   (3,    1.2053e-04)
   (4,    2.6202e-05)
   (5,    7.0580e-06)
   (6,    2.0170e-06)
   (7,    4.6300e-07)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=8,
  legend style={at={(0.50,0.01)},anchor=south},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$m$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (1, 24)
   (2, 36)
   (3, 46)
   (4, 54)
   (5, 59)
   (6, 64)
   (7, 67)
   (8, 71)};
   \pgfplotsset{cycle list shift=1}; \addlegendentry{max. TT rank};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (1, 8.7738e-04)
   (2, 1.4039e-04)
   (3, 2.1397e-05)
   (4, 2.9180e-06)
   (5, 3.9534e-07)
   (6, 5.3436e-08)
   (7, 7.6472e-09)
   (8, 1.0758e-09)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

\subsubsection{Influence of the regularization}
%for fixed $\nu_1=1/10\quad h=1/5, \quad T=30,\quad n_t=2^{10}$
% ~

% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for distributed control with $\alpha=0.$}
% \begin{tabular}{ccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (r)& mem(Y)/prod(Y.n)\\
% \hline
% $10^{3}$  &  2199.64   & 6  &  67  &  8.77e-06 \\
% $10^{2}$  &  2638.25   & 6  &  66  &  8.86e-06  \\
% $10^{1}$  &  2501.33   & 5  &  64  &  8.64e-06  \\
% $10^{0}$  &  1939.22   & 4  &  60  &  7.38e-06  \\
% $10^{-1}$ &  2262.60   & 4  &  51  &  5.40e-06   \\
% $10^{-2}$ &  2581.26   & 6  &  40  &  3.08e-06   \\
% $10^{-3}$ &  1007.56   & 4  &  29  &  1.82e-06   \\
% \hline
% \end{tabular}
% \end{table}

% \begin{tabular}{cccccccc}
% \hline
% $\beta$   &Time   & Pic iter&Max Rank (r)& mem(Y)/prod(Y.n) & vort & J & mass  \\
% \hline
% $10^{-4}$  &  695.6047   & 3  &  37  &  1.2465e-06          & 12.072155180  &   12.072180423  &  1.959117790   \\
% $10^{-3}$  &  809.5868   & 3  &  35  &  1.8179e-06          & 12.072156141  &   12.072406785  &  1.959133178   \\
% $10^{-2}$  &  1.1849e+03 & 3  &  50  &  3.0838e-06          & 12.072237403  &   12.074591819  &  1.959282297    \\
% $10^{-1}$  &  1.7693e+03 & 3  &  84  &  5.4453e-06          & 12.075644954  &   12.092059034  &  1.960514590    \\
% $10^{0}$   &  2.9298e+03 & 4  &  121 &  7.3846e-06          & 12.121485198  &   12.168364586  &  1.967048133    \\
% $10^{1}$   &  4.2312e+03 & 5  &  147 &  8.6410e-06          & 12.245154757  &   12.275779105  &  1.981490687     \\
% $10^{2}$   &  4.0478e+03 & 5  &  154 &  8.8555e-06          & 12.310294779  &   12.315566653  &  1.989307708     \\
% $10^{3}$   &  3.9007e+03 & 5  &  154 &  8.7742e-06          & 12.320554819  &   12.321121492  &  1.990540102     \\
% \hline
% \end{tabular}

In Fig. \ref{fig:beta}, we investigate how the model and the scheme depend on the regularization parameter $\beta$.
Setting $\beta \gg 1$ means solving almost an uncontrolled problem, while $\beta<1$ corresponds to a control with a certain ``power''.
We see that the CPU time and TT ranks grow significantly with $\beta$, when the model switches from effectively Stokes to an essentially nonlinear Navier-Stokes regime.
This is also reflected by a larger number of Picard iterations, needed for the system to converge.
The norm of vorticity starts growing fast when $\beta$ exceeds $1$.
On the other hand, for very large $\beta>10$ we observe a slight decrease of the CPU time due to a smaller number of the local GMRES iterations: the KKT matrix becomes anti-diagonally dominant, and the preconditioner \eqref{eq:sc2} becomes more efficient.

\begin{figure}[h!]
\caption{CPU time and the number of Picard iterations (left) and TT rank, memory reduction ratio and vorticity (right) for different regularization parameters.}
\label{fig:beta}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-4,xmax=3,
  xlabel={$\log_{10}\beta$},
  legend style={at={(0.01,0.99)},anchor=north west},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol
  ]

  \addplot+[] coordinates {
   (-4, 695.6047   )
   (-3, 809.5868   )
   (-2, 1.1849e+03 )
   (-1, 1.7693e+03 )
   (0,  2.9298e+03)
   (1,  4.2312e+03)
   (2,  4.0478e+03)
   (3,  3.9007e+03)}; \addlegendentry{CPU Time};
   \addplot+[] coordinates{(5,0)}; \addlegendentry{Picard iterations}; % from 2 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (-4, 3 )
   (-3, 3 )
   (-2, 3 )
   (-1, 3 )
   (0,  4)
   (1,  5)
   (2,  5)
   (3,  5)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-4,xmax=3,
  legend style={at={(0.99,0.01)},anchor=south east},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$\log_{10}\beta$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (-4, 37  )
   (-3, 35  )
   (-2, 50  )
   (-1, 84  )
   (0,  121)
   (1,  147)
   (2,  154)
   (3,  154)}; \addlegendentry{max. TT rank};
   \pgfplotsset{cycle list shift=1};
   \addplot+[] coordinates{(5,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$}; % from 2 axis
   \addplot+[] coordinates{(5,0)}; \addlegendentry{$\|\nabla \times \vv\|^2$}; % from 3 axis
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol,
  y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (-4, 1.2465e-06 )
   (-3, 1.8179e-06 )
   (-2, 3.0838e-06 )
   (-1, 5.4453e-06 )
   (0,  7.3846e-06)
   (1,  8.6410e-06)
   (2,  8.8555e-06)
   (3,  8.7742e-06)};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\|\nabla \times \vv\|^2$, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=green!50!black},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=3};
  \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm,color=green!50!black} }
  \addplot+[] coordinates {
   (-4, 12.072155180 )
   (-3, 12.072156141 )
   (-2, 12.072237403 )
   (-1, 12.075644954 )
   (0,  12.121485198)
   (1,  12.245154757)
   (2,  12.310294779)
   (3,  12.320554819)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}


Figures \ref{unsteady_Y_time_re_50} and \ref{unsteady_U_time_re_50} show the first two moments of the velocity and control for the final time.
For comparison, the uncontrolled flow is shown in Figure \ref{uncontrolled_Y_re_50}.
We can notice that the uncontrolled flow is more aligned to the top of the domain, while developing an eddy below the step.
The controlled minimal vorticity flow reflects the regime of a larger viscosity.
Moreover, the reattachment point shrinks significantly compared to that in Figure \ref{uncontrolled_Y_re_50}.

\begin{figure}[h!]
\centering
\caption{Plots of the mean (top) and variance (bottom) of the stream function for an unsteady distributed control flow with $\nu = 1/50$ at $t=10$}
%\vspace{-1.5mm}
\includegraphics[scale=0.6]{ns_unsteady_Y_t10_streamlines.pdf}\\
\includegraphics[scale=0.6]{ns_unsteady_varY_t10_streamlines.pdf}
%  \includegraphics[width=0.95\textwidth,height=0.65\textwidth]{ns_unsteady_Y_time_1024_re_50.eps}
\label{unsteady_Y_time_re_50}
 \vspace{-1.5mm}
\end{figure}


\begin{figure}[h!]
\centering
\caption{Plots of the mean (top) and variance (bottom) of the control for an unsteady distributed control flow with $\nu = 1/50,\; n_t=1024.$}
%\vspace{-1.5mm}
 \includegraphics[scale=0.8]{ns_unsteady_U_time_1024_re_50.pdf}
\label{unsteady_U_time_re_50}
 \vspace{-1.5mm}
\end{figure}

% \begin{figure}[h!]
% \centering
% \caption{Plots of the mean (top) and variance (bottom) of the velocity for an unsteady uncontrolled  flow with $\nu = 1/50,\; n_t=3.$}
% %\vspace{-1.5mm}
%  \includegraphics[width=0.95\textwidth,height=0.65\textwidth]{ns_unsteady_uncontrolled_partial_Y_time_3_re_50.eps}
% \label{uncontrolled_Y_re_50}
%  \vspace{-1.5mm}
% \end{figure}


\begin{figure}[h!]
\centering
\caption{Stream function for an unsteady uncontrolled flow at $t=10$.}
%\vspace{-1.5mm}
\includegraphics[scale=0.6]{ns_unsteady_uncontrolled_Y_t10_streamlines.pdf}
%  \includegraphics[width=0.95\textwidth]{ns_unsteady_uncontrolled_time_30_re_50.eps}
\label{uncontrolled_Y_re_50}
 \vspace{-1.5mm}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Time-dependent boundary control problem}
Now we switch to the boundary control, applied as a Neumann boundary condition on the step wall \eqref{kkt_cond1}.
The problem is more difficult since the partial control cannot fully steer the flow to the Stokes regime.
Strongly nonlinear effects that remain in the flow inflate the TT ranks and slow down the convergence.

In particular, in Fig. \ref{fig:nu_bnd} we show how the scheme behaves for different viscosities.
We see that, compared to the distributed control case, even the case $\nu=1/50$ is difficult to solve due to very large TT ranks.
On the other hand, changing the regularization parameter (Fig. \ref{fig:beta_bnd}) does not influence the cost and storage complexities as much as in the distributed control case.
In fact, taking small $\beta$ in the partial control case leads to a poor performance of the Schur complement preconditioner,
and hence a large CPU time due to a large number of local GMRES iterations.
Nevertheless, smaller vorticity norm and number of Picard iterations for $\beta<1$ indicates that the target functional is optimized as expected.

% \subsection{Varying  viscosity}
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for partial control with $n_t=16, \; \alpha=0.$ }
% \begin{tabular}{cccccc}
% \hline
% $\nu $&Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
%  1/10     &  4832.0   &  5    & 86  & 11.44   &1.4e-04 \\
%  1/20     &  10075.0  &  7    & 116 & 11.95   &5.8e-04 \\
%  1/40     &  37246.0  &  10   & 161 & 13.29   & 9.2e-04 \\
%  1/50     &  62334.0  &  11   & 178  & 14.02  & 1.1e-03 \\
% \hline
% \end{tabular}
% \end{table}

\begin{figure}[h!]
\caption{CPU time and the number of Picard iterations (left) and TT rank and memory reduction ratio (right) for different viscosities with the boundary control.}
\label{fig:nu_bnd}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=10,xmax=50,
  xlabel={$1/\nu$},
  ytick={00,10,20,30,40,50,60},
  yticklabels={0,10000,20000,30000,40000,50000,60000},
  ymin=0,ymax=65,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  legend style={at={(0.01,0.99)},anchor=north west},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (10,      4832.0  )
   (20,     10075.0  )
   (40,     37246.0  )
   (50,     62334.0  )}; \addlegendentry{CPU time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{Picard iterations};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (10,     5 )
   (20,     7 )
   (40,     10)
   (50,     11)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=10,xmax=50,
  legend style={at={(0.99,0.01)},anchor=south east},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$1/\nu$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (10,     86 )
   (20,     116)
   (40,     161)
   (50,     178)}; \addlegendentry{max. TT rank};
   \pgfplotsset{cycle list shift=1};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (10,      -3.8539)
   (20,      -3.2366)
   (40,      -3.0362)
   (50,      -2.9586)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

% \subsection{Varying the regularization parameter}
%  %for fixed $\nu_1=1/10\quad h=1/5, \quad T=30,\quad n_t=2^{10}$
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (Yad) &  norm(vort)&mem(Y)/prod(Y.n)\\
% \hline
%  $10^{2}$   &  11135.0   &  7    & 160   & 11.41   & 1.00e-05  \\
%  $10^{1}$   &  10157.0   &  6    & 161   & 11.33   & 1.07e-05  \\
%  $10^{-1}$  &  9884.4    &  5    & 154   & 11.12   & 1.06e-05  \\
%  $10^{-2}$  &  13342.0   &  5    & 160   & 11.11   & 1.06e-05  \\
%  $10^{-3}$  &  19509.0   &  5    & 155   & 11.11   & 1.05e-05   \\
% % \hline
% \end{tabular}
% \end{table}


\begin{figure}[h!]
\caption{CPU time and the number of Picard iterations (left) and TT rank and vorticity (right) for different regularization parameters of the boundary control.}
\label{fig:beta_bnd}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-3,xmax=2,
  xlabel={$\log_{10}\beta$},
  ytick={10,12,14,16,18,20},
  yticklabels={10000,12000,14000,16000,18000,20000},
  ymin=9,ymax=21,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (-3,        19509.0 )
   (-2,        13342.0 )
   (-1,        09884.0 )
   (1 ,        10157.0 )
   (2 ,        11135.0 )}; \addlegendentry{CPU Time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{Picard iterations};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{
   (-3, 5)
   (-2, 5)
   (-1, 5)
   (1 , 6)
   (2 , 7)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-3,xmax=2,
  ymin=80,ymax=180,
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$\log_{10}\beta$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates{
   (-3,    155)
   (-2,    160)
   (-1,    154)
   (1 ,    161)
   (2 ,    160)}; \addlegendentry{max. TT rank};
   \pgfplotsset{cycle list shift=2};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{$\|\nabla \times \vv\|^2$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\|\nabla \times \vv\|^2$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=green!50!black},every y tick label/.style={green!50!black},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=3};
  \addplot+[] coordinates{
   (-3, 11.1100 )
   (-2, 11.1100 )
   (-1, 11.1200 )
   (1 , 11.3300 )
   (2 , 11.4100 )};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

% \subsection{Varying spatial parameter}
% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $nc \;(J)$ &Time& Pic iter&Max Rank (Yad) &  norm(vort)& mem(Y)/prod(Y.n)\\
% \hline
% $3 ()$    & 410.27    &  5  &  65   & 10.39  &9.1e-03  \\
% $4 ()$    & 2623.70   &  5  &  105  & 10.83  &7.6e-04  \\
% $5 ()$    & 13342.0   &  5  &  160  & 11.11  &1.0e-05 \\
% $6 ()$    & 23800.0   &  5  &  95   & 11.59  &4.0e-06  \\
% \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying $N_t$}
% % Three levels as well.
% \begin{table}[h!]
% \centering
% \caption{The temporal parameter $N_t$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $N_t $      & Time  & Pic iter&Max Rank (Yad) &  norm(vort)& mem(Y)/prod(Y.n)\\
% \hline
% % $2^{8}$     &  5211.4  & 5     & 101  & 11.13  & 3.70e-05\\
% $2^{9}$     &  1799.1   &  5    & 94  & 10.84   & 1.3e-03\\
% $2^{10}$    &  13342.0  &  5    & 160 & 11.11   &1.06e-05\\
% % $2^{11}$    &  9333.5   & 6    & 107  & 11.11  & 5.68e-06\\
% $2^{12}$    &  13807.9  & 5    & 125  & 10.83  & 3.27e-04\\
% $2^{13}$    &  23819.0  & 5    & 196  & 11.10  & 2.05e-05\\
% \hline
% \end{tabular}
% \end{table}

% \subsection{Varying stochastic parameter $d$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $d$ for partial control with $nc=6,\; n_t=16,\; \alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $d$ &Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
% $2 $     & 16918.0     &  5   & 63     & 11.59    & 1.9e-02   \\
% $4 $     & 24720.0     &  5   & 94     & 11.59    & 3.9e-04 \\
% $6 $     & 30813.0     &  5   & 106    & 11.59    & 7.0e-06 \\
% $8 $     & 31414.0     &  5   & 119    & 11.59    & 1.2e-07 \\
%  \hline
% \end{tabular}
% \end{table}

Concerning the stochastic parameters, we see qualitatively the same behaviour as
in the distributed control case: all complexity indicators increase for smaller $\gamma$ and larger $m$,
although the CPU time and TT ranks saturate with $m$, see Fig. \ref{fig:gamma_m_bnd}.


\begin{figure}[h!]
\caption{CPU time and the TT rank for different KLE decays (left) and KLE dimensions (right), boundary control.}
\label{fig:gamma_m_bnd}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=1,xmax=4,
  xlabel={$\gamma$},
  ytick={15,20,25,30},
  yticklabels={15000,20000,25000,30000},
  ymin=11,ymax=32,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (1, 31380.0)
   (2, 20274.0)
   (3, 18453.0)
   (4, 13342.0)}; \addlegendentry{CPU Time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{max. TT rank};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=max. TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{
   (1, 285)
   (2, 230)
   (3, 188)
   (4, 160)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=2,xmax=8,
  xlabel={$m$},
  ytick={15,20,25,30},
  yticklabels={15000,20000,25000,30000},
  ymin=11,ymax=32,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  legend style={at={(0.99,0.01)},anchor=south east},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (2, 16918.0)
   (4, 24720.0)
   (6, 30813.0)
   (8, 31414.0)}; \addlegendentry{CPU Time};
   \addplot+[] coordinates{(0,0)}; \addlegendentry{max. TT rank};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=max. TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{
   (2, 63 )
   (4, 94 )
   (6, 106)
   (8, 119)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}

% \subsection{Varying stochastic parameter $\gamma$}
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  norm(vort)& mem(Y)/prod(Y.n)\\
% \hline
% $1 $    & 31380.0  &  5  & 285  & 12.70   & 1.71e-05  \\
% $2 $    & 20274.0  &  5  & 230  & 11.37   & 1.38e-05  \\
% $3 $    & 18453.0  &  5  & 188  & 11.16   & 1.19e-05   \\
% $4 $    & 13342.0  &  5  & 160  & 11.11   & 1.06e-05 \\
% \hline
% \end{tabular}
% \end{table}

%scp -r  ~/navier_stokes_uq/nsuq_partial_control_d_32.mat  pc792.mpi-magdeburg.mpg.de:/scratch/nsuq_partial_control_d_32.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection*{Partial control, time-dependent case with $\alpha=1$}
% 
% \subsection{Varying $N_t$}
% % Three levels as well.
% \begin{table}[h!]
% \centering
% \caption{The temporal parameter $N_t$}
% \begin{tabular}{cccccc}
% \hline
% $N_t $      & Time  & Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2^{8}$     &  2173.85   & 5    & 61  & 11.22  & 3.80e-05\\
% $2^{9}$     &  2372.57   & 5    & 65  & 11.21  & 2.04e-05\\
% $2^{10}$    &  3095.96   & 6    & 70  & 11.21  & 1.09e-05\\
% $2^{11}$    &  3947.89   & 6    & 75  & 11.20  & 5.94e-06\\
% $2^{12}$    &  4022.40  & 6    & 79  & 11.20  & 3.43e-06\\
% $2^{13}$    &  5511.74  & 6    & 85  & 11.20   & 2.12e-06\\
% \hline
% \end{tabular}
% \end{table}
% 
% \subsection{Varying stochastic parameter $\gamma$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for partial control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $1 $    & 4517.60   & 6    & 101  & 12.80   & 1.68e-05  \\
% $2 $    & 4221.70   & 6    & 88   & 11.47   & 1.40e-05  \\
% $3 $    & 3281.50   & 6    & 77   & 11.25   & 1.22e-05   \\
% $4 $    & 1346.10   & 6    & 70   & 11.21   & 1.09e-05  \\
% \hline
% \end{tabular}
% \end{table}
% 
% 
% \subsection{Varying  viscosity}
% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for partial control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\nu $&Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
%  20      & 769.563   &  11   & 22  & 10.88  & 2.31e-06 \\
%  10      & 1433.41   &  20   & 26  & 10.88  & 2.60e-06 \\
%   1      & 809.32    &  4    & 37  & 10.88  & 4.78e-06 \\
%  1/2     & 881.87    &  4    & 43  & 10.88  & 5.89e-06 \\
%  1/10    & 3095.96   &  6    & 70  & 11.21  & 1.09e-05 \\
%  1/20    & 8034.20   &  7    & 94  & 11.98  & 1.62e-05 \\
% \hline
% \end{tabular}
% \end{table}
% 
% 
% \subsection{Varying spatial parameter}
% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for partial control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $nc \;(J)$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2 ()$    & 226.54    &  6  & 28   &  7.94  & 4.76e-05  \\
% $3 ()$    & 310.94    &  6  & 53   & 10.27  & 2.27e-05  \\
% $4 ()$    & 805.60    &  6  & 63   & 10.91  & 1.39e-05  \\
% $5 ()$    & 3095.96   & 6   & 70   & 11.21  & 1.09e-05  \\
% $6 ()$    & 13461.00  &  5  & 74   & 11.36  & 9.63e-06 \\
% \hline
% \end{tabular}
% \end{table}
% 
% \subsection{Varying stochastic parameter $d$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $d$ for partial control with $nc=6,\; \alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $d$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2 $     & 12247.00     &  5   & 64    & 11.37    & 5.19e-04   \\
% $4 $     & 12927.00     &  5   & 75    & 11.37    & 9.63e-06 \\
% $6 $     & 14846.00     &  5   & 78    & 11.37    & 1.71e-07   \\
% $8 $     & 15098.00     &  5   & 83    & 11.37    & 2.81e-09 \\
%  \hline
% \end{tabular}
% \end{table}
% 
% \subsection{Varying the regularization parameter}
% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for  partial control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (Yad) & J &mem(Y)/prod(Y.n)\\
% \hline
% $10^{3}$   &  1657.00    & 7   &  76   & 11.42    & 9.36e-06 \\
% $10^{2}$   &  2020.00    & 7   &  79   & 11.41    & 9.36e-06  \\
% $10^{1}$   &  1979.70    & 7   &  77   & 11.38    & 1.03e-05  \\
% $10^{0}$   &  2193.70    & 7   &  71   & 11.27    & 1.02e-05  \\
% $10^{-1}$  &  2164.40    & 6   &  68   & 11.20    & 1.05e-05  \\
% $10^{-2}$  &  3095.96    & 6   &  70   & 11.21    & 1.09e-05 \\
% $10^{-3}$  &  5333.10    & 6   &  77   & 11.25    & 1.11e-05   \\
% \hline
% \end{tabular}
% \end{table}
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection*{\bf Full control, time-dependent case with $\alpha=1$}
%
% \subsection{Varying $N_t$}
% % Three levels as well.
% \begin{table}[h!]
% \centering
% \caption{The temporal parameter $N_t$  for distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $N_t $      & Time  & Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2^{8}$     &  890.50    & 4    & 45  & 12.11  & 2.45e-05\\
% $2^{9}$     &  962.46    & 4    & 45  & 12.10   & 1.27e-05\\
% $2^{10}$    &  1196.50     & 5    & 53  & 12.10   & 6.81e-06\\
% $2^{11}$    &  1292.67   & 5    & 56  & 12.09  & 3.06e-06\\
% $2^{12}$    &  1639.60   & 5    & 56  & 12.09   & 1.99e-06\\
% $2^{13}$    &  2230.21   & 5    & 56  & 12.09   &  1.18e-06\\
%
% \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying  viscosity}
% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\nu $&Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
%  20       & 450.17    &  3    & 19  & 12.07   & 2.31e-06 \\
%   10      & 560.99    &  3    & 21  & 12.07   & 2.60e-06 \\
%    1      & 641.02    &  4    & 37  & 12.07   & 3.87e-06 \\
%   1/2     & 875.36    &  4    & 42  & 12.07   & 4.48e-06 \\
%  1/10     & 1196.50     &  5    & 53  & 12.10   & 6.81e-06\\
%  1/20     & 1507.94   &  5    & 60  & 12.13   & 8.07e-06 \\
% \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying spatial parameter}
% Three levels from ifiss starting from smallest
% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $nc \;(J)$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2 ()$    & 50.87     &  3  & 12   &  9.93  & 1.27e-05  \\
% $3 ()$    & 90.19     &  4  & 28   & 11.91  & 1.03e-05  \\
% $4 ()$    & 216.87    &  4  & 39   & 12.02  & 8.19e-06  \\
% $5 ()$    & 1196.50     &  5  & 53   & 12.10   & 6.81e-06  \\
% $6 ()$    & 6658.40   &  5  & 54   & 12.11  & 6.11e-06 \\
% \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying stochastic parameter $d$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $d$ for distributed control with $nc=6,\; \alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $d$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
%  $2 $     & 5055.10     &  4   & 43    & 12.11    & 3.12e-04   \\
%  $4 $     & 6548.40     &  5   & 54    & 12.12    & 6.11e-06 \\
%  $6 $     & 6056.10     &  4   & 51    & 12.12    & 1.07e-07   \\
%  $8 $     & 7135.40     &  5   & 64    & 12.12    & 1.86e-09 \\
%  \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying the regularization parameter}
% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for  distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (Yad) & J &mem(Y)/prod(Y.n)\\
% \hline
% $10^{3}$   &  1315.20    & 7  &  76   & 12.32    & 8.45e-06 \\
% $10^{2}$   &  1545.40    & 7  &  76   & 12.32    & 8.53e-06  \\
% $10^{1}$   &  1831.90    & 7   &  75   & 12.31    & 8.85e-06  \\
% $10^{0}$   &  1357.50    & 6   &  72   & 12.27    & 8.84e-06  \\
% $10^{-1}$  &  1352.80    & 5   &  62   & 12.17    & 7.89e-06  \\
% $10^{-2}$  &  1196.50    & 5   &  53   & 12.10    & 6.81e-06 \\
% $10^{-3}$  &  939.0112   & 5   &  45   & 12.39    & 5.68e-06   \\
% \hline
% \end{tabular}
% \end{table}
%
%
% \subsection{Varying stochastic parameter $\gamma$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $1 $     & 1943.80   & 5     &  80     & 13.67    & 1.15e-05  \\
% $2 $     & 1384.50   & 5     &  69     & 12.34    & 9.42e-06  \\
% $3 $     & 1249.70   & 5     &  60     & 12.14    & 7.83e-06  \\
% $4 $     & 1196.50   & 5     &  53     & 12.10    & 6.81e-06  \\
% \hline
% \end{tabular}
% \end{table}
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% \subsection*{Distributed control, steady-state case with $\alpha=1$}
%
% \subsection{Varying  viscosity}
% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for distributed control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\nu $&Time& Pic iter     &Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
%   1/10      &  561.95    &  2    & 23   & 12.72    & 1.4e-03\\
%   1/20      &  672.13    &  3    & 25   & 12.76    & 1.6e-03\\
%   1/40      & 1110.90    &  3    & 26   & 12.85    & 2.0e-03\\
%   1/50      & 1141.10    &  3    & 28    & 12.90   & 2.1e-03 \\
%
%  \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying the regularization parameter}
% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for distributed  control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (Yad) & J &mem(Y)/prod(Y.n)\\
% \hline
%  %$10^{3}$   &  494.95    & 7   &  20   & 11.95    & 2.7e-06 \\
%  $10^{2}$   &  647.70    & 4   &  24   & 12.96    & 1.4e-03  \\
%  $10^{1}$   &  629.75   & 4   &  23   & 11.92     & 1.4e-03  \\
%  $10^{-1}$  &   614.68    & 3   &  24   & 12.80   & 1.4e-03  \\
%  $10^{-2}$  &  561.95    & 2   &  23   & 12.72    & 1.4e-03  \\
%  $10^{-3}$  &  633.25    & 3   &  22   & 13.03    & 1.4e-03 \\
% \hline
% \end{tabular}
% \end{table}
%
%
%
%
%
% \subsection{Varying spatial parameter}
% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for distributed  control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $nc \;(J)$ & Time  & Pic iter& Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $3 ()$    & 12.40     &  3    & 7    & 12.36   & 2.5e-01  \\
% $4 ()$    & 84.50     &  2    & 13   & 12.64   & 5.5e-02  \\
% $5 ()$    & 561.95    &  2    & 23   & 12.72   & 1.4e-03\\
% $6 ()$    & 4156.40   &  2    & 23   & 12.74   & 1.3e-03 \\
% \hline
% \end{tabular}
% \end{table}
%
%
% \subsection{Varying stochastic parameter $\gamma$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for distributed  control with $\alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $1 $    & 1095.0   & 3    & 50   & 14.37  & 3.1e-03  \\
% $2 $    &  623.6   & 2    & 35   & 12.99  & 2.2e-03  \\
% $3 $    &  590.7   & 2    & 26   & 12.76  & 1.7e-03   \\
% $4 $    &  561.95  &  2   & 23   & 12.72  & 1.4e-03  \\
% \hline
% \end{tabular}
% \end{table}
%
% \subsection{Varying stochastic parameter $d$ }
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $d$ for distributed  control with $nc=6,\; \alpha=1.$}
% \begin{tabular}{cccccc}
% \hline
% $d$ &Time& Pic iter&Max Rank (Yad) &  J & mem(Y)/prod(Y.n)\\
% \hline
% $2 $     & 2724.0     &  2   & 14    & 12.74    & 5.47e-02   \\
% $4 $     & 3131.7     &  2   & 23    & 12.74    & 1.30e-03 \\
% $6 $     & 3290.3     &  2   & 31    & 12.75    & 2.68e-05   \\
% $8 $     & 4799.0     &  2   & 32    & 12.74    & 4.80e-07 \\
%  \hline
% \end{tabular}
% \end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Stationary boundary control problem}
One of the sources of the rank inflation in the previous section is the simultaneous storage of
all time snapshots in a single TT representation.
In the next tests we consider the stationary Navier-Stokes equations as constraints.
In Fig. \ref{fig:nu_stat}, we vary the viscosity, and in Fig. \ref{fig:beta_stat} we investigate the influence of the regularization parameter.
We see that the results reflect qualitatively the behaviour of the time-dependent problem with the
boundary control, i.e. the stochastic components are also strongly coupled due to the nonlinearity, 
and it needs more Picard iterations and larger TT ranks for the solution to converge in a low viscosity regime.
Nonetheless, the TT ranks (and hence the computational times) are much smaller than in the time-dependent case.


% \subsection{Varying  viscosity}
%
% ~

% \begin{table}[h!]
% \centering
% \caption{Varying the values of $\nu$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\nu $&Time& Pic iter     &Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
%   1/10       &  1286.10    &  5     & 33  & 11.69    & 1.5e-03  \\
%   1/20       &  2329.90    &  7     & 39  & 12.25    & 1.9e-03 \\
%   1/40       &  6997.20    &  10     & 47  & 13.75    & 2.6e-03 \\
%   1/50       & 10261.00    &  12     & 52  & 14.60    & 2.8e-03 \\
%  \hline
% \end{tabular}
% \end{table}

\begin{figure}[h!]
\caption{CPU time and number of Picard iterations (left) and TT rank and memory reduction ratio (right) for different viscosities, boundary control, stationary problem.}
\label{fig:nu_stat}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=10,xmax=50,
  xlabel={$1/\nu$},
  ytick={2,4,6,8,10},
  yticklabels={2000,4000,6000,8000,10000},
  ymin=0,ymax=11,
  y filter/.code={\pgfmathparse{0.001*\pgfmathresult}\pgfmathresult},
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (10,    1286.10)
   (20,    2329.90)
   (40,    6997.20 )
   (50,   10261.00 )}; \addlegendentry{CPU Time};
  \addplot+[] coordinates{(0,0)}; \addlegendentry{Picard iterations};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates {
   (10,      5)
   (20,      7)
   (40,     10)
   (50,     12)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=10,xmax=50,
  ymin=29,ymax=54,
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$1/\nu$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates {
   (10,     33)
   (20,     39)
   (40,     47)
   (50,     52)}; \addlegendentry{max. TT rank};
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{(0,0)}; \addlegendentry{$\log_{10}\mbox{MemR}$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\log_{10}\mbox{MemR}$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=red},every y tick label/.style={red},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates {
   (10,   -2.8239)
   (20,   -2.7212)
   (40,   -2.5850)
   (50,   -2.5528)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}
%
% \subsection{Varying the regularization parameter}
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{The regularization parameter $\beta$ for  partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\beta$&Time& Pic iter&Max Rank (Yad) & norm(vort) &mem(Y)/prod(Y.n)\\
% \hline
% %$10^{3}$   &  853.96    & 7   &  28   & 11.96    & 2.7e-03 \\
% $10^{2}$   &  1259.60   & 7   &  30   & 12.03    & 1.4e-03  \\
% $10^{1}$   &  1325.30   & 5   &  30   & 11.95    & 1.4e-03  \\
% $10^{-1}$  &  1240.10   & 5   &  30   & 11.70    & 1.5e-03  \\
% $10^{-2}$  &  1286.10   & 5   &  33   & 11.69    & 1.5e-03 \\
% $10^{-3}$  &  1392.50   & 5   &  33   & 11.68    & 1.5e-03   \\
% \hline
% \end{tabular}
% \end{table}

\begin{figure}[h!]
\caption{CPU time and number of Picard iterations (left) and TT rank and vorticity (right) for different regularization parameters, boundary control, stationary problem.}
\label{fig:beta_stat}
\resizebox{0.98\linewidth}{!}{
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-3,xmax=2,
  xlabel={$\log_{10}\beta$},
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
  cycle list name=stokescontrol,
  ]
  \addplot+[] coordinates {
   (-3,    1392.5 )
   (-2,    1286.1 )
   (-1,    1240.1 )
   (1 ,    1325.3 )
   (2 ,    1259.6 )}; \addlegendentry{CPU Time};
  \addplot+[] coordinates{(5,0)}; \addlegendentry{Picard iterations};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=Picard iterations, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},every y tick label/.style={blue},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=1};
  \addplot+[] coordinates{
   (-3, 5)
   (-2, 5)
   (-1, 5)
   (1 , 5)
   (2 , 7)};
  \end{axis}
\end{tikzpicture}
\begin{tikzpicture}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  xmin=-3,xmax=2,
  ymin=29,ymax=54,
  legend style={at={(0.50,0.99)},anchor=north},
  ylabel=max. TT rank, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90},
  xlabel={$\log_{10}\beta$},
  cycle list name=stokescontrol
  ]
  \addplot+[] coordinates{
   (-3,    33)
   (-2,    33)
   (-1,    30)
   (1 ,    30)
   (2 ,    30)}; \addlegendentry{max. TT rank};
  \pgfplotsset{cycle list shift=2};
  \addplot+[] coordinates{(0,0)}; \addlegendentry{$\|\nabla \times \vv\|^2$};
  \end{axis}
  \begin{axis}[%
  xmode=normal,
  ymode=normal,
  legend style={at={(0.01,0.01)},anchor=south west},
  ylabel=$\|\nabla \times \vv\|^2$, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=green!50!black},every y tick label/.style={green!50!black},
  axis x line=none,
  cycle list name=stokescontrol
  ]
  \pgfplotsset{cycle list shift=3};
  \addplot+[] coordinates{
   (-3,    11.6800)
   (-2,    11.6900)
   (-1,    11.7000)
   (1 ,    11.9500)
   (2 ,    12.0300)};
  \end{axis}
\end{tikzpicture}
}
\end{figure}


Finally,  Figures \ref{stationary_partial_Y_re_50} and \ref{stationary_partial_U_re_50},
depict the moments of the velocity and control of the stationary boundary controlled flow.
The boundary controlled flow still manifests an eddy around the corner.
Nonetheless, it corresponds to a smaller vorticity than in the uncontrolled
regime (Fig. \ref{uncontrolled_Y_re_50}), where there is a rather sharp 
interface between the main stream and the corner area.
We also remark that the number of Picard iterations remains generally below 10, with just a few exceptions.

\begin{figure}[h!]
\centering
\caption{Plots of the mean (top) and variance (bottom) of the stream function for a stationary boundary control   flow with $\nu = 1/50.$}
%\vspace{-1.5mm}
\includegraphics[scale=0.6]{ns_stationary_partial_Y_streamlines.pdf} \\
\includegraphics[scale=0.6]{ns_stationary_partial_varY_streamlines.pdf}
%  \includegraphics[width=0.95\textwidth,height=0.65\textwidth]{ns_stationary_partial_Y_re_50.eps}
\label{stationary_partial_Y_re_50}
 \vspace{-1.5mm}
\end{figure}


\begin{figure}[h!]
\centering
\caption{Plots of the mean (top) and variance (bottom) of control for a  stationary boundary control  flow with $\nu = 1/50.$}
%\vspace{-1.5mm}
 \includegraphics[scale=0.8]{ns_stationary_partial_U_re_50.pdf}
\label{stationary_partial_U_re_50}
 \vspace{-1.5mm}
\end{figure}



% \begin{figure}
% \centering
% \caption{Plots of the mean (top) and variance (bottom) of the control for an uncontrolled unsteady flow with $\nu = 1/50,\; n_=16.$}
% %\vspace{-1.5mm} 
%  \includegraphics[width=0.95\textwidth,height=0.75\textwidth]{ns_unsteady_uncontrolled_partial_U_time_16_re_50.eps}
% \label{uncontrolled_U_re_50}
%  \vspace{-1.5mm} 
% \end{figure}

% \begin{figure}[h!]
% \centering
% \caption{Conceptual model of a fractured vuggy reservoir at multiple scales \cite{PEQ09}.}
% \vspace{-3.5cm} 
% \includegraphics[width=0.98\textwidth,height=0.98\textwidth]{example_SB}%{oil_reserve_SB.pdf}
%  \label{sb_example}
%  \vspace{-3.2cm} 
% \end{figure}



% \subsection{Varying spatial parameter}
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{The spatial parameter $h$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $nc \;(J)$ & Time  & Pic iter& Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
%  $3 ()$    & 20.78     &  5     & 8    & 10.94    & 2.5e-01  \\
%  $4 ()$    & 207.83    &  5     & 27   & 11.40    & 5.5e-02   \\
%  $5 ()$    & 1286.10   &  5     & 33   & 11.69    & 1.5e-03 \\
%  $6 ()$    & 9188.30   &  5     & 26   & 11.85    & 1.5e-03 \\
% \hline
% \end{tabular}
% \end{table}

%
% \subsection{Varying stochastic parameter $\gamma$ }
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $\gamma$ for partial control with $\alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $\gamma$ &Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
% $1 $    &  3971.90   & 5    & 85   & 13.36  & 3.4e-03  \\
% $2 $    &  1964.40   & 5    & 58   & 11.97  & 2.5e-03  \\
% $3 $    &  1517.30   & 5    & 43   & 11.73  & 1.7e-03   \\
% $4 $    &  1286.10   & 5    & 33   & 11.69  & 1.5e-03  \\
% \hline
% \end{tabular}
% \end{table}
%
%
% \subsection{Varying stochastic parameter $m$ }
%
% ~
%
% \begin{table}[h!]
% \centering
% \caption{The stochastic parameter $m$ for partial control with $nc=6,\; \alpha=0.$}
% \begin{tabular}{cccccc}
% \hline
% $m$ &Time& Pic iter&Max Rank (Yad) &  norm(vort) & mem(Y)/prod(Y.n)\\
% \hline
% $2 $     & 4810.10    &  5   & 15    & 11.85    & 5.7e-02   \\
% $4 $     & 6913.30    &  5   & 35    & 11.85    & 1.5e-03 \\
% $6 $     & 8777.80    &  5   & 47    & 11.85    & 2.6e-05   \\
% $8 $     & 9918.00    &  5   & 55    & 11.83    & 5.1e-07 \\
%  \hline
% \end{tabular}
% \end{table}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% \begin{figure}[h!]
% \centering
% \caption{Performance for different regularizations (left) and system times (right)}
% \label{fig:beta_T}
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={-3,-2,-1,0,1,2,3},
%   ymode=normal,
%   xlabel=$\log_{10}\beta$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=$\log_{10}\mbox{CPU Time}$, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
%   ]
%   \addplot+[] coordinates {(-3,9169.80) (-2,2263.80) (-1,1210.17) (1,887.50) (2,1027.17) (3,942.25)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={-3,-2,-1,0,1,2,3},
%   ymode=normal,ymin=40,ymax=76,
%   xlabel=$\log_{10}\beta$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(-3,61) (-2,60) (-1,57) (1,64) (2,63) (3,59)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={-3,-2,-1,0,1,2,3},
%   ymode=normal,ymin=0.9,ymax=1.32,
%   xlabel=$\log_{10}\beta$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(-3,1.20) (-2,1.18) (-1,1.11) (1,1.26) (2,1.24) (3,1.15)};
%   \end{axis}
% \end{tikzpicture}}
% %
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={20,40,60,80,100},
%   ymode=normal,
%   xlabel=$T$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol
%   ]
%   \addplot+[] coordinates {(30,1210.17) (50,1162.50) (80,945.84) (100,774.87)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={20,40,60,80,100},
%   ymode=normal,ymin=40,ymax=64,
%   xlabel=$T$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(30,57) (50,51) (80,45) (100,42)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,xtick={20,40,60,80,100},
%   ymode=normal,ymin=0.7,ymax=1.25,
%   xlabel=$T$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(30,1.11) (50,0.98) (80,0.85) (100,0.79)};
%   \end{axis}
% \end{tikzpicture}}
% \end{figure}
%
%
%
% \begin{figure}[h!]
% \centering
% \caption{Performance for different spatial (left) and temporal (right) grid sizes}
% \label{fig:h_nt}
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={2,3,4,5},
%   ymode=normal,
%   xlabel=$log_{2}(1/h)$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=$\log_{10}\mbox{CPU Time}$, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{log10(\pgfmathresult)}\pgfmathresult},
%   ]
%   \addplot+[] coordinates {(2,1210.17) (3,3221.73) (4,15004.20) (5,63520.00)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=20,ymax=82,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(2,57) (3,56) (4,54) (5,55)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=0.4,ymax=1.6,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(2,1.11) (3,0.79) (4,0.69) (5,0.68)};
%   \end{axis}
% \end{tikzpicture}}
% %
% \resizebox{0.48\linewidth}{!}{
% \begin{tikzpicture}
%   \begin{axis}[%
%   xmode=normal,xtick={8,9,10,11,12},
%   ymode=normal,
%   xlabel=$log_{2}N_t$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=CPU Time, axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol
%   ]
%   \addplot+[] coordinates {(8,488.18) (9,801.18) (10,1210.17) (11,2741.02) (12,5459.30)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=20,ymax=82,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=TT rank, axis y line*=right, y label style={at={(1.0,1.0)},anchor=south east,rotate=-90,color=blue},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(8,46) (9,51) (10,57) (11,60) (12,61)};
%   \end{axis}
%   \begin{axis}[%
%   xmode=normal,
%   ymode=normal,ymin=0.4,ymax=2.9,
%   xlabel=$log_{10}\varepsilon$,
%   legend style={at={(0.01,0.01)},anchor=south west},
%   ylabel=\% Mem, axis y line*=right, y label style={at={(1.25,1.0)},anchor=south east,rotate=-90,color=red},
%   axis x line=none,
%   cycle list name=stokescontrol
%   ]
%   \pgfplotsset{cycle list shift=2};
%   \pgfplotsset{every outer y axis line/.style={xshift=1cm}, every tick/.style={xshift=1cm}, every y tick label/.style={xshift=1cm} }
%   \addplot+[] coordinates {(8,2.81) (9,1.71) (10,1.11) (11,0.74) (12,0.53)};
%   \end{axis}
% \end{tikzpicture}}
% \end{figure}


% \begin{figure}[h!]
% \centering
% \caption{Top: squared vorticity and the total cost functional $\mathcal{J}$. Bottom: streamlines of the flow with optimized vorticity.}
% \label{fig:J}
% \begin{minipage}[t]{0.75\linewidth}
% \resizebox{\linewidth}{!}{
% \begin{tikzpicture}[yscale=0.75]
%   \begin{axis}[%
%   xmode=normal,xtick={-3,-2,-1,0,1,2,3,4,5},
%   ymode=normal,
%   xlabel=$log_{10}\beta$,
%   legend style={at={(0.01,0.99)},anchor=north west},
%   axis y line*=left, y label style={at={(-0.01,1.0)},anchor=south west,rotate=-90,color=black},
%   cycle list name=stokescontrol,
%   y filter/.code={\pgfmathparse{\pgfmathresult/2}\pgfmathresult},
%   ]
%   \addplot+[] coordinates {(-3,2.2017846e+01) (-2,2.2017358e+01) (-1,2.2028621e+01) (0,2.2234707e+01) (1,2.4400752e+01) (2,2.7861104e+01) (3,2.8909062e+01) (4,2.9041568e+01) (5,2.9055199e+01)};\addlegendentry{\LARGE $\frac{1}{2}\|\nabla \times\vv\|^2$};
%   \pgfplotsset{cycle list shift=1};
%   \addplot+[] coordinates {(-3,2.2019395e+01) (-2,2.2031279e+01) (-1,2.2143688e+01) (0,2.2946561e+01) (1,2.5742086e+01) (2,2.8326912e+01) (3,2.8969893e+01) (4,2.9047855e+01) (5,2.9055830e+01)};\addlegendentry{\LARGE $\mathcal{J}$};
%   \end{axis}
% \end{tikzpicture}}
% \end{minipage} \\
%
% \includegraphics[width=\linewidth]{streamlines.pdf}
% \end{figure}
%


\section{Conclusions and outlook}
\label{con_out}
We demonstrated the applicability of low-rank tensor decompositions to the
solution of optimal control problems constrained by unsteady Navier-Stokes equations with stochastic inputs.
This problem has a threefold challenge: a nonlinear time-dependent PDE, an optimization problem using a Lagrangian approach, and random inputs.
In a classical stochastic Galerkin discretization one needs to multiply the numbers of degrees of freedom coming from space, time and stochastic quantities.
In particular, the outer optimization implies storing all time snapshots, and random quantities are introduced as independent variables.

We never store or compute all elements of such a solution as this would quickly become infeasible even on large computers.
Instead we compress them into the Tensor Train representation.
A crucial part of the scheme is an alternating iterative solver,
which computes the TT factors directly.
Although  known in the multilinear algebra community, this idea required substantial modifications in order to be applicable to the OPNS.
We preserve the saddle-point structure in the reduced model
and accommodate components of different sizes, such as the boundary control.
Moreover, given a low-rank representation of the previous solution, we assemble the linearized operator also in a low-rank form.
The scheme provides a significant reduction of complexity, up to 9 orders of magnitude for the largest number of random variables.

% In the future, it would be interesting to consider additional parameters, coming from, for instance,
% uncertain geometry.
% As such considerations would quite likely lead to increased complexity due to larger TT ranks,
% further improvements to the scheme are desirable. Besides, these problems are difficult for other approaches,
% especially those relying on the full representation of the solution.


\bibliographystyle{siam}
\bibliography{refs}
\end{document}

paper_arxiv.bbl
\begin{thebibliography}{10}

\bibitem{AHU16}
{\sc A.~Ali, M.~Hinze, and E.~Ullmann}, {\em Multilevel {M}onte {C}arlo
  analysis for optimal control of elliptic {PDE}s with random coefficients},
  arXiv:1602.08943,  (2016).

\bibitem{babuska-hp-1988}
{\sc I.~Babu\u{s}ka and B.~Q. Guo}, {\em The hp version of the finite element
  method for domains with curved boundaries}, SIAM Journal on Numerical
  Analysis, 25 (1988), pp.~837--861.

\bibitem{babuska-regularity-1988}
\leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em Regularity of the
  solution of elliptic problems with piecewise analytic data. {P}art {I}.
  {B}oundary value problems for linear elliptic equation of second order}, SIAM
  Journal on Mathematical Analysis, 19 (1988), pp.~172--203.

\bibitem{BTZ04}
{\sc I.~Babu\u{s}ka, R.~Tempone, and G.~Zouraris}, {\em Galerkin finite element
  approximations of stochastic elliptic partial differential equations}, SIAM
  Journal on Numerical Analysis, 42 (2004), pp.~800--825.

\bibitem{BSOS215}
{\sc P.~Benner, S.~Dolgov, A.~Onwunta, and M.~Stoll}, {\em {Low-rank solvers
  for unsteady Stokes-Brinkman optimal control problem with random data}},
  Computer Methods in Applied Mechanics and Engineering, 304 (2016),
  pp.~26--54.

\bibitem{BH14}
{\sc P.~Benner and M.~W. Hess}, {\em {Reduced Basis Modeling for Uncertainty
  Quantification of Electromagnetic Problems in Stochastically Varying
  Domain}}, in Scientific Computing in Electrical Engineering, A.~Bartel,
  M.~Clemens, M.~Guenther, and E.~J.~W. ter Maten, eds., Springer, 2016,
  pp.~215 -- 222.

\bibitem{BOS14}
{\sc P.~Benner, A.~Onwunta, and M.~Stoll}, {\em {Low-rank solution of unsteady
  diffusion equations with stochastic coefficients}}, SIAM/ASA Journal on
  Uncertainty Quantification, 3 (2015), pp.~622 -- 649.

\bibitem{BOS15a}
\leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em {Block-diagonal
  preconditioning for optimal control problems constrained by {PDE}s with
  uncertain inputs}}, SIAM Journal on Matrix Analysis and Applications, 37
  (2016), pp.~491 -- 518.

\bibitem{BenGolLie05}
{\sc M.~Benzi, G.~H. Golub, and J.~Liesen}, {\em Numerical solution of saddle
  point problems}, Acta Numer, 14 (2005), pp.~1--137.

\bibitem{BPS13}
{\sc A.~Bespalov, C.~E. Powell, and D.~Silvester}, {\em Energy norm a
  posteriori error estimation for parametric operator equations}, SIAM Journal
  on Scientific Computing, 36 (2013), pp.~A339 -- A363.

\bibitem{BW09}
{\sc A.~Borzi and G.~von Winckel}, {\em Multigrid methods and sparse-grid
  collocation techniques for parabolic optimal control problems with random
  coefficients}, SIAM Journal on Scientific Computing, 31 (2009), pp.~2172 --
  2192.

\bibitem{Cheney}
{\sc W.~Cheney and D.~Kincaid}, {\em Numerical Mathematics and Computing},
  Brooks/Cole Publishing Co., Pacific Grove, CA, USA, 6th~ed., 2007.

\bibitem{ChoLu2015}
{\sc A.~J. Chorin and F.~Lu}, {\em Discrete approach to stochastic
  parametrization and dimension reduction in nonlinear dynamics}, PNAS, 112
  (2015), pp.~9804--9809.

\bibitem{SSNS17}
{\sc S.~Dolgov and M.~Stoll}, {\em Low-rank solution to an optimization problem
  constrained by the {N}avier--{S}tokes equations}, SIAM Journal on Scientific
  Computing, 39 (2017), pp.~A255--A280.

\bibitem{dkos-eigb-2014}
{\sc S.~V. Dolgov, B.~N. Khoromskij, I.~V. Oseledets, and D.~V. Savostyanov},
  {\em Computation of extreme eigenvalues in higher dimensions using block
  tensor train format}, Computer Physics Communications, 185 (2014),
  pp.~1207--1216.

\bibitem{ds-amen-2014}
{\sc S.~V. Dolgov and D.~V. Savostyanov}, {\em Alternating minimal energy
  methods for linear systems in higher dimensions}, SIAM Journal on Scientific
  Computing, 36 (2014), pp.~A2248--A2271.

\bibitem{EEOS05}
{\sc H.~Elman, O.~G. Ernst, D.~P. O'Leary, and M.~Stewart}, {\em Efficient
  iterative algorithms for the stochastic finite element method with
  applications to acoustic scattering}, Computer Methods in Applied Mechanics
  and Engineering, 194 (2005), pp.~1037--1055.

\bibitem{EF07}
{\sc H.~Elman and D.~Furnival}, {\em Solving steady-state diffusion problem
  using multigrid}, IMA Journal of Numerical Analysis, 27 (2007), pp.~675--688.

\bibitem{ESW14}
{\sc H.~Elman, D.~Silvester, and A.~Wathen}, {\em Finite {E}lements and {F}ast
  {I}terative {S}olvers}, vol.~Second Edition, Oxford University Press, 2014.

\bibitem{GS91}
{\sc R.~G. Ghanem and P.~Spanos}, {\em Stochastic Finite Elements: A Spectral
  Approach}, Springer-Verlag, New York, 1996.

\bibitem{GWZ14}
{\sc M.~D. Gunzburger, C.~G. Webster, and G.~Zhang}, {\em Stochastic finite
  element methods for equations with random input data}, Acta Numerica, 23
  (2014), pp.~521--650.

\bibitem{hackbusch-2012}
{\sc W.~Hackbusch}, {\em Tensor Spaces And Numerical Tensor Calculus},
  Springer--Verlag, Berlin, 2012.

\bibitem{HS10}
{\sc R.~Herzog and E.~Sachs}, {\em {Preconditioned conjugate gradient method
  for optimal control problems with control and state constraints}}, SIAM
  Journal on Matrix Analysis and Applications, 31 (2010), pp.~2291 -- 2317.

\bibitem{HIK02}
{\sc M.~Hinterm{\"{u}}ller, K.~Ito, and K.~Kunisch}, {\em {The primal-dual
  active set strategy as a semi-smooth {N}ewton method}}, SIAM Journal on
  Optimization, 13 (2002), pp.~865 -- 888.

\bibitem{Hinze2000}
{\sc Michael Hinze}, {\em Optimal and instantaneous control of the instationary
  Navier-Stokes equations}, Habilitation, TU Berlin, 2000.

\bibitem{holtz-ALS-DMRG-2012}
{\sc S.~Holtz, T.~Rohwedder, and R.~Schneider}, {\em The alternating linear
  scheme for tensor optimization in the tensor train format}, SIAM Journal on
  Scientific Computing, 34 (2012), pp.~A683--A713.

\bibitem{IK03}
{\sc K.~Ito and K.~Kunisch}, {\em {Semi-smooth {N}ewton methods for
  state-constrained optimal control problems}}, Systems and Control Letters, 50
  (2003), pp.~221 -- 228.

\bibitem{jeckelmann-dmrgsolve-2002}
{\sc E.~Jeckelmann}, {\em Dynamical density matrix renormalization group
  method}, Physical Review B, 66 (2002), p.~045114.

\bibitem{Kanzow04}
{\sc Ch. Kanzow}, {\em {Inexact semi-smooth {N}ewton methods for large-scale
  complementarity problems}}, Optimization Methods and Software, 19 (2004),
  pp.~309 -- 325.

\bibitem{bokh-surv-2014}
{\sc B.~N. Khoromskij}, {\em Tensor numerical methods for multidimensional
  {PDE}s: theoretical analysis and initial applications}, ESAIM Proceedings, 48
  (2015), pp.~1--28.

\bibitem{KHRW13}
{\sc D.~P. Kouri, M.~Heinkenschloss, D.~Ridzal, and B.~G. van
  Bloemen~Waanders}, {\em A trust-region algorithm with adaptive stochastic
  collocation for {PDE} optimization under uncertainty}, SIAM Journal on
  Scientific Computing, 35 (2013), pp.~A1847 -- A1879.

\bibitem{LKDNG03}
{\sc O.~P. {Le Ma\^{\i}tre}, O.~M. Knio, B.~J. Debusschere, H.~N. Najm, and
  R.~G. Ghanem}, {\em A multigrid solver for two-dimensional stochastic
  diffusion equations}, Computer Methods in Applied Mechanics and Engineering,
  192 (2003), pp.~4723--4744.

\bibitem{NT09}
{\sc F.~Nobile and R.~Tempone}, {\em Analysis and implementation issues for the
  numerical approximation of parabolic equations with random coefficients},
  International Journal for Numerical Methods in Engineering, 80 (2009),
  pp.~979 -- 1006.

\bibitem{osel-tt-2011}
{\sc I.~V. Oseledets}, {\em Tensor-train decomposition}, SIAM Journal on
  Scientific Computing, 33 (2011), pp.~2295--2317.

\bibitem{tt-toolbox}
{\sc I.~V. Oseledets, S.~Dolgov, V.~Kazeev, D.~Savostyanov, O.~Lebedeva,
  P.~Zhlobich, T.~Mach, and L.~Song}, {\em {TT-Toolbox}}, 2016.
\newblock https://github.com/oseledets/TT-Toolbox.

\bibitem{SPW11}
{\sc J.~W. Pearson, M.~Stoll, and A.~Wathen}, {\em {Preconditioners for state
  constrained optimal control problems with {M}oreau-{Y}osida penalty
  function}}, Numerical Linear Algebra with Applications,  (2011), pp.~81 --
  97.

\bibitem{PerezGarcia-mps-2007}
{\sc D.~Perez-Garcia, F.~Verstraete, M.~M. Wolf, and J.~I. Cirac}, {\em Matrix
  product state representations}, Quantum Info. Comput., 7 (2007),
  pp.~401--430.

\bibitem{PST15}
{\sc M.~Porcelli, V.~Simoncini, and M.~Tani}, {\em {Preconditioning of
  active-set {N}ewton methods for {PDE}-constrained optimal control problems}},
  SIAM Journal on Scientific Computing, 37 (2015), pp.~S472 -- S502.

\bibitem{PE09}
{\sc C.~E. Powell and H.~Elman}, {\em Block-diagonal preconditioning for
  spectral stochastic finite-element systems}, IMA Journal of Numerical
  Analysis, 29 (2009), pp.~350--375.

\bibitem{EW12}
{\sc E.~Rosseel and G.~N. Wells}, {\em Optimal control with stochastic {PDE}
  constraints and uncertain controls}, Computer Methods in Applied Mechanics
  and Engineering, 213-216 (2012), pp.~152--167.

\bibitem{schollwock-2005}
{\sc U.~Schollw\"ock}, {\em The density--matrix renormalization group}, Reviews
  of Modern Physics, 77 (2005), pp.~259--315.

\bibitem{ifiss3.4}
{\sc D.~Silvester, H.~Elman, and A.~Ramage}, {\em {I}ncompressible {F}low and
  {I}terative {S}olver {S}oftware ({IFISS}) version 3.4}, August 2015.
\newblock {\tt http://www.manchester.ac.uk/ifiss/}.

\bibitem{JS14}
{\sc J.~Sogn}, {\em Stabilized finite element methods for the {B}rinkman
  equation on fitted and fictitious domains}, Master's {T}hesis, University of
  Oslo, 2014.

\bibitem{StoB13}
{\sc M.~Stoll and T.~Breiten}, {\em {A low-rank in time approach to
  {PDE}-constrained optimization}}, SIAM Journal on Scientific Computing, 37
  (2015), pp.~B1 -- B29.

\bibitem{SW13}
{\sc M.~Stoll and A.~Wathen}, {\em {All-at-once solution of time-dependent
  Stokes control}}, Journal of Computational Physics, 232 (2013), pp.~498--515.

\bibitem{TKXP12}
{\sc H.~Tiesler, R.~M. Kirby, D.~Xiu, and T.~Preusser}, {\em Stochastic
  collocation for optimal control problems with stochastic {PDE} constraints},
  SIAM Journal on Control and Optimization, 50 (2012), pp.~2659 -- 2682.

\bibitem{book::FT2010}
{\sc F.~Tr{\"o}ltzsch}, {\em {Optimal Control of Partial Differential
  Equations: Theory, Methods and Applications}}, American Mathematical Society,
  2010.

\bibitem{white-dmrg-1993}
{\sc S.~R. White}, {\em Density matrix algorithms for quantum renormalization
  groups}, Physical Review B, 48 (1993), pp.~10345--10356.

\bibitem{NW38}
{\sc I.~Wiener}, {\em The homogeneous chaos}, American Journal of Mathematics,
  60 (1938), pp.~897 -- 936.

\bibitem{XK03}
{\sc D.~Xiu and G.~E. Karniadakis}, {\em A new stochastic approach to transient
  heat conduction modeling with uncertainty}, International Journal of Heat \&
  Mass Transfer, 46 (2003), pp.~4681--4693.

\end{thebibliography}
