Paper summary: A Linear Prolongating Coarse Mesh Finite Difference Acceleration of Discrete Ordinate Neutron Transport Calculation Based on Discontinuous Galerkin Finite Element Method

Step 4 is the crucial step which we developed in order to implement lp-CMFD for DGFEM, so we can go into more detail on how this was calculated. We will start by defining the fine mesh flux averaged over the coarse mesh geometry as $\phi_{M}$ ($M$ is the coarse mesh index) which is calculated from the fine mesh flux $\phi_{m}$ ($m$ is the fine mesh index)
$$
\phi_{M} = \text{volume weighted mean }(\phi_{m}) m \in M
$$
We then denote the coarse mesh flux solved using the finite difference method as $\Phi_{M}$. The coarse mesh difference is $\delta\Phi_{M} = \Phi_{M} – \phi_{M}$. We repeat this procedure for every coarse mesh in the domain. We can calculate the $\delta\Phi_{M}$ values at the vertex of the coarse mesh grid. On the 1D domain the vertices of the coarse mesh grid are also the boundaries between neighboring meshes. The update value between coarse mesh $M$ and $M+1$ is the average between coarse mesh $M$ and $M+1$.
$$
\delta\Phi_{M+1/2} = \frac{1}{2}(\delta\Phi_{M} + \delta\Phi_{M+1})
$$
The vertex boundary values for coarse mesh $M$ is $\delta\Phi_{M-1/2}$ on the left and $\delta\Phi_{M+1/2}$ on the right. We can therefore linearly interpolate on the interior of the coarse mesh
$$
\delta\phi(x) = \delta\Phi_{M-1/2}(1 – \frac{x}{\Delta_{M}}) + \delta\Phi_{M+1/2}\frac{x}{\Delta_{M}}
$$
where $\Delta_{M}$ is the length of the coarse mesh and $x$ is the spatial coordinate between $[0, \Delta_{M}]$ within the coarse mesh. We can employ the DGFEM weighted residual method and find the updates to the coefficients. We will define the update flux in fine mesh $m$ on the same set of basis functions which was used in the $S_{N}$ transport calculation i.e.
$$
\phi^{update}(x) = \sum_{p} s_{p} u_{p}(x),
$$
where $s_{p}$ are the update coefficients of the basis functions $u_{p}$ to be solved. The residual equation is defined as
$$
\delta\phi(x) – \phi^{update}(x) = 0
$$
We generate $P$ number of equations by multiplying the residual equation with each basis function $u_{p}$
$$
\int (\delta\phi(x) – \phi^{update}(x)) \cdot s_{p}u_{p}(x) dx = 0 \text{ for } x \in x_{m}
$$
We skip the mathematical manipulation and simply cast the $P$ number of equations in matrix form
$$
\boldsymbol{T}\vec{s} = \vec{v}
$$
and finally
$$
\vec{s} = \boldsymbol{T}^{-1}\vec{v}
$$
We finally update the coefficients in fine mesh $m$ as
$$
\vec{c}^{update}_{m} = \vec{c}^{S_{N}}_{m} + \vec{s}
$$