C5G7 — Transient example

Description

Time-dependent C5G7-TD transient driven by control-rod movements and a time-limited source. Uses the packaged MGXS library in examples/c5g7 and demonstrates moving surfaces and time-resolved tallies.

Step-by-Step Walkthrough

This example extends the C5G7 k-eigenvalue setup with time-dependent features:

  • Moving surfaces simulate control-rod insertion/withdrawal.

  • Time-resolved tallies capture the transient fission rate.

  • Time census checkpoints the particle population at specified intervals for population control.

The geometry and material setup is identical to the k-eigenvalue case. The transient-specific additions are:

  1. Surface velocities assigned via surface.move(...).

  2. A time grid added to the mesh tally.

  3. set_time_census(...) for time-step population control.

Refer to the embedded code below for the full implementation.

What to try:

  • Change the rod insertion speed to see prompt vs. delayed transient response.

  • Add more time census points for finer population control.

  • Compare power history with published C5G7-TD benchmarks.

Full Input

Click here to view the input file: examples/c5g7/transient/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
 46reflector_thickness = 21.42
 47
 48# Control rod banks fractions
 49#   All out: 0.0
 50#   All in : 1.0
 51cr1 = np.array([1.0, 1.0, 0.89, 1.0])
 52cr1_t = np.array([0.0, 10.0, 15.0, 15.0 + 1.0 - cr1[-2]])
 53
 54cr2 = np.array([1.0, 1.0, 0.0, 0.0, 0.8])
 55cr2_t = np.array([0.0, 5.0, 10.0, 15.0, 15.8])
 56
 57cr3 = np.array([0.75, 0.75, 1.0])
 58cr3_t = np.array([0.0, 15.0, 15.25])
 59
 60cr4 = np.array([1.0, 1.0, 0.5, 0.5, 1.0])
 61cr4_t = np.array(
 62    [0.0, 5.0, 5.0 + (cr4[1] - cr4[2]) / 2 * 10, 15.0, 15.0 + 1.0 - cr4[-2]]
 63)
 64
 65# Tips of the control rod banks
 66cr1_bottom = core_height * (0.5 - cr1)
 67cr2_bottom = core_height * (0.5 - cr2)
 68cr3_bottom = core_height * (0.5 - cr3)
 69cr4_bottom = core_height * (0.5 - cr4)
 70cr1_top = cr1_bottom + core_height
 71cr2_top = cr2_bottom + core_height
 72cr3_top = cr3_bottom + core_height
 73cr4_top = cr4_bottom + core_height
 74
 75# Durations of the moving tips
 76cr1_durations = cr1_t[1:] - cr1_t[:-1]
 77cr2_durations = cr2_t[1:] - cr2_t[:-1]
 78cr3_durations = cr3_t[1:] - cr3_t[:-1]
 79cr4_durations = cr4_t[1:] - cr4_t[:-1]
 80
 81# Velocities of the moving tips
 82cr1_velocities = np.zeros((len(cr1) - 1, 3))
 83cr2_velocities = np.zeros((len(cr2) - 1, 3))
 84cr3_velocities = np.zeros((len(cr3) - 1, 3))
 85cr4_velocities = np.zeros((len(cr4) - 1, 3))
 86cr1_velocities[:, 2] = (cr1_top[1:] - cr1_top[:-1]) / cr1_durations
 87cr2_velocities[:, 2] = (cr2_top[1:] - cr2_top[:-1]) / cr2_durations
 88cr3_velocities[:, 2] = (cr3_top[1:] - cr3_top[:-1]) / cr3_durations
 89cr4_velocities[:, 2] = (cr4_top[1:] - cr4_top[:-1]) / cr4_durations
 90
 91# Surfaces
 92cy = mcdc.Surface.CylinderZ(center=[0.0, 0.0], radius=radius)
 93# Control rod top and bottom tips
 94z1_top = mcdc.Surface.PlaneZ(z=cr1_top[0])
 95z1_bottom = mcdc.Surface.PlaneZ(z=cr1_bottom[0])
 96z2_top = mcdc.Surface.PlaneZ(z=cr2_top[0])
 97z2_bottom = mcdc.Surface.PlaneZ(z=cr2_bottom[0])
 98z3_top = mcdc.Surface.PlaneZ(z=cr3_top[0])
 99z3_bottom = mcdc.Surface.PlaneZ(z=cr3_bottom[0])
