Theory

Governing Equations

VERTEX-SN solves the multi-group form of the linear neutral particle transport equations with an isotropic scattering kernel,

{math} \begin{align}         \hat{\mathbf{\Omega}}_a \cdot \nabla \psi_a^g(\mathbf{r}) + \sigma^g(\mathbf{r}) \psi_a^g(\mathbf{r}) = q_e^g((\mathbf{r}),\mathbf{\Omega}_a) + \\ \sum_{g'=0}^G \sum_{l=0}^N \sigma_{sl}^{gg'}((\mathbf{r}))         \left[             \mathbf{Y}_l0^e(\mathbf{\Omega}_a)\phi_{l0}^{g'}(\mathbf{r}) + \sum_{m=1}^l              \left(                 \mathbf{Y}_lm^e(\mathbf{\Omega}_a)\phi_{lm}^{g'}(\mathbf{r}) +                 \mathbf{Y}_lm^o(\mathbf{\Omega}_a)\theta_{lm}^{g'}(\mathbf{r})             \right)         \right] \end{align}

The energy group is represented by \(g\), the discrete angle is given by \(\mathbf{\Omega}_a\), and the spatial position is given by \(\mathbf{r}\). The total number of energy groups is \(G\) and the aniostropic scattering order is \(N\). The angular flux in group \(g\) in the discrete angle \(\mathbf{\Omega}_a\) is given by \(\psi_a^g \equiv \psi^g(\mathbf{\Omega}_a)\). \(\sigma_{sl}^{gg'}\) represents the scattering cross section from group \(g'\) to group \(g\) for the \(l^{th}\) Legendre scattering moment.

The even and odd angular flux moments, \(\phi\) and \(\theta\) are defined by {math} \begin{align}     \phi_{\ell m} &= \int_{4\pi}Y^e_{\ell m}(\mathbf{\Omega})\psi(\mathbf{\Omega})\,d\mathbf{\Omega}\:,\quad m\ge0\:,\\     \vartheta_{\ell m} &= \int_{4\pi}Y^o_{\ell m}(\mathbf{\Omega})\psi(\mathbf{\Omega})\,d\mathbf{\Omega}\:,\quad m\gt0\:, \end{align} in 3D and by {math} \begin{align}     \phi_{\ell m} &= \int_{0}^{2\pi}d\varphi\int_0^1      \bigl(1+(-1)^{\ell+m}\bigr)Y^e_{\ell m}     (\theta,\varphi)\psi(\theta,\varphi)\,     d\xi\:,\quad m\ge0\:,\\     \vartheta_{\ell m} &= \int_{0}^{2\pi}d\varphi\int_0^1      \bigl(1+(-1)^{\ell+m}\bigr)Y^o_{\ell m}     (\theta,\varphi)\psi(\theta,\varphi)\,     d\xi\:,\quad m\gt0\:. \end{align} in 2D.

Operator Form

The discrete ordinates equation in operator form can be given by \(\mathbf{L}\Psi = \mathbf{MS}\Phi + Q\).

The dimensions of the operators are defined in discrete \((d)\) and moment \((m)\) space, {math} \begin{align}     d &= N_g\times N_a\times N_c\times N_b\:,\\     m &= N_g\times N_m\times N_c\times N_b\:, \end{align} so that the dimension of the matrices in the operator form are {math} (d\times d)(d\times 1) = (d\times m)(m\times m)(m\times 1)  + (d\times 1)\:.

Neglecting space (\(N_c\times N_b\)), we have {math} \mathbf{L}\begin{pmatrix} \Psi^0\\\Psi^1\\\Psi^2\\\vdots\\\Psi^G \end{pmatrix}= \begin{pmatrix} \mathbf{M}&0&0&0&0\\ 0&\mathbf{M}&0&0&0\\ 0&0&\mathbf{M}&0&0\\ 0&0&0&\ddots&0\\ 0&0&0&0&\mathbf{M} \end{pmatrix} \begin{pmatrix} \mathbf{S}^{00}&\mathbf{S}^{01}&\mathbf{S}^{02}&\cdots&\mathbf{S}^{0G}\\ \mathbf{S}^{10}&\mathbf{S}^{11}&\mathbf{S}^{12}&\cdots&\mathbf{S}^{1G}\\ \mathbf{S}^{20}&\mathbf{S}^{21}&\mathbf{S}^{22}&\cdots&\mathbf{S}^{2G}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \mathbf{S}^{G0}&\mathbf{S}^{G1}&\mathbf{S}^{G2}&\cdots&\mathbf{S}^{GG} \end{pmatrix} \begin{pmatrix} \Phi^0\\\Phi^1\\\Phi^2\\\vdots\\\Phi^G \end{pmatrix} + \begin{pmatrix} Q^0\\Q^1\\Q^2\\\vdots\\ Q^G \end{pmatrix}\:.

Here, {math} \begin{align} \mathbf{L} &= \begin{pmatrix} L^0&L^1&L^2\cdots&L^G\:\\ \end{pmatrix}^T\:.\quad (N_g\times 1)\:,\\ \Psi^g&=\begin{pmatrix} \psi^g_0&\psi^g_1&\psi^g_2&\cdots&\psi^g_A \end{pmatrix}^T\:.\quad (N_a\times 1)\:,\\ Q^g&=\begin{pmatrix} Q^g_0&Q^g_1&Q^g_2&\cdots&Q^g_A \end{pmatrix}^T\:,\quad (N_a\times 1)\:. \end{align}

The moment-to-discrete matrix \(\mathbf{M}\) has dimension \((N_a\times N_m)\) and is defined {math} \mathbf{M}=\begin{pmatrix} Y^e_{00}(\mathbf{\Omega}_0) &  Y^e_{10}(\mathbf{\Omega}_0) &  Y^o_{11}(\mathbf{\Omega}_0) &  Y^e_{11}(\mathbf{\Omega}_0) &  Y^e_{20}(\mathbf{\Omega}_0) &  \cdots & Y^o_{NN}(\mathbf{\Omega}_0) &  Y^e_{NN}(\mathbf{\Omega}_0) \\ Y^e_{00}(\mathbf{\Omega}_1) &  Y^e_{10}(\mathbf{\Omega}_1) &  Y^o_{11}(\mathbf{\Omega}_1) &  Y^e_{11}(\mathbf{\Omega}_1) &  Y^e_{20}(\mathbf{\Omega}_1) &  \cdots & Y^o_{NN}(\mathbf{\Omega}_1) &  Y^e_{NN}(\mathbf{\Omega}_1) \\ Y^e_{00}(\mathbf{\Omega}_2) &  Y^e_{10}(\mathbf{\Omega}_2) &  Y^o_{11}(\mathbf{\Omega}_2) &  Y^e_{11}(\mathbf{\Omega}_2) &  Y^e_{20}(\mathbf{\Omega}_2) &  \cdots & Y^o_{NN}(\mathbf{\Omega}_2) &  Y^e_{NN}(\mathbf{\Omega}_2) \\ \vdots&\vdots&\vdots&\vdots&\vdots&&\vdots&\vdots\\ Y^e_{00}(\mathbf{\Omega}_A) &  Y^e_{10}(\mathbf{\Omega}_A) &  Y^o_{11}(\mathbf{\Omega}_A) &  Y^e_{11}(\mathbf{\Omega}_A) &  Y^e_{20}(\mathbf{\Omega}_A) &  \cdots & Y^o_{NN}(\mathbf{\Omega}_A) &  Y^e_{NN}(\mathbf{\Omega}_A)  \end{pmatrix}\:.

The moments of the angular flux are calculated using the discrete-to-moment matrix \(\Phi = \mathbf{D}\Psi\). From the quadrature integration rules we have \(\mathbf{D} = \mathbf{M}^T\mathbf{W},\quad (N_m\times N_a)\).

where {math} \begin{align} \mathbf{W}&=\begin{pmatrix} w_0&w_1&w_2&\cdots&w_A \end{pmatrix}^T\:,\quad (N_a\times 1)\:. \end{align} Accordingly, we have

{math} \Phi^g=\begin{pmatrix} \phi^g_{00}&\phi^g_{10}&\vartheta^g_{11}&\phi^g_{11}&\phi^g_{20}& \cdots&\vartheta^g_{NN}&\phi^g_{NN} \end{pmatrix}^T\:,\quad (N_m\times1)\:. with cross sections, \((N_m\times N_m)\), defined {math} \mathbf{S}^{gg'} = \begin{pmatrix} \sigma_s{0}^{gg'}&0&0&0&0&0&0&0\\ 0&\sigma_s{1}^{gg'}&0&0&0&0&0&0\\ 0&0&\sigma_s{1}^{gg'}&0&0&0&0&0\\ 0&0&0&\sigma_s{1}^{gg'}&0&0&0&0\\ 0&0&0&0&\sigma_s{2}^{gg'}&0&0&0\\ 0&0&0&0&0&\ddots&0&0\\ 0&0&0&0&0&0&\sigma_s{N}^{gg'}&0\\ 0&0&0&0&0&0&0&\sigma_s{N}^{gg'} \end{pmatrix}\:.

FEM Equations

For any given group, angle, and source, the discrete ordinates transport equation has the following space-dependent form (where all other indices have been supressed) {math} \begin{align}     \hat{\Omega}\cdot\nabla\psi(\mathbf{r})+\sigma(\mathbf{r})\psi(\mathbf{r}) = s(\mathbf{r})\:. \end{align}

We note that this will still be true for time-dependent equations where we solve the equations using implicit sweep operations (the time-absorption get rolled into \(\sigma\)). In matrix-form, we have \(\mathbf{L}\psi = s\)

Solving this equation implicitly by direct inversion of \(\mathbf{L}\), $ \(psi`=:raw-latex:\)mathbf{L}`^{-1}s:, $ requires sweeping the grid in the direction of particle travel determined by {math} \begin{align} \Omega\cdot\mathbf{n}\lt 0\:,\quad&(\text{entering face})\:,\\ \Omega\cdot\mathbf{n}\gt 0\:,\quad&(\text{exiting face})\:.\\ \end{align}

The DG-FEM equations are: {math} \boxed{ \oint_{\partial V_e}b_i\hat{\mathbf{n}}\cdot\Omega\psi\,dA - \mathbb{T}\Psi + \sigma\mathbb{M}\Psi = \mathbb{M}S\:,\quad i=0,\ldots,B\:,} where we have expanded using basis functions \(b_j\:,\ j=0,\ldots,B\), {math} \begin{align} \psi &= \sum_{j=0}^Bb_j(\mathbf{r})\psi_e^{(j)}\:,\quad\mathbf{r}\in\partial V_e\:,\\ [\mathbb{T}]_{ij} &= \int_{V_e}\mathbf{v}\cdot\nabla b_i b_j\,dV\:,\\ [\mathbb{M}]_{ij} &= \int_{V_e}b_ib_j\,dV\:,\\ \Psi &= \begin{pmatrix} \psi_e^{(0)} & \psi_e^{(1)} & \cdots & \psi_e^{(B)} \end{pmatrix}^T\:,\\ S &= \begin{pmatrix} s_e^{(0)} & s_e^{(1)} & \cdots & s_e^{(B)} \end{pmatrix}^T\:. \end{align}

On unstructured grids, we numerically integrate to build the linear system \(\mathbf{P}\Psi=q\) in each element. Numerical integration is done in natural coordinates.

Here, we do everything for 2D quads; however this is completely generalizable to other element types and dimensionalities. Mapping between natural and physical coordinates is done via the Jacobian of the system: {math} \mathbf{J} = \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\\ \end{pmatrix}\:. For now, we are assuming isoparametric FEM representations in which the mesh and function space are represented by the same basis. And we can expand \((x,y)\) in the basis as follows $ x = \(sum\){i=0}^B b_i(:raw-latex:`xi`,:raw-latex:`eta`)x_i:,:raw-latex:`quad `y = :raw-latex:`sum`{i=0}^B b_i(\(xi\),:raw-latex:eta)y_i:. $ Taking the partial derivatives gives {math} \begin{matrix} \frac{\partial x}{\partial \xi} = \sum_i \frac{\partial b_i}{\partial \xi} x_i\:, & \frac{\partial y}{\partial \xi} = \sum_i \frac{\partial b_i}{\partial \xi} y_i\:, \\ \frac{\partial x}{\partial \eta} = \sum_i \frac{\partial b_i}{\partial \eta} x_i\:, & \frac{\partial y}{\partial \eta} = \sum_i \frac{\partial b_i}{\partial \eta} y_i\:. \end{matrix} And \(\mathbf{J}\) becomes {math} \mathbf{J} = \begin{pmatrix} \frac{\partial b_0}{\partial \xi} & \frac{\partial b_1}{\partial \xi} & \frac{\partial b_2}{\partial \xi} & \frac{\partial b_3}{\partial \xi}\\ \frac{\partial b_0}{\partial \eta} & \frac{\partial b_1}{\partial \eta} & \frac{\partial b_2}{\partial \eta} & \frac{\partial b_3}{\partial \eta} \end{pmatrix} \begin{pmatrix} x_0 & y_0 \\ x_1 & y_1 \\ x_2 & y_2 \\ x_3 & y_3  \end{pmatrix}\:. The transform of the derivative of function \(f\) is {math} \begin{align} \frac{\partial f}{\partial\xi} &=  \frac{\partial f}{\partial x}\frac{\partial x}{\partial\xi} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial\xi}\:,\\ \frac{\partial f}{\partial\eta} &=  \frac{\partial f}{\partial x}\frac{\partial x}{\partial\eta} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial\eta}\:. \end{align} Which becomes {math} \begin{pmatrix} \frac{\partial f}{\partial\xi}\\ \frac{\partial f}{\partial\eta} \end{pmatrix} = \begin{pmatrix} \frac{\partial x}{\partial\xi} & \frac{\partial y}{\partial\xi}\\ \frac{\partial x}{\partial\eta} & \frac{\partial y}{\partial\eta} \end{pmatrix} \begin{pmatrix} \frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y} \end{pmatrix}\:, or {math} \begin{pmatrix} \frac{\partial}{\partial x}\\ \frac{\partial}{\partial y} \end{pmatrix} = \mathbf{J}^{-1} \begin{pmatrix} \frac{\partial}{\partial\xi}\\ \frac{\partial}{\partial\eta} \end{pmatrix}\:.

Numerically, we can integrate \([\mathbb{M}]\), $ [\(mathbb{M}\)]{ij} = :raw-latex:`sum`{p=0}^{P} b_i(\(xi\)_p,:raw-latex:eta_p)b_j(\(xi\)_p,:raw-latex:eta_p):raw-latex:det`(:raw-latex:mathbf{J}`_p)w(\(xi\)_p)w(\(eta\)_p):, $ where \(P+1\) is the number of quadrature points, generally set to be \(O(h^{2k+2})\) accurate where \(k\) is the polynomial order of the basis. So, for \(P1\) basis functions (first-order) we need to set the quadrature to integrate \(p=2k+1\) polynomials to achieve \(2k+2\) order accuracy. For example, Gauss-Legendre quadrature will integrate \(p=2n-1\), so \(n=2\) will integrate up to 3rd order polynomials exactly and will be \(O(h^{2n})\) accurate.

