VH-1: A beginner's guide

This web page is intended to serve as a tutorial to those who are unfamiliar with VH-1: where to get it, how to get it running, how to set up a few basic problems, and what the next steps are in using VH-1 in your research. (Note: This tutorial presumes you're working in a Unix/Linux environment, or at least have all the common commands available to you.)

VH-1 is a multidimensional ideal compressible hydrodynamics code written in FORTRAN, for use on desktop workstations or supercomputers. It uses a Lagrangian remap version of the Piecewise Parabolic Method developed by Paul Woodward and Phil Colella in their 1984 paper. VH-1 comes in a variety of versions, from simple one-dimensional serial variant to multi-dimensional versions scalable to thousands of processors; this tutorial will use its most basic package, VH-0.5. This version has had most of its functionality stripped away. It can't run simulations of more than one dimension, and many of the subroutines from the full version have been streamlined or eliminated altogether. However, for our purposes, it's perfect: a simple version of the code with the key features laid out clearly and linearly.

Getting VH-1 quick-started

The first step is downloading VH-0.5 to your computer. The latest public version will always be available from John Blondin's research webpage. Follow the link at the top for "VH-1", and download VH-0.5 from the "Basic 1D code" link.

On your computer, create a clean directory with the command mkdir your_directory_name. Use either mv or cp to put VH0.5.tar into the directory. Unpack the tarball with the command tar -xvf VH0.5.tar.

If all went well you should see the following files/directories in your directory (along with VH0.5.tar):

• indat
• Makefile
• PPMLR, a directory
• vhone.f90

To compile VH-0.5, type make and the Makefile will build VH-0.5. Compiling will place create several new files in the directory: you will see some object files (*.o), a few modules (*.mod), and an executable binary named vhone. Important note: as released, Makefile uses the gfortran compiler; if your computer is using a different compiler, you should change the "F90" and "LDR" flags in Makefile to reflect the compiler you use.

VH-0.5 comes preset to evolve the Sod shock tube problem. Now that VH-1 has been compiled, open indat in your text editor of choice and verify two parameters to make sure a test run will give appropriate output: endtime should be set to 5.0e-01, and tprin should be 5.0e-03. This will result 100 output files over the course of the run, but offers sufficient temporal resolution to discern key features of the problem.

With indat adjusted, we can run VH-1. Type vhone (or ./vhone if required) to run the simulation and evolve the Sod shock tube. On modern computers VH-1 should take on the order of seconds to finish this particular computation. Once VH-1 finishes, typing ls will show you that VH-1 did indeed generate 100 output files as part of the simulation: DATAa1000.dat to DATAa1099.dat. Each of these files is a tab-delimited text file containing grid values for location, density, pressure, and x-velocity (in column order from left to right). There is also a history file, DATAahst, which contains details about the run.

Since the .dat files are formatted text, they may be plotted in whichever package you prefer. A program that will let you step through time and watch the evolution of the fluid is best, as this will allow you to track the various waves within the tube. At this point, you have a functional copy of VH-1 on your machine, and may proceed through the rest of this tutorial.

A look inside VH-0.5

The PPMLR directory can be safely ignored for the duration of this tutorial. It contains the guts of the hydrodynamics code, and should be left alone unless the problem you're working on requires changes to the subroutines there, which none of the sample problems need.

Makefile, if you haven't seen one before, contains instructions for the FORTRAN compiler on how to compile and link the subroutines of VH-1. For the moment, the contents of the Makefile may also be ignored, unless your system does not have the gfortran compiler installed; in that case you need to change the "F90" and "LDR" tags to whichever compiler you use. There are three commands associated with the Makefile that you should remember: make will compile VH-1 and output an executable that can be run; make clean removes the object files from a compiled VH-1; and make clobber removes object files, modules, and the binary executable itself. When you edit files between runs, typically make will see the modification and recompile that subroutine or module (and any others that depended on it). However, careful users may wish to use the make clobber command after modifications, just to be certain.

The indat file contains information that VH-1 will read in at the start of each run, like what prefix to use for output, how many cycles/how much simulation time to run for, or how often to output data. A full description of all the parameters in indat is beyond the scope of this tutorial, so instead the suggested changes will simply be given for each problem.

vhone.f90 is the key file of this package. It's also the file you'll be editing most frequently. In vhone.f90 you'll find the code that determines the size/spacing of the computational grid, that initializes the grid to the proper values at the start of the simulation, that evolves the simulation forward in time, and that controls the output. Instead of going line-by-line through the code, we'll only refer to segments as needed for the problem; the inline comments should be sufficient to understand how the rest of the code works.

Now that we've covered what you get in the VH-0.5 package, we'll discuss how to use it. The workflow for the vast majority of problems you'll face when using VH-1 can be distilled to four steps:

1. Determine the initial values your simulation should have, and set them in the code.
2. Edit vhone.f90 or other files to make specific adjustments for the problem you're considering.
3. Edit indat as needed to make sure your simulation evolves as long as you want it to.
4. make VH-1, and run the executable output.

The first through third steps will be different for any given problem. While you won't need to worry about the second step for any of the test cases below, we will cover the first and third steps in detail for each of the five problems in this tutorial.

Basic problems and their solutions

Now comes the meat of the tutorial: five basic hydrodynamics problems, and discussions of how to handle them using VH-1. Every problem addressed here will be run on a grid of 200 zones (i.e. imax = 200), with xmin set to 0.0 and xmax set to 1.0. We will discuss geometry, boundary conditions, and initialization within each problem.

Sod Shock Tube

The Sod shock tube is a standard test problem in computational hydrodynamics. This is because 1) it is very simple to set up, and 2) analytical solutions are known (see Sod 1978). The typical initial conditions for the Sod shock tube are to divide your grid into two halves, set the density and pressure on the left half to unity, and set the density and pressure on the right half to 0.125 and 0.1 (respectively).

