1"""
2This script builds a single PWR assembly and is a slightly more advanced
3demonstration of model building using Python. The creation of two universes for
4fuel pins and guide tube pins has been separated into functions, and then the
5overall model is built by an `assembly` function. This script also demonstrates
6the use of the `Model` class, which provides some extra convenience over using
7`Geometry`, `Materials`, and `Settings` classes directly. Finally, the script
8takes two command-line flags that indicate whether to build and/or run the
9model.
10
11"""
12
13import argparse
14from math import log10
15
16import numpy as np
17import openmc
18
19# Define surfaces
20fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
21clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
22
23# Define materials
24fuel = openmc.Material(name='Fuel')
25fuel.set_density('g/cm3', 10.29769)
26fuel.add_nuclide('U234', 4.4843e-6)
27fuel.add_nuclide('U235', 5.5815e-4)
28fuel.add_nuclide('U238', 2.2408e-2)
29fuel.add_nuclide('O16', 4.5829e-2)
30
31clad = openmc.Material(name='Cladding')
32clad.set_density('g/cm3', 6.55)
33clad.add_nuclide('Zr90', 2.1827e-2)
34clad.add_nuclide('Zr91', 4.7600e-3)
35clad.add_nuclide('Zr92', 7.2758e-3)
36clad.add_nuclide('Zr94', 7.3734e-3)
37clad.add_nuclide('Zr96', 1.1879e-3)
38
39hot_water = openmc.Material(name='Hot borated water')
40hot_water.set_density('g/cm3', 0.740582)
41hot_water.add_nuclide('H1', 4.9457e-2)
42hot_water.add_nuclide('O16', 2.4672e-2)
43hot_water.add_nuclide('B10', 8.0042e-6)
44hot_water.add_nuclide('B11', 3.2218e-5)
45hot_water.add_s_alpha_beta('c_H_in_H2O')
46
47
48def fuel_pin():
49    """Returns a fuel pin universe."""
50
51    fuel_cell = openmc.Cell(fill=fuel, region=-fuel_or)
52    clad_cell = openmc.Cell(fill=clad, region=+fuel_or & -clad_or)
53    hot_water_cell = openmc.Cell(fill=hot_water, region=+clad_or)
54
55    univ = openmc.Universe(name='Fuel Pin')
56    univ.add_cells([fuel_cell, clad_cell, hot_water_cell])
57    return univ
58
59
60def guide_tube_pin():
61    """Returns a control rod guide tube universe."""
62
63    gt_inner_cell = openmc.Cell(fill=hot_water, region=-fuel_or)
64    gt_clad_cell = openmc.Cell(fill=clad, region=+fuel_or & -clad_or)
65    gt_outer_cell = openmc.Cell(fill=hot_water, region=+clad_or)
66
67    univ = openmc.Universe(name='Guide Tube')
68    univ.add_cells([gt_inner_cell, gt_clad_cell, gt_outer_cell])
69    return univ
70
71
72def assembly_model():
73    """Returns a single PWR fuel assembly."""
74
75    model = openmc.model.Model()
76
77    # Create fuel assembly Lattice
78    pitch = 21.42
79    assembly = openmc.RectLattice(name='Fuel Assembly')
80    assembly.pitch = (pitch/17, pitch/17)
81    assembly.lower_left = (-pitch/2, -pitch/2)
82
83    # Create array indices for guide tube locations in lattice
84    gt_pos = np.array([
85                  [2, 5],   [2, 8],   [2, 11],
86             [3, 3],                        [3, 13],
87        [5, 2],   [5, 5],   [5, 8],   [5, 11],    [5, 14],
88        [8, 2],   [8, 5],   [8, 8],   [8, 11],    [8, 14],
89        [11, 2],  [11, 5],  [11, 8],  [11, 11],   [11, 14],
90             [13, 3],                       [13, 13],
91                  [14, 5],  [14, 8],  [14, 11]
92    ])
93
94    # Create 17x17 array of universes. First we create a 17x17 array all filled
95    # with the fuel pin universe. Then, we replace the guide tube positions with
96    # the guide tube pin universe (note the use of numpy fancy indexing to
97    # achieve this).
98    assembly.universes = np.full((17, 17), fuel_pin())
99    assembly.universes[gt_pos[:, 0], gt_pos[:, 1]] = guide_tube_pin()
100
101    # Create outer boundary of the geometry to surround the lattice
102    outer_boundary = openmc.model.rectangular_prism(
103        pitch, pitch, boundary_type='reflective')
104
105    # Create a cell filled with the lattice
106    main_cell = openmc.Cell(fill=assembly, region=outer_boundary)
107
108    # Finally, create geometry by providing a list of cells that fill the root
109    # universe
110    model.geometry = openmc.Geometry([main_cell])
111
112    model.settings.batches = 150
113    model.settings.inactive = 50
114    model.settings.particles = 1000
115    model.settings.source = openmc.Source(space=openmc.stats.Box(
116        (-pitch/2, -pitch/2, -1),
117        (pitch/2, pitch/2, 1),
118        only_fissionable=True
119    ))
120
121    # NOTE: We never actually created a Materials object. When you export/run
122    # using the Model object, if no materials were assigned it will look through
123    # the Geometry object and automatically export any materials that are
124    # necessary to build the model.
125    return model
126
127
128if __name__ == '__main__':
129    # Set up command-line arguments for generating/running the model
130    parser = argparse.ArgumentParser()
131    parser.add_argument('--generate', action='store_true')
132    parser.add_argument('--run', action='store_true')
133    args = parser.parse_args()
134
135    if args.generate or args.run:
136        model = assembly_model()
137        if args.generate:
138            model.export_to_xml()
139        if args.run:
140            model.run()
141