Slab Shielding#

This one-group fixed-source problem introduces the complete MC/DC workflow with two materials, two slab cells, an isotropic source, and a mesh flux tally. The First MC/DC Simulation guide explains the input step by step.

Full Input#

 1import numpy as np
 2
 3import mcdc
 4
 5simulation = mcdc.Simulation("One-group slab shielding")
 6
 7# Materials
 8source_region_material = mcdc.MaterialMG(
 9    capture=np.array([0.1]),
10    scatter=np.array([[0.9]]),
11)
12shield_material = mcdc.MaterialMG(
13    capture=np.array([0.7]),
14    scatter=np.array([[0.3]]),
15)
16
17# Geometry
18left = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum")
19interface = mcdc.Surface.PlaneZ(z=2.0)
20right = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum")
21
22source_cell = mcdc.Cell(
23    region=+left & -interface,
24    fill=source_region_material,
25)
26shield_cell = mcdc.Cell(
27    region=+interface & -right,
28    fill=shield_material,
29)
30simulation.set_model([source_cell, shield_cell])
31
32# Source
33source = mcdc.Source(
34    z=[0.0, 2.0],
35    isotropic=True,
36    energy_group=0,
37)
38simulation.set_sources([source])
39
40# Tally
41mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61))
42flux_tally = mcdc.Tally(
43    name="slab_flux",
44    mesh=mesh,
45    scores=["flux"],
46)
47simulation.set_tallies([flux_tally])
48
49# Settings
50simulation.settings.N_particle = 1_000
51simulation.settings.N_batch = 10
52simulation.settings.output_name = "slab_shielding"
53
54# Visualize model
55"""
56simulation.visualize_model(
57    vis_plane="xz",
58    x=[-1.0, 1.0],
59    y=0.0,
60    z=[0.0, 6.0],
61    pixels=(100, 300),
62    colors=None,
63    time=[0.0],
64    save_as="slab_shielding_geometry",
65)
66"""
67
68# Run simulation
69simulation.run()

Post-processing#

 1import h5py
 2import matplotlib.pyplot as plt
 3
 4with h5py.File("slab_shielding.h5", "r") as output:
 5    tally = output["tallies/slab_flux"]
 6    z = tally["grid/z"][:]
 7    flux = tally["flux/mean"][:]
 8    flux_sdev = tally["flux/sdev"][:]
 9
10dz = z[1:] - z[:-1]
11z_mid = 0.5 * (z[:-1] + z[1:])
12flux /= dz
13flux_sdev /= dz
14
15figure, axis = plt.subplots()
16axis.plot(z_mid, flux, label="Flux")
17axis.fill_between(
18    z_mid,
19    flux - flux_sdev,
20    flux + flux_sdev,
21    alpha=0.25,
22    label="Standard deviation",
23)
24axis.axvline(2.0, color="black", linestyle="--", label="Material interface")
25axis.set_xlabel("z [cm]")
26axis.set_ylabel("Flux")
27axis.grid()
28axis.legend()
29figure.tight_layout()
30figure.savefig("slab_shielding_flux.png", dpi=150)

How to Run#

From inside examples/slab_shielding:

python input.py
python process-output.py

The transport calculation writes slab_shielding.h5. The post-processing script reads its flux tally and writes slab_shielding_flux.png.