%\documentstyle[11pt]{article} 

%\input{rep12}
%\input{symbollga}
%\input{psfig}

%\begin{document}
%\setcounter{page}{22}

\newpage
 
\begin{center}
\Large{\bf IV. A BRIEF TUTORIAL}
\end{center}   
\addcontentsline{toc}{section}{\protect\numberline{IV.}{\bf ~A Brief Tutorial}}

\vspace*{5mm}
{\flushleft\large{\bf IV.A. Installing TITAN}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.A.}{ ~Installing TITAN}}
 
\noindent\rm 
In order to build TITAN one needs all its source files. They come in 2 packages:
one comprises all the subroutines as given in the flow chart, the other contains
various mathematical software in particular from BLAS and LINPACK for the 
Newton-Raphson solver. In order to obtain the first release versions of TITAN 
one has to file a registration with the LCA (laboratory for computational 
astrophysics). The user will be given access privileges to a machine where the
code resides. Begin by making the directory {\tt titan} in your local file space. 
Then signon and change into the directory {\tt /afs/ncar/projects/lca/codes/titan}. The source code is 
located there. Transfer the compressed tar-file {\tt titan.tar.Z} in binary mode from 
the directory {\tt source}. A protocol of these actions looks like:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=4.1in,width=\textwidth,silent=}
\noindent\tt 
\% mkdir titan ; cd titan\\
\% ftp landrew.ncsa.uiuc.edu\\
Connected to landrew.ncsa.uiuc.edu.\\
220 landrew.ncsa.uiuc.edu FTP server \\
Name (landrew.ncsa.uiuc.edu:username): titan\\
331 Password required for titan.\\
Password:\\
230 User titan logged in.\\
ftp> cd /afs/ncsa/projects/lca/codes/titan\\
250-Please read the file README\\
250-  it was last modified on Wed Feb  2 17:22:30 1994 - \dots\\
250 CWD command successful.\\
ftp> ls\\
200 PORT command successful.\\
150 Opening ASCII mode data connection for /bin/ls.\\
.index\\
README\\
physica\_paper\\
ref\_manual\\
user\_guide\\
source

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.7in,width=\textwidth,silent=}
\noindent\tt 
226 Transfer complete.\\
55 bytes received in 0.0057 seconds (9.4 Kbytes/s)\\
ftp> cd source\\
250 CWD command successful.\\
ftp> binary\\
200 Type set to I.\\
ftp> get titan.tar.Z\\
200 PORT command successful.\\
150 Binary data connection for titan.tar.Z \\
226 Binary Transfer complete.\\
local: titan.tar.Z remote: titan.tar.Z\\
150367 bytes received in 13 seconds (12 Kbytes/s)\\
ftp> bye\\
221 Goodbye.
\vspace{2mm}

\noindent\rm
Now the source code resides on your local directory {\tt titan}. The individual
{\tt FORTRAN} files need to be extracted. This is achieved by the following command:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.3in,width=\textwidth,silent=}
\noindent\tt 
\%ls \\
titan.tar.Z\\
\%zcat titan.tar.Z | tar xvf -\\
x advectc.f, 2419 bytes, 5 tape blocks

\noindent\tt 
\hspace{15mm} \vdots \\
x Makefile, 3158 bytes, 7 tape blocks
\vspace{2mm}

\noindent\rm
There is a {\tt Makefile} which installs TITAN on a UNIX system.  
TITAN can be compiled by issuing:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.2in,width=\textwidth,silent=}
\noindent\tt 
\%make \\
f77  -sun4 -c  advectc.f\\
advectc.f:

\noindent\tt 
\hspace{12mm} advectc:

\noindent\tt 
\hspace{15mm} \vdots \\
f77 -O advectc.o \dots ~-o titan.x
\vspace{2mm}

\noindent\rm 
TITAN is now ready to run. Its executable is called {\tt titan.x}. TITAN has 
been compiled on the following platforms:
\bct
\begin{tabular}{|lll|lr|}
    \hline
    {\tt SUN} & {\tt ~~Sparc~~2} & {\tt SunOS~4.1.1}\\
    {\tt SUN} & {\tt ~~Sparc~10} & {\tt SunOS~4.1.3}\\
    {\tt ~HP} & {\tt Apollo~735} & {\tt HP-UX~~~9.1}\\
    {\tt SGI} & {\tt IRIS~4D/35} & {\tt IRIX~~4.0.5}\\
    {\tt SGI} & {\tt Indigo~2~E} & {\tt IRIX~~4.0.5}\\
%   {\tt Cray} & {\tt Y-MP~~~~~} & {\tt UNICOS~7.0.5}\\
    \hline
\end{tabular}
\ect

\noindent\rm
The code needs the Cray vector merge function {\tt cvmgt} (and its relatives). If the user runs TITAN on 
a Cray these will automatically be provided. Otherwise the user must provide them or load {\tt craypack} 
which we supply. If the code is run on a Cray machine, the code requires some CAL-coded linear algebra 
subroutines ({\tt trslbl}, {\tt bglsdc}, {\tt bglsrd}) which are available on Cray only. Again the source 
code for these routines can be provided where necessary. The code also uses the standard LINPACK routines 
such as {\tt sgemv}, {\tt sgemm}, etc. This package is generally available on most systems. If not we can 
provide for the source code. Finally, the code utilizes BLAS (basic linear algebra routines) like 
{\tt saxpy}, {\tt sdot}, etc. extensively. Most systems will provide these automatically. If they are 
not present we have the {\tt FORTRAN} code. \vspace{2mm}

\noindent
TITAN's {\tt Makefile} is written for the Unix operating systems on workstations. For
the SUN it works straightforwardly and the command ``make'' generates the executable.
For the HP and SGI the user has to comment out line 42

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=0.25in,width=\textwidth,silent=}
\noindent\tt
\#\$(RL)                   \$(BOX) \$(LIB)
\vspace{2mm}

\noindent\rm 
to suppress the ``ranlib'' command which is automatically issued by the archive command.
%In order to compile TITAN on a Cray system all the double precision variables and 
%constants have to be changed to single precision. 


\vspace*{5mm}
{\flushleft\large{\bf B. The Input Files and How to Run a Problem}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.B.}{ ~The Input Files and How to Run a Problem}}
 
{\flushleft\large{\sl B.1. The Input Files}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.B.1.}{\sl ~The Input Files}}

\noindent\rm
TITAN is ready to run all test problems with the help of input files called 
{\tt titan.in\#\#}, the hash marks referring to the problem number {\tt lprob}. The code 
assumes the default name {\tt titan.in}. For instance, the input file for Sod's shock tube,
{\tt titan.in03}, must therefore be copied into this file before TITAN can be started. We 
have constructed the input files in a particular way. The input file {\tt titan.in03} 
looks like this:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.4in,width=\textwidth,silent=}
\noindent\tt 
\hspace{5mm} 2                                                              \hfill{\rm i5}

\noindent\tt 
\hspace{5mm} 6 \hspace{3mm} output \hspace{0mm} SOD.out \hspace{8mm} new    \hfill{\rm i5, 2x, a8, 2x, a8}
 
\noindent\tt 
\hspace{5mm} 7 \hspace{7mm} dump \hspace{0mm} SOD.dmp \hspace{8mm} new    \hfill{\rm i5, 2(2x, a8)}
 
\noindent\tt 
\hspace{5mm} 3 \hspace{3mm} 1                                   \hfill{\rm 2i5}
 
\noindent\tt 
           SOD SHOCK TUBE PROBLEM                               \hfill{\rm a80}

\noindent\tt 
\hspace{10mm}     1.0 \hspace{9mm} 1.0 \hspace{9mm} 0.0         \hfill{\rm 3f10.0}
 
\noindent\tt 
\hspace{6mm}   0.125  \hspace{9mm} 0.1 \hspace{9mm} 0.0         \hfill{\rm 3f10.0}
 
\noindent\tt 
\hspace{5mm} 1 \hspace{3mm} 1 \hspace{2mm} 10 1.000e+00  \hfill{\rm 3i5, e10.3}
 
\noindent\tt 
\hspace{10mm} ydl \hspace{9mm} ypl \hspace{9mm} yul
 
\noindent\tt 
\hspace{10mm} ydr \hspace{9mm} ypr \hspace{9mm} yur\\
c\\
c--------------------------------------------------- titan.in03
\vspace{2mm}

\noindent\rm 
The rows are read formatted in the specific format as indicated. The last line always 
identifies the name of the input file. The problem number {\tt lprob=3} appears in the 
fourth row. The following row gives a problem header which appears in the output and dump 
files. Their names are defined in lines following the first one. The rows following the 
header give particular physical parameters and job control parameters for the current run 
followed by comments eluding to the meaning of the physical ones. \vspace{2mm}

\noindent\rm
The number in the first row states how many files (aside from the {\tt tty}) should be 
opened, in this case two. In the second and third row they are defined by a unit number 
({\tt iunit = 2} or {\tt 3}), file type ({\tt output} or {\tt dump}), file name 
({\tt SOD.out} or {\tt SOD.dmp}), and file status ({\tt new}). In the fourth row the 
problem number is specified as defined in the discussion of the test problems. It is 
followed by the initialization flag ({\tt linit = 1}). The problem header is given in 
the fifth row. The eighth row lists the dump file index and time step index from where the 
computation should be started ({\tt jdump = 1} and {\tt jstep = 1}), the number of time 
steps that should be executed ({\tt nstep = 10}), and the maximum time (age) until when the
calculation should be evolved ({\tt tmax = 1.0}). The contents of these rows is 
standardized and necessary for a successful start of {\tt titan.x}. What is listed 
between the header and the job control line is optional and defined according to the 
needs of the particular problem. In our example the sixth and seventh row contain sets 
of three parameters which set the left and right density, pressure, and velocity.\\

