Science and Technology of Nuclear Installations

Volume 2017, Article ID 7072197, 10 pages

https://doi.org/10.1155/2017/7072197

## Development and Application of a New High-Efficiency Sparse Linear System Solver in the Thermal-Hydraulic System Analysis Code

^{1}Xi’an Jiaotong University, 28 W. Xianning Rd, Xi’an, China^{2}Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, China

Correspondence should be addressed to Jianqiang Shan; nc.ude.utjx.liam@nahsqj

Received 7 April 2017; Revised 3 July 2017; Accepted 22 August 2017; Published 19 September 2017

Academic Editor: Manmohan Pandey

Copyright © 2017 Li Ge et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This paper presents a faster solver named NRLU (Node Reordering Lower Upper) factorization solver to improve the solution speed for the pressure equations, which are formed by RELAP5/MOD3.3. The NRLU solver uses the oriented graph method and minimal fill-ins rule to reorder the structure of the nonsymmetry sparse pressure matrix. It solves the pressure matrix by LU factorization. Then the solver is embedded into the large scale advanced thermal-hydraulic system analysis program RELAP5/MOD3.3. The comparisons of the original solver and the NRLU solver show that the NRLU solver is faster than the original solver in RELAP5/MOD3.3, and the rate enhancement can be 44.44%. The results also show that the NRLU solver can reduce the number of fill-ins effectively. This can improve the calculation speed.

#### 1. Introduction

After entering the 21st century, the requirement of the safety and economy performance for nuclear power plant has been raised. Researchers proposed different novel concepts of advanced nuclear power systems, such as the small module reactor and the generation IV reactors. Accurate and fast simulation of these systems’ detailed behavior under transient conditions has become the key issue.

For simulating the large or complex systems using a personal computer, the existing nuclear system analysis code would spend several days to get the results. And it will take more time to modify or design a modeling system repeatedly. The future nuclear system power programs require the more fine or complex modeling, which could cause a long calculation time. So it is necessary to develop a more efficient numerical technique to improve the calculation efficiency of the nuclear system analysis code. On the other hand, the faster calculation speed is important to achieve the real-time requirement for the nuclear power plant simulators.

To improve the calculation speed, the nuclear system analysis code must be equipped with a faster unsymmetrical sparse matrix solver for the system equations, which could cost nearly half the time during the total calculation. However, the research in improving the matrix solver for nuclear system analysis code is rare. At present, the existing nuclear system analysis codes are all equipped with the matrix solvers developed at least ten years ago. RELAP5/MOD3.3 code is equipped with the matrix solver developed by Curtis and Reid in 1971. CATHARE code uses the normal Newton iteration method. These matrix solvers are universal and are not specially developed for the nuclear system analysis code. Therefore, a more efficient matrix solver for the nuclear system analysis code needs to be developed.

Over the last several decades, computer speed and memory have increased dramatically. Many methods have been developed for solving the large unsymmetric, sparse linear matrices. Standard direct methods based on the Gaussian elimination require more work than iterative schemes. As a result, such system matrices are typically solved by iterative methods with fast algorithms, such as GMRES [1] and Bi-CGSTAB [2]. For the integral equations of classical physics, iterative scheme has led to some of the fastest solvers known today. However, iterative methods still have some significant disadvantages compared with direct methods: () The iteration number for an iterative solver is highly sensitive to the conditioning of the system matrix. () Iterative methods cannot solve a linear system governed by a fixed matrix with multiple right hand sides. While direct methods can be applied to solve this kind of system at a much lower cost [3], iterative methods are very useful for a large sparse system. But for small or medium sparse matrix, iterative methods cannot do better than direct methods. In the thermal-hydraulic system analysis code, the matrices are mostly small or medium. So, to solve the system matrix of the thermal-hydraulic system analysis code, direct methods are better than interactive methods.

Standard direct methods based on Gaussian elimination require the computation work of , where is the system size. This means that this method will become infeasible while is increasing. Therefore, the coefficient matrix needs to be reordered to reduce the new nonzero elements storage and floating point operation. So far, the matrix reorder algorithms can be classified into two groups by their optimization goals. One group is to reduce the bandwidth and profile of matrix, such as the RCM algorithm, the GPS algorithm [4], the Schonauer algorithm, the King algorithm, and the Sloan algorithm. The other is to reduce the fill-ins [5], such as the MD algorithm [6], the AMD algorithm [7], and the Nested Dissection algorithm [8].

The original matrix solver in RELAP5/MOD 3.3 is a direct method based on Markowitz algorithm. The concept of the Markowitz algorithm [9] is choosing the element with minimum Markowitz factor as the pivot at each step of the factorization. If the number of nonzero elements in row is and the number of nonzero elements in column is , the Markowitz factor of element () is as follows:The solution process in RELAP5/MOD3.3 is performed in three steps:(1)First, symbolic factorization to determine the maximum memory requirements using Markowitz ordering(2)Second, numeric factorization using Markowitz ordering with threshold pivoting(3)Third, solution of the LU systems.A faster direct solver with the minimal fill-ins preprocessing will be presented in this paper, which is developed to solve sparse linear systems resulting from the application of the discretized two phase hydrodynamics equations for nuclear reactor transient problems. This solver is a direct method based on lower upper factorization with a structure reordering processing. Before numerical LU factorization, a nodes reordering is processed during symbolic factorization according to its data structure to reduce the fill-ins. This solver, named NRLU, will improve the calculation speed compared with the standard direct method and hold the advantages of the direct methods. The NRLU solver is then embedded into the large scale advanced thermal-hydraulic system analysis program RELAP5/MOD3.3. Some transients for reactor power plant systems were analyzed using RELAP5/MOD3.3 program with its original matrix solver and the NRLU solver, respectively.

#### 2. The NRLU Solver Algorithm

##### 2.1. Pressure Matrix in RELAP5/MOD3.3

In RELAP5/MOD3.3, the equations of the noncondensable density equation, the vapor/gas energy equation, the liquid energy equation, the difference density equation, and the sum density equation are firstly eliminated to obtain an equation that involves only the unknown variables , , , , and . Then substituting the velocity equations into this equation, this will result in a single equation involving only the pressure parameter. This is done for each volume, giving rise to an system of linear equations for the new pressures in a system containing control volumes.

The linear equation is formed asThe coefficient matrix will be larger as the number of the control volumes increases. And the sparsity of will increase because the problems of the nuclear system are almost 1D or 2D in RELAP5/MOD3.3.

During the long-term transient analyses of the nuclear systems, (2) will be solved more than a million times. So, every unnecessary nonzero operation and nonzero storage need to be avoided. Therefore the coefficient matrix needs to be reordered before factoring (1) to reduce the operation times and element storages during the matrix factorization. So far, there are two basic principles for reorder: () minimal fill-ins and () minimal long operations. In this paper, the minimal fill-ins will be the basic principle to renumber the nodes to improve the calculation speed of the sparse linear matrix.

Because of the existence of the time dependent volumes, the pressure coefficient matrix is almost nonsymmetric. At present, one of the popular methods is to regard the nonsymmetric structure as symmetric. However, this method needs more computer memory and wastes calculation time. In this paper, a method aimed at solving the nonsymmetric pressure matrix formed from the RELAP5/MOD3.3 code is proposed based on the oriented graph method. The key technology of this method is the node ordering optimization for the control volumes of the simulation system [10].

##### 2.2. Matrix Represented by Oriented Graph

Graph theory is a fundamental tool in sparse matrix techniques. The nonzero pattern and the relationship between the nonzero elements can be represented by the oriented graph.

Assuming that the order of the coefficient matrix is and the elements of main diagonal are nonzero, there are nodes to represent the row numbers (or column numbers) of the matrix . The nonzero elements of the matrix can be represented by the directed line between two nodes. For example, for , a directed line can be used from node to node to represent the nonzero element in row and column . The nonzero elements of the main diagonal can be represented by each node itself [11].

For example, the oriented graph for the nonsymmetric matrix (3) is shown in Figure 1.If there are two directed lines between two nodes, these two directed lines can be replaced by an undirected line. Figure 1 can be simplified as Figure 2. Obviously, the graph is undirected if matrix has symmetric structure.