The core of VH-1 is a 1D gas-dynamics code based on the PPMLR algorithm of Collela & Woodward. This core code is contained in the sub-directory PPMLR/ that comes with each version of VH-1. All versions of VH-1 use an identical copy of PPMLR/, making it straightforward to jump between versions. Here we present a description of PPMLR/, subroutine by subroutine. Equations referenced by number refer to CW (J Comp Phys, 1984).
ppmlr.f90
ppmlr.f90 is the main driver of the 1D code. It calls the various components that make up the PPMLR algorithm:- call boundary
- Fills in the ghost zones with appropriate values based on the boundary flags nleft, nright
- call paraset
- Computes the parabolic coefficients in each zone for later use in parabola.f90
- call volume
- Computes the volume of each zone vol(n) based on the geometry flag ngeom
- call flatten
- Computes a flattening coefficient in each zone flat(n) based on the presence of a shock
- call parabola
- Computes a parabolic fit in each zone for the pressure, density, and velocity
- call states
- Computes the left and right states for input to the Riemann solver by integrating the parabolic fits over the domain of the characteristics, C*dt
- call riemann
- Solves the Riemann problem using Newton iteration to find the time-averaged values of velocity and pressure at the Lagrangian zone faces.
- call evolve
- Uses the result of the Riemann problem to update the Lagrangian fluid equations
- call remap
- Remaps the updated fluid variables from the evolved Lagrangian grid back to the original Eulerian grid
boundary.f90
Based on the values of nleft, nright, this subroutine assigns values to the primary fluid variables in the ghost zones. It also assigns values to the coordinates of the ghost zones and computes the total energy based on the primary variables. The available options for boundary conditions are- nleft,nright
- = 0 Reflecting
- = 1 Zero gradient (smooth flow off the grid)
- = 2 Fixed values
- = 3 Periodic (left edge is mapped onto right edge)
- = 1 Zero gradient (smooth flow off the grid)
volume.f90
Compute the volume element of each zone based on the value of ngeom. This must be done for both the Lagrangian grid, dvol(n), and the Eulerian grid, dvol0(n). Since these grids are nominally the same at the start of the time step, they will have the same volume as well. They are calculated separately in order to allow for a moving grid where the final Eulerian coordinates are shifted from the coordinates at the beginning of the time step. Geometry options include- ngeom
- = 0 Cartesian
- = 1 Cylindrical radius
- = 2 Spherical radius
- = 3 Cylindrical polar angle
- = 4 Spherical polar angle
- = 5 Spherical azimuthal angle
- = 1 Cylindrical radius
flatten.f90
This subroutine looks for signs of strong shocks and sets the variable flat(n) between 0.0 (smooth flow) and 1.0 (strong shock) using the simplified method of CW (eqn. A.1 and A.2).
parabola.f90
This file includes two subroutines. The first is paraset, which calculates the geometric coefficients used in the second subroutine, parabola. There are 5 coefficients, stored in the array para(n,5), used in equations 1.6 and 1.7 of CW.
states.f90
This subroutine takes the parabolic profiles of pressure, density, and velocity given by (al,da,a6) and computes the left and right states of each variable for input to the Riemann solver. The ap (plus) state is the averaged value of a in the forward travelling wave (on the left of the zone boundary), and the an (negative) state is the value of a in the backward wave (on the right side of the zone boundary).First calculate the fraction of the zone to be averaged using the average sound speed in that zone times the timestep, divided by the width of the zone. A factor of 1/2 is included in the definition of Cdtdx to save multiplications later on. If the coordinate is angular, the zone width must be multiplied by the variable radius to get a physical dimension.
The states are then computed using the formulae in eq. (1.12) of CW, with the addition of force corrections to the velocity states. These corrections of delta tg/2 represent the time averaged acceleration of the fluid to first order only. At this point we also make sure that the pressure states remain positive by not letting them fall below smallp. If geometry corrections are required (currently implemented for radial cyclindrical and spherical), the velocity states are altered according to eq. (2.9) of CW . Note that other corrections may be desired at this point to account for other source terms, including forces and energy sinks and sources.
riemann.f90
The Riemann solver uses the left and right input states at each zone boundary and solves the appropriate Riemann problem for such a given configuration. The time-averaged values of pressure and velocity are found for N+1 zone interfaces; nmin -> nmax + 1.
evolve.f90
This subroutine uses umid and pmid to solve the finite difference form of the lagrangian equations of hydrodynamics given by eq. (2.10) in CW. The first step is to evolve the locations of the grid faces, xa, with the old coordinates temporarily stored in xa1. If the current coordinate direction is angular, the coordinates must be multiplied by radius to give physical dimensions. The old and new coordinates are used to compute the gravitational (and ficticious) acceleration at the beginning and end of the time step for use in evolving the velocity and energy density.Next, the new zone volumes are computed according to the geometry flag passed from ppm.f. At the same time the average area of the zone face, amid, is computed to account for non-planar geometry in the hydro equations.
Evolution of the density is particularly trivial. Because the total mass in the zone remains constant in a Lagrangian calculation, the change in the zone-averaged density is inversely proportional to the change in zone volume. The new density is protected from going negative by setting a floor value of smallr.
The velocity is updated according to the pressure gradient and the forces, both fictitious and external. The total energy is updated according to the enthalpy gradient and the external forces only. The fictitious forces do not change the total energy, but rather redistribute the momentum between the different coordinate directions. If included in the energy equation, the change due to fictitous forces in one direction should exactly cancel the change found in the orthogonal direction. The pressure is found by subtracting the kinetic energy from the new total energy, and protected from falling below the floor value, smallp. If the problem at hand is in planar coordinates and no external forces are used, one can simply comment out the lines containing grav and fict in the velocity and energy loops, 40 and 50.
remap.f90
forces.f90
This subroutine is called from evolve.f90 to calculate the zone-centered forces at the beginning and end of the timestep for use in the finite difference equations. It is also called from states.f90 to get zone-faced values of the forces to correct the input velocity states to be used in the Riemann solver. It uses the values of sweep and ngeom in global.mod to determine which loop to run, and then loops over the row calculating the appropriate grav. and fict. If no gravity (or other external force) is used in the problem, set grav = 0. If you are sweeping over the radial direction and one of the other directions is angular, then fict = centrifugal force. If you are calculating an angular row then fict = Coriolis. The quantity radius is the radius of the given row of angular zones, and is set in the calling sweep routine.
vh1mods.f90
This file contains the source code for three modules.global.mod contains variables that need to be shared outside of PPMLR/. For example, dt is computed outside of PPMLR/ but is used in states.f90 and evolve.f90. Variables are also declared here when they need to be accessed by dump.f90 for creating checkpoint files.
sweeps.mod contains variables that are only needed inside the 1D hydrodynamic code, PPMLR/.
sweepsize.mod contains only one parameter, the size of the one-dimensional arrays used in PPMLR/: maxsweep. The value of this parameter must be larger than the largest grid dimension (imax, jmax, kmax) by at least 12 to account for the 6 ghost zones on each side of the grid.