Solution of the Equations

We use a preconditioned conjugate-gradient method [7] to solve our discretized diffusion equations. The preconditioner is completely analogous to that used for the 2-D local support-operators scheme. [1] It is obtained simply by setting the off-diagonal elements of all the S -matrices, defined by Eq. (19), to zero. This makes the W -matrices diagonal, which ultimately results in a 7-point pure cell-center diffusion equation. That is to say that the balance equation for each cell in the mesh interior contains only the cell-center intensity in that cell together with the cell-center intensities in the six cells sharing a face with that cell. Each face-center intensity is expressed in terms of the cell-center intensities in the two cells that share that face-center intensity. For instance, if we set the off-diagonal elements of the S -matrices to zero, Eqs. (26) and (27) yield:

$${\frac{{2D_{i+1,j,k}}}{{\Delta^L_{i+1,j,k}}}} \left(\vphantom{ \phi_{i+\frac{1}{2},j,k} - \phi_{i+1,j,k} }\right. \phi_{{i+\frac{1}{2},j,k}}^{} \phi_{{i+1,j,k}}^{} \left.\vphantom{ \phi_{i+\frac{1}{2},j,k} - \phi_{i+1,j,k} }\right)$$ (55)

and

$${\frac{{2D_{i,j,k}}}{{\Delta^R_{i,j,k}}}} \left(\vphantom{ \phi_{i+\frac{1}{2},j,k} - \phi_{i,j,k} }\right. \phi_{{i+\frac{1}{2},j,k}}^{} \phi_{{i,j,k}}^{} \left.\vphantom{ \phi_{i+\frac{1}{2},j,k} - \phi_{i,j,k} }\right)$$ (56)

respectively, where

$$\Delta^{L}_{{i+1,j,k}} \left[\vphantom{ s_{L,L}^{LBD} V^{LBD} + s_{L,L}^{LTD} V^{LTD} + s_{L,L}^{LBU} V^{LBU} + s_{L,L}^{LTU} V^{LTU} }\right. \left.\vphantom{ s_{L,L}^{LBD} V^{LBD} + s_{L,L}^{LTD} V^{LTD} + s_{L,L}^{LBU} V^{LBU} + s_{L,L}^{LTU} V^{LTU} }\right]_{{i+1,j,k}}^{}$$ (57)

$$\Delta^{R}_{{i,j,k}} \left[\vphantom{ s_{R,R}^{RBD} V^{RBD} + s_{R,R}^{RTD} V^{RTD} + s_{R,R}^{RBU} V^{RBU} + s_{R,R}^{RTU} V^{RTU} }\right. \left.\vphantom{ s_{R,R}^{RBD} V^{RBD} + s_{R,R}^{RTD} V^{RTD} + s_{R,R}^{RBU} V^{RBU} + s_{R,R}^{RTU} V^{RTU} }\right]_{{i,j,k}}^{}$$ (58)

Substituting from Eqs. (55) and (56), into Eq. (45), we get the equation for $\phi_{{i+\frac{1}{2},j,k}}^{}$ :

$${\frac{{2D_{i,j,k}\left( \phi_{i+\frac{1}{2},j,k} - \phi_{i,j,k} \right)}}{{\Delta^R_{i,j,k}}}} {\frac{{2D_{i+1,j,k}\left( \phi_{i+\frac{1}{2},j,k} - \phi_{i+1,j,k} \right)}}{{\Delta^L_{i+1,j,k}}}}$$ (59)

Solving Eq. (59) for $\phi_{{i+\frac{1}{2},j,k}}^{}$ , we get:

$$\phi_{{i+\frac{1}{2},j,k}}^{} \left.\vphantom{\left( \phi_{i,j,k} \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \phi_{i+1,j,k}\frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} \right) }\right. \left(\vphantom{ \phi_{i,j,k} \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \phi_{i+1,j,k}\frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} }\right. \phi_{{i,j,k}}^{} {\frac{{D_{i,j,k}}}{{\Delta^R_{i,j,k}}}} \phi_{{i+1,j,k}}^{} {\frac{{D_{i+1,j,k}}}{{\Delta^L_{i+1,j,k}}}} \left.\vphantom{ \phi_{i,j,k} \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \phi_{i+1,j,k}\frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} }\right) \left.\vphantom{\left( \phi_{i,j,k} \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \phi_{i+1,j,k}\frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} \right) }\right/ \left(\vphantom{ \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} }\right. {\frac{{D_{i,j,k}}}{{\Delta^R_{i,j,k}}}} {\frac{{D_{i+1,j,k}}}{{\Delta^L_{i+1,j,k}}}} \left.\vphantom{ \frac{D_{i,j,k}}{\Delta^R_{i,j,k}} + \frac{D_{i+1,j,k}}{\Delta^L_{i+1,j,k}} }\right)$$ (60)