\noindent\rm
The following table lists the problem number, the default file name, and the
problem header for each test:
\bct
\begin{tabular}{|l|l|l|lr|}
    \hline
    {\tt lprob} & {\tt name} & {\tt ~~~~~~~~~header}\\
    \hline
    {\tt ~~3~~} & {\tt ~SOD} & {\tt SOD SHOCK TUBE PROBLEM}\\
    {\tt ~~4~~} & {\tt BLST} & {\tt WOODWARD-COLLELA BLAST WAVE PROBLEM}\\
    {\tt ~~5~~} & {\tt SDOV} & {\tt TAYLOR-SEDOV BLAST WAVE PROBLEM}\\
    {\tt ~~6~~} & {\tt HEAT} & {\tt HEATING TOWARDS RADIATIVE EQUILIBRIUM}\\
    {\tt ~~7~~} & {\tt COOL} & {\tt COOLING FROM RADIATIVE EQUILIBRIUM}\\
    {\tt ~~8~~} & {\tt SUBC} & {\tt SUBCRITICAL SHOCK PROBLEM}\\
    {\tt ~~9~~} & {\tt SUPC} & {\tt SUPERCRITICAL SHOCK PROBLEM}\\
    {\tt ~10~~} & {\tt RLST} & {\tt RADIATING BLAST WAVE PROBLEM}\\
    \hline
\end{tabular}
\ect
\vspace{5mm}
 
{\flushleft\large{\sl B.2. Starting and Restarting a Problem}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.B.2.}{\sl ~Starting and Restarting a Problem}}

\noindent\rm
The user has to select a test problem and copy the respective input file
into the default one. In our example of Sod's shock tube we type:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=0.5in,width=\textwidth,silent=}
\noindent\tt 
\%cp titan.in03 titan.in \\
overwrite titan.in? y\\

\noindent\rm
If we issued the command {\tt cat titan.in03} we would see the listing of the
previous page.\vspace{2mm}

\noindent\rm
Now we are ready to run TITAN:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.7in,width=\textwidth,silent=}
\noindent\tt 
\%titan.x\\
iter= ~2 time= 0.00E+00 d\_r/ro(42)= 1.03E-02 rhs(97)= 5.65E-03\\
iter= ~3 time= 1.00E-02 d\_r/ro(42)= 9.74E-03 rhs(97)= 3.92E-03

\noindent\tt 
\hspace{20mm} \vdots\\
iter= 77 time= 7.50E-01 d\_r/ro(10)= 1.61E-08 rhs(73)= 4.50E-11\\
iter= 78 time= 7.60E-01 d\_r/ro(10)= 9.06E-09 rhs(73)= 2.29E-11

\noindent\tt 
model jdump =     1 jstep =    ~1 archived\\
\\
\\
\\

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=3.0in,width=\textwidth,silent=}
\noindent\tt 
jstep =    ~2 dtim =  1.00E-05 tim =  1.00E-05 iter =  1\\
kmax =  \hspace{12mm} 73 \hspace{12mm} 72 \hspace{12mm} 71 \hspace{12mm} 70 \hspace{12mm} 69\\
dmax = ~3.69E-07 ~5.65E-08 ~9.43E-04 ~1.17E-02 ~7.08E-04\\
cmax = ~2.84E-03

\hspace{20mm} \vdots\\
jstep =    11 dtim =  1.90E-05 boost iteration \\
kmax =  \hspace{12mm} 74 \hspace{12mm} 73 \hspace{12mm} 72 \hspace{12mm} 71 \hspace{12mm} 70\\
dmax = -1.41E-09 -1.42E-10 ~3.08E-06 ~4.52E-06 ~1.66E-06\\
cmax = ~1.40E-04\\
\\
total energy check\\
total = 1.3749758E+00 tee = 1.3749758E+00 tes = 0.0000000E+00

\noindent\tt
\hspace{2mm} tew = 0.0000000E+00 tel = 0.0000000E+00 teq = 0.0000000E+00\\
\\
model jdump =     2 jstep =    11 archived\\

\noindent\rm
TITAN starts out by relaxing the adaptive grid from an equidistant distribution to
one that resolves the initial (density and temperature) discontinuity. This it achieves
in 78 iteration steps as indicated in the first column. The second one lists a 
non-physical time which advances by $10^{-2}$ per iteration step. In the third column
the greatest relative temporal change of the radius variable, 
$\left(r^{n+1}-r^n\right)/r^n$, at the corresponding grid index is given. We see that it 
improves from an initial $10^{-2}$ to less than $10^{-8}$ at which point we say the 
relaxation scheme to be converged. In the last column we list the magnitude of the
right hand side of the grid equation which ideally is zero when solved by the 
Newton-Raphson scheme. TITAN makes use of this preiteration to redistribute the grid
points in order to resolve the desired features for the problems with {\tt lprob} 
$\in \tt (2, 4, 5, 6)$, \ie~for all hydrodynamic problems and the first radiative
transport problem. The grid relaxation procedure is completed with writing the
resulting model into the output and dump file.\vspace{2mm}

\noindent\rm
The model with {\tt jdump = 1} and {\tt jstep = 1} serves as the initial condition from
which the partial differential (rather, finite difference) equations evolve. This is
documented in the listing by the fact that the first time step begins at {\tt jstep = 2}.
The initial time step, here {\tt dtim} $= 10^{-5}$, is set by the user. It should be 
considerably smaller than any physical time scales of the problem in order to allow the
code to compute beyond transients that could stem from the initialization procedure.
Each time step takes at least 2 iteration cycles and at most 10 to converge. The 
convergence threshold for {\tt dmax} has been set to $10^{-6}$. In our example the eighth 
and ninth time step needed only three and the eleventh time step six iteration cycles. \vspace{2mm}

\noindent\rm
Each iteration step lists the biggest value of the right hand side of each equation in 
the row beginning with {\tt dmax =} and the corresponding grid index {\tt kmax} in the 
row above. The last number {\tt cmax =} gives the maximum cell size as evaluated in 
{\tt maxdel}. When the code has converged for a particular time step it performs one 
boost iteration to improve the accuracy of the solution. The convergence control is 
described in the TITAN Code Reference Manual, chapter XVII. At the end of each time step 
the cumulative energy conservation and balancing is listed under {\tt total} {\tt energy} {\tt check}. 
This is done as described in the TITAN Code Reference Manual, chapter XI, with the 
correspondences as follows: total energy {\tt total}, sum of internal, radiation, 
kinetic, and potential energy {\tt tee} $= {\cal E}^{n+1}$, the net energy loss by 
transport through the boundary surfaces {\tt tes} $= {\cal S}^{n+1}$, the work done at 
the boundary surfaces {\tt tew} $= {\cal W}^{n+1}$, the luminosity lost through the 
boundary surfaces {\tt tel} $= {\cal L}^{n+1}$, and the energy dissipation by the 
artificial viscosity {\tt teq} $= {\cal Q}^{n+1}$. These quantities are evaluated in the
routine {\tt totnrg}.\vspace{2mm}

\noindent\rm
TITAN stops the computation either when the desired number of time steps (here: 
{\tt nstep = 10}) have been executed or when the desired model age (here: 
{\tt tmax = 1.0}) has been reached. It terminates properly by archiving the last model 
into the dump file and writing it into the output file while stating the corresponding 
{\tt jdump} and {\tt jstep}.\\

\noindent\rm
Now we want to restart TITAN at {\tt jdump = 2} and {\tt jstep = 11} and run Sod's shock
tube until the time equals unity. In order to do that we have to edit the input file
{\tt titan.in} and change the third, fourth, and eighth line in the following manner:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.3in,width=\textwidth,silent=}
\noindent\tt 
\hspace{20mm} \vdots

\noindent\tt 
\hspace{5mm} 7 \hspace{7mm} dump \hspace{0mm} SOD.dmp \hspace{8mm} old
 
\noindent\tt 
\hspace{5mm} 3 \hspace{3mm} 2

\noindent\tt 
\hspace{20mm} \vdots

\noindent\tt 
\hspace{5mm} 2 \hspace{1mm} 11 \hspace{0mm} 310 1.000e+00

\noindent\tt 
\hspace{20mm} \vdots
\vspace{2mm}

\noindent\rm
First we changed the status of the dump file from {\tt new} to {\tt old} so that TITAN
can read the last model from the previous run. By setting {\tt linit = 2}  we tell the 
code to bypass the initialization procedure and continue straight with the evolution of
the equations. In the eighth line we replace the first two integers with the appropriate
settings for {\tt jdump} and {\tt jstep} as given at the end of the previous run, and
change {\tt nstep} to {\tt 310} and set {\tt tmax = 1.0}. Now we type {\tt titan.x} again
so see the following listing at our monitor:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=4.1in,width=\textwidth,silent=}
\noindent\tt 
\\
jstep =   ~12 dtim =  2.30E-04 tim =  9.77E-04 iter =  1\\
kmax =  \hspace{12mm} 58 \hspace{12mm} 57 \hspace{12mm} 56 \hspace{12mm} 55 \hspace{12mm} 54\\
dmax = ~7.06E-05 ~3.16E-06 -5.85E-03 -2.70E-03 ~3.41E-02\\
cmax = ~1.47E-01

\noindent\tt 
\hspace{20mm} \vdots\\
jstep =   306 dtim =  1.88E-02 tim =  1.01E+00 iter =  5\\
kmax =  \hspace{12mm} 76 \hspace{12mm} 75 \hspace{12mm} 74 \hspace{12mm} 73 \hspace{12mm} 72\\
dmax = ~4.35E-06 ~5.10E-06 ~7.02E-06 ~8.66E-06 ~2.83E-06\\
cmax = ~3.00E-04\\

