# 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=\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(\xi,\eta)x_i\:,\quad y = \sum_{i=0}^B b_i(\xi,\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} = \sum_{p=0}^{P} b_i(\xi_p,\eta_p)b_j(\xi_p,\eta_p)\det(\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} = \mu\int_{y}
\int_{x}\frac{\partial b_i}{\partial x}b_j\,dxdy +
\eta\int_{y}
\int_{x}\frac{\partial b_i}{\partial y}b_j\,dxdy\:.
$
Numerically, we integrate $[\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(\psi) = \mathbb{G}^- + \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}
```