Appendix

The purpose of this appendix is to demonstrate that the coefficient matrix for our support-operators method is SPD (symmetric positive-definite) . This is achieved in the following manner. First we demonstrate that the W matrix is SPD. Next we show that the the coefficient matrix for a single-cell problem with reflective boundary conditions is SPS (symmetric positive-semidefinite) with a one-dimensional null space consisting of any set of spatially-constant intensities. At this point the demonstration becomes perfectly analogous to that given in [1] for the 2-D case. We conclude the 3-D demonstration by giving a brief description of the final steps. The full details of these steps are given in [1].

The following mathematical preliminaries are discussed in [7]. A matrix, B is symmetric if and only if

B = B^t   . (75)

A matrix, B , is SPD if and only if it is symmetric and it satisfies

$$\hat{{X}}^{t}_{} \hat{{X}} \mbox{for all vectors $\hat{X}$}$$ (76)

A matrix, B , is SPS if and only if it is symmetric and it satisfies

$$\hat{{X}}^{t}_{} \hat{{X}} \geq \mbox{for all vectors $\hat{X}$}$$ (77)

Thus every SPD matrix is also SPS. Assume that a square matrix, B , can be expressed in terms of a square matrix, K , as follows:

B = K^tK   . (78)

Then if K is not invertible, B is SPS but not SPD, and if K is invertible, B is SPD.

We begin the overall demonstration by showing that the matrix given in Eq. (35), W , is SPD. It suffices to show that its inverse, explicitly given in Eqs. (26) through (31), is SPD. We begin the construction of W^-1 by considering Eq. (25) and the S -matrices that appear in it. Each of the S -matrices is a 3 x 3 matrix that is uniquely associated with a vertex, and each of these matrices operates on a 3-vector composed of the face-area flux components associated with that vertex. We now re-express these 3 x 3 matrices as 6 x 6 matrices by having them operate on a vector composed of all six face-area flux components associated with the cell. For instance, the matrix S^LBD operates on the following vertex face-area flux vector:

$$\hat{{F}}^{{LBD}}_{} \left(\vphantom{ f^L,f^B,f^D }\right. \left.\vphantom{ f^L,f^B,f^D }\right)^{t}_{}$$ (79)

We want to redefine S^LBD so that it operates on the global vector of flux components:

$$\hat{{\mathcal F}} \left(\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right. \left.\vphantom{ f^L,f^R,f^B,f^T,f^D,f^U }\right)^{t}_{}$$ (80)

This is easily accomplished via a 3 x 6 matrix that we denote as P^LBD . In particular, the 6 x 6 version of S^LBD is given by

S^LBD_{6 x 6} = P^LBDtS^LBDP^LBD   , (81)

where

P^LBD_{L, L} = P^LBD_{B, B} = P^LBD_{D, D} = 1   , (82)

and all other elements of P^LBD are zero. The matrix S^LBD_{6 x 6} is explicitly given by

$$\left[\vphantom{ \begin{array}{cccccc} s_{L,L} & 0 & s_{L,B} & 0 ... ... s_{D,B} & 0 & s_{D,D} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} }\right. \begin{array}{cccccc} s_{L,L} & 0 & s_{L,B} & 0 & s_{L,D} & 0 \\ ... ...,L} & 0 & s_{D,B} & 0 & s_{D,D} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \left.\vphantom{ \begin{array}{cccccc} s_{L,L} & 0 & s_{L,B} & 0 ... ... s_{D,B} & 0 & s_{D,D} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} }\right]$$ (83)

For the general case, the matrix P is most easily defined with respect to the matrix S using numeric indices. To do this we simply number all vector components in the usual sequential manner, e.g.,

$$\left(\vphantom{ f^L,f^B,f^D }\right. \left.\vphantom{ f^L,f^B,f^D }\right)^{t}_{} \rightarrow \left(\vphantom{ f_1,f_2,f_3 }\right. \left.\vphantom{ f_1,f_2,f_3 }\right)^{t}_{}$$ (84)

and

$$\left(\vphantom{ f^L, f^R, f^B, f^T, f^D, f^U }\right. \left.\vphantom{ f^L,f^R,f^B,f^T,f^D,f^U }\right)^{t}_{} \rightarrow \left(\vphantom{ f_1,f_2,f_3,f_4,f_5,f_6 }\right. \left.\vphantom{ f_1,f_2,f_3,f_4,f_5,f_6 }\right)^{t}_{}$$ (85)

