I have been trying to figure how to discretize the 1D steady state monoenergetic Boltzmann Transport Equation, with a Discrete Ordinates (Sn) in angle and Continuous Galerkin Finite Element Method (CG-FEM) in space. I am confused as to how one combines the spatial and angular discretizations in order to define the linear system.
This is how I have approached the problem:
Boltzmann Transport Equation (BTE) steady state 1D, 1 energy group
\begin{align}\tag{1} \big [{\mu} \frac{d}{dx} + \Sigma_t\big ] \psi({x}, {\mu}) = \int \Sigma_s \psi({x}, {\mu}) \,d{\mu'} + S \, , \end{align}
where $\mu$ is the angular domain, with $\mu = \cos\theta$ and $\theta$ being the angle with the x-axis. $\Sigma_t \text{ and } \Sigma_s$ the total and scatter macroscopic cross-sections (constant in this case) and $S$ is the source term.
Angular discretisation: Sn
According to Lewis, E.E. and Miller, W.F., Computational Methods of Neutron Transport, p.118 or from Roberts J.A., Direct Solution of the Discrete Ordinates Equations, (see https://www.researchgate.net/publication/303471187_Direct_Solution_of_the_Discrete_Ordinates_Equations), one can discretise the angular domain to hold for a number of distinct angles $\mu_n$, Eq. 1 then becomes
$$ {\mu_n} \frac{d}{dx}\psi_n(x) + \Sigma_t \psi_n(x) = \frac{\Sigma_s}{2} \large\sum_{n=1}^{N} w_n \psi_n(x) + \frac{S}{2} \,. \tag{2} $$
Spatial Discretisation: CG-FEM
From what I understand now $\psi(x, \mu)$ has now been discretised to N directions and now our flux in solely dependant on the spatial dimension x. Therefore, by getting the weak form equation and performing integration by parts on Eq. 2 and approximating the angular flux $\psi_n(x) = \sum_{j=1}^{NB}N_j\tilde{\psi}_{n,j}(x)$, where $N_j$ is our basis functions, we get:
\begin{align} \sum_{n=1}^{N} \sum_{j=1}^{NB} \tilde{\psi}_{n,j}(x) \Biggl[ \ \int_e \big ( N_i N_j \Sigma_t - \mu_n \frac{dN_i}{dx} N_j - \frac{\Sigma_s}{2} w_n N_j\big ) dx+ \int_\Gamma \mu_n N_i N_j d\Gamma \Biggr]= \ \int_e \frac{S}{2}N_i dx \end{align} for all $i \in \{1,..., NB\}$.
Would that therefore mean, that if we used linear interpolation functions for our elements ($NB=2$) and two angular directions $N=2$ we would end up with the following linear system for an element $e$ that is not at the boundaries?
\begin{align} \left[ \ \Sigma_t h_e \left[ \begin{matrix} 1/3 & 1/6 \\ 1/6 & 1/3 \\ \end{matrix} \right] \ + \mu_1 \left[ \begin{matrix} -1/2 & -1/2 \\ 1/2 & 1/2 \\ \end{matrix} \right] \ + \Sigma_s w_1 \left( \begin{matrix} 1/4 \\ 1/4 \\ \end{matrix} \right) \ \right] \ \tilde{\psi}(x_e, \mu_1) = \ S h_e \left( \begin{matrix} 1/4 \\ 1/4 \\ \end{matrix} \right) \ \\ \left[ \ \Sigma_t h_e \left[ \begin{matrix} 1/3 & 1/6 \\ 1/6 & 1/3 \\ \end{matrix} \right] \ + \mu_2 \left[ \begin{matrix} -1/2 & -1/2 \\ 1/2 & 1/2 \\ \end{matrix} \right] \ + \Sigma_s w_2 \left( \begin{matrix} 1/4 \\ 1/4 \\ \end{matrix} \right) \ \right] \ \tilde{\psi}(x_e, \mu_2) = \ S h_e\left( \begin{matrix} 1/4 \\ 1/4 \\ \end{matrix} \right) \ \end{align}
where $h_e$ is the length of element $e$.
The way I have incorporated the scatter term (the Sn discretisation), does not seem right to me, I think somehow I messed up the sum. Any help would be greatly appreciated.