C5G7 — k-eigenvalue example#

Description#

Multigroup k-eigenvalue calculation for the C5G7 benchmark using the packaged MGXS HDF5 library in examples/c5g7. This example performs a static criticality calculation and reports \(k_{\mathrm{eff}}\) and gyration-radius diagnostics.

Step-by-Step Walkthrough#

The C5G7 benchmark uses a pre-packaged 7-group cross-section library (MGXS-C5G7-TD.h5). The input file defines the full-core geometry using MC/DC’s lattice and universe system.

Key concepts demonstrated:

  • Multigroup data loaded from an external HDF5 library and used to construct materials through mcdc.Material.multigroup(...).

  • Pin-cell universes built from cylindrical fuel pins in square moderator cells.

  • Lattice assemblies that tile pin-cell universes into fuel assemblies of different enrichments.

  • Core lattice that arranges assemblies and reflector regions.

  • k-eigenvalue mode with set_eigenmode() for criticality.

Refer to the embedded code below — comments in the source mark each section (materials, pins, assemblies, core, source, tallies, settings).

What to try:

  • Increase N_particle for better \(k_{\text{eff}}\) statistics.

  • Adjust the number of inactive/active cycles.

  • Compare \(k_{\text{eff}}\) with the published C5G7 reference value.

Full Input#

Click here to view the input file: examples/c5g7/k-eigenvalue/input.py.

The complete input used for this example is embedded below:

  1import h5py
  2import numpy as np
  3
  4import mcdc
  5
  6simulation = mcdc.Simulation("C5G7 k-eigenvalue")
  7
  8# =============================================================================
  9# Materials
 10# =============================================================================
 11
 12# Load material data
 13lib = h5py.File("../MGXS-C5G7-TD.h5", "r")
 14
 15
 16# Setter
 17def set_mat(mat):
 18    return mcdc.Material.multigroup(
 19        capture=mat["capture"][:],
 20        scatter=mat["scatter"][:],
 21        fission=mat["fission"][:],
 22        nu_p=mat["nu_p"][:],
 23        nu_d=mat["nu_d"][:],
 24        chi_p=mat["chi_p"][:],
 25        chi_d=mat["chi_d"][:],
 26        speed=mat["speed"][:],
 27        decay_rate=mat["decay"][:],
 28    )
 29
 30
 31# Materials
 32mat_uo2 = set_mat(lib["uo2"])  # Fuel: UO2
 33mat_mox43 = set_mat(lib["mox43"])  # Fuel: MOX 4.3%
 34mat_mox7 = set_mat(lib["mox7"])  # Fuel: MOX 7.0%
 35mat_mox87 = set_mat(lib["mox87"])  # Fuel: MOX 8.7%
 36mat_gt = set_mat(lib["gt"])  # Guide tube
 37mat_fc = set_mat(lib["fc"])  # Fission chamber
 38mat_cr = set_mat(lib["cr"])  # Control rod
 39mat_mod = set_mat(lib["mod"])  # Moderator
 40
 41# =============================================================================
 42# Pin cells
 43# =============================================================================
 44
 45pitch = 1.26
 46radius = 0.54
 47core_height = 128.52
 48refl_thick = 21.42
 49
 50# Control rod banks fractions
 51#   All out: 0.0
 52#   All in : 1.0
 53cr1 = 0.0
 54cr2 = 0.0
 55cr3 = 0.0
 56cr4 = 0.0
 57# Control rod banks interfaces
 58cr1 = core_height * (0.5 - cr1)
 59cr2 = core_height * (0.5 - cr2)
 60cr3 = core_height * (0.5 - cr3)
 61cr4 = core_height * (0.5 - cr4)
 62
 63# Surfaces
 64cy = mcdc.Surface.CylinderZ(center=[0.0, 0.0], radius=radius)
 65z1 = mcdc.Surface.PlaneZ(z=cr1)  # Control rod banks interfaces
 66z2 = mcdc.Surface.PlaneZ(z=cr2)
 67z3 = mcdc.Surface.PlaneZ(z=cr3)
 68z4 = mcdc.Surface.PlaneZ(z=cr4)
 69zf = mcdc.Surface.PlaneZ(z=core_height / 2)
 70
 71# Fission chamber
 72fc = mcdc.Cell(-cy, mat_fc)
 73mod = mcdc.Cell(+cy, mat_mod)
 74fission_chamber = mcdc.Universe(cells=[fc, mod])
 75
 76# Fuel rods
 77uo2 = mcdc.Cell(-cy & -zf, mat_uo2)
 78mox4 = mcdc.Cell(-cy & -zf, mat_mox43)
 79mox7 = mcdc.Cell(-cy & -zf, mat_mox7)
 80mox8 = mcdc.Cell(-cy & -zf, mat_mox87)
 81moda = mcdc.Cell(-cy & +zf, mat_mod)  # Water above pin
 82fuel_uo2 = mcdc.Universe(cells=[uo2, mod, moda])
 83fuel_mox43 = mcdc.Universe(cells=[mox4, mod, moda])
 84fuel_mox7 = mcdc.Universe(cells=[mox7, mod, moda])
 85fuel_mox87 = mcdc.Universe(cells=[mox8, mod, moda])
 86
 87# Control rods and guide tubes
 88cr1 = mcdc.Cell(-cy & +z1, mat_cr)
 89cr2 = mcdc.Cell(-cy & +z2, mat_cr)
 90cr3 = mcdc.Cell(-cy & +z3, mat_cr)
 91cr4 = mcdc.Cell(-cy & +z4, mat_cr)
 92gt1 = mcdc.Cell(-cy & -z1, mat_gt)
 93gt2 = mcdc.Cell(-cy & -z2, mat_gt)
 94gt3 = mcdc.Cell(-cy & -z3, mat_gt)
 95gt4 = mcdc.Cell(-cy & -z4, mat_gt)
 96control_rod1 = mcdc.Universe(cells=[cr1, gt1, mod])
 97control_rod2 = mcdc.Universe(cells=[cr2, gt2, mod])
 98control_rod3 = mcdc.Universe(cells=[cr3, gt3, mod])
 99control_rod4 = mcdc.Universe(cells=[cr4, gt4, mod])
