Fast diagonalization for matrix-free inversion of the nodal even-parity operator in hybrid finite element SPN transport

The hybrid finite element method stores a dense inverse of the nodal even-parity operator for every node and energy group, and that storage limits the spatial and angular expansion orders, more tightly still on a GPU, whose memory is smaller than that of a host. This work introduces fast diagonalization for the nodal even-parity operator of the hybrid finite element SPN method: the stored dense inverse is replaced by the eigendecompositions of a planar, an axial and an angular factor, and the inverse is applied without being formed explicitly. Every even-parity moment of the SPN equations obeys an equation of the same diffusion form, with no mixed derivatives, so the angular coupling is the same in every spatial direction. On a node extruded along the axis, the operator is then an exact Kronecker sum of the three factors. The factorization merely recasts the computation, so the discretization and the solution are unchanged, and the nodes are not required to be identical to one another. The method is implemented for CPU and GPU hardware and compared with the dense inverse within one solver. On a pressurized water reactor core in SP3, the two paths agree to within round-off error, the storage falls from 3.76 GB to 821 kB, and the calculation is 4.15 times faster. On a 33-group fast reactor core, the dense inverse would occupy 24.9 GB and the dense path would need more memory than both the device and the host provide, whereas the factors occupy only 1.26 MB, so the solution can be obtained. The factorization thus reduces the storage of the nodal inverse by more than three orders of magnitude and brings multigroup core calculations at their required spatial expansion orders within the memory of a single GPU.