NUREG/IA-0416 

International 

Agreement Report

Implementation of Advanced Multigroup Nodal 

and Pin Power Reconstruction Methods 

into PARCS 3.1

 

Prepared by:

Joo-il Yoon, Han Gyu Joo, Seung Hoon Ahn

 

 

 

Anthony Calvo, NRC Project Manager

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

Manuscript Completed: December 2010 

Date Published: May 2012

Prepared as part of 

The Agreement on Research Participation and Technical Exchange

Under the Thermal-Hydraulic Code Applications and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission
AVAILABILITY OF REFERENCE MATERIALS

IN NRC PUBLICATIONS

NRC Reference Material

As of November 1999, you may electronically access 

NUREG-series publications and other NRC records at 

NRC=s Public Electronic Reading Room at 

http://www.nrc.gov/reading-rm.html. Publicly released 

records include, to name a few, NUREG-series 

publications; Federal Register notices; applicant, 

licensee, and vendor documents and correspondence; 

NRC correspondence and internal memoranda; 

bulletins and information notices; inspection and 

investigative reports; licensee event reports; and 

Commission papers and their attachments.

NRC publications in the NUREG series, NRC 

regulations, and Title 10, Energy, in the Code of 

Federal Regulations may also be purchased from one 

of these two sources.

1. The Superintendent of Documents

 U.S. Government Printing Office

 Mail Stop SSOP

 Washington, DC 20402B0001

 Internet: bookstore.gpo.gov

 Telephone: 202-512-1800

 Fax: 202-512-2250

2. The National Technical Information Service

 Springfield, VA 22161B0002

 www.ntis.gov 

 1B800B553B6847 or, locally, 703B605B6000

A single copy of each NRC draft report for comment is 

available free, to the extent of supply, upon written 

request as follows:

Address: U.S. Nuclear Regulatory Commission

 Office of Administration

 Publications Branch 

 Washington, DC 20555-0001 

E-mail: DISTRIBUTION.RESOURCE@NRC.GOV 

Facsimile: 301B415B2289 

Some publications in the NUREG series that are 

posted at NRC=s Web site address

http://www.nrc.gov/reading-rm/doc-collections/nuregs 

are updated periodically and may differ from the last 

printed version. Although references to material found 

on a Web site bear the date the material was 

accessed, the material available on the date cited may 

subsequently be removed from the site.

Non-NRC Reference Material

Documents available from public and special technical 

libraries include all open literature items, such as 

books, journal articles, and transactions, Federal 

Register notices, Federal and State legislation, and 

congressional reports. Such documents as theses, 

dissertations, foreign reports and translations, and 

non-NRC conference proceedings may be purchased 

from their sponsoring organization.

Copies of industry codes and standards used in a 

substantive manner in the NRC regulatory process are 

maintained atC

The NRC Technical Library 

Two White Flint North

11545 Rockville Pike

Rockville, MD 20852B2738

These standards are available in the library for 

reference use by the public. Codes and standards are 

usually copyrighted and may be purchased from the 

originating organization or, if they are American 

National Standards, fromC

American National Standards Institute

11 West 42nd Street 

New York, NY 10036B8002

www.ansi.org 

212B642B4900

Legally binding regulatory requirements are stated only 

in laws; NRC regulations; licenses, including technical 

specifications; or orders, not in 

NUREG-series publications. The views expressed in 

contractor-prepared publications in this series are not 

necessarily those of the NRC.

The NUREG series comprises (1) technical and 

administrative reports and books prepared by the staff 

(NUREGBXXXX) or agency contractors 

(NUREG/CRBXXXX), (2) proceedings of conferences 

(NUREG/CPBXXXX), (3) reports resulting from 

international agreements (NUREG/IABXXXX), (4) 

brochures (NUREG/BRBXXXX), and (5) compilations of 

legal decisions and orders of the Commission and 

Atomic and Safety Licensing Boards and of Directors=

decisions under Section 2.206 of NRC=s regulations 

(NUREGB0750).

DISCLAIMER: This report was prepared under an 

international cooperative agreement for the exchange 

of technical information. Neither the U.S. Government 

nor any agency thereof, nor any employee, makes any 

warranty, expressed or implied, or assumes any legal 

liability or responsibility for any third party=s use, or the 

results of such use, of any information, apparatus, 

product or process disclosed in this publication, or 

represents that its use by such third party would not 

infringe privately owned rights.
NUREG/IA-0416 

International 

Agreement Report

Implementation of Advanced Multigroup Nodal 

and Pin Power Reconstruction Methods 

into PARCS 3.1

 

Prepared by:

Joo-il Yoon, Han Gyu Joo, Seung Hoon Ahn

 

 

 

Anthony Calvo, NRC Project Manager

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

Manuscript Completed: December 2010 

Date Published: May 2012

Prepared as part of 

The Agreement on Research Participation and Technical Exchange

Under the Thermal-Hydraulic Code Applications and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission

iii

ABSTRACT

The original multi-group solver of PARCS uses a nodal expansion method combined with onenode calculation scheme. But the neutronics solver was inefficient for application to the real-size 

reactor cores because of its slow performance. Moreover, the existing pin power reconstruction 

with the corner point balancing method does not support the multi-group structure. In order to 

improve the performance for multi-group problems, the semi-analytic nodal method and pin 

power reconstruction with the source expansion formulation have been incorporated into the 

PARCS code. The purpose of the nodal formulation is to solve the transverse-integrated onedimensional neutron balance equation. Also, a two-level CMFD method in which the two-group 

CMFD is used for accelerating the multi-group CMFD calculation has been constructed to 

improve the calculation efficiency. The newly implemented routines are written in the form of a 

static library of FORTRAN90 that forms link to the PARCS code during compilation. The new 

numerical method has been evaluated over benchmark problems for steady-state and transient 

calculations of real-size reactor core. The calculation results show an excellent agreement with 

the reference results of the benchmark problem. The computing time of the new routines was 

about 4 times shorter than that of the original routine.

v

CONTENTS

ABSTRACT ..................................................................................................................................iii

Page

CONTENTS .................................................................................................................................. v

LIST OF FIGURES ......................................................................................................................vii

LIST OF TABLES ........................................................................................................................ ix

EXECUTIVE SUMMARY ............................................................................................................. xi

FOREWORD...............................................................................................................................xiii

1. INTRODUCTION ...................................................................................................................1

2. NUMERICAL METHODOLOGY............................................................................................3

2.1 Source Expansion Nodal Method....................................................................................3

2.2 Two-Level Coarse Mesh Finite Difference Method .........................................................5

2.3 Pin Power Reconstruction With Source Expansion ........................................................7

2.3.1 Source Expansion .............................................................................................10

2.3.2 Determination of Corner Flux ............................................................................13

2.3.3 Corner Discontinuity Factor...............................................................................14

2.3.4 Reconstruction of Calculation Flow ...................................................................14

3. IMPLEMENTATION.............................................................................................................17

3.1 SENM Library................................................................................................................17

3.1.1 MAIN..................................................................................................................17

3.1.2 BICG..................................................................................................................19

3.1.3 CMFD ................................................................................................................20

3.1.4 SENM2N............................................................................................................20

3.1.5 PPR ...................................................................................................................20

3.1.6 TRAN.................................................................................................................20

3.2 Modification of PARCS code.........................................................................................21

4. VERIFICATION....................................................................................................................23

4.1 Part I – 2D calculation at Fixed T/H Conditions ............................................................23

4.2 Part II – 3D calculation at Hot Full Power Conditions ...................................................25

4.3 Part III – 3D calculation at Hot Zero Power Conditions.................................................27

4.4 Part IV – 3D transient calculation..................................................................................30

4.5 Computing Time Performance ......................................................................................34
vi

5. CONCLUSIONS ..................................................................................................................37

6. REFERENCES ....................................................................................................................39
vii

LIST OF FIGURES

Page

Figure 1 Two-Level CMFD Calculation Flow.............................................................................8

Figure 2 4 sets of net currents at a corner surrounded by four nodes ....................................13

Figure 3 Categories of the SENM library ................................................................................17

Figure 4 Flow chart of runsteady.............................................................................................18

Figure 5 Flow chart of runtrans ...............................................................................................19

Figure 6 Flow chart of drivesenm2n........................................................................................21

Figure 7 Axial power distribution .............................................................................................26

Figure 8 Axial power distribution error ....................................................................................26

Figure 9 Axial power distribution at HZP conditions................................................................29

Figure 10 Axial power distribution error at HZP conditions.......................................................29

Figure 11 Transient core power behavior .................................................................................31

Figure 12 Transient assembly peaking factor...........................................................................32

Figure 13 Transient point pin peaking factor.............................................................................32

Figure 14 Transient core average Doppler temperature...........................................................33

Figure 15 Transient core average moderator density ...............................................................33

Figure 16 Transient core average moderator temperature .......................................................34

ix

LIST OF TABLES

Page

Table 1 Comparison of eigenvalue and assembly power........................................................24

Table 2 Pin power %PWE at ARO...........................................................................................24

Table 3 Pin power %EWE at ARO...........................................................................................25

Table 4 Comparison of CBC, assembly power and core average T/H properties ...................25

Table 5 Pin power %PWE at HFP...........................................................................................27

Table 6 Pin power %EWE at HFP...........................................................................................27

