Nigel Davies

1. Monte Carlo Algorithm for neutronics calculations with Thermal hydraulics feedback on parallel computers Nigel Davies

2. PhD TITLE .... Whole Core PWR core-follow modelling: Coupled 3D Monte Carlo time dependent neutronics and 3D Thermal Hydraulics calculations, with comparative runtimes to coupled Deterministic methods

3. PhD In a nutshell .... Calculate k-effective, flux (in 172 groups), and power distributions for an explicitly modelled (in 3D) whole core PWR taking account of depletion and thermal hydraulics feedback – and find a way of making the calculation run in a reasonable time.

4. Presentation Outline • This presentation describes a plan for parallelising an existing 3D neutronics reactor physics and criticality code MONK. • MONK is sold commercially by the ANSWERS business – part of SERCO • Google: “ANSWERS” “SERCO” for more details • Or go to http://www.sercoassurance.com/answers/

5. Presentation Outline It covers … • A very basic description of the MONK TH-code coupling scheme • Brief review of how MONK works • A flowchart description of parallelising MONK

6. MONK-TH code coupling The basic idea is ……….. • MONK is used to calculate the power distribution at each time-step taking account of: • depletion (i.e. a burn-up calculation) • temperature distribution (temperatures from user input or TH code) • At the end of each time-step the calculated power distribution is fed into a Thermal Hydraulics calculation which in turn generates a new temperature distribution. • The new temperatures are applied to the materials in the MONK calculation at the start of the next time-step • The loop continues until all time-steps have been done.

7. MONK-TH code coupling overall scheme (a picture!) Time consuming bit is MONK – hence the need to consider the use of Parallel Computing to speed it up! ........ This is the focus of this talk

8. Capabilities of MONK • MONK is capable of modelling very complex geometries (see e.g. model of BWR). • Used widely both in UK and abroad for criticality and some reactor physics studies • Iterative solution to the Boltzmann Transport Equation (each iteration called a “stage”) 5 Initial Source Guess Solution (k-eff, flux, etc) N Stages New Source Guess

9. Tracking Techniques – WOODCOCK TRACKING • Invent a region called a HOLE MATERIAL • Make every material in this region have the SAME total cross section • To do this we invent (for every material in the HOLE) a new reaction called a “pseudo collision” (when this occurs – nothing happens!) • Cross section for a “pseudo collision” is added to the existing total cross section such that all materials have a new but equal total cross section • Well used technique in the Monte Carlo world. • Woodcock, E., Murphy, T., Hemmings, P., and Longworth, T. 1965. Techniques used in the GEM code for Monte Carlo neutronics calculations in reactors and other systems of complex geometry. In Proc. Conference on the Application of Computing Methods to Reactor Problems, ANL-7050, 557--579.

10. Tracking Techniques – WEIGHTED PARTICLES • Particles are tracked and each one has an associated weight (start weight =1) • At a “real” collision the particle weight is reduced • Absorption is modelled by reducing the weight by a proportion equal to the absorption probability (same as multiplying by scatter probability) • “Russian Roulette” technique used to determine when a particle “dies” (when below a cutoff) • Survivors continue with weight=1 W W ‘ W ‘ = W x { Σ(scat) / Σ(total) }

11. Monte Carlo Process in MONK • We have a complex geometry within the boundaries of a box • Use “Woodcock Tracking” through the ENTIRE geometry • Choose x,y,z,g,u,v,w + Wt. • Obtain mfp = 1/Σ(max) • Calculate collision point coordinates • Pseudo collision • Real collision – Wt is reduced • Roulette the sample? • Sample leaks out 1 7 2 6 5 3 4

12. Calc. SOURCE particle data TRACK to collision pt. (use Σmax) Find MATERIAL @ collision pt. REACTION Sample dead? MONK – Process flow diagram (1 of2) Let’s convert this process into a process flow diagram…. Note: on next slide I will call this part of the flowchart : SOURCE / TRACK / MATERIAL / REACTION(start) Pseudo Real Yes No Next Page

13. Accumulate Reaction Rates Calc. new g, u, v, w & Wt. Stage Over ? Finish ? MONK – Process flow diagram (2of 2) SOURCE / TRACK / MATERIAL /REACTION(start) Fill the birth store (next stage source) Next stage source ACCUMULATOR No Get next source particle Yes Calculate K-effective & Print stage result Calculate & Print Results STOP Yes Tell next source to start sampling No

14. MONK Process Flow Diagram -summarised InitialSource Birth Store Calculate Track Length Fluxes Tracking Material Search Σ(max) Reaction Processing Accumulate & Print Reaction rates Pseudo & scatter Birth Store data (Fission & Children)

15. MONK Process Flow Diagram converted to a plan S1 S2 T-1 … T-n F M-1 … M-n Σ(max) To Relevant Reaction Computer R-1 … R-n Σ(max) A Reaction rates Pseudo & scatter Birth Store data (Fission & Children)

16. Summary We have covered … • A very basic description of the MONK TH-code coupling scheme • Brief review of how MONK works • A flowchart description of parallelising MONK

17. APPENDICES The following are additional slides used in case of questions …

18. Overall method In summary : I hope I have given you a taste of how I propose to parallelise MONK. • It offers the flexibility of either the user defining how many T, M or R computers there are, or ultimately an optimisation algorithm could decide. • By using many R computers we can expand the detail that can be modelled (Domain Decomposition) • Similarly by using parallelism we can speed up the process

