Final Report Summary - HPMC (High performance Monte Carlo reactor core analysis)
Executive Summary:
The general objective of the HPMC Collaborative Project is to develop and demonstrate the application of full-core Monte Carlo calculation for time-dependent safety analysis with thermal-hydraulic feedback and burnup using high-performance computing. These goals fit into the Strategic Research Agenda of the Sustainable Nuclear Energy Technology Platform (SNETP).
The project is organised into 4 research work packages as follows
-WP1 aims at optimum Monte Carlo - thermal-hydraulics coupling; its objective is to realise efficient coupling of the Monte Carlo codes SERPENT and MCNP with the thermal-hydraulic subchannel codes SUBCHANFLOW and FLICA4, suitable for full-core application
-WP2 aims at optimum MC-burnup integration; its objective is to realise an efficient integration of burnup calculations in the Monte Carlo codes SERPENT and MCNP, suitable for full-core application
-WP3 aims at time dependence in Monte Carlo calculations; its objective is to realise an efficient time-dependence option in the Monte Carlo codes SERPENT and MCNP for safety analysis, applicable for full-core calculations
-WP4 aims at integration of high-performance parallel Monte Carlo; its objective is the demonstration of full-core Monte Carlo based dynamic safety analysis with ultimate efficiency in parallelisation
Two other work packages are defined
-WP5 aims at dissemination and training; its objective is communication of the project to potentially interested parties, including presentation of two training seminars for potential users of the developed code system with hands-on training
-WP6 covers the management of the project
Advanced and innovatives methods and code extensions were developed to be able to perform high fidelity core simulations based on Monte Carlo codes coupled with thermal hydraulics. For the demonstration of the new capabilities of the developed methods, a benchmark for a mini reactor core was defined that permitted to do code-to-code comparison. It demonstrated the efficienty and accuracy of the coupled Monte Carlo and thermal-hydraulic calculations. Integrated burnup calculations with a new and stable iteration scheme showed consistent results even for relatively large time steps. The time-dependent reactor behaviour with delayed neutron precursors after a control rod movement could be demonstrated for a mini fuel assembly. The efficiency of parallel execution of a full-size reactor criticality calculation on a massively parallel supercomputer could be increased to almost 100 %. Two successful training seminars were given to interested people from industry and research institutes.It is concluded that the project forms a solid basis for the next step forward to practical safety analysis validation of nuclear reactors based on the Monte Carlo method.
Project Context and Objectives:
The general objective of the HPMC Collaborative Project is to develop and demonstrate the application of full-core Monte Carlo calculation for time-dependent safety analysis with thermal-hydraulic feedback and burnup using high-performance computing. These goals fit into the Strategic Research Agenda of the Sustainable Nuclear Energy Technology Platform (SNETP).
The project is organised into 4 research work packages as follows
- WP1 aims at optimum Monte Carlo - thermal-hydraulics coupling; its objective is to realise efficient coupling of the Monte Carlo codes SERPENT and MCNP with the thermal-hydraulic subchannel codes SUBCHANFLOW and FLICA4, suitable for full-core application
- WP2 aims at optimum MC-burnup integration; its objective is to realise an efficient integration of burnup calculations in the Monte Carlo codes SERPENT and MCNP, suitable for full-core application
- WP3 aims at time dependence in Monte Carlo calculations; its objective is to realise an efficient time-dependence option in the Monte Carlo codes SERPENT and MCNP for safety analysis, applicable for full-core calculations
- WP4 aims at integration of high-performance parallel Monte Carlo; its objective is the demonstration of full-core Monte Carlo based dynamic safety analysis with ultimate efficiency in parallelisation
Two other work packages are defined
- WP5 aims at dissemination and training; its objective is communication of the project to potentially interested parties, including presentation of two training seminars for potential users of the developed code system with hands-on training
- WP6 covers the management of the project
Four organisations participate in the HPMC project:
o Karlsruhe Institute of Technology (KIT), Germany, also project coordinator
o Delft Nuclear Consultancy (DNC), The Netherlands, also project secretary
o VTT Technical Research Centre of Finland (VTT), Finland
o Royal Institute of Technology (KTH), Sweden
WP1 was successful in developing coupling methods for Monte Carlo (MC) and thermal-hydraulics (TH) calculations for large-scale applications. Besides the existing external coupling scheme of the MC and TH code, an internal coupling scheme was developed in which the TH code is modified to become a subroutine of the MC code. This provides higher efficiency as no communication via input and output files between the MC and TH codes is necessary. A Python based framework PIRS has been developed to define in a simple way complex geometries and to prepare the input file for either the MCNP MC code or the SCF TH code. This framework can be used to generate input files for the MC code for a complete fuel assembly and a full reactor core.
A dynamic assignment of the cross section data for the required temperature in each subzone of the reactor system was developed and implemented, which no longer requires that each subzone that can have a different temperature, has to be defined separately in the geometry input to the MC code. This facilitates the MC-TH coupling in a practical way for large fuel assemblies and full-size reactor cores.
A theory was developed and implemented to optimise the number of neutron histories to be applied in successive iteration steps of the MC-TH coupling. The number of neutron histories increases almost linearly with the MC-TH iterations. At the same time the relaxation parameter for combining the power distribution obtained from successive MC iterations as input to the TH calculation is optimised.
Other improvements are the introduction of on-the-fly temperature treatment in the SERPENT code instead of stochastic mixing of nuclide cross section data at two temperatures and the speedup of fission source convergence in the Monte Carlo calculation with the aid of an approximate deterministic calculation of the fission source distribution, taking into account the proper temperature distribution. Also new power tally estimators based on a collision estimator were developed for the detailed power distribution in a core which are much more efficient than existing estimators and hardly take more additional computing time.
To demonstrate the current capabilities for MC-TH coupled calculations a benchmark was defined for a mini core reactor with 3x3 fuel assemblies containing both UOX and MOX PWR fuel assemblies. Various codes were used showing reasonable agreement in the detailed power and temperature distributions over the core.
In WP2 current algorithms for coupling MC calculations with burnup calculations are investigated with respect to their stability. It was found that most algorithms are only conditionally stable and often instable at the relatively large time steps applied in burnup calculations.
A new algorithm, the Stochastic Implicit Euler method, was developed and successfully tested. The algorithm was implemented in the SERPENT code. Extensive calculations show stable results even at large time steps.
In WP3 time dependence was introduced in the MCNP5 Monte Carlo code, including the generation and decay of delayed neutron precursors. To reduce the statistical error in the generated reactor power in successive time intervals, a method of forced decay of precursors in each time interval was implemented. Also other variance reduction methods like the branchless collision were introduced. Nonetheless the statistical variation in the length of the fission chains and therefore the generated power per starting neutron in a time interval remains large in heterogeneous systems.
In order to simulate reactor transients, an efficient method was implemented to model the movement of control rods at any speed during a calculation.
Thermal-hydraulic feedback was implemented for the time-dependent version of MCNP to create a fully dynamic Monte Carlo code. To let the time-dependent TH calculation take into account the heating history further extensions of both codes were necessary. For an external MC-TH coupling scheme a method of exchanging signals between the independently started MC and TH codes was implemented to arrange exchange of information at the proper time instances. With further modifications an internal coupling was also realised.
The current capabilities for time-dependent calculations are demonstrated with a power transient due to control rod movement in a mini fuel assembly. Further developments are needed to accommodate safety calculations for large systems.
In WP4 various ways for parallel execution of a Monte Carlo calculation using the MPI and OpenMP application programming interfaces were investigated and their efficiency measured in terms of the speedup factor. For application on large computer clusters with different computer nodes and multiple processors per node, the optimum combination of MPI and OpenMP was determined. Application of OpenMP was introduced in the SERPENT2 code. The MCNP code was modified to use all available processor cores for neutron history simulation.
To increase the efficiency of the Monte Carlo calculation further, the amount of data communicated between computer processor cores after each cycle in a criticality calculation could be minimised. Nonetheless, the speedup of the parallel calculation using many processors was still not satisfactory. Therefore, a rigorous method was developed for application to massively parallel computer systems which avoids almost all communication between different tasks running on different computer nodes.
Access to a supercomputer via the European PRACE organisation was obtained for two successive stages of the project. It could be shown that parallel speedup for a criticality calculation of a full-size core was excellent when using the wall-clock time as a limit for criticality calculations instead of the number of cycles.
As part of the activities in WP5 a public website is setup with information about the project to bring the project and it goals to the attention of possibly interested parties. A separate website for internal use to the project was established to exchange information and documents between participants. A poster about the project was presented at the general assembly of the Sustainable Nuclear Energy Technology Platform (SNETP) in Warsaw and the FISA conference in Vilnius. Presentations about the project were given on several occasions.
With the EU project NURISP and its successor NURESAFE agreements were made about the use of benchmark problem definitions and results and about common problems to be investigated at both sides.
An extended summary of the project was prepared for publication by the European Commission. For dissemination of the results of the project a Communication Action Plan was written.
Two successful training seminars were organised in Germany and Sweden, respectively, which were attended by in total 22 external people from 11 different companies and research institutes.
Scientific results of research performed for the project were published in several articles in international scientific journals and in papers at international conferences.
WP6 covers the management of the project. The Executive Board gathered every half a year. In the board meetings the progress of work was discussed by means of presentations by all participants involved in the project. The work package leaders presented the status and outlook of the work in their work package.
All deliverables scheduled for the project were published, some of which with some delays.
It was concluded at the last Executive Board meeting that basically all tasks were executed and that important new and advanced developments were realised and successful results demonstrated for most of the tasks. Therefore the HPMC project forms a solid basis for the next step forward to practical safety analysis of nuclear reactors based on the Monte Carlo method.
Project Results:
The research and development work of the project is organised in the first four work packages. Work packages 1 to 3 cover the development of the essential items for realistic reactor safety analysis with Monte Carlo calculations, namely coupling of Monte Carlo with thermal-hydraulics (WP1), integration of burnup calculations (WP2) and time-dependent Monte Carlo methods from seconds to minutes scale (WP3).
Work package 4 is directed to optimum integration of the techniques developed in the first three work packages, to maximum efficiency in parallel execution of the Monte Carlo codes and demonstration of the developed tools for a full-size reactor core.
Work package 5 is directed to dissemination of the results of the project and to training sessions for potential users.
Work package 6 concerns management of the project and is discussed separately in Sect. 5.
3.1 Work Package 1: Optimum MC-TH coupling
In the HPMC project the main Monte Carlo codes used are MCNP and SERPENT. The main thermal-hydraulics code is SUBCHANFLOW (SCF), developed at KIT. Incidentally the FLICA4 code from CEA, Saclay, France, is used. Where necessary these codes were obtained and installed. For new codes at a partner’s institution, an extensive familiarisation process is passed through as not only the structure of the input and output files must be known, but also the structure of the code itself. Especially the SERPENT and SCF codes are under continuous development, which requires regular installation of updates.
3.1.1 General considerations and developments on MC-TH coupling
Coupled Monte Carlo (MC) and thermal-hydraulics (TH) calculations require a series of iterations as the TH calculation needs the detailed power density distribution from the MC code, while the MC code needs the temperature distribution from the TH code. Neither distributions are known on beforehand, so one has to start with a guessed power distribution for the TH code, which calculates the temperature (and coolant density) distribution. This distribution must be input to the MC calculation, from which an improved power distribution is obtained. Then successive iterations will, in general, lead to a converged solution for both the power distribution and the temperature distribution.
For the coupling of Monte Carlo calculations with thermal-hydraulics calculations two schemes are followed, external coupling and internal coupling. External coupling uses both the Monte Carlo and thermal-hydraulics code as such and provide communication between them by retrieving the necessary data from output files produced by the MC code or the TH code and integrating these data in the input file of the TH code or MC code, respectively, for the next iteration. This way of coupling has the advantage that the original codes can be used and in fact only an executable version of the codes is needed. With proper conversion programs for converting the output of one code to input to the other one, various MC and TH codes can be coupled.
With internal coupling the TH code is made a subprogram of the MC code and is called as a subroutine. It needs modification of both the MC and TH code, but provides a simpler and faster way of data exchange between the MC and TH parts.
For a coupled MC-TH calculation the MC code needs not only to deal with the detailed geometry description of a reactor with many fuel assemblies and many fuel pins in a fuel assembly, but also with the different cross section data needed for each axial subzone in each fuel pin because of the different temperatures for fuel, cladding and coolant all over the reactor. To facilitate the generation of an input file for the MC codes which includes not only the geometry description, but also the assignment of cross section data at different temperatures for each axial subzone, a sophisticated system called PIRS, based on the Python language is developed.
Coupling schemes are developed and implemented at different participating institutes, which gives the possibility to evaluate different implementations and select the most promising developments for the future. External coupling has been realised for MCNP-SCF, MCNP-FLICA, Serpent-SCF and as a side product TRIPOLI-FLICA. Internal coupling has been realised for MCNP-SCF, Serpent-SCF and Serpent-COSY (a light-weight TH code developed at VTT).
Deliverable D1.4 Parts I and II, gives a description of the MC-TH coupling schemes as developed by the various partners in the first stage of the project. Afterwards further developments took place, especially to improve the internal coupling schemes on different items.
As the cross section data for neutron interaction with material are temperature dependent, each subzone in the Monte Carlo calculation must be treated at a specific temperature according to the outcome of the TH calculation. As all temperatures in a wide range will occur in the reactor, it is not realistically possible to supply neutron cross section data at all these temperatures for all relevant nuclides. Hence the technique of stochastic mixing is further developed to solve this problem. Here for each nuclide present in a medium two references to different temperature cross section data are defined with such atom fractions that on the averaged the desired temperature data is selected in simulating histories with selecting a reacting nuclide proportional to its atom fraction. With a variation this technique can also be applied for treating S(α,β) scattering for thermal scattering to bound atoms in a molecule like H in H2O. However, also a technique for direct interpolation into the S(α,β) kernel has been developed.
The stochastic mixing technique still requires a separate composition definition for each axial subzone of each pin cell, as all temperatures may differ. For large systems like a fuel assembly and certainly a full reactor core, this requires a huge number of separate subzones that cannot be handled in practice. One solution is the development of the PIRS system (Sect. 3.1.2). Another way to handle this problem is the dynamical assignment of nuclide cross section data at the applicable temperature when a neutron enters a certain zone during its history simulation. This keeps the geometry description in the input to the MC code relatively simple and limited and provides a major step forward in MC-TH coupling for a full reactor core. The methods were fully developed and demonstrated to give good results.
Another way of treating the temperature dependence with respect to cross section data is “on-the-fly” temperature treatment of neutron cross sections, which is developed for the SERPENT code. This implies that each time a neutron at a given energy needs the cross section at the temperature of the material the neutron is in, the code samples a relative energy taking into account the thermal motion of the nucleus at the specified temperature. This method takes some additional computation time, but avoids several problems otherwise met for temperature treatment. At the moment it still has some limitations.
An important issue in optimisation of MC-TH iterations is the control over the number of neutron histories to simulate in the successive Monte Carlo runs. For the first few iterations the power profile and thus the temperature profile may be far from converged and using too many neutron histories will be a waste of computer resources. In WP2 a procedure based on a physical theory is derived for optimised control of the number of histories in successive iterations, which is not only applicable to burnup calculations, but also to MC-TH coupled iterations. The number of histories turns out to increase almost linearly with the iteration number until convergence. This theory also provides an optimised relaxation parameter for the power profile to include the results for the power profile from previous iterations. Tests show that these procedures lead to a converged result.
Besides the above-mentioned subjects, various other items are studied and methods are developed. One issue is the variance reduction in power estimation in subzones all over the reactor using the Uniform Fission Sites (UFS) method.
Another issue is the convergence of a Monte Carlo calculation with respect to the fission source distribution. To speed up the convergence it is important to start the Monte Carlo calculation with a fission source distribution as close as possible to the fundamental mode distribution, taking into account the geometry of the reactor system and the temperature distribution. To this end coupling with a deterministic nodal method is developed, taking into account the temperature distribution obtained from the TH code. The deterministic code calculates the fission source distribution with as few approximations to the geometry as possible, which is then used as the starting source for the first Monte Carlo source convergence cycle. This method shows promising results.
3.1.2 PIRS - Python Interfaces for Reactor Simulations
As the geometry description of a reactor system in a Monte Carlo input file can be rather complicated, it is important to have a tool that facilitates the generation of such an input file. Therefore, an extensive set of Python based modules with definition of various classes and class methods is developed, which allows to define geometry and code-specific data for MCNP and SCF in terms of the Python programming language. It also gives the possibility to start the codes, so iterative MC-TH calculations can be organised. Moreover, It can provide plots of the temperature and power distributions as calculated by the TH and MC codes.
The Python based framework is operational for a single pin cell geometry, a realistic fuel assembly with a large number of fuel pins and for a complete core. It can generate inputs for MCNP and Serpent, as well as for the coupling of these codes with the SCF thermal-hydraulics code.
However, for larger systems, for instance the mini core benchmark as discussed in Sect. 3.1.5 the generated input file for the MC code is extremely large as the MC input file is based on defining each subzone separately because of the varying temperatures. It takes MCNP a very long time to process the input file at initialisation of a run.
3.1.3 External MC-TH coupling
External coupling has the advantage that it can use, in principle, both the MC code and the TH code in standard form and basically only an executable version of these codes is needed. A script is needed to call successively the TH and MC code. Coupling between the codes is by output and input files. As the input for one type of calculation (either MC or TH) is totally different from the data provided in the output of the other code (TH or MC, respectively), separate conversion programs are needed to read the output file(s) from for instance the TH code and prepare from it a new input file for the MC code, using a standard input file for the MC code in which relevant data must be updated from the results of the TH code. The same must be done for preparing a new input file for the TH code.
Although the conversion programs for reworking the output from one code to a new input for the other are dependent on both the MC and TH code used, it was possible to couple in this way different MC codes (MCNP, Serpent, TRIPOLI) to different TH codes (SCF, FLICA). However, the communication via input and output files and the conversion of these files is rather inefficient and limits parallel execution to the separate runs of the MC code. For larger problems the output files (power distribution and temperature distribution, respectively) can be very large for a large number of subzones. Hence, this type of calculation is useful for comparison of codes and code combinations for smaller problems, but is not suitable for calculating large systems.
3.1.4 Advanced internal MC-TH coupling
As can be concluded from the previous sections, internal coupling of MC and TH calculations is the most efficient way and the only one to deal with large systems like a full-size reactor core. Hence, in the project, after exploitation and testing of various ways for coupling, most attention was paid to the development of efficient methods of internal MC-TH coupling.
As internal coupling requires modifications to the TH code (as well as to the MC code), the KIT code SubChanFlow (SCF) was the only TH code for which the source code was available to participants and could be adapted to become a subroutine that could be called by the MC code. Also because of the structure of the SCF code this was not a problem and various arrays with temperatures for the fuel, cladding and coolant and with the coolant density for all subzones can be made accessible in the MC code. The MC code also needs modifications to call the subroutine which is the start of the TH calculation.
For the MC-TH iterations several solutions are possible, but most efficient is to modify the MC code to include internally a loop which successively calls the TH subroutine and the subroutine starting the Monte Carlo neutron history simulation. Then the geometry, material composition and cross section data needs to be read only once. For the TH code also modifications are necessary to read the input only once and to skip certain parts which need not to be recalculated each iteration and may only reinitialise internal values that can be useful to keep for continuation the internal TH iterations when a new MC-TH iteration must be performed. The power distribution per subzone, which is needed in the TH calculation, can be put in the corresponding array in the TH code from the power tally results in the MC code after some reworking and renormalisation.
After execution of the TH calculation the optimised scheme for increasing the number of neutron histories per MC-TH iteration as mentioned in Sect. 3.1.1 is implemented in the MC code in the loop for MC-TH calculations as well as a procedure to test the convergence of the MC-TH iteration. To accommodate the different temperatures in all subzones in the reactor system, stochastic mixing is applied in the standard input to the MC code. For dynamic temperature linking to cross section data the subroutine providing the cross section data at the current neutron energy when entering (or starting in) a new medium is modified and the atom fractions of the nuclides for stochastic mixing are temporarily adapted to the temperature for that medium.
Other modifications to improve the calculation, for instance for speeding up the convergence of the fission source distribution, for improved estimators for the power distribution and for global variance reduction for estimating the detailed power distribution are discussed in Sect. 3.4.
Another important improvement concerns the introduction of a new estimator for the detailed power distribution based on a collision estimator. This reduces the additional time for tallying considerably. An advanced version did not even take noticeably more computing time.
Apart from the demonstration shown in the next section for a mini core reactor, successful coupled MC-TH calculations could be performed for a quarter core, which shows the capacities of the developed methods.
3.1.5 Demonstration and Conclusions
To demonstrate the capabilities of the developed MC-TH code coupling a computational benchmark was defined to obtain and compare numerical results. The benchmark definition and the results are detailed in deliverable D1.9. The benchmark problem consists of a mini-core of 3x3 PWR fuel assemblies, from which 6 use UOX fuel and 3 use MOX fuel with 3 different Pu concentrations. Each fuel assembly has 17x17 fuel positions from which 25 are occupied by guide tubes (without control rods) for the UOX assemblies and with WABA (Wet Annular Burnable Absorber) rods for the MOX assemblies. Hence, the geometry and composition is quite complicated. Although the project definition required the demonstration for a single fuel assembly, it was felt that the project was up to the treatment of a larger system.
Four independent solutions were obtained: two with internal coupling of MCNP-SCF (but different implementations), one with external coupling of MCNP-SCF using the PIRS system for generating the input files and one with internal coupling of Serpent-SCF.
Although the results show some differences, as discussed in the deliverable report, and further detailed analysis is desired, it can be concluded that various developments have led to coupled code systems capable of calculating a (mini) core system with realistic geometry and detailed fuel and coolant temperature and coolant density distribution throughout the core.
Fig. 1 shows from top to bottom the radial distribution through the mini core of the fuel temperature, the coolant temperature and the coolant density. The figures at the left are at the bottom of the fuel, the figures at the right are near the top of the fuel. Red colours indicate higher values, blue colours indicate lower values. The fuel assemblies with mainly light and dark blue are the MOX fuel assemblies.
Fig. 1. Radial distribution of fuel temperature, coolant temperature and coolant density at bottom, middle and top of fuel
3.2 Work Package 2: Optimum MC-burnup integration
Although there exists various examples of Monte Carlo codes that can be coupled to burnup calculations using different algorithms, the stability of the algorithms are not sufficiently investigated. Therefore, a simplified model was devised to test various coupling schemes. The model is based on the fission matrix concept and imitates reasonably well the features of standard Monte Carlo calculations, but at much lower computational cost. Moreover, Xenon feedback is implemented, which can play an important role in instabilities. As this model allows to run the simulation also in deterministic mode, it can determine exactly the statistical errors due to the Monte Carlo simulation.
Extensive calculations with this model showed that most of the applied algorithms for burnup feedback are unstable when relatively large time steps are used, as is often the case in burnup calculations.
3.2.1 A Stable MC-burnup scheme
On basis of these results a new coupling scheme was designed and successfully tested. The fully implicit scheme is robust and efficient. It is called the Stochastic Implicit Euler (SIE) method. The results of this research are detailed in deliverable D2.3 and published in two journal articles. The results are both used for burnup calculations with the SERPENT code and for coupled Monte Carlo – thermal-hydraulics calculations (see Sect. 3.1).
Fig. 2 shows the axial offset and spatial distribution of 135Xe as a function of burnup. The top rows are without equilibrium 135Xe, the bottom rows with equilibrium 135Xe calculation. For the figures from left to right the step size in the calculation is 10, 5, 2.5 and 1MWd/kgU.
Fig. 2. Results for axial offset and spatial distribution of 135Xe as a function of burnup indicated by the number of depletion time steps
3.2.2 An optimised scheme for MC-TH iterations
From the theory of stochastic approximation an optimised scheme for MC-TH iterations is developed. It prescribes the increase in number of histories (or number of cycles) in Monte Carlo runs for successive MC-TH iterations. As the inaccuracy of the temperature distribution over the reactor will decrease with the MC-TH iterations, it pays off to use more cycles and neutron histories for better statistical accuracy in later MC-TH iterations. At the same time the optimised scheme determines the relaxation of the power distribution from successive MC-TH iterations with the relaxation factor in each iteration determined by the number of active histories relative to those in the first iteration. This scheme was successfully used in most of the MC-TH coupled calculations discussed in Sect. 3.1.
3.3 Work Package 3: MC time dependence
Although time dependence is almost automatically included in Monte Carlo codes, it is normally limited to the timescale of the neutron life time, which can vary from microseconds to about 1 ms, depending on the reactor type. For reactor transient analysis with Monte Carlo the generation and decay of delayed neutron precursors must also be modeled. In long-time time-dependent Monte Carlo calculations neutron chains of successive fission neutrons, initiated by the decay of a precursor nucleus, play a dominant role. Although such chains can contain hundreds of generations of fission neutrons, their life time is normally short compared to the physical phenomena generally occurring in transient situations.
A study showed that the statistics in neutron chain length are large and therefore the generated power by a chain of fission neutrons, too. This requires special techniques when implementing time dependence in a Monte Carlo code.
In the SERPENT2 code the time-dependent transport model was implemented, but without delayed neutron precursors yet. Results for short-time time dependence are published in a conference paper. All other developments for time dependence were implemented in MCNP which is then called dynMCNP.
3.3.1 Implementation of time dependence with delayed neutron precursors
Generation and decay of delayed neutron precursors was implemented in MCNP. Like fission neutrons the precursors are stored in a bank in computer memory, but with a flag to distinguish them from prompt neutrons. To describe the produced power in a reactor system during a transient, time intervals are introduced. Internally in dynMCNP time intervals are treated like successive criticality cycles, that is, all source particles for a certain time interval are successively retrieved from the particle bank and their history simulated. In a collision the result can be a capture event, which ends the history of that neutron, a scattering event with continuation of the neutron history or fission (in a fissionable medium) in which case the history of the neutron causing fission is ended and n fission neutrons are generated from which n-1 are stored in a temporary bank and the history of one fission neutron is continued. When this history is ended, possibly after more fission events, the last fission neutron is retrieved from the temporary bank. If a neutron reaches the end of the time interval, it is stored in the particle bank as prompt neutron for the next time interval. If in a fission event a new precursor is generated, it is also stored in the particle bank, marked as a precursor. When all particles for the current interval are processed execution for the next interval will be started.
To reduce the statistical uncertainty in the produced power per time interval, existing precursors are forced to decay in each time interval with a correction for their statistical weight. Another general option for variance reduction was implemented: branchless collision. In that case always one neutron will result from a collision. Therefore, capture is treated implicitly by modifying the weight of the particle. However, in our case the treatment of fission is also modified in such a way that only one new fission neutron is generated. This requires further adaption of the neutron weight, which can possibly increase in fissionable media. A neutron history can only be terminated by escape from the system or by Russian roulette if its weight becomes too low.
In order to simulate transients due to the movement of control rods, an efficient method was implemented to determine for each neutron where the tip of the control rod is located, depending on the time variable of the neutron. It is possible to move control rods together as a “bank”, as will be the case in actual reactors. In additional input to the MC code the movement of a control rod can be described. The control rod can move linearly in time during each time interval. As the introduction of temperature dependence with an externally given temperature distribution (for instance as a result of a thermal-hydraulics calculation) per subzone with dynamic temperature linking (see Sect. 3.1.1) requires the definition of a lattice in MCNP for axial division of subzones, this prohibits the definition of cell for a control rod from the control rod tip to the top of the reactor and of a cell below the control rod tip with coolant down to the bottom of the reactor. In order to avoid multiple cell definitions of a control rod over various axial subzones, which would make modelling a movable control rod very difficult, subroutines of MCNP were modified that bypass the axial lattice in case of specific cells limited by surfaces used for determining (possibly moving) control rod tips.
Efficient parallel execution of the Monte Carlo code for the simulation of neutron histories in a single time interval requires additional measures. In order to reproduce results of a calculation with a different number of computer processors, the precursors generated during the interval and the prompt neutrons reaching the end of the time interval must be stored in a bank separate from the precursors and prompt neutrons present at the start of the time interval and must be sent separately to the master processor in order to include them in the correct order in the bank for particles starting in the new time interval.
3.3.2 Coupling with a time-dependent thermal-hydraulics code
Coupling of a time-dependent simulation with thermal-hydraulic feedback was implemented for dynMCNP to create a fully dynamic Monte Carlo code. The external MC-TH coupling scheme as developed in WP1 had to be extended, as for a time-dependent calculation the TH calculation must take into account the heating history. Therefore, both the MC and the TH codes run continuously and a method was developed of sending each other signals to indicate when they have finished the computational work for a certain time interval. When the TH code has finished the processing for a time interval the MC code receives the signal, collects the results of the TH code for temperature and coolant density distributions and starts the simulation of neutron histories for that time interval using the dynamic assignment of nuclides at the applicable temperature as developed in WP1. When the MC code has completed the work for the time interval the TH codes receives a signal from the MC code, collects the estimated power distribution from the MC code and starts the calculations for the next time interval. When one code is actively doing calculations, the other code goes asleep, waiting for a signal to awake and take over the computation. The SUBCHANFLOW (SCF) code was modified to send and receive signals for the dynamic MC-TH coupling. The coupling method was demonstrated to work properly. However, due to still consisting problems with large statistical uncertainty in the generated power in a system, it was not possible to produce a convincing demonstration of a time-dependent case due to a moving control rod with thermal-hydraulic feedback. The deliverable D3.5 describes in detail the time-dependent external MC-TH coupling mechanism, but without a demonstration calculation of a time-dependent case with thermal-hydraulic feedback.
As internal coupling of Monte Carlo and thermal-hydraulics calculations will be more efficient for large reactor systems using many processor cores on a computer cluster or supercomputer, internal coupling of dynMCNP with SCF was also developed. As in internal MC-TH coupling it requires a loop in MCNP which successively calls the routine in MCNP for processing particles for a certain time interval and then calling SCF for execution of the time step in the time-dependent version of MCNP. This requires further modification of SCF to step out of SCF after each time interval and to continue its time-dependent calculation at the next entrance, without recalculating several data. Also dynMCNP requires further modification when re-entering the simulation routine for the next time interval to continue with all data gathered at the end of the previous interval.
3.3.3 Results and Conclusions
Results for a time-dependent case (without thermal-hydraulic feedback) were obtained for a small 3x3 fuel assembly in which the central rod is replaced by a movable control rod. The outer boundaries of the system are reflective. The transient is introduced by linear movement of the control rod tip from 0.2 s to a new position at 0.5 s, after which the position was fixed. The total reactivity effect was about 0.6 $. Intervals of 0.1 s were considered. The power over the assembly was recorded per time interval during about 1.5 s. Fig. 3 shows the resulting power transient. For comparison a point-kinetics calculation was made with the estimated total reactivity from a separate static MCNP calculation. For this rather small system there is a fair agreement between both calculations.
It can be concluded that the implementation of long-time time dependence with generation and decay of delayed neutron precursors is successful. However, there are several items that need further attention and improvement. Firstly, the statistics in the recorded total power are considerable. Therefore, the calculation required simulation of as many as 106 neutrons histories per time interval. Secondly, the efficiency of the parallel execution of the Monte Carlo runs was rather low, which restricts the possibilities to scale up to massively parallel execution on a supercomputer. This is very much needed as the calculation takes a long time due to the long neutron chains per time interval because of many successive fission events. Thirdly, the implementation of internal MC-TH coupling for dynMCNP must be further demonstrated.
Although much work has been realised, for application to full-core dynamic safety calculations many further improvements are needed.
Fig. 3. Power transient in a mini fuel assembly due to control rod movement
3.4 Work Package 4: Integration of high-performance parallel MC
3.4.1 Limitations of parallel execution in present Monte Carlo codes
A study has been carried out to investigate the best way for parallel execution of criticality calculations by Monte Carlo. To that end a simplified Monte Carlo program was designed that reflects the most essential elements of a Monte Carlo criticality calculation with successive fission source cycles. Both MPI (Message Passage Interface) and OpenMP (also called threading) are applied for parallel execution, as well as in combination. OpenMP uses shared memory and can best be used to run neutron history simulation in parallel on the processor cores of one computer node in a cluster. Then MPI can be used to run the job in parallel on different nodes, each using OpenMP within the node. It was found which keywords are necessary on the command line to run a job in parallel this way on a computer cluster.
The efficiency of the parallel execution is determined by the speedup factor, which is the ratio of the computing time in sequential mode and that in parallel mode. In the ideal case the speedup factor is equal to the number of processor cores used over all computer nodes. It was found that the speedup factor severely deviates from the ideal line when using larger numbers of processor cores (20 to 50). An important cause of decrease of efficiency is the amount of data to be sent from all slaves to the master processor after each criticality cycle. This concerns mainly the variables of the newly generated fission neutrons as source for the next cycle. Hence this number of data words must be minimised. A special case can be created if the fission source data is not send to the master processor and redistributed among the slaves, but each slave continues for the next cycle with the fission neutrons generated on that slave and only sending results in terms of multiplication factors to the master at the end of the calculation. This case has the disadvantage that the results cannot be reproduced when using a different number of processor cores. Notwithstanding the advantage of this method in reducing the amount of communication, the speedup still deviates seriously from the ideal line when using tens of processor cores. This is possibly caused by the relatively slow communication between different nodes in an common computer cluster. This may be much improved on a true supercomputer.
3.4.2 Improvements for parallelisation
In order to improve the parallel execution of MCNP, considerable modifications were introduced to have the master processor also simulate neutron histories, while it is originally used only for collecting and redistributing data from and to the slaves. When using a combination of MPI and OpenMP on multiple nodes, this would mean that a full computer node is used only for communication and would remain idle for most of the time.
Other improvements concern the communication between master and MPI slaves in sending and receiving only those parts of the fission bank with all starting positions and initial energies of the fission neutrons for a new cycle (neutron generation) which are relevant for that slave task. Also sending back to the master all tally data accumulated at each slave during a single (active) cycle can be limited to sending it only at the end of the calculation after processing the required number of cycles. This restricts the possibilities of intermediate intervention during a calculation on basis of the gathered tally information (for instance, detection of a sufficient level of statistical accuracy), but will improve the parallel speedup.
Each Monte Carlo run must start from a specified source distribution. As the distribution for a criticality problem is not known in advance, one can start inactive cycles from a spatial distribution specified in the input. After running one or more criticality cycles the fission source distribution can be stored in a file (srctp file) and can be used later on as a starting source. If such a (not necessarily completely converged) source distribution is used for later runs this will help avoiding a number of inactive cycles.
During initialisation of a job considerable parallel computer time can be saved by preparing an initialisation file with all relevant data for the run (geometry, composition, cross section data, number of cycles, number of histories per cycle, etc.). Especially collecting the relevant cross section data in the initialisation phase can be time consuming. The initialisation is only done on the master processor, after which the data is sent to the slaves. In MCNP a file (runtpe file) can be produced for later continuation of the job. With some modifications this option can be used to produce the initialisation file with all relevant data, but without the array for tally bin data (which can be considerably large and only contains zeroes at initialisation) and without the fission source data, which can be made available in a separate file.
If reading and writing files on or to a user disk is relatively slow, a run can better be executed using a fast work or scratch disk on which the relevant input files must be stored before starting the Monte Carlo run.
A straight-forward way of improving the parallel efficiency is to run successive cycles on each slave without redistribution of source data. The ultimate form is typically useful for running on massively parallel supercomputers and is discussed in the next section.
Investigations into the possibilities of load balancing revealed that using OpenMP includes a form of load balancing when using the PARALLEL construct of OpenMP instead of the PARALLEL DO construct. Load balancing for different computer nodes using MPI will not be efficient in view of the data exchange needed between processors on different nodes and the small losses that will occur when using OpenMP for the processor cores of the same node.
The possibilities of hyper-threading were studied but found to be restricted to Intel processors and to give probably only limited improvements in performance. It was therefore decided not to implement hyper-threading into the Monte Carlo codes. At some computer systems or supercomputers there is a possibility of SMT (Simultaneous Multiple Threading) in which case the processor determines itself possibilities to execute parts of tasks where execution space is available. In such cases SMT can be activated by specifying for OpenMP more threads than processor cores are available. In some cases a factor of 2 in speedup could be realised on 4-way SMT processors.
In the SERPENT Monte Carlo code from version 2 on the possibility was implemented to allow for OpenMP, possibly in combination with MPI. Parallelisation was introduced in later updates for burnup calculations.
The deliverable D4.6 describes in detail the MPI and OpenMP interfaces for parallel programming, and how they can be applied in a Monte Carlo program. Also the combination of those method is discussed. Examples are shown of the speedup factor that can be obtained by using an increasing number of processor cores on nowadays standard Linux computer clusters. The modifications to MCNP and SERPENT are also discussed.
3.4.3 Use of supercomputers
As Monte Carlo calculations for a full-core reactor will be very time consuming, it was anticipated in the project to get experience in running jobs on a supercomputer and see whether the developed methods can be run efficiently on a large parallel computer.
To this end an application was submitted to the PRACE organisation for European supercomputing to get access to the new JUQUEEN supercomputer in Jülich, Germany. Although the initial application was rejected, a second application with emphasis on the large number of processors needed for the final goal of the project was granted. It took a considerable number of working days to get actual access to the supercomputer with all cross section data copied. Furthermore, running a job with MCNP on the supercomputer also took a long preparation time. Afterwards many improvements could be tested and for relatively simple Monte Carlo calculations reasonable speedups for up to thousands of processor cores could be attained. However, for more realistic cases the speedup deteriorated at higher number of processor cores.
Therefore, a new Preparatory Access proposal was later submitted, now for the SuperMUC supercomputer in Garching, Germany. Here the method of running almost independent jobs on different computer nodes was tested. This method was developed locally and consists of modifications (in this case to MCNP, but generally applicable) to start a job on many nodes with MPI but using MPI only to initiate identical runs on different nodes without further communication with each other. As each MPI task has automatically a different task number, this number can be used to advance the random number generator for each task in a different way to perform different neutron history simulations. Each MPI task will generate successive criticality cycles and accumulate power distribution data, which are written to a binary file at the end of the calculation. After all MPI tasks are finished the power tally data from each node can be accumulated, possibly by a cheap separate job in serial execution, to a single file with all accumulated tally data. It will still be advantages to use threading (OpenMP) for all available processor cores on a node, but MPI for all computer nodes.
It turned out that the running time per independent node can differ seriously due to different histories and cycles, which decrease the parallel efficiency, as the MPI job on the master will only end if all other MPI tasks are finished. This problem can be solved by not specifying a fixed number of cycles, but a certain computer time. Although this option is present in MCNP, it needs modification to account for the wall-clock time and not the CPU time. Then it works perfectly well to end all MPI tasks at the same time. MCNP checks after each history whether the maximum time is reached or exceeded and then stops the current cycle. This is no problem for the total number of (active) histories that has contributed to the power distribution.
3.4.4 Results and Conclusions
Fig. 4 shows the results for a full-core calculation with 241 PWR fuel assemblies without thermal-hydraulic feedback specifying a running time of 119 min for the Monte Carlo history simulations for each case and taking into account the extra time needed for distributing the initial data to all nodes and accumulating the tally results from all nodes with increasing number of processor cores. Each computer node comprises 16 processor cores, which are used for threading on the same node. From the figure one can see that the speedup factor relative to using one node is almost equal to the number of nodes used, which means that there are hardly any losses in parallel execution.
Using 32,768 processor cores (2048 computer nodes) for 2 hours, 12.9 x 109 histories could be simulated with the power distribution estimated for each fuel pin subdivided over 100 axial subzones, which gives a total of about 7.0 106 different tally bins (including the zero tallies from the guide tubes instead of a fuel pin). Using the UFS (Uniform Fission Sites) variance reduction technique the project target of reaching a standard deviation of less than 1 % in 24 h computing time could have easily be attained if budget was available for a longer job.
Fig. 4. Parallel speedup factor versus number of computer cores
3.5 Work Package 5: Dissemination and training
The activities in this work package are aimed at bringing the project and its goals to the attention of possibly interested parties. On the other side, activities oriented to disseminate the knowledge generated within the project to “potential users” of the HPMC codes and methods but also to the academia were developed.
3.6 Work Package 6: Management
Work Package 6 regards the project management. Progress and achievements are detailed in Sect. 5.
Potential Impact:
Potential impact: The expected impact of the the obtianed results if the HPMC project are in line with the Horizon 2020 call for EURATOM fission devoted to improved safety design and operation of fission reactors that reflects the Strategic Research Agenda/Deployment Strategy of Sustainable Nuclear Energy Technology Platform (SNETP).
Hereafter a sit of specific expected impacts is given:
- provide the knowledge basis for developing robust national and EU policies in the field of nuclear reactor safety (short term): The project will deliver improved and validated numerical simulation tools to provide reference solutions to deterministic codes and to be used for safety demonstration
-helping interacting with stakeholders and civil society on nuclear reactor safety (short term): Numerical tools and methods developed within HPMC can be used by different end-users (industry, regulators, research centres, etc.) and since they reflects the state-of-the-art it will be a basis for technical discussions between utilities and regulators aiming to keep plant safety at high level
-further improve the safety of nuclear reactors in Europe and worldwide through increased resistance of safety relevant equipment and better safety culture (medium term): HPMC tools are essential to design reactor systems with improved safety features keeping sufficient safety margins. Disseminating most advanced safety analysis tools is a key factor to improve the safety culture in Europe and worldwide.
-demonstrable progress towards competitiveness of existing nuclear installations: Because of the verification possibilities, designers can adapt their design codes and have to spend less manpower and investments in verification efforts and will therefore be more competitive. Because of acceptable smaller safety margins in the design, the reactors will also be more competitive.
Socio-economic impacts:
The HPMC project provided improved and more advanced numerical simulation tools for the design and optimisation of reactor cores as well as for the assessment of the core performance and safety assessment for any core design. SERPENT and SUBCAHNFLOW are European tools being used in a number of European institutions. For example, the SERPENT code is being used in 16 European countries by about 140 users (worldwide: 300 users in 114 organisations of 31 countries). SUBCHANFLOW is used in different research organisations (e.g. CIEMAT, JRC-IE, HZDR, Demokritos Research Center) and universities (UPM, RWTH Aachen) and Small-Medium-Enterprises(e.g. DNC). This large number of European code users will directly profit from the extended capabilities of the tools and hereby they will be able to better predict safety margins of reactor cores. By decrease of the engineering factors this can lead to better fuel exploitation and improvement of the fuel cycle economy. Finally, due to the universality of the numerical tools based on Monte Carlo codes, the codes and methods developed in the frame of the HPMC project can easily be applied to research reactors, subcritical systems and Gen-IV reactors in addition to LWRs. For this reasons, a considerable economic impact can be expected from the HPMC project outcome.
Main Dissemination and training Activities:
The activities in this work package are aimed at bringing the project and its goals to the attention of possibly interested parties. On the other side, activities oriented to disseminate the knowledge generated within the project to “potential users” of the HPMC codes and methods but also to the academia were developed.
a) Various activities
A publicly accessible website is setup with information about the project, its aims, structure and participants. A separate website for internal use to the project was established to exchange information and documents not yet ready for publication. In 2013 the organisation of the websites is changed and brought under the auspices of KIT. The current public website address is www.fp7-hpmc.eu. It contains a private part under “members Area” to which all partner can access using a password. There all deliverables, presentations at the Executive Board meetings, presentations in conferences, journal and conference papers as well as the materials of the training courses are collected and will be kept there for longer time.
As the project goals fit well in the Strategic Research Agenda of the SNETP (Sustainable Nuclear Energy Technology Platform) a poster with the project setup and goals was presented at the 3rd General Assembly of SNETP in Warsaw, November 2011. Another general poster about the HPMC project was shown at the FISA conference in Vilnius, 2013. A seminar about the project was given at KTH in Stockholm in March 2012. A presentation about the project was given during a project meeting of the EU NURISP project in Karlsruhe, April 2012.
As the EU project NURISP and its successor NURESAFE are closely related to the goals of the HPMC project, good relationships were developed with the NURESASFE project and its management. Permission was asked and obtained to use the results for various benchmark test used in the NURISP project for comparison with results to be obtained with Monte Carlo calculation in the HPMC project. Also appointments were made with the NURESAFE project to study as far as possible the same benchmark test problems for the methods to be developed in both projects. The HPMC project is directed to the Monte Carlo method for the neutron transport; the NURESAFE project uses deterministic methods. The coordinator of the HPMC project was invited to give a lecture about the status of the codes and tools of the HPMC project in a NURESAFE ExCOM Meeting. This opportunity was used to explore the common areas of interest and to intensify the exchange between both projects which are very complementary (NURESAFE ExCOM meeting 7.10.2013 in Paris).
An extended summary of the project was prepared for publication by the European Commission. The deliverable is included in Sect. 5 on Management activities as well as the Communication Action Plan.
b) Training seminars
Training seminars are scheduled for potential users of the software developed in the project and are therefore outstanding opportunities to disseminate the knowledge and tools developed in the project.
A first training seminar was conducted at KIT in Germany on 24 and 25 April 2014 with 6 speakers directly involved in the project giving 11 different presentations covering most of the subjects developed during the project. The training seminar was attended by 9 external people from 6 different companies and institutes like Areva (Germany), CEA (France), GRS (Germany), HZDR (Germany), CCFE (UK), and from KIT.
The second training seminar was conducted at KTH in Stockholm, Sweden, on 4 and 5 June 2014 with basically the same program as for the first seminar. However, various changes were made to the presentations due to feedback and experience from the first seminar. This training seminar was attended by 13 external people from 5 different companies and institutes like Westinghouse, Vattenfall, University of Uppsala, Warsaw University of Technology, and from KTH.
Deliverable 5.12 gives an extensive overview of both training seminars and the attendees. Considering the goals set for successful training seminars with at least 4 participants from at least 2 different companies/institutes for both seminars, we can conclude that the training seminars were very successful.
Apart from the training seminars, scientific results of research performed for the project were published in several articles in international scientific journals and in papers at international conferences. Various other publications are submitted or in preparation.
List of Websites:
www.fp7-hpmc.eu
Main contacts:
Dr. Victor Hugo Sanchez Espinoza (Project Coordinator)
e-mail: victor.sanchez@kit.edu
Tel: +49-721-608 22283
Prof. Dr. J. Eduard Hoogenboom (Project Secretary)
e-mail: hoogenboom@delftnuclearconsultancy.nl
Tel: +31-10-4508442
The general objective of the HPMC Collaborative Project is to develop and demonstrate the application of full-core Monte Carlo calculation for time-dependent safety analysis with thermal-hydraulic feedback and burnup using high-performance computing. These goals fit into the Strategic Research Agenda of the Sustainable Nuclear Energy Technology Platform (SNETP).
The project is organised into 4 research work packages as follows
-WP1 aims at optimum Monte Carlo - thermal-hydraulics coupling; its objective is to realise efficient coupling of the Monte Carlo codes SERPENT and MCNP with the thermal-hydraulic subchannel codes SUBCHANFLOW and FLICA4, suitable for full-core application
-WP2 aims at optimum MC-burnup integration; its objective is to realise an efficient integration of burnup calculations in the Monte Carlo codes SERPENT and MCNP, suitable for full-core application
-WP3 aims at time dependence in Monte Carlo calculations; its objective is to realise an efficient time-dependence option in the Monte Carlo codes SERPENT and MCNP for safety analysis, applicable for full-core calculations
-WP4 aims at integration of high-performance parallel Monte Carlo; its objective is the demonstration of full-core Monte Carlo based dynamic safety analysis with ultimate efficiency in parallelisation
Two other work packages are defined
-WP5 aims at dissemination and training; its objective is communication of the project to potentially interested parties, including presentation of two training seminars for potential users of the developed code system with hands-on training
-WP6 covers the management of the project
Advanced and innovatives methods and code extensions were developed to be able to perform high fidelity core simulations based on Monte Carlo codes coupled with thermal hydraulics. For the demonstration of the new capabilities of the developed methods, a benchmark for a mini reactor core was defined that permitted to do code-to-code comparison. It demonstrated the efficienty and accuracy of the coupled Monte Carlo and thermal-hydraulic calculations. Integrated burnup calculations with a new and stable iteration scheme showed consistent results even for relatively large time steps. The time-dependent reactor behaviour with delayed neutron precursors after a control rod movement could be demonstrated for a mini fuel assembly. The efficiency of parallel execution of a full-size reactor criticality calculation on a massively parallel supercomputer could be increased to almost 100 %. Two successful training seminars were given to interested people from industry and research institutes.It is concluded that the project forms a solid basis for the next step forward to practical safety analysis validation of nuclear reactors based on the Monte Carlo method.
Project Context and Objectives:
The general objective of the HPMC Collaborative Project is to develop and demonstrate the application of full-core Monte Carlo calculation for time-dependent safety analysis with thermal-hydraulic feedback and burnup using high-performance computing. These goals fit into the Strategic Research Agenda of the Sustainable Nuclear Energy Technology Platform (SNETP).
The project is organised into 4 research work packages as follows
- WP1 aims at optimum Monte Carlo - thermal-hydraulics coupling; its objective is to realise efficient coupling of the Monte Carlo codes SERPENT and MCNP with the thermal-hydraulic subchannel codes SUBCHANFLOW and FLICA4, suitable for full-core application
- WP2 aims at optimum MC-burnup integration; its objective is to realise an efficient integration of burnup calculations in the Monte Carlo codes SERPENT and MCNP, suitable for full-core application
- WP3 aims at time dependence in Monte Carlo calculations; its objective is to realise an efficient time-dependence option in the Monte Carlo codes SERPENT and MCNP for safety analysis, applicable for full-core calculations
- WP4 aims at integration of high-performance parallel Monte Carlo; its objective is the demonstration of full-core Monte Carlo based dynamic safety analysis with ultimate efficiency in parallelisation
Two other work packages are defined
- WP5 aims at dissemination and training; its objective is communication of the project to potentially interested parties, including presentation of two training seminars for potential users of the developed code system with hands-on training
- WP6 covers the management of the project
Four organisations participate in the HPMC project:
o Karlsruhe Institute of Technology (KIT), Germany, also project coordinator
o Delft Nuclear Consultancy (DNC), The Netherlands, also project secretary
o VTT Technical Research Centre of Finland (VTT), Finland
o Royal Institute of Technology (KTH), Sweden
WP1 was successful in developing coupling methods for Monte Carlo (MC) and thermal-hydraulics (TH) calculations for large-scale applications. Besides the existing external coupling scheme of the MC and TH code, an internal coupling scheme was developed in which the TH code is modified to become a subroutine of the MC code. This provides higher efficiency as no communication via input and output files between the MC and TH codes is necessary. A Python based framework PIRS has been developed to define in a simple way complex geometries and to prepare the input file for either the MCNP MC code or the SCF TH code. This framework can be used to generate input files for the MC code for a complete fuel assembly and a full reactor core.
A dynamic assignment of the cross section data for the required temperature in each subzone of the reactor system was developed and implemented, which no longer requires that each subzone that can have a different temperature, has to be defined separately in the geometry input to the MC code. This facilitates the MC-TH coupling in a practical way for large fuel assemblies and full-size reactor cores.
A theory was developed and implemented to optimise the number of neutron histories to be applied in successive iteration steps of the MC-TH coupling. The number of neutron histories increases almost linearly with the MC-TH iterations. At the same time the relaxation parameter for combining the power distribution obtained from successive MC iterations as input to the TH calculation is optimised.
Other improvements are the introduction of on-the-fly temperature treatment in the SERPENT code instead of stochastic mixing of nuclide cross section data at two temperatures and the speedup of fission source convergence in the Monte Carlo calculation with the aid of an approximate deterministic calculation of the fission source distribution, taking into account the proper temperature distribution. Also new power tally estimators based on a collision estimator were developed for the detailed power distribution in a core which are much more efficient than existing estimators and hardly take more additional computing time.
To demonstrate the current capabilities for MC-TH coupled calculations a benchmark was defined for a mini core reactor with 3x3 fuel assemblies containing both UOX and MOX PWR fuel assemblies. Various codes were used showing reasonable agreement in the detailed power and temperature distributions over the core.
In WP2 current algorithms for coupling MC calculations with burnup calculations are investigated with respect to their stability. It was found that most algorithms are only conditionally stable and often instable at the relatively large time steps applied in burnup calculations.
A new algorithm, the Stochastic Implicit Euler method, was developed and successfully tested. The algorithm was implemented in the SERPENT code. Extensive calculations show stable results even at large time steps.
In WP3 time dependence was introduced in the MCNP5 Monte Carlo code, including the generation and decay of delayed neutron precursors. To reduce the statistical error in the generated reactor power in successive time intervals, a method of forced decay of precursors in each time interval was implemented. Also other variance reduction methods like the branchless collision were introduced. Nonetheless the statistical variation in the length of the fission chains and therefore the generated power per starting neutron in a time interval remains large in heterogeneous systems.
In order to simulate reactor transients, an efficient method was implemented to model the movement of control rods at any speed during a calculation.
Thermal-hydraulic feedback was implemented for the time-dependent version of MCNP to create a fully dynamic Monte Carlo code. To let the time-dependent TH calculation take into account the heating history further extensions of both codes were necessary. For an external MC-TH coupling scheme a method of exchanging signals between the independently started MC and TH codes was implemented to arrange exchange of information at the proper time instances. With further modifications an internal coupling was also realised.
The current capabilities for time-dependent calculations are demonstrated with a power transient due to control rod movement in a mini fuel assembly. Further developments are needed to accommodate safety calculations for large systems.
In WP4 various ways for parallel execution of a Monte Carlo calculation using the MPI and OpenMP application programming interfaces were investigated and their efficiency measured in terms of the speedup factor. For application on large computer clusters with different computer nodes and multiple processors per node, the optimum combination of MPI and OpenMP was determined. Application of OpenMP was introduced in the SERPENT2 code. The MCNP code was modified to use all available processor cores for neutron history simulation.
To increase the efficiency of the Monte Carlo calculation further, the amount of data communicated between computer processor cores after each cycle in a criticality calculation could be minimised. Nonetheless, the speedup of the parallel calculation using many processors was still not satisfactory. Therefore, a rigorous method was developed for application to massively parallel computer systems which avoids almost all communication between different tasks running on different computer nodes.
Access to a supercomputer via the European PRACE organisation was obtained for two successive stages of the project. It could be shown that parallel speedup for a criticality calculation of a full-size core was excellent when using the wall-clock time as a limit for criticality calculations instead of the number of cycles.
As part of the activities in WP5 a public website is setup with information about the project to bring the project and it goals to the attention of possibly interested parties. A separate website for internal use to the project was established to exchange information and documents between participants. A poster about the project was presented at the general assembly of the Sustainable Nuclear Energy Technology Platform (SNETP) in Warsaw and the FISA conference in Vilnius. Presentations about the project were given on several occasions.
With the EU project NURISP and its successor NURESAFE agreements were made about the use of benchmark problem definitions and results and about common problems to be investigated at both sides.
An extended summary of the project was prepared for publication by the European Commission. For dissemination of the results of the project a Communication Action Plan was written.
Two successful training seminars were organised in Germany and Sweden, respectively, which were attended by in total 22 external people from 11 different companies and research institutes.
Scientific results of research performed for the project were published in several articles in international scientific journals and in papers at international conferences.
WP6 covers the management of the project. The Executive Board gathered every half a year. In the board meetings the progress of work was discussed by means of presentations by all participants involved in the project. The work package leaders presented the status and outlook of the work in their work package.
All deliverables scheduled for the project were published, some of which with some delays.
It was concluded at the last Executive Board meeting that basically all tasks were executed and that important new and advanced developments were realised and successful results demonstrated for most of the tasks. Therefore the HPMC project forms a solid basis for the next step forward to practical safety analysis of nuclear reactors based on the Monte Carlo method.
Project Results:
The research and development work of the project is organised in the first four work packages. Work packages 1 to 3 cover the development of the essential items for realistic reactor safety analysis with Monte Carlo calculations, namely coupling of Monte Carlo with thermal-hydraulics (WP1), integration of burnup calculations (WP2) and time-dependent Monte Carlo methods from seconds to minutes scale (WP3).
Work package 4 is directed to optimum integration of the techniques developed in the first three work packages, to maximum efficiency in parallel execution of the Monte Carlo codes and demonstration of the developed tools for a full-size reactor core.
Work package 5 is directed to dissemination of the results of the project and to training sessions for potential users.
Work package 6 concerns management of the project and is discussed separately in Sect. 5.
3.1 Work Package 1: Optimum MC-TH coupling
In the HPMC project the main Monte Carlo codes used are MCNP and SERPENT. The main thermal-hydraulics code is SUBCHANFLOW (SCF), developed at KIT. Incidentally the FLICA4 code from CEA, Saclay, France, is used. Where necessary these codes were obtained and installed. For new codes at a partner’s institution, an extensive familiarisation process is passed through as not only the structure of the input and output files must be known, but also the structure of the code itself. Especially the SERPENT and SCF codes are under continuous development, which requires regular installation of updates.
3.1.1 General considerations and developments on MC-TH coupling
Coupled Monte Carlo (MC) and thermal-hydraulics (TH) calculations require a series of iterations as the TH calculation needs the detailed power density distribution from the MC code, while the MC code needs the temperature distribution from the TH code. Neither distributions are known on beforehand, so one has to start with a guessed power distribution for the TH code, which calculates the temperature (and coolant density) distribution. This distribution must be input to the MC calculation, from which an improved power distribution is obtained. Then successive iterations will, in general, lead to a converged solution for both the power distribution and the temperature distribution.
For the coupling of Monte Carlo calculations with thermal-hydraulics calculations two schemes are followed, external coupling and internal coupling. External coupling uses both the Monte Carlo and thermal-hydraulics code as such and provide communication between them by retrieving the necessary data from output files produced by the MC code or the TH code and integrating these data in the input file of the TH code or MC code, respectively, for the next iteration. This way of coupling has the advantage that the original codes can be used and in fact only an executable version of the codes is needed. With proper conversion programs for converting the output of one code to input to the other one, various MC and TH codes can be coupled.
With internal coupling the TH code is made a subprogram of the MC code and is called as a subroutine. It needs modification of both the MC and TH code, but provides a simpler and faster way of data exchange between the MC and TH parts.
For a coupled MC-TH calculation the MC code needs not only to deal with the detailed geometry description of a reactor with many fuel assemblies and many fuel pins in a fuel assembly, but also with the different cross section data needed for each axial subzone in each fuel pin because of the different temperatures for fuel, cladding and coolant all over the reactor. To facilitate the generation of an input file for the MC codes which includes not only the geometry description, but also the assignment of cross section data at different temperatures for each axial subzone, a sophisticated system called PIRS, based on the Python language is developed.
Coupling schemes are developed and implemented at different participating institutes, which gives the possibility to evaluate different implementations and select the most promising developments for the future. External coupling has been realised for MCNP-SCF, MCNP-FLICA, Serpent-SCF and as a side product TRIPOLI-FLICA. Internal coupling has been realised for MCNP-SCF, Serpent-SCF and Serpent-COSY (a light-weight TH code developed at VTT).
Deliverable D1.4 Parts I and II, gives a description of the MC-TH coupling schemes as developed by the various partners in the first stage of the project. Afterwards further developments took place, especially to improve the internal coupling schemes on different items.
As the cross section data for neutron interaction with material are temperature dependent, each subzone in the Monte Carlo calculation must be treated at a specific temperature according to the outcome of the TH calculation. As all temperatures in a wide range will occur in the reactor, it is not realistically possible to supply neutron cross section data at all these temperatures for all relevant nuclides. Hence the technique of stochastic mixing is further developed to solve this problem. Here for each nuclide present in a medium two references to different temperature cross section data are defined with such atom fractions that on the averaged the desired temperature data is selected in simulating histories with selecting a reacting nuclide proportional to its atom fraction. With a variation this technique can also be applied for treating S(α,β) scattering for thermal scattering to bound atoms in a molecule like H in H2O. However, also a technique for direct interpolation into the S(α,β) kernel has been developed.
The stochastic mixing technique still requires a separate composition definition for each axial subzone of each pin cell, as all temperatures may differ. For large systems like a fuel assembly and certainly a full reactor core, this requires a huge number of separate subzones that cannot be handled in practice. One solution is the development of the PIRS system (Sect. 3.1.2). Another way to handle this problem is the dynamical assignment of nuclide cross section data at the applicable temperature when a neutron enters a certain zone during its history simulation. This keeps the geometry description in the input to the MC code relatively simple and limited and provides a major step forward in MC-TH coupling for a full reactor core. The methods were fully developed and demonstrated to give good results.
Another way of treating the temperature dependence with respect to cross section data is “on-the-fly” temperature treatment of neutron cross sections, which is developed for the SERPENT code. This implies that each time a neutron at a given energy needs the cross section at the temperature of the material the neutron is in, the code samples a relative energy taking into account the thermal motion of the nucleus at the specified temperature. This method takes some additional computation time, but avoids several problems otherwise met for temperature treatment. At the moment it still has some limitations.
An important issue in optimisation of MC-TH iterations is the control over the number of neutron histories to simulate in the successive Monte Carlo runs. For the first few iterations the power profile and thus the temperature profile may be far from converged and using too many neutron histories will be a waste of computer resources. In WP2 a procedure based on a physical theory is derived for optimised control of the number of histories in successive iterations, which is not only applicable to burnup calculations, but also to MC-TH coupled iterations. The number of histories turns out to increase almost linearly with the iteration number until convergence. This theory also provides an optimised relaxation parameter for the power profile to include the results for the power profile from previous iterations. Tests show that these procedures lead to a converged result.
Besides the above-mentioned subjects, various other items are studied and methods are developed. One issue is the variance reduction in power estimation in subzones all over the reactor using the Uniform Fission Sites (UFS) method.
Another issue is the convergence of a Monte Carlo calculation with respect to the fission source distribution. To speed up the convergence it is important to start the Monte Carlo calculation with a fission source distribution as close as possible to the fundamental mode distribution, taking into account the geometry of the reactor system and the temperature distribution. To this end coupling with a deterministic nodal method is developed, taking into account the temperature distribution obtained from the TH code. The deterministic code calculates the fission source distribution with as few approximations to the geometry as possible, which is then used as the starting source for the first Monte Carlo source convergence cycle. This method shows promising results.
3.1.2 PIRS - Python Interfaces for Reactor Simulations
As the geometry description of a reactor system in a Monte Carlo input file can be rather complicated, it is important to have a tool that facilitates the generation of such an input file. Therefore, an extensive set of Python based modules with definition of various classes and class methods is developed, which allows to define geometry and code-specific data for MCNP and SCF in terms of the Python programming language. It also gives the possibility to start the codes, so iterative MC-TH calculations can be organised. Moreover, It can provide plots of the temperature and power distributions as calculated by the TH and MC codes.
The Python based framework is operational for a single pin cell geometry, a realistic fuel assembly with a large number of fuel pins and for a complete core. It can generate inputs for MCNP and Serpent, as well as for the coupling of these codes with the SCF thermal-hydraulics code.
However, for larger systems, for instance the mini core benchmark as discussed in Sect. 3.1.5 the generated input file for the MC code is extremely large as the MC input file is based on defining each subzone separately because of the varying temperatures. It takes MCNP a very long time to process the input file at initialisation of a run.
3.1.3 External MC-TH coupling
External coupling has the advantage that it can use, in principle, both the MC code and the TH code in standard form and basically only an executable version of these codes is needed. A script is needed to call successively the TH and MC code. Coupling between the codes is by output and input files. As the input for one type of calculation (either MC or TH) is totally different from the data provided in the output of the other code (TH or MC, respectively), separate conversion programs are needed to read the output file(s) from for instance the TH code and prepare from it a new input file for the MC code, using a standard input file for the MC code in which relevant data must be updated from the results of the TH code. The same must be done for preparing a new input file for the TH code.
Although the conversion programs for reworking the output from one code to a new input for the other are dependent on both the MC and TH code used, it was possible to couple in this way different MC codes (MCNP, Serpent, TRIPOLI) to different TH codes (SCF, FLICA). However, the communication via input and output files and the conversion of these files is rather inefficient and limits parallel execution to the separate runs of the MC code. For larger problems the output files (power distribution and temperature distribution, respectively) can be very large for a large number of subzones. Hence, this type of calculation is useful for comparison of codes and code combinations for smaller problems, but is not suitable for calculating large systems.
3.1.4 Advanced internal MC-TH coupling
As can be concluded from the previous sections, internal coupling of MC and TH calculations is the most efficient way and the only one to deal with large systems like a full-size reactor core. Hence, in the project, after exploitation and testing of various ways for coupling, most attention was paid to the development of efficient methods of internal MC-TH coupling.
As internal coupling requires modifications to the TH code (as well as to the MC code), the KIT code SubChanFlow (SCF) was the only TH code for which the source code was available to participants and could be adapted to become a subroutine that could be called by the MC code. Also because of the structure of the SCF code this was not a problem and various arrays with temperatures for the fuel, cladding and coolant and with the coolant density for all subzones can be made accessible in the MC code. The MC code also needs modifications to call the subroutine which is the start of the TH calculation.
For the MC-TH iterations several solutions are possible, but most efficient is to modify the MC code to include internally a loop which successively calls the TH subroutine and the subroutine starting the Monte Carlo neutron history simulation. Then the geometry, material composition and cross section data needs to be read only once. For the TH code also modifications are necessary to read the input only once and to skip certain parts which need not to be recalculated each iteration and may only reinitialise internal values that can be useful to keep for continuation the internal TH iterations when a new MC-TH iteration must be performed. The power distribution per subzone, which is needed in the TH calculation, can be put in the corresponding array in the TH code from the power tally results in the MC code after some reworking and renormalisation.
After execution of the TH calculation the optimised scheme for increasing the number of neutron histories per MC-TH iteration as mentioned in Sect. 3.1.1 is implemented in the MC code in the loop for MC-TH calculations as well as a procedure to test the convergence of the MC-TH iteration. To accommodate the different temperatures in all subzones in the reactor system, stochastic mixing is applied in the standard input to the MC code. For dynamic temperature linking to cross section data the subroutine providing the cross section data at the current neutron energy when entering (or starting in) a new medium is modified and the atom fractions of the nuclides for stochastic mixing are temporarily adapted to the temperature for that medium.
Other modifications to improve the calculation, for instance for speeding up the convergence of the fission source distribution, for improved estimators for the power distribution and for global variance reduction for estimating the detailed power distribution are discussed in Sect. 3.4.
Another important improvement concerns the introduction of a new estimator for the detailed power distribution based on a collision estimator. This reduces the additional time for tallying considerably. An advanced version did not even take noticeably more computing time.
Apart from the demonstration shown in the next section for a mini core reactor, successful coupled MC-TH calculations could be performed for a quarter core, which shows the capacities of the developed methods.
3.1.5 Demonstration and Conclusions
To demonstrate the capabilities of the developed MC-TH code coupling a computational benchmark was defined to obtain and compare numerical results. The benchmark definition and the results are detailed in deliverable D1.9. The benchmark problem consists of a mini-core of 3x3 PWR fuel assemblies, from which 6 use UOX fuel and 3 use MOX fuel with 3 different Pu concentrations. Each fuel assembly has 17x17 fuel positions from which 25 are occupied by guide tubes (without control rods) for the UOX assemblies and with WABA (Wet Annular Burnable Absorber) rods for the MOX assemblies. Hence, the geometry and composition is quite complicated. Although the project definition required the demonstration for a single fuel assembly, it was felt that the project was up to the treatment of a larger system.
Four independent solutions were obtained: two with internal coupling of MCNP-SCF (but different implementations), one with external coupling of MCNP-SCF using the PIRS system for generating the input files and one with internal coupling of Serpent-SCF.
Although the results show some differences, as discussed in the deliverable report, and further detailed analysis is desired, it can be concluded that various developments have led to coupled code systems capable of calculating a (mini) core system with realistic geometry and detailed fuel and coolant temperature and coolant density distribution throughout the core.
Fig. 1 shows from top to bottom the radial distribution through the mini core of the fuel temperature, the coolant temperature and the coolant density. The figures at the left are at the bottom of the fuel, the figures at the right are near the top of the fuel. Red colours indicate higher values, blue colours indicate lower values. The fuel assemblies with mainly light and dark blue are the MOX fuel assemblies.
Fig. 1. Radial distribution of fuel temperature, coolant temperature and coolant density at bottom, middle and top of fuel
3.2 Work Package 2: Optimum MC-burnup integration
Although there exists various examples of Monte Carlo codes that can be coupled to burnup calculations using different algorithms, the stability of the algorithms are not sufficiently investigated. Therefore, a simplified model was devised to test various coupling schemes. The model is based on the fission matrix concept and imitates reasonably well the features of standard Monte Carlo calculations, but at much lower computational cost. Moreover, Xenon feedback is implemented, which can play an important role in instabilities. As this model allows to run the simulation also in deterministic mode, it can determine exactly the statistical errors due to the Monte Carlo simulation.
Extensive calculations with this model showed that most of the applied algorithms for burnup feedback are unstable when relatively large time steps are used, as is often the case in burnup calculations.
3.2.1 A Stable MC-burnup scheme
On basis of these results a new coupling scheme was designed and successfully tested. The fully implicit scheme is robust and efficient. It is called the Stochastic Implicit Euler (SIE) method. The results of this research are detailed in deliverable D2.3 and published in two journal articles. The results are both used for burnup calculations with the SERPENT code and for coupled Monte Carlo – thermal-hydraulics calculations (see Sect. 3.1).
Fig. 2 shows the axial offset and spatial distribution of 135Xe as a function of burnup. The top rows are without equilibrium 135Xe, the bottom rows with equilibrium 135Xe calculation. For the figures from left to right the step size in the calculation is 10, 5, 2.5 and 1MWd/kgU.
Fig. 2. Results for axial offset and spatial distribution of 135Xe as a function of burnup indicated by the number of depletion time steps
3.2.2 An optimised scheme for MC-TH iterations
From the theory of stochastic approximation an optimised scheme for MC-TH iterations is developed. It prescribes the increase in number of histories (or number of cycles) in Monte Carlo runs for successive MC-TH iterations. As the inaccuracy of the temperature distribution over the reactor will decrease with the MC-TH iterations, it pays off to use more cycles and neutron histories for better statistical accuracy in later MC-TH iterations. At the same time the optimised scheme determines the relaxation of the power distribution from successive MC-TH iterations with the relaxation factor in each iteration determined by the number of active histories relative to those in the first iteration. This scheme was successfully used in most of the MC-TH coupled calculations discussed in Sect. 3.1.
3.3 Work Package 3: MC time dependence
Although time dependence is almost automatically included in Monte Carlo codes, it is normally limited to the timescale of the neutron life time, which can vary from microseconds to about 1 ms, depending on the reactor type. For reactor transient analysis with Monte Carlo the generation and decay of delayed neutron precursors must also be modeled. In long-time time-dependent Monte Carlo calculations neutron chains of successive fission neutrons, initiated by the decay of a precursor nucleus, play a dominant role. Although such chains can contain hundreds of generations of fission neutrons, their life time is normally short compared to the physical phenomena generally occurring in transient situations.
A study showed that the statistics in neutron chain length are large and therefore the generated power by a chain of fission neutrons, too. This requires special techniques when implementing time dependence in a Monte Carlo code.
In the SERPENT2 code the time-dependent transport model was implemented, but without delayed neutron precursors yet. Results for short-time time dependence are published in a conference paper. All other developments for time dependence were implemented in MCNP which is then called dynMCNP.
3.3.1 Implementation of time dependence with delayed neutron precursors
Generation and decay of delayed neutron precursors was implemented in MCNP. Like fission neutrons the precursors are stored in a bank in computer memory, but with a flag to distinguish them from prompt neutrons. To describe the produced power in a reactor system during a transient, time intervals are introduced. Internally in dynMCNP time intervals are treated like successive criticality cycles, that is, all source particles for a certain time interval are successively retrieved from the particle bank and their history simulated. In a collision the result can be a capture event, which ends the history of that neutron, a scattering event with continuation of the neutron history or fission (in a fissionable medium) in which case the history of the neutron causing fission is ended and n fission neutrons are generated from which n-1 are stored in a temporary bank and the history of one fission neutron is continued. When this history is ended, possibly after more fission events, the last fission neutron is retrieved from the temporary bank. If a neutron reaches the end of the time interval, it is stored in the particle bank as prompt neutron for the next time interval. If in a fission event a new precursor is generated, it is also stored in the particle bank, marked as a precursor. When all particles for the current interval are processed execution for the next interval will be started.
To reduce the statistical uncertainty in the produced power per time interval, existing precursors are forced to decay in each time interval with a correction for their statistical weight. Another general option for variance reduction was implemented: branchless collision. In that case always one neutron will result from a collision. Therefore, capture is treated implicitly by modifying the weight of the particle. However, in our case the treatment of fission is also modified in such a way that only one new fission neutron is generated. This requires further adaption of the neutron weight, which can possibly increase in fissionable media. A neutron history can only be terminated by escape from the system or by Russian roulette if its weight becomes too low.
In order to simulate transients due to the movement of control rods, an efficient method was implemented to determine for each neutron where the tip of the control rod is located, depending on the time variable of the neutron. It is possible to move control rods together as a “bank”, as will be the case in actual reactors. In additional input to the MC code the movement of a control rod can be described. The control rod can move linearly in time during each time interval. As the introduction of temperature dependence with an externally given temperature distribution (for instance as a result of a thermal-hydraulics calculation) per subzone with dynamic temperature linking (see Sect. 3.1.1) requires the definition of a lattice in MCNP for axial division of subzones, this prohibits the definition of cell for a control rod from the control rod tip to the top of the reactor and of a cell below the control rod tip with coolant down to the bottom of the reactor. In order to avoid multiple cell definitions of a control rod over various axial subzones, which would make modelling a movable control rod very difficult, subroutines of MCNP were modified that bypass the axial lattice in case of specific cells limited by surfaces used for determining (possibly moving) control rod tips.
Efficient parallel execution of the Monte Carlo code for the simulation of neutron histories in a single time interval requires additional measures. In order to reproduce results of a calculation with a different number of computer processors, the precursors generated during the interval and the prompt neutrons reaching the end of the time interval must be stored in a bank separate from the precursors and prompt neutrons present at the start of the time interval and must be sent separately to the master processor in order to include them in the correct order in the bank for particles starting in the new time interval.
3.3.2 Coupling with a time-dependent thermal-hydraulics code
Coupling of a time-dependent simulation with thermal-hydraulic feedback was implemented for dynMCNP to create a fully dynamic Monte Carlo code. The external MC-TH coupling scheme as developed in WP1 had to be extended, as for a time-dependent calculation the TH calculation must take into account the heating history. Therefore, both the MC and the TH codes run continuously and a method was developed of sending each other signals to indicate when they have finished the computational work for a certain time interval. When the TH code has finished the processing for a time interval the MC code receives the signal, collects the results of the TH code for temperature and coolant density distributions and starts the simulation of neutron histories for that time interval using the dynamic assignment of nuclides at the applicable temperature as developed in WP1. When the MC code has completed the work for the time interval the TH codes receives a signal from the MC code, collects the estimated power distribution from the MC code and starts the calculations for the next time interval. When one code is actively doing calculations, the other code goes asleep, waiting for a signal to awake and take over the computation. The SUBCHANFLOW (SCF) code was modified to send and receive signals for the dynamic MC-TH coupling. The coupling method was demonstrated to work properly. However, due to still consisting problems with large statistical uncertainty in the generated power in a system, it was not possible to produce a convincing demonstration of a time-dependent case due to a moving control rod with thermal-hydraulic feedback. The deliverable D3.5 describes in detail the time-dependent external MC-TH coupling mechanism, but without a demonstration calculation of a time-dependent case with thermal-hydraulic feedback.
As internal coupling of Monte Carlo and thermal-hydraulics calculations will be more efficient for large reactor systems using many processor cores on a computer cluster or supercomputer, internal coupling of dynMCNP with SCF was also developed. As in internal MC-TH coupling it requires a loop in MCNP which successively calls the routine in MCNP for processing particles for a certain time interval and then calling SCF for execution of the time step in the time-dependent version of MCNP. This requires further modification of SCF to step out of SCF after each time interval and to continue its time-dependent calculation at the next entrance, without recalculating several data. Also dynMCNP requires further modification when re-entering the simulation routine for the next time interval to continue with all data gathered at the end of the previous interval.
3.3.3 Results and Conclusions
Results for a time-dependent case (without thermal-hydraulic feedback) were obtained for a small 3x3 fuel assembly in which the central rod is replaced by a movable control rod. The outer boundaries of the system are reflective. The transient is introduced by linear movement of the control rod tip from 0.2 s to a new position at 0.5 s, after which the position was fixed. The total reactivity effect was about 0.6 $. Intervals of 0.1 s were considered. The power over the assembly was recorded per time interval during about 1.5 s. Fig. 3 shows the resulting power transient. For comparison a point-kinetics calculation was made with the estimated total reactivity from a separate static MCNP calculation. For this rather small system there is a fair agreement between both calculations.
It can be concluded that the implementation of long-time time dependence with generation and decay of delayed neutron precursors is successful. However, there are several items that need further attention and improvement. Firstly, the statistics in the recorded total power are considerable. Therefore, the calculation required simulation of as many as 106 neutrons histories per time interval. Secondly, the efficiency of the parallel execution of the Monte Carlo runs was rather low, which restricts the possibilities to scale up to massively parallel execution on a supercomputer. This is very much needed as the calculation takes a long time due to the long neutron chains per time interval because of many successive fission events. Thirdly, the implementation of internal MC-TH coupling for dynMCNP must be further demonstrated.
Although much work has been realised, for application to full-core dynamic safety calculations many further improvements are needed.
Fig. 3. Power transient in a mini fuel assembly due to control rod movement
3.4 Work Package 4: Integration of high-performance parallel MC
3.4.1 Limitations of parallel execution in present Monte Carlo codes
A study has been carried out to investigate the best way for parallel execution of criticality calculations by Monte Carlo. To that end a simplified Monte Carlo program was designed that reflects the most essential elements of a Monte Carlo criticality calculation with successive fission source cycles. Both MPI (Message Passage Interface) and OpenMP (also called threading) are applied for parallel execution, as well as in combination. OpenMP uses shared memory and can best be used to run neutron history simulation in parallel on the processor cores of one computer node in a cluster. Then MPI can be used to run the job in parallel on different nodes, each using OpenMP within the node. It was found which keywords are necessary on the command line to run a job in parallel this way on a computer cluster.
The efficiency of the parallel execution is determined by the speedup factor, which is the ratio of the computing time in sequential mode and that in parallel mode. In the ideal case the speedup factor is equal to the number of processor cores used over all computer nodes. It was found that the speedup factor severely deviates from the ideal line when using larger numbers of processor cores (20 to 50). An important cause of decrease of efficiency is the amount of data to be sent from all slaves to the master processor after each criticality cycle. This concerns mainly the variables of the newly generated fission neutrons as source for the next cycle. Hence this number of data words must be minimised. A special case can be created if the fission source data is not send to the master processor and redistributed among the slaves, but each slave continues for the next cycle with the fission neutrons generated on that slave and only sending results in terms of multiplication factors to the master at the end of the calculation. This case has the disadvantage that the results cannot be reproduced when using a different number of processor cores. Notwithstanding the advantage of this method in reducing the amount of communication, the speedup still deviates seriously from the ideal line when using tens of processor cores. This is possibly caused by the relatively slow communication between different nodes in an common computer cluster. This may be much improved on a true supercomputer.
3.4.2 Improvements for parallelisation
In order to improve the parallel execution of MCNP, considerable modifications were introduced to have the master processor also simulate neutron histories, while it is originally used only for collecting and redistributing data from and to the slaves. When using a combination of MPI and OpenMP on multiple nodes, this would mean that a full computer node is used only for communication and would remain idle for most of the time.
Other improvements concern the communication between master and MPI slaves in sending and receiving only those parts of the fission bank with all starting positions and initial energies of the fission neutrons for a new cycle (neutron generation) which are relevant for that slave task. Also sending back to the master all tally data accumulated at each slave during a single (active) cycle can be limited to sending it only at the end of the calculation after processing the required number of cycles. This restricts the possibilities of intermediate intervention during a calculation on basis of the gathered tally information (for instance, detection of a sufficient level of statistical accuracy), but will improve the parallel speedup.
Each Monte Carlo run must start from a specified source distribution. As the distribution for a criticality problem is not known in advance, one can start inactive cycles from a spatial distribution specified in the input. After running one or more criticality cycles the fission source distribution can be stored in a file (srctp file) and can be used later on as a starting source. If such a (not necessarily completely converged) source distribution is used for later runs this will help avoiding a number of inactive cycles.
During initialisation of a job considerable parallel computer time can be saved by preparing an initialisation file with all relevant data for the run (geometry, composition, cross section data, number of cycles, number of histories per cycle, etc.). Especially collecting the relevant cross section data in the initialisation phase can be time consuming. The initialisation is only done on the master processor, after which the data is sent to the slaves. In MCNP a file (runtpe file) can be produced for later continuation of the job. With some modifications this option can be used to produce the initialisation file with all relevant data, but without the array for tally bin data (which can be considerably large and only contains zeroes at initialisation) and without the fission source data, which can be made available in a separate file.
If reading and writing files on or to a user disk is relatively slow, a run can better be executed using a fast work or scratch disk on which the relevant input files must be stored before starting the Monte Carlo run.
A straight-forward way of improving the parallel efficiency is to run successive cycles on each slave without redistribution of source data. The ultimate form is typically useful for running on massively parallel supercomputers and is discussed in the next section.
Investigations into the possibilities of load balancing revealed that using OpenMP includes a form of load balancing when using the PARALLEL construct of OpenMP instead of the PARALLEL DO construct. Load balancing for different computer nodes using MPI will not be efficient in view of the data exchange needed between processors on different nodes and the small losses that will occur when using OpenMP for the processor cores of the same node.
The possibilities of hyper-threading were studied but found to be restricted to Intel processors and to give probably only limited improvements in performance. It was therefore decided not to implement hyper-threading into the Monte Carlo codes. At some computer systems or supercomputers there is a possibility of SMT (Simultaneous Multiple Threading) in which case the processor determines itself possibilities to execute parts of tasks where execution space is available. In such cases SMT can be activated by specifying for OpenMP more threads than processor cores are available. In some cases a factor of 2 in speedup could be realised on 4-way SMT processors.
In the SERPENT Monte Carlo code from version 2 on the possibility was implemented to allow for OpenMP, possibly in combination with MPI. Parallelisation was introduced in later updates for burnup calculations.
The deliverable D4.6 describes in detail the MPI and OpenMP interfaces for parallel programming, and how they can be applied in a Monte Carlo program. Also the combination of those method is discussed. Examples are shown of the speedup factor that can be obtained by using an increasing number of processor cores on nowadays standard Linux computer clusters. The modifications to MCNP and SERPENT are also discussed.
3.4.3 Use of supercomputers
As Monte Carlo calculations for a full-core reactor will be very time consuming, it was anticipated in the project to get experience in running jobs on a supercomputer and see whether the developed methods can be run efficiently on a large parallel computer.
To this end an application was submitted to the PRACE organisation for European supercomputing to get access to the new JUQUEEN supercomputer in Jülich, Germany. Although the initial application was rejected, a second application with emphasis on the large number of processors needed for the final goal of the project was granted. It took a considerable number of working days to get actual access to the supercomputer with all cross section data copied. Furthermore, running a job with MCNP on the supercomputer also took a long preparation time. Afterwards many improvements could be tested and for relatively simple Monte Carlo calculations reasonable speedups for up to thousands of processor cores could be attained. However, for more realistic cases the speedup deteriorated at higher number of processor cores.
Therefore, a new Preparatory Access proposal was later submitted, now for the SuperMUC supercomputer in Garching, Germany. Here the method of running almost independent jobs on different computer nodes was tested. This method was developed locally and consists of modifications (in this case to MCNP, but generally applicable) to start a job on many nodes with MPI but using MPI only to initiate identical runs on different nodes without further communication with each other. As each MPI task has automatically a different task number, this number can be used to advance the random number generator for each task in a different way to perform different neutron history simulations. Each MPI task will generate successive criticality cycles and accumulate power distribution data, which are written to a binary file at the end of the calculation. After all MPI tasks are finished the power tally data from each node can be accumulated, possibly by a cheap separate job in serial execution, to a single file with all accumulated tally data. It will still be advantages to use threading (OpenMP) for all available processor cores on a node, but MPI for all computer nodes.
It turned out that the running time per independent node can differ seriously due to different histories and cycles, which decrease the parallel efficiency, as the MPI job on the master will only end if all other MPI tasks are finished. This problem can be solved by not specifying a fixed number of cycles, but a certain computer time. Although this option is present in MCNP, it needs modification to account for the wall-clock time and not the CPU time. Then it works perfectly well to end all MPI tasks at the same time. MCNP checks after each history whether the maximum time is reached or exceeded and then stops the current cycle. This is no problem for the total number of (active) histories that has contributed to the power distribution.
3.4.4 Results and Conclusions
Fig. 4 shows the results for a full-core calculation with 241 PWR fuel assemblies without thermal-hydraulic feedback specifying a running time of 119 min for the Monte Carlo history simulations for each case and taking into account the extra time needed for distributing the initial data to all nodes and accumulating the tally results from all nodes with increasing number of processor cores. Each computer node comprises 16 processor cores, which are used for threading on the same node. From the figure one can see that the speedup factor relative to using one node is almost equal to the number of nodes used, which means that there are hardly any losses in parallel execution.
Using 32,768 processor cores (2048 computer nodes) for 2 hours, 12.9 x 109 histories could be simulated with the power distribution estimated for each fuel pin subdivided over 100 axial subzones, which gives a total of about 7.0 106 different tally bins (including the zero tallies from the guide tubes instead of a fuel pin). Using the UFS (Uniform Fission Sites) variance reduction technique the project target of reaching a standard deviation of less than 1 % in 24 h computing time could have easily be attained if budget was available for a longer job.
Fig. 4. Parallel speedup factor versus number of computer cores
3.5 Work Package 5: Dissemination and training
The activities in this work package are aimed at bringing the project and its goals to the attention of possibly interested parties. On the other side, activities oriented to disseminate the knowledge generated within the project to “potential users” of the HPMC codes and methods but also to the academia were developed.
3.6 Work Package 6: Management
Work Package 6 regards the project management. Progress and achievements are detailed in Sect. 5.
Potential Impact:
Potential impact: The expected impact of the the obtianed results if the HPMC project are in line with the Horizon 2020 call for EURATOM fission devoted to improved safety design and operation of fission reactors that reflects the Strategic Research Agenda/Deployment Strategy of Sustainable Nuclear Energy Technology Platform (SNETP).
Hereafter a sit of specific expected impacts is given:
- provide the knowledge basis for developing robust national and EU policies in the field of nuclear reactor safety (short term): The project will deliver improved and validated numerical simulation tools to provide reference solutions to deterministic codes and to be used for safety demonstration
-helping interacting with stakeholders and civil society on nuclear reactor safety (short term): Numerical tools and methods developed within HPMC can be used by different end-users (industry, regulators, research centres, etc.) and since they reflects the state-of-the-art it will be a basis for technical discussions between utilities and regulators aiming to keep plant safety at high level
-further improve the safety of nuclear reactors in Europe and worldwide through increased resistance of safety relevant equipment and better safety culture (medium term): HPMC tools are essential to design reactor systems with improved safety features keeping sufficient safety margins. Disseminating most advanced safety analysis tools is a key factor to improve the safety culture in Europe and worldwide.
-demonstrable progress towards competitiveness of existing nuclear installations: Because of the verification possibilities, designers can adapt their design codes and have to spend less manpower and investments in verification efforts and will therefore be more competitive. Because of acceptable smaller safety margins in the design, the reactors will also be more competitive.
Socio-economic impacts:
The HPMC project provided improved and more advanced numerical simulation tools for the design and optimisation of reactor cores as well as for the assessment of the core performance and safety assessment for any core design. SERPENT and SUBCAHNFLOW are European tools being used in a number of European institutions. For example, the SERPENT code is being used in 16 European countries by about 140 users (worldwide: 300 users in 114 organisations of 31 countries). SUBCHANFLOW is used in different research organisations (e.g. CIEMAT, JRC-IE, HZDR, Demokritos Research Center) and universities (UPM, RWTH Aachen) and Small-Medium-Enterprises(e.g. DNC). This large number of European code users will directly profit from the extended capabilities of the tools and hereby they will be able to better predict safety margins of reactor cores. By decrease of the engineering factors this can lead to better fuel exploitation and improvement of the fuel cycle economy. Finally, due to the universality of the numerical tools based on Monte Carlo codes, the codes and methods developed in the frame of the HPMC project can easily be applied to research reactors, subcritical systems and Gen-IV reactors in addition to LWRs. For this reasons, a considerable economic impact can be expected from the HPMC project outcome.
Main Dissemination and training Activities:
The activities in this work package are aimed at bringing the project and its goals to the attention of possibly interested parties. On the other side, activities oriented to disseminate the knowledge generated within the project to “potential users” of the HPMC codes and methods but also to the academia were developed.
a) Various activities
A publicly accessible website is setup with information about the project, its aims, structure and participants. A separate website for internal use to the project was established to exchange information and documents not yet ready for publication. In 2013 the organisation of the websites is changed and brought under the auspices of KIT. The current public website address is www.fp7-hpmc.eu. It contains a private part under “members Area” to which all partner can access using a password. There all deliverables, presentations at the Executive Board meetings, presentations in conferences, journal and conference papers as well as the materials of the training courses are collected and will be kept there for longer time.
As the project goals fit well in the Strategic Research Agenda of the SNETP (Sustainable Nuclear Energy Technology Platform) a poster with the project setup and goals was presented at the 3rd General Assembly of SNETP in Warsaw, November 2011. Another general poster about the HPMC project was shown at the FISA conference in Vilnius, 2013. A seminar about the project was given at KTH in Stockholm in March 2012. A presentation about the project was given during a project meeting of the EU NURISP project in Karlsruhe, April 2012.
As the EU project NURISP and its successor NURESAFE are closely related to the goals of the HPMC project, good relationships were developed with the NURESASFE project and its management. Permission was asked and obtained to use the results for various benchmark test used in the NURISP project for comparison with results to be obtained with Monte Carlo calculation in the HPMC project. Also appointments were made with the NURESAFE project to study as far as possible the same benchmark test problems for the methods to be developed in both projects. The HPMC project is directed to the Monte Carlo method for the neutron transport; the NURESAFE project uses deterministic methods. The coordinator of the HPMC project was invited to give a lecture about the status of the codes and tools of the HPMC project in a NURESAFE ExCOM Meeting. This opportunity was used to explore the common areas of interest and to intensify the exchange between both projects which are very complementary (NURESAFE ExCOM meeting 7.10.2013 in Paris).
An extended summary of the project was prepared for publication by the European Commission. The deliverable is included in Sect. 5 on Management activities as well as the Communication Action Plan.
b) Training seminars
Training seminars are scheduled for potential users of the software developed in the project and are therefore outstanding opportunities to disseminate the knowledge and tools developed in the project.
A first training seminar was conducted at KIT in Germany on 24 and 25 April 2014 with 6 speakers directly involved in the project giving 11 different presentations covering most of the subjects developed during the project. The training seminar was attended by 9 external people from 6 different companies and institutes like Areva (Germany), CEA (France), GRS (Germany), HZDR (Germany), CCFE (UK), and from KIT.
The second training seminar was conducted at KTH in Stockholm, Sweden, on 4 and 5 June 2014 with basically the same program as for the first seminar. However, various changes were made to the presentations due to feedback and experience from the first seminar. This training seminar was attended by 13 external people from 5 different companies and institutes like Westinghouse, Vattenfall, University of Uppsala, Warsaw University of Technology, and from KTH.
Deliverable 5.12 gives an extensive overview of both training seminars and the attendees. Considering the goals set for successful training seminars with at least 4 participants from at least 2 different companies/institutes for both seminars, we can conclude that the training seminars were very successful.
Apart from the training seminars, scientific results of research performed for the project were published in several articles in international scientific journals and in papers at international conferences. Various other publications are submitted or in preparation.
List of Websites:
www.fp7-hpmc.eu
Main contacts:
Dr. Victor Hugo Sanchez Espinoza (Project Coordinator)
e-mail: victor.sanchez@kit.edu
Tel: +49-721-608 22283
Prof. Dr. J. Eduard Hoogenboom (Project Secretary)
e-mail: hoogenboom@delftnuclearconsultancy.nl
Tel: +31-10-4508442