Table 7 Comparison of critical boron concentration and assembly power ..............................28

Table 8 Pin power %PWE at HZP conditions ..........................................................................30

Table 9 Pin power %EWE at HZP conditions ..........................................................................30

Table 10 Transient peak power and power integral ..................................................................31

Table 11 Comparison of computing times for the steady-state calculation (Part II) ..................35

Table 12 Comparison of computing time for the transient calculation (Part IV) ........................35

xi

EXECUTIVE SUMMARY

In order to improve the performance of PARCS code for multi-group problems, a new semianalytic nodal method and a pin power reconstruction method have been developed by 

introducing the source expansion method. These are named as the source expansion nodal and 

pin power reconstruction methods. The neutron balance equations are coupled with each other 

through the fissioning and scattering terms. Therefore, the coupled equations need to be solved 

simultaneously to obtain the group-wise solutions for all the groups. However, this is 

burdensome because the calculating performance deteriorates as the number of energy groups 

increases. The source expansion method enables to avoid simultaneous calculations for solving 

the group-wise neutron balance equations. This is the reason why the source expansion method 

is used in this work. The PARCS code modified with the source expansion method has 

improved its calculating performance without compromising its accuracy. 

The calculation procedure of the source expansion nodal method is as follows.

1. The right hand side of the neutron balance equation which consists of fissioning, scattering 

and transverse leakage terms is approximated by a fourth order polynomial function based 

on the Legendre polynomial. In the initial calculation, the shape of the polynomial function 

can be flattened by using the average flux value which is determined by finite difference 

calculation.

2. Once the right hand side becomes a known value, the neutron balance equation can be 

solved with the proper boundary conditions. The average flux of two neighbored nodes and 

the flux and current continuities at the interface of those nodes are used as the constraints 

in determining the homogeneous solutions. 

3. The steps 1 and 2 are repeated for all energy groups to determine the group-wise general 

solutions.

4. The general solution is expanded to a fourth order polynomial function with the orthogonal

expansion using the least square method.

5. The right hand side is updated using all of the group-wise expanded solutions.

6. The steps from 1 to 5 are repeated until the intra-nodal shape of fission source is converged. 

Then, the group-wise exact solutions of the neutron balance equation are obtained.

The final solution by the source expansion nodal method with sufficient iterations becomes

essentially identical to the solution by simultaneous calculation. The source expansion pin 

power reconstruction method also follows the same procedure of the source expansion nodal 

method except determining the homogeneous solution with the 8 boundary conditions; 4 surface 

currents and 4 corner fluxes of each node.

A PWR MOX/UO2 core benchmark problem has been used as the reference in order to evaluate 

the new source expansion kernel of the PARCS code. In all parts of the problem including 

steady-state and transient cases, the new source expansion kernel provided quite good 

accuracy compared to the original nodal expansion kernel. Especially, the multi-group pin power 

results which are not obtainable with the original kernel have become available with the new 

source expansion kernel and agrees well with the reference. The source expansion kernel 

shows a calculation performance about 4 times faster than the nodal expansion kernel for the 

3D real-size core problems. The newly implemented source expansion kernel has provided a 

great accuracy, yet very fast speed regardless structure of energy groups chosen.

xiii

FOREWORD

This report represents one of the code improvement submitted to fulfill the bilateral agreement 

for cooperation in thermal-hydraulic activities between Korean Institute of Nuclear Safety (KINS) 

and the U.S. Nuclear Regulatory Commission (NRC) in the form of a Korean contribution to the 

NRC’s Code Assessment and Maintenance Program (CAMP), the main purpose of which is to 

validate the TRAC/RELAP Advanced Computational Engine (TRACE) Code.

KINS and other participants in Korea have established a coordinated frame called NuSTEP with 

three objectives: to fulfill the formal CAMP requirements, to keep maintenance and assessment 

activities for safety analysis codes including the MARS code, and to improve the quality of the 

technical support groups that provide services to Korean Utilities, KINS, research institutes, and

engineering companies.

The NuSTEP Review Committee has reviewed this report as the contribution by KEPCO 

Nuclear Fuel and Seoul National University to NuSTEP.

1

1. INTRODUCTION

Although there has been growing interest in the transport based core neutronics analysis 

methods for more accurate calculation with high performance computers, it is yet impractical to 

apply them in the real core design activities because their performance is not so great on 

ordinary desktop or server machines. For this reason, most of the neutronics codes for reactor 

core calculation are still subject to the two-step calculation procedure which consists of 

homogenized group constant generation and core diffusion calculation. In the core calculation 

step which is the main concern of this work, nodal codes based on the diffusion theory have

been used to determine the multiplication factor and flux distribution. Practically, almost all nodal 

codes employ the two-group formulation involving only fast and thermal neutron energy groups 

for the applications to light water reactors. However, numerical calculations with the two-group 

structure is not appropriate in the analysis of the cores loaded with mixed oxide fuels or of the 

fast breed reactors, since the neutron spectrum is influenced more by the core environment 

requiring much more energy groups than two groups.

The PARCS code [1] is one of the very efficient nodal codes in which the analytic nodal method

is implemented within the CMFD formulation. It also employs a unique and accurate twodimensional analytic solution based pin power reconstruction method. The analytic nodal 

method and pin power reconstruction solvers are, however, efficient only for two group solution. 

As an optional solver, PARCS has a nodal expansion method based multi-group nodal solver 

employing the one-node formulation. But this one-node NEM multi-group solver was not 

optimized and is too slow to calculate the real-size reactor cores. Furthermore, the pin power 

reconstruction solver does not support multi-group solutions. It is therefore now planned to add 

an efficient multi-group calculation capability for both the nodal solution and pin power 

reconstruction.

Among the several nodal methods, there have been renewed interests in the semi-analytic 

nodal method [2-6] in the last decade because of the ease of multi-group extension while 

retaining the level of solution accuracy provided by the analytic nodal method. The source 

expansion nodal method [6] developed at Korea Nuclear Fuel and Seoul National University is a

semi-analytic nodal method and optimized for multi-group calculation such that the source 

expansion form provides an advantage of decoupling group dependence and the ease of 

increasing the source expansion order which can reduce calculation time. 

As an acceleration method for the nodal calculation, the coarse mesh finite difference (CMFD) 

formulation [7, 8] is widely used. A two-group CMFD solver was already implemented in PARCS 

to accelerate the two-group analytic nodal kernel and also the multi-group nodal expansion 

kernel. In the multi-group nodal expansion kernel which is based on the one-node scheme, the 

incoming partial currents are specified at the surfaces of a node of interest as the boundary 

condition and the one-dimensional flux shape in each direction is determined so that the node 

average flux and the outgoing partial currents can be newly determined. Therefore, the multigroup spectrum within a node can be updated by the one-node nodal calculation and the CMFD 

calculation is only used as an acceleration scheme. On the contrary, the two-node scheme has 

been chosen for the multi-group source expansion nodal method primarily because the twonode scheme provides better convergence characteristics than the one-node scheme primarily. 

It is due to that better current vs. node average flux relations can be obtainable from the twonode scheme in which the properties of the both nodes of an interface are reflected 

simultaneously in the determination of the interface current. However, in the two-node scheme, 

the node average fluxes of the two nodes are specified as the constraints and the flux shape 

satisfying the current and flux continuity at the interface of the two nodes as well as the node 
2

average flux constraints is determined. Therefore, a multi-group CMFD solution is inevitable 

when using the two-node source expansion nodal method. Since the multi-group CMFD 

calculation converges more slowly than the two-group CMFD calculation because of the source 

iteration, it would be beneficial to establish a two-group CMFD from the multi-group CMFD to 

accelerate the multi-group calculation. Then the efficient two-group solution method such as the 

Wielandt eigenvalue shift method and the simultaneous group solution scheme (as opposed to 

the groupwise solution scheme) can be utilized.

In this document, the multi-group numerical methodologies are presented in Chapter 2 and their 

implementation is introduced in Chapter 3. Then, the computational performance of the new 

implementation is evaluated with a MOX/UO2 fueled benchmark problem [9] for steady-state 

and transient cases in Chapter 4.
3

2. NUMERICAL METHODOLOGY

2.1

The source expansion nodal method (SENM) is based on the transverse-integrated neutron 

diffusion equation which can be written as follows in terms of the normalized coordinate variable 

Source Expansion Nodal Method

2 / u u ξ = u h , where hu is the node width in the u-direction, with the omission of the subscript u

for brevity:

2

2 2

4 () () () () () g g g

rg g g g

eff

D

S L

h k

φ ξ χ φ ξ ψξ ξ ξ ξ∂ − +Σ = + −

∂ , (1)

where the right hand side terms representing the fission source ( ), scattering source (Sg), and 

the transverse leakage(Lg), respectively. The origin of the coordinate is the center of the node 

and thus the range of the normalized coordinate is [-1,1]. This choice of normalization is dictated 

by the choice of Legendre polynomial ( ) Pi ξ for the source expansion. It is supposed that the 

transverse leakage is given by a quadratic polynomial.

In the fourth order source expansion nodal method, the fission and scattering sources are 

approximated by a quartic polynomial so that the entire right hand side (RHS) of Eq. (1) can be 

represented as:

4

0

() () () () () g

g g g gi i

eff i

Q S L qP k

χ ξ ψξ ξ ξ ξ

