Skip to content

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:

\[ \begin{align} \mu \frac{\partial\psi(x, \mu)}{\partial x} + \Sigma_t \psi(x, \mu) = \\ \frac{\Sigma_s}{2} \int_{-1}^{1} \psi(x, \mu') d\mu' + \frac{\nu \Sigma_f}{2k_{\mathrm{eff}}} \int_{-1}^{1} \psi(x, \mu') d\mu' \end{align} \]

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.