In many applications, one may wish to solve fluid flow problems in a non-rectangular domain. One way to do this is to map a logically Cartesian two- or three-dimensional grid to a domain which is non-rectangular. Classic examples of such grids include polar grids, spherical grids, mappings to toroidal shaped regions, or the constricted tube at the right. In each of these cases, the grid coordinates no longer coorrespond to the usual (x,y,z) coordinates of a Cartesian framework,but rather to other geometric terms such radius or angle.
There are two basic philosophies behind solving partial differential equations (PDEs) on a general mesh.
The first approach has the advantage in that once the PDE has been transformed to the new coordinate system, it is generally obvious how to approximate the derivatives in the new set of equations, since these derivatives are taken with respect to a uniform Cartesian mesh. The disadvantage, however, is that the transformed PDE will in general contain geometric terms which were not present in the original PDE. In the context of hyperbolic problems, these terms may act to modify wave speeds, possibly leading to stiffness problems, or appear as source terms, potentially requiring splitting techniques. Furthermore, the analytic mapping function, which is generally required, may have singularities which are not present in the original PDE. Despite the disadvantages, this approach is used widely, at least for simple mapped grids. The developers of Overture at the Lawrence Berkeley Laboratory make extensive use of this technique, and achieve some very impressive results on very complicated composite meshes.
The second approach has the advantage that the original PDE is solved in its original form. However, one must derive accurate approximations to the derivatives in the physical coordinates. Generally, the stencils involved in such an approach will vary in space. This can be a problem for implicit discretization schemes which take advantage of symmetric matrices, but for explicit algorithms, one wouldn't expect to see any performance degradation due to the spatially varying stencils.
This second approach is the one we use to develop two- and three-dimensional finite volume methods for quadrilateral and hexahedral meshes. One notable feature of our approach is that we do not require an analytic mapping function. The user only need to be able to provide the four or eight corners, in physical coordinates, of the mesh cell corresponding to the logical cell C_{ijk}.
Below, we provide details for the two- and three-dimensional implementations.
Here is an example, taken from the online examples list above, of a shallow water wave calculation on an annular shaped region.
Initial conditions for a shallow water wave calculation on an annular shaped region
Although the above grid was generated using an analytic mapping function, this function is not used in computations using Clawpack. Only the corners of mesh cells are required. We then assume that cell edges are straight lines, and use simple formulas for quadrilaterals to compute the volume, lengths of edges, and normal and tangent vectors to each face of the quadrilateral mesh cell. The volumes and lengths are then used automatically by the Clawpack software to accurately discretize the PDE on the mapped grid. The normal and tangent information is used to solve Riemann problems in directions aligned with the mesh cell edges.
Here is an example of a supersonic nozzle, initialized with a shock Mach number of 1.25:
Nozzle problem. Contours of the gas Mach number are plotted.
Hexahedral mesh cells from a torus (left) and a general mesh (right) with normal and two tangent vectors computed at center of faces. The faces on the left mesh cell are planar, whereas those on the right are not.
Our current focus is to make the existing 2d functionality available in the three-dimensional versions of the Clawpack and AMRClaw software. To this end, we have addressed many of the technical issues related to using the wave propagation algorithm of LeVeque on a general logically Cartesian, hexahedral mesh.
One of the main issues we had to address was how to obtain necessary geometric information from the general hexahedral mesh cell. In keeping with the spirit of the two-dimensional implementation, we only require that the user to specify, in physical coordinates, the eight corners of a mesh cell. From this corner information, we then compute the volume of the mesh cell, the area of the faces, and normal and tangent vectors at each face. Since in general, the faces of a hexahedral mesh cell are not planar, we use a trilinear map, obtained from the eight corners, to obtain all the necessary information. The trilinear map approximates the surfaces of each face of each mesh cell by a ruled surface. We compute the volume of the hexahedral mesh cell using a an exact formula, compute the area of the faces using an approximate formula (which is exact for planar faces) and finally use the Jacobian of the trilinear map and a Gram-Schmid procedure to compute an orthogonal set of normal and tangential vectors at the center of each face. This information is then used in the Riemann solvers to solve Riemann problems in directions aligned with these vectors.
We have some preliminary results on our algorithm, and so far the results look promising. The following is an example of an Euler calculation in an enclosed nozzle shaped region, with an initial shock speed Mach number > 1.
Evolution of the shock front in an enclosed, constricted tube. Contours of the gas Mach number are shown.
