炉物理の研究 第74号 (2022年3月)

Comparative Review of Coarse Mesh Acceleration Methods and the Two-Level p-CMFD
Acceleration in the Whole-Core Transport Calculation
Nam Zin Cho
Professor Emeritus
Korea Advanced Institute of Science and Technology (KAIST)
291 Daehak-ro, Yuseong-gu, Daejeon, Korea 34141
nzcho@kaist.ac.kr
March 2022
KEY ACRONYMS
MOC : Method of Characteristics
CMR : Coarse Mesh Rebalance
GCMR : Generalized Coarse Mesh Rebalance
CMFD : Coarse Mesh Finite Difference
p-CMFD : Partial current-based Coarse Mesh Finite Difference
p-IMFD : Partial current-based Intermediate Mesh Finite Difference
PCDF : Partial Current Discontinuity Factor
TL(N)/p-CMFD : Two-Level p-CMFD with N global/local iterations
PI : Power Iteration
G-S : Gauss-Seidel
WS : Wielandt Shift
BiCGSTAB : Bi-Conjugate Gradient Stabilized
JFNK : Jacobian-Free Newton Krylov
INTRODUCTION
Due to the structural design characteristics in the current generation of nuclear reactor cores, theneutron transport theory methods are now called for in the neutronics design and analysis [1]. In
typical reactor types such as the thermal and fast reactors, the 2D (method of characteristics)/1D(discrete ordinates method) fusion method [2] is a novel approach and getting increasing attention[3,4]. Still, the 2D method of characteristics calculation takes most of the computing time and thus
requires good acceleration schemes [1].
The coarse mesh based acceleration methods [1–11] in neutron transport calculation are widelyused, because they are easily applied to the original transport calculation with various geometries.
However, they exhibit slow convergence or divergent behavior (in some of the methods) for
optically thick coarse mesh cells. This drawback limits the size of coarse mesh cells, incurring
significant computational burden per iteration. It is also very important to use converging schemesin the coupling of the neutronics code with the thermal-hydraulics code for feedback [12].

1炉物理の研究 第74号 (2022年3月)
炉物理の研究 第 74 号 (2022 年 3 月)

The present article revisits the partial current-based coarse mesh finite difference (p-CMFD)
method [7–10] that was developed at KAIST, and that is unconditionally stable and more stable
than other coarse mesh based acceleration methods. A recent extension with a two-level (two-grid)convergence speedup scheme [13–16] was made by using the Krylov subspace method. Numericalresults of its application to the standard OECD benchmark problem and an enlarged problem showthat the two-level scheme wrapped around by the Krylov subspace enhances the convergence rateof p-CMFD significantly, especially for optically thick coarse mesh cells (of assembly size).
COARSE MESH BASED ACCELERATION METHODS AND LIMITATIONS
Coarse mesh based acceleration methods usually consist of two parts: a high-order calculation anda low-order calculation. The high-order calculation employs transport methods (usually based on
the discrete ordinates method or the method of characteristics in multigroup transport equations)with fixed fission source. The low-order calculation uses the balance equation, in which the highorder calculation provides parameters over each coarse mesh cell. The low-order calculation gives
the multiplication factor (for the whole reactor) and the coarse mesh cell averaged scalar fluxes.
The coarse mesh cell averaged scalar fluxes are then “modulated” or prolongated to be used in thefission source of the high-order equation for the next iteration.
To study the efficiency of an acceleration method, Fourier convergence analysis was performedon coarse mesh rebalance (CMR) [5,10], coarse mesh finite difference (CMFD) [10], partial
current-based CMFD (p-CMFD) [7–10, 13], and generalized coarse mesh rebalance (GCMR) [11].
The results show that p-CMFD is always stable and more stable than the other methods. CMRshows divergent behavior for optically thin coarse mesh cells. CMFD, GCMR, and p-CMFD arevery efficient for optically thin coarse mesh cells but not efficient for optically thick coarse meshcells; See Fig. 1 (CMFD even diverges for optically thick coarse mesh cells, if D is used from thestandard diffusion theory ( D = 1/3 t )). Note that a recent work on CMFD with linearlyinterpolating prolongation lp-CMFD [17] stabilizes this divergence).

Fig. 1. Results of Fourier stability analysis of CMR, CMFD, and p-CMFD with SC and S16 forthe eigenvalue problem (from Ref. 13).

2

炉物理の研究 第74号 (2022年3月)

