VH-1: A beginner's guide

Quick Start

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

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

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

For this beginner's guide we will be working with the one-dimensional starter code in src/Starter. Note that all versions (Starter, Serial, Parallel) use the same core hydrodynamics code in src/PPMLR/.

The workflow for the vast majority of problems you'll face when using VH-1 can be distilled into five 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 vh1-starter
5. Plot the output for analysis

VH-1 comes preset to initialize the grid with a boring, steady-state solution of uniform density and pressure.

To compile VH-1, type make in the src/Starter directory. Compiling will create several new files in the directory: you will see some object files (*.o) and a few modules (*.mod). The executable binary named vh1-starter is moved up to the top-level directory. 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. First, give your simulation a unique name (prefix = ). You can ignore the rstrt when using the Starter code. Next, set controls for how long to run the simulation and how frequently to output data. Since this is a very boring simulation, lets set it to run for a total of 20 time steps (ncycend = 20) and output data every 10 time steps (nprin = 10). This will result in 3 output files over the course of the run (at time steps 0, 10, and 20). For now, leave the time controls (endtime and tprin) at large numbers so they are never used.

With indat configured, you can run VH-1. Type vh1-starter (or ./vh1-starter if required) to run the simulation. Once VH-1 finishes, go to the directory containing the output: cd output/prefix where prefix is the name you put into the indat file. Type ls to see all the output files generated by your simulation. The history file records the logistics of the simulation (type cat prefix.hst to view contents). You will also see three output files from the simulation: prefix_1001.dat to prefix_1003.dat. Each of these files is a text file containing grid values for location, density, pressure, and velocity (in column order from left to right). In this boring example all the numbers should bethe same. Since the .dat files are formatted text, they may be plotted with whatever package you prefer (e.g., gnuplot).

At this point, you have a functional copy of VH-1 on your machine. The following is a step-by-step introduction to the basics of Computational Fluid Dynamics. If you are already familiar with the basic concepts you can jump right to the Sod Shocktube problem, a standard test problem for CFD codes.

Getting Started With CFD

Advection

Let's begin with a simple advection problem where we advect a density bump downstream. Use a code editor to open up vhone.f90 in src/Starter.

CREATE GRID: For this section we will keep the grid unchanged: One hundred uniform zones spanning from 0 to 1. Use cartesian geometry (ngeom = 0), meaning a constant cross-sectional area - think flow down a uniform pipe. Set the boundary conditions to be periodic (nleftx = 3) so that whatever flows off the right edge of the grid flows onto the grid from the left side.

INITIAL CONDITIONS: The initial set up has a uniform density (damb = 1.0) and uniform pressure (pamb = 1.0/gam). What is the speed of sound given these initial conditions? We want to advect something downstream, so set a non-zero (but subsonic) value for the velocity (uamb = 0.1). Now create something to advect - a density perturbation that is a function of position. Here we use a Gaussian profile centered on X = 0.5 with a width of 0.1 and an amplitude of 0.1 (relative to the ambient density). In VH-1, xa0(n) is the position of the left edge of zone n and dx0(n) is the width of zone n. The position of the center of zone n (for uniform zone spacing) is therefore given by xa0(n)+0.5*dx0(n).

! initialize every zone on the grid with density, pressure, and velocity
do n = 1, nmax          ! loop over the whole grid
  r(n) = damb*(1.0 + 0.1*exp(-0.5*((xa0(n)+0.5*dx0(n)-0.5)/0.1)**2))
  p(n) = pamb
  u(n) = uamb
  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
  f(n) = 0.0            ! set initial flattening to zero (used for smoothing shocks)
enddo
 

We can use the same indat file as above. That only evolves for 20 cycles, but that is enough to see if our initial setup is working. Compile the code by typing make in src/Starter and then running the executable file vh1-starter in the top-level directory. The result should look like the following...

The first thing to do is convince yourself that your initial conditions are correct (the 1001 file). Then ask yourself if the subsequent evolution makes sense. Here you expect the initial profile to move to the right with whatever speed you set for uamb.

Now lets see how accurate this code is (or isn't). With the periodic boundary conditions this density profile will advect off the right edge of the grid and re-enter from the left side. Since the width of the grid is set to unity, it should take a simulation time of 1.0/uamb to go all the way around the grid and reach the initial starting point. Configure the indat file to do just that. Set ncycend and nprin to large numbers (99999) so they are never reached, and then set endtime and tprin to 1.0/uamb so that the run ends (and outputs) at the time when the density profile reaches the starting point.

The final result looks pretty similar to the initial profile, except near the peak of the profile.

Pressure perturbation

Now try a similar experiment but perturbing the pressure. Here we set the background velocity to zero (uamb = 0.0) and perturb the pressure rather than the density. Here it is important to use a small perturbation, say 1% or less. Compile the code and edit indat. Perhaps try a different name so that you keep the density advection data and start a new subdirectory in output. Decrease tprin to something like 0.1 so you get more frequent outputs and can track the early evolution. Below we plot the pressure (column #3) for the first few time outputs.

Sound wave

The above problem generated two sound waves moving in opposite directions. Now modify your initial conditions to produce a single wave propogating to the right. Start with a perturbed velocity profile using the same Gaussian function and an amplitude that is highly subsonic (maybe 1% the speed of sound). A little bit of acoustic theory tells us that if the velocity profile of a linear wave is given by u^* (normalized to the sound speed, but we intialized the background solution such that the sound speed is unity) then the density perturbation is given by 1 + u^* and the pressure perturbation is given by 1 + γu^*. Code up these solutions in the initial conditions section of vhone.f90, compile, configure indat, and run the code.

Boundary conditions

Use your propagating wave to test different boundary conditions. Replace the periodic boundaries with reflecting boundaries (nrightx = 0). What do you expect to happen in this case? Compile, configure indat (maybe change the name?), and run the code. Here is the resulting evolution of the velocity, ending at a time of unity when the reflected wave has reached the center of the grid. What do the other variables (density and pressure) look like?