A sweep-based low-rank method for the discrete ordinate transport equation

We define the orthonormal bases $X_{i}$ and $V_{j}$ as

where $Z_{m}(x)=\frac{1}{\sqrt{\Delta x}}$ with $x\in[x_{m-\frac{1}{2}},x_{m+\frac{1}{2}}]$ are based on a finite volume discretization in space with a constant mesh spacing $\Delta x$ and $M$ zones; $m$ is the cell number. Here $u^{k}_{mi}$ are components of the time dependent matrix $U^{k}(t)\in\mathbb{R}^{M\times r}$, and $v^{k}_{*,j}$ refers to the $j_{th}$ column of the $V^{k}(t)\in\mathbb{R}^{N\times r}$. Note that the ${N}\times 1$ vector $v^{k}_{*,j}$ is a group of linear factors of all positive or negative ordinates, so there are ${r}$ combinations for each set of angles.

Next we describe procedures to solve Eq. (24) using source iteration and transport sweeps. By applying the standard upwinding technique, we write the spatial derivative term using upwind finite differences as

Note that $V^{k}_{\ell}\bm{\mu}^{k}V^{k}_{j}$ forms a matrix $R^{k}=(r^{k}_{\ell j})\in\mathbb{R}^{r\times r}$ that is precomputed. Thus we can discretize Eq. (24) using backward Euler method for the time integration

where $(s^{k}_{ij})^{\pm}=\frac{1}{2}\sigma_{\mathrm{s}}\phi V^{\pm}+\frac{1}{2}QV^{\pm}+\frac{1}{h}K^{\pm,0}\in\mathbb{R}^{M\times r}$ is the source term and the superscript 1 denotes the current time step and 0 denotes the previous time step. For maximal clarity we write out Eq. (32) in full, matrix form:

Here $L^{+}_{i}=\frac{1}{\Delta x}R^{+}+\left(\sigma_{t}(i)+\frac{1}{h}\right)I_{r}$, $L^{-}_{i}=-\frac{1}{\Delta x}R^{-}+\left(\sigma_{t}(i)+\frac{1}{h}\right)I_{r}$, $I_{r}$ is the $r\times r$ identity matrix, $\bar{R}^{k}=\frac{1}{\Delta x}R^{k}$, $\sigma_{t}(i)$ denotes the value of total cross-section in cell $i$, $b^{k}_{L}$ and $b^{k}_{R}$ are the left and the right boundary values. We set $b^{k}_{L}=b^{k}_{R}=0$ to correspond to vacuum boundary conditions. A compact form of the marching scheme for Eq.(31) is

Eqs. (25) and (26) are matrix ordinary differential equations and can be solved by an implicit Runge–Kutta method.

This work focuses on reducing the computer memory cost for radiative transfer simulations. The memory consumption for solving Eq. (3) using the classical iteration solution procedure is $2\times M\times(N_{\Omega}+1)$ during each time step to store the angular and scalar intensities. Here we do not include the simulation parameters such as the scattering/absorption cross-sections and the prescribed source based on the assumption that they can be stored efficiently during the computation.

In the low-rank algorithm the memory requirements are determined by the size of matrices $U$, $S$ and $V$. The initial conditions $U^{k,0}$, $V^{k,0}$, and $S^{k,0}$ take $(M\times r+N\times r+r^{2})\times n_{o}$ floating point numbers. Solving Eq. (24) requires $M\times r\times n_{o}+M$ more memory to store the scattering source and the updated matrix $U$ (factorized from the updated matrix $K$). Eq.(25) is solved next and uses $N\times r\times n_{o}$ memory for the updated $V$ for a total of $(M\times r+N\times r+r^{2})\times n_{o}$ for the updated basis. Lastly, the $n_{o}$ equations in Eq. (26) are $r\times r$ Sylvester matrix equations which require $n_{o}\times\mathcal{O}(r^{3})$ of memory when using the algorithm proposed in [24]. In this step the memory usage is $(M\times r+N\times r+r^{2})\times n_{o}+n_{o}\times\mathcal{O}(r^{3})$. In practice, we choose the rank $r\ll\min(M,N)$, so we can assume $\mathcal{O}(r^{3})+r^{2}<M\times r+N\times r$. To summarize we use $2\times(M\times r\times n_{o}+N\times r\times n_{o}+M)$ to approximate the largest memory usage during the low-rank calculations. We can see that the low-rank solution requires much smaller memory than the full solution, when $r\ll\min(M,N)$.

In numerical experiments we measure the low-rank memory occupied by double precision floating point numbers in megabytes (MB) as

and for the full-rank solution the memory required is