Ref. 18 adjusted (optimized) the diffusion coefficient in CMFD to improve the convergence rate.
The diffusion coefficient in p-CMFD can be also adjusted to optimize the spectral radius. Figure
2 shows the results in comparison. Even with the optimized D, CMFD and p-CMFD show slow
convergence in optically thick coarse mesh cells. This is a common property for all coarse meshbased acceleration methods; see Ref. 18 for CMFD, Ref. 11 for GCMR, Refs. 10 and 13 for pCMFD, and Ref. 17 for lp-CMFD. Therefore, the need arises to improve convergence for thickcoarse mesh cells.

Fig. 2. Spectral radii of the D optimized CMFD and p-CMFD. (See additional resultsin Ref. 13.)
TWO-LEVEL p-CMFD ACCELERATION
W-Cycle Scheme
The p-CMFD acceleration is a modification of the CMFD acceleration, in that p-CMFD is basedon the use of partial currents (instead of net currents used in CMFD), resulting in:i) p-CMFD is unconditionally stable,ii) p-CMFD provides additional information, that is, transport partial currents on theinterface of two coarse mesh cells, and thus,iii) the incoming partial current allows a fixed-source problem formulation for a coarse
mesh cell, whose solution provides in turn improved flux distribution and outgoing partial
current (to the neighboring coarse mesh cell).
Note that CMFD [6] uses one correction coefficient to preserve the surface net current. On theother hand, p-CMFD uses two correction coefficients to preserve the two surface partial currents,
respectively (The surface net current is automatically preserved by the two surface partial currents).
In p-CMFD, there is little incentive to deal with the optimized D (the improvement is only
marginal in thin coarse mesh cells, where it is already very fast), instead we focus on thick coarse
mesh cells, utilizing transport partial currents (see Fig. 2). This section describes a two-level spatialcoarse grid p-CMFD technique to speedup p-CMFD acceleration, particularly in optically thickcoarse mesh cells. To reduce the spectral radius further, we also introduce a global/local iterative3炉物理の研究 第74号 (2022年3月)

procedure to the low-level calculation in the two-level p-CMFD, as shown in the multi-gridliterature notation in Fig. 3 (W-cycle). In Fig. 3, a quarter pin forms an intermediate-mesh cell in
p-IMFD and each assembly forms a coarse-mesh cell in p-CMFD. Note that the p-IMFD rebalancecalculations are used in restriction and prolongation for the scalar flux and partial currents.

Fig. 3. Two-grid representation of the p-CMFD computational flow for N = 2 in W-cycle. (*withupdated transport incoming partial currents)

Figure 4 is the convergence result of global/local iterations in the two-level p-CMFD acceleration
for a 1D version OECD/NEA benchmark problem. As the number of global/local iterations
increases (N goes large), the spectral radius for optically thick coarse mesh cells reduces drastically,
approaching that of the fine mesh limit.
1.2

TL(1)/p-CMFD
TL(2)/p-CMFD

Spectral radius (ρ)

1.0

TL(10)/p-CMFD

0.8

TL(100)/p-CMFD
p-CMFD

0.6
0.4
0.2
0.0

0.01

0.10
1.00
10.00
100.00Coarse-mesh cell optical thickness (σtΔ)

Fig. 4. Spectral radius of global/local iteration in two-level p-CMFD acceleration (TL(N)/pCMFD: two-level p-CMFD with N global/local iterations).

4炉物理の研究 第74号 (2022年3月)

Krylov Solution Algorithms
In the W-cycle scheme in Fig. 3, Restriction I results in p-IMFD of intermediate-mesh cells andleads to the form,Aф = b,(1)that is a very large problem. It is then followed by the p-CMFD eigenvalue problem in coarsemesh cells of the form through Restriction II,
Mф = 1/keff Fф.

(2)

They are both whole-core problems. Eq. (1) is solved by the BiCGSTAB algorithm [19] and Eq.
(2) by the Jacobian-Free Newton-Krylov (JFNK) algorithm [20]. In Eq. (2), ф and keff areconsidered as unknowns simultaneously, rendering them as roots of a nonlinear equation system.
They are compared to the traditional methods of Gauss-Seidel (G-S) and power iteration (PI),respectively. The Wielandt shift (WS) acceleration of power iteration (PI) is to be avoided, because
it may render the resulting matrix near singular.
NUMERICAL RESULTS
3 x 3 Problem
The two-level p-CMFD is tested on a two-dimensional problem with seven-group cross section
data that are specified in the OECD/NEA C5G7 benchmark report [21], and is shown in Figure 5.