=

≡ +−=∑ (2)

where Lg(ξ) is given in the form of quadratic polynomial and will be obtained from the leakages 

over the three adjacent nodes. The source expansion coefficient (qi

) is predetermined with an 

initial guess or a previous solution which has been approximated by a fourth order polynomial 

given in Eq. (7). 

Given the RHS as a quartic polynomial, the analytic solution of Eq. (1) is obtained as the 

following:

4

0

( ) ( ) ( ) sinh( ) cosh( ) ( ) H P

g g g g g g g gi i

i

φ ξ φ ξ φ ξ κξ κξ ξ A B c P

=

= += + +∑ , (3)

where the superscripts H and P stand for the homogeneous and particular solutions, 

respectively, and 2

rg

g

g

h

D

κ

Σ = .

The coefficients of the particular solution (cgi) are all determined uniquely by the source term 

coefficients as follow:
4

0 0 24 4 2 4

1 1 32 2 4 2 2

3 4

3 4

11 1 (3 10 ) 105 ,

1 15 1 35 , ,

, and .

g g gg g

rg g g

g g gg g g

rg g g rg

g g

g g

rg rg

c q qq q

c q qc q q

q q

c c

κ κ

κ κ

 

= + ++   Σ  

  

=+ = +    Σ Σ   

= = Σ Σ

(4)

Next, the coefficients of the homogeneous solution should be determined by imposing 

constraints or boundary conditions. In the source expansion nodal method, the even term 

coefficient of the homogeneous solution, Bg, which is obtained from the node average flux

constraint as follows.

( 0 ) sinh( )

g

g g g

g

B c φ

κ

κ = − (5)

Now the only unknown coefficient is Ag and it needs to be determined for each node leading to 

two unknowns per node. These two are determined by solving the following coupled equations 

which are given by the following flux and current continuity at the interface:

sinh( ) sinh( )

2 cosh( ) 2 cosh( )

l r l

g g g

ll l rr r r

gg g gg g g J

A

A

φ κ κ γ

βκ κ βκ κ γ

        =    −       

(6)

where

( )

( )

( )

4

0

123 4

123 4

cosh( ) cosh( ) ( 1)

2 sinh( ) 3 6 10

2 sinh( ) 3 6 10 .

r r l l ir l

g g g g gi gi

i

l ll l l l l l

J g gg g g g g g

r rr r r r r r

g gg g g g g g

B B c c

B aa a a

B aa a a

φ γκ κ

γ βκ κ

βκ κ

=

= − +− −

= − ++ + +

− −+ − +

∑

With the five particular coefficients and the two homogeneous coefficients, the solution given by 

Eq. (3) is uniquely determined and then the flux shape involving the sinh and cosh functions can 

be approximated by a quartic polynomial given by Eq. (7) using the orthogonal property of 

Legendre polynomials as follows:

4

1

( ) ( ) g g gi i

i

φξ φ ξ a P

=

≈ +∑ (7)
5

2 2

2 2

1 1

2 2

3 3

4 4

,

3 cosh( ) ,

,

, and

( ), ( )

( ), ( )

sinh( )

3cosh( ) sinh( ) 3 5 (1 )

7 15 15 sinh( ) (1+ )cosh( ) (6 )

5 9

g

g

g g

g g

g g

g

g g g

g

g g

g g g g

g

i

gi

i i

g g

g g

g g

g g

A

B

A

P

a

P P

a c

a c

a c

a c

κ

κ

κ

κ

κ

φξ ξ

ξ ξ

κ

κ

κ κ

κ κ

κ

κ κ κ

κ

 

+ −  

 

    −

 

 

+  

 

< > = < >

=

= − +

= − +

= − 2 2 2 4 . 21 45 105 sinh( ) (2 ) cosh( ) (1 ) g

g g

g g g g

κ B

κ

κ

κ κ κ

   

 

+ −+ +

(8)

Here the bracket stands for inner product of the functions involving an integral over [-1,1]. Note 

that the orthogonal expansion corresponds to the least square fitting of a function based on the 

weighted residual procedure.

The source expansion nodal method involves a scheme of source iteration. It is thus necessary 

to establish a proper condition to check the convergence of the multi-group flux solution. Since 

the change in the flux shape is directly reflected into the fission source shape, the fission source 

shape within the local nodes can be monitored to check the convergence. Instead of checking 

fission source convergence, it is also possible to check convergence of individual group fluxes. 

But that kind of strict convergence check is not necessary because the source level and shape 

keep changing during the outer iteration. Specifically, the following convergence parameter is 

used in the source expansion nodal method:

( )

1/2 2 1 ( 1) ( 1) ( )

( 1) 1

1 1 () () 2

n n n

n δ ψ ξψ ξ ξd

ψ

+ +

+ −

  =   −

  ∫ (9)

with n being the source iteration index.

By expanding the fission source into a quartic polynomial using Legendre polynomials and using 

the orthogonal property of the Legendre polynomials, the above integral can be reduced as 

follows:

1/2 4

( 1) ( 1) ( ) 2

( 1)

0

1 1 ( ) 2 1

n n n

n i i

i

q q i ψ ψ δ

ψ

+ +

+

=

  =   − +   ∑ , (10)

where ( ) n

i qψ is the i-th order fission source expansion coefficient at iteration step n. In case of 

non-fuel nodes, the fission source is replaced with the lowest energy group flux.

2.2

There can be two views of the coarse mesh finite difference (CMFD) formulation which employs 

a correction term in the current vs. node average flux relation. The first view is to regard the 

Two-Level Coarse Mesh Finite Difference Method
6

CMFD as an acceleration scheme for the higher order methods whereas the second view is to 

regard it as a complete nodal balance formulation employing higher order current relations. The 

latter view enables more extensive and efficient application of the CMFD formulation to various 

problems. For instance, a CMFD standalone spatial kinetics calculation is possible by applying a 

suitable temporal differencing to the CMFD balance equations. In this case, only occasional 

updates of the higher order correction factors are required and the correction factor appearing in 

the current vs. node average flux relation can be obtained by solving the local problem with the 

source expansion nodal method.

Suppose that multi-group surface currents are specified for the local problem like the source 

expansion nodal calculation. It is then possible to define the current correction factor such as:

( ) ˆ r l

gg g g

g r l

g g

D J

D φ φ

φ φ

− + = − +



(11)

where the superscripts, l and r, denote the left and right nodes, respectively, g J is the surface 

current, and D

gˆ  is the finite difference method (FDM) based proportional constant relating the 

current to the flux difference defined as:

2

( )

r l

g g

g r l

g g

D D

D

hD D = +

 .

The multi-group CMFD calculation can be performed to update globally the node average fluxes 

with the correction factor. However, the convergence of the multi-group calculation is generally 

slower than the corresponding two-group problem because of the ineffectiveness of the source 

iteration scheme. For instance, the Chebyshev acceleration scheme employed in the multigroup eigenvalue problem shows an inferior performance to the Wielandt eigenvalue shift 

method which is readily employed in a two-group problem. This motivates a dynamic two-group 

formulation during the multi-group CMFD calculation.

Once the multi-group node average fluxes and surface currents are available from the multigroup CMFD calculation, it is then possible to define the corresponding two-group cross section 

and the current correction factor as follows:

ˆ

ˆ ,

g

g g g

g S

α α φ ∈

Σ= Σ ∑ (12)

( ) ˆ r l

gg g g

g r l

g g

D J

D φ φ

φ φ

− + = − +



(13)

where gˆ is the coarse group index, α is the reaction type, and

ˆˆ

ˆˆ

ˆ

, , ( 1,2). ˆ

g g

g

g g gg g

g g S g S

J Jg φ φ φφ φ ∈ ∈= = = = ∑ ∑

Here Sgˆ designates the set of the fine group number belonging to coarse group gˆ . The coarse 
7

group current correction is used in the two-group CMFD formulation to relate the interface 

current with the two node average fluxes of the interface as follows:

ˆˆˆˆˆˆˆ ˆ ( )( ) rl rl

g gg g gg g JD D =− − − +  φφ φφ (14)

Once the two-group CMFD calculation is performed, the two-group node average flux 

distribution becomes available and is used to adjust the multi-group fluxes with the multi-group 

intra-node spectrum stored in the previous multi-group calculation. The subsequent multi-group 

CMFD calculation is done once again to update the solution satisfying the nodal balance 

equation with the expanded multi-group fluxes. Then, the source expansion nodal calculation is 

repeated with the updated average fluxes.

This overall flow of two level CMFD calculation involving two-group and multi-group CMFD 

solutions is depicted in Fig. 1. As seen in the flow chart, there are several decision points in the 

two-level CMFD calculation. First, the overall convergence check is made only once per nodal 

update cycle. It is done at the first MG CMFD group sweep performed right after each nodal 

update. This is to make sure that the overall convergence is achieved even after the nodal 

update. The exit from the MG CMFD to 2G CMFD is made when sufficient residual error 

reduction is achieved. The exit from the 2G CMFD to MG CMFD is also made after certain 

residual reduction during the 2G calculation. The exit from the MG to the nodal update is made 

at the first MG CMFD after a 2G CMFD calculation.

2.3

After integrating the three-dimensional neutron diffusion equation axially over a plane of 

thickness,

Pin Power Reconstruction With Source Expansion