100z4_top = mcdc.Surface.PlaneZ(z=cr4_top[0])
101z4_bottom = mcdc.Surface.PlaneZ(z=cr4_bottom[0])
102# Fuel top
103#   (Bottom is bounded by the universe cell)
104zf = mcdc.Surface.PlaneZ(z=0.5 * core_height)
105
106# Move the control tips
107z1_top.move(cr1_velocities, cr1_durations)
108z1_bottom.move(cr1_velocities, cr1_durations)
109z2_top.move(cr2_velocities, cr2_durations)
110z2_bottom.move(cr2_velocities, cr2_durations)
111z3_top.move(cr3_velocities, cr3_durations)
112z3_bottom.move(cr3_velocities, cr3_durations)
113z4_top.move(cr4_velocities, cr4_durations)
114z4_bottom.move(cr4_velocities, cr4_durations)
115
116# Fission chamber pin
117fc = mcdc.Cell(-cy, mat_fc)
118mod = mcdc.Cell(+cy, mat_mod)
119fission_chamber = mcdc.Universe(cells=[fc, mod])
120
121# Fuel rods
122uo2 = mcdc.Cell(-cy & -zf, mat_uo2)
123mox4 = mcdc.Cell(-cy & -zf, mat_mox43)
124mox7 = mcdc.Cell(-cy & -zf, mat_mox7)
125mox8 = mcdc.Cell(-cy & -zf, mat_mox87)
126moda = mcdc.Cell(-cy & +zf, mat_mod)  # Water above pin
127fuel_uo2 = mcdc.Universe(cells=[uo2, mod, moda])
128fuel_mox43 = mcdc.Universe(cells=[mox4, mod, moda])
129fuel_mox7 = mcdc.Universe(cells=[mox7, mod, moda])
130fuel_mox87 = mcdc.Universe(cells=[mox8, mod, moda])
131
132# Control rods and guide tubes
133cr1 = mcdc.Cell(-cy & +z1_bottom & -z1_top, mat_cr)
134gt1_lower = mcdc.Cell(-cy & -z1_bottom, mat_gt)
135gt1_upper = mcdc.Cell(-cy & +z1_top, mat_gt)
136#
137cr2 = mcdc.Cell(-cy & +z2_bottom & -z2_top, mat_cr)
138gt2_lower = mcdc.Cell(-cy & -z2_bottom, mat_gt)
139gt2_upper = mcdc.Cell(-cy & +z2_top, mat_gt)
140#
141cr3 = mcdc.Cell(-cy & +z3_bottom & -z3_top, mat_cr)
142gt3_lower = mcdc.Cell(-cy & -z3_bottom, mat_gt)
143gt3_upper = mcdc.Cell(-cy & +z3_top, mat_gt)
144#
145cr4 = mcdc.Cell(-cy & +z4_bottom & -z4_top, mat_cr)
146gt4_lower = mcdc.Cell(-cy & -z4_bottom, mat_gt)
147gt4_upper = mcdc.Cell(-cy & +z4_top, mat_gt)
148#
149control_rod1 = mcdc.Universe(cells=[cr1, gt1_lower, gt1_upper, mod])
150control_rod2 = mcdc.Universe(cells=[cr2, gt2_lower, gt2_upper, mod])
151control_rod3 = mcdc.Universe(cells=[cr3, gt3_lower, gt3_upper, mod])
152control_rod4 = mcdc.Universe(cells=[cr4, gt4_lower, gt4_upper, mod])
153
154# =============================================================================
155# Fuel lattices
156# =============================================================================
157
158# UO2 lattice 1
159u = fuel_uo2
160c = control_rod1
161f = fission_chamber
162lattice_1 = mcdc.Lattice(
163    x=[-pitch * 17 / 2, pitch, 17],
164    y=[-pitch * 17 / 2, pitch, 17],
165    universes=[
166        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
167        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
168        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
169        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
170        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
171        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
172        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
173        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
174        [u, u, c, u, u, c, u, u, f, u, u, c, u, u, c, u, u],
175        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
176        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
177        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
178        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
179        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
180        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
181        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
182        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
183    ],
184)
185
186# MOX lattice 2
187l = fuel_mox43
188m = fuel_mox7
189n = fuel_mox87
190c = control_rod2
191f = fission_chamber
192lattice_2 = mcdc.Lattice(
193    x=[-pitch * 17 / 2, pitch, 17],
194    y=[-pitch * 17 / 2, pitch, 17],
195    universes=[
196        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
197        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
198        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
199        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
200        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
201        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
202        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
203        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
204        [l, m, c, n, n, c, n, n, f, n, n, c, n, n, c, m, l],
205        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
206        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
207        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
208        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
209        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
210        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
211        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
212        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
213    ],
214)
215
216# MOX lattice 3
217l = fuel_mox43
218m = fuel_mox7
219n = fuel_mox87
220c = control_rod3
221f = fission_chamber
222lattice_3 = mcdc.Lattice(
223    x=[-pitch * 17 / 2, pitch, 17],
224    y=[-pitch * 17 / 2, pitch, 17],
225    universes=[
226        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
227        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
228        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
229        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
230        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
231        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
232        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
233        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
234        [l, m, c, n, n, c, n, n, f, n, n, c, n, n, c, m, l],
235        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
236        [l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
237        [l, m, c, n, n, c, n, n, c, n, n, c, n, n, c, m, l],
238        [l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
239        [l, m, m, c, m, n, n, n, n, n, n, n, m, c, m, m, l],
240        [l, m, m, m, m, c, m, m, c, m, m, c, m, m, m, m, l],
241        [l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
242        [l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
243    ],
244)
245
246# UO2 lattice 4
247u = fuel_uo2
248c = control_rod4
249f = fission_chamber
250lattice_4 = mcdc.Lattice(
251    x=[-pitch * 17 / 2, pitch, 17],
252    y=[-pitch * 17 / 2, pitch, 17],
253    universes=[
254        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
255        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
256        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
257        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
258        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
259        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
260        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
261        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
262        [u, u, c, u, u, c, u, u, f, u, u, c, u, u, c, u, u],
263        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
264        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
265        [u, u, c, u, u, c, u, u, c, u, u, c, u, u, c, u, u],
266        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
267        [u, u, u, c, u, u, u, u, u, u, u, u, u, c, u, u, u],
268        [u, u, u, u, u, c, u, u, c, u, u, c, u, u, u, u, u],
269        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
270        [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
271    ],
272)
273
274# =============================================================================
275# Assemblies and core
276# =============================================================================
277
278# Surfaces
279x0 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="reflective")
280x1 = mcdc.Surface.PlaneX(x=pitch * 17)
281x2 = mcdc.Surface.PlaneX(x=pitch * 17 * 2)
282x3 = mcdc.Surface.PlaneX(x=pitch * 17 * 3, boundary_condition="vacuum")
283
284y0 = mcdc.Surface.PlaneY(y=-pitch * 17 * 3, boundary_condition="vacuum")
285y1 = mcdc.Surface.PlaneY(y=-pitch * 17 * 2)
286y2 = mcdc.Surface.PlaneY(y=-pitch * 17)
287y3 = mcdc.Surface.PlaneY(y=0.0, boundary_condition="reflective")
288
289z0 = mcdc.Surface.PlaneZ(
290    z=-(core_height / 2 + reflector_thickness), boundary_condition="vacuum"
291)
292z1 = mcdc.Surface.PlaneZ(z=-(core_height / 2))
293z2 = mcdc.Surface.PlaneZ(
294    z=(core_height / 2 + reflector_thickness), boundary_condition="vacuum"
295)
296
297# Assembly cells
298center = np.array([pitch * 17 / 2, -pitch * 17 / 2, 0.0])
299assembly_1 = mcdc.Cell(+x0 & -x1 & +y2 & -y3 & +z1 & -z2, lattice_1, translation=center)
300
301center += np.array([pitch * 17, 0.0, 0.0])
302assembly_2 = mcdc.Cell(+x1 & -x2 & +y2 & -y3 & +z1 & -z2, lattice_2, translation=center)
303
304center += np.array([-pitch * 17, -pitch * 17, 0.0])
305assembly_3 = mcdc.Cell(+x0 & -x1 & +y1 & -y2 & +z1 & -z2, lattice_3, translation=center)
306
307center += np.array([pitch * 17, 0.0, 0.0])
308assembly_4 = mcdc.Cell(+x1 & -x2 & +y1 & -y2 & +z1 & -z2, lattice_4, translation=center)
309
310# Bottom reflector cell
311reflector_bottom = mcdc.Cell(+x0 & -x3 & +y0 & -y3 & +z0 & -z1, mat_mod)
312
313# Side reflectors
314reflector_south = mcdc.Cell(+x0 & -x3 & +y0 & -y1 & +z1 & -z2, mat_mod)
315reflector_east = mcdc.Cell(+x2 & -x3 & +y1 & -y3 & +z1 & -z2, mat_mod)
316
317# Root universe
318mcdc.simulation.set_root_universe(
319    cells=[
320        assembly_1,
321        assembly_2,
322        assembly_3,
323        assembly_4,
324        reflector_bottom,
325        reflector_south,
326        reflector_east,
327    ],
328)
329
330# =============================================================================
331# Set source
332# =============================================================================
333# Throughout the active center pin of Assembly four, at highest energy,
334# for the first 15 seconds
335
336source = mcdc.Source(
337    x=np.array([pitch * 17 * 3 / 2] * 2) + np.array([-pitch / 2, +pitch / 2]),
338    y=np.array([-pitch * 17 * 3 / 2] * 2) + np.array([-pitch / 2, +pitch / 2]),
339    z=[-core_height / 2, core_height / 2],
340    isotropic=True,
341    energy_group=0,  # Highest energy
342    time=[0.0, 15.0],
343)
344
345# =============================================================================
346# Set tallies, settings, techniques and run MC/DC
347# =============================================================================
348
349# Tallies
350Nt = 100
351Nx = 17 * 2
352Ny = 17 * 2
353Nz = 17 * 6
354t = np.linspace(0.0, 20.0, Nt + 1)
355x = np.linspace(0.0, pitch * 17 * 2, Nx + 1)
356y = np.linspace(-pitch * 17 * 2, 0.0, Ny + 1)
357z = np.linspace(-core_height / 2, core_height / 2, Nz + 1)
358mesh = mcdc.MeshStructured(x=x, y=y, z=z)
359mcdc.Tally(mesh=mesh, scores=["fission"], time=t)
360
361# Settings
362mcdc.settings.N_particle = 10000
363mcdc.settings.N_batch = 2
364mcdc.settings.active_bank_buffer = 1000
365
366# Run
367mcdc.run()

How to Run

From the repository root run:

python examples/c5g7/transient/input.py

Expected Output

HDF5 tallies with time-resolved fission rates and PNG visualisations for fission and relative standard deviation per time step produced by the companion plotting scripts.