100
101# =============================================================================
102# Fuel lattices
103# =============================================================================
104
105# UO2 lattice 1
106u = fuel_uo2
107c = control_rod1
108f = fission_chamber
109lattice_1 = mcdc.Lattice(
110    x=[-pitch * 17 / 2, pitch, 17],
111    y=[-pitch * 17 / 2, pitch, 17],
112    universes=[
113        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
114        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
115        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
116        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
117        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
118        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
119        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
120        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
121        [u, u, c, u, u, c, u, u, f, u, u, c, u, u, c, u, u],
122        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
123        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
124        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
125        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
126        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
127        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
128        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
129        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
130    ],
131)
132
133# MOX lattice 2
134l = fuel_mox43
135m = fuel_mox7
136n = fuel_mox87
137c = control_rod2
138f = fission_chamber
139lattice_2 = mcdc.Lattice(
140    x=[-pitch * 17 / 2, pitch, 17],
141    y=[-pitch * 17 / 2, pitch, 17],
142    universes=[
143        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
144        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
145        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
146        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
147        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
148        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
149        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
150        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
151        [l, m, c, n, n, c, n, n, f, n, n, c, n, n, c, m, l],
152        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
153        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
154        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
155        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
156        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
157        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
158        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
159        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
160    ],
161)
162
163# MOX lattice 3
164l = fuel_mox43
165m = fuel_mox7
166n = fuel_mox87
167c = control_rod3
168f = fission_chamber
169lattice_3 = mcdc.Lattice(
170    x=[-pitch * 17 / 2, pitch, 17],
171    y=[-pitch * 17 / 2, pitch, 17],
172    universes=[
173        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
174        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
175        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
176        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
177        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
178        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
179        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
180        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
181        [l, m, c, n, n, c, n, n, f, n, n, c, n, n, c, m, l],
182        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
183        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
184        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
185        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
186        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
187        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
188        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
189        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
190    ],
191)
192
193# UO2 lattice 4
194u = fuel_uo2
195c = control_rod4
196f = fission_chamber
197lattice_4 = mcdc.Lattice(
198    x=[-pitch * 17 / 2, pitch, 17],
199    y=[-pitch * 17 / 2, pitch, 17],
200    universes=[
201        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
202        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
203        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
204        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
205        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
206        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
207        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
208        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
209        [u, u, c, u, u, c, u, u, f, u, u, c, u, u, c, u, u],
210        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
211        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
212        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
213        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
214        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
215        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
216        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
217        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
218    ],
219)
220
221# =============================================================================
222# Assemblies and core
223# =============================================================================
224
225# Surfaces
226x0 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="reflective")
227x1 = mcdc.Surface.PlaneX(x=pitch * 17)
228x2 = mcdc.Surface.PlaneX(x=pitch * 17 * 2)
229x3 = mcdc.Surface.PlaneX(x=pitch * 17 * 3, boundary_condition="vacuum")
230
231y0 = mcdc.Surface.PlaneY(y=-pitch * 17 * 3, boundary_condition="vacuum")
232y1 = mcdc.Surface.PlaneY(y=-pitch * 17 * 2)
233y2 = mcdc.Surface.PlaneY(y=-pitch * 17)
234y3 = mcdc.Surface.PlaneY(y=0.0, boundary_condition="reflective")
235
236z0 = mcdc.Surface.PlaneZ(z=-(core_height / 2 + refl_thick), boundary_condition="vacuum")
237z1 = mcdc.Surface.PlaneZ(z=-(core_height / 2))
238z2 = mcdc.Surface.PlaneZ(z=(core_height / 2 + refl_thick), boundary_condition="vacuum")
239
240# Assembly cells
241center = np.array([pitch * 17 / 2, -pitch * 17 / 2, 0.0])
242assembly_1 = mcdc.Cell(+x0 & -x1 & +y2 & -y3 & +z1 & -z2, lattice_1, translation=center)
243
244center += np.array([pitch * 17, 0.0, 0.0])
245assembly_2 = mcdc.Cell(+x1 & -x2 & +y2 & -y3 & +z1 & -z2, lattice_2, translation=center)
246
247center += np.array([-pitch * 17, -pitch * 17, 0.0])
248assembly_3 = mcdc.Cell(+x0 & -x1 & +y1 & -y2 & +z1 & -z2, lattice_3, translation=center)
249
250center += np.array([pitch * 17, 0.0, 0.0])
251assembly_4 = mcdc.Cell(+x1 & -x2 & +y1 & -y2 & +z1 & -z2, lattice_4, translation=center)
252
253# Bottom reflector cell
254reflector_bottom = mcdc.Cell(+x0 & -x3 & +y0 & -y3 & +z0 & -z1, mat_mod)
255
256# Side reflectors
257reflector_south = mcdc.Cell(+x0 & -x3 & +y0 & -y1 & +z1 & -z2, mat_mod)
258reflector_east = mcdc.Cell(+x2 & -x3 & +y1 & -y3 & +z1 & -z2, mat_mod)
259
260# Set model
261simulation.set_model(
262    [
263        assembly_1,
264        assembly_2,
265        assembly_3,
266        assembly_4,
267        reflector_bottom,
268        reflector_south,
269        reflector_east,
270    ]
271)
272
273# =============================================================================
274# Set source
275# =============================================================================
276
277source = mcdc.Source(
278    x=[0.0, pitch * 17 * 2],
279    y=[-pitch * 17 * 2, 0.0],
280    z=[-core_height / 2, core_height / 2],
281    isotropic=True,
282    energy=0,  # Highest energy
283)
284simulation.set_sources([source])
285
286# =============================================================================
287# Set tallies, settings, techniques and run MC/DC
288# =============================================================================
289
290# Tally
291x_grid = np.linspace(0.0, pitch * 17 * 3, 17 * 3 + 1)
292y_grid = np.linspace(-pitch * 17 * 3, 0.0, 17 * 3 + 1)
293z_grid = np.linspace(
294    -(core_height / 2 + refl_thick), (core_height / 2 + refl_thick), 102 + 17 * 2 + 1
295)
296g_grid = np.array([-0.5, 3.5, 6.5])  # Collapsing to fast (1-4) and slow (5-7)
297mesh = mcdc.MeshStructured(x=x_grid, y=y_grid, z=z_grid)
298tally = mcdc.Tally(mesh=mesh, scores=["flux"], energy=g_grid)
299simulation.set_tallies([tally])
300
301# Settings
302simulation.settings.N_particle = 50
303simulation.settings.census_bank_buffer_ratio = 4.0
304simulation.settings.set_eigenmode(N_inactive=5, N_active=10, gyration_radius="all")
305
306# Techniques
307simulation.technique.population_control()
308
309# Run
310simulation.run()

How to Run#

From inside examples/c5g7/k-eigenvalue run:

python input.py

Expected Output#

Eigenvalue history printed to stdout and HDF5 tally data for flux and gyration-radius diagnostics saved in the build artifacts folder.