\noindent\tt 
jstep =   306 dtim =  1.88E-02 boost iteration\\
kmax =  \hspace{12mm} 76 \hspace{12mm} 75 \hspace{12mm} 74 \hspace{12mm} 73 \hspace{12mm} 72\\
dmax = -9.34E-08 -1.09E-07 -1.32E-07 -1.78E-07 -6.18E-08\\
cmax = ~5.17E-06\\
\\
total energy check\\
total = 1.3593973E+00 tee = 1.3593973E+00 tes = 0.0000000E+00

\noindent\tt
\hspace{2mm} tew = 0.0000000E+00 tel = 0.0000000E+00 teq = 1.6190176E-19\\
\\
model jdump =     8 jstep =   306 archived\\

\noindent\rm 
TITAN completed the calculation in 306 time steps and made 8 dumps. The last time
step has converged after 5 iteration cycles followed by the boost iteration. Whenever
convergence is reached quickly it is achieved quadratically as expected from the 
Newton-Raphson algorithm. This is always true for a fixed grid. For the adaptive grid this
cannot happen all the time especially if the correct grid distribution is not found right
away. A look at the total energies at {\tt jstep = 11} and {\tt 306} tells that we
have conserved it to within one percent. We save the output file to {\tt SOD.out1}.
TITAN wrote a model every five time steps in the format as given in {\tt writout}. There
are eight columns: the grid index {\tt k}, the radius, grid resolution, density, artificial
viscous stress, gas energy, and gas pressure. At four different times we extract the 
density profile and show them as snap shots in figure 1 ($+$'s indicate grid points).\vspace{2mm}

\psfig{figure=sod1.ps}
\centerline{\small Fig. 1: Density profiles of classical shock tube with adaptive grid}
\vspace{2mm}

\noindent\rm 
The user can design new shock tube problems by changing the physical parameters
in the input file. Let the density jumps be twice as large as in Sod's case. We have to 
change the first number in the sixth line of {\tt titan.in} to {\tt ydl = 2}. 
The following figure illustrates the result in form of density profiles at four different
times ($+$'s indicating again individual grid points). Comparing this with figure 1 we see 
that the new run completes at a later time. The sound speed, proportional to 
$\sqrt{\Pg/\rho}$, is a factor $\sqrt{2}$ larger thus making the shock to propagate slower 
by the inverse factor.\\
 
{\flushleft\large{\sl B.3. Hints for Running the Test Problems}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.B.2.}{\sl ~Hints for Running the Test Problems}}

\noindent\rm
In the spirit of the previous section the user can run all the other test problems. It is 
useful to break up the computation into several runs in order to avoid the output files 
{\tt \#\#\#\#.out} to become too big. This is especially important if one wishes to write 
out many models, every fifth time step say, to produce space-time diagrams or movies. The 
dump and output frequencies are set in the {\tt start\#\#} files.\vspace{2mm}

\psfig{figure=sod4.ps}
\centerline{\small Fig. 2: Density profiles of modified shock tube with adaptive grid}
\vspace{2mm}

\noindent\rm
The input files of the problems {\tt lrpob = 8} and {\tt 9} have only one difference
in the seventh line which states the physical parameters. For the subcritical test we
have set {\tt upstn} $= 6\times 10^5$ whereas it is $2\times 10^6$ for the supercritical
test. The start files for both tests are identical and a duplicate exists only for 
completeness reasons. Either test could be run with the problem number of the 
other.\vspace{2mm}

\noindent\rm
Some care must be taken when running the two radiative transfer tests. The heating towards
radiative equilibrium needs only {\tt nstep = 150} for the relaxation computation. The 
results for the different opacities are obtained by varying the last physical parameter
{\tt ratio} in the input file. The cooling tests require the final model as an input. Since 
this is the same for all heating runs the user can take any such dump file called
{\tt HEAT.dmp} and copy it into {\tt COOL.dmp}. The cooling test then must start at the 
last dump ({\tt jdump = 31}) and time step index ({\tt jstep = 151}). It proceeds for a 
default {\tt nstep = 400} and can be restarted by changing the {\tt linit} flag and the
job control parameters {\tt jdump} and {\tt jstep}.


\vspace*{5mm}
{\flushleft\large{\bf C. I/O-Files}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.C.}{ ~I/O-Files}}

\noindent\rm
TITAN permits access to the nine input/output files.  All of the unit numbers and file 
types are defined in {\tt opnfil}. Only {\tt itty} and {\tt iin} are assigned to non-zero 
and hardwired to {\tt 4} and {\tt 5} there. The other unit numbers can be assigned to any 
different integer. The file names are all eight character long with the default 
{\tt ********}. The file type (format {\tt a8}) is used to identify the unit number and 
file status with the file name. All files are formatted, with file status being either {\tt new} 
or {\tt old}, except the dump file which is unformatted. The available files are given in the 
following table:\vspace{2mm}
\bct
\begin{tabular}{|l|l|l|l|lr|}
    \hline
    {\tt ~unit} & {\tt file type} & {\tt file name}     & {\tt ~~~~access}\\
    \hline
    {\tt ~itty} & {\tt ~~monitor} & {\tt ~/dev/tty}     & \\
    {\tt ~~iin} & {\tt ~~~~input} & {\tt ~titan.in}     & \\
    {\tt ~iout} & {\tt ~~~output} & {\tt ~\#\#\#\#.out} & {\tt sequential}\\
    {\tt idump} & {\tt ~~~~~dump} & {\tt ~\#\#\#\#.dmp} & {\tt ~~~~direct}\\
    {\tt ~ieos} & {\tt ~~~~~~eos} &                     & {\tt sequential}\\
    {\tt iopac} & {\tt ~~~~~opac} &                     & {\tt sequential}\\
    {\tt ihist} & {\tt ~~history} & {\tt ~\#\#\#\#.sty} & {\tt sequential}\\
    {\tt iinit} & {\tt ~~initial} & {\tt ~\#\#\#\#.ini} & {\tt sequential}\\
    {\tt ~idoc} & {\tt ~~~~~~doc} & {\tt ~\#\#\#\#.doc} & {\tt sequential}\\
    \hline
\end{tabular}
\ect
\vspace{2mm}

\noindent\rm
The files for {\tt eos} and {\tt opac} define the constitutive relations in tabular form 
and come either from OP or OPAL. They are opened in {\tt opnfil} and called by their 
respective subroutines. In a future release some will be available and carry the file names 
as follows:\vspace{2mm}
\bct
\begin{tabular}{|l|l|lr|}
    \hline
    {\tt      file name} & {\tt contents} \\
    \hline
    {\tt ~\#\#\#\#\#.pg} & {\tt gas pressure} \\
    {\tt ~\#\#\#\#\#.eg} & {\tt internal energy} \\
    {\tt ~\#\#\#\#\#.pe} & {\tt electron pressure} \\
    {\tt ~\#\#\#\#\#.cf} & {\tt Rosseland mean opacity} \\
    {\tt ~\#\#\#\#\#.ke} & {\tt Planck mean opacity} \\
    \hline
\end{tabular}
\ect
\vspace{2mm}

\noindent\rm
The user can add further files if so desired by editing {\tt opnfil} and {\tt clzfil}
and copying the structure of these routines.
 
%\vspace*{5mm}
{\flushleft\large{\bf D. The Start Routines and How to Modify a Problem}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.D.}{ ~The Start Routines and How to Modify a Problem}}
 
{\flushleft\large{\sl D.1. The Start Routines}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.D.1.}{\sl ~The Start Routines}}

\noindent\rm
All start routines are constructed in the same fashion. At the top are comment cards
containing the name, purpose, and some additional information. These are followed by the 
three major segments:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.5in,width=\textwidth,silent=}
\noindent\tt
c=============================================================\\
c\hspace{9mm} general setup\\
c=============================================================

\noindent\tt
\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} generate initial solution\\
c=============================================================

\noindent\tt
\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} recover model from dump file\\
c=============================================================

\noindent\tt
\hspace{20mm} \vdots\\


\noindent\rm
In {\tt general setup} all logical switches and all job control and physical parameters 
are specified as defined in the TITAN Code Reference Manual, chapter XIV and XV. 
The general setup has the following layout:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.2in,width=\textwidth,silent=}
\noindent\tt
c-------------------------------------------------------------\\
c\hspace{9mm} read in header and parameters\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} geometry: lgeom = 0  => planar geometry\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} gravity: lgrav = 0  => no gravity\\
c-------------------------------------------------------------

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=7.2in,width=\textwidth,silent=}
\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} radiation transport: lrad  = 0  => no radiation\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} radiation boundary conditions:\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} declare equations to be solved:\\
c-------------------------------------------------------------
\\
\hspace{20mm} \vdots\\
\noindent\tt
c-------------------------------------------------------------\\
c\hspace{9mm} numerics:\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} set hydro parameters: lhydr = 1  =>    hydro\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} hydrodynamic boundary conditions:\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} grid specification: lgrid = 1  => adaptive   grid\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} iteration control \& integration control:\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} output control:\\
c-------------------------------------------------------------\\


\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.9in,width=\textwidth,silent=}
\noindent\tt
\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} eos: leos = 2  => perfect gas\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} opacity: lopac = 4  => constant opacity\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\