19. Function of each computer • S1 = Generates the source particle data in the first stage and every alternate subsequent stage. In the second stage (and every alternate stage subsequent to the second stage) it acts as the birth store. • S2 = similar to S1 but only starts generating source in the second stage. This is because it acts as the birth store in the first stage. • T-1:T-n = These are the TRACKING computers. They receive data from S1/S2, calculate collision point data and passes this on to the next available MATERIAL SEARCH computer. (These can also pass track length data to the FLUX computer) • M-1:M1n = These are the MATERIAL SEARCH computers. They search for the material number at the collision point and pass the data on to the relevant REACTION computer.

20. Function of each computer • R-1:R-n = These are the REACTION computers. Each one processes collisions only in a given range of Artificial Materials. This enables them to hold only a limited amount of the nuclear data (which is the main source of the large RAM requirement). • They also send pseudo collision and real scattered particles back to the TRACKING computers • They also send the reaction rates to the ACCUMULATOR computer • They also send birth store data to the relevant SOURCE computer S1 or S2, as well as information to S1/S2 when a particle dies (through roulette) or leakage – this way both the ACCUMULATOR and S1/S2 know when a stage completes and subsequently when the cycle completes. • A and F = These are the ACCUMULATOR and FLUX computers respectively. These process reaction rates and fluxes and control when results are printed.

21. Function of each computer • We have only looked at processing and output here but the other aspect is input and pre-processing • However using MPI we can differentiate between computers (by rank) and thus extract only what input is required. • For example only the REACTION computers need the nuclear data -hence the need for them to pass Σ(max) values to the ACCUMULATOR which processes them to find the true maximum in each group. This data is then passed to the TRACKING computers prior to them tracking (other than that the TRACKING computers require none of the nuclear data and none of the geometry data!)

22. MONK-TH code coupling: BU Mesh and TH Mesh • Enables the user to superimpose a X-Y-Z mesh over the MONK geometry (defines the user-required mesh for a burn-up calculation - called the BU-mesh) • In addition the user superimposes another X-Y-Z mesh to define the mesh for calculating the power distribution (for the TH code) and for receiving new temperatures and densities (from the TH code) • Such a X-Y-Z mesh (called the TH-mesh) is easier to translate to and from the TH code mesh than MONK geometry regions are

23. BUTH Mesh and Artificial Materials • X-Y-Z meshing is defined by superimposed a BOX body over the MONK geometry. • The 2 meshes (BU-mesh and TH-mesh) are combined by the code into a so-called BUTH-mesh. • Reaction rates are calculated for each material inside a BUTH mesh element. • Each material inside a BUTH mesh element is reassigned to a so-called “Artificial Material”

24. “Artificial Materials” • How does the code find out what materials there are inside a BUTH mesh element ? • Before any calculation is done a set of tracks are passed through each BUTH mesh element • Each track stops at regular locations along its length • The code determines (just as it does during its normal tracking) what material is at each position • An “Artificial Material number” is assigned to each “User defined” material number found in there.

25. WOODCOCK TRACKING • The following slides give a brief overview of Woodcock tracking …

26. Hole Material Tracking • Tracking body boundaries is necessary only because the mean free path varies between materials. • If the mean free path was constant within a multi-material region, the tracking of boundaries would not be required. • Since the mean free path is not naturally constant between materials, artificial extra cross-sections would have to be added if we want the total cross-section, and hence the mean free path, to be constant for all materials in a region.

27. Hole Material Tracking • Consider a region containing materials M1, M2, .... Mn with macroscopic total cross-sections 1,  2, ....  n. • In order to achieve a constant mean free path within the region the largest cross-section is located ( max) and artificial cross-sections are added to the other values to make them equal to  max: i.e.  max = max( 1,  2, ....  n)  max =  1 +  1’ =  2 +  2’ = .... =  n +  n’ where  1’,  2’, ....  n’ are additional artificial cross-sections that must be dealt with.

28. Hole Material Tracking • Tracking within the hole material now proceeds using the minimum mean free path (1/  max). • At each collision site the material present is identified. • If the material present is the one with the largest cross-section then a normal (real) collision occurs. • If the material is one of the others (say material i), then there is a probability of  i/  max that the collision is a real collision. Otherwise the collision is due to the artificial cross-section  i’.

29. Hole Material Tracking • The device of increasing the total cross-sections and introducing additional artificial cross-sections therefore produces additional artificial collisions. These must clearly not disturb the passage of the particle, so that the distribution of real particle collisions is not disturbed in any way. B 3 2 4 1 A 5 C Artificial collisions (1- 5) do not disturb the passage of the particle in its direction of travel between real collisions (A - C)

30. Hole Material Tracking • It can be shown that although the total number of collisions is increased by the addition of the artificial cross-sections, the distribution of real collisions is unchanged, and hence the process is unbiased. • That it is correct to do this may be seen from the following consideration. The equation from which the distribution of collision points is derived is: I(x+dx) = I(x) -  i I(x) dx where  i has the meaning above, and the function I(x) is the intensity of neutrons at a point x along an axis. The second term on the right hand side relates to the removal of neutrons due to collisions.

31. Hole Material Tracking • If one now uses a cross-section , a term must be created relating to the production of neutrons. The equation becomes: I(x+dx) = I(x) -  I(x) dx + ( -  i) I(x) dx • The two equations are therefore equivalent and hence the distribution of real collisions will be given correctly