# A toy Monte Carlo criticality code

Here's a quick implementation of a Monte Carlo neutron transport code for calculating the criticality eigenvalue (i.e. \(k_{\mathrm{eff}}\)) for a fissile system (in slab geometry).

This complements other 1D toy codes for radiation transport I've written:
`snes`

, which calculates
\(k_{\mathrm{eff}}\) using the (deterministic) discrete-ordinates method,
and `fcimc`

, an
implicit Monte Carlo thermal
radiation code that models a Marshak wave in an optically thick medium.

## Transport equation

A highly simplified form of the transport equation in planar geometry is:

where the total, fission, and scattering macroscopic cross-sections, \(\Sigma_t\), \(\Sigma_f\) and \(\Sigma_s\) respectively, and the mean number of neutrons per fission, \(\nu\), are all constant. This is equivalent to assuming that the system is homegeneous, scattering is isotropic, and there is no energy dependence.

The system is also made time-independent by including a scaling factor in the fission source that enforces an exact balance between neutron loss due to leakage and collision (on the left-hand side) and production (on the right-hand side, if we think of scattering as both the loss and the production of a single neutron). This scaling factor is by definition equal to \(k_{\mathrm{eff}}\), which is what we want to determine.

As a result of all this, the angular flux, \(\psi\), is a function only of the single space variable \(x\) and the direction cosine \(\mu\).

## Monte Carlo algorithm

The outline of the Monte Carlo solution algorithm is given by the following Python code.

```
def run(num_particles):
"""
A single independent run with a fixed number of particles.
"""
# Code to get object `cfg` that contains the paramaters of the problem
# has been omitted
cfg = ...
# Get a uniformly distributed set of starting positions
start_positions = [
sample_position(cfg.slab_thickness_cm)
for _ in range(num_particles)
]
# Reset all the tallies to zero
tallies = initialise_tallies()
# Main loop over particles
for _ in range(num_particles):
# Particle history
tallies = simulate_single_history(cfg, tallies, start_positions)
# Can happen for small numbers of starting particles
if tallies["collision"] == 0:
sys.exit("Zero collisions")
# Estimate k_eff
keff = cfg.nu * tallies["fission"] / num_particles
c = tallies["secondary"] / tallies["collision"]
if verbose:
print(tallies)
print(
"Fission",
tallies["fission"] / tallies["collision"],
cfg.fission_prob
)
print(
"Scatter",
tallies["scatter"] / tallies["collision"],
cfg.scatter_prob,
)
print(
"Capture",
tallies["capture"] / tallies["collision"],
1 - cfg.fission_prob - cfg.scatter_prob,
)
print("c", c)
print(cfg.nu * cfg.fission_prob / (c - cfg.scatter_prob))
print(f"keff = {keff}")
# Sanity checks
assert tallies["history"] == num_particles
assert (
tallies["capture"] + tallies["leakage"] + tallies["fission"]
== num_particles
)
assert (
tallies["scatter"] + tallies["fission"] + tallies["capture"]
== tallies["collision"]
)
return keff
```

The starting positions of the neutrons are randomly distributed throughout the width of the slab. The various tallies we want are set up, such as the total number of each type of interaction, which are needed later for estimating \(k_{\mathrm{eff}}\). All the particle histories are looped over, and a collision esimator for \(k_{\mathrm{eff}}\) is computed. Various print and sanity-checks (that can be omitted from the final code) are shown.

If \(k_{\mathrm{eff}}\) is unfamiliar, it's directly analagous to the reproduction number or "R number" we heard about in the Covid pandemic. The reproduction number is the number of secondary infections per infected person; \(k_{\mathrm{eff}}\) is the number of secondary neutrons (created by fission) per neutron. We estimate it by finding the ratio of the total number of new neutrons created (\(\nu N_{\mathrm{fissions}}\)) to those removed (\(N_{\mathrm{capture}} + N_{\mathrm{leakage}} + N_{\mathrm{fission}}\), which is equivalent to the size of the initial population).

