First MC/DC Simulation#

This tutorial constructs and runs a one-group shielding calculation. It assumes familiarity with the basic concepts of Monte Carlo radiation transport.

Monte Carlo transport inputs are generally assembled from the same core components: materials, geometry, particle sources, tallies, and simulation settings. In MC/DC, these components are represented by Python objects and assembled into a mcdc.Simulation.

This tutorial applies the general Simulation Lifecycle to one complete problem.

The example uses multigroup data defined directly in the input, so it does not require an external nuclear-data library. Its complete, executable source is available in examples/slab_shielding.

MC/DC Workflow#

An MC/DC calculation follows five main steps:

  1. Construct the materials, geometry, sources, and tallies.

  2. Create a simulation and attach the model objects to it.

  3. Configure settings and compile the connected object graph.

  4. Visualize or run the prepared model.

  5. Generate and post-process the output.

Important

Validate the geometry, source distribution, tally definitions, and settings before interpreting simulation results. A successfully executed calculation is not necessarily a correctly specified physical model.

Problem Description#

The model contains two adjacent slab regions:

  • A mostly scattering source region over \(0 < z < 2\) cm.

  • A more strongly absorbing shield over \(2 < z < 6\) cm.

Both outer boundaries are vacuum. Particles are emitted isotropically throughout the source region, and a mesh tally records the flux across the entire domain.

One-group material data#

Region

Range (cm)

\(\Sigma_c\) (cm-1)

\(\Sigma_s\) (cm-1)

Source region

\(0 < z < 2\)

0.1

0.9

Shield

\(2 < z < 6\)

0.7

0.3

The total cross section is \(1.0\ \text{cm}^{-1}\) in both regions. Changing the capture-to-scatter ratio isolates the effect of the shield on the flux distribution.

Building the Input#

Imports and Simulation#

NumPy provides the numerical arrays used for cross sections and tally grids. The Simulation instance collects the model and controls its execution:

import numpy as np

import mcdc


simulation = mcdc.Simulation("One-group slab shielding")

Materials#

MaterialMG represents multigroup interaction data. A one-element capture array and a \(1 \times 1\) scattering matrix define a one-group material:

source_region_material = mcdc.MaterialMG(
    capture=np.array([0.1]),
    scatter=np.array([[0.9]]),
)
shield_material = mcdc.MaterialMG(
    capture=np.array([0.7]),
    scatter=np.array([[0.3]]),
)

Continuous-energy calculations instead use mcdc.Material and require an MC/DC nuclear-data library. See Generating a Data Library from ACE Files for configuration instructions.

Geometry#

Three z-planes define the two slab regions. The outer planes use vacuum boundary conditions; the plane at \(z=2\) cm is an internal interface:

left = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum")
interface = mcdc.Surface.PlaneZ(z=2.0)
right = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum")

A cell combines a region with the material that fills it. A positive half-space selects points above a PlaneZ, while a negative half-space selects points below it:

source_cell = mcdc.Cell(
    region=+left & -interface,
    fill=source_region_material,
)
shield_cell = mcdc.Cell(
    region=+interface & -right,
    fill=shield_material,
)
simulation.set_model([source_cell, shield_cell])

The cells are the roots of the geometry. MC/DC reaches their materials and surfaces when it compiles the simulation, so those objects do not require separate setter calls.

Source#

The source emits group-0 particles isotropically and uniformly throughout the first cell:

source = mcdc.Source(
    z=[0.0, 2.0],
    isotropic=True,
    energy_group=0,
)
simulation.set_sources([source])

Tallies#

The structured mesh divides the domain into 60 equal spatial bins. The tally scores the track-length estimate of flux in each bin:

mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61))
flux_tally = mcdc.Tally(
    name="slab_flux",
    mesh=mesh,
    scores=["flux"],
)
simulation.set_tallies([flux_tally])

Naming the tally makes its location in the output file predictable: tallies/slab_flux.

Settings and Execution#

This example runs 1,000 particle histories in each of 10 statistically independent batches. Multiple batches allow MC/DC to estimate the standard deviation of each tally bin:

simulation.settings.N_particle = 1_000
simulation.settings.N_batch = 10
simulation.settings.output_name = "slab_shielding"

simulation.run()

The calculation writes its results to slab_shielding.h5.

Complete Input#

The complete runnable input is embedded directly from examples/slab_shielding/input.py:

 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()

Visualizing the Model#

Before running transport, insert the following call immediately before simulation.run() to render an x-z slice of the material geometry:

simulation.visualize_model(
    vis_plane="xz",
    x=[-1.0, 1.0],
    y=0.0,
    z=[0.0, 6.0],
    pixels=(100, 300),
    colors=None,
    time=[0.0],
    save_as="slab_shielding_geometry",
)

The image is saved as slab_shielding_geometry.png. Visualization compiles the current model when necessary.

Running the Example#

Enter the problem directory, then run the input in pure Python mode:

cd examples/slab_shielding
python input.py

Pure Python mode avoids compilation overhead and is suitable for checking a small model. For accelerated or parallel calculations, see Running MC/DC on CPUs, Running MC/DC on GPUs, and Batch Job Scripts.

Post-processing#

MC/DC writes tally results and runtime information to HDF5. The companion script reads the spatial grid, normalizes the flux and standard deviation by the mesh-bin widths, and plots the result:

 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)

After the transport calculation finishes, run the companion script from the same problem directory:

python process-output.py

The script writes slab_shielding_flux.png. The dashed line marks the material interface at \(z=2\) cm. The flux is expected to decrease more rapidly in the shield because capture accounts for a larger fraction of its total cross section.

Next Steps#

After running the original problem, useful variations include:

  • Increase N_particle and compare the reported standard deviation.

  • Change the shield capture and scattering cross sections.

  • Move the material interface and observe the change in attenuation.

  • Add an energy group or another spatial region.

  • Add a surface-crossing tally at the material interface.

See Example Problems for examples involving lattices, moving geometry, time-dependent transport, and reactor benchmarks.