\noindent\rm
First the header and the physical parameters are read from the input file {\tt titan.in}.
In {\tt geometry} we specify slab (0), cylindrical (1), or spherical (2) symmetry. The
quantity $\mu =$ {\tt lgeom} and its relatives as used in {\tt peq\#\#} and the definition 
of the equations is set. In {\tt radiation transport} we set the flag {\tt lrad} either on
or off. In the first case we need to specify what transfer scheme we would like. That is
done by setting {\tt ltran} along with {\tt lam}. We further have to give the number of
rays that intersect the core {\tt ncor} for the Eddington factor calculation, \cf~section
X. Finally the radiant boundary conditions have to be defined via {\tt lribc} and 
{\tt lrobc}. In the case of them being equal to unity the outgoing/incoming specific
intensity incident on the inner/outer boundary, {\tt xipl} $= I_L^+$ to {\tt ximr} $= I_R^-$,
need to be given. In {\tt declare equations to be solved} the number of grid points in the
computational domain {\tt ncell} must be set along with the number of equations at each
grid point {\tt neqn}. The include file {\tt titan.par} sets upper limits on these numbers,
{\tt ncell} $\le$ {\tt mgr - 9} and {\tt neqn} $\le$ {\tt meqn}, as well as for {\tt ncor} 
$\le$ {\tt mcor}. Each equation is assigned a dummy index. For adaptive grid computations
it is meaningful to set {\tt ir = 1} because the code has to find the coordinate system.
In {\tt set hydro parameters} there is a series of hydrodynamical control parameters:
\bct
\begin{tabular}{|l|l|l|lr|}
\hline
{\tt  parameter name}      & {\tt symbols}                              & {\tt equation}\\
\hline
{\tt time-centering parameter}          & {\tt theta} $= \theta$                 & {\rm (14)} \\
{\tt pseudoviscous length} & {\tt ql0} $= \ell_0$, {\tt ql1} $= \ell_1$ & {\rm (25),}~{\tt FE70}\\
{\tt pseudoviscous coefficients}        & {\tt cq1} $= C_1$, {\tt cq2} $= C_2$   & {\rm (25),}~{\tt FE83}\\
{\tt pseudoviscous pressure ratio}      & {\tt q0} $= q_0$ & {\tt AG87} \\
{\tt artificial diffusion coeff.}       & {\tt sigd} $= \sr$, {\tt side} $= \se$ & {\tt C42} \\
    \hline
\end{tabular}
\ect
\vspace{2mm}

\noindent\rm
Furthermore, the order of the advection scheme has to be defined: for {\tt cadv = 1} we
choose Donor-cell and for {\tt cadv = 2} we choose Van Leer advection. When we decide for
a specific hydrodynamic boundary condition we have to set:
\bct
\begin{tabular}{|l|l|l|lr|}
    \hline
    {\tt boundary}   & {\tt switch}       & {\tt constraint}\\
    \hline
    {\tt Lagrangean} & {\tt llibc}        & {\tt uextl} = $U_L$ \\
		     & {\tt llobc}        & {\tt uextr} = $U_R$, {\tt pextr} $= \Pi_R$\\
    {\tt Eulerian}   & {\tt leibc}        & {\tt phil\{0,1,2\}} = $\Phi_L^{0,1,2}$ \\
                     & {\tt leobc}        & {\tt phir\{0,1,2\}} = $\Phi_R^{0,1,2}$ \\
    \hline
\end{tabular}
\ect
\vspace{2mm}

\noindent\rm
In {\tt grid specification} the user has to set the grid flag {\tt lgrid}: the choice is
among an adaptive grid (1), a fixed Eulerian grid (2), or a fixed Lagrangean grid (3). In
the first case a series of grid control parameters need to be specified:
\bct
\begin{tabular}{|l|l|l|lr|}
    \hline
    {\tt parameter name}   & {\tt symbols}       & {\tt equation}\\
    \hline
    {\tt diffusion coefficient} & {\tt alph} $= \alpha$ & {\tt AG5}\\
    {\tt time constant}         & {\tt tau} $= \tau$    & {\tt AG4}\\
    {\tt delay exponent}        & {\tt ibet} $= \beta$  & {\tt AG4}\\
    {\tt spatial scale}         & {\tt xscale} $= R_{scale}$  & {\tt AG7}\\
    {\tt ordinate scales}       & {\tt yscl(l)} $= \left(F_{scale}\right)_l$  & {\tt AG9}\\
    {\tt grid weights}          & {\tt wt(l)} $\in (W_m, W_\rho, W_T, W_E,$   & \\
				& ~~~~~~~~~~~~~~~$W_p, W_e, W_\kappa, W_q)$   & \\
    \hline
\end{tabular}
\ect
\vspace{2mm}

\noindent\rm
The user has a choice between three settings for the scalings of the grid quantities: for
{\tt ladx = 1} together with {\tt xscale} $= R_{scale}$ one choses a linear spatial 
resolution ({\tt AG7}), for {\tt ladx = 2} a logarithmic spatial resolution ({\tt AG8}), 
for {\tt lady(l) = 1} together with {\tt yscl(l)} $= \left(F_{scale}\right)_l$ one choses
a linear ordinate resolution ({\tt AG9}), for {\tt lady(l) = 2} a logarithmic ordinate 
resolution ({\tt AG10}), and for {\tt lady(l) = 3} a harmonic ordinate resolution 
({\tt AG11}). In general one can reduce the pseudoviscous length scale by at least one
order of magnitude in the adaptive grid over the fixed grid case.\vspace{2mm}

\noindent\rm
For {\tt iteration} {\tt control} {\tt \&} {\tt integration} {\tt control} we recommend the user to stay 
with the parameters as preset. On the other hand, the user should choose her/his own settings for 
{\tt output control}: {\tt ndump} is the number of timesteps between dumps of models and to 
{\tt \#\#\#\#.dmp} and {\tt nout} is the number of timesteps between printouts of models to 
{\tt \#\#\#\#.out}. The dump frequency can be chosen generously since it is used to store 
possible locations for a job restart. The write frequency depends on whether the user 
wishes to get a quick overview over the solution or wants to do a detailed study.\vspace{2mm}


\noindent\rm
Finally there are a few parameters to set in {\tt eos} and {\tt opac}. In general the 
switch settings {\tt leos = 1} and {\tt lopac = 1} give access to the constitutive 
relations in tabular form. The user has to provide for her/his own tables until we give out 
some OP data. For the equation of state one has also the option of an ideal gas 
({\tt leos = 2}), in case of which the polytropic exponent {\tt gam} $= \gamma$ and the 
mean molecular weight {\tt gmu} $= \mu_m$ have to be given, or of the Stellingwerf fit 
formula for variable star envelopes {\tt leos = 3}, in case of which the chemical 
abundancies {\tt xabun} $= X$, {\tt yabun} $= Y$, {\tt zabun} $= Z$ have to be set 
(X, Y, Z = 1 - X - Y are not allowed to be arbitrary numbers). For the opacities one also 
has the option of a constant law {\tt lopac} $\le$ {\tt 3}, in case of which the ratio 
between the Rosseland and Planck mean opacity {\tt ratio} $= \fr$ must be defined (see 
II.E.3), or of the Stellingwerf fit formula for variable star envelopes {\tt lopac = 2}, 
in case of which chemical abundancies together with the electron number density 
{\tt xnen(k)} $= \Ne(\rho, T)$ must be specified.\\

\noindent\rm
The part {\tt generate} {\tt initial} {\tt solution} is executed for the setting {\tt linit = 1}.
It is particular to the individual test problem. Accordingly the various sections differ 
with {\tt lprob}. In them the initial model for the time dependent equations is constructed 
in a consistent fashion. For our example {\tt start03} this is achieved in the following 
steps under the heading {\tt initialize} {\tt variables}: we begin with setting up the 
{\tt radial} {\tt grid}, then define the $r$-dependent physical quantities {\tt density} {\tt and} 
{\tt temperature}, and {\tt mass} {\tt and} {\tt velocity}, along with the {\tt phantom} {\tt zones}, 
followed by the evaluation of the {\tt eos}, after which we are ready to {\tt initiate} {\tt old} {\tt time} 
{\tt solution}, to {\tt initialize} {\tt terms} {\tt in} {\tt total} {\tt energy} {\tt equation}, to 
{\tt zero} {\tt unneeded} {\tt vari- ables}, to {\tt zero} {\tt initial} {\tt errors}, to {\tt initialize} 
{\tt time} {\tt and} {\tt timestep}, and to {\tt initialize} {\tt grid} {\tt equation} optionally for 
{\tt lgrid = 1}. In the case that we opted for the adaptive grid we have to reset various quantities which 
is done in {\tt zero} {\tt unneeded} {\tt variables} again, {\tt initialize} {\tt terms} {\tt in} {\tt total}
{\tt energy} {\tt equation}, {\tt mass} {\tt and} {\tt velocity}, and finally {\tt initialize} {\tt time} 
{\tt} {\tt and} {\tt timestep} again. The initial model is now archived in {\tt write} {\tt initial} 
{\tt model} {\tt on} {\tt dump} {\tt file} and assigned the numbers {\tt jstep = 1} and {\tt jdump = 1}. 
The user is asked to study the various initialization procedures for the test problems.\\

\noindent\rm
The start routines are always concluded by {\tt recover} {\tt model} {\tt from} {\tt dump} {\tt file}. 
This section is executed for {\tt linit} $\ne 1$. It begins with reading the problem header
and parameters from the dump file and comparing them with those specified in the input
file. If they agree (otherwise resulting in an error message) the new dump number and
timestep for the starting model are read in from {\tt titan.in}. It is checked that the 
specified starting model exists. The model is read in from the dump file and it is
verified that the code starts from the correct one. Here is the generic layout:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=3.5in,width=\textwidth,silent=}
\noindent\tt
c-------------------------------------------------------------\\
c\hspace{9mm} compare model header\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} read dump number and timestep number of initial model;\\
c\hspace{9mm} check that specified starting model exists\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} read model\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\
c-------------------------------------------------------------\\
c\hspace{9mm} check that we got the right one\\
c-------------------------------------------------------------

