%
\documentstyle [11pt,epsf] {article} 
%
\pagestyle{empty}
%
\topmargin=-0.75in
\hoffset=-0.75in
\textheight=9.0in \vsize\textheight
\textwidth=6.5in \hsize\textwidth
%
\newcommand{\ad}{{\tt ad}}
\newcommand{\NR}{{\tt NR}}
\newcommand{\NT}{{\tt NT}}
\newcommand{\tabspace}{0.5cm}
\newcommand{\T}{\hspace*{\tabspace}}
\newcommand{\hilo}[2]{{}^{\raise2pt\hbox{$\scriptstyle{#1}$}}
                _{\lower2pt\hbox{$\scriptstyle{#2}$}}}
\newcommand{\ind}[2]{{}\hilo{#1}{#2}}
\newcommand{\dt}{\,\triangle t\,}
\newcommand{\dr}{\,\triangle r\,}
\newcommand{\lb}{\left(}
\newcommand{\rb}{\right)}
%
\title{
   A Summary of Choptuik's Implementation of Mesh Refinement for the
      Wave Equation in Spherical Symmetry
}
%
\date{October 1994}
%
\author{  }
%
\begin {document}
%
\maketitle
%
\section{Introduction}
%
In the past few years, Matthew Choptuik has reported on some exciting new
results in numerical relativity, namely, the discovery of critical phenomena
in the black hole formation via. imploding scalar fields \cite{MattPRL}.
The success
of his code depends crucially on the use of a mesh refinement algorithm due
to Berger and Oliger \cite{BergerOliger}.  Realizing that his code could
be used as a benchmark for mesh refinement applied to numerical relativity,
Choptuik presented a ``stripped down'' version of his code (\ad) in March 
of 1994.
This version takes out all of the physics complications by simply evolving
a scalar field via. the wave equation in spherical symmetry.  What is left is
the ``bare essentials'' of the mesh refinement algorithm.  This paper is 
a summary of the ideas behind Choptuik's code and the Berger and Oliger mesh
refinement algorithm.  The first two sections (Getting the Source Code and
The Model Problem) are taken word for word from Choptuik's documentation
that accompanies the code.  The next section(s) describe the mesh refinement
technique used.
%
\section{Getting the Source Code(Choptuik)}
%
Source code, documentation, sample control files and some utilities are 
available via anonymous {\tt ftp} from {\tt helmholtz.ph.utexas.edu}.
It is assumed that the code will be installed on a ``typical'' scientific 
workstation (Sun, SGI, IBM RS6000, Dec Alpha, $\, \ldots$) running 
Unix, and that {\tt C} and {\tt FORTRAN 77} compilers and a {\tt perl}
interpreter are available.

To retrieve \ad, {\tt cd} to the directory where you want \ad\ installed,
then enter:
%
\tt
\begin{verbatim}
          ftp helmholtz.ph.utexas.edu
          anonymous
          <user name>@
\end{verbatim}
\rm
%
Once the {\tt ftp} connection has been established, type
%
\tt
\begin{verbatim}
          bin
          cd pub
          get ad.tar.Z
          quit
\end{verbatim}
\rm
%
Now execute
%
\tt
\begin{verbatim}
          zcat ad.tar.Z | tar xvf -
\end{verbatim}
\rm
%
to reconstruct the \ad\ directory hierarchy.  The file {\tt ad/README} 
will tell you how to proceed.
%
\section{The Model Problem(Choptuik)}
%
\ad\ solves the wave equation in spherical symmetry using second-order 
finite difference techniques and adaptive mesh refinement.  The 
fundamental dynamical variable is a scalar field, $\phi(r,t)$,
whose evolution is governed by 
%
\begin {equation} 
\frac{\partial^2 \phi} {\partial t^2} =
	\frac{1}{r^2} \frac{\partial}{\partial r} \lb
		r^2 \frac{\partial \phi} {\partial r} \rb =
   3 \frac{\partial}{\partial r^3} \lb
		r^2 \frac{\partial \phi} {\partial r} \rb
\label {MODEL1}
\end   {equation}
%
\noindent
I rewrite this equation in first-order (in time) form, by introducing
auxiliary variables $\Phi$ and $\Pi$ defined by:
%
\begin {equation}
\Phi \equiv \frac{\partial \phi}{\partial r} \quad\quad\quad\quad\quad
\Pi  \equiv \frac{\partial \phi}{\partial t}
\label {MODEL2}
\end   {equation}
%
With these definitions, equation~(\ref{MODEL1}) may be rewritten as:
%
\begin {equation}
\frac{\partial \Phi}{\partial t} = \frac{\partial \Pi}{\partial r}
\quad\quad\quad\quad\quad
\frac{\partial \Pi}{\partial t} = 3 \frac{\partial}{\partial r^3} \lb
	r^2 \Phi \rb
\label {MODEL3}
\end   {equation}
%
The wave equation is solved on the finite domain $0 \le r \le R$, 
$0 \le t \le T$.
Introducing a standard finite difference notation for grid functions
(with a {\tt FORTRAN} convention for the spatial grid index): 
%
\begin {equation}
f\ind{n}{j} \equiv f \lb r\ind{}{j} , t \ind{n}{} \rb 
            \equiv f \lb \lb j  - 1 \rb \dr , n \dt \rb
\label {MODEL4}
\end   {equation}
%
where $\dr$ and $\dt$ are the {\sl uniform\/} spatial (radial) and 
temporal mesh-spacings
respectively, I use the following second-order finite difference approximations 
of~(\ref{MODEL3}):
%
\begin {eqnarray}
	\frac{\Phi\ind{n+1}{j} - \Phi\ind{n-1}{j}}{2 \dt} & = &
	\frac{\Pi \ind{n}{j+1} - \Pi \ind{n}{j-1}}{2 \dr}
\label{MODEL5}
\\
	\frac{\Pi \ind{n+1}{j} - \Pi \ind{n-1}{j}}{2 \dt} & = &
	3 \, \frac{r\ind{2}{j+1} \, \Phi\ind{n}{j+1} - 
           r\ind{2}{j-1} \, \Phi\ind{n}{j-1}}
          {r\ind{3}{j+1} - r\ind{3}{j-1}}
\label{MODEL6}
\end   {eqnarray}
%
Note that the spatial index, $j$, runs over the values 
$1, \, 2, \, \ldots \, , \, \NR$, where $\NR \equiv R / \dr + 1$, and 
the temporal index, $n$, runs over the values $0, \, 1 \, \ldots \, , \,
\NT$, where $\NT \equiv T / \dt + 1$.
(\ref{MODEL5}) and (\ref{MODEL6}) are trivially solved for 
$\Phi\ind{n+1}{j}$ and $\Pi\ind{n+1}{j}$, yielding the basic algebraic 
equations used to advance the solution from one time-step to the next:
%
\begin {eqnarray}
   \Phi\ind{n+1}{j} & := & \Phi\ind{n-1}{j} + 
   \frac{\dt}{\dr} \lb \Pi \ind{n}{j+1} - \Pi \ind{n}{j-1} \rb
\label{MODEL7}
\\
   \Pi \ind{n+1}{j}  & := & \Pi \ind{n-1}{j} + 
   \frac{6 \dt}{r\ind{3}{j+1} - r\ind{3}{j-1}} \, 
      \lb r\ind{2}{j+1} \, \Phi\ind{n}{j+1} -
           r\ind{2}{j-1} \, \Phi\ind{n}{j-1} \rb
\label{MODEL8}
\end   {eqnarray}
%
(see {\tt subroutine elfppi()} and {\tt subroutine elfpii()} in
{\tt unigrid.f}).  Equations (\ref{MODEL7}--\ref{MODEL8}) are used to 
update {\sl interior\/} grid function values: i.e. those with spatial 
grid indexes $2 \le j \le \NR$.  At the boundaries $j = 1$, $j = \NR$,
the following additional conditions are used 
(see {\tt subroutine elfpp()} and {\tt subroutine elfpi()} in
{\tt unigrid.f}):
%
\begin {eqnarray}
   \Phi\ind{n+1}{1} & := & 0
\label{MODEL9}
\\
   \Pi\ind{n+1}{1} & := & \frac{4 \Pi\ind{n+1}{2} - \Pi\ind{n+1}{3}}{3}
\label{MODEL10}
\\
   \Phi\ind{n+1}{\NR} & := & \lb \frac{3}{\dt} + \frac{3}{\dr} + 
                              \frac{2}{r\ind{}{\NR}} \rb^{-1}
		 \lb \frac{4 \Phi\ind{n}{\NR} - \Phi\ind{n-1}{\NR}}{\dt} +
           \frac{4 \Phi\ind{n+1}{\NR-1} - \Phi\ind{n+1}{\NR-2}}{\dr} \rb
\label{MODEL11}
\\
   \Pi\ind{n}{\NR} & := & \lb \frac{3}{\dt} + \frac{3}{\dr} + 
                              \frac{2}{r\ind{}{\NR}} \rb^{-1}
		 \lb \frac{4 \Pi\ind{n}{\NR} - \Pi\ind{n-1}{\NR}}{\dt} +
           \frac{4 \Pi\ind{n+1}{\NR-1} - \Pi\ind{n+1}{\NR-2}}{\dr} \rb
\label{MODEL12}
\end   {eqnarray}
%
Equations (\ref{MODEL9}--\ref{MODEL10}) are based on regularity of the scalar 
field $\phi$ at the origin:
%
\begin {equation}
	\lim_{r\to0} \,\, \phi \lb r , t \rb = \phi_0 \lb t \rb + 
       r^2 \, \phi_2 \lb t \rb \, + \, \cdots 
\label{MODEL13}
\end   {equation}
%
while (\ref{MODEL11}--\ref{MODEL12}) are based on the (presumed) outgoing-only
nature of solutions of~(\ref{MODEL1}) at large $r$
%
\begin {equation}
	\lim_{r\to\infty} \,\, r \phi \lb r , t \rb  = f \lb r - t \rb
\label{MODEL14}
\end   {equation}
%
which implies 
%
\begin {eqnarray}
	\lim_{r\to\infty} \,\,\, \frac{\partial \Phi}{\partial t} \,\, + \,\,
         \frac{\partial \Phi}{\partial r} \,\, + \,\, \frac{\Phi}{r} & = &  
         O \lb r^{-3} \rb \approx 0
\label{MODEL15}
\\
	\lim_{r\to\infty} \,\,\, \frac{\partial \Pi}{\partial t} \,\, + \,\,
         \frac{\partial \Pi}{\partial r} \,\, + \,\, \frac{\Pi}{r} & = & 0
\label{MODEL16}
\end   {eqnarray}
%
Equations~(\ref{MODEL11}--\ref{MODEL12}) follow 
from~(\ref{MODEL15}--\ref{MODEL16})
by replacement of first derivatives with second-order {\sl backwards}
differences applied at $(r\ind{}{\NR},t\ind{n+1}{})$.

In order to regularize (damp high frequency components) the numerical solution,
\ad\ applies a form of dissipation advocated by Kreiss \& Oliger~\cite{KO}
for use with wave-equation discretizations such as~(\ref{MODEL7}--\ref{MODEL8}).
For example, following the computation of the advanced values 
$\Phi\ind{n+1}{j}$ from~(\ref{MODEL7}), dissipation is added 
via:
%
\begin {equation}  
	\Phi\ind{n+1}{j} := \Phi\ind{n+1}{j} - \frac{\epsilon}{16} \lb
       \Phi\ind{n-1}{j+2} - 4 \Phi\ind{n-1}{j+1} + 6 \Phi\ind{n-1}{j}
                          - 4 \Phi\ind{n-1}{j-1} + \Phi\ind{n-1}{j-2} \rb
	\label{MODEL17}
\end   {equation}
%
where $\epsilon > 0$ is an adjustable parameter which
must satisfy $\epsilon < 1$
for stability.  A precisely analogous update is applied to $\Pi\ind{n+1}{j}$. 
The value of $\epsilon$ is implemented as a run-time parameter, {\tt EPSDIS},
which is read from the control file {\tt ad.control}.
(See {\tt subroutine adddis()} in {\tt unigrid.f} and calls to 
{\tt adddis()} in {\tt ad.f} and {\tt ctres.f}).
 
\ad\ also uses the following discretization of~(\ref{MODEL2}b) to maintain 
an approximation of the scalar field $\phi$ itself:
%
\begin {equation}
	\Pi\ind{n}{j} = \frac{\phi\ind{n+1}{j} - \phi\ind{n-1}{j}}{2 \dt}
\label{MODEL18}
\end   {equation}
%
This yields algebraic equations for the update of the mesh quantities 
$\phi\ind{n+1}{j}$:
%
\begin {equation}
	\phi\ind{n+1}{j} := \phi\ind{n-1}{j} + 2 \dt \Pi\ind{n}{j}
\label{MODEL19}
\end   {equation}
%
which are implemented in {\tt subroutine elfph()} in {\tt unigrid.f}.

Finally, chiefly for diagnostic purposes, \ad\ computes approximations to 
(i) the radial energy density, $dm/dr$, (ii) the spatial integral, $m$ 
(known as the {\sl mass aspect\/} function)  of $dm/dr$ and (iii) 
$m /r$.  The following discrete forms are used (see {\tt subroutine 
clcmas()} in {\tt unigrid.f}):
%
\begin {eqnarray}
%
	{\frac{dm}{dr}}^{n+1}_{j} & := & 2 \pi r\ind{2}{j} \, \lb
		{\Phi\ind{n+1}{j}\,}^2 + {\Pi\ind{n+1}{j}\,}^2 \rb
\label{MODEL20}
\\
	m\ind{n+1}{j}  &  := & m\ind{n+1}{j-1} + \frac{1}{2}    
      \lb r\ind{}{j} - r\ind{}{j-1} \rb \lb {\frac{dm}{dr}}^{n+1}_{j} -
      {\frac{dm}{dr}}^{n+1}_{j-1} \rb
\label{MODEL21}
\\
	{\frac{m}{r}} ^{n+1}_{j}  & := &  \frac{m\ind{n+1}{j}}{r\ind{}{j}}
\label{MODEL22}
%
\end   {eqnarray}
%
\section{Grid Description and Data Structures}
%
The program \ad maintains a tree structure of nested, aligned, uniform meshes.
They are uniform in the sense of scale:  all meshes have the same ratio
$\Delta t / \Delta r$.  Each grid is characterized by the level, L, of
discretization where $L = 1$ is the coarsest level of discretization. An 
example grid with three levels of refinement ($L = 1,2,3$) is:

\vspace{1.0 cm}

\epsfxsize=5in \epsfbox{graphics.dir/threerefine.eps}

\vspace{1.0 cm}

The data structure used to keep track of all of the subgrids is an ordered
tree data structure, where each node represents a grid.  Each node can 
also have multiple descendents.  In this representation, we see that
a node will have multiple descendents if the coarse grid corresponding to 
that node has several fine grids containded in it.  

\section{Integration Algorithm}

The evolution scheme is executed by starting at the coarsest grid level and
working down through to the finest grid level.  The algorithm is
best demonstrated by example.  Lets say we wanted to evolve a grid with
refinement shown in Fig. 1.  The initial spatial surface has the
following structure:

\vspace{1.0cm}

\epsfxsize=5in \epsfbox{graphics.dir/init.eps}

\vspace{1.0cm}

The first level grid ($L = 1$) is evolved, using the finite difference
updating equations:

\vspace{1.0cm}

\epsfxsize=5in \epsfbox{graphics.dir/level1.eps}

\vspace{1.0cm}

The next level ($L = 2$) is evolved, using extrapolated values of the field
on the coarse grid as boundary values for the fine grid:

\vspace{1.0cm}

\epsfxsize=5in \epsfbox{graphics.dir/level2.eps}

\vspace{1.0cm}

At any time during the updating procedure where a fine grid evolves the same
point as a coarser grid, the value found during the fine grid evolution
takes precedence and overwrites the field value found during the
coarser grid evolution.  Finally the last level ($L = 3$) is evolved, again
using the extrapolated values of the field on the coarser grid it is contained
in as boundary values:

\vspace{1.0cm}

\epsfxsize=5in \epsfbox{graphics.dir/level3.eps}

\vspace{1.0cm}

\section{Grid Refinement: Flagging Points that Need Refinement}

Deciding which points are in need of refinement is done by approximating
the local truncation error using Richardson extrapolation.  The basic idea
is to integrate a grid using only half of the spatial resolution as
would normally be used (i.e. ignoring every other spatial grid point).
 The result is then compared to the result obtained
by integrating using the normal spatial resolution (using all grid points).
The results are pointwise compared,
and if the difference between the two is larger
than some predecided tolerence level, the point is flagged for regridding.

\vspace{1.0cm}

\epsfxsize=5in \epsfbox{graphics.dir/regrid.eps}

\vspace{1.0cm}

%
\begin{thebibliography}{99}
%
\bibitem{BergerOliger}
M.~J.~Berger and J.~Oliger,  {\sl J. Comp. Phys}, {\bf 53}, 484 (1984)
%
\bibitem{MattPRL}
M.~W.~Choptuik, {\sl Phys. Rev. Lett.,} {\bf 70}, 9--12 (1993)
%
\bibitem{KO}
H.-O. Kreiss and J. Oliger, {\bf Methods for the Approximate Solution
of Time Dependent Problems}, GARP Publications Series No. 10, (1973)
%
\end  {thebibliography}
%
\end {document}