Set your geometry and boundary flags to 0, if they aren't already. This tells VH-1 that our grid is Cartesian (as opposed to some variety of polar), and that both edges of the grid are reflecting. You'll need to declare four variables at the top of vhone.f90: the pressures and densities for both halves of the grid. You could choose any name you want for those variables, but we suggest pleft, pright, dleft, and dright. In addition to being easy to remember, they are also easy for others to interpret; it's good programming practice. After the comment "Set up parameters from the problem" you may wish to add the phrase "(Sod shock tube)" in case you come back to VH-1 after a long break and can't remember where you left off. If the ratio of specific heats gam isn't set to 5/3 (or 5.0 / 3.0 in fortran) you should do it here. Also, specify values for the four variables you declared earlier.

It's now time to initialize the grid. There are five fluid variables that VH-1 uses in its hydrodynamics calculations. In VH-0.5, they are p(n), the pressure; r(n), the density (from the Greek \rho  ); u(n), the component of fluid velocity along our grid dimension; and v(n) and w(n), the transverse components of the fluid velocity. We're running a one-dimensional simulation, so VH-1 knows not to update the transverse velocities, but they must be set and carried around anyway. Each one of the five fluid variables must be set in every zone when VH-1 initializes, or you will get some very strange results.

The problem requires two regions, and there are two ways to go about locating the midpoint to separate them. We could either look at the zone number (n = imax/2 + 6), or we could use the physical position (xa0 = 0.5). Either one would work, and we arbitrarily choose the latter option. The grid initialization takes place as a loop over all the zones in the grid; here it's a single one-dimensional loop, though more complex problems require nested loops for multiple dimensions. So within the loop over n, check to see if the zone in question is before xa0 = 0.5. If it is, set the pressure and density to pleft and dleft; otherwise use pright and dright. (Watchful readers might note that this avoids the issue of a zone falling exactly on the xa0 = 0.5 boundary. Instead of needing to average over the contact interface, we just displace the interface one grid zone to the left.) All three components of fluid velocities should be zero at the start.