\noindent\tt
\hspace{20mm} \vdots\\

\vspace*{5mm}
{\flushleft\large{\sl D.2. Modifying a Problem}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.D.2.}{\sl ~Modifying a Problem}}

\noindent\rm
We can change the default job control and physical parameter settings of all test
problems in the respective {\tt start\#\#} routines. Again the hash marks refer to the
problem number. For illustrative purposes we employ again Sod's shock tube as our example
and {\tt start03} is the start routine under consideration. For instance, we might be
interested in repeating the complete run with a fixed (Eulerian) grid to compare it with
the adaptive result. We edit {\tt start03} and forward to line 180 where the section of
the grid specification begins. There we see that the logical switch {\tt lgrid} is used
to choose between the adaptive ({\tt 1}) and Eulerian ({\tt 2}) grid. Since we are now
interested in the last one we have to modify line 206 to {\tt lgrid = 2}. We recompile
this new setting by typing {\tt make titan}. Before we can run the Eulerian version of
the problem we first have to reset the input file {\tt titan.in}. First the status of
the dump file is changed back to {\tt new}, then {\tt linit = 1}, and finally we set
{\tt jdump = 1} and {\tt jstep = 1} in the correct formatting. We leave the number of
time steps {\tt nstep = 310} because we wish to run the evolution until {\tt tmax = 1}.
Now we can carry out the computation:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=5.8in,width=\textwidth,silent=}
\noindent\tt 

\noindent\tt 
model jdump =     1 jstep =    ~1 archived\\
\\
jstep =   ~~2 dtim =  1.00E-05 tim =  1.00E-05 iter =  1\\
kmax =  \hspace{12mm} 99 \hspace{12mm} 98 \hspace{12mm} 97 \hspace{12mm} 96 \hspace{12mm} 95\\
dmax = ~0.00E+00 ~8.01E-16 -5.37E-21 ~1.95E-17 ~4.44E-14\\
cmax = ~1.00-300\\
\\
jstep =   ~~2 dtim =  1.00E-05 boost iteration\\
kmax =  \hspace{12mm} 99 \hspace{12mm} 98 \hspace{12mm} 97 \hspace{12mm} 96 \hspace{12mm} 95\\
dmax = ~0.00E+00 -5.34E-17 -5.37E-21 ~4.16E-20 ~3.47E-16\\
cmax = ~1.00-300\\
\\
total energy check\\
total =  1.3750000E+00 tee =  1.3750000E+00 tes = ~0.0000000E+00

\noindent\tt 
\hspace{2mm} tew =    0.0000000E+00 tel =  0.0000000E+00 teq = -9.5764228E-34

\noindent\tt
\hspace{20mm} \vdots\\
\noindent\tt
jstep =   311 dtim =  1.89E-03 tim =  3.49E-01 iter =  1\\
kmax =  \hspace{10mm} 100 \hspace{12mm} 99 \hspace{12mm} 98 \hspace{12mm} 97 \hspace{12mm} 96\\
dmax = ~0.00E+00 -1.93E-05 ~1.57E-02 ~1.01E-01 -6.62E-03\\
cmax = ~1.00-300\\
\\
jstep =   311 dtim =  1.89E-03 tim =  3.49E-01 iter =  2\\
kmax =  \hspace{10mm} 100 \hspace{12mm} 99 \hspace{12mm} 98 \hspace{12mm} 97 \hspace{12mm} 96\\
dmax = ~0.00E+00 ~4.43E-08 ~7.39E-05 ~8.21E-04 -1.87E-03\\
cmax = ~1.00-300\\
\\
jstep =   311 dtim =  1.89E-03 tim =  3.49E-01 iter =  3\\
kmax =  \hspace{12mm} 85 \hspace{12mm} 84 \hspace{12mm} 83 \hspace{12mm} 82 \hspace{12mm} 81\\
dmax = ~0.00E+00 -7.61E-13 -3.12E-10 ~2.34E-12 -1.83E-11\\
cmax = ~1.00-300
\\

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.0in,width=\textwidth,silent=}
\noindent\tt
jstep =   311 dtim =  1.89E-03 boost iteration\\
kmax =  \hspace{12mm} 77 \hspace{12mm} 76 \hspace{12mm} 75 \hspace{12mm} 74 \hspace{12mm} 73\\
dmax = ~0.00E+00 ~5.73E-17 ~3.26E-17 ~3.62E-16 ~4.66E-15\\
cmax = ~1.00-300\\
\\
\noindent\tt
total energy check\\
total =  1.3687542E+00 tee =  1.3687548E+00 tes =  0.0000000E+00

\noindent\tt
\hspace{2mm} tew =    0.0000000E+00 tel =  0.0000000E+00 teq = 1.3583605E-22\\
\\
model jdump =     8 jstep =   311 archived\\

\noindent\rm 
At the beginning TITAN takes only two iteration steps per time step until it increases to 
about $10^{-3}$. The Newton-Raphson then converges strictly quadratically in three to four
iterations. We see that {\tt cmax =} $10^{-300}$ at every time step which is the default 
setting in {\tt maxdel}. This comes from the fact that the right hand side of the grid 
equation now is always identical to zero as seen in the first entry of the {\tt dmax} row. 
The Eulerian version has reached the time 0.349 in {\tt 311} time steps. We therefore want 
to continue the calculation. For this we have to change {\tt titan.in} again in the way 
as described in section IV.B. but this time we set {\tt nstep = 999} to make sure that
the run will be completed. We execute {\tt titan.x} and after another 354 time steps the
code reaches {\tt tmax = 1} and terminates with {\tt jdump = 16} and {\tt jstep = 665}.
We now can look at {\tt SOD.out} containing the Eulerian models and compare them with
{\tt SOD.out1} where we saved the adaptive grid models.\vspace{2mm}

\noindent\rm 
In figure 3 we show the density profiles of this fixed grid computation. The four times 
were chosen to be close to the ones in figure 1. We observe that the density jumps lack 
the sharpness of this figure (or figure 2 for that matter). Already at early times the 
shock, contact discontinuity, and rarefaction fan are not well resolved. About 2 grid 
points define the shock and about 6 the contact discontinuity. In contrast to that the
propagating shock front contains between 10 and 20 grid points and the contact 
discontinuity about 15 in the adaptive grid computation. The ``extra'' grid points were 
drawn from the other regions and from the rarefaction fan in particular, which is 
represented by considerably fewer grid points in the adaptive grid computation. Our 
Eulerian calculation proceeds satisfactory (upper right plot) until the shock reflects 
from the wall. At this instant (lower left plot) the numerical effects from wall heating 
are visible. The corresponding spike stays

\psfig{figure=sode.ps}
\centerline{\small Fig. 3: Density profiles of classical shock tube with fixed grid}
\vspace{5mm}

\noindent\rm 
through the interaction with the contact discontinuity. The matter could be helped by 
increasing the coefficients for artificial mass and energy diffusion $\sr$ and $\se$ in 
the lines before the definition of the hydrodynamic boundary conditions.\\
 
\vspace*{5mm}
{\flushleft\large{\bf E. Rules for Writing Your Own Problem}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.E.}{ ~Rules for Writing Your Own
Problem}}

\noindent\rm
The radiation hydrodynamics equations in TITAN can be modified, augmented, or completely
replaced with your own set of equations. Examples for this are: addition of conductive
terms for heat transport, a complete $O(v/c)$ formulation of the radiation hydrodynamical 
equations allowing fluid flows $v \approx 0.1 c$, a general relativistic formulation, 
addition of species equations for multi fluid problems for chemical or nuclear reaction
networks, addition of transfer equations for non-photonic energy fluxes, replacement by
models for cosmic ray acceleration or any other set of nonlinear partial differential
equations. Corresponding to the changes is the amount of work that goes into adapting
TITAN for your problem of choice. In the simplest case,\ie~using the equations as provided,
one has be concerned with the initial model only. In the case of adding terms or equations
one also needs to do coding of the right hand sides and entries into the Jacobian matrix. 
The modular nature of TITAN's architecture allows one to achieve this with the smallest
effort possible. In general, the code itself should always be refered to for guidance.\\
 
{\flushleft\large{\sl E.1. Writing Your Own Start Routine}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.E.1.}{\sl ~Writing Your Own Start Routine}}

\noindent\rm
A new initial model has to be generated in different ways depending on whether it can
be formulated analytically via formulae or whether it can be defined by a set of static
equations or whether is has to be obtained from mapping a provided model onto TITAN's
finite difference scheme. Examples for the first are all test problems, for the second
are the stellar structure equations, and for the third is a machine readable model from
a colleague.\\

\noindent\rm
We have provided for two dummy start routines, {\tt start01} and {\tt start02}, where
the user can state her/his problem of choice. They contain a template for the first segment,
{\tt general setup}. The user has to fill in the appropriate numbers. For the proper choices 
one can refer to the examples in the other start routines and the TITAN Code Reference Manual, 
chapters XIV and XV. The user also has to make a decision on which parameters should be read in. 
This is done in the section {\tt read in header and parameters}. There the problem header must 
be read and it is written onto the monitor and output file. \vspace{2mm} 

\noindent\rm
The section where the initial model is archived carrying the heading {\tt write} {\tt initial} {\tt model} 
{\tt on} {\tt dump file} and the last segment {\tt recover} {\tt model} {\tt from} {\tt dump} {\tt file} 
exist already. The user can optionally write to/read from the dump file as the first record ({\tt irec = 1}) 
a desired selection of parameters together with the {\tt header}. In the test problems the parameters are 
those which are read in from the input file. Under the section {\tt compare} {\tt model} {\tt header} the 
parameters are again user defined. Those from the restart model are compared with those from the inital 
model to check the correctness of the restart model. Suppose that we have read in the variables {\tt xmass}, 
{\tt xlum}, {\tt xrad}, {\tt teff}, from {\tt titan.in}. In the start routine we therefore have to define:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=7.0in,width=\textwidth,silent=}
\noindent\tt 

