MCNP criticality search Python script

One of the occasional tasks of a nuclear engineering student is to perform a criticality search using some parameter. That is, by varying (for example) the radius r of a uranium sphere, at some point the neutron production (from fission) and loss (from absorption and leakage) terms exactly balance, and the reactor is critical. The issue is complicated a little when using a Monte Carlo simulation, since the k-eigenvalue (unity when critical) is only reported with a certain probability or tolerance; running more particles will decrease the uncertainty.

Undergraduates, being unthinking animals, usually try this arduous task by hand, laboriously changing the parameter in an input file, running it, and repeating. As a graduate student, such tediousness is beneath me, as I would rather spend an equal amount of time writing a script to do the task for me. (Experience dictates that most of the time, some error or revision will require re-doing the task at least once, which means that a script will often save time by a factor of two or more. Besides, scripts are reusable and an educational experience.)

These Python scripts (gzipped archive) will do just that. It's a fairly simple binary search (no Newton–Krylov iterations) partly because of the complicating factor of the uncertainty in each answer. When criticality is reported to within the experimental standard deviation, the script increases the number of particles until such a parameter is determined that results in criticality within the user's specified tolerance. The script will run MCNP, process the output to determine k-effective, and iterate to your heart's content. If you have matplotlib installed, it will also produce a plot of how k varies with the parameter:

The MCNP search script is by no means perfect, and it assumes you have mcnp5 in your path. You have to know some Python to create a subclass of McnpSolver for your particular problem.

Included are scripts for finding the critical parameters of HEU cubes, cylinders, hemispheres, etc. There's also a script for generating a spherical shell volume tally to calculate the radial flux of the converged critical sphere.

Future updates of the MCNP criticality script will be posted on my projects page, although you should not expect there to be any. If you use this script, please cite its use if you report or present the results.