The box below contains the code we used to initialize our Sod shock tube problem. If you get results that diverge greatly from the images and movies included at the end of this section, you can compare your code against ours for differences. Note that we included a description of the initial conditions in the comments for that section; while this isn't necessary, it makes your code more intelligible to others.

! Set up parameters for the problem (Sod shock tube)

pright = 0.1
dright = 0.125
pleft  = 1.0
dleft  = 1.0

! We also need a ratio of specific heats, gamma...
gam  = 5.0 / 3.0
gamm = gam - 1.0

! initialize grid for Sod shock problem:
!  left half of grid has density, pressure given by dleft, pleft
!  right half of grid has density, pressure given by dright, pright
!  velocity is zero everywhere

do n = nmin, nmax
 if (xa0(n) > 0.5) then
   r(n) = dright
   p(n) = pright
 else
   r(n) = dleft
   p(n) = pleft
 endif
 u(n) = 0.0
 v(n) = 0.0           ! note that we have to carry around the transverse
 w(n) = 0.0           ! velocities even though this is a 1D code
enddo

Having completed the initialization, it's important (read: critical) to keep records of what was done, and when. Virtually (probably absolutely) everyone in academia has forgotten where a key data set went, or what a certain figure is supposed to show. VH-1 outputs a history file for every run it does, and it is highly recommended that these history files be complete and kept with the data from the runs in question. In the history file section of vhone.f90 you can add whatever you feel is necessary to explain the run. This should include the version of VH-1 being used, the grid size and geometry, and a summary of the initial conditions/key parameters. For the Sod shock tube, this summary need only refer to the ratios in pressures and densities, since those are what drive the dynamics.

With vhone.f90 set up properly, it's time to edit indat. Though in your own work this step will involve a lot of intuition and guesswork, in the tutorial we'll make things easier for you. Edit the value of prefix to whatever you want; VH-1 is coded to use 5-character prefixes, so make sure yours is 5 characters long. Change the values of endtime, tprin and tmovie to be much larger than 1.0; everything we'll be watching for will be finished by a simulation time of 1, but it's better to include a safety margin of several orders of magnitude. The way VH-1 calculates timesteps, the denser your grid is the more steps it takes to reach a particular time; with an imax of 200, setting ncycend to 400 will capture most of the dynamics. The value of nprin should be below 10; if it is set higher we won't have the temporal resolution to see a few important effects. Don't worry about the nmovie variable, as it isn't used in VH-0.5. Included below is the indat used to generate the sample images and movies for the Sod shock tube.

&hinput
 rstrt   = 'no'
 prefix  = 'Sod1d'
 ncycend = 400
 ndump   = 5000000
 nprin   = 4
 nmovie  = 100
 endtime = 1.8e+11
 tprin   = 2.0e+12
 tmovie  = 5.0e+10 /

If you've set up the grid correctly, you'll see two major things happen as you evolve the simulation. A shock will propagate to the right of the initial contact discontinuity, and a rarefaction wave will travel to the left. STILL TO TALK ABOUT: STRONG STANDING SHOCK AT END, SQUIGGLES AT BOUNDARIES, CHANGING RESOLUTION TO LOOK AT EFFECTS (ESP. SHOCK SMEARING AND COMPUTATION TIME)

Sedov-Taylor Blast Wave

! Set up parameters from the problem (Sedov-Taylor blast wave)
pamb = 1.0
damb = 1.0
pstb = 1.0e9

! We also need a ratio of specific heats, gamma...
gam  = 5.0 / 3.0
gamm = gam - 1.0

! initialize grid for Sedov-Taylor blast wave:
!  grid everywhere has low density, pressure
!  first two zones of grid have extremely high pressure
!  velocity is zero everywhere