\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} write initial model on dump file\\
c=============================================================\\
c

\hspace{5mm} if ( idump .le. 0 ) return\\
c

\hspace{5mm} irec = 1

\hspace{5mm} write(idump, rec = irec, iostat = ier)

\hspace{3mm} . \hfill header, xmass, xlum, xrad, teff

\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} recover model from dump file\\
c=============================================================\\
c

\hspace{5mm} if ( linit .ne. 1 ) then\\
c\\
c-------------------------------------------------------------\\
c\hspace{9mm} compare model header\\
c-------------------------------------------------------------\\
c

\hspace{5mm} irec = 1

\hspace{5mm} read~(idump, rec = irec, iostat = ier)

\hspace{3mm} . \hfill yheader, ymass, ylum, yrad, yteff

\hspace{20mm} \vdots

\hspace{5mm}   if ( yheader .ne. header .or. ymass .ne. xmass .or.

\hspace{3mm} .   ~~~ylum~~~ .ne. xlum~~ .or. yrad~ .ne. xrad~ .or.

\hspace{3mm} .   ~~~yteff~~ .ne. teff~~~~~~) ~then\\
c

\hspace{5mm}   ~~~write(itty, '(a80/a80/1p4e15.7/1p4e15.7)')

\hspace{3mm} . \hspace{30mm} header, yheader,

\hspace{3mm} . \hspace{30mm} xmass, xlum, xrad, ~teff,

\hspace{3mm} . \hspace{30mm} ymass, ylum, yrad, yteff,

\hspace{7mm}   ~~~stop 'start21'

\hspace{7mm}  end if\\

\noindent\rm
In the main section, {\tt generate} {\tt inital} {\tt solution}, the part of the coding already exists 
which is necessary for a proper start of the evolution equations. This encompasses the pieces {\tt zero} 
{\tt unneeded} {\tt variables}, {\tt radial grid}, {\tt mass and velocity}, {\tt radiation field}, 
{\tt phantom} {\tt zones}, {\tt eos}, {\tt opacity}, {\tt edding- ton} {\tt factors}, {\tt initialize} 
{\tt terms} {\tt in} {\tt total} {\tt energy} {\tt equation}, {\tt initialize} {\tt old} {\tt time} 
{\tt solution}, {\tt initialize} {\tt time} {\tt and} {\tt timestep}, and {\tt zero} {\tt initial} 
{\tt errors}. These sections are intended to follow the user's specification of the inital model
including the grid relaxation. Their main purpose is to reset the physical quantities properly after
an optional grid relaxation and therefore must be executed in the order given. We recommend to follow this 
procedure in any case. In {\tt radial} {\tt grid} and {\tt mass} {\tt and} {\tt velocity} some important
quantities like the cell volume {\tt dvoln(k)} and cell mass {\tt xmn(k)} are re-evaluated and the
velocity field is set to zero, {\tt u(k) =} 0. For radiative problems the {\tt radiation} {\tt field}
is re-evaluated. So are the {\tt eos}, {\tt opacity}, and {\tt eddington} {\tt factors} which are also
needed for the evaluation of the {\tt total} {\tt energy} of the initial model. In {\tt initialize} 
{\tt time} {\tt and} {\tt timestep} the beginning timestep {\tt dtime} for the evolution needs to be 
set. In general we advise to make it a very small number because the code quickly relaxes to the
physical time scale. \\
 
{\flushleft\large{\sl E.2. Hints For Generating the Initial Model}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.E.2.}{\sl ~Hints For Generating the Initial Model}}
 
\begin{center} \it Reference \end{center}
 
{\flushleft\rm
\hspace{10mm}M. Steffen, \it Astron. Astrophys., \bf 239, \rm 443, 1990}\\
\vspace*{2mm}

\noindent\rm
Most of the main section of the start routine is devoted to the definition of the initial model. 
For all practical purposes it can be described either analytically by formulae or by a set of 
ordinary differential equations or is given in numerical form. All test problems belong to the first 
category whereby the last one ({\tt lprob = 10}) employs a model relaxation. In general the user 
should formulate the initial model first on a fixed (Eulerian or Lagrangean) grid and relax it there. In 
the second step, if necessary or desirable, this solution is mapped onto the adaptive grid and relaxed 
on that. In all hydrodynamic test cases we had to relax the adaptive grid to resolve the prescribed 
discontinuities (see V.B). On the other hand, in all radiation hydrodynamic tests we stated the initial 
model on the Eulerian grid and used that to start up the evolution. This was possible because the 
initial models were smooth. Their evolution induced steep features which were resolved through the 
parabolic nature of the adaptive grid equation. \vspace{2mm}

\noindent\rm
The following flow diagram summarizes the initialization procedure:\\
\flowchart 
\vspace{4.2in}

\begin{description}
\item Some recommendations on how to proceed if the initial model is given

\item``by formulae'':\\
Define either the radius or the mass variable equidistantly. Express as many physical quantities as
possible in terms of the radius or mass, and all others as functions of the first. As an example may 
serve $\rho(r) = \rho_0 \left(r/r_0\right)^7$ which translates immediately into
$\rho_k^n = \rho_0 \left(r_k^n/r_0\right)^7$. The user should pay close attention to the functional 
dependencies if the same formula occurs in the evolution equations. For instance, the mass in a shell 
is written in the form 
$\Delta m = 4\pi \rho \Delta\left(r^{\mu+1}/{\mu+1}\right)$, \cf~chapter III.B of the TITAN Code
Reference Manual, which has a different truncation error than 
$\Delta m = 4\pi \rho r^\mu \Delta(r)~$! Coding the proper formulae constitutes already a good starting
model for a fixed grid computation. If it contains sharp features then the adaptive mesh should be 
allowed to resolve them in a relaxation procedure with a selection of grid parameters. As described in the
previous section (E.1) some important quantities need to be re-evaluated on the new non-uniform grid 
distribution. The model is ready to evolve with that set of grid parameters.
       
\item``by a set of equations'':\\
Define either the radius or the mass variable equidistantly. Write down the equations with one of them
as the independent variable. In the case of ordinary differential equations it is crucial to formulate 
them in a way that is consistent with the formulation of the partial differential equations. {\em In
praxi} this means that one uses the identical finite difference approximations. The resulting expressions
might look somewhat awkward but they admit exactly the same truncation errors as the evolution equations.
Being extremely attentive to this point saves a lot of debugging hassle later on. As an example may serve
the radiative flux in diffusion approximation. Its formula is given in equation (7), see chapter II.A. Our 
finite difference approximation, as found in chapter IX.A of the TITAN Code Reference Manual, is derived 
from 
$$
     \Fr = -{c\over3} {1\over\rho\kF} {\p\er\overp\left({r^{\mu+1}\over{\mu+1}}\right)}   \eqno{(54)}$$ 
which is analytically equivalent to (7) but not numerically! Our rule of thumb in formulating the finite 
difference expressions is to leave them as compact as possible. Relax or integrate the equations to find the 
fixed grid solution of the initial model. The integration is best carried out with the Newton-Raphson scheme.
If one needs to resolve physical structures then the adaptive grid has to be relaxed. This means that now 
the independent variable is viewed as a dependent one and all quantities become functions of the grid index 
$k$. The grid equation with a set of grid parameters must be solved together with the given equations. As 
soon as a desired grid distribution is obtained the initial model can serve as an input into the evolution
equations.\\

\item``in numerical form'':\\
The physical quantities are given in dependency on the grid index. They must be transformed into the set 
of fundamental variables ($r$, $m$, $\rho$, $u$, $\Fr$, $T$, $\er$) that is used by TITAN. Further, they 
must be mapped onto the finite difference approximations employed by TITAN (see the examples above). 
Finally, the need to map between a different amount of grid points in the computational domain can arise. 
For this purpose we provide for an interpolation routine which utilized monotonized piecewise cubic 
polynomials according to M. Steffen. If the given numerical description of the initial model is steady-state 
then the mapped model should be relaxed first on a fixed grid and then, if necessary, together with the grid 
equation. If the given numerical description of the initial model is a snap shot from an evolution process, 
\ex~containing a particular velocity field, then a relaxation of the mapped model might not be possible and 
it has to be fed into TITAN's evolution equations ``as is''. In this case one should try to find the best 
possible mapping satisfying the above criteria. A test for the quality of the mapping is how well the model 
evolves physically as well as numerically.
\end{description}
\vspace{5mm}
 
{\flushleft\large{\sl E.3. Relaxing the Adaptive Grid}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.E.3.}{\sl ~Relaxing the Adaptive Grid}}

\noindent\rm
There are several ways to relax the adaptive grid if the initial model contains sharp features 
or discontinuities. Here we shall describe a method that works reliable and well. This relaxation
scheme is employed in our test problems.\vspace{2mm}

\noindent\rm
The basic principle is as follows: Decide by which physical quantities the discontinuity can be described
most easily. They will become part of our selection of the grid parameters. Express the discontinuity in 
terms of them as a function of the radius or mass variable. Write this down in form of one or several 
steady-state equations, possibly by introducing dummy quantities. These equations are to be coupled with 
the grid equation. Since it possesses a time derivative one can use the standard Newton-Raphson iteration 
to perform a relaxation of the adaptive grid. Once the grid distribution is satisfactory the relaxation
is terminated.\vspace{2mm}