Fig. 5. Geometry of 2D C5G7 benchmark problem.
The computational conditions are as follows: The method of characteristics (MOC) is used, with40 mesh cells (5 ring divisions and 8 azimuthal sector divisions) per fuel pin cell and 64 mesh cells
5炉物理の研究 第74号 (2022年3月)(8 by 8 rectangular mesh cells) per reflector (pseudo) pin cell. An intermediate-mesh cell is aquarter of a fuel pin or a reflector pin. Table I shows the various results in comparison. Note thatthe number of transport sweeps (that is most time-consuming among the various phases of thecomputation) reduces significantly from No Acceleration to p-CMFD, and more significantly to
TL(N)/p-CMFD, as N increases. This effect outweighs well enough the increase in the number of
p-CMFD power (or JFNK) iterations.
Table Ⅰ. Results of 2D C5G7 (3 by 3) Benchmark Problem
Methods

No Acc.
a*p-CMFDa*a*

TL(1)/p-CMFD
b*

c*

a*

TL(4)/p-CMFD
b*

c*

keff
1.18658 1.18658 1.18658 1.18658 1.18658 1.18658 1.18658 1.18658
Number oftransport sweep
1017
473
28
28
27
18
18
18
iterations
Number of pCMFD Power
0
4316
436
436
217
723
723
386
(or JFNK)
Iterations
Transport
sweep
86.188 82.766 55.467 56.747 56.993
calculation time 3148.618 1468.528 87.509(sec)†
Whole-core p0
0
5.311
3.034
2.06
2.8
1.859
1.376
IMFD time (sec)
Local p-IMFD
0
0
4.558
2.834
2.162
8.819
6.602
4.859
time (sec)Whole-core p0
15.802
0.937
0.925
0.944
2.368
2.379
2.581
CMFD time (sec)
Total calculation
3149.01 1484.67 98.412 93.081 88.017 69.552 67.687
65.9
time (sec)
Speedup
1
2.12
32.00
33.83
35.78
45.28
46.52
47.78
Percentage oftransport sweep
99.99
98.91
88.92
92.59
94.03
79.75
83.84
86.48
time in total
calculation timea*: Gauss-Seidel +Power Iteration, b*: BiCGSTAB +Power Iteration, c*: BiCGSTAB +JFNK†: a single thread of a CPU (Intel i7-7700K)

5 x 5 Problem
Additional results on an enlarged OECD/NEA C5G7 benchmark problem [22] are available andincluded below in Fig. 6 and Table II. As expected, the speedup becomes larger as the problem
becomes larger.

6炉物理の研究 第74号 (2022年3月)

Fig. 6. Geometry of enlarged 2D C5G7 benchmark problem.
Table Ⅱ. Results of Enlarged 2D C5G7 (5 by 5) Benchmark Problem
Methods

No
Acc.
a*

p-CMFDa*

TL(1)/p-CMFD

a*

b*

c*

a*

TL(4)/p-CMFD
b*

c*

keff
1.23141 1.23141 1.23141 1.23141 1.23141 1.23141 1.23141 1.23141
Number of transport
2241
438
48
48
49
18
18
18
sweep iterations
Number of p-CMFD
Power
0
4584
895
893
646
1023
1023
742
(or JFNK) Iterations
Transport sweep
19874.19
3852.734 424.896 433.586 432.902 159.347 159.241 158.804
calculation time (sec)†
4
Whole-core p-IMFD
0
0
23.143 13.831
8.467
8.203
5.515
3.536
time (sec)
Local p-IMFD time
00
19.212 12.682
9.307
25.41
18.688 11.893(sec)
Whole-core p-CMFD
038.363
4.177
4.308
4.765
6.353
6.333
6.91
time (sec)
Total calculation time
19874.33 3891.21 471.662 464.529 455.566 199.425 189.887 181.254(sec)
Speedup
1
5.11
42.14
42.78
43.63
99.66
104.66 109.65
Percentage of transport
sweep time in total
100.0
99.01
90.08
93.34
95.03
79.90
83.86
87.61
calculation timea*: Gauss-Seidel +Power Iteration, b*: BiCGSTAB +Power Iteration, c*: BiCGSTAB +JFNK†: a single thread of a CPU (Intel i7-7700K)

7炉物理の研究 第74号 (2022年3月)

