The system of partial differential is discretized by a finite volume
approach--see for example [7] or [11].
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
Kutta scheme.
A single time-step involves four stages like:

Here is the volume a cell and , are coordinate differences between end-points of the face. Since the dependent variables and fluxes are defined at cell centers, their values at a cell face in equation 3 is are approximated as the average of the values from the two cells meeting at the face.

So at its most basic level the program for integrating the Euler equations consists of a series of steps like:

- Calculate , from current (via equations of state).
- Calculate from , , .
- Calculate from , , .
- Calculate from , .
- Update .

The four fields

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
arrays).
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:

The actual code in the body is a fairly complicated interpolation based on Riemann invariants. In general access to