Using a Fixed Source
In this example, we’ll modify the previous pin-cell model to use a fixed source instead of a criticality source.
Generating the Model
The pin-cell can be generated from the MCNPy example classes.
import mcnpy as mp
# Pin-cell example class
pincell = mp.example.Pincell()
# Deck object for the pin-cell
deck = pincell.deck
Removing Criticality Source Cards
First, we’ll remove the cards associated with the criticality source.
cards_to_remove = []
for card in deck.src_settings:
if isinstance(card, mp.CriticalitySourcePoints):
cards_to_remove.append(card)
elif isinstance(card, mp.CriticalitySource):
cards_to_remove.append(card)
deck -= cards_to_remove
Note
deck.src_settings returns a list and not a dict since many source cards do not have IDs to use as keys.
Adding a Fixed Source
With the old cards removed, we can create a new source definition.
sdef = mp.Source(particle=mp.Particle.NEUTRON, position=[0,0,0])
deck += sdef
Important
Due to the lengthy list of options for mcnpy.Source, the class only accepts keyword arguments.
This will produce an isotropic neutron source at (0, 0, 0) using MCNP’s default source particle energy. To make out source more interesting, let’s add an energy distribution. This begins by creating Source information and probability cards.
# Source information - the energy values
energies = [0.0253e-6, 1e-6, 1, 10] # arbitrary distribution
si = mp.SourceInfo(name=10, option=mp.SourceInfoOption.DISCRETE, values=energies)
# Source probability - the probability of each value
probability = [1, 1, 1, 1]
sp = mp.SourceProbability(name=10, option=mp.SourceProbabilityOption.PROBABILITIES,
values=probability)
deck += [si, sp]
The distribution can now be added to our source.
sdef.energies = [si]
Important
If an mcnpy.Distributions object is not explicitly created, a list of objects of the type mcnpy.Distribution must be used.
The deck now has a source definition with an energy distribution. To use MCNP’s built-in distributions, one could omit si and substitute sp with something like
# Watt fission spectrum for thermal neutron-induced U235 fission.
sp = mp.SourceProbability.Function(name=10).watt(a=0.988 b=2.249)
Setting the Simulation Cutoff
To control when the MCNP simulation terminates, we need to add a cutoff card like NPS or CTME. Here we’ll use an mcnpy.Cutoff.History object to end the simulation after a set number of particle histories.
# Terminate after 1e5 histories
deck += mp.Cutoff.History(histories=1e5)
Note
Choose a number of particle histories that is appropriate for your system and target uncertainty. More histories result in lower uncertainties, but proportionally longer simulations. E.g. 1E+06 histories will take approximately 10 times as long as 1E+05 histories.
Adding Tallies
Tallies tell MCNP which quantities to measure or “tally up” and report in the simulation outputs. One of the most common quantities to measure is flux. Here, we’ll add a cell flux tally for neutrons to our model. Defining a tally usually begins with defining its cell or surface bins.
In their simplest form, a bin will be the cell or surface a tally measures in/on. Additional bins (energy, time, etc.) are defined using their own cards. Much like in native MCNP syntax, bins can be defined with operators.
Operation |
Symbol (Python) |
Symbol (MCNP) |
Semantics |
|---|---|---|---|
Union |
|
|
OR |
Level |
|
space |
Grouping of cells |
Levels |
|
|
Within repeated structure levels |
Index |
|
|
Elements of lattice cell |
For our fuel pin, let’s tally the flux inside the fuel region. At a minimum, we must provide the tally with a particles list and bins. MCNPy will automatically number any unnamed tallies using MCNP’s tally ID conventions.
Important
Manually defined tally names must follow MCNP conventions. E.g. 4, 14, 34, etc. are all valid names for an “F4” cell flux tally.
fuel_flux = mp.Tally.CellFlux(particles=[mp.Particle.NEUTRON], bins=[deck.cells[1]])
deck += fuel_flux
Important
Even though particles and bins use only 1 item, they must be provided as lists.
It is also straightforward to add some energy bins to our tally. This will produce a tally with thermal and fast energy bins. Of course, Python makes generating or reading in lengthy bins structures simple as well.
ergs = [1e-6, 15]
deck += mp.Tally.Bins.Energies(tally=fuel_flux, max_energies=ergs,
particles=[mp.Particle.NEUTRON])
Calling print(deck.serialize)) will then produce the following deck.
$ This file was written with mcnpy
1 1 -10.0 -1 -4 IMP:N=1.0
2 3 -6.6 2 -3 -4 IMP:N=1.0
3 2 -1.0 3 -4 IMP:N=1.0
4 4 -1.78E-4 1 -2 -4 IMP:N=1.0
5 0 4 IMP:N=0.0
1 CZ 0.39
2 CZ 0.4
3 CZ 0.46
*4 RPP -0.63 0.63 -0.63 0.63 -0.63 0.63
M1
92238 0.9696278737791809
92235 0.030372126220819046
8016 2.0
M2
1001 2.0
8016 1.0
M3
40000 -0.982
50000 -0.015
26000 -0.002
24000 -0.001
M4
1001 0.5
2004 0.5
MT1 O2-U U-O2
MT2 LWTR
PRDMP 0.0 -60.0 1.0
SDEF ERG (D10) POS 0.0 0.0 0.0 PAR N
SI10 L 2.53E-8 1.0E-6 1.0 10.0
SP10 1.0 1.0 1.0 1.0
NPS 100000.0
F4:N 1
E4:N 1.0E-6 15.0