CONCLUDING REMARKSThis article presented a further speedup technique for p-CMFD acceleration in the whole-coretransport calculation, that is also effective for optically thick coarse mesh cells (of assembly size).
It is based on two-level spatial coarse grid p-CMFD with global/local iterations in the low-order
calculation. The results show that we obtain fast convergence (even with large coarse-mesh cells),
if we use a number of global/local iterations in TL(N)/p-CMFD with Krylov algorithms. This isowing to the availability of transport partial currents on the interface of two coarse mesh cells,
allowing rebalance calculation via p-IMFD in each coarse-mesh problem (that can be performedin parallel and thus can reduce the local p-IMFD time further). More importantly, the most timeconsuming transport sweep calculation footprint can be reduced further in a major way, if thetransport sweep operations are performed characteristics by characteristics in parallel.
As a future work in another direction, the two-level p-CMFD method can be applied to the fast
reactor analysis. Since the neutron mean free path is quite long in a fast reactor, the computational
cells of super hexagonal- and rhombic-stencils could be considered, following the procedures
proposed in Ref. 23, based on the HIRE-theoretic “multigroup” transport equations [24–26] withpartial current discontinuity factors (PCDFs).
ACKNOWLEDGMENTS
The author expresses special thanks to Professor T. Kitada of Osaka University, who invited me
to write this article for the Reactor Physics Division of the Atomic Energy Society of Japan (AESJ).
The author also expresses sincere thanks for special associations and friendships of many years toProfessor T. Takeda of Osaka University (now of Fukui University), Professors K. Kobayashi and
C. Pyeon of Kyoto University, Professor A. Yamamoto of Nagoya University, Professors H.
Sekimoto and T. Obara of Tokyo Institute of Technology, Dr. Y. Tahara of Mitsubishi Heavy
Industries, Ltd., and Dr. M. Tatsumi of Nuclear Engineering International.
The author also expresses thanks to Dr. Seungsu Yuk of Korea Atomic Energy Research Instituteand Dr. Seongdong Jang of Korea Advanced Institute of Science and Technology, who provided
valuable assistance during the preparation of this article.
REFERENCES
1. M. L. ADAMS and E. W. LARSEN, “Fast Iterative Methods for Discrete Ordinates Particle
Transport Calculations,” Prog. Nucl. Energy, 40, 3 (2002).
2. N. Z. CHO, G. S. LEE, and C. J. PARK, “Fusion of Method of Characteristics and Nodal
Method for a 3-D Whole-Core Transport Calculation,” Trans. Am. Nucl. Soc., 86, 322 (2002).
3. B. Collins, et. al., “Stability and Accuracy of 3D Neutron Transport Simulations Using the
2D/1D Method in MvPACT,” J. of Computational Physics,” 326, 612 (2016).
4. G. Zhang, A. Shieh, W. S. Yang, and Y. S. Jung, “Consistent pCMFD Acceleration Schemesof the Three-Dimensional Transport Code PROTEUS-MOC,” Nucl.Sci. Eng., 193, 828 (2019).

8炉物理の研究 第74号 (2022年3月)