Using this numeric indexing, the matrix P is defined for the general case as follows: If the i 'th component of the local vector $\hat{{F}}^{{vertex}}_{}$ associated with S^vertex is the j 'th component of the global vector $\hat{{\mathcal{F}}}$ , then

p_{i, j} = 1   , (86)

otherwise

p_{i, j} = 0   . (87)

It is convenient at this point to assign the vertices with the indices LBD, RBD, LTD, RTD, LBU, RBU, LTU, RTU, to the respective numeric indices 1, 2, 3, 4, 5, 6, 7, 8. This enables us to re-express Eq. (25) as follows:

$$\hat{{\mathcal{H}}}^{t}_{} \hat{{\Phi}} \sum_{{n=1}}^{{8}} \hat{{\mathcal{H}}}^{t}_{} \hat{{\mathcal{F}}} \hat{{\mathcal{H}}}^{t}_{} \left(\vphantom{ \phi^C \hat{1} }\right. \phi^{C}_{} \hat{{1}} \left.\vphantom{ \phi^C \hat{1} }\right)$$ (88)

where n is the numeric vertex index, and where

$$\hat{{1}} \left(\vphantom{ 1,1,1,1,1,1 }\right. \left.\vphantom{ 1,1,1,1,1,1 }\right)^{t}_{}$$ (89)

$$\hat{{\Phi}} \left(\vphantom{ \phi^L, \phi^R, \phi^B, \phi^T, \phi^D, \phi^U }\right. \phi^{L}_{} \phi^{R}_{} \phi^{B}_{} \phi^{T}_{} \phi^{D}_{} \phi^{U}_{} \left.\vphantom{ \phi^L, \phi^R, \phi^B, \phi^T, \phi^D, \phi^U }\right)^{t}_{}$$ (90)

$$\hat{{\mathcal{H}}} \left(\vphantom{ h^L, h^R, h^B, h^T, h^D, h^U }\right. \left.\vphantom{ h^L, h^R, h^B, h^T, h^D, h^U }\right)^{t}_{}$$ (91)

Since Eq. (88) must hold for all possible $\hat{{\mathcal{H}}}$ , it follows that

$$\hat{{\Phi}} \left[\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right. \sum_{{n=1}}^{{8}} \left.\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right] \hat{{\mathcal{F}}} \phi_{C}^{} \hat{{1}}$$ (92)

Further manipulating Eq. (92), we obtain

$$\left[\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right. \sum_{{n=1}}^{{8}} \left.\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right] \hat{{\mathcal{F}}} \Delta \hat{{\Phi}}$$ (93)

where $\Delta$$\hat{{\Phi}}$ is defined by Eq. (34). Comparing Eqs. (32) and (93) it follows that

$$\left[\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right. \sum_{{n=1}}^{{8}} \left.\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{S}_n \mathbf{P}_n }\right]$$ (94)

From Eq. (19) it follows that each 3 x 3 S -matrix is the product of a matrix A and its transpose. Substituting from Eq. (19) into Eq. (94), we get,

$$\left[\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{A}_n^t \mathbf{A}_n \mathbf{P}_n }\right. \sum_{{n=1}}^{{8}} \left.\vphantom{ \sum_{n=1}^{8} V_n \mathbf{P}_n^t \mathbf{A}_n^t \mathbf{A}_n \mathbf{P}_n }\right]$$ W^-1 =

$$\left[\vphantom{ \sum_{n=1}^{8} V_n \left( \mathbf{A}_n \mathbf{P}_n \right)^t \left( \mathbf{A}_n \mathbf{P}_n \right) }\right. \sum_{{n=1}}^{{8}} \left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right. \left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right)^{t}_{} \left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right. \left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right) \left.\vphantom{ \sum_{n=1}^{8} V_n \left( \mathbf{A}_n \mathbf{P}_n \right)^t \left( \mathbf{A}_n \mathbf{P}_n \right) }\right]$$ = (95)

Since

• the matrix, $\left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right.$A_nP_n$\left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right)^{t}_{}$$\left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right.$A_nP_n$\left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right)$ , must be SPS for each value of n ,
• an SPS matrix multiplied by a positive scalar remains SPS,
• the diffusion coefficient will always be positive,
• the vertex volumes will be positive with any reasonably well-formed mesh,
• the A -matrices will be invertible with any well-formed mesh,
• the P -matrices are not invertible,

