IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes

In the extreme violence of merger and mass accretion, compact objects like black holes and neutron stars are thought to launch some of the most luminous outbursts of electromagnetic and gravitational wave energy in the Universe. Modeling these systems realistically is a central problem in theoretical astrophysics, but has proven extremely challenging, requiring the development of numerical relativity codes that solve Einstein's equations for the spacetime, coupled to the equations of general relativistic (ideal) magnetohydrodynamics (GRMHD) for the magnetized fluids. Over the past decade, the Illinois Numerical Relativity (ILNR) Group's dynamical spacetime GRMHD code has proven itself as a robust and reliable tool for theoretical modeling of such GRMHD phenomena. However, the code was written "by experts and for experts" of the code, with a steep learning curve that would severely hinder community adoption if it were open-sourced. Here we present IllinoisGRMHD, which is an open-source, highly-extensible rewrite of the original closed-source GRMHD code of the ILNR Group. Reducing the learning curve was the primary focus of this rewrite, with the goal of facilitating community involvement in the code's use and development, as well as the minimization of human effort in generating new science. IllinoisGRMHD also saves computer time, generating roundoff-precision identical output to the original code on adaptive-mesh grids, but nearly twice as fast at scales of hundreds to thousands of cores.Comment: 37 pages, 6 figures, single column. Matches published versio

Numerical Models of Binary Neutron Star System Mergers. I.: Numerical Methods and Equilibrium Data for Newtonian Models

The numerical modeling of binary neutron star mergers has become a subject of much interest in recent years. While a full and accurate model of this phenomenon would require the evolution of the equations of relativistic hydrodynamics along with the Einstein field equations, a qualitative study of the early stages on inspiral can be accomplished by either Newtonian or post-Newtonian models, which are more tractable. In this paper we offer a comparison of results from both rotating and non-rotating (inertial) frame Newtonian calculations. We find that the rotating frame calculations offer significantly improved accuracy as compared with the inertial frame models. Furthermore, we show that inertial frame models exhibit significant and erroneous angular momentum loss during the simulations that leads to an unphysical inspiral of the two neutron stars. We also examine the dependence of the models on initial conditions by considering initial configurations that consist of spherical neutron stars as well as stars that are in equilibrium and which are tidally distorted. We compare our models those of Rasio & Shapiro (1992,1994a) and New & Tohline (1997). Finally, we investigate the use of the isolated star approximation for the construction of initial data.Comment: 32 pages, 19 gif figures, manuscript with postscript figures available at http://www.astro.sunysb.edu/dswesty/docs/nspap1.p

Parareal in time 3D numerical solver for the LWR Benchmark neutron diffusion transient model

We present a parareal in time algorithm for the simulation of neutron diffusion transient model. The method is made efficient by means of a coarse solver defined with large time steps and steady control rods model. Using finite element for the space discretization, our implementation provides a good scalability of the algorithm. Numerical results show the efficiency of the parareal method on large light water reactor transient model corresponding to the Langenbuch-Maurer-Werner (LMW) benchmark [1]

Asymptotic Neutronic Solutions for Fast Burst Reactor Design

Deterministic numerical methodologies for solving time-eigenvalue problems are valuable in characterizing the inherent rapid transient neutron behavior of a Fast Burst Reactor (FBR). New nonlinear solution techniques used to solve eigenvalue problems show great promise in modeling the neutronics of reactors. This research utilizes nonlinear solution techniques to solve for the dominant time-eigenvalue associated with the asymptotic (exponential) solution to the neutron diffusion and even-parity form of the neutron transport equation, and lays the foundation for coupling with other physics phenomena associated with FBRs. High security costs and proliferation risks associated with Highly Enriched Uranium (HEU) fueled FBRs are the motivation for this research. Use of Low Enriched Uranium (LEU) as fuel reduces these risks to acceptable levels. However, the use of LEU fuel introduces complexities such as, increased volume, and longer neutron lifetimes. Numerical techniques are sought to explore these complexities and determine the limitations and potential of a LEU fueled FBR. A combination of deterministic and stochastic computational modeling techniques are tools used to investigate the effects these complexities have on reactor design and performance. Monte Carlo N-Particle (MCNP) code is useful to determine criticality and calculate reactor kinetics parameters of current and proposed designs. New deterministic methods are developed to directly calculate the fundamental time-eigenvalue in a way that will support multi-physics coupling. The methods incorporate Jacobian Free Newton Krylov solution techniques to address the nonlinear nature of the neutronics equations. These new deterministic models produce data to determine LEU designs that may meet the performance requirements of proven HEU FBRs in terms of neutron burst yield and burst duration (pulse width) based on the Nordheim-Fuchs model. This computational data and measured performance characteristics of historical LEU FBRs show that LEU designs can generate pulses that are beneficial for meeting Research and Development (R&D) requirements. These modern computational neutronic results indicate that a LEU fueled FBR is a plausible alternative to current HEU fueled reactors

직접 전노심 과도해석능 고도화

학위논문(박사) -- 서울대학교대학원 : 공과대학 에너지시스템공학부, 2022.2. 주한규.원자로 규제 기준이 강화되고, 고정밀도 다물리 연계계산에 대한 수요가 증가하면서 고신뢰도 전노심 직접 과도해석이 요구되는 상황이다. 그러나 전노심 직접 계산의 막대한 계산요구량 때문에 실제적인 노심 문제 해석에 사용하는 경우 많은 계산 시간을 필요로 하거나, 수천 코어 수준의 대규모 컴퓨팅 시설에 의존해야 한다는 한계를 보였다. 본 연구는 GPU 컴퓨팅 기술 적용 및 과도해석 방법론 개선을 통해서 효율적인 과도해석능을 전노심 직접해석 코드 nTRACER에 구현하는 것을 목표로 한다. nTRACER의 삼차원 직접 전노심 수송해석은 이차원 층별 특성곡선법 (MOC) 삼차원 소격격자 유한 차분법 (CMFD), 일차원 축방향 특성곡선법 등의 계산요소들의 연계를 통해서 이루어진다. 본 연구에서는 기존 해법에 점근사 동특성 방정식 (PKE)을 도입하여 노심 전체 중성자속의 거동을 해석하는 준정적 해법 (Quasi-static method)을 도입하였다. 이를 통해 MOC/CMFD/PKE로 이루어진 3단계 해석체계를 구현하고 각 단계마다 다른 시구간 크기를 적용하였다. 요구되는 계산량이 비교적 크고 계산을 통해서 결정하는 변수의 크기의 변화율이 작은 MOC와 CMFD 계산에는 비교적 큰 시구간을 사용하고, 계산량이 작은 PKE 계산에는 작은 시구간을 사용함으로써 계산량 대비 높은 정확도를 얻을 수 있다. 과도 상황에서 시간에 따라 각 변수의 변화율 또한 변화하기 때문에 적응형 해법을 구현하여 불필요한 MOC 및 CMFD 계산을 줄이고 정확도 대비 요구되는 계산량을 최소화하였다. MOC 계산의 경우 노심 조건이 크게 변화하는 시구간에서만 계산을 수행하는 조건적 수송계산 해법을 통해서 적응형 해법을 구현하였다. 특히 본 연구에서는 조건적 MOC 발동 기준을 기존의 단순한 일군 반응단면적 변화 기준 대신 세부 격자 잔차항 기준을 도입하여 MOC 계산의 영향을 더 정교하게 평가하였다. CMFD 계산의 시구간 크기는 적응형 시구간 조정 알고리즘을 도입하여 각 시구간에서 발생하는 에러 값이 주어진 허용치 이하로 유지되도록 조정하였다. 해당 알고리즘을 위해서 각 시구간에서 발생할 수 있는 오차 모델을 유도했고, 유도된 오차 모델을 통해 과도 계산 중 발생하는 오차를 실시간으로 추정하여 허용치 기준을 만족하는 시구간 크기를 산출한다. 5×5 핵연료 집합체 문제 해석을 통해서 적응형 해법들을 검증하였다. 새로운 MOC 발동 기준은 이전의 MOC 발동 기준보다 최대 출력 상대 오차 값을 약 80 % 감소시켰다. 적응형 시구간 조정 알고리즘 도입 결과, 검증문제 계산 시 출력이 빠르게 변화하는 구간에서 발생하는 오차가 주어진 허용치 이하로 유지되었고, 같은 수의 시구간을 사용한 고정 시구간 결과 대비 최대 출력 상대 오차가 약 80 % 감소하였다. nTRACER 과도해석 요소 중 MOC 및 CMFD의 선형계 해법 등 연산 집약적인 요소들에 GPU 컴퓨팅이 적용되었다. GPU의 특성은 기존의 CPU와 다르기 때문에 이에 맞게 최적화가 이루어졌다. 특히 본 연구에서는 과도 CMFD 선형계 해법의 최적화가 중점적으로 수행되었다. 우선, 대규모 데이터 처리에 적합한 GPU 특성에 맞지 않는 그룹 우선 배치 방식의 선형계 해법 대신 다중 그룹 직접 해법을 적용하여 수렴 안정성과 계산 속도를 향상시켰다. 또한, CMFD 선형계 해법에서 사용되는 선조건자를 기존의 불완전 LU 분해 기반 선조건자 대신 대규모 병렬 실행이 가능한 희소 근사 역행렬 (SPAI) 선조건자로 대체하였다. SPAI 선조건자는 CMFD 행렬의 역행렬을 미리 정해진 희소행렬 구조에 따라서 근사한다. 희소행렬 구조에 따라서 SPAI 선조건자의 생성 비용과 선형계 해법 반복계산수가 결정되기 때문에 최적의 희소행렬 구조가 요구된다. 본 연구에서는 SPAI 추정 알고리즘을 도입하여 CMFD 행렬에 따른 희소행렬 구조 최적화를 수행하였다. 추정 알고리즘 도입 결과 2 %의 추정 탈락 기준치를 사용했을 때, 기존의 고정 희소행렬 구조 대비 요구되는 선형계 해법 반복계산수가 13 % 감소하였다. 새로운 직접 전노심 과도해석능의 유효성을 다양한 실제적인 노심 문제 해석을 통해서 입증하였다. SPERT III E-core 반응도 사고 실험을 해석하고 이를 실험치와 비교하여 새로운 과도해석능의 중성자 거동 해석, 동적 데이터 처리, 열궤환 모델의 정확성을 확인하였다. 서로 다른 조건의 5가지 대표 문제에 대해서 최대 출력, 노심 주기, 에너지 방출량 등 주요 실험치와 실험 불확실도 범위 이내에서 일치하였다. 같은 조건에서 수행한 이단계 해석방법을 적용 시 노심 주기에서 실험 불확실도보다 큰 오차를 보여, 직접 전노심 과도해석능을 정확도를 확인하였다. 각 계산은 20개의 상용 GPU를 장착한 클러스터에서 7시간 이내에 수행되었고, 320개의 CPU를 장착한 클러스터에서 CPU 기반 nTRACER 과도해석능을 사용한 계산보다 약 7배 빠른 속도를 보였다. 각 계산에 사용된 컴퓨팅 시설의 가격과 시간을 고려했을 때, 새로운 과도해석능은 기존 과도해석능 대비 약 13 배 높은 가격대비성능을 가진다고 볼 수 있다. APR1400에서 발생하는 가상의 제어봉 이탈 사고를 24개의 GPU를 장착한 클러스터를 사용하여 계산하였고, 125 개의 시구간을 통한 1초 동안의 과도해석을 19시간 이내에 수행하여 직접 전노심 과도해석의 실용적 활용에 대한 가능성을 확인하였다.As regulations for nuclear reactors have been strengthened, and demand for high-fidelity multi-physics simulation has increased, high-fidelity direct whole core transient calculation is required. However, due to its heavy computational burden, application of direct whole core transient calculation to realistic core problems requires too long computing time or leadership class computing facilities containing thousands of CPU cores. In this work, the efficiency of direct whole core transient capability of nTRACER is enhanced through applying GPU computing technique and improvement of methodology. The 3D direct whole core transport calculation of nTRACER is composed of 2D planar method of characteristics (MOC), 3D coarse mesh finite difference (CMFD), and 1D axial MOC calculations. In this work, quasi-static method is employed which solve the point kinetics equation (PKE) to simulate the temporal variation of overall amplitude of neutron flux. Consequently, 3-level method composed of MOC/CMFD/PKE which uses different time step size for each level is implemented in nTRACER. Relatively large time step sizes are used for MOC and CMFD calculations, whose computational burden is heavy and results vary slowly. On the other hand, small time step size is used for PKE calculation. As a consequence, the calculation burden of 3D direct whole core calculation can be alleviated without significant loss of accuracy. Since the temporal change rate of each variation vary over time, adaptive solution method is implemented to avoid unnecessary MOC and CMFD calculation and minimize the required cost to the accuracy. For MOC calculation, adaptive solution is implemented by performing MOC calculation only when the core conditions are change significantly. In this work, rather than using simple cross section change, fine mesh residual norm is used as a MOC invoking criterion, which can evaluate the effect of MOC update more precisely. The time step size for CMFD calculation is controlled by the algorithm that is designed to maintain the error occur at each step below prescribed tolerance. To implement the adaptive time step control algorithm, new error model is derived. By using this model, error occur at each time step is calculated during transient calculation and the calculated error is used to determine the time step size. 5 ×5 fuel assembly problem is used to evaluate the new adaptive solution methods. The new MOC invoking criteria reduces the relative peak power error by 80 % level when compared to previous criteria. The adaptive time step size control algorithm effectively controls the error below prescribed tolerance for the interval where the core power level changes rapidly. When compared to the fixed time step size case which has similar number of time step, the adaptive time step size control algorithm decreases the relative peak power error by 80 % level. GPU computing technique is applied for computationally intensive components of nTRACER such as MOC or CMFD. Since the characteristics of GPU are different from CPU, optimization is made accordingly. Especially, optimization of transient CMFD linear system solution is performed intensively in this work. First, rather than the group-major ordered linear system solution which is inappropriate for GPU computing, multi-group direct solution method is employed to enhance the stability and speed of the convergence. And, instead of ILU (Incomplete LU) preconditioner, SPAI (Sparse approximate inverse) preconditioner is used to utilize the massive parallelism of GPU. SPAI preconditioner is constructed by approximating inverse of the CMFD matrix using prescribed sparsity structure. Since the sparsity structure used for the construction determines the construction cost and the quality of precondition- ing, it is important to find the optimal sparsity structure. In this work, SPAI prediction algorithm is devised for optimization of sparsity structure. By using the SPAI prediction algorithm with 2 % drop criteria, the required iteration number is reduced by 13 % when compared to existing fixed sparsity structure. The effectiveness of new direct whole core transient capability is verified through various realistic core problems. Verification of neutron kinetics, kinetics data treatment, and thermal feedback model is performed using SPERT III E-core RIA (Reactivity initiated accident) experiments. For 5 representative tests, the calculated values for experimental data such as peak power level, reactor period, and released energy show good agreement within uncertainty range. When using conventional two step method to analyze same experiment, large difference between calculated value and experimental data occurs. All calculations are run on the cluster containing 20 commercial GPUs and are finished within 7 hours. It is 7 times faster than the calculation using CPU version of nTRACER on the cluster containing 320 GPUs. Considering the price of computing facilities for each calculation, the new transient capability is about 13 times cost effective than the previous CPU version. To check the possibility of actual use of direct whole core transient calculation for realistic core problem, hypothetical RIA in APR1400 core is used. The RUA is simulated up to 1 s using 125 time steps. The calculation is run on the cluster containing 24 GPUs and is finished within 19 hours.Contents Abstract i Contents v List of Figures viii List of Tables xi 1 Introduction 1 1.1 Study Background and Motivation 1 1.2 Objectives and Scopes 6 2 Direct Whole Core Transient Calculation Methodology 10 2.1 Time-dependent Neutron Transport Solutions 13 2.1.1 Time-dependent Planar MOC Solution 13 2.1.2 Time-dependent CMFD Solution 18 2.2 Effective Cross Section Generation 21 2.3 Kinetics Parameters Treatment 23 2.4 CMFD-based Adjoint Capability 28 2.5 Approximate Flux Weighting Method 31 3 The Multi-level Method 33 3.1 Intermittent Transport Update 35 3.2 Neutron Flux Factorization Methods 36 3.2.1 Improved Quasi-Static Method 36 3.2.2 Predictor Corrector Quasi-Static Method 39 3.2.3 Exponential Transform Method 39 3.3 Delayed Neutron Precursor Treatment 41 3.4 Examination of Flux Factorization Methods 43 3.4.1 C5G7-TD Results 43 3.4.2 SPERT III E-core Results 46 4 Adaptive time step Control 51 4.1 Conditional Transport Update 51 4.1.1 MOC Invoking Criteria 52 4.1.2 Evaluation of Flux Shape Change Estimator 55 4.2 Time Step Control of CMFD 59 4.2.1 Error Analysis of Multi-level Method 59 4.2.2 Estimation of Error 70 4.2.3 Determination of Time Step 73 4.2.4 Evaluation of Adaptive Time Step Size Control 75 4.3 Employment of Multi T/H Steps 82 5 Enhancement of CMFD Solution 87 5.1 Formulations of Transient CMFD 87 5.1.1 Group Major Ordering 87 5.1.2 Multi-group Direct Solution 91 5.1.3 Numerical Calculation Results 91 5.2 Preconditioner for Node Major Transient CMFD 95 5.2.1 Sparse Approximate Inverse Preconditioner 95 5.2.2 A Priori Sparsity Structure for SPAI preconditioner 97 6 Numerical Analyses 104 6.1 SPERT III E-core RIA Experiments 104 6.1.1 Calculation Options and Basic Information 104 6.1.2 Core Property Calculation at Zero Power Conditions 105 6.1.3 Analysis of the RIA simulation results 107 6.1.4 Computing Time Results 122 6.2 APR1400 Full Core Analysis 125 7 Conclusion 131 Bibliography 136 Appendix A SPERT III E-core Modelling 139 초록 149박

Low-Order Multiphysics Coupling Techniques for Nuclear Reactor Applications

The accurate modeling and simulation of nuclear reactor designs depends greatly on the ability to couple differing sets of physics together. Current coupling techniques most often use a fixed-point, or Picard, iteration scheme in which each set of physics is solved separately, and the resulting solutions are passed between each solver. In the work presented here, two different coupling techniques are investigated: a Jacobian-Free Newton-Krylov (JFNK) approach and a new methodology called Coarse Mesh Finite Difference Coupling (CMFD-Coupling). What both of these techniques have in common is that they are applied to the low-order CMFD system of equations. This allows for the multiphysics feedback effects to be captured on the low-order system without having to perform a neutron transport solve.The JFNK and CMFD-Coupling approaches were implemented in the MPACT (Michigan Parallel Analysis based on Characteristic Tracing) neutron transport code, which is being developed for the Consortium for Advanced Simulation of Light Water Reactors (CASL). These methods were tested on a wide range of practical reactor physics problems, from a 2D pin cell to a massively parallel 3D full core problem. Initially, JFNK was implemented only as an eigenvalue solver without any feedback enabled. However this led to greatly increased runtimes without any obvious benefit. When multiphysics problems were investigated with both JFNK and CMFD-Coupling, it was concluded that CMFD-Coupling outperformed JFNK in terms of both accuracy and runtime for every problem. When applied to large full core problems with multiple sources of strong feedback enabled, CMFD-Coupling reduced the overall number of transport sweeps required for convergence