The outline of a single particle history is:

```
def simulate_single_history(
cfg,
tallies,
start_positions,
):
"""
Function to simulate a single neutron history in a 1D slab.
"""
tallies["history"] += 1
current_position = start_positions.pop()
while True:
# Free flight to next reaction/collision
direction_cosine = sample_direction_cosine()
scatter_distance = sample_scattering_distance(cfg.mean_free_path)
current_position = update_neutron_position(
current_position,
cfg.slab_thickness_cm,
direction_cosine,
cfg.left_boundary_condition,
scatter_distance,
)
# Leakage
if current_position < 0:
tallies["leakage"] += 1
return tallies
# Collisions
tallies["collision"] += 1
interaction_type = sample_interaction_type(
cfg.scatter_prob,
cfg.fission_prob
)
tallies[interaction_type] += 1
# Capture
if interaction_type == "capture":
return tallies
# Fission
if interaction_type == "fission":
num_secondaries = sample_neutrons_emitted(cfg.nu)
tallies["secondary"] += num_secondaries
return tallies
# Scattering
tallies["secondary"] += 1
```

First a starting location is popped off the pre-generated list. (At this stage,
there's no reason to pre-generate the list and pass it to the particle-history
function in this way; the `sample_initial_position()`

function could be called
directly here instead. But doing it this way will come in handy later.)

The main loop then teleports the neutron from interaction site to interaction
site, until it either leaks out of the boundaries or gets swallowed up in a
fission or capture event. For each leg, the direction cosine is randomly chosen
from a uniform distribution between -1 and 1, and the position is updated by
sampling, within the `update_neutron_position`

function, the distance to the
next collision from an exponential distribution, given a mean-free-path of
\(1/\Sigma_t\). This routine also handles the boundary conditions such that the
returned position is set to a negative marker value to indicate the neutron
has left the system.

For neutrons that remain, the interaction type is selected in
`sample_interaction_type`

using roulette sampling, and the corresponding
tally is incremented. These tallies are the main output, for the estimator of
\(k_{\mathrm{eff}}\), as seen above.

## Basic results

We'll see soon that this algorithm isn't sufficient, but let's look at how it performs.

This shows the mean and standard deviation of ten runs at each particle number, for two different estimators of \(k_{\mathrm{eff}}\), and although at very high particle numbers the statistics are good, the code is converging to the wrong answer.

The solution to this is to run multiple generations inside the code. The first generation of neutrons refers to the initial batch that are created directly. But each fission caused by a neutron in this first generation gives rise to further neutrons, which collectively make up the second generation, and so on.

If we make a record of the location of and the number of neutrons emitted from each fission in the first generation, then we have a second generation which we can also simulate, and so on.

## Simulating multiple generations

Here's how the code for this looks. Overall:

```
def simulate_single_history(
cfg,
tallies,
start_positions,
):
"""
Function to simulate a single neutron history in a 1D slab.
"""
tallies["history"] += 1
current_position = start_positions.pop()
new_positions = []
while True:
# Free flight to next reaction/collision
direction_cosine = sample_direction_cosine()
scatter_distance = sample_scattering_distance(cfg.mean_free_path)
current_position = update_neutron_position(
current_position,
cfg.slab_thickness_cm,
direction_cosine,
cfg.left_boundary_condition,
scatter_distance,
)
# Leakage
if current_position < 0:
tallies["leakage"] += 1
return tallies, new_positions
# Collisions
tallies["collision"] += 1
interaction_type = sample_interaction_type(
cfg.scatter_prob,
cfg.fission_prob
)
tallies[interaction_type] += 1
# Capture
if interaction_type == "capture":
return tallies, new_positions
# Fission
if interaction_type == "fission":
num_secondaries = sample_neutrons_emitted(cfg.nu)
tallies["secondary"] += num_secondaries
for fission_neutron in range(num_secondaries):
new_positions.append(current_position)
return tallies, new_positions
# Scattering
tallies["secondary"] += 1
def run(num_generations, num_particles):
"""
A single independent run with a fixed number of generations and particles.
"""
# Code to get object `cfg` that contains the paramaters of the problem
# has been omitted
cfg = ...
# Initialise the per-generation number of particles
num_particles_in_generation = cfg.num_particles
# Get a uniformly distributed set of starting positions for the initial
# generation of particles
start_positions = [
sample_position(cfg.slab_thickness_cm)
for _ in range(num_particles_in_generation)
]
# Lists for storing the estimates of k_effective across generations
# k1 and k2 are two different estimators for k_effective
k1 = []
k2 = []
for gen in range(cfg.num_generations):
# Reset all the tallies to zero for this generation
tallies = initialise_tallies()
# Initialise a list for the start positions of the next generation,
# arising from fission in this generation
next_start_positions = []
# Main loop over particles in this generation
for _ in range(num_particles_in_generation):
# Particle history
tallies, new_start_positions = simulate_single_history(
cfg,
tallies,
start_positions,
)
# Add the new start positions from this history to the global list
# of start positions for the next generation
next_start_positions += new_start_positions
if tallies["collision"] == 0:
sys.exit("Zero collisions")
# Estimate k_eff
k1.append(
cfg.nu
* tallies["fission"]
/ (tallies["capture"] + tallies["leakage"] + tallies["fission"])
)
k2.append(len(next_start_positions) / num_particles_in_generation)
c = tallies["secondary"] / tallies["collision"]
verbose = False
if verbose:
print(tallies)
print(
"Fission",
tallies["fission"] / tallies["collision"],
cfg.fission_prob
)
print(
"Scatter",
tallies["scatter"] / tallies["collision"],
cfg.scatter_prob,
)
print(
"Absorbs",
tallies["capture"] / tallies["collision"],
1 - cfg.fission_prob - cfg.scatter_prob,
)
print("c", c)
print(cfg.nu * cfg.fission_prob / (c - cfg.scatter_prob))
print(f"k1 = {k1}")
print(f"k2 = {k2}")
# Sanity checks
assert tallies["history"] == num_particles_in_generation
assert (
tallies["capture"] + tallies["leakage"] + tallies["fission"]
== num_particles_in_generation
)
assert (
tallies["scatter"] + tallies["fission"] + tallies["capture"]
== tallies["collision"]
)
num_particles_in_generation = len(next_start_positions)
if num_particles_in_generation == 0:
sys.exit("Zero particles")
start_positions = next_start_positions
return k1, k2
```

A loop over generations has been added, and in each generation a list of
starting positions for the next generation is constructed, by having
`simulate_single_history`

return the starting positions of the fission neutrons
produced in each history. (Explicitly recording every fission neutron location
is a very simplistic way of seeding the next generation, which is only
practical because of the overall small number of neutrons here. This is only
for demonstration.)

Note also that the second estimator for \(k_{\mathrm{eff}}\) is now shown; \(k_{\mathrm{eff}}\) is also, by definition, the ratio of the number of neutrons between subsequent generations, which we can compute directly.

## Multi-generation results

Let's look at how our estimate of \(k_{\mathrm{eff}}\) behaves as we simulate more generations. The next figures show the results of the same particle convergence study but where for each run we simulate two consecutive generations, taking the estimate of \(k_{\mathrm{eff}}\) from the final one only:

Then three generations:

which looks pretty good. If we look a the variation with generation for a fixed number of particles (128k), we see again that by the third generation things look good:

## Fission rate

Plotting a histogram of the start position of all the neutrons in each generations gives a representation of the fission rate, which is itself representative of the scalar flux profile, since the cross-sections are constant. Here we can see the initial flat distribution in the zeroth generation, followed by a rapid convergence to the expected flux profile by around the third generation:

Getting this distribution right is what allows the correct estimation of \(k_{\mathrm{eff}}\).

## Conclusion

The full code is available at https://github.com/msleigh/mccc.