it follows from Eq. (95) that M_n must be SPS but not SPD for each value of n , where

$$\left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right. \left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right)^{t}_{} \left(\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right. \left.\vphantom{ \mathbf{A}_n \mathbf{P}_n }\right)$$ (96)

Substituting from Eq. (96) into Eq. (95) we find that W^-1 is a sum of matrices with each constituent matrix, M_n , being SPS:

$$\sum_{{n=1}}^{{8}}$$ (97)

It is shown in [1] that if a matrix is a sum of SPS matrices, it is SPS, and its null space is the intersection of the null spaces of the constituent matrices. From the definitions of the A -matrices and the P -matrices (see Eqs. (17), (86), and (87)), it follows that each M -matrix has a three-dimensional null space. For instance, the null space of M₁ (corresponding to the LBD corner) consists of any vector of the form

$$\hat{{\mathcal{F}}} \left(\vphantom{ 0,f^R,0,f^T,0,f^U }\right. \left.\vphantom{ 0,f^R,0,f^T,0,f^U }\right)^{t}_{}$$ (98)

where f^R , f^T , and f^U are free to take on any values. There is no one face-area flux component that is common to the null spaces of all eight M -matrices, so the intersection of their null spaces is the null set. This implies that W^-1 has an empty null space. Since it is also SPS, it follows that W^-1 is SPD. Finally, if W^-1 is SPD, then W must be SPD.

The next step in the demonstration is to construct the discrete diffusion equations for a single cell with reflective boundary conditions. We neglect the time-derivative term in Eq. (1) and consider only the diffusion operator. Let us assume a solution vector, $\hat{{\Phi}}$ , of the form given in Eq. (90). In order to use numeric indices for the coefficient matrix of the single-cell system, we number this vector in the usual manner, i.e.,

$$\left(\vphantom{ \phi^L,\phi^R,\phi^B,\phi^T,\phi^D,\phi^U,\phi^C }\right. \phi^{L}_{} \phi^{R}_{} \phi^{B}_{} \phi^{T}_{} \phi^{D}_{} \phi^{U}_{} \phi^{C}_{} \left.\vphantom{ \phi^L,\phi^R,\phi^B,\phi^T,\phi^D,\phi^U,\phi^C }\right)^{t}_{} \rightarrow \left(\vphantom{ \phi_1,\phi_2,\phi_3,\phi_4,\phi_5,\phi_6,\phi_7 }\right. \phi_{1}^{} \phi_{2}^{} \phi_{3}^{} \phi_{4}^{} \phi_{5}^{} \phi_{6}^{} \phi_{7}^{} \left.\vphantom{ \phi_1,\phi_2,\phi_3,\phi_4,\phi_5,\phi_6,\phi_7 }\right)^{t}_{}$$ (99)

The first 6 equations for a single cell are the equations for the face-center intensities. For a reflective boundary condition, these equations simply state that the face-area flux component on each face is zero. However, in analogy with Eqs. (45) through (47), we equivalently require that the negative of each component be zero. The W -matrix relates the face-area flux components to the differences between the cell-center intensity and the face-center intensities in accordance with Eq. (35). Thus the first 6 equations can be expressed in terms of the matrix W as follows:

$$\Delta \hat{{\Phi}}$$ (100)

where in accordance with Eqs. (34) and (99):

$$\Delta \hat{{\Phi}} \left(\vphantom{ \phi_7 -\phi_1, \phi_7 - \phi_2, \phi_7 - \phi_3, \phi_7 - \phi_4, \phi_7 - \phi_5, \phi_7 - \phi_6 }\right. \phi_{7}^{} \phi_{1}^{} \phi_{7}^{} \phi_{2}^{} \phi_{7}^{} \phi_{3}^{} \phi_{7}^{} \phi_{4}^{} \phi_{7}^{} \phi_{5}^{} \phi_{7}^{} \phi_{6}^{} \left.\vphantom{ \phi_7 -\phi_1, \phi_7 - \phi_2, \phi_7 - \phi_3, \phi_7 - \phi_4, \phi_7 - \phi_5, \phi_7 - \phi_6 }\right)^{{t}}_{}$$ (101)

Using Eqs. (100), and (101), one can easily construct the first six rows of the single-cell coefficient matrix, C , as follows:

$$\mbox{$i=1,6$} \mbox{$j=1,6$}$$ c_{i, j} = (102)

$$\sum_{{j=1}}^{{6}} \mbox{$i=1,6$}$$ c_{i, 7} = (103)

The seventh and last row of C corresponds to the steady-state balance equation, i.e., Eq. (38) with the time-derivative set to zero:

f^L + f^R + f^B + f^T + f^D + f^U = Q^CV   . (104)

Using Eqs. (35) and (101) through (104), we define the last row of the coefficient matrix:

$$\sum_{{i=1}}^{{6}} \mbox{$i=1,6$}$$ c_{7, j} = (105)

$$\sum_{{i=1}}^{{6}} \sum_{{j=1}}^{{6}}$$ c_{7, 7} = (106)

To summarize, the coefficient matrix takes the following block form:

$$\left[\vphantom{ \begin{array}{cc} \mathbf{W} & \mathbf{W}_r\\ \mathbf{W}_c & \mathbf{W}_{rc} \end{array} }\right. \begin{array}{cc} \mathbf{W} & \mathbf{W}_r\\ \mathbf{W}_c & \mathbf{W}_{rc} \end{array} \left.\vphantom{ \begin{array}{cc} \mathbf{W} & \mathbf{W}_r\\ \mathbf{W}_c & \mathbf{W}_{rc} \end{array} }\right]$$ (107)

where W_r is a 6 x 1 matrix obtained by summing the rows of W , W_c is a 1 x 6 matrix obtained by summing the columns of W , and W_rc is a 1 x 1 matrix obtained by summing all of the elements of W . Note that W_c is the transpose of W_r because W is symmetric. Thus C is symmetric. To prove that it is SPS, we need only show that it is positive-semidefinite. Towards this end we note that any vector $\hat{{\Phi}}$ can clearly be re-expressed as follows:

$$\hat{{\Phi}} \left(\vphantom{ \phi_1,\phi_2,\phi_3,\phi_4,\phi_5,\phi_6,\phi_7 }\right. \phi_{1}^{} \phi_{2}^{} \phi_{3}^{} \phi_{4}^{} \phi_{5}^{} \phi_{6}^{} \phi_{7}^{} \left.\vphantom{ \phi_1,\phi_2,\phi_3,\phi_4,\phi_5,\phi_6,\phi_7 }\right)^{t}_{} \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{}$$ (108)

where

$$\hat{{\Phi}}_{f}^{} \left(\vphantom{ \phi_1-\phi_7,\phi_2-\phi_7,\phi_3-\phi_7,\phi_4-\phi_7,\phi_5-\phi_7, \phi_6-\phi_7,0 }\right. \phi_{1}^{} \phi_{7}^{} \phi_{2}^{} \phi_{7}^{} \phi_{3}^{} \phi_{7}^{} \phi_{4}^{} \phi_{7}^{} \phi_{5}^{} \phi_{7}^{} \phi_{6}^{} \phi_{7}^{} \left.\vphantom{ \phi_1-\phi_7,\phi_2-\phi_7,\phi_3-\phi_7,\phi_4-\phi_7,\phi_5-\phi_7, \phi_6-\phi_7,0 }\right)^{t}_{}$$ (109)

and

$$\hat{{\Phi}}_{c}^{} \left(\vphantom{ \phi_7,\phi_7,\phi_7,\phi_7,\phi_7, \phi_7,\phi_7 }\right. \phi_{7}^{} \phi_{7}^{} \phi_{7}^{} \phi_{7}^{} \phi_{7}^{} \phi_{7}^{} \phi_{7}^{} \left.\vphantom{ \phi_7,\phi_7,\phi_7,\phi_7,\phi_7, \phi_7,\phi_7 }\right)^{t}_{}$$ (110)

Taking the inner product of $\hat{{\Phi}}$ with C $\hat{{\Phi}}$ , we get

$$\left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right.$$$$\hat{{\Phi}}_{f}^{}$$ + $$\hat{{\Phi}}_{c}^{}$$$$\left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right)^{t}_{}$$C$$\left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right.$$$$\hat{{\Phi}}_{f}^{}$$ + $$\hat{{\Phi}}_{c}^{}$$$$\left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right)$$ =