z h , the following two-dimensional balance equation is obtained for each group with all 

the source terms placed on the right hand side (RHS).

2 1 (, ) (, ) (, ) (, ) (, ) g g rg g g fg g gg g gz

eff g g g

D xy xy x y xy L xy k φ φ χ νφ′ ′ ′ ′ φ ′ ′≠

− ∇ +Σ = Σ + Σ − ∑ ∑ (15)

where the axial transverse leakage is defined as ( , ) ( ( , ) ( , )) / T B L xy J xy J xy h gz = − gz gz z in term of 

the currents at the top and bottom surfaces of the plane and keff is the multiplication factor. The 

entire RHS terms can be approximated by a polynomial of two spatial coordinate variables. The 

polynomial for a two-dimensional shape is iteratively updated. With such a polynomial

approximation, the solution of Eq. (15) would be obtained straight forwardly for a square node 

whose width is h. 
8

Figure 1 Two-Level CMFD Calculation Flow

Since it is advantageous to use the Legendre polynomial for the polynomial approximation 

owing to its orthogonality property, we first normalize the independent variable such that it varies 

from -1.0 to 1.0 in the node. This leads to the following equation with the group index g omitted:

(16)

where ξ η = = 2/, 2/ xh yh , and Q(,) ξ η represents the distribution of the entire sources. As 

an approximation for the source distribution, we use a quartic polynomial given in terms of 

Legendre polynomials:

4 4

,

0 0

4

(,) () () ij i j

i j

i j

Q qP P ξ η ξ η = =

+ ≤

= ∑∑ (17)

where ( ) Pi ξ is the i-th order Legendre polynomial. Note that this 15-term polynomial contains

2 2

22 2

4 (,) (,) (,) rD Q h φξη φξη ξη ξ η

  ∂ ∂

− +   +Σ =

  ∂ ∂
9

fourth order cross terms such as 1 3 P P () () ξ η and 2 2 P P () () ξ η representing a quartic variation 

in each direction. In the following, it is assumed that the source term coefficients are known and 

so are the 8 boundary conditions for each group.

The solution of the non-homogeneous second order differential equation, Eq. (16), consists of 

the homogeneous and particular solutions. The particular solution has the same functional form 

as the source polynomial, Eq. (17). The coefficients, pi,j, of the polynomial particular solution, 

( , ), φ ξη p can be determined by the method of undetermined coefficient as follows.

2

0,4 4,0 2,2 0,2 2,0 0,4 4,0 0,0

0,0 4 3 2 2

48 (35( ) 6 ) 4 (3( ) 10( ))

r r r

D q q q Dq q q q q

p h h

+ + ++ +

=+ +

Σ Σ Σ ,

0,3 2,1 0,1 3,0 1,2 1,0

0,1 2 2 1,0 2 2

12 (5 ) 12 (5 ) ,

r r r r

Dq q q Dq q q

p p h h

+ +

= += +

ΣΣ ΣΣ , (18)

2,2 0,4 0,2 2,2 4,0 2,0

0,2 2 2 2,0 2 2

(12 140 ) (12 140 ) ,

r r r r

Dq q q Dq q q

p p h h

+ +

= += +

ΣΣ ΣΣ ,

,

, for other (i,j) combinations. i j

i j

r

q

p = Σ

The homogeneous equation of Eq. (16) can be rewritten by dividing by 2 4 / D h as follows:

2 2

2

2 2 (,) (,) φξη φξη Bt ξ η

  ∂ ∂   + =

  ∂ ∂ (19)

where the dimensionless buckling is defined as:

2 2

2

2 or 

44 2

r

t t

hh h B B

DL L

Σ = = = (20)

with L being the diffusion length. Eq. (19) can be solved by separation of variables after splitting 

the total buckling as follows:

22 2 BB B tt t = + , , ξ η (21)

where = cos and = sin . BB BB t t , , ξ η θ θ t t There are infinite numbers of splittings by Eq. (21) and 

the solution satisfying the boundary condition at all the boundary points can be obtained in 

terms of an infinite series. However, we choose only eight pairs with the following angles which 

would require only 8 boundary conditions: 

, 0..7 4 k k k π θ = = (22)
10

This leads to the following form of the homogeneous equation:

7

( cos sin )

0

(,) ˆ Btkk

H k

k

a e ξ θη θ φ ξη +

=

= ∑ (23)

which can be rewritten in terms of sinh and cosh as:

12 34

5 6

7 8

( , ) sinh( ) cosh( ) sinh( ) cosh( )

 sinh cosh( ) sinh( )sinh 22 22

 cosh sinh( ) cosh( )cosh 22 2 2

H t tt t

tt tt

tt t t

a Ba Ba Ba B

BB BB

a a

BB B B

a a

φ ξη ξ ξ η η

ηη ξη

ηη ξ η

=+ ++

    + +        

    + +        

(24)

The four surface average currents and four corner fluxes are chosen as the boundary conditions 

to determine the eight coefficients of Eq. (24). The reason for choosing the four surface average 

current as the boundary condition is to force that the average of the 2D flux solution over the 

node be exactly the same as the one determined in the 3D nodal calculation. This is possible 

because the nodal solution and the two-dimensional solution satisfy the nodal balance equation 

involving the same surface average currents.

2.3.1 Source Expansion

In deriving the eight homogeneous solution coefficients, it is convenient to make the boundary 

condition apply only to the homogeneous solution by subtracting the particular solution part from 

the whole boundary condition. This is because the particular solution is always determined 

uniquely once the source polynomial is specified. For example, the homogenous part of the flux 

at the northwest corner ( 1, 1) ξ η =− = is obtained as follows:

0,0 0,1 0,2 0,3 0,4 1,0

2,0 3,0 4,0 1,1 1,2 2,1 2,2

 (

)

H P

NW NW NW

NW ppp ppp

p p p pp p p

φ φφ

φ

= −

= − ++++−

+ − + −−+ +

. (25)

And the surface average current at the west surface is given by:

1

, 1

1,0 2,0 3,0 4,0

1 ( 1, ) 2

2 = ( 3 6 10 )

Hl P

xl x x

l

x

JJ J d

D J pp p p h

η η − =− −

+ −+−

∫

(26)

where the superscript l designates the left surface and the current for the particular solution part 

is given by:

2 (,) (,) P P

x

D J

h

φ ξη ξ η ξ∂ = − ∂ . (27)
11

The corner flux at other corner points and the surface average current at other surfaces for the 

homogeneous solution component can be obtained similarly.

In terms of these boundary conditions for the homogeneous solution, the coefficients are 

determined as follows:

2

, 1 , 8

1 2

,

2

, 3 , 8

3 4

, 1 2

5 6 2

, 3

7

,

8

2 2 , , 2( ) 2

2 2 , , 2( ) 2

( 2) , , 2( )

( 2), 2( )

( (

HH H

t xs t xd

t g g tg

HH H

t ys t yd

t t

H H H

t xs

t

H H

t ys

t

H

t xd y

B hJ D B hJ DS a

a a

D S BC D S B

B hJ D B hJ DS a

a a

D S B C DSB

B hSJ CD

a a

D BC S X S

B hSJ CD

a

D BC S X

B Ch J J

a

φ

φ

φ φ

φ

+ + = = − −

+ + = = − −

+ = = −

+ = −

+ =





, 4

2 2

)2 ) , 2( 2 )

H H

d

t

DS

D B SC CS

+ φ

−

(28)

where

2 2

sinh( ), sinh , 2

sinh( ), cosh , 2

,

t

t

t

t

B S BS

B C BC

XSC

  = =    

  = =    

= ×







, ,,

, ,,

, ,,

, ,,

1

2

H HH

x s xr xl

H H H

x d xr xl

H HH

y s t yr yl

H HH

y d yr yl

J JJ

J JJ

J JJ B

J JJ

   +    − = +

−   

and

,

1

2

3

4

1

4

H H HH H

NE NW SE SW

H H H HH

NE NW SE SW

H H H HH

NE NW SE SW

H H H HH

NE NW SE SW

φ φ φ φφ

φ φφ φφ

φ φ φ φφ

φ φ φ φφ

    − +−     − −+   =   + −−         + ++

.
12

The coefficients given by Eqs. (18) and (28) for the particular and homogeneous solutions, 

respectively, complete the final solution of:

(,) (,) (,) φξη φ ξη φ ξη = + H P . (29)

In the above derivation, it was assumed that the source distribution is known and describable by 

a quartic polynomial. With the new 2D solution flux solution available, the primary source 

distribution which consists of the fission and scattering sources should be updated for use in the 

next step of the iterative solution sequence. Since the 2D flux distribution for each group 

contains the exponential function components as well as the polynomial function components 

whereas the source is to be approximated by only a polynomial function, it is necessary to 

obtain first the 2D quartic polynomial approximation of the flux distribution for each group. Then 

the source term coefficients can be obtained by multiplying the flux coefficients by the proper 

cross section and then by summing over the groups.

As a method for approximating the 2D flux distribution, the least square fitting which makes use 

of the orthogonal expansion of an analytic function in terms of the Legendre function is 

employed. The expansion coefficients can be given by: 

1 1

1 1

, 1 1 2 2

1 1

1 1

1 1

(,) () ()

() ()

1 (2 1)(2 1) ( , ) ( ) ( )

4

i j

i j

i j

i j

P P dd

c

