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:

  • Multi-group materials loaded from an external HDF5 library via mcdc.MaterialMG(library=...).

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

How to Run

From the repository root run:

python examples/c5g7/k-eigenvalue/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.