% lides for presentation at Casa colloquium
%
%\input {slide-macros-v4.tex}
\documentclass [
%lignesbleues %mooi, maar zwarte achtergrond, dus figuur onzichtbaar
%frames %niet zo mooi
azure %heel mooi
%contemporain
%nuancegris %niet
%troispoints % kermis
%darkblue %should take this with slideColor and colorBG
%alienglow %green on black
%autumn %mooi, maar dan TU logo wit
,pdf,slideColor,colorBG
%,nototal
] {prosper}
\usepackage {algorithmic,amsmath,amssymb}
\newcommand \bm [1] {{\bf #1}}
\newcommand \eps {\varepsilon}
\renewcommand \vec [1] {{\bf #1 }}
\renewcommand \emph [1] {\textcolor {lightcyan} {#1}}
\newrgbcolor{lightcyan}{0.0 0.9 0.9}
%\renewcommand{\bottomtext}{MG presentation, Jos van Kan,
%  DIAM, Delft, November 2003 (slide \thepage)}
\title {Under pressure!}
\subtitle {{\small Fast pressure calculations in the Deft package}}
\author {Jos van Kan}
\email {j.vankan@math.tudelft.nl}
\institution {Delft Institute of Applied mathematics (DIAM)}
\slideCaption {Jos van Kan (DIAM), CASA talk February 15, 2006}
%\Logo (0,-1.1){\includegraphics[width=1.5cm] {tud_01logobl-zw.eps}}
\Logo (0,-1.1){\includegraphics[width=1.5cm] {tulogo.eps}}

\begin{document}
\DefaultTransition{Box}

% Title Slide
\maketitle
% Slide Deft package
\begin{slide}{Deft package}
\emph{De}lft \emph{F}low and \emph{T}ransport
\begin{itemize}
\item Navier Stokes equations for \emph{incompressible} flow on general domains
\item Offshoot of ISNaS (\emph{I}nformation \emph{S}ystem \emph{Na}vier \emph{S}tokes)
%\item \emph {Use only proven technology}
\end{itemize}
\end{slide}
%Slide Deft package design decisions
\begin{slide}{Design decisions}
\begin{itemize}
\item Finite Volume Method
\item Rectangular blocks of curvilinear coordinates
\item Staggered grid
\item \emph {Time-dependent} algorithm
\item Pressure correction for the incompressibility condition
\end{itemize}
\end{slide}
%Slide Navier Stokes equations
\begin{slide}{Navier Stokes equations}
\begin{align*}
\frac{\partial \vec u}{\partial t} +  \nabla \cdot (\vec u\vec u^T) + \nabla p&=
\nabla \cdot T + \vec f\\
\nabla \cdot \vec u &=0
\end{align*}
\begin{itemize}
\item $\vec u$ velocities
\item $p$ pressure
\item $T$ stress tensor
\item $\vec f$ body forces like gravity
\end{itemize}
Variations in density are not taken into account. $p$ and $T$ are
scaled quantities.
\end{slide}
%Slide Boundary conditions
%\begin {slide} {Boundary Conditions}
%\begin{itemize}
%\item Essential: Velocity given
%\item Natural: Viscous stress given
%\end{itemize}
%\begin{itemize}
%\item [2D] Two boundary conditions: normal and tangential. 
%\item [3D] Three boundary conditions: normal and \emph {two} tangential.
%\end{itemize}

%\end {slide}
%Slide Space discretization
\begin{slide}{Method of lines}
\begin {align*}
\frac{d\vec u_h}{dt} + G\vec p_h &= N(\vec u_h) + L(\vec u_h) + \vec f_h\\
D\vec u_h &= 0
\end {align*}
\end {slide}
%Slide Time integration
\begin{slide}{Time Integration}
Implicit time integration, like for instance Crank-Nicolson
\begin {align*}
\vec u_h^{n+1}+&\Delta t G\vec p_h^{n+1/2}=\\
\vec u^n_h+&\frac12\Delta t\left(N(\vec u_h^{n+1})+
L(\vec u_h^{n+1})+\vec f_h^{n+1}\right) + \\
&\frac12\Delta t\left(N(\vec u_h^{n})+ L(\vec u_h^{n})+\vec f_h^n\right)\\
D\vec u_h^{n+1}&=0
\end {align*}
\end{slide}
%Slide Matrix structure
\begin {slide}{Matrix Structure}
\begin{center}
\includegraphics [height=6cm]{matrix_struc.eps}
\end{center}
\end{slide}
%Slide Pressure correction
\overlays {2}{
\begin {slide}
{Pressure Correction}
\begin {align*}
\vec u_h^{*}+\Delta t G\vec p_h^{n-1/2}&=\\
\vec u^n_h+&\frac12\Delta t\left(N(\vec u_h^{*})+L(\vec u_h^{*})+\vec f_h^{n+1}\right) + \\
&\frac12\Delta t\left(N(\vec u_h^{n})+ L(\vec u_h^{n})+\vec f_h^n\right)\\
\vec u_h^{n+1} - \vec u_h^* + &\Delta t G \Delta \vec p=0\\
-D\vec u_h^* + &\Delta t DG \Delta \vec p =0\\
\vec p^{n+1/2}&=\vec p^{n-1/2} + \Delta \vec p\\
%\onlySlide*{1}{-D\vec u_h^* + &\Delta t DG \Delta \vec p =0}%
%\onlySlide*{2}{\color {red} -D\vec u_h^*+ &\color {red}  \Delta t DG \Delta \vec p =0}%
\end{align*}
\end {slide}}

%Slide pressure equation
%\begin{slide} {Pressure equation}
%Solve
%\begin{equation*}
%DG \Delta \vec p = D\vec u_h^*/\Delta t
%\end{equation*}
%\begin {itemize}
%\item $D$ discretization of the divergence
%\item $G$ discretization of the gradient
%\item \emph {Right hand side is $O(\Delta t).$ }
%\end {itemize}
%\end{slide}

%Slide Structure pressure equation
%\begin {slide} {Molecule structure}
%\begin {itemize}
%\item Cartesian: 5 (7) point molecule in 2D (3D)
%\item General: 9 (19) point molecule in 2D (3D)
%\end{itemize}
%\end {slide}
%
\begin {slide} {$p$ Matrix structure 2D}
Example: $n \times n$ block, $N=n^2.$\\
\begin{itemize}
\item Tridiagonal block matrix of tridiagonal matrices
\item Bandwidth: $O(\sqrt N)$
\item Flops $LU$ decomposition: $O(N^2).$ (One time only on fixed domains)
\item Flops $LU$ backsubstitution: $O(N^{3/2})$ 
\item Flops matrix vector multiplication: $O(N).$
\end{itemize}
\end {slide}

\begin {slide} {$p$ Matrix structure (3D)}
Example: $n \times n \times n$ block, $N=n^3.$\\
\begin{itemize}
\item Tridiagonal block matrix of tridiagonal blockmatrices of tridiagonal matrices
\item Bandwidth $O(N^{2/3})$
\item Flops $LU$ decomposition: $O(N^{7/3}).$ (Only once on fixed domains)
\item Flops $LU$ backsubstitution $O(N^{5/3})$
\item Flops matrix vector multiplication: $O(N)$
\end{itemize}
\end {slide}

\begin {slide} {A Classic: Defect Correction}
Solve $A \vec x = \vec b$\\
\begin{algorithmic}
\STATE Presets: $\vec x^0 =0, \vec r^0 = \vec b - A\vec x^0=\vec b, k=0$
\WHILE {$\|r^k\|_\infty > \eps \|b\|_\infty$} 
\STATE Solve $P\vec c^k = \vec r^k$ \COMMENT {P is preconditioner}
\STATE $\vec x^{k+1}= \vec x^k + \vec c^k$
\STATE $\vec r^{k+1}= \vec r^k - A\vec c^k$
\STATE $k=k+1$
\ENDWHILE
\end{algorithmic}
\bigskip
Each iteration requires $O(N)$ flops.
\end{slide}

\begin {slide} {Preconditioners}
\begin{large}\emph{Classic}\end{large}
With $A=D-L-U$
\begin{itemize}
\item Jacobi: $P=D$
\item Gauss-Seidel: $P=D-L$
\item Successive overrelaxation: $P=(D/\omega-L)$
\end{itemize}
\begin{large}\emph{Modern}\end{large}
%\end {slide}
%\begin {slide} {Preconditioners II}
With $A=LU$
\begin {itemize}
\item Incomplete LU (ILU): \\
$A= \tilde L \tilde U + E$, $P = \tilde L \tilde U$. $\tilde L$ and $\tilde U$ sparse, usually the same sparsity pattern as $A$.
\item Incomplete Block LU (IBLU).
\end {itemize}
\end {slide}

\begin {slide} {DC Error Reduction}
\begin{align*}
\vec r^{k+1}=\vec r^k-A \vec c^k&=(I-AP^{-1})\vec r^k\\
\vec A^{-1}\vec b - \vec x^{k+1} &= A^{-1}\vec r^{k+1}=\eps^{k+1}\\
\eps^{k+1}=A^{-1}(I-AP^{-1})A \eps^k&=(I-P^{-1}A)\eps^k
\end{align*}
Reduction governed by spectral radius of $(I-AP^{-1})$.\\
For the Laplacian:
\begin{itemize}
\item Jacobi and Gauss-Seidel: $1-O(h^2)$
\item SOR with optimal $\omega$ and a whole slew of other conditions: $1-O(h)$.
\end{itemize}.
\end {slide}
\overlays {2}{
\begin {slide} {Effectiveness of DC}
\onlySlide*{1}{How many iterations to gain a decimal digit?
\begin{align*}
\eps^{n+k}&=\lambda_1^k\eps^n\\
\lambda_1^k&=0.1\\
k \log \lambda_1 &= - \log 10\\
k&=-\frac {2.3}{\log \lambda_1}\\
k&=O(\frac{1}{1-\lambda_1})
\end{align*}}
\fromSlide{2}{
Jacobi and Gauss-Seidel $O(h^{-2})$ iterations.\\
SOR $O(h^{-1})$ iterations.\\
In 2D:\\
\begin {itemize}
\item Jacobi and GS $O(N^2)$ flops, worse than $LU$ 
\item SOR $(O(N^{3/2})$ flops, order equal to $LU$ 
\end{itemize}
In 3D:\\
\begin{itemize}
\item Jacobi and GS $O(N^{5/3})$ flops, order equal to $LU$ 
\item SOR $O(N^{4/3})$ flops, \color{lightcyan} better than $LU$ 
\end{itemize}}

\end {slide}}

\overlays {4} {
\begin{slide} {Convergence properties}
\begin{center}
\onlySlide*{1}{Damped Jacobi, 10 iterations\\\includegraphics[height=7cm]{conv_prop01.eps}}%
\onlySlide*{2}{Chebyshev10, $\lambda_0=0.1$\\\includegraphics[height=7cm]{conv_prop02.eps}}%
\onlySlide*{3}{Chebyshev10, $\lambda_0=0.03$\\\includegraphics[height=7cm]{conv_prop03.eps}}%
\onlySlide*{4}{Chebyshev10, $\lambda_0=0.01$\\\includegraphics[height=7cm]{conv_prop04.eps}}%
\end{center}
\end{slide}}


\begin {slide}
{Gradient Methods}
\begin{center}
\includegraphics[height=4cm]{clsspc_cyan.eps}
\end{center}
\begin {itemize}
\item "Best" polynomial on the fine structure of the spectrum
\item \emph{Gradient Methods are always better than Defect Correction}
\item {\begin{footnotesize}Irregular convergence behaviour\end{footnotesize}}
\end {itemize}
\end {slide}

\begin {slide}
{Multigrid}
\emph{Defect Correction: very effective on a \textcolor{white}{part} of the spectrum.}\\
\includegraphics[height=5cm]{conv_prop01.eps}\\
The eigenspace of the spectral interval $(0.2,1)$ is virtually reduced to 0 in a few iterations.
\end {slide}
\begin {slide} {1D example}
Consider $-\frac{d^2u}{dx^2}=f, \quad u(0)=u(1)=0$\\
Discretize  into $N$ intervals):
$A\vec u = \vec f$\\
Eigenvalues of $1-P^{-1}A$ are $$\lambda_k = 1-\sin^2 \frac {k\pi}{2 N}, k=1,\dots,N-1.$$ 
Corresponding eigenvectors $$v_{kj}=\sin \frac{kj\pi}{N},k,j = 1,\dots,N-1.$$ 
Eigenvalues \emph{close to 1} correspond to {smooth} eigenvectors, \emph{also for the Laplacian in 2 and 3D.}\\
\end {slide}
\begin {slide}
{Rough and smooth spectrum}
\begin {itemize}
\item \emph {Rough} part of the spectrum: defect correction,  \emph{smoother} in MG speak.
\item \emph {Smooth} part of the spectrum: solve problem on coarser grid and interpolate. \emph {Coarse grid correction} in MG speak.
\end {itemize}
\end {slide}
\begin{slide}
{Restriction and prolongation}
Fine grid correction: $A_h \vec c_h = \vec r_h$\\
Coarse grid correction $A_H \vec c_H = \vec r_H$\\
Transfer operators:
\begin{itemize}
\item $P_{hH}$: prolongation from coarse to fine grid. Interpolation usually.
\item $R_{Hh}$: restriction from fine grid to coarse grid. $R=P^T$ in symmetric problems.
\end{itemize}
The coarse grid correction becomes:
$R_{Hh}A_hP_{hH} \vec c_H = R_{Hh} \vec r_h$\\
\end{slide}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\overlays{3}{
\begin {slide}{Two Grid Algorithm}
\begin{algorithmic}
\fromSlide{1}
{\STATE Presets: $\bm u_h^0$, $\bm r_h^0 =\bm f_h - A \bm u_h^0$
\STATE $\bm u_h^{\text{prs}} = S(\bm u_h^0 , \bm b, A, n_0)$\COMMENT 
{Presmoothing}\STATE $\bm r_H = R_{Hh}\bm r_h$}
\fromSlide{2}{
\STATE Solve $A_H \bm  c_H = \bm r_H$
\STATE $\bm u_h^\text{cgc} = \bm u_h^{\text{prs}} + P_{hH} \bm c_H$ \COMMENT{Coarse Grid 
Correction}}
\fromSlide{3}{ 
 \STATE  $\bm u_h^{\text{pos}} = S(\bm u_h^ \text{cgc}, \bm b, A, 
n_1)$\COMMENT {Postsmoothing}}
\end{algorithmic}
\end {slide}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\overlays{4}{
\begin {slide}{Multi Grid Algorithm}
\begin {algorithmic}
\REQUIRE $A_{\ell+1}=R_{\ell+1,\ell} A_{\ell} P_{\ell,\ell+1}$ have been calculated on all levels
\fromSlide{1}{
\STATE {\textbf{MGRecursive 
$(A_\ell, \bm r_\ell, \bm c_\ell, \ell)$}}
\IF {$\ell<p$}}
\fromSlide{2}{\STATE $\bm c_\ell = S(\bm 0 , \bm r_\ell, A_\ell, n_0)$\COMMENT {Presmoothing}
\STATE $\bm r_{\ell+1}=R_{\ell+1,\ell} (\bm r_\ell - A_\ell \bm c_\ell) $ 
\COMMENT {Calculate coarse grid residual}}
\fromSlide{3}{
\STATE call MGRecursive $(A_{\ell+1}, \bm r_{\ell+1}, \bm c_{\ell+1}, 
\ell+1) $
\STATE $\bm c_\ell = \bm c_\ell + P_{\ell,\ell+1} \bm c_{\ell+1}$ \COMMENT 
{Coarse grid correction}}
\fromSlide{4}{
\STATE $\bm c_{\ell}=S(\bm c_\ell, \bm r_\ell, A_\ell, n_1)$ 
\COMMENT{Postsmoothing} }
\ELSE 
\STATE Solve $A_{p} \bm c_{p}= \bm r_{p}$\COMMENT {Direct solution on 
coarsest level}
\ENDIF
\end {algorithmic}
\end{slide}
}
\begin {slide} {The MG miracle}
\begin {itemize}
\item Spectrum of $I-P^{-1}A$ bounded away from $1$ \emph{uniformly} in $h$. 
\item Number of iterations does not depend on $h$.
\item The workload is theoretically $O(N)$ flops.
\end {itemize}
\emph{But how big is the multiplicative constant going to be?}
\end {slide}

\begin {slide} {Robust Blackbox}
Wishlist:
\begin {itemize}
\item Good smoother under various circumstances (anisotropy, stretched and skew cells)
\item Arbitrary number of points in either direction
\end {itemize}
\end {slide}
\overlays{2}{
\begin{slide}{Smoothers tested}
\begin{itemize}
\item (Alternating) damped line Jacobi, 1 or 2 postsmoothing steps, no presmoothing
\onlySlide*{1}{\includegraphics[width=3cm]{jac_x.eps}\quad\quad
\includegraphics[width=3cm]{jac_y.eps}%
}
\onlySlide*{2}{\item Incomplete Block $LU$ decomposition, 1 postsmoothing step, no presmoothing
\includegraphics[width=3cm]{blockL.eps}\quad\quad
\includegraphics[width=3cm]{blockU.eps}%
}
\end{itemize}
\end{slide}}

\begin {slide} {Comparison}
\begin{center}
\includegraphics[width=9cm]{compar.eps}
\end{center}
\end {slide}

\begin {slide} {Recursion}
\begin{itemize}
\item Line Jacobi is recursive \emph {per line} and can be massively parallelized, especially in 3D.
\item What about IBLU? Classic IBLU is fully recursive.
\end{itemize}
\end{slide}

\begin {slide} {Divide and Conquer}
The inversion of an $n\times n$ tridiagonal matrix can be executed in ${}^2\log n$ non recursive steps.

\begin {multline*}
(I+LD^{-1}+UD^{-1})(D-L-U)=\\
D-LD^{-1}U-UD^{-1}L-LD^{-1}L-UD^{-1}U
\end {multline*}

$$\quad\quad=D_1-L_1-U_1$$

Bandwith is doubled in this operation. 
\end {slide}

%\overlays {5}{
\begin{slide}
{Incomplete Block Div and Conq}
\begin {itemize}
\item Use the same formula, \emph {interpreted als blocks}
\item Use incomplete versions of $LD^{-1}U$ etc.
\item $D$, $L$ and $U$ are (block)diagonals, consisting of tridiagonal blocks. 
\item We need 7 diagonals of $D^{-1}$
\item Calculate from the productform of the inverse, keeping only 7 diagonals in ${}^2\log n$ steps
\end {itemize}
Recursion uses $O({}^2\log n)$ steps as claimed.
\end{slide}%}

\end {document}
  