P P dd

i j P P dd

φξη ξ η ξ η

ξ η ξη

φξη ξ η ξ η

− −

− −

− −

=

=++

∫ ∫

∫ ∫

∫ ∫

(30)

In Eq. (16), the semi-analytic form of the flux solution, Eq. (29) is to be inserted in the integrand 

and the integral is trivial. Since it is cumbersome to present all the 15 coefficients here, only the

fourth order term is given below:

4,0 4 32 2

24 3 2 8 4,0

105 105 45 10 9

18 420 210 2 90 2 10 2 .

t t tt

tt t t t

S C SC

c S a

B B BB

SS C S C Sa p BB B B B

 

= − +−+    

 

+ − + − ++  

 

    

(31)

Once the coefficients of the polynomial flux are updated, the source coefficients can be 

determined by summing the groupwise source terms:

, , , ,

g g g g g

i j fg i j g g i j i j

eff g g g

q c cl k

χ

ν ′ ′

′ ′

′ ′≠

= Σ +Σ − ∑ ∑ . (32)

The updated source coefficients will be used to determine the new particular solution at the next 

iteration step.
13

2.3.2 Determination of Corner Flux

Suppose that the two outgoing currents are known at the interface of two nodes in a onedimensional problem as follows:

1 1

1

4 2

net

s J J φ + 1

= + and 2 2

1

4 2

net

s J J φ + 1 = − . (33)

Here the numeric subscripts indicate the node number. Node 1 is for left node and 2 is for the 

right node. It is assumed that the net currents determined at the two nodes are not the same 

because the current continuity condition is not imposed. With the two outgoing currents, we can 

determine the new surface flux as follows:

12 1 2 2( ) ( ) new net net

s s φ φ JJ J J + + = + =+ − . (34)

The equation above states that the new surface flux increases by the difference in the net 

current. 

Now consider a corner surrounded by four nodes as shown in Fig. 2. 

Figure 2 4 sets of net currents at a corner surrounded by four nodes

From the 2D flux solution of each node, x- and y-directional net currents at the corner become 

available. For the four nodes, there are 4 sets of x- and y-directional currents. In principle, the 

four x-directional currents must be the same. But since no corner current continuity condition is 

imposed, there can be difference in the net currents. Using the four sets of net current difference, 

four estimates of the new corner flux can be obtained by using Eq. (34). By taking the average 

of the four, the new corner flux is obtained as:

( 1, 2, 3, 4, 1, 3, 2, 4, ) 1

4

new net net net net net net net net

c c x xx xy y yy φ φ =+ − + − + − + − JJJJJJJJ . (35)

This iterative update scheme will converge if the summation of the net current differences 

vanishes, which is the corner point balance condition. 

1 2

3 4
14

2.3.3 Corner Discontinuity Factor

In the nodal calculation, the assembly discontinuity factor is used to allow discontinuity of 

homogeneous flux shapes in order to preserve the current of the heterogeneous configuration. 

This discontinuity must be incorporated when determining the corner fluxes primarily because 

the nodal flux solution is discontinuous at the interfaces. However, there should be given

another consideration for the discontinuity. That is related with the subsequent multiplication of 

heterogeneous form factors to the homogeneous flux shape to reflect intra-assembly 

heterogeneity. Let the discontinuity factor which would correspond to the form factor at the 

corner be i ζ for assembly i. Then the heterogeneous flux at the corner is obtained by 

multiplying the corner discontinuity factor (CDF) to the homogeneous flux at the corner. Namely,

,

het hom φ ζφ c i ci = . (36)

Since the heterogeneous flux is supposed to be continuous, the homogeneous flux should be 

discontinuous. 

Once the common corner flux is determined by Eq. (35), we can take that as the heterogeneous 

corner flux and the homogeneous corner flux for each assembly is determined by diving it by the 

CDF. The homogeneous corner flux is then used in Eq. (28).

2.3.4 Reconstruction of Calculation Flow

The pin power reconstruction is to be performed after the multi-group transverse-integrated 

nodal solution is obtained. The basic information available for each node from the nodal 

calculation consists of one node average flux, four surface average fluxes, and four surface 

average currents. The reconstruction starts with this information to determine the homogeneous 

flux shape with 2D source expansion nodal method (SENM) presented in the previous section. 

Repeated 2D SENM calculations are needed to update the corner fluxes and the source 

expansion iteratively. The heterogeneous pin power distribution is finally constructed by 

multiplying the heterogeneous form factor to the homogeneous flux distribution. This 

reconstruction procedure is explained with a little more details below.

The initial corner flux can be derived by employing the method of successive smoothing

(MSS)[10] which yields the following expression for the corner flux:

x y φφ φ φ cs s =+− (37)

where u φs is the surface average flux for the u-directional edge from the corner and φ is the 

node average flux. Another approximation of the corner flux is possible by simply expressing the 

2-D flux distribution with a product of the two one-dimensional fluxes (x and y). This yields with 

the proper normalization:

x y

s s

c

φ φ φ φ = . (38)

Note that Eq. (37) can potentially yield a negative flux while Eq. (38) cannot. Since the MSS 

estimate would give a better estimate than the latter one in which the cross term effect is totally 
15

neglected, we use Eq. (37) as the default formula for the initial corner flux, but switch to Eq. (38)

if the MSS estimate becomes negative. For each corner, normally four estimates of the corner 

flux value are available, one from each node surrounding the corner. The final corner flux is 

determined by merely averaging the four available estimates.

With the four corner flux information, there are now 13 pieces of information available for each 

node. The initial source shape can be determined by fitting these values with a 13 term 

polynomial. This 13 term polynomial is constructed by omitting the two fourth order cross terms

of 3 1 P P () () ξ η and 1 3 P P () () ξ η .

Once the source expansion and corner fluxes are determined for each group, it is now possible 

to determine newly the 15 particular solution coefficients and then the 8 homogeneous solution 

coefficients for each node and each group. The iterative calculation to update the corner flux 

and the source expansion coefficients proceed with the iteration scheme shown in Fig. 2. In the 

practical applications, it turned out that the initial guess of the corner fluxes given by Eq. (37) or 

(38) are sufficiently accurate so that only about 5 iterations are enough to obtain converged 

corner fluxes as well as the converged source expansion coefficients.

Once the converged homogeneous flux distribution is obtained, the pinwise flux distribution can 

be obtained for the homogenized node. The heterogeneity within the node is then reflected by 

multiplying the groupwise power form factor to the groupwise flux.

17

3. IMPLEMENTATION

The multi-group solution is implemented by modification of the PARCS code and new addition of 

routines regarding the multi-group solutions in the form of a static library, namely source 

expansion nodal method (SENM) library. The SENM library contains various Fortran90 modules 

for determining node-wise flux distribution during two-group and multi-group steady-state and 

transient calculations with fixed group constants by use of the source expansion nodal method 

based on coarse method finite difference (CMFD) formulation. The Wieland shift method is

employed to achieve fast fission source convergence by accelerating two-group CMFD 

calculation. And the multi-group CMFD calculation is accelerated by two-group CMFD 

calculation with condensed two-group constants. In addition, the bi-conjugate gradient stabilized 

(BICGSTAB) method with blockwise and pointwise incomplete LU precondition technologies is

implemented to solve fixed source problems. Modification of PARCS to integrate the library was 

confined to input handling and calling driver routines of the library.

3.1

Variables and subroutines of the SENM library are categorized into six sets shown in the 

following figure in terms of functionalities and each set is composed of several modules in 

FORTRAN90. 

SENM Library

Figure 3 Categories of the SENM library

3.1.1 MAIN

This first category contains eight Fortran90 modules in which numerous subroutines and 

variables corresponding to the geometry, cross-sections and solutions exist. Several 

subroutines in this category are opened by calls from PARCS while the other categories are 

encapsulated. For managing the memory of the computer and avoiding inadvertent memory 

copies, most of the variables are declared as POINTER that can be point variables belonging to 

PARCS without duplicating memory allocations and initialized when the subroutine, initvar, is 

called.

Among all of the modules, the main module senm serves as the main driver routine for steady 

and transient calculations and calls runsteady and runtrans. Once the cross-sections are 

SENM Library

BICG

PPR

MAIN

TRAN SENM2N

CMFD
18

calculated incorporating the changes in the boron and thermal-hydraulic conditions in PARCS, 

the driving subroutines, runsteady and runtrans, are called for one-step of steady and transient 

calculations, respectively. These determine the average flux distribution, surface current and flux 

as the results of the CMFD calculations. The following flow charts show the calling sequence of 

the subroutines in the drivers. In case of two-group problems, only the two-group CMFD 

calculation is executed for the efficient calculation and an argument to determine to drive the 

nodal scheme is provided. Convergence check is made through the relative l2 norm of fission 

source.

FLAG2G SETLSMG

UPDLS

DRIVECMFD2G

FLAG2G

IFNODAL

NO UPDLS YES FACILUMG DRIVECMFDMG

NO

YES

NO

NO EXPPHI DRIVECMFDMG

YES

YES DRIVENODAL

START

END

NO

YES SETLS2G FACILU2G

UPDDHAT2G UPDDTIL2G COLXS COLPHI

Figure 4 Flow chart of runsteady
19

FLAG2G SETLSMG

DRIVECMFD2G

FLAG2G

IFNODAL