Usually, the initial conditions $U^{k,0}$, $V^{k,0}$, and $S^{k,0}$ for the low-rank evolution are obtained by taking the SVD of the given initial intensity and then truncating the singular values smaller than the $r$th largest in matrix $S$ and removing the corresponding rows and columns in $U$ and $V$. Nevertheless, setting the initial condition is not always straightforward. If the intensities are zero or a constant initially (as in many common test problems), the initial condition will be of rank 0 or 1. One way to deal with this issue is to add $r$ randomly generated orthogonal basis to $U$ and $V$ and set $S$ to a $r\times r$ matrix of zeros [12]. We find that this approach does not work well for implicit schemes with large time steps. In this work, we choose to let the low-rank basis evolve for a very small time step, such as CFL$=0.01$, to obtain initial conditions with a basis that is in the range of the low-rank operator. After this treatment, we can set time steps as intended.

To demonstrate the accuracy of our low-rank algorithm, we consider a multi-material 2-D problem with an asymmetric layout originally presented in [14]. As detailed in Figure 2, this problem has purely scattering zones with $\sigma_{\mathrm{t}}=\sigma_{\mathrm{s}}=0.01\,\mathrm{cm}^{-1}$, absorbing zones $\sigma_{\mathrm{t}}=100\,\mathrm{cm}^{-1},\sigma_{\mathrm{s}}=0.1\,\mathrm{cm}^{-1}$, and an isotropic source at the bottom with $Q=1$ of thickness 0.1 $\mathrm{mm}$. We solve this problem using the spatial grid of size $90\times 90$ for the computational domain $[0,9\,\mathrm{mm}]\times[0,9\,\mathrm{mm}]$, and the simulation time $t=0.9$ s. The CFL number is set to be 3.

Figure 3 compares the low-rank solutions for the scalar intensity with S₃₂ and varying rank to the full rank S₃₂ solution and the S₁₀₀ benchmark solution. As we can see S₃₂ with rank 36 is close to the full rank S₃₂ solution but it is not enough to resolve the particle distribution behind the second obstacle and the negative results are observed. The low-rank solution with rank 25 is nearly identical to the full rank S₃₂ solution, which means rank 25 is sufficient to capture the important features in this problem.

Figure 4 presents a comparison for scalar intensities along $x=0.45$ cm and $x=0.6$ cm. The left plot emphasizes their different behavior behind the second chevron, and the right plot shows the particle distribution along the gap that particles can travel through. We find that the low-rank solution with rank 36 has a similar error to the full rank S₃₂ solution when compared to the benchmark. We also notice that S₃₂ is not enough to capture the particle stream along the gap and it requires high angular resolution to obtain.

The computational efficiency of our low-rank algorithm is presented in Figure 5. Our solutions’ accuracy is measured by the root mean square (RMS) error to the benchmark solution. In Figure 5(a) we plot the error of our solutions versus the running memory in the sweeping process. As we can see, the third point in the red dotted line, which indicates the solution with rank 16, has already achieved the same error level as the full rank S₃₂ solution, but with only $10\%$ of the memory. It also can be seen that the low-rank and the full rank solution use a similar amount of memory for the same rank. As explained by our memory formula (33), the memory usage depends on the rank when $M$ is fixed, and $N$ is much smaller than $M$.

Figure 5(b) makes the similar comparison with the running time per one time step ($\Delta t=0.03$ s for this test). It appears that our rank 16 solution could save $95\%$ of the computational time compared to the full rank solution. We also point out that our low-rank method will take more time than the full rank solution with the same rank. For example, the rightmost dot in the red line is the solution with rank 100, the same as the full rank S₂₀ solution, but it runs twice as long. But when rank is properly chosen and relatively small, we can still save a large amount of memory and time.

The computer memory occupied in MATLAB during the calculation and the theoretical values are shown in Figure 6. The coefficients of determination ($R^{2}$) for the linear fit are approximately one, which indicates a strong linear relationship between the theoretical and actual memory. We conclude that Eqs. (33) and (34) are valid estimations for the actual memory usage, and the low-rank algorithm reduces the memory requirement proportionally to the reduction in rank.

Next, we consider the lattice problem: a $7\times 7$ cm checkerboard consisting of pure absorbers, purely scattering regions, and a strong source turned on at $t=0$ with a zero initial condition. We use a computational domain of $[0,7]\times[0,7]$ with a $210$ spatial grid. We solve this problem with a single, large time step using CFL $=10^{4}$. The layout and the reference solution calculated with S₆₄ for this problem are shown in Figure 7.