Thus we see from Eq. (60) that neglecting the off-diagonal elements of the S -matrices makes each interior-mesh face-center intensity a weighted-average of the two cell-center intensities adjacent to it. Substituting from Eq. (60) into Eqs. (55) and (56) we find that the face-area fluxes on the right and left faces of cells (i, j, k) and (i + 1, j, k) , respectively, can be expressed in terms of a difference between the cell-center intensities in those two cells:

$${\frac{{ D_{i+\frac{1}{2},j,k}}}{{\Delta_{i+\frac{1}{2},j,k}}}} \left(\vphantom{ \phi_{i+1,j,k} - \phi_{i,j,k} }\right. \phi_{{i+1,j,k}}^{} \phi_{{i,j,k}}^{} \left.\vphantom{ \phi_{i+1,j,k} - \phi_{i,j,k} }\right)$$ (61)

where

$${\frac{{1}}{{2}}} \left[\vphantom{ \left. \left( \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k... ...ht)\right / \left( \Delta^{R}_{i,j,k} + \Delta^{L}_{i+1,j,k} \right) }\right. \left.\vphantom{ \left( \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k}} + \frac{\Delta^{L}_{i+1,j,k}}{D_{i+1,j,k}} \right)}\right. \left(\vphantom{ \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k}} + \frac{\Delta^{L}_{i+1,j,k}}{D_{i+1,j,k}} }\right. {\frac{{\Delta^{R}_{i,j,k}}}{{D_{i,j,k}}}} {\frac{{\Delta^{L}_{i+1,j,k}}}{{D_{i+1,j,k}}}} \left.\vphantom{ \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k}} + \frac{\Delta^{L}_{i+1,j,k}}{D_{i+1,j,k}} }\right) \left.\vphantom{ \left( \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k}} + \frac{\Delta^{L}_{i+1,j,k}}{D_{i+1,j,k}} \right)}\right/ \left(\vphantom{ \Delta^{R}_{i,j,k} + \Delta^{L}_{i+1,j,k} }\right. \Delta^{{R}}_{{i,j,k}} \Delta^{{L}}_{{i+1,j,k}} \left.\vphantom{ \Delta^{R}_{i,j,k} + \Delta^{L}_{i+1,j,k} }\right) \left.\vphantom{ \left. \left( \frac{\Delta^{R}_{i,j,k}}{D_{i,j,k... .../ \left( \Delta^{R}_{i,j,k} + \Delta^{L}_{i+1,j,k} \right) }\right]^{{-1}}_{}$$ (62)

and

$$\Delta_{{i+\frac{1}{2},j,k}}^{} {\frac{{\Delta^{R}_{i,j,k} + \Delta^{L}_{i+1,j,k}}}{{2}}}$$ (63)

Thus we see that each interior-mesh face-area flux can be expressed in terms of a difference between the two adjacent cell-center intensities. Given Eq. (61), it is not difficult to see that the balance equation, Eq. (38), can be constructed from the cell-center intensities alone, resulting in a 7-point cell-center diffusion discretization for each cell on the mesh interior. In particular, the balance equation for cell (i, j, k) (and the equation for $\phi_{{i,j,k}}^{}$ ) is

$${\frac{{\partial \phi_{i,j,k}}}{{\partial t}}}$$Vi, j, k +	- $${\frac{{ D_{i+\frac{1}{2},j,k}}}{{\Delta_{i+\frac{1}{2},j,k}}}}$$$$\left(\vphantom{ \phi_{i+1,j,k} - \phi_{i,j,k} }\right.$$$$\phi_{{i+1,j,k}}^{}$$ - $$\phi_{{i,j,k}}^{}$$$$\left.\vphantom{ \phi_{i+1,j,k} - \phi_{i,j,k} }\right)$$ + $${\frac{{D_{i-\frac{1}{2},j,k}}}{{\Delta_{i-\frac{1}{2},j,k}}}}$$$$\left(\vphantom{ \phi_{i,j,k} - \phi_{i-1,j,k} }\right.$$$$\phi_{{i,j,k}}^{}$$ - $$\phi_{{i-1,j,k}}^{}$$$$\left.\vphantom{ \phi_{i,j,k} - \phi_{i-1,j,k} }\right)$$
	- $${\frac{{D_{i,j+\frac{1}{2},k}}}{{\Delta_{i,j+\frac{1}{2},k}}}}$$$$\left(\vphantom{ \phi_{i,j+1,k} - \phi_{i,j,k} }\right.$$$$\phi_{{i,j+1,k}}^{}$$ - $$\phi_{{i,j,k}}^{}$$$$\left.\vphantom{ \phi_{i,j+1,k} - \phi_{i,j,k} }\right)$$ + $${\frac{{D_{i,j-\frac{1}{2},k}}}{{\Delta_{i,j-\frac{1}{2},k}}}}$$$$\left(\vphantom{ \phi_{i,j,k} - \phi_{i,j-1,k} }\right.$$$$\phi_{{i,j,k}}^{}$$ - $$\phi_{{i,j-1,k}}^{}$$$$\left.\vphantom{ \phi_{i,j,k} - \phi_{i,j-1,k} }\right)$$
	- $${\frac{{D_{i,j,k+\frac{1}{2}}}}{{\Delta_{i,j,k+\frac{1}{2}}}}}$$$$\left(\vphantom{ \phi_{i,j,k+1} - \phi_{i,j,k} }\right.$$$$\phi_{{i,j,k+1}}^{}$$ - $$\phi_{{i,j,k}}^{}$$$$\left.\vphantom{ \phi_{i,j,k+1} - \phi_{i,j,k} }\right)$$ + $${\frac{{D_{i,j,k-\frac{1}{2}}}}{{\Delta_{i,j,k-\frac{1}{2}}}}}$$$$\left(\vphantom{ \phi_{i,j,k} - \phi_{i,j,k-1} }\right.$$$$\phi_{{i,j,k}}^{}$$ - $$\phi_{{i,j,k-1}}^{}$$$$\left.\vphantom{ \phi_{i,j,k} - \phi_{i,j,k-1} }\right)$$	= Qi, j, kVi, j, k   .

(64)

To obtain the analog of Eq. (61) for a cell face on the outer mesh boundary, we again consider a cell (1, j, k) , whose left face is on the boundary with its other faces in the mesh interior. Substituting from Eqs. (53) and (55) into Eq. (54), we obtain the equation for $\phi_{{\frac{1}{2},j,k}}^{}$ :

$${\frac{{2 D_{1,j,k}}}{{\Delta_{0,j,k}^{R}}}} \left(\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{\frac{1}{2},j,k} }\right. \phi_{{\frac{1}{2},j,k}}^{} \phi^{e}_{{\frac{1}{2},j,k}} \left.\vphantom{ \phi_{\frac{1}{2},j,k} - \phi^e_{\frac{1}{2},j,k} }\right) {\frac{{2 D_{1,j,k}}}{{\Delta^L_{1,j,k}}}} \left(\vphantom{ \phi_{\frac{1}{2},j,k} - \phi_{1,j,k} }\right. \phi_{{\frac{1}{2},j,k}}^{} \phi_{{1,j,k}}^{} \left.\vphantom{ \phi_{\frac{1}{2},j,k} - \phi_{1,j,k} }\right)$$ (65)

where

$$\Delta_{{0,j,k}}^{R} {\frac{{2 d^e_{0,j,k}}}{{\Vert \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{A}$}^L_{1,j,k} \Vert}}}$$ (66)