YES FACILUMG DRIVECMFDMG

NO

YES EXPPHI DRIVECMFDMG

NO

YES DRIVENODAL

START

END

NO

UPDDHAT2G UPDDTIL2G COLXS COLPHI

SETLS2G

Figure 5 Flow chart of runtrans

3.1.2 BICG

The BICG method stands for Bi-Conjugate Gradient Stabilized method and this method was 

implemented for solving two-group and multi-group fixed source problems with fully optimized 

routines in this category. There are two modules; bicg2g and bicgmg. In case of two-group 

problems, the fixed source problem is represented by the blockwise matrix-vector form in which 

the matrix has N by N block elements and each block element is a 2 by 2 matrix representing 

group coupling. The module, bicg2g, provides several subroutines to solve the fixed source 

problem with BICGSTAB preconditioned by blockwise incomplete LU technique. The other 

module, bicgmg, is also a linear system solver but has point elements, not block elements, since 

the linear system is constructed for only one energy group and a scheme of source iterations 

are employed to consider group-coupling by fission and up-scattering sources. Moreover, the 

point-wise incomplete LU factorization is used for accelerating BICGSTAB in case of the multigroup solver. 
20

Among the subroutines of bicg2g, facilu2g and solbicg2g called from the module cmfd2g, the 

subroutine facilu2g factorize the matrix into incomplete LU with block elements. By looping over 

rows, the incomplete LU factorization is performed for matrices according to each row and 

inverse matrices for which the variables del, al and au were used are obtained.

3.1.3 CMFD

This category provides all functionalities solving steady-state and transient problems with the 

CMFD formulation, in which there are two modules, bicg2g and bicgmg. These modules consist 

of four main subroutines: upddtil2g(mg), upddhat2g(mg), setls2g(mg) and drivecmf2g(mg). 

Additionally, the subroutine wiel is provided for two-group CMFD calculation for accelerating the 

calculation with the Wielandt Shift method.

3.1.4 SENM2N

The only module senm2n in this category is a Fortran90 module providing subroutines and local 

variables for calculating the transverse integrated one-dimensional neutron balance equation 

with the source expansion nodal method. The driving subroutine drivesenm2ncalls a set of local 

subroutines in the module sequentially as shown in the following flow chart and determines net 

current and flux at each surface of all nodes.

3.1.5 PPR

The source expansion pin power reconstruction with multi-group structure was implemented in 

this category and there is a main module ppr that contains a set of subroutines. The module 

does determine the homogeneous pin power distribution only and not heterogeneous pin power 

distribution which is calculated by multiplying the heterogeneous form function to the 

homogeneous pin power distribution. Instead, PARCS determines the heterogeneous one using 

the homogeneous one provided by the SENM library. It is a proper programming approach to do 

so because the SENM library does not need to know the heterogeneous form function. The 

subroutine ppr is composed of initppr, calhomo and calpinpow.

3.1.6 TRAN

All modules of the category, TRAN, are materials for the multi-group transient calculation and 

used in the subroutine runtrans. Several subroutines for initializing, subsets of transient 

calculations and numerous variables represent kinetic parameters and pre-calculated values in 

the modules; tranxsec and tran. The first module is tranxsec which provides variables

corresponding to additional group constants needed for the transient calculation such as

delayed neutron parameters. And the subroutine, malloctranxsec, has a role of allocating 

memories for the variables. This subroutine uses the memory sharing technique to avoid 

duplicate memory allocation by PARCS. The other module tran is a main module to provide 

subroutines to calculate transient parameters.
21

Figure 6 Flow chart of drivesenm2n

3.2

In PARCS, several subroutines and variables were added or modified to integrate the SENM 

library. 

Modification of PARCS Code

 SenmIntgM (added)

The module senmintgm was added to PARCS to integrate the SENM library with PARCS. 

The module provides several subroutines to initialize and calculate steady-state, transient 

and pin power reconstruction calculations as an agent for the SENM library which

encapsulates the SENM library from PARCS and allows modification of the PARCS code as 

little as possible. 

 MemAllocM (modified)

Variables limited to the two-group structure were modified for multi-group calculation and 

the module moxtrm was allocated. 

 InitializedM (modified) 

This was modified to handle the benchmark table-sets and to initialize multi-group fission 

source terms. 

START

RESETSENM2N RESETJNET CALTRL

CALSENM2N

BOUNDARY

CALSENM2NBND

YES

NO

END

UPDSRCCFF CALBY1N

CALBY2N

CONVERGE UPDPHICFF

CALJNET

NO

OVER ALL SURFACES
22

 InpProcM (modified)

In this module, some input variables were added to choose SENM as a solver and solve 

OECD/NEA MOX Transient benchmark problem with the module moxtrm.

 NormalizedM (modified) 

The subroutine norm of the module senm is called when the main solver is configured as 

the source expansion nodal method.

 XsecFdbkM (modified)

If the input 'FUNC_TYPE' = 20, the option ioptxsf is 20 and the xsecmox which calculates 

cross-sections for the MOX transient benchmark problem is called.

 SSeigM (modified)

This was modified to call outerss of senmintgm when SENM is chosen as the steady-state 

solver.

 TransDriverM (modified)

The subroutine outertr of senmintgm was added in the subroutine transient of this module.

 DepMainM (modified)

The subroutine driveppr of senmintgm is called when the solver is set as SENM.

 Parcsexe (modified)

The subroutine senmintg of senmintgm is called after finishing initialization which is 

executed with p_init in this subroutine when SENM is assigned for 'NODAL_KERN' in the 

input file.

 ItrcntlM (modified)

The logical variable lsenm was added to check whether the solver is SENM or not.
23

4. VERIFICATION

The OECD/NEA AND U.S. NRC PWR MOX/UO2 CORE TRANSIENT BENCHMARK Problem [9] published 

by Purdue University in 2003 was used for verification of the newly implemented multi-group 

solver. The benchmark problem report contains four parts representing calculation results for 

steady-states in various conditions and a fast transient as follows:

 Part I: Multiplication factor, rod worth, assembly and pin power at fixed T/H conditions

 Part II: Critical boron concentration, assembly and pin power at HFP 

 Part III: Critical boron concentration, assembly and pin power at HZP

 Part IV: Simulation of control rod ejection accident from HZP conditions

Fundamental functionalities of the multi-group routines, such as implementation of numerical 

formulations and basic integration with PARCS, are evaluated in Part I. And then, the steadystate routines including integration with T/H routines of PARCS are tested and the performance 

for solving the steady-state problem in various conditions is assessed in Parts II and III. Finally, 

the transient routines are evaluated in Part IV and computational performance of the new multigroup routines is compared to the existing routines, additionally. Since the heterogeneous

solutions calculated by DeCART in the final benchmark report are used as the references in 

Parts I and III, accuracy of the new multi-group kernel is evaluated by comparing the calculation

results quantitatively to those of the existing nodal and pin power routines. On the other hand, it 

is tried to verify implementation of the new routines by comparing results qualitatively in Parts II 

and IV since any heterogeneous solution by DeCART was not provided. 

4.1

This part was to solve steady-state problem with two-dimensional core configuration with fixed 

T/H conditions of hot zero power. The basic capabilities of the multi-group nodal and pin power 

routines are assessed in this part.

Part I – 2D Calculation at Fixed T/H Conditions

Table 1 shows a comparison of multiplication factors and assembly power errors against Part I

of the benchmark report. All solutions determined by the SENM kernel match well to the results 

determined by the ANM and NEM kernels at ARO and ARI conditions. Comparing two-group 

results indicated by ANM 2G and SENM 2G, it is verified that the source expansion nodal 

method has an accuracy equivalent to the ANM solution. The multi-group results of the SENM 

kernel are also comparable to those of the existing NEM kernel. 

Like the existing kernels, assembly power distributions of the SENM kernel were in a very good 

agreement to the heterogeneous solution by the DeCART code at ARI as well as ARO 

conditions in all group structures within 1 ~ 2% of PWE and EWE. It was noticeable that the twogroup results agree well with the reference contrary to the expectation that the assembly power 

distribution makes a larger distortion in the MOX loaded core. The reason could be found in the 

spatial discretization. Since the cross-section is not interpolated for the fixed T/H condition, the 

difference in reaction rates between the homogeneous and heterogeneous solutions is only 

ascribable to the spatial discretization which is unrelated to the group configuration.
24

Table 1 Comparison of eigenvalue and assembly power

Nodal

Kernel

Eigenvalue TOTAL ROD

WORTH

[dk/k, pcm]

Assembly Power Error

ARO ARI

ARO ARI %PWE %EWE %PWE %EWE

ANM 2G 1.06379 0.99154 6850 0.96 1.63 1.67 2.18

SENM 2G 1.06379 0.99154 6850 0.83 1.63 1.57 2.18

NEM 4G 1.06376 0.99136 6865 0.90 1.42 1.61 2.26

SENM 4G 1.06378 0.99136 6868 0.76 1.47 1.59 2.20

NEM 8G 1.06354 0.99114 6868 0.86 1.25 1.65 2.49

SENM 8G 1.06358 0.99114 6871 0.84 1.29 1.49 2.37

*DeCART 1.05852 0.98743 6801 ref ref ref ref

*DeCART: Reference

Table 2 and Table 3 show the pin power results at the ARO condition. Like the eigenvalue and 

