diff options
Diffstat (limited to 'Carpet/doc/documentation.tex')
-rw-r--r-- | Carpet/doc/documentation.tex | 892 |
1 files changed, 892 insertions, 0 deletions
diff --git a/Carpet/doc/documentation.tex b/Carpet/doc/documentation.tex new file mode 100644 index 000000000..546b93d4d --- /dev/null +++ b/Carpet/doc/documentation.tex @@ -0,0 +1,892 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header: /home/eschnett/C/carpet/Carpet/Carpet/doc/documentation.tex,v 1.10 2003/07/21 18:48:08 schnetter Exp $ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevent thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be sparated with a \\ or a comma +% - You can define your own macros are OK, but they must appear after +% the START CACTUS THORNGUIDE line, and do not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphix package. +% More specifically, with the "includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us (later) to create a PDF version of the ThornGuide +% via pdflatex. | +% - References should be included with the latex "bibitem" command. +% - use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header: /home/eschnett/C/carpet/Carpet/Carpet/doc/documentation.tex,v 1.10 2003/07/21 18:48:08 schnetter Exp $ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../doc/latex/cactus} +\usepackage{hyperref} + +\begin{document} + +% The author of the documentation +\author{Erik Schnetter \textless schnetter@uni-tuebingen.de\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{Carpet} + +% the date your document was last changed, if your document is in CVS, +% please use: +\date{$ $Date: 2003/07/21 18:48:08 $ $} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} +This text describes the Carpet arrangement. Carpet is a mesh +refinement driver for Cactus that can replace PUGH, the standard +unigrid driver. Carpet supports multiple refinement levels and +multiple grid patches. Carpet can run in parallel, but not yet very +efficiently so. Carpet does not yet support multiple grid +hierarchies, i.e.\ shadow hierarchies or automatic convergence tests. +\end{abstract} + + + +\section{Overview} + +\subsection{Fixed Mesh Refinement, aka Box-in-Box} + +Fixed Mesh Refinement (FMR), also known as box-in-box, is a way to +increase the local resolution in unigrid applications, while retaining +the basic unigrid character of an application. A small number (maybe +two or three) of grids with varying resolution overlay each other, +where the coarsest grid has the largest extent. This allows the +application to benefit from the higher resolution of the smaller grids +while keeping the outer boundary far out at the same time. The main +advantage of FMR are that it needs far less resources than globally +increasing the resolution. + +\subsection{Carpet} + +Carpet is the name of an FMR driver, i.e.\ the back end that handles +storage allocation for the grid functions, parallelism, I/O, and the +various inter-grid operations. Carpet was developed in early summer +of 2000 by Erik Schnetter \cite{Carpet__erik-schnetter}, then a +research scholar in the Department for Astronomy and Astrophysics +\cite{Carpet__astro-psu-edu} of Penn State University +\cite{Carpet__psu-edu}. In spring 2001, Carpet was coupled to Cactus +as a drop-in enhancement for the standard unigrid Cactus driver PUGH. + +\subsection{Cactus} + +From the main Cactus web pages \cite{Carpet__cactuscode-org}: +\begin{quote} +Cactus is an open source problem solving environment designed for +scientests and engineers. Its modular structure easily enables +parallel computation across different architectures and collaborative +code development between different groups. Cactus originated in the +academic research community, where it was developed and used over many +years by a large international collaboration of physicists and +computational scientists. +\end{quote} + + + +\section{Introduction} + +\subsection{Fixed Mesh Refinement} + +A standard way of solving partial differential equations are finite +differences on a regular grid. This is also called \emph{unigrid}. +Such an application discretises its problem space onto a single, +rectangular grid which has everywhere the same grid spacing. This +grid might be broken up into several parts for parallelisation +purposes, but parallelisation should be transparent to the physics +part of the application. + +Increasing the resolution in a unigrid application is somewhat +expensive. For example, increasing the resolution by a factor of two +requires a factor of eight more storage in three dimensions. Given a +constant Courant factor, the calculation time will even go up by a +factor of sixteen. This behaviour makes it easy to find problems that +cannot be solved on contemporary supercomputers, no matter how big and +fast those computers are. + +Apart from physical insight, which often has to be used to decrease +the problem size until it fits the current hardware, there are also +numerical and algorithmic methods to decrease the resource +requirements of the application. Most applications need the high +resolution only in a part of the simulation domain. Discretisation +methods that don't require a uniform resolution, such as finite +elements, can implement non-uniform resolutions very naturally. One +problem with finite elements is that many physicists today are not +familiar with finite elements, or shy away from their perceived +complexity, or are not willing to adapt existing finite difference +code. + +Fixed Mesh Refinement (FMR) is a poor man's way of implementing a +non-uniform resolution into a unigrid application with minimal changes +to its structure. Instead of only one grid, there are several grids +or grid patches with different resolutions. The coarsest grid usually +encloses the whole simulation domain. Successively finer grids +overlay the coarse grid at those locations where a higher resolutions +is needed. The coarser grids provide boundary conditions to the finer +grid through interpolation. + +Instead of updating only one grid, the application has to update all +grids. The usual approach is to first take a step on the coarsest +grid, and then recursively take several smaller steps on the finer +grids. The Courant criterion requires that the step sizes on the +finer grids be smaller than on the coarse grid. The boundary values +for the finer grids are found through interpolation in space and time +from the coarser grid. In the end, the information on the finer grids +is injected into the coarse grids. + +Strictly speaking there is no need for a coarse grid on the regions +covered by the finer grids. But as stated above, the resources +required for treating the overlapping region on the coarse grid are +only minimal compared to treating the finer grids. And because a +coarse grid with a hole often creates complications, this obvious +optimisation is often left out. + +\subsection{Carpet} + +Carpet is a C++ library that provides infrastructure to describe +regions of varying resolution in a convenient and efficient way. +Carpet contains routines to manage grid hierarchies, containing the +relationships between the components of the grid on the different +refinement and convergence levels. Carpet has a notion of simulation +time and grid spacing, which are necessary for interpolation, and +contains efficient interpolators. + +Carpet can run on several processors in parallel using MPI for +communication. Each grid can be broken down into several components, +and every component has a home processor. Carpet also contains +operators to move certain regions to a different processor, or to +synchronise all components of a grid. + +Carpet is also an arrangement of thorns for Cactus, implementing a +driver and associated I/O routines for both ASCII and binary I/O. It +should be possible to substitute Carpet for the standard Cactus driver +PUGH without changes to the application thorns and thus use Carpet as +a unigrid driver. Making use of the FMR capabilities of Carpet +usually requires some rearranging of the application, comparable in +general to the changes necessary for a uniprocessor application to run +on multiple processors. + +The driver section of Carpet contains the logic to manage storage for +the grid functions, to traverse the grid hierarchy for all scheduled +routines, and to automatically apply the necessary inter-grid +operators for prolongation (interpolation of the fine grid boundaries) +and restriction (injecting the fine grid information back into the +coarse grid). + +The ASCII I/O routines use the quasi-standard gnuplot +\cite{Carpet__gnuplot-info} format. The binary I/O routines use the +FlexIO library \cite{Carpet__FlexIO} written by John Shalf. It allows +efficient and platform independent I/O. The FlexIO format is based on +HDF \cite{Carpet__HDF} and also supported by several visualisation +packages. + +Carpet is copyrighted by Erik Schnetter, and is available under the +GPL licence from a CVS \cite{Carpet__CVS} repository. + +\subsection{WaveToy} + +Cactus comes with a sample application called \emph{WaveToy}, which +solves the scalar wave equation with various initial data and boundary +conditions. An an example, I have extended WaveToy so that is uses +Carpet's FMR capabilities. WaveToy serves both as a test case for +Carpet, and as example of how to convert an application to using FMR. + +The equation solved by WaveToy is the well known scalar wave equation, +discretised using the Leapfrog method with three time levels, yielding +second order accuracy in space and time. A typical set of initial +data are a plane wave, and a typical boundary condition is +periodicity. Those allow long term simulations as well as easy and +meaningful comparisons to the analytic solution. + + + +\section{Compiling Cactus With Carpet} + +Carpet has been written in C++, using templates and the STL (Standard +Template Library). Both templates and the STL make writing and +debugging code a lot easier. Without templates, I would have had to +put much effort into making Carpet support all of Cactus' data types. +Without the STL, I would have had to spend quite some time +implementing basic containers such as lists or sets. I still had to +implement a custom vector type, because STL's vector type is optimised +for large vectors only, and I needed threedimensional vectors of +integers. + +The inner loops of Carpet are the inter-grid operators, that is the +routines that copy, restrict, and prolongate between grids. Due to +Cactus it was rather easy to write these in \textsc{Fortran 77}, which +makes them both fast and portable. + +Carpet is an arrangement in Cactus. It can theoretically be compiled +without any external library, if you omit the binary I/O support which +requires the FlexIO library. FlexIO is already part of Cactus in the +thorn CactusExternal/FlexIO. I suggest that you enable support for +the HDF format in the FlexIO library, although this is not necessary. +For that, you have to install the HDF5 libraries first. + +\subsection{Hurdle 1: STL} + +Some operating systems do not have a compliant STL (Standard Template +Library) installed. If not, then you are in trouble. Carpet does +make use of the STL, and there is no way around that. + +\subsection{Hurdle 2: Templates} + +Some compilers contain switches to instantiate some or all templates +automatically. This usually does not work when files are put into +libraries, which is what Cactus does. The scheme that I found working +on all machines is to instantiate most templates by hand, and have the +compiler instantiate the missing templates for every object file. +This is the default for gcc. On SGIs, you have to pass the options +\texttt{-no\_auto\_include -ptused} to the C++ compiler. + +The C++ standard specifies a limit when using templates as template +parameters. Carpet's use of the GNU STL exceeds this limit. Gcc +requires the option \texttt{-ftemplate-depth-30} to enable this. + +\subsection{WaveToy} + +Unfortunately, PUGH and Carpet cannot yet be both compiled into a +single application. (This will be fixed soon.) That means that you +will have separate executables for unigrid and for mesh refinement +applications. + +Configuring Carpet is not quite trivial, because Cactus provides +currently no way to autodetect the settings for Carpet. Hence you +will have to set the settings manually. I propose that you start with +on of the pre-made options files in the directory +\texttt{Carpet/Carpet/options}. Try e.g.\ \texttt{carpet-harpo-sgi} +for an SGI, or \texttt{carpet-lilypond} for Linux with gcc, or +\texttt{carpet-lilypond-ic} for Linux with the Intel compilers. Once +you have a working options file for your machine, send it to me, so +that I can include it. + +As for the thorn list: Carpet has its own ASCII output thorn, which +outputs more information than CactusBase/IOASCII. The thorn list that +I use is + +\begin{verbatim} +CactusBase/Boundary # boundary (grid) [ ] { } +CactusBase/CartGrid3D # grid ( ) [ ] {driver} +#CactusBase/IOASCII # IOASCII (IO,Hyperslab) [ ] {IO} +CactusBase/IOBasic # IOBasic (IO) [ ] {IO} +CactusBase/IOUtil # IO ( ) [ ] { } +CactusBase/LocalInterp # LocalInterp ( ) [ ] { } +CactusBase/Time # time ( ) [ ] { } +CactusConnect/HTTPD # HTTPD (Socket) [ ] {Cactus} +CactusConnect/HTTPDExtra # http_utils (httpd,IO) [ ] { } +CactusConnect/Socket # Socket ( ) [ ] { } +CactusExternal/FlexIO # FlexIO ( ) [ ] { } +CactusExternal/jpeg6b # jpeg6b ( ) [ ] { } +CactusIO/IOJpeg # IOJpeg (IO,Hyperslab,jpeg6b) [ ] {IO} +CactusUtils/NaNChecker # NaNChecker ( ) [ ] { } +CactusWave/IDScalarWave # idscalarwave (wavetoy,grid) [ ] {grid} +CactusWave/IDScalarWaveC # idscalarwave (wavetoy,grid) [ ] {grid} +CactusWave/IDScalarWaveCXX # idscalarwave (wavetoy,grid) [ ] {grid} +#CactusWave/IDScalarWaveElliptic # idscalarwaveelliptic (grid,wavetoy,ellbase) [ ] {idscalarwave} +CactusWave/WaveBinarySource # binarysource (wavetoy,grid,idscalarwave) [ ] { } +CactusWave/WaveToyC # wavetoy (Grid,Boundary) [ ] { } +CactusWave/WaveToyCXX # wavetoy (Grid,Boundary) [ ] { } +CactusWave/WaveToyF77 # wavetoy (Grid,Boundary) [ ] { } +#CactusWave/WaveToyF90 # wavetoy (Grid,Boundary) [ ] { } +#CactusWave/WaveToyFreeF90 # wavetoy (Grid,Boundary) [ ] { } +Carpet/Carpet # driver (CarpetLib) [ ] {Cactus,IO} +Carpet/CarpetIOASCII # IOASCII (CarpetLib,driver,Hyperslab) [ ] {IO} +Carpet/CarpetIOFlexIO # IOFlexIO (CarpetLib,driver,Hyperslab,FlexIO) [ ] {IO} +#Carpet/CarpetIOHDF5 # IOHDF5 (CarpetLib,driver,Hyperslab) [ ] {IO} +#Carpet/CarpetIOSer # IOSer (CarpetLib,driver,Hyperslab) [ ] {IO} +Carpet/CarpetLib # CarpetLib ( ) [ ] { } +Carpet/CarpetReduce # reduce (CarpetLib,driver) [ ] { } +Carpet/CarpetRegrid # CarpetRegrid (CarpetLib,driver) [ ] { } +Carpet/CarpetSlab # Hyperslab (CarpetLib,driver) [ ] { } +\end{verbatim} + +The thorns prefixed with \texttt{\#} are disabled. IOASCII conflicts +with CarpetIOASCII. I disabled IDScalarWaveElliptic because there is +no elliptic solver for mesh refinement, and I disabled WaveToyF90 and +WaveToyFreeF90 because gcc does not yet contain a Fortran 90 compiler. +CarpetIOHDF5 is not yet finished, and CarpetIOSer needs the Ser +library which is not publically available. + +The CactusConnect, CactusIO, and CactusUtils thorns are not necessary, +but are nice to have around. You can safely omit these. + + + +\section{Running The Example Applications} + +Although Carpet works fine with the standard WaveToy thorns, all the +example parameter files in the CactusWave arrangement use PUGH, and +can therefore not be directly used. + +The coordinate thorn CactusBase/CartGrid3D does not provide periodic +boundary conditions. These are normally provided by the driver PUGH. +However, Carpet does not contain any boundary conditions. If you want +to apply periodic boundaries, you will therefore have to use the +AlphaThorns/Cart3d coordinate thorn instead, which does provide +periodicity. Unfortunately, AlphaThorns/Cart3d is incompatible with +CactusBase/CartGrid3D. There is a version of WaveToy in the Carpet +arrangement that has been adapted to AlphaThorns/Cart3d. I suggest +that you use this version of WaveToy instead of CactusWave to run test +problems, because periodicity makes for nice testing setups. + +You can find quite a few example parameter files in the directory +\texttt{Carpet/WaveToyF77/par}. I especially recommend the +\texttt{wavetoyf77\_periodic\_*} set, which comes in two sizes +(\texttt{coarse} and \texttt{fine}, corresponding to a small and a +large simulation domain) and three different refinement hierarchies +(with one, two, and three level altogether, respectively). This set +thus forms a convergence test, which you can run and test yourself. +The set \texttt{wavetoyf77\_rad\_full\_*} uses radiative instead of +periodic boundaries and should also be nice to look at. The file +\texttt{wavetoyf77\_rad\_automatic.par} is an attempt at adaptive mesh +refinement, which may or may not work, depending on the current status +of Carpet. + +Second order convergence requires second order interpolation in time, +which requires that at least three time levels are present. + + + +\section{Fold Your Own FMR Application} + +There are three steps to take from a simple unigrid uniprocessor toy +application to a full-blown FMR multiprocessor production application. +Those steps are almost independent, and I would like to explain them +and their implications in some detail below. + +\subsection{Multiple Processors} + +The probably best known of these is the step from using one to using +several processors, also known as parallelisation. Because many +people are already familiar with this step, I will describe it first. + +In a uniprocessor application, it is possible to access every grid +point in arbitrary manners. In order to allow multiple processors to +run efficiently in parallel, the grid is broken down into several +rectangular components, and each processor is assigned one of these +components. + +The components will usually overlap by a few grid points, so as to +allow the processors to e.g.\ calculate spatial derivatives (which +require neighbouring grid points) without having to communicate for +every grid point. From time to time it is then necessary to +synchronise the overlapping region, which is the only time at which +communication happens. This allows the application to run almost +unchanged, i.e.\ without invoking communication itself. The +synchronisation routine is provided by the driver and not by the +application. + +Of course a serial applicate usually will have to be changed to +support multiple processors. In order to do so, all the operations +that the application performs have to be classified into one of two +categories: + +One category contains the so-called \emph{local} operations. These +are operations that are applied to each and every grid point +individually, and that do not depend on any other grid point except +nearby neighbours. Each local operation will thus involve a loop over +all grid points, and in order to run on multiple processors, after +each such loop the synchronisation routine has to be called. An +example of a local operation would be calculating a spatial +derivative. + +The other category contains so-called \emph{global} operations. These +operations do not depend on individual grid points, and thus do not +involve loops over grid points. The result of a global operation is +the same on all processors; therefore global operations don't involve +communication and don't require synchronisation. An example of a +global operation would be to check how many time steps have been +taken, and decide whether the simulation should be terminated. + +Typically most operations can be classified or rewritten to be either +local or global. But often there are operations that fit neither +category, and these parts of an application are hardest to +parallelise. Applying the boundary conditions, to give another +example, might seem at first to be neither local nor global. But in a +slight (yet completely correct) stretch of the term "applied to all +grid points", boundary conditions can be classified as local; they are +a local operation that just does nothing to most grid points. + +To give one more example, calculating an error norm does not fit these +categories. It is neither local nor global. It is not local because +the results involved all grid points (and not only nearby neighbours), +and it is not global because it does involve the grid points. All +operations that do not fit the two category require typically special +handling, and often require hand-coded communication in the +application. Luckily calculating various norms is such a common case +that there are special routines for that already present, called +\emph{reduction operators}. + +\subsection{Multiple Resolution Levels} + +There are several reasons why an application might want to incorporate +more than one grid, overlapping and each with a different resolution. + +The most commonly known reason is probably a convergence test, where +the very same problem is treated in different resolutions. +Differences in the result are then likely caused by insufficient +resolution on the coarser (or on all) grids. For a convergence test, +the grids are completely independent, and it does not matter whether +the simulation runs on all grids simultaneously or sequentially. In +order to treat the grid sequentially, the application does not have to +be changed at all. + +The reason of interest here is of course FMR. For FMR, the order in +which the grids are treated is fixed. As described above, there is +first a time step on the coarse grid, and then recursively several +smaller steps on the finer grids. This order does require certain +changes in the application. The sequence of operations that form a +single time step have to be identified and isolated. (Which is to say +that there has to be a routine that calculates a time step, that is, a +complete time step, and nothing else.) It is then the task of the FMR +driver to call this routine for the correct grids in the correct +order. + +Other reasons for multiple resolution levels are e.g.\ multigrid +algorithms for elliptic equations, which I do not want to mention +here, or shadow hierarchies to determine truncation errors, which I +also want to skip here. Shadow hierarchies are very similar to the +convergence levels described above. + +Apart from this order in which the operations are performed on the +grids, there is one more complication for FMR. The boundary values of +the finer grids have to be calculated from the coarser grids through +interpolation. An because the time steps on the finer grids are +smaller, there is not always a corresponding value on the coarser +grids available. This makes it necessary to interpolate in time +between time steps on the coarser grids. The alternative would be to +take smaller steps on the coarser grids, and this would be very +expensive. + +These interpolations in time make it necessary that the driver knows +which grid function contains values corresponding to what time. The +usual way to achieve this is to have several time levels per grid +function; three time levels allow for a second order interpolation in +time. Only grid functions with enough time levels can be +interpolated, i.e.\ boundary conditions can be calculated only for +those. + +Fortunately time levels are rather widespread in applications, so they +are no new concept to introduce. Unfortunately they are often abused, +so that values corresponding to the wrong time are stored in a time +level, usually with the excuse of saving storage. This will in +general not work with FMR, because the driver then cannot interpolate +in time, leading to incorrect values on the boundaries of the finer +grids. + +\subsection{Multiple Grid Components} + +Sometimes it is convenient to have a simulation domain that is not a +rectangle. It might instead be an L-shaped simulation domain, or a +domain that consists of two disconnected rectangular regions. This +issue becomes more important with FMR, because there it is often +convenient to have several disconnected refined regions. As long as +there are enough processors available, each processor can be assigned +a region or a part thereof, and no new concept need be introduced. +If, however, there are fewer processors than regions, then a new +problem arises. + +A common case for that problem might be a simulation containing just +two refined regions, and running on a single processor. The refined +grid the consists of two component. The problem then is that the two +components cannot be treated sequentially: Imagine the time evolution +routine working on (say) the first component. It will at some time +call the synchronisation routine. At that time there are no values +from the second component available, because the second component has +not been treated yet. Therefore the synchronisation routine cannot +complete. That means in turn that the time evolution routine cannot +complete working on the first component, leading to a deadlock. Work +on neither component can be completed before work on the other +component. + +The solution is to break up the time evolution routine into several +smaller routines, each consisting of a single either local or global +operation. (``Local'' and ``global'' have here the exact same +meanings that were defined above for parallelisation.) A local +operation works, by definition, on individual grid points. Hence the +local routines have to be called once for every grid component. A +global operation, by definition, does not depend on individual grid +points. Hence it has to be called only once per processor, and not +once per component. That means that the driver has to be told the +category individual routine is in. + +\subsection{Example} + +Let me finish this section with an detailed example. Suppose you want +to solve the equation +\begin{eqnarray} + \frac{d}{dt} u & = & f(u) \quad , +\end{eqnarray} +integrating using the midpoint rule, i.e.\ the simplemost second-order +time integration scheme. Given values at the previous time $u^{n-1}$, +one first calculates a first order solution using an Euler step, +leading to the intermediate result +\begin{eqnarray} + v^n & = & u^{n-1} + dt\; f(u^{n-1}) \quad . +\end{eqnarray} +The second and final step is then calculated via +\begin{eqnarray} + u^n & = & u^{n-1} + dt\; f(\frac{1}{2} [u^{n-1} + v^n]) \quad . +\end{eqnarray} + +The corresponding pseudo code would look like +\begin{enumerate} +\item +Calculate Euler step, storing the result into $u^n$ +\item +Apply boundary conditions to $u^n$ +\item +Synchronise $u^n$ +\item +Calculate average of $u^{n-1}$ and $u^n$, storing the result into +$v^n$ +\item +Calculate second step, storing the result again into $u^n$ +\item +Apply boundary conditions again to $u^n$ +\item +Synchronise again $u^n$ +\end{enumerate} + +The above algorithm looks a bit different from a naive implementation +of the midpoint rule. One difference is that both the first and the +second step store their result into $u^n$. This is necessary because +it would be inconvenient to apply boundary conditions to the +intermediate value $v^n$. Remember, in order to apply boundary +conditions on the finer grids, there have to be several time levels +present. With the above scheme, only $u$ needs several time levels. +$v$ is used only as a temporary (and could conceivably be completely +eliminated). + +Note also that the first step goes all the way from time level $n-1$ +to time level $n$. The midpoint rule can be rewritten (in fact, is +usually written) so that the first step is only a half step, leading +to the time level $n - \frac{1}{2}$. This is not possible for FMR, +because interpolating to the time $n - \frac{1}{2}$ is not possible, +and thus there could be no boundary conditions applied after the first +step. + +The second thing to note is that the application of the boundary +condition and the synchronisation have been separated rather +artificially. Normally synchronisation would be considered part of +the boundary condition. In this case, however, the applying the +boundary condition is a local operation, whereas synchronisation +counts as global operation. (It is not obvious that synchronisation +should be global, but as the synchronisation routine is a part of +Carpet, it was up to me to decide this.) As explained above, local +and global operations have to be separated. + +Separating the evolution steps and the boundary condition routines is, +on the other hand, just a notational convenience. There could well be +a single routine implementing both. + +For Cactus, the order in which to call the individual parts of the +time evolution routines is described in the schedule routines, i.e.\ +in the files called \texttt{schedule.ccl}. By default a routine is +assumed to be local; global routines have to be tagged with +\texttt{OPTIONS: GLOBAL}. + +The tag \texttt{SYNC: groupname} indicates that the group +\texttt{groupname} should be synchronised after the scheduled routine +has been called for all grid components. This obviously makes sense +only for local routines. Using the \texttt{SYNC:} tag is preferred +over calling the synchronisation routine \texttt{CCTK\_SyncGroup} +directly. + +The example thorn WaveToy in Carpet's arrangement is a bit simpler +than what is described here, because it uses the Leapfrog scheme which +consists of only a single step. I would suggest looking at WaveToy as +an initial FMR example. + +The thorn SpaceToy is implemented very close to the way described +here. It evolves two variables phi and psi, but it is also coupled to +the thorn HydroToy. This coupling introduces some additional +complications. The thorn HydroToy, on the other hand uses a +predictor-corrector scheme, which is also a two step scheme and thus +more complex that WaveToy. All the coupling between SpaceToy and +HydroToy is contained in SpaceToy. I would thus suggest looking at +HydroToy first. + +I assume that converting an application to FMR is straightforward +after handling the time levels has been straightened out. + + + +\section{Further documentation} + +The individual thorns in the Carpet arrangement might contain further +documentation, which is also available in the thorn guide. +Additionally, there is a document \texttt{internals.tex} in the +arrangement's doc directory, and a document +\texttt{threelev\_initdata.tex} in thorn \texttt{Carpet}'s doc +directory. + + +\section{Frequently Asked Questions} +\label{sec:faq} + +Here are a few of the more frequently asked questions with some +answers. +\begin{enumerate} +\item \textbf{If I run without any refined grids, why don't I get the + same results as with PUGH?} + + There are two possible reasons. The most common is that the you are + not comparing exactly the same output. It used to be the case that + norms would disagree (this is no longer the case). If it is the + ASCII output that disagress, then you should note that the default + output format for CarpetIOASCII gives more digits than + CactusBase/IOASCII. If you want to get ``identical'' results for + this output, try setting \texttt{IOASCII::out\_format = ".14f"}). + + The second reason is subtle differences are bugs in the + implementation. Good luck finding these... +\item \textbf{I switch on a refined grid. Why do I not see it output? + Why is the output strange?} + +\begin{figure}[htbp] + \begin{center} + \includegraphics[scale=0.5]{Grid1.eps} + \caption{How the grids are indexed in Carpet. This is an + artificial three level example using C-style numbering (0 + origin). Note that the numbering is with respect to the finest + grid.} + \label{fig:Grid1} + \end{center} +\end{figure} + As soon as you switch on refinement the way the grids are numbered + by index changes. The numbering is done with respect to the + \textit{finest} grid but covers the entire domain. An example of how + the numbering works is given in figure~\ref{fig:Grid1}. It is + important to note that this also applies to the numbering in + time. So with the grid structure of figure~\ref{fig:Grid1} output + for the coarsest grid only occurs on iterations $0,4,8,\dots$, for + the medium grid only on iterations $0,2,4,\dots$, and for the finest + grid on iterations $0,1,2,\dots$. Note that here the finest grid is + not the finest \textit{existing} grid, but the finest + \textit{possible} grid. This is controlled by the + \texttt{Carpet::max\_refinement\_levels} parameter. + + So, there are plenty of reasons why the output might be strange: + \begin{itemize} + \item You are requesting output on iterations when not all grids are + output. For example, requesting output every $5^{\text{}th}$ + iteration with the above grid structure would only output the + coarse grid every 20 iterations. + \item You are requesting output along an index that does not + intersect with any grid points. For example, the line defined by + $j = 6$ in the example above corresponds to the center of the box, + but does not intersect the coarse grid at all! + \item Requesting output along a line defined by a coordinate value + will give you the index closest to it. This may not agree on the + different refinement levels. In the example above the coordinate + value $y=5.1$ is closest to $j=5$ on the fine grid, $j=6$ on the + medium grid, and $j=4$ on the coarse grid. All the different lines + will be output but you should not expect points that appear to + overlap in the output to agree as they're actually not at the same + point. + \item CarpetRegrid (which sets up the refined boxes) knows nothing + about symmetries. So if you have a simulation in, for example, + octant mode with $x,y,z\in[0,10]$ and you leave all the parameters + to be the defaults, the following will happen: + \begin{itemize} + \item CarpetRegrid creates a refined box at the center of the + \textit{index space}. This might cover something like + $x,y,z\in[3,7]$. + \item When the IO thorn requests the output lines and planes it + does know the symmetries, so tries to put the lines and planes + as close to the origin $x=y=z=0$ as possible. + \item When output occurs the lines and planes don't intersect the + fine grid and so you get no output. + \end{itemize} + \end{itemize} + + Morals: Comparing 1D output on different refinement levels can be + very frustrating. 2D output is usually much more informative. Using + symmetry conditions with Carpet is tricky. + +\item {\bf I want to run with periodic boundaries. Why aren't things + working correctly?} + + You thought symmetry boundaries were bad? Periodic boundaries are + even worse. + + Firstly, Carpet does not itself implement periodic boundaries. The + thorn {\tt TAT/Periodic} is ``more or less'' driver independent and + does. This should be used to implement the actual boundary + conditions. You should not need to change your code - just activate + the thorn with the appropriate parameters. + + Secondly, periodic boundaries do {\bf not} work the same way as + symmetry boundaries. This is because you cannot specify a point in + coordinate space where the boundary actually lies; it really lies in + the index space. The following example will hopefully help. + + Take a 1D slice through the grid. There are 7 points with 2 boundary + (ghost) zones (0,2 and 10,12), so only 3 points are actually being + evolved (4, 6, 8). Periodic boundaries means that the boundary points + are identified with certain evolved points. For example, point 2 is + to the left of the first evolved point and so must be identified + with the \textit{last} evolved point (8). The identifications are + shown in figure~\ref{fig:Periodic1}. + \begin{figure}[htbp] + \begin{center} + \includegraphics[scale=0.5]{Periodic1.eps} + \caption{Periodic grids identify boundary points and interior + points. The interior points are given by circles and the + boundary points by squares. The identifications are shown by the + arrows.} + \label{fig:Periodic1} + \end{center} + \end{figure} + + We then want to place a refined region across the entire width of + the domain but also have the correct periodic boundaries. The + crucial point is to ensure that points that are identified on the + coarse grid are identified in the same way on the fine grid. For + example, point 2 must still be identified with point 8. Therefore + point 2 must remain a boundary point and point 8 an interior + point. Point 4 must also be identified with point 10. There are + therefore 2 possibilities: + \begin{itemize} + \item Point 3 is the first interior point on the refined grid and + point 8 the last. Therefore the point to the ``left'' of point 3, + point 2, is still identified with point 8. + \item Point 4 is the first interior point on the refined grid and + point 9 the last. This possibility is illustrated in + figure~\ref{fig:Periodic2}. + \end{itemize} + \begin{figure}[htbp] + \begin{center} + \includegraphics[scale=0.5]{Periodic2.eps} + \caption{A periodic refined grid. The boundary zones are blue + plus signs, the interior blue crosses. Note that the interior + points on the refined grid extend \textit{outside} the + interior on the base grid. However, equivalent points on both + grids coincide.} + \label{fig:Periodic2} + \end{center} + \end{figure} + + So to specify the particular refined grid shown in + figure~\ref{fig:Periodic2} you would specify a lower bound of 2, an + upper bound of 11, and that both boundaries are outer boundaries. An + example for a $44 \times 7 \times 7$ grid where the ``centre half'' + of the grid in the $x$ direction is refined and the refined region + covers the entirety of the $y$ and $z$ directions, you could use +\begin{verbatim} +carpet::max_refinement_levels = 2 +carpetregrid::refinement_levels = 2 +carpetregrid::refined_regions = "manual-gridpoint-list" +carpetregrid::gridpoints = "[ [ ([22,2,2]:[62,11,11]:[1,1,1]) ] ]" +carpetregrid::outerbounds = "[ [ [[0,0],[1,1],[1,1]] ] ]" +\end{verbatim} + +\end{enumerate} + +%% \bibliographystyle{amsalpha} % initials + year +%% \bibliography{carpet} + +\begin{thebibliography}{{Pen}} + +\bibitem[AA]{Carpet__astro-psu-edu} +{Department for} Astronomy and Astrophysics, + \emph{{http://www.astro.psu.edu/}}. + +\bibitem[{Cac}]{Carpet__cactuscode-org} +{Cactus web pages}, \emph{{http://www.cactuscode.org/}}. + +\bibitem[CVS]{Carpet__CVS} +CVS, \emph{{http://www.cvshome.org/}}. + +\bibitem[gnu]{Carpet__gnuplot-info} +gnuplot, \emph{{http://www.gnuplot.info/}}. + +\bibitem[HDF]{Carpet__HDF} +HDF, \emph{{http://hdf.ncsa.uiuc.edu/}}. + +\bibitem[{Pen}]{Carpet__psu-edu} +{Penn State University}, \emph{{http://www.psu.edu/}}. + +\bibitem[Sch]{Carpet__erik-schnetter} +Erik Schnetter, \emph{{\textless + schnetter@uni-tuebingen.de\textgreater}}. + +\bibitem[Sha]{Carpet__FlexIO} +John Shalf, \emph{{FlexIO} library: + {http://zeus.ncsa.uiuc.edu/\textasciitilde jshalf/FlexIO/}}. + +\bibitem[TAT]{Carpet__tat-physik-uni-tuebingen-de} +Theoretische Astrophysik T\"ubingen, + \emph{{http://www.tat.physik.uni-tuebingen.de/}}. + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} |