Homec4science
Diffusion SPClibs (master) Browse

spclibs/multigrid/docs/multigrid.tex
03540012fedemaster

Last Change
Raw File

multigrid.tex
View Options

%
% @file multigrid.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}
%\usepackage{verbatim}
%\usepackage[notref]{showkeys}
\title{\tt Multigrid Formulation for Finite Elements.}
\author{Trach-Minh Tran, Stephan Brunner}
\date{v0.2, October 2012}
\abstract{A multigrid formulation for finite elements is
derived, using variational principles. More specifically the grid
transfer operators will be derived and tested in 1D Cartesian,
cylindrical and spherical geometry for arbitrary order B-Splines.}
\begin{document}
\maketitle
\tableofcontents
\section{The discretized problem}
Consider the one-dimensional linear integro-differential problem
\begin{equation}
\label{eq:oned_prob}
\mathcal{L}(u) = f, \qquad 0\le x\le L,
\end{equation}
with suitable boundary conditions. On an \emph{equidistant} mesh with
interval $h=L/N$ and using the \emph{weak form} of Eq.~(\ref{eq:oned_prob}),
the linear system to be solved on this grid (which will be
referred as the \emph{fine} grid) can be written as (see
\cite{SOLVERS}, \cite{BSPLINES}):
\begin{equation}
\label{eq:fine}
\sum_{i'=1}^{N+p}A_{ii'}^hu^h_{i'} = b^h_i, \qquad
A^h_{ii'}=\int_{0}^{L}\Lambda^h_i
\mathcal{L}(\Lambda^h_{i'})\,x^{\alpha}dx, \qquad
b^h_i = \int_{0}^{L}f\Lambda^h_i\,x^{\alpha}dx,
\end{equation}
where $p$ is the order the Splines $\Lambda^h_i$ and
$\alpha=0,1,2$ for Cartesian, cylindrical and spherical coordinates
respectively. It should be noted that the unknowns $u^h_i$ of this
linear system are the \emph{expansion coefficients} of the discretized
solution of the problem $u^h(x)$
\begin{equation}
u^h(x) = \sum_{i'=1}^{N+p}u^h_{i'}\Lambda^h_{i'}(x)
\end{equation}
and the right hand sides $b^h_i$ are defined as the \emph{projection} of
$f(x)$ on the same basis functions, in contrast with the Finite
Differences (FD) or Finite Volume (FV) formulations where $u^h_i$ and
$b^h_i$ are the \emph{nodal values} of $u$ and $f$.
On the \emph{coarser} mesh with interval $2h=2L/N$, the discretized
linear system can be written as
\begin{equation}
\label{eq:coarse}
\sum_{i'=1}^{N/2+p}A_{ii'}^{2h}u^{2h}_{i'} = b^{2h}_i, \qquad
A^{2h}_{ii'}=\int_{0}^{L}\Lambda^{2h}_i
\mathcal{L}(\Lambda^{2h}_{i'})\,x^{\alpha}dx, \qquad
b^{2h}_i = \int_{0}^{L}f\Lambda^{2h}_i\,x^{\alpha}dx.
\end{equation}
\section{Transfer operators}
\label{sec:twogrid}
For simplicity let consider the two-grid procedure \cite{Briggs} which can be
summarized as follow:
\begin{enumerate}
\item Obtain an approximation $\mathbf{u}^h$ on the \emph{fine} grid, using a
Gauss-Seidel (GS) or a weighted Jacobi scheme. This procedure is
also called \emph{smoothing} or \emph{relaxation}.
\item Compute the \emph{residuals}:
$\mathbf{r}^h=\mathbf{b}^h-\mathbf{A}^h\mathbf{u}^h$.
\item Obtain the residuals on the coarse mesh $\mathbf{r^{2h}}$ by
\emph{restriction} of $\mathbf{r}^h$.
\item Direct solve $\mathbf{A}^{2h}\mathbf{e}^{2h}=\mathbf{r^{2h}}$ to obtain
the error $\mathbf{e}^{2h}=\mathbf{u}-\mathbf{u}^{2h}$
\item Interpolate (\emph{prolong}) the error to obtain $\mathbf{e}^h$.
\item Correct the approximation obtained on the fine grid:
$\mathbf{u}^h\leftarrow \mathbf{u}^h+\mathbf{e}^h$.
\item Relax on $\mathbf{A}^{h}\mathbf{u}^{h}=\mathbf{f^{h}}$, using
the previously computed $\mathbf{u}^h$ as a guess.
\end{enumerate}
Steps 3 and 5 are called \emph{grid transfers} and are detailed in the
following. It should be noted that the fine to coarse transfer
(restriction) applies to the right hand side $\mathbf{b}^h$ while the
prolongation applies to the expansion coefficients $\mathbf{u}^{2h}$.
\subsection{Fine to coarse grid transfer (restriction)}
The right hand side on the fine and coarse grid can be written as
\begin{equation*}
\begin{split}
b^{h}_i &= \int_0^Lf\Lambda^{h}_i\,x^{\alpha}dx =
\sum_{i'=1}^{N+p}f^h_{i'}\underbrace{\int_0^L\Lambda^{h}_{i}\Lambda^{h}_{i'}\,x^{\alpha}dx}_
{M^{h,h}_{ii'}}, \\
b^{2h}_i &= \int_0^Lf\Lambda^{2h}_i\,x^{\alpha}dx = \sum_{i'=1}^{N+p}f^h_{i'}
\underbrace{\int_0^L\Lambda^{2h}_{i}\Lambda^{h}_{i'}\,x^{\alpha}dx}_{M^{2h,h}_{ii'}} \\
\end{split}
\end{equation*}
where the expansion $f(x)=\sum_{i=1}^{N+p}f^h_i\Lambda^h_i(x)$ has been used.
Elimination of $\mathbf{f^h}$ leads to the definition of the
\emph{restriction} matrix:
\begin{equation}
\label{eq:restriction}
\mathbf{b}^{2h} = \mathbf{R}^{2h}_{h}\mathbf{b}^h,
\qquad \boxed{\mathbf{R}^{2h}_{h}=\mathbf{M}^{2h,h}(\mathbf{M}^{h,h})^{-1}}.
\end{equation}
Note that the computation of the \emph{mass matrices}
$\mathbf{M}^{h,h}$ and $\mathbf{M}^{2h,h}$ can be done \emph{exactly}
using a Gauss integration with $N_G=\lceil p+(\alpha+1)/2 \rceil$ points.
Another way to derive the restriction operator $\mathbf{R}^{2h}_{h}$
is by noting that the basis functions $\Lambda^{2h}_{i}$ are
\emph{piecewise} $C^{p-1}_h$ \emph{polynomials} with \emph{breaks} on
the fine grid points $x_i=ih$, and thus can be expressed
\emph{uniquely} as
\begin{equation}
\Lambda^{2h}_{i}(x) = \sum_{i'=1}^{N+p}c_{ii'}\Lambda^{h}_{i'}(x),
\quad i=1\ldots N/2+p.
\end{equation}
Projecting this equation on the basis $\Lambda^{h}_j$ then leads to
\begin{equation*}
\begin{split}
\sum_{i'=1}^{N+p}c_{ii'}\int_0^L\Lambda^{h}_{i'}\Lambda^{h}_{j}\,x^{\alpha}dx
&= \int_0^L\Lambda^{2h}_{i}\Lambda^{h}_{j}\,x^{\alpha}dx, \qquad
i=1\ldots N/2+p, \quad j=1,\ldots N+p \\
\Longrightarrow \mathbf{c}\cdot\mathbf{M}^{h,h} &=
\mathbf{M}^{2h,h} \Longrightarrow \mathbf{c} =
\mathbf{M}^{2h,h}(\mathbf{M}^{h,h})^{-1} = \mathbf{R}^{2h}_{h} \\
\end{split}
\end{equation*}
and finally
\begin{equation}
\label{eq:restrict_gen}
\boxed{\Lambda^{2h}_{i}(x) =
\sum_{i'=1}^{N+p}\left(\mathbf{R}^{2h}_{h}\right)_{ii'}\Lambda^{h}_{i'}(x),
\quad i=1\ldots N/2+p}
\end{equation}
Because the expansion coefficients $c_{ii'}$ of $\Lambda^{2h}_{i}(x)$
(rows of the restriction matrix $\mathbf{R}^{2h}_{h}$) on the
fine mesh basis are \emph{unique}, $\mathbf{R}^{2h}_{h}$ should be
independent of the geometry exponent $\alpha$ or more generally, of the
definition of the \emph{projection} (or scalar product) used to
calculate the restriction matrix. Furthermore, since the supports of
both $\Lambda^{h}_i$ and $\Lambda^{2h}_i$ are \emph{compact}, the
matrix $\mathbf{R}^{2h}_{h}$ should be \emph{sparse}.
One can show that, using (\ref{eq:restrict_gen}), the
restriction of the fine mesh FE matrix $\mathbf{A}^h$ is given by
\begin{equation}
\label{eq_coarse_mat}
\mathbf{A}^{2h} = \mathbf{R}^{2h}_h\mathbf{A}^{h}\left(\mathbf{R}^{2h}_h\right)^{T}.
\end{equation}
\subsection{Coarse to fine grid transfer (prolongation)}
Let denote the discretized solution on the coarse mesh of
$\mathbf{A}^{2h}\mathbf{u}^{2h}=\mathbf{R}^{2h}_{h}\mathbf{b}^h$
by
\begin{equation*}
u^{2h}(x) = \sum_{i=1}^{N/2+p}u^{2h}_{i}\Lambda^{2h}_{i}(x),
\end{equation*}
and seek for an approximated solution on the fine mesh $\mathbf{u}^{h}$
\begin{equation*}
u^{h}(x) = \sum_{i=1}^{N+p}u^{h}_{i}\Lambda^{h}_{i}(x).
\end{equation*}
by \emph{prolongation} of $\mathbf{u}^{2h}$ (instead of solving
$\mathbf{A}^{h}\mathbf{u}^{h}=\mathbf{b}^h$). A reasonable solution is
to \emph{minimize} the square of the error norm defined as
\begin{equation*}
\begin{split}
\epsilon^2 &= \|u^{h}(x)-u^{2h}(x)\|^2 \equiv \int_0^L
[u^{h}(x)-u^{2h}(x)]^2\,x^\alpha dx, \\
\frac{\partial\epsilon^2}{\partial u^h_i} &=0 \Longrightarrow
\sum_{i'=1}^{N+p}u^{h}_{i}\int_0^L
\Lambda^{h}_{i}\Lambda^{h}_{i'}\,x^\alpha dx =
\sum_{i'=1}^{N/2+p}u^{2h}_{i}\int_0^L
\Lambda^{h}_{i}\Lambda^{2h}_{i'}\,x^\alpha dx. \\
\end{split}
\end{equation*}
This yields the prolonged (or interpolated) \emph{coarse grid}
solution on the \emph{fine grid}
\begin{equation}
\label{eq:prolongation}
\mathbf{u}^h = \mathbf{P}^h_{2h}\mathbf{u}^{2h}, \qquad
\boxed{\mathbf{P}^h_{2h} =
(\mathbf{M}^{h,h})^{-1}\mathbf{M}^{h,2h}=(\mathbf{R}^{2h}_{h})^T}
\end{equation}
and the coarse FE matrix can be finally expressed as\begin{equation}
\label{eq:coarse_mat}
\boxed{\mathbf{A}^{2h} = \mathbf{R}^{2h}_{h}\mathbf{A}^{h}\mathbf{P}^h_{2h}}
\end{equation}
\subsection{An alternative derivation of grid transfer operators}
Starting from the inter grid transformation of the basis functions
Eq.(\ref{eq:restrict_gen}), the
restriction of $\mathbf{b}^h$ and the prolongation of
$\mathbf{u}^{2h}$ can be derived as follow
\begin{gather*}
b^{2h}_i = \int_0^L f\Lambda^{2h}_i\,x^{\alpha}dx = \sum_{i'=1}^{N+p}
\left(\mathbf{R}^{2h}_{h}\right)_{ii'}\int_0^Lf\Lambda^{h}_{i'}\,x^{\alpha}dx
=\sum_{i'=1}^{N+p} \left(\mathbf{R}^{2h}_{h}\right)_{ii'}b^h_{i'}, \\
u^{2h}(x) = \sum_{i=1}^{N/2+p} u^{2h}_{i} \Lambda^{2h}_{i}=
\sum_{i'=1}^{N+p}\underbrace{\left[\sum_{i=1}^{N/2+p}
\left(\mathbf{R}^{2h}_{h}\right)_{ii'}
u^{2h}_{i}\right]}_{u^h_{i'}}\Lambda^{h}_{i'}(x) \Longrightarrow
\mathbf{u}^h = \left(\mathbf{R}^{2h}_{h}\right)^T\mathbf{u}^{2h} =
\mathbf{P}^h_{2h}\mathbf{u}^{2h}.\\
\end{gather*}
\section{Numerical results for the transfer operators}
The prolongation matrix as defined in Eq.~(\ref{eq:prolongation}) was
calculated using the BSPLINES module. A Gauss integration with
$N_G=\lceil p+(\alpha+1)/2 \rceil$ points is used to carry out the
numerical integrations.
In the following, the results are presented
for linear, quadratic and cubic Splines. Since the
restriction matrix is just the transpose of the prolongation matrix,
only the latter is shown. As expected, all the obtained matrices
are found to be \emph{independent} of $\alpha$ and \emph{sparse}.
During the calculations, it was checked that
\begin{itemize}
\item The coarse matrix computed using Eq.~(\ref{eq:coarse_mat}) and
the transfer matrix, is identical to the matrix assembled directly
on the coarse grid.
\item The sum of each row of the prolongation matrix is 1, since a
constant function ($\mathbf{u}^{2h}=1$) should remain constant after
the grid transfer.
\end{itemize}
\subsection{Linear Splines}
For $N=8$, the prolongation is a $9\times 5$ matrix given by
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
1 & 0 & 0 & 0 & 0 \\
1/2 & 1/2 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 \\
0 & 1/2 & 1/2 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 1/2 & 1/2 & 0 \\
0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1/2 & 1/2\\
0 & 0 & 0 & 0 & 1 \\
\end{matrix}\right)
\end{equation}
As expected, the prolongation matrix for linear Splines is
identical to the one obtained for first order FD discretization, where
a linear interpolation is used. One can easily check that
\begin{equation*}
\begin{split}
\Lambda^{2h}_1(x) &= \Lambda^h_1(x) + \frac{1}{2}\Lambda^h_2(x), \\
\Lambda^{2h}_2(x) &= \frac{1}{2}\Lambda^h_2(x) + \Lambda^h_3(x) +
\frac{1}{2}\Lambda^h_4(x), \\
\end{split}
\end{equation*}
as expected from (\ref{eq:restrict_gen}).
\subsection{Quadratic Splines}
For $N=8$, the prolongation is a $10\times 6$ matrix given by
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
1 & 0 & 0 & 0 & 0 & 0 \\
1/2 & 1/2 & 0 & 0 & 0 & 0 \\
0 & 3/4 & 1/4 & 0 & 0 & 0 \\
0 & 1/4 & 3/4 & 0 & 0 & 0 \\
0 & 0 & 3/4 & 1/4 & 0 & 0 \\
0 & 0 & 1/4 & 3/4 & 0 & 0 \\
0 & 0 & 0 & 3/4 & 1/4 & 0 \\
0 & 0 & 0 & 1/4 & 3/4 & 0 \\
0 & 0 & 0 & 0 & 1/2 & 1/2\\
0 & 0 & 0 & 0 & 0 & 1 \\
\end{matrix}\right)
\end{equation}
\subsection{Cubic Splines}
For $N=10$, the prolongation is a $13\times 8$ matrix given by
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1/2 & 1/2 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 3/4 & 1/4 & 0 & 0 & 0 & 0 & 0 \\
0 & 3/16 & 11/16 & 1/8 & 0 & 0 & 0 & 0 \\
0 & 0 & 1/2 & 1/2 & 0 & 0 & 0 & 0 \\
0 & 0 & 1/8 & 3/4 & 1/8 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/2 & 1/2 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/8 & 3/4 & 1/8 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/2 & 1/2 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/8 & 11/16 & 3/16 & 0 \\
0 & 0 & 0 & 0 & 0 & 1/4 & 3/4 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1/2 & 1/2\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
\end{matrix}\right)
\end{equation}
Note that from the results shown above, it is straightforward to derive the
prolongation matrix for other number of intervals $N$.
\section{Practical Considerations}
\subsection{Boundary conditions}
The \emph{essential Dirichlet boundary conditions} are imposed by zeroing the column
and row (first column and first row for the left boundary and last column
and last row for the right boundary) of the FE matrix $\mathbf{A}^h$ and
putting 1 on the diagonal. The same operation should be also performed
on the prolongation matrix, preserving thus the relation
(\ref{eq:coarse_mat}). For non-homogeneous Dirichlet boundary
conditions, the elements of the column should be saved before the
zeroing operation (for example $A^h_{21}, A^h_{31}, \ldots$ for the left
boundary condition). They will be used later to modify the right hand
side:
\begin{equation*}
b^h_i \leftarrow b^h_i-A^h_{i1}u^h_1, \quad i=2,\ldots
\end{equation*}
Nothing has to be done for \emph{natural boundary conditions}.
\subsection{Residual norm and error}
The residual norm is simply defined as the Euclidean norm of the
residue:
\begin{equation}
\label{eq:resid}
\|r\|_2 = \|\mathbf{b}-\mathbf{A}\mathbf{u}\|_2 =
\sqrt{\sum_i\left(b_i-\sum_{i'}A_{ii'}u_{i'}\right)^2}.
\end{equation}
When the \emph{exact} solution $u(x)$ is known, the \emph{discretization error}
can defined as
\begin{equation}
\label{eq:discerr}
\|e\|_2 = \sqrt{\int x^\alpha dx\left[\sum_{i}u_{i}\Lambda_i(x)-u(x)\right]^2}
\end{equation}
and computed using a Gauss quadrature. Note that for Splines of order
$p$, $\|e\|_2(h)$ converges to zero as $O(h^{p+1})$.
\section{The Model Problems}
\subsection{Cartesian geometry}
The following second-order boundary value problem is considered:
\begin{equation}
\label{eq:cartesian_problem}
\begin{split}
-\frac{d^2}{dx^2} u(x) + \sigma u(x) &= \sin (\pi kx), \qquad 0\le x\le 1 \\
u(0)=u(1) &= 0 \\
\Rightarrow u(x) = \frac{\sin(\pi kx)}{\pi^2k^2+\sigma}.& \\
\end{split}
\end{equation}
Using the weak form, the FE discretized matrix and right hand side can
be computed as
\begin{equation}
A_{ii'} = \int_0^1dx\left[\Lambda'_i(x)\Lambda'_{i'}(x) +
\sigma\Lambda_i(x)\Lambda_{i'}(x)\right], \qquad
b_i = \int_0^1dx \sin (\pi kx)\Lambda_i(x).
\end{equation}
For Splines of order $p$, the integration is done with a $\lceil
p+1/2\rceil$ point Gauss quadrature which is \emph{exact}
for the matrix $\mathbf{A}$ if $\sigma$ is constant.
The boundary conditions are simply imposed by setting
\begin{equation*}
A_{ki}=A_{ik}=\delta_{ik} \qquad\mbox{and} \qquad b_k=0
\end{equation*}
for $k=1$ (the first equation) and $N+p$ (the last equation).
\subsection{Cylindrical geometry}
The following second-order boundary value problem is considered:
\begin{equation}
\begin{split}
-\frac{1}{r}\frac{d}{dr}r\frac{d}{dr}u(r) + \frac{m^2}{r^2}u(r) &=
j^2_{ms}J_{m}(j_{ms}r), \qquad 0\le r\le 1, \quad j_{ms} =
s^{th}\mbox{ zero of }J_{m}, \\
u(1) &= 0 \\
\Rightarrow u(r) = J_{m}(j_{ms}r).& \\
\end{split}
\end{equation}
Using the weak form, the FE discretized matrix and right hand side can
be computed as
\begin{equation}
A_{ii'} = \int_0^1rdr\left[\Lambda'_i(r)\Lambda'_{i'}(r) +
\frac{m^2}{r^2}\Lambda_i(r)\Lambda_{i'}(r)\right], \qquad
b_i = \int_0^1rdr j^2_{ms}\,J_{m}(j_{ms}r)\Lambda_i(r).
\end{equation}
The boundary condition has only to be imposed on the last equation,
using the same procedure described for the Cartesian geometry.
It should be noted here that for $m\neq 0$, the matrix elements
$A_{1i}$ and $A_{i1}$ \emph{diverge} since $\Lambda_1(r)$ is not
equal to zero at $r=0$. However, using a \emph{direct solver}, one can observe
that the resulting \emph{discretization errors} as defined by
Eq.(\ref{eq:discerr}) converge for number of Gauss points $N_G$ slightly
larger than $p+1$, as shown in Table~\ref{tab:gauss_conv}. Then, using
$N_G=4$ and $6$ for the linear and cubic splines respectively, the
discretization error as a function of the number of grid intervals
(Fig~\ref{fig:cyl_conv}) show the expected quadratic and quartic
scaling respectively for the linear and cubic Splines.
\begin{figure}
\centering
\includegraphics[angle=0,width=0.8\hsize]{cyl_conv}
\caption{Discretization errors $\|e\|_2$ obtained by a \emph{direct
solver} versus the number of grid intervals $N$. A linear fit
yields a quadratic scaling ($\sim N^{-2.0}$) for the linear Splines
and a quartic convergence ($\sim N^{-4.3}$) for the cubic Splines.}
\label{fig:cyl_conv}
\end{figure}
\begin{table}
\centering
\begin{tabular}{|c|c|c|}\hline
Number of Gauss points & $p=1$ & $p=3$ \\ \hline
2 & 8.319E-04 & \\
4 & 9.277E-04 & 5.799E-07 \\
6 & 9.276E-04 & 5.936E-07 \\
8 & 9.276E-04 & 5.936E-07 \\\hline
\end{tabular}
\caption{Convergence of the \emph{discretization error} with respect
to the number of Gauss points for the cylindrical problem with $m=1$,
$s=10$ on a $128$ interval grid.}
\label{tab:gauss_conv}
\end{table}
\section{The Multigrid Schemes}
The two grid procedure described in section (\ref{sec:twogrid}) can be
generalized as follow.
Let $\nu_1$, $\nu_2$ and $\mu$ be three iteration parameters.
Given a guess $\mathbf{u}^h$ and right hand side $\mathbf{b}^h$ at the
\emph{finest} level, a MG cycle represented by
\begin{equation*}
\boxed{\mathbf{u}^h \leftarrow MG^h(\mathbf{u}^h,\mathbf{b}^h)}
\end{equation*}
will compute a \emph{new} $\mathbf{u}^h$ and is defined recursively by the
following steps:
\begin{enumerate}
\item If $h$ is the coarsest mesh size, direct solve
$\mathbf{A}^h\mathbf{u}^h=\mathbf{b}^h$ and return.
\item Else
\begin{itemize}
\item Relax $\mathbf{u}^h$ $\nu_1$ times.
\item $\mathbf{b}^{2h} \leftarrow
\mathbf{R}^{2h}_h(\mathbf{b}^h-\mathbf{A}^h\mathbf{u}^h), \quad
\mathbf{u}^{2h}\leftarrow 0$.
\item $\mathbf{u}^{2h} \leftarrow MG^{2h}(\mathbf{u}^{2h},\mathbf{b}^{2h})$ $\mu$ times.
\item $\mathbf{u}^h\leftarrow
\mathbf{u}^h+\mathbf{P}^{h}_{2h}\mathbf{u}^{2h}$.
\item Relax $\mathbf{u}^h$ $\nu_2$ times.
\end{itemize}
\end{enumerate}
The standard $V$-cycle is obtained for $\mu=1$ while $\mu=2$ results in
the $W$-cycle. Usually the number of \emph{pre-smooth} and
\emph{post-smooth} sweeps $\nu_1$ and $\nu_2$ is limited to 1 or 2. In the
following a $V$-cycle will be denoted by $V(\nu_1,\nu_2)$.
Another multigrid algorithm called \emph{Full Multigrid} or FMG does
not require an input guess $\mathbf{u}^h$ but solves first the
problem on coarser grids and uses one or many MG cycles to obtain the problem
solution. It can be represented by
\begin{equation*}
\boxed{\mathbf{u}^h \leftarrow FMG^h(\mathbf{b}^h)}
\end{equation*}
and defined recursively by the following steps:
\begin{enumerate}
\item If $h$ is the coarsest mesh size, direct solve
$\mathbf{A}^h\mathbf{u}^h=\mathbf{b}^h$ and return.
\item Else
\begin{itemize}
\item $\mathbf{b}^{2h} \leftarrow \mathbf{R}^{2h}_h(\mathbf{b}^h)$.
\item $\mathbf{u}^{2h} \leftarrow FMG^{2h}(\mathbf{b}^{2h})$.
\item $\mathbf{u}^h \leftarrow \mathbf{P}^{h}_{2h}\mathbf{u}^{2h}$.
\item $\mathbf{u}^h \leftarrow MG^h(\mathbf{u}^h,\mathbf{b}^h)$ $\nu_0$ times.
\end{itemize}
\end{enumerate}
Note that while the MG process is an iterative process (started by
setting for example the initial guess $\mathbf{u}^h=0$), the FMG is
more like a \emph{direct solver} with appropriate values of $\nu_1$,
$\nu_2$ and $\nu_0$ determined experimentally.
\section{Numerical Experiments}
The residual norm $\|r\|_2$ and error $\|e\|_2$ defined
previously are reported after each $V$-cycle in
Table~\ref{tab:cartesian1} for the Cartesian model problem and in
Table~\ref{tab:cylindrical1} for the cylindrical one . The ratio
between successive cycle $\|r\|_2$
and $\|e\|_2$ are shown in columns labeled \emph{ratio} and measure
the rate of iteration convergence. The \emph{asymptotic} ratio of
$\|r\|_2$ is called the \emph{convergence factor}.
In all the cases shown, one can note that $\|e\|_2$ level off quickly to the
discretization error obtained by using the \emph{direct solver} on the
finest grid, while the residual norms $\|r\|_2$ continue to decrease until
the machine zero is eventually reached. One can also verify that the \emph{final}
discretization errors scale approximately as $8^2$ and $8^4$
respectively for linear and cubic Splines, as $N$ is increased from
$128$ to $1024$.
Most interestingly, the \emph{iterative performance} depends very weakly
on the problem size $N$, for both the Cartesian and the cylindrical
cases. Moreover, the multigrid seems to be less efficient when linear Splines are used for
the problem discretization. This iterative performance can be further improved by
increasing the \emph{iteration parameters} $\nu_1$, $\nu_2$ and $\mu$,
as shown in Table~\ref{tab:improv}. One can also observe in the same
table that the Jacobi relaxation is systematically less efficient than Gauss
Seidel relaxation.
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline
\multicolumn{9}{|c|}{ Linear B-Splines $p=1$} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 6.219E-02 & & 7.164E-04 & & 2.210E-02 & & 7.164E-04 & \\
1 & 2.169E-02 & 0.35 & 5.880E-05 & 0.08 & 9.699E-03 & 0.44 & 3.622E-05 & 0.05 \\
2 & 3.801E-03 & 0.18 & 7.806E-06 & 0.13 & 1.790E-03 & 0.18 & 1.965E-06 & 0.05 \\
3 & 5.061E-04 & 0.13 & 3.666E-06 & 0.47 & 2.923E-04 & 0.16 & 1.583E-07 & 0.08 \\
4 & 6.762E-05 & 0.13 & 3.564E-06 & 0.97 & 4.055E-05 & 0.14 & 6.197E-08 & 0.39 \\
5 & 8.902E-06 & 0.13 & 3.585E-06 & 1.01 & 5.586E-06 & 0.14 & 5.655E-08 & 0.91 \\
6 & 1.199E-06 & 0.13 & 3.589E-06 & 1.00 & 7.122E-07 & 0.13 & 5.622E-08 & 0.99 \\
7 & 1.585E-07 & 0.13 & 3.590E-06 & 1.00 & 9.815E-08 & 0.14 & 5.620E-08 & 1.00 \\
8 & 2.089E-08 & 0.13 & 3.590E-06 & 1.00 & 1.320E-08 & 0.13 & 5.619E-08 & 1.00 \\
9 & 2.746E-09 & 0.13 & 3.590E-06 & 1.00 & 1.887E-09 & 0.14 & 5.619E-08 & 1.00 \\
10 & 3.741E-10 & 0.14 & 3.590E-06 & 1.00 & 2.533E-10 & 0.13 & 5.619E-08 & 1.00 \\
\hline \hline
\multicolumn{9}{|c|}{ Cubic B-Splines $p=3$} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 6.187E-02 & & 7.164E-04 & & 2.209E-02 & & 7.164E-04 & \\
1 & 1.948E-04 & 0.00 & 1.893E-06 & 0.00 & 1.685E-05 & 0.00 & 4.292E-08 & 0.00 \\
2 & 4.316E-06 & 0.02 & 3.927E-09 & 0.00 & 1.241E-07 & 0.01 & 7.156E-11 & 0.00 \\
3 & 1.554E-07 & 0.04 & 2.374E-09 & 0.60 & 4.184E-09 & 0.03 & 6.198E-13 & 0.01 \\
4 & 5.750E-09 & 0.04 & 2.373E-09 & 1.00 & 1.560E-10 & 0.04 & 5.635E-13 & 0.91 \\
5 & 2.153E-10 & 0.04 & 2.373E-09 & 1.00 & 5.912E-12 & 0.04 & 5.635E-13 & 1.00 \\
6 & 8.122E-12 & 0.04 & 2.373E-09 & 1.00 & 2.258E-13 & 0.04 & 5.635E-13 & 1.00 \\
7 & 3.079E-13 & 0.04 & 2.373E-09 & 1.00 & 8.777E-15 & 0.04 & 5.635E-13 & 1.00 \\
8 & 1.173E-14 & 0.04 & 2.373E-09 & 1.00 & 1.758E-15 & 0.20 & 5.635E-13 & 1.00 \\
9 & 4.489E-16 & 0.04 & 2.373E-09 & 1.00 & 1.709E-15 & 0.97 & 5.635E-13 & 1.00 \\
10 & 9.571E-17 & 0.21 & 2.373E-09 & 1.00 & 1.761E-15 & 1.03 &
5.635E-13 & 1.00 \\ \hline
\end{tabular}
\caption{The multigrid $V(1,1)$ performance with Gauss-Seidel relation
for a \emph{Cartesian} problem with $k=10$ and $\sigma=0$,
discretized on a grid with $N=128$ and $1024$ intervals, using
linear and cubic B-splines. For
both grid sizes, a total of 6 grid levels were considered.}
\label{tab:cartesian1}
\end{table}
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline
\multicolumn{9}{|c|}{ Linear B-Splines $p=1$} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 1.789E+01 & & 9.354E-02 & & 6.400E+00 & & 9.354E-02 & \\
1 & 3.373E+00 & 0.19 & 3.068E-03 & 0.03 & 1.826E+00 & 0.29 & 3.036E-03 & 0.03 \\
2 & 4.895E-01 & 0.15 & 8.064E-04 & 0.26 & 3.133E-01 & 0.17 & 1.624E-04 & 0.05 \\
3 & 6.160E-02 & 0.13 & 6.704E-04 & 0.83 & 4.581E-02 & 0.15 & 1.411E-05 & 0.09 \\
4 & 8.013E-03 & 0.13 & 6.811E-04 & 1.02 & 5.959E-03 & 0.13 & 1.062E-05 & 0.75 \\
5 & 9.871E-04 & 0.12 & 6.844E-04 & 1.00 & 8.098E-04 & 0.14 & 1.069E-05 & 1.01 \\
6 & 1.283E-04 & 0.13 & 6.847E-04 & 1.00 & 1.048E-04 & 0.13 & 1.070E-05 & 1.00 \\
7 & 1.613E-05 & 0.13 & 6.847E-04 & 1.00 & 1.504E-05 & 0.14 & 1.070E-05 & 1.00 \\
8 & 2.097E-06 & 0.13 & 6.847E-04 & 1.00 & 2.050E-06 & 0.14 & 1.070E-05 & 1.00 \\
9 & 2.639E-07 & 0.13 & 6.847E-04 & 1.00 & 3.008E-07 & 0.15 & 1.070E-05 & 1.00 \\
10 & 3.500E-08 & 0.13 & 6.847E-04 & 1.00 & 4.074E-08 & 0.14 & 1.070E-05 & 1.00 \\
\hline \hline
\multicolumn{9}{|c|}{ Cubic B-Splines $p=3$} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 1.768E+01 & & 9.354E-02 & & 6.399E+00 & & 9.354E-02 & \\
1 & 4.243E-02 & 0.00 & 4.727E-05 & 0.00 & 4.975E-03 & 0.00 & 6.588E-06 & 0.00 \\
2 & 1.378E-03 & 0.03 & 1.897E-06 & 0.04 & 7.835E-05 & 0.02 & 6.578E-09 & 0.00 \\
3 & 4.773E-05 & 0.03 & 1.814E-06 & 0.96 & 2.797E-06 & 0.04 & 4.125E-10 & 0.06 \\
4 & 2.174E-06 & 0.05 & 1.814E-06 & 1.00 & 1.041E-07 & 0.04 & 4.092E-10 & 0.99 \\
5 & 4.816E-07 & 0.22 & 1.814E-06 & 1.00 & 3.935E-09 & 0.04 & 4.092E-10 & 1.00 \\
6 & 1.942E-07 & 0.40 & 1.814E-06 & 1.00 & 1.499E-10 & 0.04 & 4.092E-10 & 1.00 \\
7 & 8.887E-08 & 0.46 & 1.814E-06 & 1.00 & 5.757E-12 & 0.04 & 4.092E-10 & 1.00 \\
8 & 4.449E-08 & 0.50 & 1.814E-06 & 1.00 & 2.517E-13 & 0.04 & 4.092E-10 & 1.00 \\
9 & 2.377E-08 & 0.53 & 1.814E-06 & 1.00 & 1.360E-13 & 0.54 & 4.092E-10 & 1.00 \\
10 & 1.328E-08 & 0.56 & 1.814E-06 & 1.00 & 1.384E-13 & 1.02 & 4.092E-10 & 1.00 \\ \hline
\end{tabular}
\caption{The multigrid $V(1,1)$ performance with Gauss-Seidel relation
for a one-dimensional \emph{cylindrical} problem with $m=22$ and
$s=10$, discretized on a grid with
$N=128$ and $1024$ intervals, using linear and cubic B-splines. For
both grid sizes, a total of 6 grid levels were considered.}
\label{tab:cylindrical1}
\end{table}
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|}\hline
& \multicolumn{2}{c|}{ Cartesian problem}
& \multicolumn{2}{c|}{ Cylindrical problem} \\ \cline{2-5}
& $N=128$ & $N=1024$ & $N=128$ & $N=1024$\\ \hline
$\nu_1=1, \nu_2=1, \mu=1$ & 0.13 & 0.14 & 0.13 & 0.14 \\
$\nu_1=1, \nu_2=2, \mu=1$ & 0.08 & 0.08 (\emph{0.10}) & 0.08 & 0.08 (\emph{0.09}) \\
$\nu_1=2, \nu_2=1, \mu=1$ & 0.08 & 0.08 & 0.08 & 0.08 \\
$\nu_1=2, \nu_2=2, \mu=1$ & 0.04 & 0.04 (\emph{0.08}) & 0.02 & 0.03 (\emph{0.08})\\
$\nu_1=1, \nu_2=1, \mu=2$ & 0.12 & 0.11 & 0.12 & 0.11 \\ \hline
\end{tabular}
\caption{The \emph{convergence factor} (averaged over
the last 5 cycles) for different iteration parameters $\nu_1$, $\nu_2$ and
$\mu$, using the linear Splines for both Cartesian ($k=10$,
$\sigma=0$) and cylindrical ($m=22$, $s=10$) problems. The last
entry is usually called
a $W$-cycle while the first four designate a $V(\nu_1,\nu_2)$
cycle. Gauss Seidel relaxation is used except for the results enclosed
in parenthesis which are obtained with the Jacobi
(weighted with $\omega=2/3$) relaxation.}
\label{tab:improv}
\end{table}
The next experiment is shown on Table~\ref{tab:fmg}, where two FMG$(\nu_1,\nu_2)$
schemes are applied to the $m=22$, $s=10$ cylindrical problem with
grid sizes up to $N=2048$. Note that the problem is solved to the level of
discretization for $N\ge 128$ with FMG$(2,1)$ but not with FMG$(1,1)$. Solving the same
problem with the $V(2,1)$ cycle required 3 iterations for all the values of
$N$ shown. Since the cost of one FMG(2,1) is $\sim 2$ the cost
of one $V(2,1)$ (see Appendix \ref{sec:cost}), it appears that FMG is
more efficient for $N\ge 128$.
Finally, in all the cases shown here, the equality (\ref{eq:coarse_mat})
is verified numerically, except for the cylindrical case with $m\neq 0$. This is
expected since as noted earlier, the matrix elements $A_{i1}$ and
$A_{1i}$ diverge unless $m=0$ in the cylindrical problem.
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|}\hline
& \multicolumn{2}{c|}{FMG(1,1)}
& \multicolumn{2}{c|}{FMG(2,1)} \\ \cline{2-5}
$N$ & $\|e\|_2$ & $\|e\|_2/\|e\|_d$ & $\|e\|_2$ & $\|e\|_2/\|e\|_d$ \\ \hline
4 & 1.011E-01 & 0.968 & 1.012E-01 & 0.969 \\
8 & 7.781E-02 & 1.031 & 7.679E-02 & 1.018 \\
16 & 3.332E-02 & 1.310 & 2.808E-02 & 1.104 \\
32 & 1.516E-03 & 1.421 & 1.098E-03 & 1.030 \\
64 & 5.168E-05 & 1.443 & 3.652E-05 & 1.019 \\
128 & 2.012E-06 & 1.109 & 1.818E-06 & 1.002 \\
256 & 1.125E-07 & 1.053 & 1.069E-07 & 1.001 \\
512 & 6.819E-09 & 1.037 & 6.576E-09 & 1.000 \\
1024 & 4.224E-10 & 1.032 & 4.093E-10 & 1.000 \\
2048 & 2.634E-11 & 1.031 & 2.556E-11 & 1.000 \\
\hline
\end{tabular}
\caption{The discretization errors $\|e\|_2$ obtained from a
FMG$(\nu_1,\nu_2)$ sweep with $\nu_0=1$ for different grid sizes $N$. The columns
$\|e\|_2/\|e\|_d$ display their ratio with the discretization errors
obtained from a \emph{direct} solver. The cylindrical problem with
$m=22$ and $s=10$ using cubic Splines is considered here.}
\label{tab:fmg}
\end{table}
\section{Periodic Case}
\subsection{Transfer operators}
For periodic problems, we use \emph{periodic} Splines \cite{BSPLINES} which
satisfy the periodic boundary condition
$\Lambda^h_{i+N}(x)=\Lambda^h_i(x-Nh)$. As a result, both the expansion
coefficients and the right hand sides are periodic
with periodicity $N$ ($u^h_{i+N}=u^h_i$, $b^h_{i+N}=b^h_i$) and he rank
of all matrices should be $N$ instead of $N+p$ as in the non-periodic case.
The \emph{prolongation} matrix $\mathbf{P}^h_{2h}$ as given by
(\ref{eq:prolongation}) are computed numerically and the results
for $N=8$ are given below for linear, quadratic and cubic Splines.
\begin{itemize}
\item Linear Splines
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
1 & 0 & 0 & 0 \\
1/2 & 1/2 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 1/2 & 1/2 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 1/2 & 1/2 \\
0 & 0 & 0 & 1 \\
1/2 & 0 & 0 & 1/2 \\
\end{matrix}\right)
\end{equation}
\item Quadratic Splines
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
3/4 & 1/4 & 0 & 0 \\
1/4 & 3/4 & 0 & 0 \\
0 & 3/4 & 1/4 & 0 \\
0 & 1/4 & 3/4 & 0 \\
0 & 0 & 3/4 & 1/4 \\
0 & 0 & 1/4 & 3/4 \\
1/4 & 0 & 0 & 3/4 \\
3/4 & 0 & 0 & 1/4 \\
\end{matrix}\right)
\end{equation}
\item Cubic Splines
\begin{equation}
\mathbf{P}^{h}_{2h} =
\left(
\begin{matrix}
1/2 & 1/2 & 0 & 0 \\
1/8 & 3/4 & 1/8 & 0 \\
0 & 1/2 & 1/2 & 0 \\
0 & 1/8 & 3/4 & 1/8 \\
0 & 0 & 1/2 & 1/2 \\
1/8 & 0 & 1/8 & 3/4 \\
1/2 & 0 & 0 & 1/2 \\
3/4 & 1/8 & 0 & 1/8 \\
\end{matrix}\right)
\end{equation}
\end{itemize}
The restriction matrix is simply $\mathbf{R}^{2h}_h=
(\mathbf{P}^{h}_{2h})^T$. Generalization for any other number of intervals $N$ should be
straightforward.
\subsection{Numerical Experiments}
In order to test the grid transfer operators obtained above, the same
second-order problem (\ref{eq:cartesian_problem}) but with the
periodic boundary condition $u(x+1)=u(x)$ is considered. It should be
noted that in that case, if $\sigma=0$, the problem is singular since
the solution is not \emph{unique}! But we have observed that this
problem can be avoided for a slightly non zero $\sigma$,
With $\sigma=0.01$ and $k=10$ and using linear and cubic Splines, we
recover the same multigrid iterative performances shown in Table
\ref{tab:cartesian1} obtained previously for non-periodic Dirichlet
boundary conditions. Table \ref{tab:quad_splines} also shows similar
iterative efficiencies for \emph{quadratic} non-periodic and periodic problems.
The identity (\ref{eq:coarse_mat}) is numerically verified in all the
cases considered.
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline
\multicolumn{9}{|c|}{ Cartesian problem with quadratic splines} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 6.203E-02 & & 7.164E-04 & & 2.209E-02 & & 7.164E-04 & \\
1 & 8.114E-04 & 0.01 & 6.375E-06 & 0.01 & 1.003E-04 & 0.00 & 4.509E-07 & 0.00 \\
2 & 1.891E-05 & 0.02 & 6.079E-08 & 0.01 & 1.769E-06 & 0.02 & 8.061E-10 & 0.00 \\
3 & 1.103E-06 & 0.06 & 5.220E-08 & 0.86 & 7.018E-08 & 0.04 & 9.970E-11 & 0.12 \\
4 & 8.148E-08 & 0.07 & 5.220E-08 & 1.00 & 5.620E-09 & 0.08 & 9.958E-11 & 1.00 \\
5 & 6.368E-09 & 0.08 & 5.220E-08 & 1.00 & 4.772E-10 & 0.08 & 9.958E-11 & 1.00 \\
6 & 4.969E-10 & 0.08 & 5.220E-08 & 1.00 & 4.101E-11 & 0.09 & 9.958E-11 & 1.00 \\
7 & 3.874E-11 & 0.08 & 5.220E-08 & 1.00 & 3.548E-12 & 0.09 & 9.958E-11 & 1.00 \\
8 & 3.081E-12 & 0.08 & 5.220E-08 & 1.00 & 3.081E-13 & 0.09 & 9.958E-11 & 1.00 \\
9 & 2.489E-13 & 0.08 & 5.220E-08 & 1.00 & 2.690E-14 & 0.09 & 9.958E-11 & 1.00 \\
10 & 1.986E-14 & 0.08 & 5.220E-08 & 1.00 & 3.212E-15 & 0.12 & 9.958E-11 & 1.00 \\
\hline \hline
\multicolumn{9}{|c|}{ Periodic problem with quadratic splines} \\ \hline
& \multicolumn{4}{c|}{ $N=128$}
& \multicolumn{4}{c|}{ $N=1024$} \\ \cline{2-9}
$V$-cycle & $\|r\|_2$ & ratio & $\|e\|_2$ & ratio &
$\|r\|_2$ & ratio & $\|e\|_2$ & ratio \\ \hline
0 & 6.203E-02 & & 7.164E-04 & & 2.209E-02 & & 7.164E-04 & \\
1 & 1.285E-03 & 0.02 & 1.294E-05 & 0.02 & 3.116E-04 & 0.01 & 5.862E-06 & 0.01 \\
2 & 7.878E-05 & 0.06 & 6.569E-07 & 0.05 & 2.893E-05 & 0.09 & 1.626E-07 & 0.03 \\
3 & 6.573E-06 & 0.08 & 6.511E-08 & 0.10 & 2.691E-06 & 0.09 & 5.631E-09 & 0.03 \\
4 & 5.681E-07 & 0.09 & 5.224E-08 & 0.80 & 2.385E-07 & 0.09 & 2.400E-10 & 0.04 \\
5 & 4.890E-08 & 0.09 & 5.219E-08 & 1.00 & 2.097E-08 & 0.09 & 9.997E-11 & 0.42 \\
6 & 4.198E-09 & 0.09 & 5.219E-08 & 1.00 & 1.828E-09 & 0.09 & 9.958E-11 & 1.00 \\
7 & 3.607E-10 & 0.09 & 5.219E-08 & 1.00 & 1.584E-10 & 0.09 & 9.958E-11 & 1.00 \\
8 & 3.103E-11 & 0.09 & 5.219E-08 & 1.00 & 1.370E-11 & 0.09 & 9.958E-11 & 1.00 \\
9 & 2.674E-12 & 0.09 & 5.219E-08 & 1.00 & 1.184E-12 & 0.09 & 9.958E-11 & 1.00 \\
10 & 2.307E-13 & 0.09 & 5.219E-08 & 1.00 & 1.025E-13 & 0.09 & 9.958E-11 & 1.00 \\\hline
\end{tabular}
\caption{The multigrid $V(1,1)$ performance with Gauss-Seidel relation
for \emph{Cartesian} problem ($k=10$, $\sigma=0$) and
\emph{periodic} problem ($k=10$, $\sigma=0.01$), discretized on a grid with
$N=128$ and $1024$ intervals, using quadratic B-splines. For
both grid sizes, a total of 6 grid levels were considered.}
\label{tab:quad_splines}
\end{table}
\section{Conclusion}
Using the variational principle, we have derived the expressions of
the grid transfer matrices for Finite Elements using Splines of any
order. It is found that:
\begin{itemize}
\item The grid transfer matrices do not depend of the geometries
characterized by the Jacobian as defined in $dV = x^\alpha dx$.
\item The standard grid transfer operator used for first order finite
difference (FD) discretization for Cartesian geometry is recovered when
linear Spline finite elements (FE) are used.
\item Applying these transfer matrices, we have solved Cartesian,
cylindrical as well as periodic one dimensional
problems, and obtained essentially the same multigrid iterative performances
as found for standard first order FD Cartesian problems.
\item No performance \emph{degradation} is observed when the order of Splines FE
is increased from 1 to 3, or when the cylindrical geometry is considered.
\end{itemize}
For two dimensional problems, notice that for both Cartesian ($dV=dxdy$) and
standard curvilinear geometries ($dV=r^\alpha drd\theta$), the
Jacobian is \emph{separable}. Using this property, one can show that
the two dimensional grid transfer consists of simply applying
successively one dimensional grid transfer on each of the $x$ and
$y$ (or $r$ and $\theta$) grids. With the solution
$\mathbf{u}^h=[u^h_{ij}]$ and right hand side
$\mathbf{b}^h=[b^h_{ij}]$ defined by
\begin{equation}
u(x,y) = \sum_{ij} u^h_{ij}\Lambda^h_i(x)\Lambda^h_j(y), \qquad
b^h_{ij} = \int dx\Lambda^h_i(x)\int dy\Lambda^h_j(y)f(x,y),
\end{equation}
the two dimension grid transfers can be expressed as (see Appendix \ref{sec:twod})
\begin{equation}
\begin{split}
\mathbf{u}^h &= {_{x}\mathbf{P}^h_{2h}}\; \mathbf{u}^{2h}
\left(_{y}\mathbf{P}^h_{2h}\right)^T, \\
\mathbf{b}^{2h} &={_{x}\mathbf{R}^{2h}_{h}}\; \mathbf{b}^{h}
\left(_{y}\mathbf{R}^{2h}_{h}\right)^T. \\
\end{split}
\end{equation}
For more general curvilinear coordinates such as found in tokamak
magnetic coordinates defined by $dV=J(s,\theta)dsd\theta$, we will
assume that the grid transfer operators derived above are still
applicable. The validity of this assumption will be the object of the
next task.
\appendix
\section{Multigrid Cost Estimation}
\label{sec:cost}
Assuming that the \emph{coarsest} grid is fixed to $2$, the total
number of grid levels $L$ is given by $N/2^{L-1}=2$ or $L=\log_2(N)$, where $N$
is the number of intervals in the \emph{finest grid}. Since both
\emph{relaxation} and intergrid transfer are proportional to the
number of problem unknowns, the cost of the
$V$-cycle can be estimated as:
\begin{equation}
\begin{split}
\mbox{MG}(N) &= c\left[ (N+p)+(N/2+p) +\ldots + (N/2^{L-2}+p) \right] \\
&= c \left[ 2N-4 +(L-1)p \right] ,\\
\end{split}
\end{equation}
where $p$ is the order of Splines used for the discretization. The FMG
can then be deduced, assuming $\nu_0=1$ as
\begin{equation}
\begin{split}
\mbox{FMG}(N) &= \mbox{MG}(N) +\mbox{MG}(N/2) + \ldots + \mbox{MG}(N/2^{L-2})\\
&= c\left[ 4N-8 +(L-1)(pL/2-4) \right] .\\
\end{split}
\end{equation}
As expected a single FMG cycle (with $\nu_0=1$) costs about two $V$-cycles.
\section{Two dimensional Grid Transfer}
\label{sec:twod}
On the fine and the coarse grids, the problem solution $u(x,y)$ can be written
as:
\begin{equation*}
u(x,y) = \sum_{i'j'} u^h_{i'j'}\Lambda^h_{i'}(x)\Lambda^h_{j'}(y)
= \sum_{i'j'} u^{2h}_{i'j'}\Lambda^{2h}_{i'}(x)\Lambda^{2h}_{j'}(y).
\end{equation*}
Projecting these two expansions on the two dimensional basis
functions $\Lambda^h_{i}(x)\Lambda^h_{j}(y)$ yields
\begin{gather*}
\sum_{i'j'} u^h_{i'j'}
\underbrace{\int dx\Lambda^h_{i}(x)\Lambda^h_{i'}(x)}_{M^{h,h}_{ii'}}
\underbrace{\int dy\Lambda^h_{j}(x)\Lambda^h_{j'}(y)}_{N^{h,h}_{jj'}}
= \sum_{i'j'} u^{2h}_{i'j'}
\underbrace{\int dx\Lambda^h_{i}(x)\Lambda^{2h}_{i'}(x)}_{M^{h,2h}_{ii'}}
\underbrace{\int
dy\Lambda^h_{j}(x)\Lambda^{2h}_{j'}(y)}_{N^{h,2h}_{jj'}} \\
\Longrightarrow \quad
\mathbf{M}^{h,h}\;\mathbf{u}^h\;\left(\mathbf{N}^{h,h}\right)^T =
\mathbf{M}^{h,2h}\mathbf{u}^{2h}\left(\mathbf{N}^{h,2h}\right)^T \\
\Longrightarrow \quad
\mathbf{u}^h = \left(\mathbf{M}^{h,h}\right)^{-1}
\mathbf{M}^{h,2h}\;\mathbf{u}^{2h}\;\left[\left(\mathbf{N}^{h,h}\right)^{-1}\mathbf{N}^{h,2h}\right]^T.
\end{gather*}
The right hand side can be written on the fine and coarse grids as
\begin{equation*}
\begin{split}
b^{h}_{ij} = \int dx\Lambda^{h}_i(x)\int dy\Lambda^{h}_j(y)\;f(x,y)
= \sum_{i'j'} M^{h,h}_{ii'}f^{h}_{i'j'}N^{h,h}_{jj'} \quad\Longrightarrow\quad
\mathbf{b}^h &= \mathbf{M}^{h,h}\;\mathbf{f}^h\left(\mathbf{N}^{h,h}\right)^{T}, \\
b^{2h}_{ij} = \int dx\Lambda^{2h}_i(x)\int dy\Lambda^{2h}_j(y)\;f(x,y)
= \sum_{i'j'} M^{2h,h}_{ii'}f^{h}_{i'j'}N^{2h,h}_{jj'} \quad\Longrightarrow\quad
\mathbf{b}^{2h} &= \mathbf{M}^{2h,h}\;\mathbf{f}^h\left(\mathbf{N}^{2h,h}\right)^{T}, \\
\end{split}
\end{equation*}
where the expansion of $f(x,y)$ on the \emph{fine} mesh has been
used. Elimination of $\mathbf{f}^h$ then yields
\begin{equation*}
\mathbf{b}^{2h} =
\mathbf{M}^{2h,h}\left(\mathbf{M}^{h,h}\right)^{-1} \;
\mathbf{b}^h\;
\left[\mathbf{N}^{2h,h}\left(\mathbf{N}^{h,h}\right)^{-1} \right]^T.
\end{equation*}
\begin{thebibliography}{99}
\bibitem{SOLVERS} {\tt The Solvers in BSPLINES},
\url{https://crppsvn.epfl.ch/repos/bsplines/trunk/docs/solvers.pdf}
\bibitem{BSPLINES} {\tt BSPLINES} Reference Guide,
\url{https://crppsvn.epfl.ch/repos/bsplines/trunk/docs/bsplines.pdf}
\bibitem{Briggs} {W.L.~Briggs, V.E.~Henson and S.F.~McCormick, A
Multigrid Tutorial, Second Edition, Siam (2000)}.
\end{thebibliography}
\end{document}
c4science · Help