# 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 is the main program 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.

# vhone.f90

The basic structure of the vhone.f90 looks like- Read indat file for job control parameters
- Open history file for recording metadata
- Create the simulation grid
- Set initial conditions
- Loop over time (time = time + dt)
- Hydrodynamic update
- Compute maximum stable time step
- 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.

Example: read indat file |
---|

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.

Example grid initialization code |
---|

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

Example grid initialization code |
---|

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

Example: initialization code for the Sod shock tube problem |
---|

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

Example grid initialization code |
---|

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

Example: Record initial data and parameters |
---|

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

Example: Main Computational Loop |
---|

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

Example: Calculate time step for next loop |
---|

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

Example: Data output |
---|

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

**THE END! THAT'S THE ENTIRE CODE!**

# Makefile

The file Makefile is a text file that contains instructions for the operating system to compile and link the subroutines of VH-1. The contents of the Makefile must reflect your specific computer system. For the purposes of VH-0.5, the only machine specific parameter in the Makefile is the name of the FORTRAN compiler, which is assumed to be gfortran. If your system does not have the gfortran compiler installed 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.

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

Example indat for the Sod shock tube problem |
---|

&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]