%%%%%%%%%%%%%%%%%%%%%%%%% -*-LaTeX-*- %%% multigrid.tex %%% %%%%%%%%%%%%%%%%%%%%%%%%% %% %% Date: 21-Jun-2007 %% Author: wd (Wolfgang.Dobler@ucalgary.ca) %% Description: %% %% Choose appropriate driver: \RequirePackage{ifpdf} \ifpdf \def\mydriver{pdftex} % writing PDF \else \def\mydriver{dvips} % writing DVI \fi \documentclass[\mydriver,12pt,twoside,notitlepage,letterpaper]{article} %\usepackage{german,a4} %\usepackage[german,british]{babel} %\usepackage[latin1]{inputenc} %\usepackage[T1]{fontenc}\usepackage{lmodern} \usepackage{amsmath,amssymb,bm,wasysym} %\usepackage{latexsym,exscale}%\usepackage{newcent,helvet}%\usepackage{pxfonts}%\usepackage{txfonts} %\usepackage[it,footnotesize]{caption2} %\setlength{\abovecaptionskip}{5pt} % Space before caption %\setlength{\belowcaptionskip}{5pt} % Space after caption %% Load titlesec before hyperref or \section links will be off with dvipdfm %\usepackage[bf,sf,small,nonindentfirst]{titlesec} %\newcommand{\sectionbreak}{\clearpage} %\titleformat{\subsubsection}{\normalfont\itshape}{\thesubsubsection}{.5em}{} %\titlespacing{\subsubsection}{0pt}{*1{*-1} %\usepackage{graphicx} \usepackage{multicol} \usepackage{parskip} \usepackage[hmargin=25mm,top=20mm,bottom=20mm]{geometry} \usepackage{makeidx} \usepackage[bookmarksopen=true,bookmarksopenlevel=2]{hyperref} % should come late/last \ifpdf \usepackage{pdfpages} \else \newcommand{\includepdf}[2][]{% \centerline{% \fbox{% \parbox[c][0.5\textheight]{0.5\textwidth}{% \vspace*{\stretch{1}}% \hfill Include file \path{#2}\hfill\mbox{}% \vspace*{\stretch{2}}% }% }% }% }% \fi \frenchspacing \sloppy %%% Multiple, three-column indexes \makeindex %% The following is adapted from hyperref.sty and fixes hyperrefs in the %% index after all of our nasty manipulations: % \makeatletter % \@ifpackageloaded{hyperref}{% % \let\HyInd@org@wrindex\@wrindex % \def\@wrindex#1#2{\HyInd@@wrindex{#1}#2||\\}% % \def\HyInd@@wrindex#1#2|#3|#4\\{% % \ifx\\#3\\% % \HyInd@org@wrindex{#1}{#2|hyperpage}% % \else % \def\Hy@temp@A{#3}% % \ifx\Hy@temp@A\HyInd@ParenLeft % HyInd@org@wrindex{#1}{#2|#3hyperpage}% % \else % \HyInd@org@wrindex{#1}{#2|#3}% % \fi % \fi % }% % }{} % \makeatother %% Redefine index to be in three columns (adapted from `index.sty'): \makeatletter \renewenvironment{theindex}{% %\edef\indexname{\the\@nameuse{idxtitle@\@indextype}}% \if@twocolumn\@restonecolfalse \else\@restonecoltrue \fi \columnseprule \z@ \columnsep 35\p@ \begin{multicols}{3}[\section*{\indexname}% %\ifx\index@prologue\@empty% %\else\index@prologue\bigskip %\fi ]% \@mkboth{\MakeUppercase\indexname}% {\MakeUppercase\indexname}% \thispagestyle{plain}% \parindent\z@ \parskip\z@ \@plus .3\p@\relax \let\item\@idxitem } { \end{multicols} \if@restonecol\onecolumn\else\clearpage\fi } \makeatother %% Layout macros \newcommand{\righteqn}[1]{\llap{$\displaystyle #1$}} \newcommand{\centereqn}[1]{\clap{$\displaystyle #1$}} % Like \rlap, \llap, but centered % (http://math.arizona.edu/~aprl/publications/mathclap/perlis_mathclap_24Jun2003.pdf) \newcommand{\clap}[1]{\hbox to 0pt{\hss#1\hss}} %% Matrices and tensors \newcommand{\matx}[1]{\mbox{\boldmath${\cal #1}$\unboldmath}} \newcommand{\tensor}[1]{\matx{#1}} %% Math macros \newcommand{\Av} {\mathbf{A}} \newcommand{\Bv} {\mathbf{B}} \newcommand{\Ev} {\mathbf{E}} \newcommand{\fv} {\mathbf{f}} \newcommand{\gv} {\mathbf{g}} \newcommand{\jv} {\mathbf{j}} \newcommand{\Laplace} { \mathop{\Delta}\nolimits} \newcommand{\mtildei} {\widetilde{m}_{\rm i}} \newcommand{\uv} {\mathbf{u}} \newcommand{\uve} {\uv_{\rm e}} \newcommand{\uvi} {\uv_{\rm i}} \newcommand{\vv} {\mathbf{v}} \newcommand{\vve} {\vv_{\rm e}} \newcommand{\xv}{\mathbf{x}} \title{Multigrid in the Pencil Code} \author{Wolfgang Dobler} % ---------------------------------------------------------------------- % \begin{document} \thispagestyle{empty} \maketitle \tableofcontents %%%%%%%% \section{Introduction} %%%%%%%% The main purpose of this set of notes is to document the implementation of multigrid solvers in the Pencil Code. At the same time, I use it as a reference for the many confusing terms, a pointer to literature, and as a place for ideas and things that have not been implemented so far. %%%%%%%% \section{Glossary and more} %%%%%%%% To me, one of the most difficult aspects of multigrid methods is the jargon that is used in the literature. For some straight-forward concepts and operations, there are several opaque synonyms, which makes my mind spin, precess, and nutate at the same time. \footnote{ I admit to adding yet my own jargon here by introducing the terms ``coarse-graining'' and ``fine-graining''. But I think \emph{coarse-graining} is a way more intuitive term than \emph{restricting}, \emph{relaxing} [???], etc.} So, as a reference for myself and others in the same situation, we start this document (after a very superficial outline of the multigrid methodology) with a glossary that will also introduce all relevant concepts of multigrid. %%%%%%%%%%% \subsection{A bird's view of multigrid} A multigrid solver for an elliptic PDE includes the following steps: \begin{enumerate} \item Apply one or more Gauss--Seidel (or weighted Jacobi) iterations to solve the discretized equations (converges quickly on the smallest scales). \item Calculate the residual (which measures to what extent the current approximation does not satisfy the equations and should thus decrease in the process). \item Coarse-grain the operator and residual to construct a coarse-grained equation for the correction to the current approximation. \item Apply multigrid on the coarser grid (the algorithm is recursive). \item Fine-grain the correction and apply it on the finer grid. \item Repeat this a number of times, typically with varying depths of recursion. \end{enumerate} To be somewhat specific, let us consider the inhomogeneous nonuniform Helmholtz equation (or generalized Poisson equation) \begin{equation} [\Delta - h(\xv)] f(\xv) = g(\xv) \; , \end{equation} which can be discretized as \begin{eqnarray} &&\centereqn{ \dfrac{f_{i{-}1,j,k} - 2 f_{ijk} + f_{i{+}1,j,k}}{\delta x^2} + \dfrac{f_{i,j{-}1,k} - 2 f_{ijk} + f_{i,j{+}1,k}}{\delta y^2} + \dfrac{f_{i,j,k{-}1} - 2 f_{ijk} + f_{i,j,k{+}1}}{\delta z^2} } \nonumber \\ && \label{Eq-elliptic-discrete} {} - h_{ijk} f_{ijk} = g_{ijk} \; . \end{eqnarray} Equation~(\ref{Eq-elliptic-discrete}) can be solved for $f_{ijk}$ to give \begin{equation} \label{Eq-discrete-explicit} f_{ijk} = \mathcal{L}_{ijk} f \; , \end{equation} with the linear operator \begin{equation} \mathcal{L}_{ijk} f = \dfrac{\dfrac{f_{i{-}1,j,k} + f_{i{+}1,j,k}}{\delta x^2} + \dfrac{f_{i,j{-}1,k} + f_{i,j{+}1,k}}{\delta y^2} + \dfrac{f_{i,j,k{-}1} + f_{i,j,k{+}1}}{\delta y^2} - g_{ijk} } {\dfrac{2}{\delta x^2} + \dfrac{2}{\delta y^2} + \dfrac{2}{\delta z^2} + h_{ijk} } \; . \end{equation} %%%%%%%%%%% \subsection{Glossary} \begin{description} \item[Gauss--Seidel iteration:] \index{Gauss--Seidel} Method to solve the discretized system (\ref{Eq-elliptic-discrete}) iteratively. One \emph{Gauss--Seidel} smoothing step loops over $i,j,k$ and applies the iteration \begin{eqnarray} \label{Eq-Gauss-Seidel} && \mathtt{for\ (i,j,k)=}(N_x,N_y,N_z) \nonumber\\ && \quad f_{ijk} \leftarrow \mathcal{L}_{ijk} f\\ && \mathtt{endfor} \nonumber \end{eqnarray} where the right-hand-side is always evaluated using the latest values of $f$ available. Gauss--Seidel is faster than (unweighted) Jacobi and uses less memory. [$\rightarrow$ Matrix representation!] \item[Red-black Gauss--Seidel iteration:] \index{red-black Gauss--Seidel}% \index{Gauss--Seidel!red-black} A variant of Gauss--Seidel where the grid is divided into two (or more) ``colours'', based on whether $i+j+k \mod 2$ is $0$ or $1$ (think of a Go board or chequerboard). May be cargo cult, as I haven't seen clear evidence in the literature that this is more efficient than plain Gauss--Seidel. [$\rightarrow$ Matrix representation!] \item[Jacobi iteration:] \index{Jacobi} Similar to Gauss--Seidel iteration, but does not use the continuously updated values on the right-hand side: \begin{eqnarray} \label{Eq-Jacobi} && f^{\rm(old)} \leftarrow f \nonumber\\ && \mathtt{for\ (i,j,k)=}(N_x,N_y,N_z) \nonumber\\ && \quad f_{ijk} \leftarrow \mathcal{L}_{ijk} f^{\rm(old)} \\ && \mathtt{endfor} \nonumber\\ \end{eqnarray} \begin{itemize} \item Slower than Gauss-Seidel (need to iterate about twice as long) and uses more memory. \item Does not damp Nyquist signal. \item Popular in the literature, mainly because it is easier to analyze than Gauss--Seidel (symmetric iteration matrix). \end{itemize} \item[Weighted Jacobi iteration:] \index{weighted Jacobi} \index{Jacobi!weighted} A modification of Jacobi iteration that damps at the Nyquist scale: \begin{eqnarray} \label{Eq-weighted-Jacobi} && f^{\rm(old)} \leftarrow f \nonumber\\ && \mathtt{for\ (i,j,k)=}(N_x,N_y,N_z) \nonumber\\ && \quad f^{\rm(new)}_{ijk} \leftarrow \mathcal{L}_{ijk} f^{\rm(old)} \\ && \mathtt{endfor} \nonumber\\ && f = (1{-}w) f^{\rm(old)} + w f^{\rm(new)} \nonumber \end{eqnarray} Typically, the weighting coefficient is chosen to be \begin{equation} w = \dfrac{2}{3} \; . \end{equation} \item[Restricting:] Same as \emph{coarse-graining} \index{restricting} \item[Coarse-graining:] \index{coarse-graining} \hypertarget{coarse-graining} Mapping values from the given grid to the next coarser one. For face-centred multigrid, the coarser grid consists of points that are also on the finer grid, \begin{equation} \tilde{x}_i = x_{2i} \end{equation} Mathematically, coarse-graining needs to strongly damp any wavenumbers between the coarser-grid Nyquist wavenumber $\tilde{k}_{\rm Ny}=k_{\rm Ny}/2$ and the fine-grid Nyquist wavenumber $k_{\rm Ny}$, because those will otherwise get aliased as spurious low-wavenumber signals. For cell-centred multigrid, the coarse-grid points are \emph{not} points on the fine grid (except for the boundary points). One can generalize the approach that leads to full weighting (see below) and construct e.g.~a filter that uses two neighbouring points on either side and interpolates using a Nyquist-removing second-order polynomial. However, my experiments indicate that using fully-weighted interpolation to every point of the fine grid and then trilinear interpolation to the coarse grid points leads to faster convergence. \item[Fully weighted restriction:] \index{fully weighted restriction} \hypertarget{fully-weighted} A popular coarse-graining filter: \begin{equation} \label{Eq-fully-weighted} \tilde{f}_{i} = \dfrac{f_{2i-1} + 2 f_{2i} + f_{2i+1}} {4} \end{equation} Fully weighted restriction completely filters out $k_{\rm Ny}$ and damps $\tilde{k}_{\rm Ny}=k_{\rm Ny}/2$ by a factor $1/2$: \centerline{\includegraphics[width=0.6\textwidth]{fully_weighted}} \item[Refining:] Same as \emph{fine-graining} \index{refining} \item[Fine-graining:] \index{fine-graining} Interpolation of values (typically the correction to the current approximation) from the coarser to the finer grid. Typically implemented as trilinear interpolation. I have experimented with splines (in 1D) and did not find good evidence for an improvement. \item[Smoothing:] \index{smoothing} Damping the small-scale errors by applying a Gauss--Seidel (or similar) iteration step. \item[Relaxation:] \index{relaxation} Damping the large-scale errors by going to the coarser grid (and recursively higher up) and smoothing there. \item[Vertex-centred vs. cell-centred multigrid:] Classical discussions of multigrid restriction (=\hyperlink{coarse-graining}{coarse-graining}), and in particular of \hyperlink{fully-weighted}{fully weighted restriction} assume that every point of the coarse grid is also a point of all finer grids (``vertex-centered multigrid''). For the Pencil Code, this condition holds for a periodic grid, but not if we have boundaries (which implies ``cell-centered multigrid'') \addvspace{1\baselineskip} \centerline{\includegraphics[width=0.8\textwidth]{vertex-cell-centered}} To deal with this, the current (August 2007) implementation first smoothes using Eq.~(\ref{Eq-fully-weighted}) on the original (finer) grid, and then interpolates that smoothed (i.e.~low-pass filtered) version to the coarser grid points using tri-linear interpolation. \bigskip Suggestion by Chao-Chin Yang: Instead of insisting that the coarser-grid boundary points are identical with the main-grid ones, we could use coarser grids like this: \addvspace{1\baselineskip} \centerline{\includegraphics[width=0.8\textwidth]{vertex-cell-centered2}} (The left is Chao-Chin's suggestion, the right is yet another option). Advantage of the currently implemented scheme: will also work for $N_x = 390$. Disadvantage: currently uses significantly more memory (the fully-weighted field is stored on the finest grid). Could be improved by `chunking', i.e.~doing the ffully-weighted restriction only on a data sub-cube that then gets recycled. It might also be possible to directly formulate an expression for the result of full weighting and interpolating together. \item[Trilinear interpolation:] \index{trilinear interpolation} [Formula\ldots] \item[Adaptive multigrid:] \item[Algebraic multigrid:] \item[Semi-coarsening:] \index{semi-coarsening} If $\delta x \gg \delta y$, the operator is effectively treating the different coordinate directions very differently. [The same can also happen for equal grid spacing when the operator itself is highly anisotropic.] This can make the multigrid convergence very slow, and one way out is to only coarsen the grid in the finer direction (i.e.~$y$ in the example above). This should be pretty easy to implement: at every level, compare $\delta x$ and $\delta y$, and if $\delta x > 2\delta y$, do not coarse-grain in $x$. Another approach is \emph{line relaxation} or \emph{plane relaxation}. \item[Line relaxation:] \index{line relaxation} \item[Plane relaxation:] \index{plane relaxation} Another solution to the problem of a strongly anisotropic operator (see \emph{semi-coarsening} above). [Cite paper, discuss\ldots] \item[FDA(?):] \item[V-cycle:] \index{V cycle} \item[W-cycle:] \index{W cycle} \item[FMG cycle:] \index{FMG cycle} \item[Variational property:] \end{description} %%%%%%%% \section{Implementation details} %%%%%%%% %%%%%%%% \section{Inverse Laplacian} %%%%%%%% [As discussed above] %%%%%%%% \section{Inverse curl\,curl} %%%%%%%% The equations of inertial MHD (i.e. ``MHD'' including electron inertia) are \begin{eqnarray} \dfrac{D\ln\varrho}{Dt} &=& - \nabla\cdot\uv \; , \\[1.5ex] % \dfrac{D\uv}{Dt} &=& - \dfrac{m_{\rm e}} {\mtildei} (\vve\cdot\nabla)\vve - \dfrac{\nabla P_{\rm e}} {\varrho} + \dfrac{\jv \times \Bv} {\varrho} \; , \\[1.5ex] % \left[ 1 +\lambda_{\rm e}^2 \nabla\times (\nabla\times) \right] \; \dfrac{\partial \Av}{\partial t} &=& \Bigl[ \uv\times\Bv - \eta\mu_0\jv - \Ev_{\rm Hall+\mathit{P}_e} \nonumber \\ & & \quad{} - \mu_0 \lambda_{\rm e}^2\, \nabla\cdot(\uv\jv {+} \jv\uv) + \dfrac{m_{\rm e}}{e} (\vve\cdot\nabla)\vve - \nabla\psi \biggr] \; , \end{eqnarray} % where \emph{(link to hydro:)} \begin{eqnarray} n &=& n_{\rm e} \; , \\ \varrho &=& \mtildei n \; , \\ \mtildei &=& \dfrac{m_{\rm i}}{Z} \; , \\ \varrho\uv &=& n ( \mtildei \uvi + m_{\rm e} \uve ) \; , \\ \uv &=& \uvi + \dfrac{m_{\rm e}}{\mtildei} \uve \; , \\ \vve &=& \uv_{\rm e} - \uv_{\rm i} \ = \ \dfrac{\jv}{en} \; . \end{eqnarray} \emph{(magnetic field, etc:)} \begin{eqnarray} \jv &=& - e n \uve \; , \\ \Bv &=& \nabla\times\Av \; , \\ \mu_0 \jv &=& \nabla\times\Bv \; , \\ \nabla\cdot(\uv\jv {+} \jv\uv) &=& (\uv\cdot\nabla)\jv + (\jv\cdot\nabla)\uv + \jv(\nabla\cdot\uv)\; , \\ (\uv\cdot\nabla)\jv &=& \dfrac{1}{\mu_0} \uv\cdot \tensor{D} \; , \\ % \Ev_{\rm Hall} &=& \dfrac{\jv\times\Bv}{n e c} \; , \\ \Ev_{p_{\rm e}} &=& - \dfrac{\nabla P_{\rm e}}{n e c} \; , \\ \Ev_{\rm inertia} &=& \mu_0 \lambda_{\rm e}^2 \left[ \partial_t \jv + \nabla\cdot(\uv \jv + \jv \uv) \right] \; , \\ % (\tensor{D})_{ik} &=& \partial_i\partial_k\, \partial_l A_l + \partial_i \partial_l^2 A_k \; . \\ \end{eqnarray} %%%%%%%%%%% \subsection{The elliptic operator} At every time step, we need to invert the elliptic operator \begin{equation} \mathcal{L} \equiv 1 +\lambda_{\rm e}^2 \nabla\times (\nabla\times) = 1 - \lambda_{\rm e}^2 \widetilde{\Laplace} \; , \end{equation} where \begin{equation} \widetilde{\Laplace} \equiv -\nabla\times (\nabla\times) = \Laplace - \nabla (\nabla \cdot) \end{equation} is another elliptic operator, closely related to the Laplacian. Note that we cannot replace $\widetilde{\Laplace}$ with $\Laplace$ under the pretext of gauge invariance, because we cannot add an arbitrary gradient to $\Av$ any more -- we can only add it to the right-hand side of Eq.~(\ref{Eq-Alfven-Induction-Av})\footnote{Of course we still have the freedom to add a time-independent gradient to $\Av$, but with $\lambda_{\rm e}$ being a function of time, this does not help.} In component notation, we have \begin{equation} \widetilde{\Laplace} \fv \equiv \widetilde{\Laplace}_{ik} f_k \end{equation} with \begin{equation} \widetilde{\Laplace}_{ik} = \delta_{ik} \partial_j^2 - \partial_i \partial_k \end{equation} In matrix notation, \begin{equation} \widetilde{\Laplace} = \begin{pmatrix} \partial_y^2{+}\partial_z^2 & -\partial_x \partial_y & -\partial_x \partial_z \\ -\partial_x \partial_y & \partial_x^2{+}\partial_z^2 & -\partial_y \partial_z \\ -\partial_x \partial_z & -\partial_y \partial_z & \partial_x^2{+}\partial_y^2 \end{pmatrix} \end{equation} %%% \subsection{Discretization} To invert $\mathcal{L}$, we need to find a vector function $\fv$ that satisfies \begin{equation} (1 - \lambda_{\rm e}^2 \widetilde{\Laplace})^{-1} \gv = \fv \; . \end{equation} Thus $(1- \lambda_{\rm e}^2 \widetilde{\Laplace}) \fv = \gv$, or \begin{equation} \label{Eq-Alfven-generalized-Laplace} (\widetilde{\Laplace} - h) \fv = \tilde{\gv} \; , \end{equation} with \begin{equation} \tilde{\gv} = -\dfrac{1}{\lambda_{\rm e}^2} \gv \; \qquad h = \dfrac{1}{\lambda_{\rm e}^2} \; . \end{equation} Discretizing the first component equation of Eq.~(\ref{Eq-Alfven-generalized-Laplace}), we get \begin{eqnarray} \dfrac{\sum\limits_{\pm} f^{x}_{i,j\pm1,k} - 2 f^{x}_{i,j,k}} {\delta y^2} + \dfrac{\sum\limits_{\pm} f^{x}_{i,j,k\pm1} - 2 f^{x}_{i,j,k}} {\delta z^2} \hspace*{7em} \nonumber \\[1ex] {} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{y}_{i\pm1,j\pm'1,k}} {4\,\delta x\,\delta y} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{z}_{i\pm1,j,k\pm'1}} {4\,\delta x\,\delta z} - h_{i,j,k} f^{x}_{i,j,k} &=& \tilde{g}^{x}_{i,j,k} \; . \end{eqnarray} Isolating $f^{x}_{i,j,k}$, we get\footnote{ Compare this to the iteration \begin{equation} f_{ijk} = \dfrac{\dfrac{\sum\limits_{\pm} f_{i\pm1,j,k}} {\delta x^2} + \dfrac{\sum\limits_{\pm} f_{i,j\pm1,k}} {\delta y^2} + \dfrac{\sum\limits_{\pm} f_{i,j,k\pm1}} {\delta z^2} - g_{ijk} } {\dfrac{2}{\delta x^2} + \dfrac{2}{\delta y^2} + \dfrac{2}{\delta z^2} + h_{ijk} } \end{equation} obtained for the scalar Laplace (or rather inhomogeneous Helmholtz) equation \begin{equation} (\Laplace - h) f = g \; . \end{equation} } \begin{eqnarray} f^{x}_{i,j,k} &=& \dfrac{\dfrac{\sum\limits_{\pm} f^{x}_{i,j\pm1,k}} {\delta y^2} + \dfrac{\sum\limits_{\pm} f^{x}_{i,j,k+\sigma}} {\delta z^2} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{y}_{i\pm1,j\pm'1,k}} {4\,\delta x\,\delta y} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{z}_{i\pm1,j,k\pm'1}} {4\,\delta x\,\delta z} - \tilde{g}^{x}_{i,j,k} } {\dfrac{2}{\delta y^2} + \dfrac{2}{\delta z^2} + h_{i,j,k}} \; , \end{eqnarray} and similarly \begin{eqnarray} f^{y}_{i,j,k} &=& \dfrac{\dfrac{\sum\limits_{\pm} f^{y}_{i\pm1,j,k}} {\delta x^2} + \dfrac{\sum\limits_{\pm} f^{y}_{i,j,k\pm1}} {\delta z^2} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{x}_{i\pm1,j\pm'1,k}} {4\,\delta x\,\delta y} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{z}_{i,j\pm1,k\pm'1}} {4\,\delta y\,\delta z} - \tilde{g}^{y}_{i,j,k} } {\dfrac{2}{\delta x^2} + \dfrac{2}{\delta z^2} + h_{i,j,k}} \; , \end{eqnarray} and \begin{eqnarray} f^{z}_{i,j,k} &=& \dfrac{\dfrac{\sum\limits_{\pm} f^{z}_{i\pm1,j,k}} {\delta x^2} + \dfrac{\sum\limits_{\pm} f^{z}_{i,j\pm1,k}} {\delta y^2} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{x}_{i\pm1,j\pm'1,k}} {4\,\delta x\,\delta z} - \dfrac{\sum\limits_{\pm,\pm'} (\pm) (\pm') f^{y}_{i\pm1,j,k\pm'1}} {4\,\delta y\,\delta z} - \tilde{g}^{z}_{i,j,k} } {\dfrac{2}{\delta x^2} + \dfrac{2}{\delta y^2} + h_{i,j,k}} \; . \end{eqnarray} This equations can be used to iteratively invert $\widetilde{\Delta}$, which naturally lends itself as the soothing process of a multigrid method. %%%%%%%% \subsection{Memory management} %%%%%%%% \section{To Do} %%%%%%%% %%%%%%%% \subsection{What needs to be done} \begin{itemize} \item Parallelization \end{itemize} %%%%%%%% \subsection{Potential extensions} \begin{itemize} \item line/plane relaxation \item semi-coarsening \item FDA? \item Full multigrid cycle, W cycles, \ldots \end{itemize} %%%%%%%%% \appendix %%%%%%%% \section{Literature} %%%%%%%% \printindex %%%%%%%%% \section*{An excerpt from my lecture notes} See \url{http://www.capca.ucalgary.ca/~wdobler/teaching/phys535/notes2.pdf} \includepdf[pages=-,scale=0.9,frame,pagecommand={}]{lecture_notes.pdf} \end{document} % End of file multigrid.tex %%% Please leave this for Emacs [wd]: %% Local Variables: %% ispell-check-comments: t %% Local IspellDict: canadian %% End: % LocalWords: ifpdf pdftex dvips multigrid PDE Seidel ijk endfor Ny wavenumber % LocalWords: wavenumbers FMG MHD Dt ik Eq pagecommand pdf