Expected Value Estimator Anlaysis

[sslattery]

Okten's expected value estimator for his 2005 paper has been added to MCLS in addition to Gelbard's collision estimator which has been used to this point. Based on the literature and my parallel implementation, I expect the new estimator to significantly outperform the collision estimator. This new estimator does not solve the problem of long stochastic histories, but it should significantly reduce the number of histories need for convergence; with or without MCSA.

The purpose of this task is to determine whether or not the above statements about the new estimator are true within the context of the SPN challenge problems I am working on. This is done when I have run a series of problems that compares iteration and timing performance for MCSA with the two estimators.

[sslattery]

A couple of upates on this:

1. Okten's 2005 paper has an error in the expected value estimator formulation. He would likely have never discovered this due to the test problems used in his paper. It took me a few not so fun hours to find it. I have corrected this error and will describe later.
2. The expected value implementation fits right into the parallel framework used for the MCLS collision estimator at the cost of effectively storing an extra matrix to implement the tallies. It should not effect the scaling of the algorithm.
3. The memory costs are complemented by an improvement in the iterative performance. For the 1-group SP1 problem 3 system preconditioned with ILUT, MCSA with the collision estimator converged in ~35 iterations per eigenvalue iteration while MCSA with the expected value estimator converged in ~10 iterations per eigenvalue iteration. Both runs used the same ILUT and MCSA parameters. As a reference, Aztec BiCGSTAB preconditioned with the same ILUT parameters converged in ~17 iterations per eigenvalue iteration.
4. It looks like I can relax the number of histories per iteration and the weight cutoff when using the expected value estimator which will not only shorten the Markov chains but also reduce the number required for convergence. This may bring Jacobi preconditioning back in to play.
5. MCSA with the different estimators responds differently to the Neumann relaxation parameter. For the collision estimator, an optimal relaxation parameter was around 0.9 for problem 3 while the expected value estimator performed best with a value of around 1.0. This makes sense because the expected value estimator is giving a better correction for the same number of histories and therefore the Richardson extrapolation from which the Monte Carlo method is derived can be pushed a little harder.

More detailed numerical studies will be tomorrow. Some code refactoring will be completed as well.

[sslattery]

Just wanted to give an update here that there was not actually an error in Okten's paper but in fact a simple error in my implementation of the history weight calculation. The above results are still valid with this correction as my modification to Okten's estimator was purely to overcome my programming mistake.

[sslattery]

Here are some nice results from this comparison. Using a 1-group SP1 problem 3 with 20,088 DOFs, ILUT preconditioning was applied with a drop tolerance of 1.0e-3 and a fill level of 2.5. The independent variable in this test was the number of stochastic histories per MCSA iteration.

This is what the takeaways are from this data:

1. MCSA always converges for this problem and preconditioning with any number of stochastic histories using the expected value estimator, all the way down to 1. This makes sense as at a very minimum, even with no histories or the statistical noise from a few, we are always at least getting the source vector back out for the correction, effectively having done a single Richardson iteration.
2. MCSA does not always converge with the collision estimator. If too few histories are used loss of convergence is observed as the stochastic error is too large relative to the size of the correction.
3. In general with these ILUT values, strict monotonicity of the residual is not maintained although when lost it was not by much. I expect a higher fill level or smaller drop tolerance would alleviate this.
4. The dense composite operator means that the expected value estimator is really slow. At each row reached in the system, the expected value estimator tallies into all columns in that row. There are lots of columns if the matrix is dense. This gives more information and therefore the better convergence over the collision estimator that only tallies into the row state.

Overall, if the density of the iteration matrix can be reduced the expected value estimator is clear winner. This will reduce the tally time and significantly improve robustness by guaranteeing a properly conditioned problem will always converge regardless of the number of histories chosen. Next I'm going to explore how these estimators respond to the relaxation parameters.

[sslattery]

The evidence is growing against using the collision estimator for the SP1 1-group problem. Here are the residual convergence plots for both collision estimators with 5,000 histories per MCSA iteration:

The difference is even clearer at 1,000 histories per MCSA iteration:

Robustness is a serious issue here. On top of that, we get a major reduction in the amount of Monte Carlo transport required to get a good solution.

[sslattery]

This should do it for the collision estimator. Here's the same problem with SP3 discretization and 1-group:

It takes even more histories to get convergence out of the collision estimator while performance is about the same. As a final reason to only use the expected value estimator for SPN problems, here's a SP1 2-group problem with a fast and thermal group:

The collision estimator flat out won't converge with a reasonable number of histories. Given the results in this issue, the superior convergence results of the expected value estimator have been proven for the SPN equations. For all future SPN work, I will be using this estimator. Obviously, its performance as well as that of the rest of the algorithm is dependent on the sparsity of the preconditioned domain.

I am putting this issue up for review by Paul.

[sslattery]

Based on these results, I'll be making another issue that looks at the Richardson and Neumann relaxation parameters only for the expected value estimator.