This test aims to show the accuracy improvement by adding more angular directions while keeping the rank fixed. As shown in Figure 8, the full rank S₆ solution has ray effects in the solution near the bottom, right and left sides. These can be alleviated by using more discrete ordinates as observed in solutions with S₆₄ and rank 9. Similar phenomena are found in solutions with rank 16, where the full rank solution suffered from ray effects while the low-rank solution is closer to our reference. We also plot the solution along $x=3.5$ cm, where we can see that low-rank solutions are closer to the benchmark solution than the full rank solution.

Figure 9 gives a quantitative comparison between low-rank and full rank solutions in terms of their memory usage and error measured by RMS deviations to the reference. We notice that the low-rank solutions converge faster than the full rank solutions, where the low-rank solution with rank 49 achieves the same accuracy as the full rank S₅₀ solution, but only requires $10\%$ of the memory.

The line source problem describes a beam of radiation is spreading out in a purely scattering plane. In this problem we have the initial condition $\psi(x,z,\mu,\varphi,t)=\delta(x)\delta(z)$, and the purely scattering medium $\sigma_{\mathrm{s}}=\sigma_{\mathrm{t}}=1$ with no source $Q=0$. We use a computational domain of $[-1.5,1.5]\times[-1.5,1.5]$ for the simulation time $t=1\,$s, while the spatial grid is set to be $150\times 150$. Figure 10 shows the analytic solution and our benchmark solution.

We use the line source test to demonstrate the benefits of the low-rank method with high angular resolutions. Figure 11 compares the full rank solution to low-rank solutions with the same rank. We notice that all the full-rank solutions have remarkable ray effects, which means the number of angular directions are insufficiently dense to move particles to all the regions . The low-rank method could significantly improve the full rank solution by using a small amount of extra memory. As we can see, the ray effects are alleviated in the solution with rank 16, and the solution with rank 36 is comparable to the benchmark solution.

The last 2-D problem is a modified Hohlraum problem as shown in Figure 12. The hohlraum is surrounded by vacuums, except for the incoming source on the left. In this problem, we are interested in the particle distribution behind the block that particles cannot reach directly, which will require more angular direction samples to resolve. For this problem we use a computational domain of $[0,1.3]\times[0,1.3]$ for the simulation time $t=0.9\,$s, and set the spatial grid to $150\times 150$.

We make a similar comparison between the low-rank solutions to the full rank solutions with the same rank in Figure 13. We observe beam-shaped shadows in the left part of the full rank S₈ solution, but not in the low-rank solution with S₁₀₀ and rank 16. There are negative regions in the low-rank solutions brought by the truncation from low-rank algorithms.

Figure 14 presents the cut along $z=0.65$ cm to compare the particle distribution behind the obstacle. The left plot shows that the solution with rank 36 is nearly identical to the full rank solution. The right plot shows that the low-rank solutions are closer to the benchmark than the full rank solution with the same rank.

As shown in Figure 15, we introduce the 3-D wires problem, inspired by the 2-D problem in [8], has an emitting dense material embedded vacuum. Particles are emitted from the central wire with $Q=1$. All the five wires are made by the same material with $\sigma_{\mathrm{s}}=0.1$ and $\sigma_{\mathrm{t}}=5$. The streaming regions has $\sigma_{\mathrm{s}}=\sigma_{\mathrm{t}}=0$ and we use vacuum boundary conditions. The computational domain for this problem is set to $[-0.01,0.01]\times[-0.01,0.01]\times[0,0.02]$ with a $40\times 40\times 40$ mesh grid. We simulate this problem to $t=0.6$s using a time step $\Delta t=0.2$ (CFL$=1200$).

The following results compare the low-rank solution and the full rank solution with the same rank. Figure 16 shows the solution in the plane of $x-y$ along $z=0.1$. First, we notice strong ray effects in the full rank S₄ solution, especially along $y=0$. The low-rank solution with S₆₄ and rank 4 looks better in this aspect. The full rank solution with S₈ achieves better distributions than S₄ but is still unable to avoid the ray effects along $y=0$. But we can see that the low-rank solution with rank 16 is close to the full rank solution. If we further increase the rank to 36, the low-rank solution is identical to the benchmark solution.

We also plot the solution in the y-z plane along $x=0$ shown in Figure 17. We can see that the solution with rank 4 is has large, qualitative errors, especially in the center source region. The low-rank solution significantly improves the situation with the errors greatly reduced.

The running memory to carry out these simulations in MATLAB is shown in Figure 18. The rank 64 solution with S₆₄, represented by the rightmost dot in the solid red line, uses 25% of the full rank $S_{32}$ solution but achieves the same accuracy. It also reaches five times lower error than the full rank $S_{16}$ solution with slightly more memory.