Solving Eq. (65) for $\phi_{{\frac{1}{2},j,k}}^{}$ . we get

$$\phi_{{\frac{1}{2},j,k}}^{} \left.\vphantom{\left( \phi_{0,j,k} \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \phi_{1,j,k}\frac{D_{1,j,k}}{\Delta^L_{1,j,k}} \right) }\right. \left(\vphantom{ \phi_{0,j,k} \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \phi_{1,j,k}\frac{D_{1,j,k}}{\Delta^L_{1,j,k}} }\right. \phi_{{0,j,k}}^{} {\frac{{D_{1,j,k}}}{{\Delta^R_{0,j,k}}}} \phi_{{1,j,k}}^{} {\frac{{D_{1,j,k}}}{{\Delta^L_{1,j,k}}}} \left.\vphantom{ \phi_{0,j,k} \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \phi_{1,j,k}\frac{D_{1,j,k}}{\Delta^L_{1,j,k}} }\right) \left.\vphantom{\left( \phi_{0,j,k} \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \phi_{1,j,k}\frac{D_{1,j,k}}{\Delta^L_{1,j,k}} \right) }\right/ \left(\vphantom{ \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \frac{D_{1,j,k}}{\Delta^L_{1,j,k}} }\right. {\frac{{D_{1,j,k}}}{{\Delta^R_{0,j,k}}}} {\frac{{D_{1,j,k}}}{{\Delta^L_{1,j,k}}}} \left.\vphantom{ \frac{D_{1,j,k}}{\Delta^R_{0,j,k}} + \frac{D_{1,j,k}}{\Delta^L_{1,j,k}} }\right)$$ (67)

