%\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} .