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

How to Run#

From inside examples/c5g7/transient run:

python 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.