Substituting from Eq. (67) into Eqs. (53) and (55), we obtain the desired expression for the face-area flux component on a boundary face:

$${\frac{{D_{\frac{1}{2},j,k}}}{{\Delta_{\frac{1}{2},j,k}}}} \left(\vphantom{ \phi_{\frac{1}{2},j,k} - \phi_{0,j,k}^e }\right. \phi_{{\frac{1}{2},j,k}}^{} \phi_{{0,j,k}}^{e} \left.\vphantom{ \phi_{\frac{1}{2},j,k} - \phi_{0,j,k}^e }\right)$$ (68)

where

$${\frac{{1}}{{2}}} \left[\vphantom{ \left. \left( \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k... ...ight)\right / \left( \Delta^{R}_{0,j,k} + \Delta^{L}_{1,j,k} \right) }\right. \left.\vphantom{ \left( \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k}} + \frac{\Delta^{L}_{1,j,k}}{D_{1,j,k}} \right)}\right. \left(\vphantom{ \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k}} + \frac{\Delta^{L}_{1,j,k}}{D_{1,j,k}} }\right. {\frac{{\Delta^{R}_{0,j,k}}}{{D_{1,j,k}}}} {\frac{{\Delta^{L}_{1,j,k}}}{{D_{1,j,k}}}} \left.\vphantom{ \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k}} + \frac{\Delta^{L}_{1,j,k}}{D_{1,j,k}} }\right) \left.\vphantom{ \left( \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k}} + \frac{\Delta^{L}_{1,j,k}}{D_{1,j,k}} \right)}\right/ \left(\vphantom{ \Delta^{R}_{0,j,k} + \Delta^{L}_{1,j,k} }\right. \Delta^{{R}}_{{0,j,k}} \Delta^{{L}}_{{1,j,k}} \left.\vphantom{ \Delta^{R}_{0,j,k} + \Delta^{L}_{1,j,k} }\right) \left.\vphantom{ \left. \left( \frac{\Delta^{R}_{0,j,k}}{D_{1,j,k... ...t / \left( \Delta^{R}_{0,j,k} + \Delta^{L}_{1,j,k} \right) }\right]^{{-1}}_{}$$ (69)

and $\Delta_{{\frac{1}{2},j,k}}^{}$ is given by Eq. (63) evaluated with i = 0 and Eq. (66). This completes the derivation of the approximate cell-center diffusion scheme used to precondition the full cell-center/cell-face scheme.

Once the 7-point discretization equations have been solved for the cell-center intensities, the face-center intensities can be directly calculated. In particular, Eq. (60) and its analogs for the Bottom/Top and Down/Up faces are used to calculate the face-center intensities on the mesh interior, while Eq. (65) and its analogs for the Bottom/Top and Down/Up faces are used to calculate the face-center intensities on the outer mesh boundary. The 7-point cell-center discretization yields an SPD coefficient matrix. Most importantly, since the S -matrices are diagonal when the mesh is orthogonal, it follows that the 7-point cell-center discretization is equivalent to the full cell-center/face-center discretization whenever the mesh is orthogonal. Thus our preconditioner can be expected to be very effective if the mesh is not too skewed. Our preconditioning system costs much less to solve than the full system because the coefficient matrix of the preconditioning system has roughly one-fourth as many rows and one-sixth as many elements as the full-system coefficient matrix. Computational results presented in the next section confirm this expectation.

In closing this section we note that when the 7-point system is used for preconditioning purposes, an inhomogeneous source term will generally appear in both the cell-center and face-center intensity equations. We did not include such a source in our derivation of the face-center intensity equations because they do not appear in standard calculations. One must remember to include these sources before the face-center intensities are eliminated to obtain the 7-point cell-center system. This matter is extensively discussed for the 2-D case in [1].