VH-1: A beginner's guide

This web page is intended to serve as a tutorial for using the Beginner's version of VH-1. If you are unfamiliar with VH-1 this is the place to start. Here you will find answers to 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.)

To use the Beginner's version of VH-1 you will only need a few basic tools common to any unix operating system. Here we list common open-source solutions available on most linux systems.

F90 COMPILER

gfortran is available for most unix systems. Any FORTRAN 90 compiler will suffice. Commercial compilers exist for virtually every computing platform.

TEXT EDITOR

gedit is a common open-source choice, but many options exist. Do not use any editor that inserts unseen formatting tags into the file. (Don't use WORD!) Many commercial compilers come bundled with an editor.

PLOT PROGRAM

gnuplot is a command-line program. It is very easy to get a simple plot, but a bit tedious to tailor a plot to your liking. GUI-driven plotting programs are also available.

Quick Start

The first step is downloading the Beginner VH-1 code to your computer, which you can do from this link -> VH-0.5.

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 tar file 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

The workflow for the vast majority of problems you'll face when using VH-1 can be distilled to four steps:

1. Edit vhone.f90 (and other files) to reflect the specific problem you're simulating.
2. Compile the code using make
3. Edit indat as needed to make sure your simulation evolves as long as you want it to.
4. Run the executable file vhone
5. Plot the output for analysis

VH-0.5 comes preset to evolve the Sod shock tube problem, so you can bypass step 1 for a quick start.

To compile VH-0.5, type make. Compiling will 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.

The next step is to edit indat to control how many time steps the simulation will evolve and how often to output data. The copy that comes with VH-0.5 should have endtime set to 2.1e-01, and tprin set to 2.0e-02. This will result in 10 output files over the course of the run, but you can edit it to output more or fewer files.

With indat configured, you can run VH-1. Type vhone (or ./vhone if required) to run the simulation and evolve the Sod shock tube. This simple problem should take less than a second to finish. Once VH-1 finishes, typing ls will show you that VH-1 did indeed generate 10 output files as part of the simulation: DATAa1000.dat to DATAa1009.dat. Each of these files is text file containing grid values for location, density, pressure, and 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 with whatever package you prefer. As an example, type gnuplot at the command prompt, and then type
gnuplot> plot 'DATAa1005.dat' w lp, 'DATAa1007.dat' w lp, 'DATAa1009.dat' w lp
to get something that looks similar (but not quite) like

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

VH-0.5 comes with three files ( Makefile, indat, vhone.f90) and a subdirectory PPMLR. The PPMLR directory contains the guts of the hydrodynamic code. In a nutshell, PPMLR does one time-update on one-dimensional arrays of density, pressure and velocity using the PPM Lagrange-Remap algorithm. The files in PPMLR should be left alone unless the problem you're working on requires changes to the subroutines there. It can be safely ignored for the duration of this tutorial. Makefile tells your computer how to compile VH-1 and indat is a text file defining job control parameters.

1. Read indat file for job control parameters
2. Open history file for recording metadata
3. Create the simulation grid
4. Set initial conditions
5. Loop over time (time = time + dt)
   1. Hydrodynamic update
   2. Compute maximum stable time step
   3. Output data to disk

READ INDAT
VH-1 uses a FORTRAN namelist to read in a set of job control parameters from the file indat.

NAMELIST / hinput / rstrt,prefix,ncycend,ndump,nprin,nmovie,endtime,tprin,tmovie
!-------------------------------------------------------------------------------
! Begin by reading input file (namelist) for job control values

open (unit=15,file='indat',status='old',form='formatted')
 read (15,nml=hinput)
close(15)
 

WRITE HISTORY
Open up a history file for recording simulation metadata. Anything written to FORTRAN unit 8 will be written into this file. Start by recording today's date.

! Name and open a file for recording the simulation metadata

hstfile  = prefix // 'hst'
open (unit=8,file=hstfile,form='formatted')
call date_and_time(todayis)
write (8,*) 'VH-1 run on ',todayis(5:6),' / ',todayis(7:8),' / ',todayis(1:4)
write (8,*) 
 

CREATE GRID
We start by defining the simulation grid, including the geometry, boundary conditions, number of grid zones and spatial coordinates of each zone.

! CREATE GRID

! Set some flags for geometry and boundary conditions

sweep  = 'x'   ! default is x-sweep (only option for 1D)
ngeom  = 0     ! Cartesian geometry
nleft  = 0     ! Reflecting at xmin
nright = 0     ! Reflecting at xmax

! Create a grid of imax zones, making room for 6 'ghost zones' on each end

imax = 60            ! total number of real zones on grid
nmin = 7             ! first real zone
nmax = imax + 6      ! last real zone
xmin = 0.            ! x value at left edge of grid
xmax = 1.0           ! x value at right edge of grid

! Initialize grid coordinates: xa0(n) is coordinate of left edge of zone n

dxfac = (xmax - xmin) / float(imax)  ! width of each zone
do n = nmin, nmax
  xa0(n) = xmin + (n-nmin)*dxfac 
  dx0(n) = dxfac
enddo

 

Geometry
The first step is to choose one of the predefined grid geometries using the integer flag ngeom. Options include Cartesian (0), Cylindrical radius (1), Spherical radius (2), Cylindrical polar angle (3), Spherical polar angle (4), and Spherical azimuthal angle (5).

Boundary Conditions
The integer flags nleft and nright set the boundary conditions on the left and right sides of the grid, respectively. Options include reflecting (0), zero gradients to facilitate smooth flow off the grid (1), fixed values (2), and periodic (3). Using option 2 requires setting values for fluid variables at the boundary. For the left side these are dinflo, pinflo, uinflo, vinflo, winflo, and for the right side these are dotflo, potflo, uotflo, votflo, wotflo.

Simulation Grid
The integer imax sets the number of zones in the grid. In order to implement boundary conditions, the code adds 6 ghost zones on each side of the grid. The integers nmin, nmax are the array indicies of the first and last real grid zones, 7 and imax+6. The variables xmin, xmax set the coordinates of the left and right side of the grid. With these parameters set, the code computes the coordinates of each grid zone. xa0(n) is the coordinate of the left edge of zone n, and dx0 is the width of zone n.

INITIAL CONDITIONS
It's now time to initialize the grid. Each one of the five fluid variables used in VH-1 must be set in every zone (from nmin to nmax) when VH-1 initializes.

! INITIAL CONDITIONS

! 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

Problem Parameters
In many cases the initial conditions will depend on some parameters. In the original version of the code set up for the Sod shock problem, these variables are pleft, pright, dleft, and dright. You will need to declare any such variables at the top of vhone.f90, and then initialize their values. The ratio of specific heats gam must also be initialized.

Initialize Variables
Loop over all real zones, from zone nmin to zone nmax, and set initial values for all five variables. The one-dimensional arrays representing the fluid variables are:

 p(n), the pressure

 r(n), the density (from the Greek \rho )

 u(n), the component of fluid velocity along our grid dimension

 v(n), the transverse components of the fluid velocity

 w(n), the transverse components of the fluid velocity

Initialize Time Step
Finally, a value for the initial time step is computed based on the Courant condition.

! Compute initial time step based on Courant condition to ensure stability

ridt = 0.
do n = nmin, nmax
 svel = sqrt(gam*p(n)/r(n))/dx0(n)
 xvel = abs(u(n)) / dx0(n)
 ridt = max(xvel,svel,ridt)
enddo
dt = courant / ridt 
  

Record Initial Conditions
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.

! Write out initial data to a file

filename = prefix // 'init.dat'
open(unit=3,file=filename,form='formatted')
do n = nmin, nmax
  write(3, 1003) xa0(n), r(n), p(n), u(n)
enddo
close(3)

! Log parameters of problem in history file

write (8,*) 'Sod shock tube in 1 dimension'
write (8,"('Grid dimensions: ',i4)") imax
write (8,*) 
write (8,*) 'Adiabatic index, gamma = ', gam
write (8,*) 'Pressure ratio is ', pright/pleft
write (8,*) 'Density ratio is ', dright/dleft
write (8,*) 
  

TIME EVOLUTION
After initializing the simulation time and some counters, we are ready to begin the hydrodynamic evolution. The time evolution is computed by looping over time in steps of duration dt. This main loop is controlled by the variable ncycle, and the evolution continues until the value of ncycle reaches ncycend.

!  MAIN COMPUTATIONAL LOOP

do while (ncycle < ncycend)

  if ( time + dt  >  endtime ) then 
    dt = endtime - time     ! set dt so time lands on endtime
    ncycend = ncycle-1      ! set ncycend so simulation stops after this step
  endif

  ! construct total energy; reset Lagrangian coordinates to Eulerian grid
  !    The Eulerian coordinates, xa0(n), will stay fixed; 
  !    The Lagrangian coordinates, xa(n), will move with the flow each time step
  do n = nmin, nmax
    e (n) = p(n)/(r(n)*gamm) + 0.5*(u(n)**2+v(n)**2+w(n)**2)
    xa(n) = xa0(n)
    dx(n) = dx0(n)
  enddo
  
  ! perform 1D hydro update with PPMLR algorithm
  call ppmlr

  ! update time and cycle counters
  time   = time  + dt
  timep  = timep + dt
  ncycle = ncycle + 1
  ncycp  = ncycp  + 1
  
  

The hydrodynamic update begins by calculating the total energy from the five primitive variables and resetting the Lagrangian coordinates to the original Eulerian values. The update itself is done by the call to subroutine ppmlr. This subroutine calculates the updated values (at time + dt) of the primitive fluid variables at each zone.

Time Step
At each time step a new value of dt is computed to ensure stability. The value of dt is controlled by the parameter courant, which must be less than one (typically courant=0.5). The size of the time step is limited to increase by no more than 10%.

  ! TIME STEP CONTROL
  
  ridt = 0.         ! searching for maximum, so start out with zero
  do n = nmin, nmax
   svel = sqrt(gam*p(n)/r(n))/dx0(n)
   xvel = abs(u(n)) / dx0(n)
   ridt = max(xvel,svel,ridt)
  enddo
  dtx  = courant / ridt     ! global time constraint for given courant parameter
  dt3  = 1.1 * dt           ! limiting constraint on rate of increase of dt 
  dt   = min( dt3, dtx )    ! use smallest required timestep 
  
  

Data Output
When the cycle counter ncycp reaches the value of nprin or the time counter timep reaches the value of tprin a copy of the simulation variables are written to a new file.

  ! DATA OUTPUT
  
  if ( ncycp >= nprin .or. timep >= tprin ) then    ! Print out a data file

    write(tmp1,"(i4)") nfile + 1000 ; nfile = nfile + 1
    filename = prefix // tmp1 // '.dat'
    write(8,6001) filename, time, ncycle

    open(unit=3,file=filename,form='formatted')
     do n = nmin, nmax
       write(3, 1003) xa0(n), r(n), p(n), u(n)
     enddo
    close(3)
  
    timep = 0.
    ncycp = 0
    
  endif
  
  

Filename
A file name is constructed using the five-letter prefix prefix read in from indat and the four-digit number nfile. The output file contains one line for each numerical zone. Each line has four numbers corresponding to the grid coordinate, density, pressure and velocity. After the data is written to the file the output counters timep and ncycp are reset to zero, and the value of nfile is incremented by one so that the next output file has a unique name.

Indat

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, and how often to output data.

&hinput
 rstrt   = 'no'
 prefix  = 'DATA0'
 ncycend = 400
 ndump   = 5000000
 nprin   = 40
 nmovie  = 1000000
 endtime = 9.9e+19
 tprin   = 9.9e+19
 tmovie  = 9.9e+19 /
    

A simulation can be controlled by cycle numbers using ncycend, ndump, nprin, nmovie or by simulation time using endtime, tprin, tmovie. For example, the simulation will stop as soon as it reaches a cycle number of ncycend or a simulation time of endtime, whichever comes first. In practice one sets these values to use cycles or time exclusively. In the case shown here the time values are set to very large numbers so that they are never reached and the control is determined by cycle numbers.

Several of these parameters are ignored in VH-0.5 since 1D runs take mere seconds to compute and there is no need for 2D animations.

rstrt

This parameter directs VH-1 to either start a new simulation (rstrt = 'no') or restart from an old simulation by reading in a restart file with a name of the form DATA0dxx where the last two charachters are the value of rstrt. [Ignored in VH-0.5]

prefix

This character string will be used to name all the output files. VH-1 is coded to use 5-character prefixes, so make sure yours is 5 characters long.

ncycend

The total number of time steps to evolve the simulation.

ndump

The number of time steps between 'restart dumps'. If running a simulation takes many days, you will want to checkpoint your simulation every few hours so that you don't have to start from the beginning if something goes wrong. [Ignored in VH-0.5]

nprin

The number of time steps between writing output files. In the example shown here, there will be 400/40 = 10 output files written to disk (plus the initial data).

nmovie

The number of time steps between writing a frame to a movie file. In 2D/3D the movie file(s) are 3D files (two spatial dimensions and one time) that can be stepped through to animate the simulation results. In 1D the output is a space-time 2D array. [Ignored in VH-0.5]

endtime

The value of the simulation time at which the code stops.

tprin

The interval of simulation time between writing output files.

tmovie

The interval of simulation time between writing animation frames. [Ignored in VH-0.5]