Homec4science
Diffusion SPClibs (master) Browse

spclibs/multigrid/docs/multigrid_2d.tex
93a1e9712d04master

Last Change
Raw File

multigrid_2d.tex
View Options

%
% @file multigrid_2d.tex
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% SPClibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
\documentclass[a4paper]{article}
\usepackage{linuxdoc-sgml}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{amsmath}
\title{Multigrid for Finite Elements using Splines.}
\author{Trach-Minh Tran, Stephan Brunner}
\date{v0.1, January 2013}
\abstract{A multigrid formulation for finite elements is
derived, using variational principles. More specifically the grid
transfer operators will be derived and tested in 2D Cartesian and
cylindrical geometry for arbitrary order B-Splines.}
\begin{document}
\maketitle
\tableofcontents
\section{The Model Problems}
\subsection{Cartesian Geometry}
The following second-order boundary value problem is considered
\begin{equation}
\label{eq:cartesian_problem}
\begin{split}
-\left[\frac{\partial^2}{\partial x^2}
+\frac{\partial^2}{\partial y^2} \right] u(x,y) &= f(x,y) \qquad 0\le x\le 1,\quad 0\le y\le 1 \\
u(0,y) = u(1,y) &= u(x,0) = u(x,1) = 0.
\end{split}
\end{equation}
By choosing
\begin{equation*}
f(x,y) = \sin (\pi k_xx + \pi k_yy),
\end{equation*}
where $k_x$ and $k_y$ are integers, the solution of the BVP is simply
\begin{equation*}
u(x,y) = \frac{\sin(\pi k_xx+\pi k_yy)}{\pi^2(k_x^2+k_y^2)}.
\end{equation*}
Using a weak formulation on Eq.(\ref{eq:cartesian_problem}) and a grid
of $N_x\times N_y$ intervals, one obtains the following discretized
linear system
\begin{equation}
\sum_{i'=1}^{N_x+p}\sum_{j'=1}^{N_y+p}A_{iji'j'}u_{i'j'} = b_{ij},
\qquad i=1,\ldots,N_x+p,\quad j=1,\ldots,N_y+p,
\end{equation}
where the unknowns $u_{ij}$ are the Spline (of order $p$) expansion
coefficients of the solution
\begin{equation}
u(x,y) = \sum_{i=1}^{N_x+p}\sum_{j=1}^{N_y+p}u_{ij}\Lambda_i(x)\Lambda_j(y),
\end{equation}
and the matrix $A$ and right hand side $b$ are determined from
\begin{align}
A_{iji'j'} &= \int^1_0\int^1_0 dxdy
\left[\Lambda'_{i'}(x)\Lambda_{j'}(y)\Lambda'_i(x)\Lambda_j(y) +
\Lambda_{i'}(x)\Lambda'_{j'}(y)\Lambda_i(x)\Lambda'_j(y) \right], \\
b_{ij} &= \int^1_0\int^1_0 dxdy\Lambda_i(x)\Lambda_j(y)f(x,y).
\end{align}
Note that using a Gauss quadrature with $\lceil(2p+1)/2\rceil$ points
per interval to calculate the matrix $A$ would yield an exact integration.
\subsection{Cylindrical Geometry}
The following second-order boundary value problem is considered:
\begin{equation}
\label{eq:cylindrical_problem}
\begin{split}
-\left[\frac{1}{r}\frac{\partial}{\partial r}r\frac{\partial}{\partial r}
+ \frac{1}{r^2} \frac{\partial^2}{\partial\theta^2}
\right]u(r,\theta) &= f(r,\theta) \qquad 0\le r\le 1,\quad
0 \le \theta < 2\pi \\
u(1,\theta) &= 0, \\
\end{split}
\end{equation}
By choosing
\begin{equation*}
f(r,\theta) = j^2_{ms} J_{m}(j_{ms}r)\cos(m\theta),
\end{equation*}
where $m$ is an integer and $j_{ms}$, the $s^{th}$ zero of $J_{m}$,
the solution of this BVP is
\begin{equation*}
u(r,\theta) = J_{m}(j_{ms}r)\cos(m\theta).
\end{equation*}
Using a weak formulation on Eq.(\ref{eq:cylindrical_problem})and a grid
of $N_r\times N_\theta$ intervals, one obtains the following discretized
linear system
\begin{equation}
\sum_{i'=1}^{N_r+p}\sum_{j'=1}^{N_\theta}A_{iji'j'}u_{i'j'} = b_{ij},
\qquad i=1,\ldots,N_r+p,\quad j=1,\ldots,N_\theta,
\end{equation}
where the unknowns $u_{ij}$ are the Spline (of order $p$) expansion
coefficients of the solution
\begin{equation}
u(r,\theta) = \sum_{i=1}^{N_r+p}\sum_{j=1}^{N_\theta}u_{ij}\Lambda_i(r)\Lambda_j(\theta),
\end{equation}
and the matrix $A$ and right hand side $b$ are determined from
\begin{align}
A_{iji'j'} &= \int^1_0\int^{2\pi}_0 rdrd\theta
\left[\Lambda'_{i'}(r)\Lambda_{j'}(\theta)\Lambda'_i(r)\Lambda_j(\theta) +
\frac{1}{r^2}
\Lambda_{i'}(r)\Lambda'_{j'}(\theta)\Lambda_i(r)\Lambda'_j(\theta)
\right], \\
b_{ij} &= \int^1_0\int^{2\pi}_0 rdrd\theta\Lambda_i(r)\Lambda_j(\theta)f(r,\theta).
\end{align}
Note that $A$ has an $1/r$ singularity in the
integrand. For $m\neq0$, this should not be problematic since the
converged solution behaves as $\sim r^m$ near $r=0$. The case $m=0$
will be investigated numerically latter in this report, together withe
the $m\neq 0$ case.
\section{Restriction Operator}
In the following, let us use the superscripts $h$ and $2h$ to denote
quantities defined respectively on a \emph{fine} ($N_x\times N_y$
or $N_r\times N_\theta$) and a \emph{coarser} ($N_x/2\times N_y/2$
or $N_r/2\times N_\theta/2$) grid.
The two grid transfers required in the standard \emph{multigrid}
\cite{MG1D,Briggs} are:
\begin{enumerate}
\item the \emph{restriction} of the right hand side: $\mathbf{b}^{h}
\longrightarrow \mathbf{b}^{2h}$ and
\item the \emph{prolongation} of the solution: $\mathbf{u}^{2h}
\longrightarrow \mathbf{u}^{h}$.
\end{enumerate}
Noting that the basis functions $\Lambda^{2h}_i(x)$, which are \emph{piecewise}
$C^{p-1}$ polynomials with \emph{breaks} on the \emph{coarse} grid
points $x^{2h}_k=(2h)k$ can be also considered as \emph{piecewise}
$C^{p-1}$ polynomials with \emph{breaks} on the \emph{fine} grid
$x^h_k=kh$, they can be expressed \emph{uniquely} as a linear
combination of the \emph{fine} grid basis functions:
\begin{equation}
\label{eq:basis_transf}
\Lambda^{2h}_i(x) = \sum^{N+p}_{i'=1}c_{ii'}\Lambda^h_{i'}(x), \quad
i=1,\ldots,N/2+p.
\end{equation}
The (rectangular) matrix $c_{ii'}$ can be identified as the
one-dimensional \emph{restriction} $\mathbf{R}$ since
\begin{equation*}
b^{2h}_i = \int_0^1 dx f(x)\Lambda^{2h}_i(x) =
\sum^{N+p}_{i'=1}c_{ii'}\;b^h_{i'} =
\sum^{N+p}_{i'=1}R_{ii'}\;b^h_{i'}.
\end{equation*}
It can be computed by simply projecting Eq.(\ref{eq:basis_transf})
on the fine grid basis function $\Lambda^h_{j}(x)$ \cite{MG1D}:
\begin{equation}
\sum^{N+p}_{i'=1}R_{ii'}\underbrace{\int_0^1 dx
\Lambda^h_{i'}(x)\Lambda^h_{j}(x)}_{M^h_{i'j}} = \underbrace{\int_0^1 dx
\Lambda^{2h}_i(x)\Lambda^h_{j}(x)}_{M^{2h,h}_{i'j}} \Longrightarrow
\mathbf{R}=\mathbf{M}^{2h,h}\cdot(\mathbf{M}^{h})^{-1}.
\end{equation}
It should be stressed that the representation for $\Lambda^{2h}_i(x)$
in Eq.(\ref{eq:basis_transf}) is \emph {unique}. This is
checked by verifying that the same matrix $R_{ii'}$ is obtained using
for example the \emph{collocation} methods. One such method, which is used for
this check is detailed in Appendix \ref{sec:colloc}. The calculated
grid transfer matrices for linear, quadratic and cubic periodic and
non-periodic Splines are given in \cite{MG1D}.
Denoting the restriction on $x$ and $y$ respectively
by $\mathbf{R}^x$ and $\mathbf{R}^y$, the two-dimensional restriction of $b^h_{ij}$ is
defined as
\begin{equation*}
b^{2h}_{ij} = \int_0^1\int_0^1 dxdy f(x,y)\Lambda^{2h}_i(x)\Lambda^{2h}_j(y) =
\sum^{N+p}_{i'=1}\sum^{N+p}_{j'=1}R^x_{ii'}R^y_{jj'}b^h_{i'j'},
\end{equation*}
and thus
\begin{equation}
\label{eq:restriction}
\boxed{\mathbf{b}^{2h} = \mathbf{R}^x \cdot \mathbf{b}^{h}
\cdot (\mathbf{R}^y)^T.}
\end{equation}
\section{Prolongation Operator}
Using Eq.(\ref{eq:basis_transf}) (with $c_{ii'}=R_{ii'}$), the
solution at the coarse grid can be expressed as
\begin{equation*}
u^{2h}(x) = \sum_{i=1}^{N/2+p}u^{2h}_{i}\Lambda^{2h}_{i}(x) =
\sum_{i'=1}^{N+p}\left[\sum_{i=1}^{N/2+p} R_{ii'}u^{2h}_{i}\right]
\Lambda^h_{i'}(x) =
\sum_{i'=1}^{N+p}\underbrace{\left[\sum_{i=1}^{N/2+p} (R)^T_{i'i}u^{2h}_{i}\right]}_{\tilde{u}^h_{i'}}
\Lambda^h_{i'}(x),
\end{equation*}
from which one obvious choice for the \emph{prolongation} operator would be
\begin{equation}
\mathbf{P} = \mathbf{R}^T = (\mathbf{M}^{h})^{-1}\cdot\mathbf{M}^{h,2h}.
\end{equation}
Generalization to a two-dimensional prolongation is obtained as
follows, where summation over repeated indices is assumed:
\begin{equation*}
u^{2h}(x,y) = u^{2h}_{ij}\Lambda^{2h}_{i}(x)\Lambda^{2h}_{j}(y) =
\left[R^x_{ii'}u^{2h}_{ij}R^y_{jj'}\right]\Lambda^{h}_{i'}(x)\Lambda^{h}_{j'}(y)
\end{equation*}
which leads to the prolonged solution $\tilde{\mathbf{u}}^{h}$ given by
\begin{equation}
\label{eq:prolongation}
\boxed{\tilde{\mathbf{u}}^{h} = \mathbf{P}^x \cdot \mathbf{u}^{2h} \cdot
(\mathbf{P}^y)^T .}
\end{equation}
It should be noted here that, while the restricted right hand side
$\mathbf{b}^{2h}$ as defined in Eq.(\ref{eq:restriction}) is
\emph{exactly identical} to the assembled right hand side, the
prolonged solution $\tilde{\mathbf{u}}^{h}$ defined in
Eq.(\ref{eq:prolongation}) is just a
representation of $u^{2h}(x,y)$ on the fine mesh and \emph{not} the
solution $u^h(x,y)$ which can only be obtained by solving the problem
on the fine mesh!
\section{Numerical Experiments}
The multigrid performance can be characterized by looking at the
convergence of the residual Euclidean norm for the linear system
$\mathbf{A}\mathbf{u}=\mathbf{b}$:
\begin{equation}
\label{eq:resid}
\|\mathbf{r}\|_2 = \|\mathbf{b}-\mathbf{A}\mathbf{u}\|_2.
\end{equation}
When the \emph{exact} solution $u(x,y)$ is known, the \emph{discretization error}
can defined as
\begin{equation}
\label{eq:discerr}
\|e\|_2 = \sqrt{\int dV\left[\sum_{ij}u_{ij}\Lambda_{ij}(x,y)-u(x,y)\right]^2}
\end{equation}
and computed using a Gauss quadrature. Note that for Splines of order
$p$, $\|e(x,y)\|_2(h)$ converges to zero as $O(h^{p+1})$.
\subsection{Cartesian Geometry}
The multigrid performances for varying problem sizes are displayed in
Fig.(\ref{fig:linear_mg2d}) for linear Splines and
Fig.(\ref{fig:cubic_mg2d}) for cubic Splines. They show that the number
of iterations required for convergence (abount 3 for both linear and cubic
Splines) is insensitive to the problem sizes. Compared to direct
methods, the multigrid should scale much better for large problem
sizes, as indicated in Table~\ref{tab:comparison1}. For this model
problem, using cubic Splines seems to converge
slightly faster than linear Splines!
\begin{figure}
\centering
\includegraphics[angle=0,width=\hsize]{linear_mg2d}
\caption{Performance of the multigrid $V(2,1)$ scheme using a
Gauss-Seidel relaxation and \emph{linear Splines} for different
problem sizes. The size of the \emph{coarsest} grid is $2\times 2$.}
\label{fig:linear_mg2d}
\end{figure}
\begin{figure}
\centering
\includegraphics[angle=0,width=\hsize]{cubic_mg2d}
\caption{Performance of the multigrid $V(2,1)$ scheme using a
Gauss-Seidel relaxation and \emph{cubic Splines} for different
problem sizes. The size of the \emph{coarsest} grid is $2\times 2$.}
\label{fig:cubic_mg2d}
\end{figure}
\begin{table}[htb]
\centering
\begin{tabular}{|c|c|c|c|c|c|}\hline
& \multicolumn{2}{c|}{Linear Splines}
& \multicolumn{2}{c|}{Cubic Splines} \\ \hline
$N$ & $V(2,1)$ & Direct & $V(2,1)$ & Direct \\ \hline
16 & 8.844E-04 & 2.051E-03 & 2.653E-03 & 3.970E-03 \\
32 & 1.661E-03 & 5.345E-03 & 4.983E-03 & 1.540E-02 \\
64 & 5.766E-03 & 2.054E-02 & 1.730E-02 & 7.492E-02 \\
128 & 2.347E-02 & 3.288E-01 & 7.042E-02 & 1.060E+00 \\
\hline
\end{tabular}
\caption{Times (in seconds) used by a the \emph{direct sparse} solver MUMPS-4.10.0
for different problem sizes versus the times used by \emph{three}
multigrid $V(2,1)$ cycles. The Intel Fortran-13.0 compiler is used on an Intel
i7 platform.}
\label{tab:comparison1}
\end{table}
The effects of the relaxation parameters $\nu_1,\nu_2$ on the
multigrid performnace (Fig.(\ref{fig:cubic_mg2d_relax})) indicates
that only a few relaxations are sufficient to achieve a good multigrid
performance. Further analysis of the computational cost is required
however to determine the \emph{optimal} $\nu_1,\nu_2$.
Finally, the effects of the number of grid levels are analyzed in
Fig.(\ref{fig:cubic_mg2d_levels}). In addition to the computational
cost (see Table~\ref{tab:comparison2}), the memory required for the
\emph{direct solver} at the coarsest grid level should be taken into
account for the choice of the optimal number of grid levels,
especially for very large problems.
\begin{figure}
\centering
\includegraphics[angle=0,width=\hsize]{cubic_mg2d_relax}
\caption{Effect of the number of the relaxation sweeps $\nu_1,\nu_2$
on the performance of the multigrid $V(\nu_1,\nu_2)$-cycle for
\emph{Cubic Splines}. The finest grid has $128\times 128$
intervals.}
\label{fig:cubic_mg2d_relax}
\end{figure}
\begin{figure}
\centering
\includegraphics[angle=0,width=\hsize]{cubic_mg2d_levels}
\caption{Effect of the number grid levels
on the performance of the multigrid $V(2,1)$-cycle for \emph{Cubic
Splines}. The finest grid has $128\times 128$ intervals.}
\label{fig:cubic_mg2d_levels}
\end{figure}
\begin{table}[htb]
\centering
\begin{tabular}{|c|c|c|c|}\hline
Number of levels & $V(1,0)$ & $V(1,1)$ & $V(2,1)$ \\ \hline
2 & 3.386E-02 & 3.881E-02 & 4.031E-02 \\
3 & 2.923E-02 & 3.398E-02 & 3.605E-02 \\
4 & 2.880E-02 & 3.275E-02 & 3.595E-02 \\
7 & 2.912E-02 & 3.236E-02 & 3.566E-02 \\
\hline
\end{tabular}
\caption{Effects of the times in seconds used per $V$-cyclefor
different number of grid levels and relaxation paramters for a
$128\times 128$ problem. The Intel Fortran-13.0 compiler is used on an Intel
i7 platform.}
\label{tab:comparison2}
\end{table}
\subsection{Cylindrical Geometry}
\newpage
\appendix
\section{Grid transfer matrix by collocation}
\label{sec:colloc}
Let first consider the \emph{periodic case}. Denoting $N$ as the
number of intervals of the fine grid, the \emph{periodic} Spline basis
functions on the \emph{coarse} grid $\Lambda^{2h}_i$ can be expressed
as linear combinations of the \emph{fine} grid Spline basis functions
as:
\begin{equation}
\Lambda^{2h}_i(x) = \sum^{N}_{i'=1}R_{ii'}\Lambda^h_{i'}(x), \quad
i=1,\ldots,N/2.
\end{equation}
For any given $i$, the coefficients $R_{ii'}$ can be calculated by
expressing the relation above on exactly $N$ points on the $x$-grid. For
\emph{odd} Spline order $p$, these \emph{collocation} (or
interpolating) points can be chosen as the \emph{break} points of the fine
grid $x^h_k,\quad k=0,\ldots,N-1$. For \emph{even} values of $p$, the
collocation points should be $x^h_{k+1/2}=(x^h_{k}+x^h_{k+1})/2$ in
order to obtain a non-singular linear system of equations
\cite{BSPLINES}. The resulting system of equations to solve for
$R_{ii'}$ are given below:
\begin{equation}
\begin{split}
p\mbox{ odd}: \qquad &\sum^{N}_{i'=1}\Lambda^h_{i'}(x^h_k)\,R_{ii'} =
\Lambda^{2h}_i(x^h_k), \qquad k=0,\ldots,N-1,\quad i=1,\ldots,N/2, \\
p\mbox{ even}: \qquad &\sum^{N}_{i'=1}\Lambda^h_{i'}(x^h_{k+1/2})\,R_{ii'} =
\Lambda^{2h}_i(x^h_{k+1/2}), \qquad k=0,\ldots,N-1,\quad i=1,\ldots,N/2.\\
\end{split}
\end{equation}
For \emph{non-periodic} Splines, there are $N+p$ and
$N/2+p$ basis functions respectively on the fine and coarse grid:
\begin{equation}
\Lambda^{2h}_i(x) = \sum^{N+p}_{i'=1}R_{ii'}\Lambda^h_{i'}(x), \quad
i=1,\ldots,N/2+p.
\end{equation}
This implies that for any given $\Lambda^{2h}_i$, $N+p$ conditions
are required to determined the $N+p$ terms of row $i$ of the matrix $R_{ii'}$. For odd $p$,
$N+1$ collocation points $x_k,\quad k=0,\ldots,N$ can be used with the
missing $p-1$ equations obtained by expressing all the $(p-1)/2$ derivatives
of $\Lambda^{2h}_i(x)$ at the end points $x_0$ and $x_N$:
\begin{equation}
\frac{d^\alpha}{dx^\alpha}\Lambda^{2h}_i(x) =
\sum^{N+p}_{i'=1}R_{ii'}\frac{d^\alpha}{dx^\alpha} \Lambda^h_{i'}(x)
, \quad \alpha=1,\dots,\frac{p-1}{2} \quad (\mbox{$p$ odd}).
\end{equation}
For \emph{even} $p$, in addition to the $N$ relations obtained with
the collocation points $x_{k+1/2}$ (as in the \emph{periodic} case),
the missing $p$ conditions can be obtained by expressing $\Lambda^{2h}_i$
and its derivatives up to $p/2-1$ at the end points $x_0$ and $x_N$:
\begin{equation}
\frac{d^\alpha}{dx^\alpha}\Lambda^{2h}_i(x) =
\sum^{N+p}_{i'=1}R_{ii'}\frac{d^\alpha}{dx^\alpha} \Lambda^h_{i'}(x)
, \quad \alpha=0,\dots,\frac{p}{2}-1 \quad (\mbox{$p$ even}).
\end{equation}
\begin{thebibliography}{99}
\bibitem{MG1D} {\tt Multigrid Formulation for Finite Elements},\\
\url{https://crppsvn.epfl.ch/repos/bsplines/trunk/multigrid/docs/multigrid.pdf}
\bibitem{Briggs} {W.L.~Briggs, V.E.~Henson and S.F.~McCormick, A
Multigrid Tutorial, Second Edition, Siam (2000)}.
\bibitem{BSPLINES} {\tt BSPLINES} Reference Guide,
\url{https://crppsvn.epfl.ch/repos/bsplines/trunk/docs/bsplines.pdf}
\bibitem{SOLVERS} {\tt The Solvers in BSPLINES},
\url{https://crppsvn.epfl.ch/repos/bsplines/trunk/docs/solvers.pdf}
\end{thebibliography}
\end{document}
c4science · Help