the assembly power, the pin power results also agree very well with the reference results. The 

existing pin power kernel of PARCS was based on the corner point balancing method with 

corner discontinuity and provided accurate results because it solved the two-dimensional 

neutron diffusion equation analytically, but it was only applicable to the two-group problems. On 

the other hand, the new source expansion pin power reconstruction kernel provides accurate 

results with any group structures as shown in the tables. There is slight difference between the 

two-group and multi-group results but they are consistent with respect to the assembly power 

distributions.

Table 2 Pin power %PWE at ARO

Nodal

Kernel

Assembly Position

(A, 1) (B,2) (C,3) (D,4) (E,5) (F,6)

ANM 2G 0.28 0.21 0.29 0.54 0.32 0.51

SENM 2G 0.23 0.25 0.33 0.41 0.28 0.77

NEM 4G - - - - - -

SENM 4G 0.23 0.21 0.35 0.36 0.30 0.77

NEM 8G - - - - - -

SENM 8G 0.28 0.18 0.38 0.41 0.32 0.57
25

Table 3 Pin power %EWE at ARO

Nodal

Kernel

Assembly Position

(A,1) (B,2) (C,3) (D,4) (E,5) (F,6)

ANM 2G 0.54 0.45 0.41 0.79 0.60 1.05

SENM 2G 0.73 0.52 0.60 1.06 0.78 6.63

NEM 4G - - - - - -

SENM 4G 0.75 0.48 0.60 1.03 0.82 6.10

NEM 8G - - - - - -

SENM 8G 0.94 0.50 0.67 1.07 0.84 5.68

4.2

Part II describes three-dimensional steady-state calculations with thermal-hydraulic feedback at 

hot full power condition. Unlike Part I, the PARCS 2G calculation results are used as the basis 

for comparison of assembly and pin-power distribution, since there exist no reference results

calculated by DeCART in the final report of the benchmark problem. Critical boron 

concentrations, assembly power and core-averaged T/H properties are compared for different

nodal kernels and group structures in Table 4.

Part II – 3D calculation at Hot Full Power Conditions

The critical boron concentrations calculated by the SENM 2G kernel are very close to the ANM 

2G results within one ppm and the assembly power distributions and core-averaged T/H 

properties are essentially identical. The SENM kernel shows good agreement with the existing 

NEM kernel in four- and eight-group calculation results.

Table 4 Comparison of CBC, assembly power and core average T/H properties

Nodal

Kernel

Critical 

Boron 

Concent.

[ppm]

Ass. Power 

Error Core Average T/H Properties

%PWE %EWE

Doppler

Temp.

[K]

Mod.

Density

[kg/m3]

Mod.

Temp.

[K]

Outlet 

Mod. 

Density

[kg/m3]

Outlet 

Mod.

Temp.

[K]

ANM 2G 1679 ref ref 836.0 706.1 581.3 662.1 598.8

SENM 2G 1678 0.00 0.02 836.2 706.1 581.3 662.1 598.8

NEM 4G 1674 0.31 0.50 836.1 706.1 581.3 662.1 598.8

SENM 4G 1675 0.03 0.34 836.2 706.1 581.3 662.1 598.8

NEM 8G 1672 0.55 0.86 836.2 706.1 581.3 662.1 598.8

SENM 8G 1673 0.39 0.69 836.2 706.1 581.3 662.1 598.8

The SENM kernel appears to determine properly the axial power distributions in Figure 7 but 

noticeable differences are found near the top and bottom axial reflecting regions as shown in
26

Figure 8. The axial power distributions of the SENM kernel show about 7 ~ 8 % differences from 

those of the ANM result near the reflectors, while they are quite similar near the center. This is 

because the assembly discontinuity factor is considered along the axial direction as well as the

radial direction in the SENM kernel but not in the existing kernels. From these results, it can be 

judged that the source expansion nodal routines have been integrated well with routines for T/H 

feedbacks and cross-section handling with interpolation in PARCS. The multi-group results are

also similar to the two-group results.

Figure 7 Axial power distribution

Figure 8 Axial power distribution error
27

Comparison of pin-power distributions are shown in Table 5 and Table 6. Comparing with the 

results of the ANM 2G kernel, the SENM kernel gave reasonable results with the two- and multigroup structures.

Table 5 Pin power %PWE at HFP

Nodal

Kernel

Assembly Position

(A,1) (B,2) (C,3) (D,4) (E,5) (F,6)

ANM 2G ref ref ref ref ref ref

SENM 2G 0.00 0.00 0.00 0.00 0.02 0.25

NEM 4G - - - - - -

SENM 4G 0.00 0.00 0.02 0.03 0.14 0.41

NEM 8G - - - - - -

SENM 8G 0.10 0.01 0.03 0.33 0.18 0.39

Table 6 Pin power %EWE at HFP

Nodal

Kernel

Assembly Position

(A,1) (B,2) (C,3) (D,4) (E,5) (F,6)

ANM 2G ref ref ref ref ref ref

SENM 2G 0.35 0.35 0.35 0.35 0.35 5.00

NEM 4G - - - - - -

SENM 4G 0.36 0.35 0.36 0.38 0.40 4.60

NEM 8G - - - - - -

SENM 8G 0.37 0.37 0.36 0.47 0.45 4.05

4.3

Part III of the benchmark problem report describes the steady-state at hot zero power condition. 

The critical boron concentration and assembly power distributions results are shown in Table 7. 

Since the errors due to the cross-section interpolation and feedback handling are negligible at 

HZP conditions, the error are mainly ascribable to the spatial discretization. The assembly 

power errors become slightly larger than those of Part I because the effect of spatial 

discretization increases in three-dimensional calculations. But all of the nodal kernels show

good agreements with the reference solution. 

Part III – 3D Calculation at Hot Zero Power Conditions
28

Table 7 Comparison of critical boron concentration and assembly power

Nodal

Kernel

CBC

[ppm]

*Assembly Power Error

%PWE %EWE

ANM 2G 1341 1.05 3.49

SENM 2G 1340 1.04 3.49

NEM 4G 1337 1.11 3.06

SENM 4G 1337 1.09 3.22

NEM 8G 1334 1.20 2.85

SENM 8G 1335 1.14 2.91

*Assembly Power Error: the DeCART result is used as the reference

Axial power distributions and the relative errors at each axial location are shown in Figures 9

and 10, respectively. Since the axial reflector model has a large effect on determining the axial

power distribution, the same axial reflector model is used for all cases of calculations in this 

benchmark problem. The axial power distributions generated by SENM kernel are in a good 

agreement with the reference results in all axial regions, whereas the existing nodal kernels of 

PARCS show large relative errors of about 6 ~ 10 % in the top and bottom fuel regions. The 

reason why the existing nodal kernels show such large errors near the reflector regions seems

that the assembly discontinuity factor is not considered in the axial direction. From these results, 

it can be concluded that the assembly discontinuity factor should be also considered in the axial 

direction.

Table 8 and Table 9 show the results of pin power calculations. In the three-dimensional pin 

power calculation, two-dimensional pin power calculation is performed several times for each 

plane and the radial pin power distribution is obtained by averaging the plane-wise pin powers 

over the axial direction. The pin power results show similar trends with the results of Part 1 for 

different nodal kernels and group structures. In this calculation, it was worthwhile to note that 

the source expansion pin power reconstruction works well for three-dimensional calculations as 

well as two-dimensional calculations.
29

Figure 9 Axial power distribution at HZP conditions

Figure 10 Axial power distribution error at HZP conditions
30

Table 8 Pin power %PWE at HZP conditions

Nodal

Kernel

Assembly Position

(B,2) (D,4) (F,6)

ANM 2G 0.82 0.62 0.48

SENM 2G 0.78 0.51 0.71

SENM 4G 0.78 0.46 0.72

SENM 8G 0.76 0.44 0.60

Table 9 Pin power %EWE at HZP conditions

Nodal

Kernel

Assembly Position

(B,2) (D,4) (F,6)

ANM 2G 1.83 1.39 0.90

SENM 2G 1.83 1.71 5.76

SENM 4G 1.88 1.85 5.39

SENM 8G 1.85 1.88 5.06

4.4

In Parts I ~ III, capabilities of the source expansion nodal and pin power kernels have been

evaluated and their performance was found superior for steady-state problems. In this section, 

the transient capability is evaluated by simulating a fast transient accident in which a single rod 

is ejected in 0.1 second. Like Part II, the transient results are not compared to the reference

results, but the results of SENM kernel are compared with those of existing nodal kernels 

qualitatively because any reference results are not given in Part IV of the benchmark report. The 

ejected rod worth is about 1.13$ so that it causes a super-prompt critical transient involving a 

sharp power pulse from the hot-zero-power condition.

Part IV – 3D Transient Calculation

In Table 10, peak time, peak power and power integral are compared and Figure 11 shows

transient core power behaviors. From these comparisons, one can conclude that the two-group 

results calculated by the SENM kernel are essentially identical to the results of the ANM kernel

and the SENM kernel provides the same accuracy as the analytic solution in transient

calculations as well as steady-state calculations. In case of the multi-group results, the SENM 

and the NEM kernels show a slight difference of about 10% in the peak power whereas they 

show coincidence in their peak times. It was found in a previous work regarding the source 

expansion methodology that the results of NEM kernel approach to those of the SENM kernel 

