The system of partial differential is discretized by a finite volume
approach--see for example  or .
Space is divided into a series of quadrilateral (but not necessarily
rectangular) cells labelled .
This reduces the PDEs
to a large coupled system of ordinary differential equations.
These are integrated by a variant of the well-known 4th order Runge
A single time-step involves four stages like:
So at its most basic level the program for integrating the Euler equations consists of a series of steps like:
Another noteworthy thing is that these overall statements work naturally with aligned data--no communication is needed here. Out of the five stages enumerated above, only computation of R involves non-local terms (formally because of the use of averages across adjacent cells for the flow values at the faces). The code can be written easily using ghost regions, shifted indices, and the writeHalo() operation. Again it involves a single overall nest with a long body. A much-ellided outline is given in Figure 5. The optional arguments wlo, whi to Adlib.writeHalo() define the widths of the parts ghost regions that need updating (the default is to update the whole of the ghost regions of the array, whatever their width). In the current case these vectors both have value [1, 1]--because shifted indices displace one site in positive and negative x and y directions.
The arrays xnode and ynode hold coordinates of the cell vertices. Because these arrays are constant through the computation, the ghost regions of these arrays are initialized once during startup.
We will briefly discuss two other interesting complications: handling of so-called artificial viscosity, and imposition of boundary conditions.
Artificial viscosity (or artifical smoothing) is added to damp out a numerical instability in the Runge Kutta time-stepping scheme, which otherwise causes unphysical oscillatory modes associated with the discretization to grow. An accepted scheme adds small terms proportional to 2nd and 4th order finite difference operators to the update of . From the point of view of HPJava programming one interesting issue is that 4th order damping implies an update stencil requiring source elements offset two places from the destination element (unlike Figure 5, for example, where the maximum offset is one). This is handled by creating the array with ghost regions of width 2.
Implementing numerically stable boundary conditions for the Euler
equations is non-trivial.
In our implementation the domain of cells is rectangular, though the grid
is distorted into an irregular pipe profile by the choice of physical
coordinates attached to grid points (xnode, ynode distributed
HPJava has an additional control construct called at,
which can be used to update edges (it has other uses). The at
statement is a degenerate form of the overall statement. It only
``enumerates'' a single location in its specified range. To work
along the line , for example, one may write code like: