A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0

JMCT is a large-scale, high-fidelity, three-dimensional general neutron–photon–electron–proton transport Monte Carlo software system. It was developed based on the combinatorial geometry parallel infrastructure JCOGIN and the adaptive structured mesh infrastructure JASMIN. JMCT is equipped with CAD modeling and visualizes the image output. It supports the geometry of the body and the structured/unstructured mesh. JMCT has most functions, variance reduction techniques, and tallies of the traditional Monte Carlo particle transport codes. Two energy models, multi-group and continuous, are provided. In recent years, some new functions and algorithms have been developed, such as Doppler broadening on-the-fly (OTF), uniform tally density (UTD), consistent adjoint driven importance sampling (CADIS), fast criticality search of boron concentration (FCSBC) domain decomposition (DD), adaptive control rod moving (ACRM), and random geometry (RG) etc. The JMCT is also coupled with the discrete ordinate SN code JSNT to generate source-biasing factors and weight-window parameters. At present, the number of geometric bodies, materials, tallies, depletion zones, and parallel processors are sufficiently large to simulate extremely complicated device problems. JMCT can be used to simulate reactor physics, criticality safety analysis, radiation shielding, detector response, nuclear well logging, and dosimetry calculations etc. In particular, JMCT can be coupled with depletion and thermal-hydraulics for the simulation of reactor nuclear-hot feedback effects. This paper describes the progress in advanced modeling, high-performance numerical simulation of particle transport, multiphysics coupled calculations, and large-scale parallel computing.


Introduction
Today, numerical simulations allow us to obtain as much detailed information as possible for any nuclear system based on the Monte Carlo (MC) method and software. The physical quantities can easily obtain both measurable and immeasurable values, especially for some extreme conditions inaccessible by experiments. Monte Carlo methods and software have been widely used in nuclear science engineering, statistical physics, biomedicine, quantum mechanics, molecular dynamics, petroleum geophysical exploration, finance, information, operational research, and polymer chemistry.
Joint Monte Carlo Transport (JMCT) has been developed by the IAPCM and CAEP-SCNS. It is designed to simulate neutrons, photons, electrons, proton, light radiation, atmosphere and molecule transport. This passes the validation and verification of a large number of benchmarks and experiments. The simulation problems cover the standard, critical fission, and adjoint sources. The data libraries use the continuous point-wise cross section of the ACE format and multigroup cross section of the ANISN This work was supported by the National Natural Science Foundation of China (Nos. 11805017 and 12001050).
format. Abundant variance reduction techniques have been incorporated into the JMCT. The JMCT software has been developed based on the parallel infrastructure JCOGIN [1] and JASMIN [2], where JCOGIN is a parallel combinatorial geometry infrastructure and JASMIN is an adaptive structured mesh infrastructure. The input uses CAD modeling, and the calculated results are outputted as a visual format or text file. Some problems in which the memory exceeds the limit of a single core or node can be simulated using the DD algorithm. At present, the number of geometric bodies, materials, tallies, depletion zones, sample numbers, period of random numbers, and parallel processors are large enough and expandable. In particular, the JMCT is coupled depletion and thermal hydraulic, where the depletion supports the TTA and CRAM methods, and thermal hydraulic supports the sub-channel method [3].
In this paper, first, the basic features and key algorithms are introduced. Then, the numerical results of the BEAVRS [4], VERA [5], AP1000 [6] and CAP1400 [7] reactor models are presented. For the BEAVRS, VERA, and AP1000 models, the numerical results of the JMCT are in good agreement with the results of MC21 [8], OpenMC [9], KENO-VI [10] and experiments. Finally, the initial physical parameters of the CAP1400 first core were predicted. The high accuracy of the JMCT was demonstrated in comparison with other well-known Monte Carlo programs.

JMCT system
JMCT is a program of the Joint Particle Transport System (JPTS) package. The JPTS package includes three floors, where the top floor contains the data library NuDa produced by NJOY [11]. The middle floor contains four application programs: JSNT (S N ) [12], JMCT [3] (MC), JBURN [13] (deletion), and JTH (thermal-hydraulic). The lower floor comprises the two support frameworks of JASMIN [2] and JCOGIN [1]. The unified CAD preprocessor JLAMT [14] and view postprocessor TeraVAP [15] are equipped for the JPTS (see Fig. 1). Figure 2 shows the structure flow of the JMCT, where the JMCT is coupled with the JSNT for radiation shielding and the JMCT is coupled with the JBURN and JTH for the reactor core.

JMCT modeling and geometry
The JMCT input is performed with the help of preprocess software JLAMT, which provides CAD modeling. Figure 3 shows the input interface. The input parameters include geometry, material, density, temperature, importance, tally types, source information, sample numbers, a wide variety of variance reduction techniques, energy cutoff, weight cutoff, and time cutoff, etc. The material in the geometric zone consists of multiple nuclei, including isotopes. Specifically, JMCT supports repetitive structures of the same geometry zone with different materials (MCNP [16] only supports repetitive structures of the same geometry zone with the same material). This function is available for the tightly coupled of neutron transport and depletion. The JMCT tally assorts a global tally and user-specified tally. By default, a global tally is performed for all zones according to the track length estimator. In addition, a mesh tally was added to the structured mesh.
The JMCT was developed on the basis of the JCOGIN infrastructure. JCOGIN supports the real combination geometry body and it integrates the common Monte Carlo particle transport methods, such as geometry descriptions, Boolean operations, track length calculations, random number generators, sample of arbitrary distributions, calculation of volumes and areas, domain decomposition (DD), parallel computations, and load balance technologies etc. At present, the JCOGIN infrastructure has supported the development of several Monte Carlo programs, and the JMCT is one of these programs.
JMCT uses real combinatorial geometry. The basic geometric bodies include cube, cylinder, sphere, cone, ellipsoid, and torus, etc. Some special bodies, such as reactor components, 1/n pillars, 1/n spheres, and special fuel assemblies, can be customized (see Fig. 4).

JMCT output
The JMCT supports visual output and text file output. The postprocessor TeraVAP supports image output in TBs scale. Figure 5 shows some examples.

Nuclear data
The nuclear data were from the ENDF/B-VII [17] and CENDL-3.2 [18] evaluation libraries. The cross section parameters were generated by using NJOY code [11]. The energy ranges are 10 -11 -20 MeV for neutrons, 1 eV-100 GeV for photons, and 10 eV-1 GeV for electrons, especially for 1 eV-1 keV of photons and 10 eV-1 keV of electrons are newly increased. It makes that JMCT has a simulation capability for visible light problems. The continuous point-wise cross section accounts for all the neutron reactions in the evaluation library. Thermal neutron processes in free-gas mode or S(a,b) mode. For photons, coherent and incoherent scattering, fluorescence emission, and electron pair generation are considered to follow photoelectric absorption. Nuclear data tables exist for neutron interactions, neutron-induced photons, photon interactions, neutron dosimetry or activation, and thermal scattering, S(a,b). Photon and electron data are atomic rather than nuclear in nature. The functions of proton, atmosphere and molecule transport are newly added in JMCT3.0.
Simultaneously, JMCT provides the function of multigroup calculation, where a 47-group neutron and 20-group photon P 5 cross section library is mainly used for radiation shielding, which its energy structure is same as the Bugle library [19]. Another 172-group neutrons, 40-group upper scatter neutrons, and 30-group photon P 5 library is mainly used for the reactor core.

Source specification
JMCT provides several standard sources. Users can specify a wide variety of source conditions. Independent probability distributions may be specified for the source variables of the position, energy, direction, time, and zone or surface. Information regarding the geometrical extent of the source can also be provided. In addition, the source variables may depend on other source variables. The user can bias all of the input distributions. For the energy, it includes various analytic functions for fission and fusion energy spectra, such as Watt, Maxwellian and Gaussian spectra. Biasing can also be accomplished by using special building functions. Simultaneously, the user source subroutine also provides a standard port.

Estimation of errors
Monte Carlo error e satisfies the center limit theorem as following x i is the mean, N is the sample number (or histories), a ¼ Eðx i Þ; i ¼ 1; 2; . . . is a mathematical expectation of various random x i , e ¼ kr ffiffi ffi N p x is an error, r is a variance and /ðkÞ ¼ 1 ffiffiffiffi Àk e Àt 2 =2 dt is the degree of confidence.
The FOM (a figure of merit) is calculated for one tally bin of each tally as a function of the sample numbers. It is defined as where t is the computation time (min). FOM value was used to determine the Monte Carlo calculation efficiency.   The larger FOM means the lesser the computer time t and a small variance r.

Variance reduction techniques
Reducing the variance r is important for studying the Monte Carlo method. In JMCT, the main techniques include weight window, importance, mesh window, and source biasing. Biasing parameters may be produced based on adjoint calculations. At present, JMCT has a function of adjoint calculation. Other variance reduction techniques include geometry splitting and Russian roulette, energy splitting/roulette, implicit capture, forced collisions, source variable biasing, correlated sampling, exponential

Tallies
All tallies can be expressed as where P ¼ ðr; E; X; tÞ, /ðPÞ is a flux, gðPÞ is a responding function. /ðPÞ satisfies the Boltzmann differential-integral equation as follows [20]: In the JMCT, the standard tallies include k-effective, current, surface flux, volume flux, point flux, deposition energy, detector response, and various reactive rates. In addition, the JMCT also provides a standard port of tally subroutine for user.

Algorithms
In recent years, a series of special algorithms for different problems have been developed for JMCT. Here, only a few key algorithms are introduced.

Domain decomposition and parallelization [21]
Currently, most Monte Carlo particle transport programs do not have a DD function which includes the MCNP program. The Mercury Monte Carlo particle transport program was the first program with a DD function [22]. If the case happens when one zone is divided into two or more zones (see Fig. 6a), there are some differences between a DD case and the no DD case. For JMCT, another DD algorithm has been developed, in which the nature of the geometry interface is determined by random sampling (see Fig. 6b). Due to the topology structure maintaining, the same results are obtained in the DD case and no DD case. Of cause, it is necessary for designing a clear chain of random numbers (see Fig. 7).
JMCT supports two-level parallelization of MPI and OpenMP. Figure 8a shows the communication between different domains in the DD case. Figure 8b shows domain replication, domain decomposition, and parallelization. Dr. Gang Li developed several algorithms to solve data exchange, dynamic load balance, and variance lower questions [21]. Figure 8c shows a comparison of the parallel efficiency in the DD and no DD case. A high parallel efficiency was still obtained even in the DD case.

Uniform tally density (UTD)
Monte Carlo methods are widely used for criticality calculations. For this purpose, the power iteration method is one of the most important techniques for the Monte Carlo program. To obtain an accurate multiplication factor and global tallies, inactive iteration cycles must be performed until the source distribution converges. Generally, tallying should only be invoked in the subsequent active cycles after completion of the inactive cycles. If the directed simulation uses, the errors in some margin assemblies may exceed the convergence rule of 95%. Therefore, the MC21 Monte Carlo program proposes a uniform fission site (UFS) algorithm for computing the expected number of fission secondary neutrons. UFS algorithm significantly improves the convergence of margin assemblies [23]. For JMCT, an improved algorithm based on the UFS algorithm has been proposed. This algorithm is called the uniform tally density (UTD) algorithm, and the modified formula is as follows: where m is the expected number of fission neutrons.k eff is the multiplication factor estimated in the last cycle, v k is a fraction of V occupied by zone k, in which the collision occurs, V is the volume of the problem domain that comprises all fissionable materials, s k is a fraction of the fission source contained in zone k. For same problem, a large FOM value is obtained based on UTD [24].

Doppler broadening on-the-fly (OTF)
The on-the-fly (OTF) Doppler method was used to calculate the temperature-dependent cross section of any nuclide at any temperature in the range of 300-3000 K based on a cross section library of 300 K. For the OTF method, first, a series of temperature-dependent cross section is produced by NJOY [11]. Second, a uniform energy grid was evaluated using temperature-dependent cross section. Third, a polynomial was used to fit the temperature-dependent cross section of each energy grid. The coefficient of the polynomial was obtained using a single-value decomposition algorithm. The basic formula of OTF cross-section is as follows: where t; a; f express the total, absorption, and fission. Finally, the coefficients of the polynomial in all energy grids and the energy grids themselves were written in a text file for interpolation [25]. The algorithms are based on the neutron balance equation. When considering significant reactions, such as (n, a), (n, f), (n, c), and (n, 2n), the following equation is obtained: where r boron is a microscopic cross section of soluble boron, R are the related macroscopic cross sections, F presents the neutron produced by fission. Assuming that the ratio of soluble boron density to water density is the same in the entire reactor, the boron concentration can be iterated together with k eff during the cycles. The new algorithm is only one cycle iteration for convergence and the traditional method usually needs four to five time iterations [26].

Coupled MC and S N
It is well known that consistent adjoint driven importance sampling (CADIS, or FW-CADIS) is an ideal method for solving the deep penetration problem. The importance values of Monte Carlo transport are produced by S N code. However, it needs two modeling for S N and MC calculation, respectively. It is good luck that JMCT (MC) and JSNT (S N ) use same modeling, it makes the CADIS method become easy and highly efficient. The mesh weight windows, importance values, and source biasing coefficients were produced by solving the adjoint equation based on the discrete ordinate S N program JSNT [12]. Figure 9 shows the coupled flow of MC and S N . Good efficiency has been achieved for some deep-penetration shielding problems [27,28].

Multi-physics coupled calculation
Tightly coupled neutron transport (JMCT), depletion (JBURN), and thermal hydraulics (JTH) have been developed for numerical reactor simulations. Recently, the JMCT has added the analysis of fuel properties and makes the simulation closer to the true. Figure 10 shows the iteration flow.

Collision mechanism based on material [29]
For most Monte Carlo multigroup neutron transport calculations, the collision mechanism is based on nuclides. A new collision mechanism has been developed based on the material. This mechanism is suitable for the case of a large amount of fission productions. The simulation can save a significant amount of memory and computational time. The basic principle is that suppose the material k being consists of m(k) types of nuclides. The combined cross sections of material k are defined as follows: where n i (i = 1,2, …, m(k)) is the percentage of each nuclide, and g, t, a, f express the energy group, total, absorption, and fission, respectively. The scattering cross section is defined as follows: If r g f 6 ¼ 0, then the number of fission neutrons of material k is defined as follows: The fission spectrum of material k is defined as follows: where v g;i is the fission spectrum of i fission nuclide; g = 1, 2, …, G, where g = 1 corresponds to the highest energy group and g = G corresponds to the lowest energy group. The scattering transfer matrix of material k is defined as follows:  If upper scattering is not considered, that is, nup = 0, then the scattering transfer matrix becomes a triangular matrix.

Treatment of multigroup angular distribution
The multigroup scattering cross section can be expressed as follows: where l ¼ X 0 Á X;X 0 , and X are the incident and emission directions, respectively, and f ðlÞ is an angle distribution function that satisfies normalization, that is, In general,f ðlÞ is expanded in a Legendre series as following where Legendre coefficients are The approximate distribution is obtained by making Lorder cutoff: Due to this cutoff, f L ðlÞ may appear negative in [-1, 1]. To avoid this case, the generalized Gaussian quadrature is used [30] and an approximate discrete distribution is obtained as follows: where n ¼ Lþ1 is the probability of the sample l i ( P n i¼1 p i ¼ 1). Q i ðlÞ f g n i¼1 is an orthogonal polynomial system with respect to f L ðlÞ, l i f g n i¼1 is the root of Q n ðlÞ and satisfies where d ij is the Kronecker delta function and N i is a normalization constant [29]. f Ã ðlÞ replaces f L ðlÞ and l is produced through sampling.

Random geometry modeling
The JMCT is used to simulate a high-temperature reactor (HTR). The core of the HTR consists of many fuel pebbles, with each fuel pebble consisting of many individual pellets [31]. According to the traditional method, these pellets are uniformly mixed into a single material. Today, precision physical research is needed to study the effect of pellet numbers and different distributions on pebble criticality. Therefore, the geometry must describe according to the real distribution of pellet. So, a sample method has been developed based on the real distribution of pellet. Figure 11 shows the uniform pellet distribution in the pebble.

Numerical tests 3.1 BEAVRS model
The BEAVRS model is thought to be a significant challenge for computers and software. The model was released by the MIT Computational Reactor Physics Group on July 7, 2013 (www.crpg.mit.edu) [4]. It includes detailed specifications of the operating 4-loop Westinghouse PWR (3411 MW), two cycles of measured data, hot  zero power (HZP) and hot full power (HFP) data, fuel loads by assembling as built, and three enriched fuels (1.6%, 2.4%, and 3.1%). Two cycles of measured data were used to validate the high-fidelity core analysis program. The simulation results of MC21 and OpenMC in the HZP condition were presented in PHYSOR2014 [32]. The basic parameters are as follows: The core contained eight types of rods, as shown in Fig. 12a, and nine types of assemblies, as shown in Fig. 12b. Each assembly had six segment grids in the axial   Temperature coefficient under 560°F, T 2 = 570°F, T 1 = 550°F direction, as shown in Fig. 12c. Because the model memory is too large, 2 9 2 9 2 DD are used, as shown in Fig. 12d. The simulation tracked 4 million neutrons per cycle, 1000 cycles in total, and discarded 400 cycles. The simulation consumed 200 processors and 5.3 h in the Tianhe-II computer. Figure 13 shows a comparison of the radial integral assembly power between the MC21 and JMCT. The maximal difference is 3.17% in the E2 assembly, and the minimum difference is 0 in the E6 assembly (MC21 appears in the same assemblies). Figure 14 shows a comparison of the detector tallies and measured data in the instrument tubes. The maximal deviation of the power assemble was -14.77% in the B13 assembly (MC21 was -17.05%). The minimum deviation of the power assemble was -5.648% in the L15 assembly. Therefore, comparisons of the axial detector power signals were performed for the B13 and L15 assemblies. The JMCT results were located between MC21 and the measured values (see Fig. 15). Figure 16 shows a comparison of the radial pin power in the axial plane between the MC21 and JMCT. The standard deviation of 99% of the fuel pin power was less than 1% (see Fig. 16c). Table 1 shows the k eff comparison of JMCT, OpenMC, and MC21 at the different locations of the control rods and boron concentrations. Table 2 lists the reactivity worth of the control rods at 560°F. Table 3 presents a comparison of the temperature coefficients. The integrated axial power distribution of the experiment had a slight asymmetry. We guess that this was caused by the instrument tubes. Generally, the comparison shows good agreement among the JMCT, MC21, and OpenMC models [33].

HZP condition
As a marked achievement of the CASL plan [34], VERA core physics benchmark progression problem specifications were published on March 29, 2013 [5]. Figure 17 shows the VERA modeling by the JLAMT. Table 4 compares the k eff and control rod worth between the JMCT and KENO-VI [10]. Figure 18 shows a comparison of the radial integral power distributions of the JMCT and KENO-VI. The maximum difference was 0.85%. Good agreement was achieved between the JMCT and the KENO-VI.   Figure 19 shows a comparison of the radial integral assembly power distribution between the JMCT and MC21. The maximum difference was 2.27%. Figure 20 shows the CAD modeling using JLAMT. The start of the physical experiment of the first core in units A(b) and B(b) was simulated by JMCT, NECP-X [35] and SCAP [36]. Table 5 lists the deviations of the simulation results relative to the experimental results for unit A(b).    Figure 21 shows the deviation of the control rod worth of units A(b) and B(b).

Reactor core
For CAP1400, the initial physical parameters of the first core were calculated in a similar fashion to that of AP1000. The JMCT simulation results were selected as the reference values. Table 6 lists the deviations of NECP-X and SCAP relative to JMCT [37].

Shielding
The radiation shielding of the CAP1400 includes eight operating cases (i.e., eight models). Here, only the simulation results for a typical case are presented. First, a complex source subroutine is designed according to the distribution of position, energy, and direction. The fixed source and mesh are tally adopted, where the importance values are produced by the S N program JSNT. A total of 26,000 million neutron histories were simulated using 260,000 processors. The total CPU time was 17.7 min. The calculation was completed using a Tianhe-II computer. Figure 22 shows the neutron dose distribution and secondary photon dose distribution of the full factory.  Control bank worth of MA/10 -5 205.6 -5.2 9 10 -5 -7.6 9 10 -5 Control bank worth of MB/10 -5 172.8 4.0 9 10 -5 2.8 9 10 -5 Control bank worth of MV/10 -5 210.9 0.4 9 10 -5 -2.0 9 10 -5 Control bank worth of MD/10 -5 187.5 -4.8 9 10 -5 -6.9 9 10 -5 Control bank worth of M1/10 -5 498.

Summary
A high-fidelity general-purpose 3-D Monte Carlo particle transport program, JMCT3.0, was developed for integrated simulation of nuclear systems, where the geometry zones break through ten million, the depletion regions break through millions, and parallel processors break through one hundred thousand. Advanced computer technologies, automatic CAD modeling, and visualization make the program with a high performance. The JMCT has carried out a number of simulations in reactor full-core, radiation shielding [38], nuclear detection [39] and nuclear medicine [40]. A strong simulation capability was shown. In recent years, JMCT has completed the simulations of many large and complex nuclear devices, such as CFR600, ACP100, Hualong, Qingshan-I, VVER, CFETR, and HTR-10 etc (see Fig. 23). The necessary technology supports ware provided for the optimization of nuclear device designs. Future efforts will be forced to enhance computing efficiency. Depletion, complicating the uncertainty quantification and propagation of errors is an essential area to consider in the future. Furthermore, new algorithms need to be developed to reduce computing fees. JMCT is still evolving toward this goal, and all algorithms are being actively developed at present.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.