Preface
back to Table of Contents
VH-1 was written by the numerical astrophysics group at the University of Virginia in 1990-1991, and is based on the description of the piece-wise parabolic method (PPM) given in Colella and Woodward (1984, hereafter "CW"). While many individuals at UVa contibuted to the development of VH-1, the bulk of the code writing and testing was done by John Hawley, John Blondin, Greg Lindahl, and Eric Lufkin. This User's Guide was written by J.B., and he assumes all responsibility for its content. This text is intended to answer most questions regarding the organization of VH-1 and its use in solving hydrodynamics problems. The User's Guide should therefore eliminate the necessity of contacting one of the code authors for further information.The purpose of VH-1 is to provide the astrophysics community with a working version of PPM that can be adapted by workers familiar with finite differencing to their own research. VH-1 is NOT idiot proof. No attempt has been made to provide a universal code that can be blindly applied to any problem in hydrodynamics. To that end we expect the VH-1 user to be able to change the code to suit his/her own problem. In addition, we assume that the reader of this Guide and user of this code is FORTRAN literate (yes, fortran!), and that he/she will know how to compile, link, and run VH-1 on a local machine. A UNIX makefile is included with the code. If the code was aquired as a concatenation of all the subroutines, it is recommended that they be separated into individual files and compiled with the aid of a Makefile.
This guide to VH-1 will not repeat the detailed description of PPM given in CW, but will merely provide a brief overview of the algorithm in chapter 2. It is imperative that the VH-1 user read and understand CW in addition to reading this manual.
The following chapters provide a brief overview of PPM, a detailed description of VH-1, and instructions on how to run it. We also discuss some general features of the code and algorithm that we have found necessary to pay careful attention to, and that therefore should be of interest to the VH-1 user. These include dissipation techniques and coordinate singularities. A few test problems are included with the code. A description of each of these problems and an example of VH-1 output for each problem is given in S3.
1. PPM: A Brief Review
back to Table of Contents
VH-1 uses finite difference techniques to solve the equations of ideal inviscid compressible fluid flow. These equations, written in conservative, Eulerian form are:
where 
is the
total specific energy, rho is the mass density, p the pressure,
and F and G are momentum and energy source terms (e.g., gravity, radiative
cooling).
VH-1 does not include an energy source term, but such a term can
be readily added using an operator splitting method. Because a force
term is already present in the form of fictitous forces for a curvilinear grid,
we have included an arbitrary body force term in VH-1.
VH-1 is written as a Lagrangian hydrodynamics code coupled
with a remap onto the original Eulerian grid. This implementation of
PPM is referred to as PPMLR. Our primary reason for using PPMLR rather
than the direct Eulerian method is to avoid the subtleties associated with
constructing the correct input states for the Riemann problem in PPMDE
(see pages 190-191 in CW ). While PPMLR requires a second call to the
interpolation routines for the remap step, the computation of the input
states is much simplified, and the Courant constraint on the timestep
is less restrictive, being the minimum of 
In practice we have found PPMLR to
be less dissipative (not always a good thing) than PPMDE, and better at
maintaining contact discontinuities without the aid of a contact
steepener.
The following is an outline of the PPMLR algorithm. For specific
details see CW. As an illustration, let us consider the equations
of ideal hydrodynamics in Lagrangian coordinates in
one dimension, planar geometry and no source terms. Written in
conservative form these equations are:
where V, u, p, and E are the specific
volume, velocity,
pressure, and total energy of the gas. The gas density rho and internal
energy e are related to these quantities via the relations
These conservation equations can be finite differenced as:
where the subscript j refers to zone averaged values, and subscripts j-1/2 and j+1/2 refer to values at the left and right hand sides of the zone, and the superscript is the timestep.
The variables u-bar and p-bar are the time-averaged values of the velocity and pressure at the (Lagrangian) zone interfaces. The trick is clearly to obtain accurate, stable estimates of these quantities.
The approach of Godunov's method is to obtain these time-averaged quantities by approximating the flow at each zone interface during each timestep with a Riemann shock tube problem. At the beginning of the timestep, the zone interface is modelled as a discontinuity separating two uniform states given by the zone averages on the left and right side of the zone interface (zones i-1 and i). This constructed Riemann problem is then solved to find the time-averaged value of the velocity and pressure at this discontinuity (zone interface), u-bar and p-bar. The solution of the non-linear Riemann problem typically requires some kind of iterative procedure, as described in van Leer (1979) and in Woodward (1986).
PPM improves upon this method by using more accurate guesses for the input states to the Riemann problem (the values on either side of the interface). Using a quadratic interpolation of the fluid variables in each zone, the Riemann input states are taken to be the average over that part of the zone that can be reached by a sound wave in a time dt, i.e., the characteristic domain of dependence.
Once the hydrodynamic equations have been differenced to obtain the values at t+dt, the fluid variables can be instantaneously remapped from the Lagrangian coordinate system to the stationary Eulerian grid. This remap step uses the same quadratic interpolation method that was used in the hydrodynamics step.
This description makes PPM sound relatively simple, but there are many key ingredients that we have skipped over, including the parabolic interpolation and monotonicity constraints, the iterative Riemann solver, and additional dissipation mechanisms. We will give a brief discussion of these first two topics in the appropriate subsections of Chapter 2 and postpone the discusion of dissipation until Chapter 4 (see also CW).
2. VH-1
back to Table of Contents
The following description of the code assumes the reader is familiar with finite differencing in general and with PPM in particular. For specific details we will occasionaly refer the reader the various references on PPM (CW and Woodward 1982). In the following two subsections we provide a detailed description of the code, subroutine by subroutine. But first, let us give a general introduction to VH-1. File names (e.g., subroutines) are referred to in bold type, vhone.f, while variable names are given in italic typeface, ncycle. The file names in this HTML file are also WWW-linked to the corresponding FTP file, for speed and convenience.VH-1 is currently written as both a 1D and a 2D code, both versions employing the same subroutines for the hydrodynamics calculation. Essentially, the only difference is that the 2D code calls the 1D hydro sweep for both the X and the Y directions, while the 1D code only calls the hydro sweep for the one dimension. Extension of the code to 3D is relatively trivial given the 2D code, but we recommend becoming familiar with the 1D and 2D versions before attempting 3D computations. Both the 1D and 2D versions are driven by vhone.f and call all of the same subroutines. The dimension of the problem is set in init.f by a global variable, ndim that is used by the subroutines that must know of the dimension of the problem, such as output routines and dtcon.f.
The variables needed to run this code (ignoring for the moment the local variables in the lower subroutines) are stored in three common blocks. global.h contains all the ``global'' quantities needed throughout the code, including boundary conditions, physical parameters, and loop variables. Sweep.h contains all the 1D arrays needed for the one dimensional sweeps, and is included in all subroutines below sweepx.f . In addition, sweepsize.h contains the parameter maxsweep, used to dimension local working arrays in the 1D subroutines. Zone.h contains the multidimensional arrays for the fluid variables (density, pressure, and velocities), grid coordinates, and image interpolation coefficients (if needed for making movie images from a curvilinear grid). A description of each variable and its location in the code is listed in appendix A of the paper manual.
The input used to run the code is written in a namelist file, indat . This file currently provides the prefix for the output file names, the adiabatic index, information on the grid size, the Courant parameter, the number of timesteps to evolve, when to write a dump file, a data file, and a movie file, and the maximum time to evolve the flow. It is expected that the user will adapt this input file for his/her own use.
The output file convention implemented in VH-1 uses a five letter name (e.g., filea) followed by a three letter suffix or a consecutively increasing number. A separate data file is written out each time prin.f is called, and the number used as the filename suffix, nfile, is incremented by one for each data print. Thus the first data file will be filea.1000, and the second data file will be filea.1001. The binary dump files that are written periodically to allow a restart if something goes wrong have the suffix daa with the last two letters incremented by one character each time dump.f is called. The frequency of these outputs are determined by the inputs ndump and nprin. The history file that contains all the pertinent information about the run has the suffix hst, and a time series of HDF raster images is stored in a file with suffix mvy.
A schematic of the code is shown in Figure 1. The code naturally separates into two components: the main driver and the 1D hydro sweep. The main driver performs all of the peripheral activities such as input, output, and timestep control, while the hydro is computed within the 1D hydro sweeps. For multidimensional problems this is done by calling the 1D hydro sweep in alternating directions. The following two subsections will describe these two primary components of VH-1.
2.1 ONE DIMENSIONAL HYDRODYNAMICS
back to Table of Contents
In this subsection we will describe how VH-1 performs the 1D hydro as called by sweepx.f (or sweepy.f). We will defer discussion of the peripheral routines, such as grid.f, init.f, outputs, and drivers to the next section, and focus on the actual hydrodynamics calculation in this section. The convention in this code is to use the index n for the 1D hydro, and (i,j) for the multidimensional arrays, so the following discussion will only use index n.
sweepx.f
In loading the 1D arrays we must create two ghost zones on each end of the 1D array in order to handle boundary conditions. Thus the first real zone is nmin=3, and the last real zone is nmax=imax+2, where imax is the number of zones in the x direction. These boundary conditions must be set for the coordinates before the call to paraset.f, which calculates the geometric coefficients used in the interpolation scheme. These boundary conditions for the Eulerian grid coordinates will be independent of the type of BC, and will normally just assume a constant grid spacing.
Finally the call to paraset.f is made and the loop over each row can be run. To accomodate curvilinear coordinates, the radius of the angular columns is given by radius. For cartesian coordinates, this variable is set to 1 before the loop over rows or columns, and left there.
Each loop over a row begins by loading the 1D arrays from a given row in the 2 or 3D arrays. Note again the padding of two ghost zones on each end of the 1D array. The coordinates must be loaded each time as the previous calculation will have evolved the grid points in the lagrangian step. At the same time we put a floor value on the pressure and calculate the total energy in each zone. The zone volume is then calculated based on the value of the geometry flag ngeomx (ngeomy), and the variable radius is set to the radius of the current vector if the current coordinate is angular.
This procedure of loading up 1D arrays from the 2D arrays is not needed in the 1D simulations, but we retain this structure for the 1D code to maintain compatibility and similarity between the 1D and 2D versions of this code. If extensive work is to be done in 1D, the user may wish to "strip" out the 1D version and remove the unneccesary 2D structure.
Before calculating the hydrodynamics the BCs must be set in the four ghost zones. These are set by a call to sweepbc.f with the current boundary condition flags passed as parameters. Boundary conditions are set for all fluid variables and coordinates in all four ghost zones, except that xa(nmax+1) is not set in sweepbc.f because it is a real quantity rather than a boundary condition (i.e., it is evolved with the real flow). The transverse velocities are not needed in this step (they just advect with the flow), but are nonetheless set in sweepbc.f for simplicity.
The call to ppm.f performs the Lagrangian hydrodynamics step (see following sections), updating the 1D array with new fluid variables and coordinates. The call to ppm.f passes the BC flags, nleft and nright, the geometry, ngeom, and the array of parabolic coefficients, paracoff.
The ghost zones are reset according to the new updated variables with another call to sweepbc.f. This time the transverse velocities must also be included, as the next step will remap all fluid quantities to the old grid. Note that xa(nmax+1) is updated within ppm, and should not be reset here! remap.f does just what its name implies, it remaps the new fluid quantities from the new, evolved lagrangian grid given by the xa's, to the old, Eulerian grid, given by the xa0's. Again, remap is passed the BC and geometry flags.
Finally, the new 1D arrays are loaded back into the multidimensional arrays. sweepy.f is similar to sweepx.f except that it loops over a different index and uses the appropriate coordinate (e.g., xa = zya).
sweepbc.f
forces.f
ppm.f
The input states on the left side of the left boundary, alft(nmin), and the right side of the right boundary, argh(nmax+1), must be set according to the BC flags passed from the calling sweep routine. These input states are passed to riemann.f which calculates the time averaged value of the pressure and velocity at each zone face, from nmin to nmax + 1: umid and pmid. Note that the density is not needed lagrangian evolution. Currently we force umid = 0 on a grid boundary if it uses reflecting boundary conditions. Ideally this will be the case automatically, but round off errors can creep in.
Finally, evolve.f is called to update the fluid variables using the time-averaged quantities, umid and pmid.
flatten.f
The variable shock is unity if a shock is in the zone and zero if there is no shock, based on the pressure gradient and the sign of the velocity gradient (to distinguish from rarefaction waves). The steepness of the pressure gradient is stored in steep, and is set to zero if there is no shock in the zone (shock=0). The flatten parameter is set to the highest value of steep in that zone and the two neighboring zones, ie, a zone is flattened if a strong shock is in that zone or only one zone away. The parameters used in this routine are chosen based on extensive testing and the advice of CW and W.
paraset.f
The last four vecors are used in the more detailed parabolic interpolation used when dealing with coordinate singularites. This curvilinear version is considerably more complex and slower, and is only advantageous near coordinate singularities. For that reason we have left this additional coding commented out of the standard code. To reincorporate these geometry additions one must uncomment the lines beginning with CG, and follow any specific information written into the code and marked with CGEOMETRY. For more information on this alternative see S5 on curvilinear coordinates.
parabola.f
If iskip = 1 (only when calling with total energy from remap) then the values of al and deltaa are already known and we only need calculate a6 and apply the monotonicity constraints.
Otherwise calculate the ar's (and al(n + 1) = ar(n)) according to CW. In calculating the ar's we apply a different limitation than that given by eq. (1.8) in CW . We simply require ar(n) to be in the range of a(n) and a(n + 1). This is important in the vicinity of shocks where a steep pressure profile might lead to an overshoot in ar(n), giving a negative pressure just in front of the shock. While our constraint will not steepen sharp gradients as would eq. (1.8) of CW , it can be applied to a correct curvilinear interpolation scheme as described in chapter 5.
If iflat = 1, (when calling from {\bf ppm.f}, not from {\bf remap.f}) then use the passed array of flattening coefficients, flat, to smear out steep shocks.
Finally, apply the monotonicity contstraints. These are very important in providing a stable algorithm! We apply two different constraints in this loop (510). If the zone average value a is a local extremum, flatten the profile in that zone: ar(n) -> a(n). If the derived parabola introduces a new extremum, ie, if the parabola exceeds ar or al, then reset ar or al so that the slope -> 0 at the side closest to where the parabola exceeded the boundary value (see CW). The order of these operations in loop 510 is; (1) calculate the deltaa's and a6's, (2) compute a flag, scrch, that is negative if a is a local extremum, (3) flatten ar and al if scrch < 0, (4), compute the monotized values of al and ar for the second monotonize test, (5) compute a flag, scrch, to flag the left boundary for monotonizing, and reset al if scrch < 0, (6) compute scrch to flag the right boundary for monotonizing, and reset ar if scrch < 0, and finally, (7) recompute deltaa and a6. While it is a little cumbersome to pack all this into one loop, it provides for considerable speed up on vector machines such as the Cray.
states.f
First calculate the fraction of the zone to be averaged using the average sound speed in that zone times the timestep, divided by the width of the zone. A factor of 1/2 is included in the definition of Cdtdx to save multiplications later on. If the coordinate is angular, the zone width must be multiplied by the variable radius to get a physical dimension.
The states are then computed using the formulae in eq. (1.12) of CW, with the addition of force corrections to the velocity states. These corrections of delta tg/2 represent the time averaged acceleration of the fluid to first order only. At this point we also make sure that the pressure states remain positive by not letting them fall below smallp. If geometry corrections are required (currently implemented for radial cyclindrical and spherical), the velocity states are altered according to eq. (2.9) of CW . Note that other corrections may be desired at this point to account for other source terms, including forces and energy sinks and sources.
riemann.f
This particular subroutine can be replaced by other riemann solvers developed in the future that may be more accurate or more efficient.
The nonlinear equations to be solved at each zone interface
look like (CW, van Leer 1979):
Where it is assumed that the equation for the lagrangian wave speed of a shock
is a sufficiently valid approximation for a rarefaction wave so that only
one wave speed need be calculated.
To start off the iteration, the wave speeds W_L,W_R are approximated by
the adiabatic lagrangian sound speeds, i.e., C_L^2 = gamma p_L rho_L ,
and the
above equations are solved for p-bar. In smooth flow this is all that
is required. In the vicinity of discontinuities, we must iterate further
to converge on the correct value of p-bar.
First, we find the current guesses for u-bar based on the current guess
for p-bar and the left and right wave equations above. Let us call these
u_L^* and u_R^*. We then use the
slope of the left and right adiabat (Z_L, Z_R)
in the u - p plane at these two points
(p-bar,u_L^*) and (p-bar,u_R^*),
to extrapolate the tangents to the point of intersection,
from which we find our next guess for p-bar:
This new guess can then be used to start another iteration (compute wave
speeds, solve for p-bar) until the
desired accuracy is obtained. In practice we automatically iterate
5 times on all zones, although more efficient algorithms are certainly
possible.
Once a converged value for p-bar has been found, the value of
u-bar can be found by combining the above equations into an
expression for u-bar:
evolve.f
Next, the new zone volumes are computed according to the geometry flag passed from ppm.f. At the same time the average area of the zone face, amid, is computed to account for non-planar geometry in the hydro equations.
Evolution of the density is particularly trivial. Because the total mass in the zone remains constant in a Lagrangian calculation, the change in the zone-averaged density is inversely proportional to the change in zone volume. The new density is protected from going negative by setting a floor value of smallr.
The velocity is updated according to the pressure gradient and the forces, both fictitious and external. The total energy is updated according to the enthalpy gradient and the external forces only. The fictitious forces do not change the total energy, but rather redistribute the momentum between the different coordinate directions. If included in the energy equation, the change due to fictitous forces in one direction should exactly cancel the change found in the orthogonal direction. The pressure is found by subtracting the kinetic energy from the new total energy, and protected from falling below the floor value, smallp. If the problem at hand is in planar coordinates and no external forces are used, one can simply comment out the lines containing grav and fict in the velocity and energy loops, 40 and 50.
remap.f
The first order of business is to compute the new mass elements dmu and the width of the overlapping subshells, delta. Note that delta can be positive or negative, depending on the sign of the velocity at the zone interface. As in the hydro step, we call paraset to get the parabolic coefficients, paracoff, for the evolved Lagrangian grid. This time no flattening is desired, so we pass 0 for the flattening flag along with a dummy array in place of flat. Next we make individual calls to parabola to get the profiles of the variables to be remapped. The values of these variables (and the coordinates) in the ghost zones have already been set in the calling sweep routine. As discussed in CW , a ``fix'' must be implemented for the remap of the total energy to try and account for the fact that total energy contains the square of the velocity. A parabola fit to the total energy will not match the square of the parabola fit to the velocity, ie, although we have E = P/rho + 0.5u^2 for the zone averages, the same will not be true for each point on the respective parabolic interpolations. To compensate for this possible discrepancy, CW suggest using the left and right values of the velocity and internal energy to construct the total energy profile. In this way the different parabolas will at least satisfy the above energy equation at the zone boundaries, which is where the calculation of the subfluxes originate. (In the limit of small grid motion, we will only need the quantities very near the zone interfaces). The loop 12 does the construction of el and er as described above, and the subsequent call to parabola passes a flag to tell parabola that the left and right values are already determined.
Next, the ugly face of BCs raises its ugly head! If there is any motion of the grid faces at the boundaries of the grid, we need to know the variable profiles in the first ghost zones, nmin - 1 and nmax + 1. We set these values based on the boundary condition flags passed down to remap.f: nleft and nright. Note that we reset the boundaries even in the case of a reflecting boundary in order to avoid problems with uninitialized variables.
We remap the conserved quantities in two steps: First, we compute the integral of the conserved quantities in the overlapping subshell using the parabolic profiles calculated above, and second, we update the zone averages by adding (or subracting) the appropriate quantities in the subshells.
The total quantity of a given variable in a subshell, or ``flux'', is stored in an array; fluxr (total mass), fluxu, fluxv (velocities), and fluxe (total energy) equal to the integral of that variable over the volume of the subshell: If delta > 0 the flux is positive and equal to the integral from the right side of zone n - 1 back to the original zone boundary, xa0(n). If delta < 0 the flux is negative and equal to the integral from the left side of zone $n$ forward to the original zone boundary, xa0(n). The advection is then just a matter of redistributing these subshell fluxes to the appropriate zone:
The quantities to be remapped are (1) the total mass, fluxr, (2) the momentum, fluxr*fluxu, and (3) the total energy, fluxr*fluxe. We use the roundoff protection for the pressure and zone mass, and enforce the floor pressure, smallp.
2.2 THE REST OF THE CODE
back to Table of Contents
The main driver for VH-1 is vhone.f. This program reads the input, sets up the grid ( grid.f) and initializes it (init.f ), handles the output of data ( prin.f and images.f), constrains the timestep (dtcon.f ), calls the appropriate sweep routines, and other general managing tasks. In this section we will describe the driver and each of the peripheral subroutines.
vhone.f
The main loop consists of calling the appropriate sweep subroutines and dtcon, incrementing the loop variables, (e.g., time, ncycle), and calling the appropriate peripheral routines based on the values of nmovie ( images.f), nprin (prin.f), and ndump (dump.f). If the time or cycle number is not yet at the final value, loop back to label 1 and do it again. For 2D runs the cycles are split into even or odd values of ncycle so that the order of sweep calls is XYYX (to approach second order accuracy).
When the code has reached the maximum number of iterations (ncycle = ncyclend), the code closes the output files and exits. If the value of time at the next timestep will be greater than the maximum time stated in indat (time + dt > endtime), the value of dt is adjusted to end exactly at endtime and ncycp, ncycd are set to nprin and ndump so that a final data file and dump file will be written after the last timestep.
grid.f
init.f
dtcon.f
where
This can be a big advantage over Direct
Eulerian codes that must use the sum of the
fluid velocity and the sound speed, 
An
additional constraint is implemented that does not allow the timestep
to grow too fast, in this case limiting the new timestep to 1.1 times
the old timestep. If source terms are included in the calculation, they
may impose additional constraints on the timestep.
A value of courant of 0.7 or less is recommended to insure stability.
The 1D Sod shock tube test problem provides a good example of the error
introduced at a shock front when a larger value of courant is used.
prin.f
images.f
If non-cartesian computational grids are used (ngeom > 2), the data must be interpolated onto a cartesian grid for output by images.f. This is accomplished using the arrays of interpolation indices ipa and jpa created in init.f. This feature is automatically invoked when ngeom > 2, but the user must be careful to set ig, jg in zone.h to appropriate values before compiling the code (e.g., ig = 128, jg = 128). These arrays ipa, jpa store the number of the computational grid point closest to (but less than) the grid point on the cartesian raster grid. i.e., ipa and jpa are of dimensions equal to the dimension of the raster grid, and contain the grid coordinates from the computational grid that fall closest to the given point on the rastar grid. If a cartesian grid is used there is no need for declaring ipa and jpa so these can be declared length 1 in zone.h in order to save on memory allocation.
dump.f
2.3 RUNNING VH-1
back to Table of Contents
In order to run VH-1 properly, the user must first tailor the code to the specific proplem that he/she desires. Normally this involves (1) including the proper input data in the file indat, (2) setting up the initialization subroutine, init.f, and (3) possibly altering the output subroutines to write out the data in the user's desired format. This is the minimum list of required changes, and more difficult problems may involve more changes to the code for proper results. Such changes might involve additional body forces, switching to Lagrangian mode in 1D, adding new boundary conditions or incorporating additional source terms. Often the user will add parameters to global.h for use on a specific problem. In order to maintain the ability to restart the code from a dumpfile, the user must be sure to increase the size of BDUMP accordingly.The input deck containing the information needed to define a specific run is written in the file indat. The user is encouraged to adapt this file to suit his/her needs. The release version of VH-1 provides the following inputs in indat:
restrt suffix of dumpfile for restart, else if='no' then
start new problem.
prefx 5 character prefix for naming output files.
gam Adiabatic index gamma
imax Number of zones in the X direction
jmax Number of zones in the Y direction
xmin Physical location of left side of grid
xmax Physical location of right side of grid
ymin Physical location of bottom side of grid
ymax Physical location of top side of grid
courant Courant number, usually less than 0.8
ncycend Ending cycle number
ndump Number of cycles at which to dump a restart file
nprin Number of cycles at which to print out a data file
nmovie Number of cycles at which to write out a movie file
tprin Interval of time between writing out a data file
endtime Time at which to end the run
The number of zones in a given direction must be equal to or smaller than
the array sizes set in
zone.h for each direction (imax <= ii,
jmax <= jj),
and at least 4 zones smaller than the variable maxsweep set in
sweepsize.h,
which defines the length of the 1D sweep arrays (imax,jmax
< maxsweep + 4).
These extra 4 zones are to accomodate the 2 ghost zones on each side of
the 1D arrays used for boundary conditions.The initialization subroutine init.f contains several required sections. First, several flags are set that define the problem (geometry, boundary conditions, dimensions). Then the computational grid is initialized with the desired initial conditions. If inflow boundary conditions are used, the inflow values must be specified. If movie images are to be created, the normalization is set here. Finally, a guess for the initial timestep is made based on the initial data.
The geometry for each dimension is set in the data statements
data ngeomx / 0 /
data ngeomy / 0 /
where the value of ngeomx, ngeomy corresponds to:
0 planar
1 radial cylindrical
2 radial spherical
3 angular cylindrical
4 angular spherical polar
The boundary conditions for each dimension on the left and right
side are set in similar data statements, with the data values
corresponding to the following available boundary conditions:
0 reflecting
1 zero gradient in/outflow
2 fixed inflow (values must be set in init.f)
3 periodic
5 personalized boundary conditions
These flags are used in
sweepbc.f, ppm.f and
remap.f to implement
boundary conditions within the 1D hydro sweep. The last option, 5,
is left open for the user to include his own problem specific boundary
condition. To use this, the user must write his/her new boundary
condition into each of the three subroutines requiring knowledge of
the boundary conditions:
sweepbc.f, ppm.f and
remap.f.
Occasionally more complicated boundary conditions will be needed
that are dependent upon the position within the computational grid.
This problem arises in two of the test problems provided with VH-1:
the oblique shock and the Mach 3 wind tunnel with step. In this case,
sweepx.f and
sweepy.f are altered to set the value of
the BC flag as a function of grid location.The dimension of the problem, either 1 or 2, is set with ndim.
The entire computational grid is then initialized according to the problem at hand. An initial guess at a timestep is also required, but this will be changed by a call to dtcon.f before the first timestep. Because dtcon.f will not increase dt too fast, it is best to give an overestimate of dt and let dtcon.f reduce it according to the courant condition. If fixed inflow boundaries are used in the problem, they must also be set in init.f (or possibly set for each sweep in the appropriate sweep subroutine if they are different for each dimension or row). The details of the initialization can be written to the history file as a record of the run.
If movie images (e.g., of density) are to be output from
images.f,
the value of cmin and cmax must be set. These are them minimum
and maximum values of the imaged variable expected in the simulation.
If the range in values is large, you may want to consider outputing
the logarithm of the variable, in which case cmin and cmax must
be logarithms and you must add a line taking the log of the imaged
variable to the subroutine
images.f, i.e., 
If the problem is not run on a cartesian grid, and you wish to output
movie files on a cartesian grid, it saves time to do the interpolation
in init.f and store the interpolation
values in memory. (see
subsection on init.f and
images.f.)
It is also likely that the output subroutines will have to be modified for the desired problem. If a format other than HDF is desired, the output routines will have to be changed. The user is encouraged to adapt these routines to his/her specific needs.
3. TEST PROBLEMS
back to Table of Contents
Several test problems are included with the distributed version of VH-1. These problems serve two specific purposes. First, they provide the new user with immediate, working examples with which he/she can begin to learn how to operate the code. It is highly recommended that the new user begin by running all of the enclosed test problems to ensure that the current code is working properly for a variety of boundary conditions and geometries. The second purpose of these sample problems is to provide a series of tests that can be rerun each time the code is altered. In addition to the problems presented here, it may prove useful to add problems that are similar to the types of flow problems the user intends to work on with VH-1. In general, the closer the test problem is to the actual working problem, the less chance for errors to arise through ``insignificant" changes in the code.To run these test problems one need only copy the relevant files from the test suite directory into a directory with a working copy of VH-1. Normally this just involves the init.f subroutine and the indat input file, but more complicated problems will also require copying over new sweepx.f and sweepy.f subroutines. The prefix of these test subroutines must be removed (e.g., ms_init.f -> init.f) and the entire code compiled. The size of the declared arrays may also have to be changed prior to compiling the code. As a caution, all of the *.o files should be removed by typing ``make clobber" and the code recompiled anytime any of the common blocks are altered. To run the test problem one then types the magic word, ``vhone".
Output files from these test problems are provided in the Output directory. Your test results should look identical to the corresponding images stored in /Output. All of these runs were run at relatively low resolution so that the cpu time needed to perform these tests would remain minimal.
3.1 SOD SHOCK TUBE
s_init.fback to Table of Contents
The Sod shock tube has become a standard test problem in computational hydrodynamics. The initial conditions are very simple, a contact discontinuity separating gasses with different pressures and densities. In the standard case the density and pressure on the left are unity, and the density on the right side of the contact is 0.125 and the pressure is 0.1. As the evolution begins, a shock propogates to the right while a rarefaction wave travels to the left. In the standard case these waves are relatively weak, and most hydro codes produce good results. A more demanding test is imposed if the initial density and pressure ratios are increased by an order of magnitude!Figure 2 shows the flow variables for the standard Sod shock tube problem at t=23.0, along with the analytic solution drawn as a solid line. The shock and contact are both very sharp.
This problem can be expanded to multiple dimensions by placing the discontinuity across a 2D or 3D grid. We have chosen to position the contact at a 45-degree angle so that the shock and rarefaction wave propagate diagonally across the grid. Figure 3 shows the density of this 2D shock tube problem where the initial high density gas was in the lower left corner. The image is very close to symmetric about the diagonal, but some small differences can be seen. The image is at time = 0.871. The shock wave has reflected off of the far walls and is converging back towards the orginal corner. The contact discontinuity is beginning to ripple after the recent passage of the reflected shock. This is the Richtmyer-Meshkov instablity.
3.2 DOUBLE MACH REFLECTION
o_init.fback to Table of Contents
This problem is taken from Woodward and Colella (1984), and is used extensively as a test problem because of the complicated structures generated and the availability of experimental solutions to compare with numerical results. In this case problem the shock is set at an angle of 30-degrees and is allowed to run into the lower reflecting boundary. The shock is anchored at the lower left corner by applying outflow boundary conditions for the first few zones on the lower boundary (see Woodward and Colella 1984 for details). Figure 4 shows the density structure at t = 0.2. This simulation can be compared with Figures (9) of Woodward and Colella.
3.3 MACH 3 WIND TUNNEL WITH A STEP
m_init.fback to Table of Contents
This test problem is again pulled from Woodward and Colella (1984). In this case a uniform grid is initialized with planar supersonic flow from the left. A step is put into the flow downstream from the left boundary by adjusting the location of the reflecting boundary with the sweepx.f and sweepy.f subroutines. In Figure 5 the simulation is stopped at t = 4.0. This simulation can also be compared with the results of Woodward and Colella (see their figures 7). We have not implemented the corrections to the grid zones around the step corner as described in Woodward and Colella. This affects the flow solution along the surface of the step, and as such this region is different from the solutions published in Woodward and Colella. We have run this problem with such corrections and obtained good results, but have not bothered to incorporate the added complexity into the distributed test version. The contact discontinuity along the top of the flow (y ~ 0.8) is Kelvin-Helmholtz unstable, and slight wiggles can be seen even in this low resolution run. The seed perturbations for this instability are attributed to postshock noise (see Chap. 4).
3.4 SUPERSONIC FLOW PAST A CYLINDER/SPHERE
b_init.fback to Table of Contents
We have included a generic test problem of a plane parallel flow past a rigid sphere or cylinder (for ngeom = 2 or 1, respectively). This test problem verifies the geometry dependent terms and fictitious fources by evolving a plane parallel flow on a curvilinear grid. The inner radial boundary condition is reflecting (the solid sphere/cylinder), and the outer radial boundary is set for free inflow/outflow. The angular boundaries (theta = 0, pi) are set for reflecting conditions. The grid is initialized with a uniform flow at Mach number umach set in init.f. The resulting simulations compare favorably with shadow-graph images of similar experiments in fluid flow. In particular, the presence of the slip line and a shock reflecting off the axis behind the sphere agree well with supersonic flow experiments. Figure 6 shows a Mach 3 flow past a sphere. Similar results are obtained for flow past a cylinder.
4. DISSIPATION
back to Table of Contents
PPM is a very nice method for compressible hydrodynamics, but it does have an Achilles heel. When strong shocks propagate slowly across the numerical grid low frequency noise is generated in the post shock region. This noise is a consequence of the narrow shock structure generated by PPM. Because the shock is resolved into only one or two zones, the exact profile of the shock is dependent on the relative location of this narrow profile with respect to the zone faces. When the shock front crosses a zone face the profile necessarily changes. If the shock is moving relatively fast with respect to the grid, these changes result in high frequency noise that is quickly damped. Slowly moving shocks on the other hand produce low frequency noise that the highly accurate PPM algorithm is able to accurately follow without much dissipation. By page count, CW spend 1/3 of their paper discussing dissipation mechanisms to try and minimize this problem!This problem can be seen in the creeping shock test problem (see Chapter 3), where a strong (Mach 10) shock is propagating at 1% the shock speed with respect to the numerical grid. Intuitively one could imaging trying to suppress this instability by (a) broadening the shock structure so that it does not change appreciably while crossing the grid, (b) increase the relative velocity of the shock so that the generated noise is of sufficiently high frequency that it is effectively damped, or (c) add some sort of artificial dissipation at shock fronts to dampen any high frequency component. These are in fact the three options suggested by CW, flattening, grid wiggling, and artificial viscosity. The first of these options, shock flattening, has proven simple, robust, and unobtrusive. The second of these, grid wiggling, is found to be beneficial in some circumstances, but can be rather cumbersome. The last option, artificial viscosity, has not been shown to work with VH-1, despite the fact that it is beneficial in many other versions of PPM (this remains a mystery). We will discuss these three techniques below, but we have only implemented the first of these options in VH-1.
4.1 SHOCK FLATTENING
back to Table of Contents
This technique involves a reshaping of the interpolated parabolas in the vicinity of shock waves so that the shock is effectively broadened over a couple zones. This is accomplished in two steps, (1) the calculation of a flattening parameter based on the presence and strength of a shock (for VH-1 this is done in flatten.f), and (2) the mixing of the interpolated parabola and the zone average (for VH-1 this is done in parabola.f). In the extreme limit of everything being flattened, this method would transform to the original Godunov scheme, using only zone averages to compute the input states for the Riemann problem. In Figure 7 we also show the same problem computed with flat(n) = 1 in both the hydro and in the remap step, i.e., a Godunov scheme. The postshock oscillations are indeed removed, but the shock is now spread out over many zones (cf. CW who find some noise remains with complete flattening). CW describe several stages of complexity in computing the flattening coefficient. We have found that the more complex methods do not show significantly better results in our slowly moving strong shock problem than the simplest method. We have therefore only included this simple approach. Figure 7 also shows the shock structure when flattening is applied as it is implemented in the release version of VH-1. While there is still some postshock noise, it is of substantially smaller amplitude than in the unflattened solution.One variation that we have not included is a sensing of the presence of a shock transverse to the current direction of updating. This could concievably make a difference in cases where a strong shock is sitting nearly parallel (but not exactly) to the numerical grid. In this case the postshock structure is flattened perpendicular to the shock front but not along the direction of the shock front. Instituting a multidimensional flattening routine is much more cumbersome than in 1D and we have not found the results sufficiently encouraging to warrant such an addition to the standard code. In these difficult cases we find that grid wiggling provides a more effective suppression of postshock noise.
4.2 GRID WIGGLING
back to Table of Contents
The easiest way to avoid this whole problem is to only work on problems where the shocks are rapidly moving across the numerical grid. Unfortunately it is often the case that the interesting thing one is studying is associated with the shock and it is desireable to have the shock in one location on the grid throughout the simulation. If you can't move the shock with respect to the grid, move the grid with respect to the shock! This is the essence of the grid wiggling technique described in CW. At each timestep one can move the grid back and forth, effectively smearing out the steep shock profile over more than one zone. One can imagine going about this in several ways:(1) As in CW, wiggle the grid at each zone location according to the amount of dissipation needed and then wiggle back to the original grid on the next timestep.
(2) Wiggle the entire grid a fixed fraction of a zone each timestep, alternating directions of the wiggle each time.
(3) Wiggle the grid on the remap step and then re-remap the variables back onto the original grid at the same timestep.
In more than one dimension the grid must be orthogonal in between timesteps, so that if option (1) were used there must be two X sweeps (wiggle forward, wiggle backward) followed by two Y sweeps. In practice this is done anyway to minimize cross-term errors in the directional splitting. In practice we have found this technique to be most effective when there is a strong shock almost parallel to the zone gridding. In this case small changes in the shock/zone orientation from row to row can produce large differences in the postshock entropy, generating substantial postshock noise along the inside surface of the shock. If the grid is wiggled along the shock front (rather than across it), these postshock variations from zone to zone are effectively smoothed out.
Option (3) requires an extra remap step and hence is less efficient, but we include it here because there may be circumstances where it is the only feasable option.
While using a fixed length wiggle will add dissipation where it is not needed (away from shocks), this added dissipation is not significantly more than what is already present in the algorithm (which is very little). We have not found any significant degradation when a fixed length wiggle is used rather than a selective wiggle as described in CW, although we have not performed any quantitative testing of these different possibilities.
4.3 ARTIFICAL VISCOSITY
back to Table of Contents
Artificial viscosity has been used effectively to dampen short wavelength oscillations in several versions of PPM (e.g., codes written by Fryxell, Balsara, Stone). We have not been able to formulate an artificial viscosity term that improves the performance of VH-1, and so have not included it with this code. The reason why VH-1 behaves in this manner compared with other versions of PPMLR is not clear.
5. CURVILINEAR GEOMETRY
back to Table of Contents
Curvilinear geometry is treated on two different levels in VH-1. The first level of complexity is just computing the appropriate volume (dvol(n)) and area (amid(n)) elements and including the fictitious forces (centrifugal and coriolis) on angular grids. This is sufficient for most applications, but is prone to large numerical errors near coordinate singularities. A correct derivation of PPMLR in curvilinear geometry is presented in Blondin and Lufkin (1992). The modifications to the code required by this technique are included in the released version of VH-1 but have been commented out because they are needed only for specific problems. If the user encounters a problem requiring accurate solutions near a coordinate singularity, the necessary coding can be uncommented in paraset.f and parabola.f by following the directions in those subroutines.
GLOSSARY OF VARIABLES USED IN VH-1
A to M
M to R
R to Z
back to Table of Contents
The following is an alphabetical listing of all variables used in VH-1. The location of the declaration of each variable is listed. This will be either one of the include blocks (global.h, zone.h, or sweep.h), a specific subroutine, or local if the same variable is used locally in more than one subroutine. All of the variable arrays correspond to the 1D sweep array used in the individual sweeps, except for the arrays beginning with z, which are the 2D arrays of the fluid variables and coordinates.
VAR LOC DESCRIPTION
a parabola.f generic fluid variable p,rho, u, etc. adump zone.h vector to hold as much memory as needed to for restart al parabola.f value of interpolation at left zone face (returned as pl, rl, etc.) almon parabola.f monotonized value of interpolation at left face amid evolve.f time averaged area of left zone face ar parabola.f value of interpolation at right zone face armon parabola.f monotonized value of interpolation at right face avgmass global.h average mass per particle (not needed for adiabatic hydro) a6 parabola.f parabolic coefficient for interpolated variable bdump global.h vector to hold as much memory as needed to for restart (holds /global/) boltzman global.h Boltzmann's constant Cdtdx states.f fraction of the zone contained by the domain of dependence cleft riemann.f lagrangian sound speed in the left input state courant global.h courant parameter crght riemann.f lagrangian sound speed in the right input state da parabola.f functional variable used to calculate interpolation delta remap.f dvol - dvol0 deltaa parabola.f linear coefficient for interpolated variable diffa parabola.f zone difference of variable, diffa = a(n+1) - a(n) dinflo global.h value of density at inflow boundary (used if nleft, nright = 2) dmpfile vhone.f name of restart dump file dmu local lagrangian mass element before remap (used in evolve, remap) dmu0 remap.f lagrangian mass element after remap dp local linear coefficient of interpolated pressure (calculated in parabola) dt global.h timestep dr local linear coefficient of interpolated density du local linear coefficient of interpolated velocity dvol sweep.h zone volume dvol0 sweep.h dvol on fixed, Eulerian grid dvol1 evolve.f dvol, prior to lagrangian update dx sweep.h xa(j+1) - xa(j) dx0 sweep.h dx on the fixed, Eulerian grid dxfac grid.f fixed size of dx for making grid coordinates einflo global.h value of internal energy at inflow boundary enrg sweep.h total energy: zone average fCdtdx states.f 1. - Cdtdx fict local fictitious forces (i.e., centrifugal, coriolis) fict1 evolve.f fictitious forces (i.e., centrifugal, coriolis) at t+dt filename prin.f output file for HDF scientific data sets flat local array of flattening coefficients calculated in flatten fluxe remap.f total (total) energy in overlapping subshell fluxr remap.f total density in overlapping subshell fluxu remap.f total velocity (direct) in overlapping subshell fluxv remap.f total velocity (tangential) in overlapping subshell gam global.h equation of state gamma gamm global.h gam - 1 gamma riemann.f =gam grav local gravitational (or other) force grav0 evolve.f gravitational (or other) force at t grav1 evolve.f gravitational (or other) force at t+dt hstfile vhone.f name of run history file i local zone number in X direction iflat local if iflat = 1, use flattening ig zone.h number of zones in image array (ig X jg) ii zone.h number of zones in x direction imax zone.h number of physical zones in the x direction ipa zone.h 2D array of i coords for nearest neighbor on image grid j local zone number in Y direction jg zone.h number of zones in image array (ig X jg) jj zone.h number of zones in y direction jmax zone.h number of physical zones in the y direction jpa zone.h 2D array of j coords for nearest neighbor on image grid maxsweep sweep.h maximum size of any sweep set in sweepsize.h mvyfile vhone.f name of file containing HDF raster images n local zone number in 1D sweep array (n=i+2 or j+2) ncycend vhone.f maximum number of cycles to run program ncycd vhone.f number of cycles from last restart dump ncycle global.h cycle number ncycm vhone.f number of cycles from last imagefile write ncycp vhone.f number of cycles from last dataset write ndim global.h number of spatial dimensions in problem ndump vhone.f number of cycles between writing out restart dump nfile global.h incremental number used as suffix for output filename ngeom local geometry flag for current sweep direction ngeomx global.h geometry flag for X direction ngeomy global.h geometry flag for Y direction nleft local BC flag for left side of current sweep direction nleftx global.h BC flag for left side of X direction nlefty global.h BC flag for right side of Y direction nmax sweep.h last real zone in 1D arrays nmin sweep.h first real zone in 1D arrays nmovie vhone.f number of cycles between writing out imagefile nprin vhone.f number of cycles between writing out dataset nright local BC flag for right side of current sweep direction nrightx global.h BC flag for right side of X direction nrighty global.h BC flag for right side of Y direction ntotal zone.h total length of stuff to be stored in restart file nzones grid.f number of physical zones in current direction p sweep.h pressure: zone average paracoff local array of parabola interpolation coefficients pi global.h pi = 3.141592654 pinflo global.h value of pressure at inflow boundary pl local interpolated pressure at left zone face plft local average pressure in the + wave (left side of interface) pmid local time averaged pressure at left zone face prgh local average pressure in the - wave (right side of interface) prefx vhone.f prefix for naming output data files p6 local parabola coefficient for interpolated pressure radius sweep.h radius of current column (for curvilinear geometry) rho sweep.h density: zone average rndoff global.h magnitude of limited roundoff error rl local interpolated density at left zone face rlft local average density in the + wave (left side of interface) rrgh local average density in the - wave (right side of interface) r6 local parabola coefficient for interpolated density scrch* local scratch arrays shock flatten.f flag presence of a shock, = 1 if shock, else = 0 smallp global.h floor value of pressure smallr global.h floor value of density steep flatten.f array of steepness parameter for calculating flat sweep sweep.h direction of current sweep, = 'x' or 'y' time global.h time timep vhone.f time from last dataset write tprin vhone.f time interval between writing out dataset u sweep.h zone average velocity uinflo global.h value of velocity at inflow boundary ul local interpolated velocity at left zone interface ulft local average velocity in the + wave (left side of interface) umid local time averaged velocity at left zone face uold evolve.f velocity in zone before lagrangian update upmid evolve.f umid * pmid \cr urgh local average velocity in the - wave (right side of interface) u6 local parabola coefficient for interpolated velocity v sweep.h zone average transverse velocity vinflo global.h value of transverse velocity at inflow boundary wlft riemann.f wave speed in left state of Riemann shock tube problem wrgh riemann.f wave speed in right state of Riemann shock tube problem xa sweep.h grid zone face locations for 1D sweep array xa0 sweep.h grid zone face locations on fixed 1D sweep Eulerian grid xa1 evolve.f grid zone face locations prior to lagrangian update xa2 evolve.f zone center coordinate prior to lagrangian update xa3 evolve.f zone center coordinate after lagrangian update xmax vhone.f physical location of end of X coordinate xmin vhone.f physical location of beginning of X coordinate xsweep local sweep flag: xsweep = 1 if sweeping in X direction, else = 0 ymax vhone.f physical location of end of Y coordinate ymin vhone.f physical location of beginning of Y coordinate z*** arrays that begin with z correspond to the fixed, 2D grid zrho zone.h density: zone average zp zone.h pressure: zone average zux zone.h zone average velocity x direction zuy zone.h zone average velocity y direction zxa zone.h x grid zone face locations zya zone.h y grid zone face locations zdx zone.h xa(i+1) - xa(i) zdy zone.h ya(i+1) - ya(i)