Expanding gradient term in the \([\mathbb{T}]\) Matrix, $ [\(mathbb{T}\)]{ij} = :raw-latex:`mu`:raw-latex:`int`{y} \(int\){x}:raw-latex:`frac{partial b_i}{partial x}`b_j,dxdy + :raw-latex:`eta`:raw-latex:`int`{y} \(int\)_{x}:raw-latex:frac{partial b_i}{partial y}`b_j,dxdy:. $ Numerically, we integrate :math:`[mathbb{T}] using {math} \begin{multline}     [\mathbb{T}]_{ij} = \sum_{p=0}^P     \Biggl(\mu\biggl([\mathbf{J}^{-1}]_{00}\frac{\partial b_i(\xi_p,\eta_p)}{\partial\xi}     + [\mathbf{J}^{-1}]_{01}\frac{\partial b_i(\xi_p,\eta_p)}{\partial\eta}\biggr)\\     + \eta\biggl([\mathbf{J}^{-1}]_{10}\frac{\partial b_i(\xi_p,\eta_p)}{\partial\xi}     + [\mathbf{J}^{-1}]_{11}\frac{\partial b_i(\xi_p,\eta_p)}{\partial\eta}\biggr)\Biggr)     b_j(\xi_p,\eta_p)\det(\mathbf{J}_p)w(\xi_p)w(\eta_p)\:. \end{multline} The leftmost term in the DG-FEM equation represents streaming inflow/outflow in a cell. Expanding as a sum over faces we define {math} \Gamma(\psi) = \oint_{\partial V_e} b_i\hat{\mathbf{n}}\cdot\Omega\psi\,dA =  \sum_{\mathrm{f}} \Omega\cdot\hat{\mathbf{n}}_\mathrm{f}\int_{A_\mathrm{f}} b_i \psi\,dA\:,\quad i=0,\cdots,B To calculate \(\Gamma\) we define $ \(Gamma`(:raw-latex:\)psi`) = \(mathbb{G}`^- + :raw-latex:\)mathbb{G}`^+:, $ where the \(\pm\) superscripts indicate incoming and outgoing terms, respectively, and \(\mathbb{G}^{-} = (N\times 1)\) and \(\mathbb{G}^{+}=(N\times N)\). Now for all incoming faces we have {math} \begin{align} \mathbb{G}^- &= \sum_{\mathrm{f}} \Omega\cdot\hat{\mathbf{n}}_\mathrm{f} \sum_{p=0}^{F} b_i(\beta_{p})\Bigl( b_m(\beta_{p})\psi_{mn}^{(m)} + b_n(\beta_{p})\psi_{mn}^{(n)} \Bigr) \det(\mathbf{J}_{p}) w(\beta_{p}) \:,\\ \mathrm{f}&\in\{\mathrm{f}_{mn} |\ \forall \Omega\cdot\hat{\mathbf{n}}_\mathrm{f} < 0\}\:,\\ i&=0,\ldots,B\:, \end{align} where \(\psi_{mn}^{(m/n)}\) are the upwind fluxes on face \(\mathrm{f}\) at node \(m\) or \(n\), and \(F+1\) is the number of quadrature points used to integrate face \(\mathrm{f}\). For the outgoing faces, the terms are {math} \begin{align} \left[\mathbb{G}^{+}\right]_{ij} &= \sum_{\mathrm{f}} \Omega\cdot\hat{\mathbf{n}}_\mathrm{f} \sum_{p=0}^{F} b_i(\beta_{p})b_j(\beta_p) \det(\mathbf{J}_p)w(\beta_{p})\:,\\ \mathrm{f}&\in\{\mathrm{f}_{mn} |\ \forall \Omega\cdot\hat{\mathbf{n}}_\mathrm{f} > 0\}\:,\\ j&=m,n\:. \end{align} And, the \(\beta_{p}\) are the quadrature points on face \(\mathrm{f}\). Each edge, \(\mathbf{f}_{mn}\), of the 2D quadrilateral is a vector defined by {math} \begin{align} \mathbf{f}_{01} &= \mathbf{f}_{23} = \biggl(\frac{\partial x}{\partial\xi}, \frac{\partial y}{\partial\xi}\biggr) = \mathbf{J}_\xi\:,\\ \mathbf{f}_{12} &= \mathbf{f}_{03} = \biggl(\frac{\partial x}{\partial\eta},\frac{\partial y}{\partial\eta}\biggr) = \mathbf{J}_\eta\:, \end{align} and \(\mathbf{J} = \begin{pmatrix}\mathbf{J}_\xi & \mathbf{J}_\eta\end{pmatrix}^T\).

Integrating along the edge in natural coordinates requires a transformation of the differential line segment, {math} \mathrm{f} = |\mathbf{J}_{\beta}|d\beta = \sqrt{\biggl(\frac{\partial x}{\partial\beta}\biggr)^2 + \biggl(\frac{\partial y}{\partial\beta}\biggr)^2}d\beta\:,\quad \beta = \xi,\eta\:. The coordinate transformations are defined {math} \det(\mathbf{J}_p)=|\mathbf{J}_p| = \begin{cases} |\mathbf{J}_\xi(\beta_{p})| & \mathrm{f}\in\{\mathrm{f}_{01},\mathrm{f}_{23}\}\:,\\ |\mathbf{J}_\eta(\beta_{p})| & \mathrm{f}\in\{\mathrm{f}_{12},\mathrm{f}_{03}\}\:. \end{cases}

Finally, we note that the normals are defined by rotating the edge vectors through \(\theta = \pi/2\) at each edge center point in natural coordinates, {math} \hat{\mathbf{n}}_\mathrm{f} = \begin{pmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{pmatrix}\mathbf{J}_\beta\:,\quad \beta=\xi,\eta\:.

Thecomplete system that needs to be solved in each cell after upwinding is {math} (\mathbb{G}^+-\mathbb{T} + \sigma\mathbb{M})\Psi = \mathbb{M}S-\mathbb{G}^-\:.

Reference Discretizations

Here we give reference DG-FEM equations for a regular, Cartesian grid. On a regular grid, the basis functions can be integrated analytically. For example, in 2D the basis functions are: {math} \begin{align} b_0 &= \Bigl(\frac{x_R-x}{\Delta_x}\Bigr)     \Bigl(\frac{y_T-y}{\Delta_y}\Bigr)\:,\\ b_1 &= \Bigl(\frac{x-x_L}{\Delta_x}\Bigr)     \Bigl(\frac{y_T-y}{\Delta_y}\Bigr)\:,\\ b_2 &= \Bigl(\frac{x-x_L}{\Delta_x}\Bigr)     \Bigl(\frac{y-y_B}{\Delta_y}\Bigr)\:,\\ b_3 &= \Bigl(\frac{x_R-x}{\Delta_x}\Bigr)     \Bigl(\frac{y-y_B}{\Delta_y}\Bigr)\:. \end{align} and {math} b_n(\mathbf{r}) = \begin{cases} 1 & \mathbf{r} = (x_n,y_n)\\ 0 & \mathbf{r} = (x_k,y_k)\ \forall\ k\ne n \end{cases}\:.

{math} \begin{align} \Gamma(\psi) &= \frac{\Delta_{x} \eta}{6} \begin{pmatrix}- 2 \psi_{01}^{(0)} - \psi_{01}^{(1)}\\- \psi_{01}^{(0)} - 2 \psi_{01}^{(1)}\\ 2 \psi_{23}^{(2)} + \psi_{23}^{(3)}\\\psi_{23}^{(2)} + 2 \psi_{23}^{(3)} \end{pmatrix} +  \frac{\Delta_{y} \mu}{6} \begin{pmatrix} - 2 \psi_{03}^{(0)} - \psi_{03}^{(3)}\\ 2 \psi_{12}^{(1)} + \psi_{12}^{(2)}\\ \psi_{12}^{(1)} + 2 \psi_{12}^{(2)}\\ - \psi_{03}^{(0)} - 2 \psi_{03}^{(3)} \end{pmatrix}\:,\\ \mathbb{T} &= \frac{\Delta_{y} \mu}{12} \begin{pmatrix}-2 & -2 & -1 & -1\\2 & 2 & 1 & 1\\1 & 1 & 2 & 2\\-1 & -1 & -2 & -2\end{pmatrix}+ \frac{\Delta_{x} \eta}{12} \begin{pmatrix}-2 & -1 & -1 & -2\\-1 & -2 & -2 & -1\\1 & 2 & 2 & 1\\2 & 1 & 1 & 2\end{pmatrix}\:,\\ \mathbb{M} &= \frac{\Delta_{x} \Delta_{y}}{36} \begin{pmatrix}4 & 2 & 1 & 2\\2 & 4 & 2 & 1\\1 & 2 & 4 & 2\\2 & 1 & 2 & 4\end{pmatrix}\:. \end{align}