\noindent\rm
We illustrate this method using Sod's shock tube as an example. The complete relaxation scheme has been
coded outside of {\tt start03} in the file {\tt gridinit}. The first routine there, {\tt grid3nit}, is
called by {\tt start03}. It contains the parts:\\

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.5in,width=\textwidth,silent=}
\noindent\tt 
c=============================================================\\
c\hspace{9mm} modify original code parameters\\
c=============================================================

\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} newton-raphson iteration\\ 
c=============================================================

\hspace{20mm} \vdots\\
c=============================================================\\
c\hspace{9mm} restore original code parameters\\
c=============================================================

\hspace{20mm} \vdots\\

\noindent\rm
In the first part we select temperature and density to express the initial discontinuity in terms of the
radius variable. Together with the grid equation we therefore want to iterate three equations. Since the 
grid equation is implicit we have to solve for $\tt r(k)$ via a matrix inversion. TITAN provides already
for a complete linear solution scheme which we can employ here as well. To do this we set {\tt neqn = 3}
and assign three equation indices {\tt ir = 1}, {\tt id = 2}, {\tt it = 3}. Then we choose as grid 
parameters the density, pressure, and internal energy in logarithmic scaling with weights equal unity.
The grid diffusion coefficient $\alpha$ is left as set in {\tt start03} but the grid time constant $\tau$
is arbitrarily set to {\tt tau = }$10^{-2}$. So is the time step, {\tt dtime =}$10^{-2}$, which serves
only as a counter in the relaxation method. Their ratio $\tau/\delta t$ thus becomes equal to unity
which forces the grid equation to redistribute the zones. In the last part the particular selection
of code parameters are nullified by resetting them to the original ones.\vspace{2mm}

\noindent\rm
In the middle part, {\tt newton-raphson} {\tt iteration}, the grid relaxation with our selection of grid 
parameters is carried out. It follows the general solution procedure as outlined in {\tt step}. Under the 
heading {\tt calculate} {\tt derivatives} we describe the density and temperature discontinuity by a very 
steep, yet smooth function. We choose the tanh-function because it allows us to control the location and 
the slope of the step. Having the left and right density prescribed as {\tt ydl} and {\tt ydr} we can 
express the density step as a function of the radius in one compact formula:\\
$$
     \rho(r) = {\tt ydr} + {1\over2}({\tt ydl}-{\tt ydr})
	       \left[ 1 - {\rm tanh}\left( {r - {\tt rc}\over{\tt rw} }\right)\right]    \eqno{(55)}$$
The quantities {\tt rc} and {\tt rw} give the location and the steepness of the step. From the formula 
(55) we see that the density approaches a step function with the discontinuity at $r =$ {\tt rc} for 
{\tt rw} $\rightarrow\infty$. We use an artificial viscosity to describe the emerging shock. It is 
therefore meaningful to relate {\tt rw} to the shock width which is given by the pseudoviscous length; in
{\tt start03} we simply set {\tt rw} $=\ell_0$. The temperature step is also expressed in the form (55).
Furthermore, (55) admits a straightforward derivative through the cosh-function:
$$
     {d\rho(r)\over{d r}} = - {1\over2}({\tt ydl}-{\tt ydr})
               {1\over{\tt rw}} {\rm cosh}^{-2}\left( {r - {\tt rc}\over{\tt rw} }\right) \eqno{(56)}$$
The next section, {\tt initialize} {\tt all} {\tt matrices} {\tt to} {\tt zero}, the relaxation iteraion 
begins with resetting all entries of the Jacobian matrix and the auxiliary solution arrays to zero. This is 
identical to the beginning of {\tt matgen}. This is followed by defining two equations through their right 
hand sides, \ex~{\tt rhs(k) = dn(k) - }$\rho${\tt (rn(k))}, and elements of the Jacobian matrix. Then the 
grid equation is called. The system of three equations is solved on the computational domain by {\tt matslv}.
We use only the radius array from the solution, contained in {\tt rhs(k)}, and correction vector, $\delta r$.
We evaluate its largest element, called {\tt dmax}, and compute the largest grid-spacing, called {\tt cmax}. 
If {\tt dmax} or {\tt cmax} exceed the preset values {\tt dtol} or {\tt ctol} we scale $\delta r$ by the 
factor {\tt min(1., dtol/dmax, ctol/cmax)}. The new radial distribution is computed from the solution vector,
\ie~{\tt rhs(k)}, and the old radial distribution {\tt rn(k)}. The dummy time step is advanced by {\tt dtime}
and our two equations are re-evaluated. We employ as the convergence criterion the relative change in the 
grid distribution between two relaxation steps: $\max_k\left(r_k^{n+1}/r_k^n-1\right) < 10^{-8}$. By 
non-fulfillment the relaxation cycle is continued with {\tt initialize} {\tt all} {\tt matrices} {\tt to} 
{\tt zero}.\vspace{2mm} 

\noindent\rm
This grid relaxation appears to be somewhat cumbersome because it necessitates at least one more equation
in addition to the grid equation. All we care about is to redistribute the grid points to resolve the
discontinuity. Thus we actually only have to relax the grid equation by itself. It, however, needs 
input from the discontinuity via a the physical variables, such as density and temperature and their spatial
derivatives. The description of the discontinuity, however, depends on the specific problem that one
wishes to investigate. Iterating only the grid equation thus means that one has to code it differently
for each problem. We thought that this is way too complicated and prone to admit many mistakes. In order
to avoid this we decided to keep the grid equation as programmed and add a description of the discontinuity
in equation form. Thus we only have to specify that for a given problem while having the solution procedure 
already at hand.\\
 
{\flushleft\large{\sl E.4. Coding Your Own Equations}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.E.4.}{\sl ~Coding Your Own Equations}}

\noindent\rm
In principle, any number of new terms, and any number of new equations can be implemented
in TITAN. An example for the first case is the addition of electronic heat conduction, and
for the latter is the addition of an equation for the turbulent kinetic energy. In both
cases the user has to define new global quantities and list them as a common block in the
include file {\tt titan.com}. Since adding a new term is a special case of adding a new
equation we shall discuss some rules that need to be observed in implementing the latter.\\

\noindent\rm
Suppose that we wish to add an equation for the turbulent kinetic energy that is of the form
$$
     \rho\Dt(\et) + {\p (r^\mu \Ft)\over{r^\mu\p r}}
                + \Pt{\p(r^\mu u)\over r^\mu\p r}
                + \rho \Bigl[\Kt - \St\Bigr] = 0    \eqno{(57)} 
$$
with the notation $\et$ for the turbulent kinetic energy, $\Ft$ for the turbulent flux, 
$\Pt$ for the turbulent pressure, and $\St$ and $\Kt$ for sink and source terms. Before 
we can code equation (57) we have to transform it to arbitrary coordinates in the spirit 
of chapter II B. The adaptive mesh version reads:
$$
   {d\over d t} \int_V \et~\rho dV - \oint_{\p V}~\et \rho\ur~dS
       + \int_V d[r^\mu~\Ft]   
       + \int_V \Pt~d[r^\mu\rho~u] $$
$$     + \int_V~\Bigl[\Kt-\St\Bigr]~\rho dV = 0     \eqno{(58)}
$$
Expressed in symbolic notation:
$$
   \frac{d}{dt} \left( \et \rho\Delta V \right) 
                   - \Delta\left[\frac{dm}{dt}~\et - r^\mu\Ft \right] 
                   + \Pt\Delta(r^{\mu}u) 
                   = \Bigl[\St - \Kt]~\rho\Delta V  \eqno{(59)} 
$$
The translation from this to the actual finite difference representation is completely
analogous to that of any equation (17)--(21) with the help of the TITAN Code Reference Manual.
Let the {\tt FORTRAN} names of the new quantities be {\tt et} $= \et$, {\tt ft} $= \Ft$, {\tt pt}
$= \Pt$, {\tt source} $= \St$, and {\tt sink} $= \Kt$.\\

\noindent\rm
The new fundamental quantity is $\et$ since we determine it from equation (59). All the other 
quantities are algebraic combinations of $\et$ with the existing quantities $\rho$, $\eg$, 
$\er$, $u$, $\Fr$, $r$, and the auxiliary quantity $m$. In our example we have the formulae
$\Pt = {2\over3}\rho \et$ and $\Ft = - \Lambda {\p(\et^{3/2})\overp r}$. We won't need the 
definitions of the other quantities for the purpose of our discussion. Thus, instead of seven 
equations we now deal with eight, and correspondingly we have eight fundamental quantities that 
are iterated in the Newton-Raphson algorithm. Those and their related quantities are specified
in {\tt titan.com}. Most of them appear in the first columns of the common blocks with the names 
{\tt common /mid/}, {\tt common /new/}, and {\tt common /old/}, referring to the time-centered 
quantities ($x$), the quantities at the new ($x^{n+1}$) and old ($x^n$) time level, respectively. In
our case we have listed the time-centered {\tt et} in the row below {\tt er} and above {\tt fr} 
(in line 130 to be precise). Similarly, we add the turbulent kinetic energy at the new time level
{\tt etn} $= \et^{n+1}$ to {\tt common /new/}, and {\tt eto} $= \et^n$ to {\tt common /old/}
(lines 145 and 157). Since $\et$ is also advected we denote the advected {\tt et} by {\tt etb}.
Turbulent kinetic energy is a scalar quantity and therefore cell-centered. The appropriate advection
routine is {\tt advectc}. It employs the advected quantity at both time levels and its monotonized 
slopes at both time level, see chapter IV, equations $C22$--$C36$, in the TITAN Code Reference Manual. 
We have to provide for storage space for the monotonized slopes of the advected $\et$ which we call 
{\tt etsn} and {\tt etso}. We enter them also in the common blocks {\tt /new/} and {\tt /old/}. 
Furthermore, we list {\tt ft}, {\tt pt}, {\tt sink}, and {\tt source} in {\tt /mid/}. The new common 
block {\tt /mid/} now reads:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.7in,width=\textwidth,silent=}
\noindent\tt 
\hspace{6mm} common /mid/