do n = nmin, nmax
  if (n < 3) then
    p(n) = pstb
  else
    p(n) = pamb
  endif
  r(n) = damb
  u(n) = 0.0
  v(n) = 0.0           ! note that we have to carry around the transverse
  w(n) = 0.0           ! velocities even though this is a 1D code
enddo

Colliding Blast Waves

! Set up parameters from the problem (colliding blast waves)
pleft  = 1000.
pmid   = 0.01
pright = 100.
damb   = 1.0

! We also need a ratio of specific heats, gamma...
gam  = 5.0 / 3.0
gamm = gam - 1.0

! initialize grid for colliding blast wave problem:
!  density everywhere unity, initial velocities zero
!  left tenth of grid is very high pressure
!  right tenth of grid is high pressure
!  middle of grid is very low pressure

do n = nmin, nmax
  if (xa0(n) < 0.1) then
    p(n) = pleft
  elseif (xa0(n) > 0.9) then
    p(n) = pright
  else
    p(n) = pmid
  endif
    r(n) = damb
  endif
  u(n) = 0.0
  v(n) = 0.0           ! note that we have to carry around the transverse
  w(n) = 0.0           ! velocities even though this is a 1D code
enddo

Strong Standing Shock

! Set up parameters from the problem (strong standing shock)
mach = 10.  ! mach number of incoming shock
pamb = 1.
damb = 1.

! We also need a ratio of specific heats, gamma...
gam  = 5.0 / 3.0
gamm = gam - 1.0

! More parameters from the problem, placed here because they depend on gam
ishock = (nmax-nmin)/2
pratio = (2*gam*(mach**2) - (gam-1.)) / (gam+1.)
pshock = pamb * pratio
dratio = (gam+1.)*(mach**2) / (2. + (mach**2)*(gam-1.))
dshock = damb * dratio
ushock = mach * sqrt(gam*pamb/damb) * (1. - 1./dratio)
upwind = mach * sqrt(gam*pamb/damb) * 0.9

! set inflow/outflow parameters for boundary conditions
dinflo = dshock
uinflo = ushock - upwind
vinflo = 0.
winflo = 0.
pinflo = pshock
!--
dotflo = damb
uotflo = -upwind
votflo = 0.
wotflo = 0.
potflo = pamb

! initialize grid for strong standing shock problem:
!  shock front at center of grid
!  left half is downstream (shocked material), right half is upstream (unshocked)
!  entire grid in moving frame with velocity 'upwind' to the left

do n = nmin, ishock
  r(n) = dshock
  p(n) = pshock
  u(n) = ushock - upwind
  v(n) = 0.
  w(n) = 0.
enddo
  r(ishock) = dshock * 0.5
  p(ishock) = pshock * 0.2
  u(ishock) = ushock * 0.5 - upwind
  v(ishock) = 0.
  w(ishock) = 0.
do n = ishock+1, nmax
  r(n) = damb
  p(n) = pamb
  u(n) = -upwind
  v(n) = 0.
  w(n) = 0.
enddo

Stellar Wind Bubble

! Set up parameters from the problem (stellar wind bubble)
pamb = 1.0
damb = 1.0
vwnd = 20.0

! We also need a ratio of specific heats, gamma...
gam  = 5.0 / 3.0
gamm = gam - 1.0

! set inflow parameters for boundary conditions
dinflo = damb
pinflo = pamb
uinflo = vwnd
vinflo = 0.0
winflo = 0.0

! initialize grid for stellar wind bubble problem:
!  supersonic wind sweeping into grid from left side (note that
!   c_s = sqrt(gam*pressure/density) = sqrt(5/3), and v_wind is 20;
!   our wind is indeed supersonic, with a mach number of ~15)
!  density, pressure held constant
!  velocity is zero everywhere

do n = nmin, nmax
  r(n) = damb
  p(n) = pamb
  u(n) = 0.0
  v(n) = 0.0           ! note that we have to carry around the transverse
  w(n) = 0.0           ! velocities even though this is a 1D code
enddo