$$\hat{{\Phi}}^{t}_{f} \hat{{\Phi}}_{f}^{} \hat{{\Phi}}^{t}_{f} \hat{{\Phi}}_{c}^{} \hat{{\Phi}}^{t}_{c} \hat{{\Phi}}_{f}^{} \hat{{\Phi}}^{t}_{c} \hat{{\Phi}}_{c}^{}$$ (111)

It is easily verified that

$$\hat{{\Phi}}_{c}^{} \hat{{0}} \mbox{for all $\hat{\Phi}_c$}$$ (112)

Substituting from Eq. (112) into Eq. (111), we get

$$\left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right)^{t}_{} \left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right) \hat{{\Phi}}^{t}_{f} \hat{{\Phi}}_{f}^{} \hat{{\Phi^t}}_{c}^{} \hat{{\Phi}}_{f}^{}$$ (113)

Since

$$\hat{{\Phi}}^{t}_{c} \hat{{\Phi}}_{f}^{} \hat{{\Phi}}^{t}_{f} \hat{{\Phi}}_{c}^{}$$ (114)

Eq. (113) reduces to

$$\left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right)^{t}_{} \left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right) \hat{{\Phi}}_{f}^{t} \hat{{\Phi}}_{f}^{}$$ (115)

Using Eq. (107), it is easily shown that

$$\hat{{\Phi}}_{f}^{t} \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{{f6}}^{t} \hat{{\Phi}}_{{f6}}^{}$$ (116)

where

$$\hat{{\Phi}}_{{f6}}^{} \left(\vphantom{ \phi_1-\phi_7,\phi_2-\phi_7,\phi_3-\phi_7,\phi_4-\phi_7,\phi_5-\phi_7, \phi_6-\phi_7 }\right. \phi_{1}^{} \phi_{7}^{} \phi_{2}^{} \phi_{7}^{} \phi_{3}^{} \phi_{7}^{} \phi_{4}^{} \phi_{7}^{} \phi_{5}^{} \phi_{7}^{} \phi_{6}^{} \phi_{7}^{} \left.\vphantom{ \phi_1-\phi_7,\phi_2-\phi_7,\phi_3-\phi_7,\phi_4-\phi_7,\phi_5-\phi_7, \phi_6-\phi_7 }\right)^{t}_{}$$ (117)

Since W is SPD, it follows from Eqs. (114) and (117) that

$$\left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right)^{t}_{} \left(\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right. \hat{{\Phi}}_{f}^{} \hat{{\Phi}}_{c}^{} \left.\vphantom{ \hat{\Phi}_f + \hat{\Phi}_c }\right) \mbox{if $\hat{\Phi}_{f}=\hat{0}$\ }$$ =

> 0   ,otherwise. (118)

Thus C is positive-semidefinite. Since it is also symmetric, C is SPS. Note from Eq. (118) that the null space of C is spanned by all vectors $\hat{{\phi}}^{c}_{}$ . Following Eq. (110), it is clear that the null space of C is spanned by all vectors of constant intensity.

The remainder of the demonstration is identical to that given for the 2-D case in [1]. The final steps can be briefly described as follows:

1. Given a multi-cell mesh with N cells, the C -matrices for each cell are expanded to operate on the global vector of intensities for the entire mesh. This step is conceptually analogous to the expansion of the S^LBD matrix given in Eq. (83). Since the C -matrices are SPS, their expansions must be SPS.

2. It is shown that the sum of the expanded C -matrices represents the coefficient matrix for entire mesh with reflective conditions on the outer boundary faces. Since the global coefficient matrix is the sum of SPS matrices, it must be SPS. Furthermore, the null space of the full coefficient matrix must be equal to the intersection of the null spaces of the expanded C -matrices.

3. It is shown that the null space of the full coefficient matrix is spanned by all vectors of constant intensity. This is the correct result because the analytic diffusion operator has a null space spanned by all constant intensity functions if the reflective condition is imposed on the entire outer boundary. The analytic diffusion operator becomes invertible if the reflective condition is replaced with an extrapolated boundary condition on any portion of the outer boundary surface.

4. Finally, it is shown that if the reflective boundary condition is replaced with an extrapolated condition on any outer-boundary cell face, the expanded C -matrix for the cell containing the boundary face has a null space that is disjoint from the null spaces of all the other expanded C -matrices. Thus the intersection of the null spaces of all the expanded C -matrices is the null set. Since the global coefficient matrix is the sum of the expanded C -matrices, and the expanded C -matrices are SPS, it follows that the global coefficient matrix is SPD.