VH-1: A beginner's guide

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.

vhone.f90

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

ngeom

= 0 Cartesian

= 1 Cylindrical radius

= 2 Spherical radius

= 3 Cylindrical polar angle

= 4 Spherical polar angle

= 5 Spherical azimuthal angle

Boundary Conditions
The integer flags nleft and nright set the boundary conditions on the left and right sides of the grid, respectively. Options include

nleft,nright

= 0 Reflecting

= 1 Zero gradient (smooth flow off the grid)

= 2 Fixed values

= 3 Periodic (left edge is mapped onto right edge)

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. To ensure numerical stability, the time step must be smaller than both the flow time across any zone and the sound-crossing time for any zone. The variable courant is normally set to 0.5 in vh1mods.f90.

! 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.

Makefile

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.

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]