Code Basics

Interface Structure

A block face is the entire side of a block. A subface is a subset of a particular face and is defined by the block-side and the tangential index ranges. For example, a block with 20 × 20 × 20 nodes might have a subface on side is = 2 (j max side), and have tangential index ranges of k = [1, 10] and m = [11, 20]. An interface is a particular type of subface that has blocks on either side. An interface indicates how the two blocks are oriented. Hence, the interface must provide several pieces of information:

• • which blocks are adjacent;

  • the sides (1–6) of the blocks that meet at the interface, and;

  • tangential index ranges for both blocks.

There are several ways of storing the tangential index ranges; We use ICEMCFD’s convention. In this approach, an array of the form dit(1:2,1:3) is stored for each block on either side of the interface. The first tangential coordinate direction, i.e. one of (j, k, m) → (1, 2, 3), is stored in dit(1,1). The beginning of the index range is given by dit(1,2), and the number of nodes in the first tangential direction is dit(1,3). Note that the coordinate direction, dit(1,1), may be positive or negative; the sign indicates if the tangent indices should be traversed in increasing (positive) or decreasing (negative) order. The second tangential direction information is stored similarly, with the first array index in dit replaced with 2. For example, consider the interface shown in the figure above. For block 1 we might have dit(1,1) = 1 and dit(2,1) = 3; the first tangential direction is j and the second is m. I say “we might” only because the order and sign of the tangential indices may be different. Similarly, for block 2 we might have dit(1,1) = 2 and dit(2,1) = 3.

Boundary Structure

A boundary structure is a one-sided subface and must indicate:

• which block it belongs to;

• the type of boundary (far-field, symmetry, slip, no-slip, etc.);

• the side (1–6) of the block that the boundary lies on, and;

• tangential index ranges for the boundary subface.

Again, a dit(1:2,1:3) array is used to indicate the tangential coordinate information.

Structure Storage in the Code

The blocks are partitioned among the processors according to some algorithm. This partitioning gives rise to the Partition_Mod module. You can think of the Partition_Mod as a container for all the data needed by a individual processor. The block list array, blk, is the set of blocks assigned to a processor together with the set of necessary ghost blocks. The local blocks are listed first and the last local block has the index nBlk. Including the ghost blocks there are totBlk blocks in total. Similarly, the iface and bcface arrays are lists of interfaces and boundary faces. The boundary list is a straightforward array with nbcf boundary subface structures. The interface list is a bit more complicated, since an interface may be between two local blocks, or between a local block and a ghost block. We will call the former a local interface and the latter an external interface. The external interfaces are listed first in the iface array, and they are grouped according to the neighbouring processor responsible for the corresponding ghost block. The ifptr array can be used to access the interfaces corresponding to a particular neighbouring processor. For example, the interfaces for neighbour p are given by iface(ifptr(p):ifptr(p+1)-1). Suppose there are nBPnbr neighbouring processors. Then iface(ifptr(nBPnbr+1)) is the first internal interface and iface(nface) is the last internal interface, where nface is the total number of interfaces stored locally. The partition module contains more than lists for the basic structures: this module must coordinate the entire solution process. It is responsible for the highest level subroutines including initialization, outer (Newton) iterations, and residual evaluation. This module also holds the Jacobian and precondition matrix data types of the pKISlib library.

Additional Information

Stencil Numbers

I define the stencil of a node to be the set of neighbouring nodes used to construct the discrete flow equations. Consider, in 2-D, a second-order accurate approximation to the Laplace equation on a uniform Cartesian grid. The discrete equation at node (j, k) has the stencil {(j, k), (j −1, k), (j +1, k), (j, k−1), (j, k+1)}. Each node in the full stencil is assigned a binary bit; 0, 1, 2, 4, 8 in this case. If a node is absent from the stencil, say at a boundary, then the bit corresponding to that neighbouring node is turned off, i.e. set to zero. Hence, adding all the binary bits produces a number which encodes the stencil information. This may not seem very useful for the second-order accurate Laplace equation, but it becomes very handy in 3 dimensions and higher-order approximations. Why? Once you learn a bit about the compressed-sparse-row and block-compressed-row matrix storage schemes, you may start to appreciate the usefulness of a stencil number. Basically, it tells you (quickly) where to store neighbour information in the compressed, 1D array for the sparse matrix. The best way to appreciate how it works is to look at the code. Warning: the default integer in most implementations of the Fortran compiler will only be 32 bits. Thus, if you plan on using the stencil integer, your stencil had best not include more than 31 positions!

Finding Tangential Coordinate Directions

Often, when dealing with the sides of a block, it is necessary to find the coordinate indices corresponding to the normal and tangential directions. The calculation for finding the normal coordinate direction was described earlier; just add 1 to the side index and divide by two. Suppose this produces the coordinate direction index di. The two tangential coordinate directions can be found easily as follows: it1 = mod(di,1)+1 and it2 = mod(di+1,1)+1. Note, that the coordinate system (di,it1,it2) obeys a right-hand-rule. This is occasionally very important. The indices it1 and it2 are also used frequently to describe the tangential coordinate directions on an interface or boundary, i.e. it1=dit(1,1) and it2=dit(2,1). In these cases, the system (di,it1,it2) may not obey a right-hand-rule.