5. G. R. CEFUS and E. W. LARSEN, “Stability Analysis of Coarse-Mesh Rebalance,” Nucl. Sci.
Eng., 105, 31 (1990).
6. K. S. SMITH and J. D. RHODES III, “Full-Core, 2-D, LWR Core Calculations with CASMO4E,” Proc. PHYSOR 2002, Seoul, Korea, October 7–10, 2002, American Nuclear Society (2002)(CD-ROM).
7. N. Z. CHO, G. S. LEE, and C. J. PARK, “On a New Acceleration Method for 3D Whole-CoreTransport Calculations,” Annual Meeting, Sasebo, Japan, March 27–29, 2003, Atomic Energy
Society of Japan (2003).
8. N. Z. CHO, G. S. LEE, and C. J. PARK, “Partial Current-Based CMFD Acceleration of the
2D/1D Fusion Method for 3D Whole-Core Transport Calculations,” Trans. Am. Nucl. Soc.,
88, 594 (2003).
9. N. Z. CHO, “The Partial Current-Based CMFD (p-CMFD) Method Revisited,” Trans. Kor.
Nucl. Soc., Gyeongju, Korea, October 25–26, 2012, Korean Nuclear Society (2012).
10. N. Z. CHO and C. J. PARK, “A Comparison of Coarse Mesh Rebalance and Coarse Mesh
Finite Difference Accelerations for the Neutron Transport Calculations,” Proc. M&C 2003,
Gatlinburg, Tennessee, April 6–11, 2003, American Nuclear Society (2003); see also KAIST
internal report NURAPT-2002-02, (with addendum on p-CMFD);
https://github.com/nzcho/Nurapt-Archives/blob/master/NurapT2002_2rev2.pdf
11. A. YAMAMOTO, “Generalized Coarse-Mesh Rebalance Method for Acceleration of Neutron
Transport Calculations,” Nucl. Sci. Eng., 151, 274 (2005).
12. B. KOCHUNAS, A. FITZGERALD, and E. LARSEN, “Fourier Analysis of Iteration Schemes
for k-Eigenvalue Transport Problems with Flux-Dependent Cross Sections,” J. of Comp. Phys.,345, 294 (2017).
13. S. YUK and N. Z. CHO, “Two-Level Convergence Speedup Schemes for p-CMFD
Acceleration in Neutron Transport Calculation,” Nucl. Sci. Eng., 188, 1 (2017); andCorrigendum, 188, 303 (2017).
14. N. Z. CHO, “An Overview of p-CMFD Acceleration and Its Applications to Reactor Physics
Transport Calculation,” Trans. Am. Nucl. Soc., 117, 1247-1250 (2017).
15. N. Z. CHO, “The Roles of Transport Partial Current Information in Two-Level p-CMFD
Acceleration in the Whole-Core Transport Calculation,” Trans. Am. Nucl. Soc., 121, 13231326 (2019).
16. N. Z. CHO, “Krylov Subspace Wraps Around the Two-Level p-CMFD Acceleration in the
Whole-Core Transport Calculation,” Trans. Am. Nucl. Soc., 123, 1327 (2020).

9炉物理の研究 第74号 (2022年3月)

17. D. WANG and S. XIAO, “A Linear Prolongation Approach to Stabilizing CMFD,” Nucl. Sci.
Eng., 190, 45 (2018).
18. A. ZHU, M. JARRETT, Y. XU, B. KOCHUNAS, E. W. LARSEN, and T. DOWNAR, “AnOptimally Diffusive Coarse Mesh Finite Difference Method to Accelerate Neutron Transport
Calculations,” Ann. Nucl. Energy, 95, 116 (2016).
19. Y. SAAD, “Iterative Methods for Sparse Linear Systems,” SIAM, 2003, pp. 231-234(BiCGSTAB algorithm for Krylov method).
20. H. PARK, D. A. KNOLL, and C. K. NEWMAN, “Nonlinear Acceleration of Transport
Criticality Problems,” Nucl. Sci, Eng., 172, 52 (2012).
21. M. A. SMITH, E. E. LEWIS, and B.-C. NA, “Benchmark on Deterministic Transport
Calculations Without Spatial Homogenization: A 2-D/3-D MOX Fuel Assembly 3-D
Benchmark,” NEA/NSC/DOC(2003)16, Organisation for Economic Co-operation andDevelopment, Nuclear Energy Agency (2003).
22. N. Z. CHO, “An Enlarged 2-D OECD/NEA Benchmark Problem and Two-Level p-CMFD
Solutions,” Korea Advanced Institute of Science and Technology, June 2020;
https://github.com/nzcho/NuraptArchives/blob/master/NZCho_ANS_W2020_Addendum.4.pdf
23. N. Z. CHO, “A Proposed New Framework for Reactor Physics Analysis with Transport
Calculation,” Trans. Kor. Nucl. Soc., Yeosu, Korea, October 24-26, 2018, Korean Nuclear
Society (2018); https://www.kns.org/files/pre_paper/40/18A-132 조남진 I.pdf
24. N. Z. CHO, Y. G. JO, and S. YUK, “Multigroup Transport Equations Derived via
Homogeneity and Isotropy Restoration Theory,” Trans. Am. Nucl. Soc., 115, 592-595 (2016).
25. N. Z. CHO, Y. G. JO, and S. YUK, “A New Derivation of Multigroup Transport Equations
via Homogeneity and Isotropy Restoration Theory,” Ann. Nucl. Energy, 110, 798 (2017).
26. N. Z. CHO, S. Yuk, and Y. G. JO, “Partial Current Discontinuity Factors Determined via
JFNK in HIRE-Theoretic Multigroup Transport Equations,” Fall Meeting, Sapporo, Japan,
September 13-15, 2017, Atomic Energy Society of Japan (2017).

10