as number of nodes per fuel assembly is increased. The SENM kernel simulates onedimensional flux shape using the trigonometric functions and therefore provides better results

than the NEM kernel.
31

Table 10 Transient peak power and power integral

Figure 11 Transient core power behavior

Behaviors of transient assembly peaking and three-dimensional point pin peaking factors are

shown in Figure 12 and 13, respectively. While the existing routines were only able to calculate 

the two-group transient pin power distribution, the new routines are able to calculate it with all 

group structures as shown in the figures. The assembly and pin peakings increase quickly 

during the rod ejection period and then decreases slowly with T/H feedback. Especially, the 

multi-group transient calculation of the source expansion pin power reconstruction routines has

been verified by the fact that the pin peaking behavior of the multi-group calculation agrees well

with that of the existing two-group calculation.

Nodal

Kernel

Peak Time

[sec]

Peak Power

[%]

Power Integral

[%-sec]

ANM 2G 0.34 142 27.2

SENM 2G 0.34 146 27.3

NEM 4G 0.33 152 27.8

SENM 4G 0.33 163 28.4

NEM 8G 0.32 172 29.1

SENM 8G 0.32 180 29.8
32

Figure 12 Transient assembly peaking factor

Figure 13 Transient point pin peaking factor

The transient behavior of T/H properties such as the Doppler temperature, the moderator 

density and moderator temperature are shown in Figure 14, 15 and 16, respectively. Since the

same T/H routines are used with similar power behaviors, these are also in good agreements.
33

Figure 14 Transient core average Doppler temperature

Figure 15 Transient core average moderator density
34

Figure 16 Transient core average moderator temperature

4.5

One of the main reasons to develop efficient multi-group nodal method is to enhance the 

computational performance of the existing NEM kernel in which one-node nodal expansion 

method is accelerated by two-group CMFD formulation. The computing times taken for 

calculations of Parts II and IV are compared to assess the computational performance of the 

SENM kernel. The computing environment is as follows. 

Computing Time Performance

- CPU: Intel Xeon CPU E5520 2.27Ghz X 2

- MEMORY: DDR3 10600 12.0GB

- OPERATING SYSTEM: Windows 7 Professional 64 Bit Edition

Table 11 shows the iteration information of each multi-group nodal kernel. Compared with the 

existing kernels, the SENM kernel shows much faster computing performance than the NEM 

kernel. The largest difference in the computing performance is found in the nodal calculation. 

The NEM kernel takes computing time more than 10 times that of the SENM kernel. Since

SENM 8G takes about 1.04 second per iteration which is about two times 0.46 second of SENM 

4G, the computing time seems to be proportional to the number of energy groups in case of 

SENM kernel. However, the computing time of the NEM kernel is proportional to the square of 

energy groups as 3.15 seconds of NEM 4G is increased to 11.27 seconds for NEM 8G. By 

these results, the SENM kernel shows a superior computing performance especially for multigroup calculations. 

Computing performance of the nodal kernels becomes more important in transient problems

because they take much more iterations and time. 

Table 12 compares elapsed times taken for Part IV calculations. The calculation time differ by

about 5000 seconds in case of the eight-group structured problem which should be a 

considerable advantage in analyzing reactor core transients. 
35

Table 11 Comparison of computing times for the steady-state calculation (Part II)

Nodal

Kernel

NODAL CMFD XSEC MISC. TOTAL

Iter. # Time (sec) Iter. # Time (sec) Time (sec) Time (sec) Time (sec)

NEM 4G 15 47.31 44 6.06 3.94 0.22 57.53

SENM 4G 9 4.11 49 3.36 1.59 0.23 9.30

NEM 8G 12 135.29 - 5.01 5.89 0.24 146.44

SENM 8G 8 8.33 47 6.72 4.19 0.19 19.42

Table 12 Comparison of computing time for the transient calculation (Part IV)

Nodal

Kernel

NODAL+CMFD OTHERS TOTAL

Time (sec) Time (sec) Time (sec)

NEM 4G 1696.6 221.0 1917.6

SENM 4G 473.1 178.4 651.5

NEM 8G 6306.8 390.0 6696.8

SENM 8G 957.4 497.9 1455.3

37

5. CONCLUSIONS

An efficient multi-group solution based on source expansion nodal method consisting of nodal 

and pin power calculation schemes has been formulated successfully for multi-group problems. 

The source expansion formulation is to expand source terms located on the right hand side of 

the neutron balance equation with a set of polynomial functions and provides the decoupled 

groupwise neutron balance equations which can be solved more easily than the coupled 

equations. Source iteration is essentially needed to achieve more exact solutions. 

The source expansion nodal scheme solving the transverse-integrated one-dimensional neutron

balance equation is coupled with the multi-group CMFD formulation to find the global neutron 

flux distribution. And the two-level CMFD scheme in which two-group CMFD calculation 

accelerates the multi-group CMFD calculation is also employed. All of the modules in cmfd

category provide subroutines and variables for the two-level CMFD calculation and the source 

expansion nodal method is implemented in the senm category. In addition, the BICG 

preconditioned by incomplete LU factorization is implemented in the bicg2g and bicgmg 

categories to solve the CMFD linear systems. The pin power reconstruction with the source 

expansion approach is implemented into the ppr category with consideration of the corner 

discontinuity factor.

The numerical coding is written in the form of a static library of FORTRAN90 and all subroutines 

of the library except the main category is encapsulated from PARCS. This programming strategy 

can make the library independent from the PARCS code and allows easy maintenance of 

routines. Also, it minimizes modification of the PARCS code.

In order to evaluate the newly implemented routines in terms of accuracy and speed, the 

OECD/NEA AND U.S. NRC PWR MOX/UO2 CORE TRANSIENT BENCHMARK Problems was

chosen as the reference. Throughout all parts of the benchmark problem containing steadystate and transient cases, the new source expansion kernel provides enhanced accuracy 

compared to the existing nodal expansion kernel. Especially, the multi-group pin power which is 

not obtainable with the original kernel has become available with the new source expansion 

kernel and agrees well with the reference. The source expansion kernel is about 4 times faster 

than the nodal expansion kernel in the 3D real-size core problems. From the results, it is 

concluded that the newly implemented source expansion kernel provides a good accuracy, yet 

very fast calculation speed regardless the structure of energy groups chosen.

39

6. REFERENCES

1. H. G. Joo, D. A. Barber, G. Jiang and T. J. Downar, “PARCS, A Multi-Dimensional Two

Group Reactor Kinetics Code Based on the Nonlinear Analytic Nodal Method,” PU/NE-98-

26, Sept. 1998.

2. P.D. Esser and K.S. Smith, "A Semi-Analytic Two-Group Nodal Model for SIMULATE-3," 

Trans. Am. Nucl. Soc., 68, 220, 1993.

3. V. G. Zimin, H. Ninokata and L. R. Pogosbekyan, "Polynomial and Semi-Analytic Nodal 

Methods for Nonlinear Iteration Procedure," Proc. Int. Conf. Phys. Nucl. Sci. Tech. 

(PHSOR-98), Long Island, NY, October 5-8, 1998, 2, 994, 1998.

4. V. G. Zimin and H. Ninokata, "Nodal Neutron Kinetics Model Based on Nonlinear Iteration 

Procedure for LWR Analysis," Ann. Nucl. Energy, 25, 507, 1998.

5. Y. I. Kim, Y. J. Kim, S. J. Kim and T. K. Kim, "A Semi-Analytic Multigroup Nodal Method," 

Ann. Nucl. Energy, 26, 699, 1999.

6. X. D. Fu and N. Z. Cho, “Nonlinear Analytic and Semi-Analytic Nodal Methods for 

Multigroup Neutron Diffusion Calculations,” J. Nucl. Sci. Tech., 39, 1015, 2002.

7. D. J. Lee, T. J. Downar and Y. H. Kim, " Convergence Analysis of the Nonlinear CoarseMesh Finite Difference Method for One-Dimensional Fixed-Source Neutron Diffusion 

Problem," Nucl. Sci. Eng., 147, 127, 2004.

8. T. Y. Han, H. G. Joo and C. H. Kim, "Two-Group CMFD Accelerated Multi-Group 

Calculation with a Semi-Analytic Nodal Kernel," Proc. PHYSOR 2006, ANS Topical Mtg. on 

Reactor Physics, Vancouver, B.C., Canada, Sept. 10-14, 2006.

9. T. Kozlowski and T. Downar, "Pressurized Water Reactor MOX/UO2 Core Transient 

Benchmark, Final Report," OECD Nuclear Energy Agency, 2006.

10. Böer, R and Finnemann, H., “Fast Analytical Flux Reconstruction Method for Nodal SpaceTime Nuclear Reactor Analysis,” Annl. Nucl. Energy, 19, 617, 1992.

11. J. I. Yoon and H. G. Joo, “Two-Level Finite Difference Formulation with Multigorup Source 

Expansion Nodal Kernel,”Journal of Nuc. Sci. and Tech., 45, 668, 2008.




UNITED STATES 

NUCLEAR REGULATORY COMMISSION

WASHINGTON, DC 20555-0001

--------------------

OFFICIAL BUSINESS
NUREG/IA

-0416 Implementation of Advanced Multigroup Nodal and Pin Power 

Reconstruction Methods into PARCS 3.1