\noindent\tt 
\hspace{5mm} . \hspace{5mm} r~~~~ (mgr), rmu~~ (mgr), rmup1 (mgr), rmum1 (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} dvol~ (mgr), xm~~~ (mgr), d~~~~ (mgr), db~~~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} u~~~~ (mgr), ub~~~ (mgr), unom~ (mgr), as~~~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} eg~~~ (mgr), egb~~ (mgr), pg~~~ (mgr), xne~~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} er~~~ (mgr), erb~~ (mgr), egr~~ (mgr), egrb~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} et~~~ (mgr), etb~~ (mgr), egt~~ (mgr), egtb~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} fr~~~ (mgr), frb~~ (mgr), egrt~ (mgr), egrtb (mgr),

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=1.3in,width=\textwidth,silent=}
\noindent\tt 
\hspace{5mm} . \hspace{5mm} frnom (mgr), avchi (mgr), pt~~~ (mgr), ambda (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} ft~~~ (mgr), fc~~~ (mgr), sink~ (mgr), source(mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} ftnom (mgr), fcnom (mgr), floor (mgr), acous (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} chif~ (mgr), xke~~ (mgr), xkp~~ (mgr), plf~~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} beta~ (mgr), dels~ (mgr), cp~~~ (mgr), fedd~ (mgr),

\noindent\tt 
\hspace{5mm} . \hspace{5mm} t~~~~ (mgr), geddl~~~~~~, geddr\\

\noindent\rm
It might be necessary to dedicate a separate common block for the definition of turbulent quantities.
In our case a candidate for this is the turbulent flux {\tt ft} $= \Psi$ since we need this quantity
in the turbulent energy routine, in the internal energy routine which contains {\tt egrt} 
$= \eg+\er/\rho+\et$, and in the total energy routine {\tt totnrg}. In {\tt titan.com} we demonstrate this
in {\tt common /turb1/} which contains all derivatives of {\tt ft} with respect to all fundamental
quantities at various grid points, \ex~{\tt dftdld00} $= \p${\tt ft}$/\p\ln(${\tt dn}$)|_k$ or {\tt dftdlrp1}
$= \p${\tt ft}$/\p\ln(${\tt rn}$)|_{k+1}$. Finally we need to provide for additional equation indices 
{\tt ic} and {\tt jc} in {\tt common /index/}. \vspace{2mm}

\noindent\rm
Now we are ready to code the turbulent energy equation in the new file {\tt turnrg}. An easy way of 
getting most terms correctly is to copy the contents of {\tt gasnrgh}. We begin by changing all equation
indices {\tt it} $\rightarrow$ {\tt ic}, {\tt jt} $\rightarrow$ {\tt jc}. Then we code (59) in the section
{\tt right} {\tt hand} {\tt side} in the following manner:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.1in,width=\textwidth,silent=}
\noindent\tt 
\hspace{6mm} rhs(ic,k) =

\noindent\tt 
\hspace{5mm} . \hfill ( dn(k) * etn(k) * dvoln(k)~~~~~~~~

\noindent\tt 
\hspace{5mm} . \hfill - do(k) * eto(k) * dvolo(k) )/dtime

\noindent\tt 
\hspace{5mm} . \hfill - dmdt(k+1) *  etb(k+1)~

\noindent\tt 
\hspace{5mm} . \hfill + dmdt(k~~) *  etb(k~~)~

\noindent\tt 
\hspace{5mm} . \hfill + rmu(k+1) *  ~ft(k+1)~

\noindent\tt 
\hspace{5mm} . \hfill - rmu(k~~) *  ~ft(k~~)~

\noindent\tt 
\hspace{5mm} . \hfill + pt(k) * ( rmu(k+1) *  ~~u(k+1)~

\noindent\tt 
\hspace{5mm} . \hfill ~ ~~~~~ ~ - rmu(k~~) *  ~~u(k~~))

\noindent\tt 
\hspace{5mm} . \hfill + ( sink(k) - source(k) ) * d(k)~ * dvol(k~~)~
\\

\noindent\rm
In the section {\tt set} 
{\tt interface} {\tt energy} {\tt densities} {\tt for} {\tt advection} we have to rename the character
string {\tt eg} $\rightarrow$ {\tt et}. In the section {\tt time} {\tt derivative} we repeat the renaming
and delete the derivative with respect to the density, \ie~{\tt e00(ic,jd,k)}, since $\p\et/\p\rho = 0$.
We proceed analogously in the section {\tt work}: {\tt pg} $\rightarrow$ {\tt pt}. We can keep its
derivatives {\tt dptdldn} $= \p${\tt pt}$/\p\ln(${\tt dn}$)$ and {\tt dptdlcn} 
$= \p${\tt pt}$/\p\ln(${\tt et}$)$ if we define them beforehand. It is advisable to include them in the
common block {\tt common /turb1/} so that they can be referred to in the internal energy equation. Now
we have to compute the Jacobian matrix entries from the flux divergence. We could copy the section
{\tt flux} {\tt divergence} from the file {\tt radnrgr} and change {\tt ie} $\rightarrow$ {\tt ic}.
Since {\tt ft} is a derived variable we have to right down all of its derivatives with respect to all
fundamental quantities at all grid points as listed in {\tt titan.com} under {\tt common /turb1/}. Again
these are best evaluated beforehand at the beginning of {\tt turnrg}. We conclude the main programming
part with the sinks and sources. The section {\tt source}-{\tt sink} {\tt terms} in {\tt radnrgr} is
a good guide to follow. In the section {\tt interior} {\tt zones} of our file {\tt turnrg} we already
have defined the right hand side. The advection of {\tt et} is coded by changing {\tt egb} into
{\tt etb}, the other derivatives are correct. We conclude the coding of {\tt turnrg} by setting all
inner and outer boundary conditions and the phantom zones. \vspace{2mm}

\noindent\rm
If the user introduces a new fundamental quantity then it and its relatives must be listed in the following 
routines:
\begin{description}
\item{\tt start01}: Add {\tt etn}, {\tt etsn} to the write statement to the dump file with record number 
{\tt irec = jdump + 2} in the section {\tt write} {\tt initial} {\tt model} {\tt on} {\tt dump} {\tt file} 
and in the corresponding read statement in the section {\tt read} {\tt model}. Assign the additional 
equation index {\tt ic = 8} and the numbers of equations {\tt neqn = 8} in the section {\tt declare} 
{\tt equations} {\tt to} {\tt be} {\tt solved}. Finally we have to define the flag {\tt ltur = 1} to
admit the turbulence equation.

\item{\tt archive}: Add {\tt etn}, {\tt etsn} to the write statement to the dump file with record number 
{\tt irec = jdump + 2}.

\item{\tt reset}: Add {\tt etn}, {\tt etsn} to the loop {\tt do 2}.

\item{\tt totnrg}: Add {\tt et(k)} = 0.5 ({\tt etn(k)} + {\tt eto(k)}) to the loop {\tt do 2}. The formula 
for {\tt ft(k)} is already available through the common block {\tt /mid/} and need not to be re-computed in 
loop {\tt do 1}. The turbulent flux needs also to be included in the computation of the luminosity.
\end{description}

\noindent\rm
The following two new routines should be written:
\begin{description}
\item{\tt peqtrh}: Copy the routine {\tt peqrh} which prepares the radiation hydrodynamics equations for the 
Newton-Raphson iteration into this new file. Add {\tt eto}, {\tt et}, {\tt etn}, {\tt etso} to the loop
{\tt do 2}. Define inner and outer boundary conditions for {\tt etn}. Evaluate the time-centered 
quantity {\tt et} = $\theta$ {\tt etn} + (1 - $\theta$) {\tt eto} in loop {\tt do 12}.

\item{\tt updatetrh}: Copy the routine {\tt updaterh} which updates the solution of the radiation 
hydrodynamics equations after a Newton-Raphson iteration into this new file. Add the line {\tt etn(k)} = 
(1 + {\tt rhs(ic,k)}) {\tt etn(k)} to the loop {\tt do 12}. Just like in {\tt peqtrh} the inner
and outer boundary conditions for {\tt etn} must be defined and the time-centered {\tt et} must be evaluated.
\end{description}

\noindent\rm
Before we can test the new equation we have to call it from {\tt matgen} along with the other equations:

\psfig{file=box1.ps,rheight=0bp,rwidth=0bp,height=2.4in,width=\textwidth,silent=}
\noindent\tt 
\hspace{6mm} if (lhydr.eq.1~.and.~lrad.eq.1~.and.~ltur.ge.1) then 

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call mass

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call contin

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call viscous

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call turnrg

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call gasmomtrh

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call gasnrgtrh

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call radnrgrh

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call radmomrh

\noindent\tt 
\hspace{5mm} . \hspace{10mm} call grid

\noindent\tt 
\hspace{5mm} . \hspace{10mm} return

\noindent\tt 
\hspace{6mm} end if\\

\noindent\rm
The internal energy equation {\tt gasnrgtrh} is the routine {\tt gasnrgrh} with all the new terms from 
{\tt turnrg}. The turbulent sink-source term is omitted just like the turbulent one because we regard the 
sum {\tt egrt} in this file. Because of the turbulent pressure term we also get a corresponding {\tt pt}
gradient in the momentum equation which is included in {\tt gasmomtrh}. The routines {\tt peqtrh} and 
{\tt updatetrh} have to be called from {\tt peq} and {\tt update} with the logic as shown in the example
above. The implementation of the new equation (57) is now complete.

%\end{document}
.
