Spectral Analysis of SPn Equations

[sslattery]

A spectral analysis of the SPn equations for multi-group, higher-order expansions are required. This spectra will give us an idea of convergence, preconditioning requirements, and permit a hypothesis for performance vs. subspace methods.

[sslattery]

I'm going to start today by looking at the iteration matrix spectra generated by single and multigroup P1, P3, P5, and P7 equations to get a feel for their structure. Analytically, I may see if I can perform the same spectral analysis on the general SPn equations that I did for the diffusion equation (although the group coupling may make that more difficult).

[sslattery]

Here's the eigenvalues of the Jacobi-preconditioned iteration matrix for an SPn7, P3, 3-group problem with upscatter. Note that all eigenvalues are on the real number line (with roundoff error) but Jacobi preconditioning DOES NOT reduce the spectral radius to less than one.

[sslattery]

Heres the eigenvalues for the same problem with SPn5 instead of SPn7. The largest positive eigenvalue is the same for both cases (0.9945), but reducing the SPn order brings the spectral radius down to 1.024. I've also double checked the matlab code used to generate these plots and it looks correct (the resulting iteration matrix appears correctly scaled with 0's on the diagonal).

[sslattery]

Finally, heres the problem with SPn3 with the spectral radius now < 1:

This appears to be related to the symmetry of the problem. For the symmetric (a single energy group), Jacobi preconditioning should always be enough for these equations to reduce the spectral radius ( a function of the diagonal dominance of the matrix). The SPn3 equations with 3 groups are fairly symmetric as shown by this zoom-in of the diagonal region of a sparsity plot for this iteration matrix:

And so the symmetric component of this matrix system is responding to the scaling by always reducing the largest positive eigenvalue to less than one for the 3 cases given here while not always doing so for the negative eigenvalues. This makes sense as the M&C paper analysis I did gave positive-real eigenvalues for the diffusion equation (to which these SPn equations ultimately reduce) and Jacobi preconditioning was enough in this case to reduce spectral radius for the iteration matrix. So system parameters like mesh size, cross sections, etc. that are present in the symmetric problem seem to dictate the positive eigenvalues as they were constants in all of these problems. The SPn order in this case was primarily affecting the negative eigenvalues produced.

[sslattery]

One last piece of data here for discussion. I generated the Gauss-Seidel iteration matrix for the first case above (SPn7, P3, 3-groups). The eigenvalues of this iteration are plotted here:

This is exactly what we expect here with the eigenvalues in a well-defined domain of convergence inside of the unit circle, giving an absolutely convergent stationary method. There is a paper that implements a Gauss-Seidel type of splitting for Monte Carlo methods:

A. Srinivasan, "Monte Carlo linear solvers with non-diagonal splitting", Mathematics and Computers in Simulation, Vol 80, 2010, pp. 1133-1143.

I never though I'd need to implement it, but we may have to. The only other option I see right now is to generate a spectral mapping as I've discussed previously. This shouldn't be too tough from a numerical analysis perspective because we've shown above that for Jacobi splitting the eigenvalues are real. This gives some simple functional forms that we can use for the spectral mapping.

[tmdelellis]

Stuart,

Just FYI, what does a 1 group SP7 spectral plot look like?

Tom

On Tue, Feb 26, 2013 at 3:05 PM, Stuart Slattery notifications@github.comwrote:

One last piece of data here for discussion. I generated the Gauss-Seidel iteration matrix for the first case above (SPn7, P3, 3-groups). The eigenvalues of this iteration are plotted here:

[image: gsSPn7P3G3]https://f.cloud.github.com/assets/755901/198091/f935a730-804e-11e2-810f-7de4b2f58adb.png

This is exactly what we expect here with the eigenvalues in a well-defined domain of convergence inside of the unit circle, giving an absolutely convergent stationary method. There is a paper that implements a Gauss-Seidel type of splitting for Monte Carlo methods:

A. Srinivasan, "Monte Carlo linear solvers with non-diagonal splitting", Mathematics and Computers in Simulation, Vol 80, 2010, pp. 1133-1143.

I never though I'd need to implement it, but we may have to. The only other option I see right now is to generate a spectral mapping as I've discussed previously. This shouldn't be too tough from a numerical analysis perspective because we've shown above that for Jacobi splitting the eigenvalues are real. This gives some simple functional forms that we can use for the spectral mapping.

— Reply to this email directly or view it on GitHubhttps://github.com/sslattery/Chimera/issues/25#issuecomment-14136596 .

[sslattery]

Here's a Jacobi iteration matrix Eigenvalue plot for SPn7, P3, 1-group:

I'll have to do some more of these to find out whether the legendre expansions for angle or scattering or the multigroup formulation have the greatest effect. So far, it looks like the SPn order has the greatest impact.

[sslattery]

I ran the iteration matrix from Srinivasan's method in the reference I inserted above and found it was not better than Jacobi. In fact, it asymptotically performs equivalently to the Jacobi method in his paper, so we don't expect any better results for the eigenvalue scaling. Check that one off the list as a no.

[sslattery]

Per the fixes to the SPN equations in Denovo I reran a parameter study for the spectral radii generated by Jacobi preconditioning. The SPN order, PN order, number of energy groups, and upscatter/downscatter were varied to yield an array of matrices (some of them very large). The spectral radius results are tabulated below. The red results required a relaxed eigenvalue solver in tolerance compared to the 1.0e-8 used for the rest of the computations.

As you can see, the SPN order has the greatest effect here and for all cases SPN=5,7 will not converge with Jacobi preconditioning.

Next, I tried block Jacobi preconditioning (which is still pretty easy to implement with a dense matrix inversion function). The block size is N_g * (N+1)/2 where N_g is the number of energy groups and N is the SPN order. This correlates to the block sizes observed in the formulation of the equations and the actual sparse matrix constructed in matlab. Here are the spectral radius results for the same set of matrices:

This works really well! In fact, these spectral radii are so small that the previous work I did for the M&C paper tells me that the MC solver should be really quick, and likely competitive with preconditioned Krylov methods. The next step is to decide whether to implement the preconditioning in Denovo or MCLS. I vote internal point and block preconditioning in MCLS as it is easy to implement and I can directly access the LAPACK routines to invert the dense blocks through Trilinos. That way, if block Jacobi is good enough for Navier Stokes (which it just might be considering these results) then that problem becomes plug and play.

[sslattery]

I'm going to put this ticket into review as this was everything I was looking for.