aparler/openmc-designs
1
0
Fork 0
forked from crp/openmc-designs
Code Pull requests Activity

Added new stress test, 'validation'

Browse source
This commit is contained in:
Adam Parler 2023-12-16 19:41:24 -08:00
parent fe2b97b7dc
commit 05285cd188
5 changed files with 2722 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

509
TRIGA/model_creation.py Normal file
View file

@ -0,0 +1,509 @@
#!/usr/bin/env python3
# # A TRIGA geometry
# This notebook can be used as a template for modeling TRIGA reactors.
import openmc
import numpy as np
# Materials definitions
# Borated water
water = openmc.Material(name='Borated Water')
water.set_density('g/cm3', 0.740582)
water.add_nuclide('H1', 4.9457e-2)
water.add_nuclide('O16', 2.4732e-2)
water.add_nuclide('B10', 8.0042e-6)
# 20% enriched uranium zirconium hydride fuel
uzrh = openmc.Material(name='UZrH')
uzrh.set_density('g/cm3', 6.128)
uzrh.add_nuclide('U235', .02376, 'wo')
uzrh.add_nuclide('U238', .09619, 'wo')
uzrh.add_element('H', .03, 'wo')
uzrh.add_element('Zr', .85, 'wo')
molybdenum = openmc.Material(name='Molybdenum')
molybdenum.add_element('Mo', 1.0)
molybdenum.set_density('g/cm3', 10.22)
graphite = openmc.Material(name='Graphite')
graphite.set_density('g/cm3', 1.70)
graphite.add_element('C', 1.0)
graphite.add_s_alpha_beta('c_Graphite')
# Stainless steel
ss304 = openmc.Material(name='Stainless Steel 304')
ss304.set_density('g/cm3', 8.0)
ss304.add_element('C',.002,'wo')
ss304.add_element('Si',.004,'wo')
ss304.add_element('P',.0003,'wo')
ss304.add_element('S',.0002,'wo')
ss304.add_element('V',.003,'wo')
ss304.add_element('Cr',.115,'wo')
ss304.add_element('Mn',.006,'wo')
ss304.add_element('Fe',.8495,'wo')
ss304.add_element('Ni',.005,'wo')
ss304.add_element('Mo',.01,'wo')
ss304.add_element('W',.005,'wo')
# Boron carbide
b4c = openmc.Material(name='Boron Carbide')
b4c.set_density('g/cm3', 2.52)
b4c.add_element('B', 4)
b4c.add_element('C', 1)
zirconium = openmc.Material(name='Zirconium')
zirconium.add_element('Zr', 1.0)
zirconium.set_density('g/cm3', 6.506)
void = openmc.Material(name='Void')
void.set_density('g/cm3', 0.001205)
void.add_element('Ni', 0.755268, 'wo')
void.add_element('C', 0.000124, 'wo')
void.add_element('O', 0.231781, 'wo')
void.add_element('Ar', 0.012827, 'wo')
aluminum = openmc.Material(name='Aluminum')
aluminum.add_element('Al', 1.0)
aluminum.set_density('g/cm3', 2.6)
# Instantiate a Materials collection and export to xml
materials_file = openmc.Materials([aluminum, zirconium, b4c, ss304, graphite, molybdenum, uzrh, water, void])
materials_file.export_to_xml()
# Geometry definitions for the fuel rod
rod_outer_radius = openmc.ZCylinder(r=0.635)
uzrh_outer_radius = openmc.ZCylinder(r=3.6449)
molybdenum_outer_radius = openmc.ZCylinder(r=3.6449)
graphite_outer_radius = openmc.ZCylinder(r=3.6449)
ss304_outer_radius = openmc.ZCylinder(r=3.6449)
clad_outer_radius = openmc.ZCylinder(r=3.75412)
empty_space_min = openmc.ZPlane(z0=-114.3)
empty_space_max = openmc.ZPlane(z0=-55.079265)
rod_min = openmc.ZPlane(z0=-55.079265)
rod_max = openmc.ZPlane(z0=-16.979265)
uzrh_min = openmc.ZPlane(z0=-55.079265)
uzrh_max = openmc.ZPlane(z0=-16.979265)
molybdenum_min = openmc.ZPlane(z0=-55.15864)
molybdenum_max = openmc.ZPlane(z0=-55.079265)
graphite_upper_min = openmc.ZPlane(z0=-16.979265)
graphite_upper_max = openmc.ZPlane(z0=-10.375265)
graphite_lower_min = openmc.ZPlane(z0=-64.621664)
graphite_lower_max = openmc.ZPlane(z0=-55.15864)
ss304_upper_min = openmc.ZPlane(z0=-10.375265)
ss304_upper_max = openmc.ZPlane(z0=+0)
ss304_lower_min = openmc.ZPlane(z0=-72.0598)
ss304_lower_max = openmc.ZPlane(z0=-64.621664)
clad_min = openmc.ZPlane(z0=-72.0598)
clad_max = openmc.ZPlane(z0=0)
# Create a Universe to encapsulate the fuel rod
fuel_universe = openmc.Universe(name='UZrH Fuel Universe')
# Create rod cell
rod_cell = openmc.Cell(name='Zr Rod')
rod_cell.fill = zirconium
rod_cell.region = -rod_outer_radius & +rod_min & -rod_max
fuel_universe.add_cell(rod_cell)
# Create uzrh cell
uzrh_cell = openmc.Cell(name='UZrH')
uzrh_cell.fill = uzrh
uzrh_cell.region = +rod_outer_radius & -uzrh_outer_radius & +uzrh_min & -uzrh_max
fuel_universe.add_cell(uzrh_cell)
# Create molybdenum disk
molybdenum_cell = openmc.Cell(name='Molybdenum')
molybdenum_cell.fill = molybdenum
molybdenum_cell.region = -molybdenum_outer_radius & +molybdenum_min & -molybdenum_max
fuel_universe.add_cell(molybdenum_cell)
# Create upper graphite cell
graphite_upper_cell = openmc.Cell(name='Upper Graphite')
graphite_upper_cell.fill = graphite
graphite_upper_cell.region = -graphite_outer_radius & +graphite_upper_min & -graphite_upper_max
fuel_universe.add_cell(graphite_upper_cell)
# Create lower graphite cell
graphite_lower_cell = openmc.Cell(name='Lower Graphite')
graphite_lower_cell.fill = graphite
graphite_lower_cell.region = -graphite_outer_radius & +graphite_lower_min & -graphite_lower_max
fuel_universe.add_cell(graphite_lower_cell)
# Create upper ss304 cell
ss304_upper_cell = openmc.Cell(name='Upper Stainless Steel 304')
ss304_upper_cell.fill = ss304
ss304_upper_cell.region = -ss304_outer_radius & +ss304_upper_min & -ss304_upper_max
fuel_universe.add_cell(ss304_upper_cell)
# Create lower ss304 cell
ss304_lower_cell = openmc.Cell(name='Lower Stainless Steel 304')
ss304_lower_cell.fill = ss304
ss304_lower_cell.region = -ss304_outer_radius & +ss304_lower_min & -ss304_lower_max
fuel_universe.add_cell(ss304_lower_cell)
# Create clad cell
clad_cell = openmc.Cell(name='Stainless Steel 304 Cladding')
clad_cell.fill = ss304
clad_cell.region = -clad_outer_radius & +ss304_outer_radius & +clad_min & -clad_max
fuel_universe.add_cell(clad_cell)
# Create empty space cell
empty_space_cell = openmc.Cell(name='Empty space before fuel rod')
empty_space_cell.fill = water
empty_space_cell.region = -clad_outer_radius & +empty_space_min & -empty_space_max
fuel_universe.add_cell(empty_space_cell)
# Geometry definitions for the transient rod
void_outer_radius = openmc.ZCylinder(r=3.03276)
b4c_outer_radius = openmc.ZCylinder(r=3.03276)
clad_outer_radius = openmc.ZCylinder(r=3.38455)
aluminum_outer_radius = openmc.ZCylinder(r=3.03276)
aluminum_1_min = openmc.ZPlane(z0=-114.3)
aluminum_1_max = openmc.ZPlane(z0=-113.03)
void_1_min = openmc.ZPlane(z0=-113.03)
void_1_max = openmc.ZPlane(z0=-57.785)
aluminum_2_min = openmc.ZPlane(z0=-57.785)
aluminum_2_max = openmc.ZPlane(z0=-56.515)
b4c_min = openmc.ZPlane(z0=-56.515)
b4c_max = openmc.ZPlane(z0=-18.415)
void_2_min = openmc.ZPlane(z0=-18.415)
void_2_max = openmc.ZPlane(z0=-18.0975)
aluminum_3_min = openmc.ZPlane(z0=-18.0975)
aluminum_3_max = openmc.ZPlane(z0=-16.8275)
void_3_min = openmc.ZPlane(z0=-16.8275)
void_3_max = openmc.ZPlane(z0=-7.3025)
aluminum_4_min = openmc.ZPlane(z0=-7.3025)
aluminum_4_max = openmc.ZPlane(z0=0)
clad_min = openmc.ZPlane(z0=-114.3)
clad_max = openmc.ZPlane(z0=0)
# Create a Universe to encapsulate the transient rod
transient_universe = openmc.Universe(name='Transient Universe')
# Create void 1 cell
void_1_cell = openmc.Cell(name= 'Void 1')
void_1_cell.fill = void
void_1_cell.region = -void_outer_radius & +void_1_min & -void_1_max
transient_universe.add_cell(void_1_cell)
# Create void 2 cell
void_2_cell = openmc.Cell(name= 'Void 2')
void_2_cell.fill = void
void_2_cell.region = -void_outer_radius & +void_2_min & -void_2_max
transient_universe.add_cell(void_2_cell)
# Create void 3 cell
void_3_cell = openmc.Cell(name= 'Void 3')
void_3_cell.fill = void
void_3_cell.region = -void_outer_radius & +void_3_min & -void_3_max
transient_universe.add_cell(void_3_cell)
# Create b4c cell
b4c_cell = openmc.Cell(name='Boron Carbide')
b4c_cell.fill = b4c
b4c_cell.region = -b4c_outer_radius & -b4c_max & +b4c_min
transient_universe.add_cell(b4c_cell)
# Create aluminum 1 cell
aluminum_1_cell = openmc.Cell(name='Aluminum')
aluminum_1_cell.fill = aluminum
aluminum_1_cell.region = -aluminum_outer_radius & +aluminum_1_min & -aluminum_1_max
transient_universe.add_cell(aluminum_1_cell)
# Create aluminum 2 cell
aluminum_2_cell = openmc.Cell(name='Aluminum')
aluminum_2_cell.fill = aluminum
aluminum_2_cell.region = -aluminum_outer_radius & +aluminum_2_min & -aluminum_2_max
transient_universe.add_cell(aluminum_2_cell)
# Create aluminum 3 cell
aluminum_3_cell = openmc.Cell(name='Aluminum')
aluminum_3_cell.fill = aluminum
aluminum_3_cell.region = -aluminum_outer_radius & +aluminum_3_min & -aluminum_3_max
transient_universe.add_cell(aluminum_3_cell)
# Create aluminum 4 cell
aluminum_4_cell = openmc.Cell(name='Aluminum')
aluminum_4_cell.fill = aluminum
aluminum_4_cell.region = -aluminum_outer_radius & +aluminum_4_min & -aluminum_4_max
transient_universe.add_cell(aluminum_4_cell)
# Create a clad cell
clad_cell = openmc.Cell(name='Aluminum Cladding')
clad_cell.fill = aluminum
clad_cell.region = -clad_outer_radius & +b4c_outer_radius & +clad_min & -clad_max
transient_universe.add_cell(clad_cell)
# Geometry definitions for the control rod
rod_outer_radius = openmc.ZCylinder(r=0.635)
uzrh_outer_radius = openmc.ZCylinder(r=3.33375)
void_outer_radius = openmc.ZCylinder(r=3.33375)
b4c_outer_radius = openmc.ZCylinder(r=3.33375)
ss304_outer_radius = openmc.ZCylinder(r=3.33375)
clad_outer_radius = openmc.ZCylinder(r=3.38455)
ss304_5_min = openmc.ZPlane(z0=-114.3)
ss304_5_max = openmc.ZPlane(z0=-113.03)
void_4_min = openmc.ZPlane(z0=-113.03)
void_4_max = openmc.ZPlane(z0=-99.06)
ss304_4_min = openmc.ZPlane(z0=-99.06)
ss304_4_max = openmc.ZPlane(z0=-96.52)
rod_min = openmc.ZPlane(z0=-96.52)
rod_max = openmc.ZPlane(z0=-58.42)
uzrh_min = openmc.ZPlane(z0=-96.52)
uzrh_max = openmc.ZPlane(z0=-58.42)
void_1_min = openmc.ZPlane(z0=-58.42)
void_1_max = openmc.ZPlane(z0=-57.785)
ss304_1_min = openmc.ZPlane(z0=-57.785)
ss304_1_max = openmc.ZPlane(z0=-56.515)
b4c_1_min = openmc.ZPlane(z0=-56.515)
b4c_1_max = openmc.ZPlane(z0=-18.415)
void_2_min = openmc.ZPlane(z0=-18.415)
void_2_max = openmc.ZPlane(z0=-18.0975)
ss304_2_min = openmc.ZPlane(z0=-18.0975)
ss304_2_max = openmc.ZPlane(z0=-16.8275)
void_3_min = openmc.ZPlane(z0=-16.8275)
void_3_max = openmc.ZPlane(z0=-7.3025)
ss304_3_min = openmc.ZPlane(z0=-7.3025)
ss304_3_max = openmc.ZPlane(z0=0.0)
clad_min = openmc.ZPlane(z0=-114.3)
clad_max = openmc.ZPlane(z0=0.0)
# Create a Universe to encapsulate the control rod
control_universe = openmc.Universe(name='Control Universe')
# Create rod cell
rod_cell = openmc.Cell(name='Zr Rod')
rod_cell.fill = zirconium
rod_cell.region = -rod_outer_radius & +rod_min & -rod_max
control_universe.add_cell(rod_cell)
# Create uzrh cell
uzrh_cell = openmc.Cell(name='UZrH')
uzrh_cell.fill = uzrh
uzrh_cell.region = +rod_outer_radius & -uzrh_outer_radius & +uzrh_min & -uzrh_max
control_universe.add_cell(uzrh_cell)
# Create void 1 cell
void_1_cell = openmc.Cell(name= 'Void 1')
void_1_cell.fill = void
void_1_cell.region = -void_outer_radius & +void_1_min & -void_1_max
control_universe.add_cell(void_1_cell)
# Create void 2 cell
void_2_cell = openmc.Cell(name= 'Void 2')
void_2_cell.fill = void
void_2_cell.region = -void_outer_radius & +void_2_min & -void_2_max
control_universe.add_cell(void_2_cell)
# Create void 3 cell
void_3_cell = openmc.Cell(name= 'Void 3')
void_3_cell.fill = void
void_3_cell.region = -void_outer_radius & +void_3_min & -void_3_max
control_universe.add_cell(void_3_cell)
# Create void 4 cell
void_4_cell = openmc.Cell(name= 'Void 3')
void_4_cell.fill = void
void_4_cell.region = -void_outer_radius & +void_4_min & -void_4_max
control_universe.add_cell(void_4_cell)
# Create ss304 1 cell
ss304_1_cell = openmc.Cell(name='Stainless Steel 304 Cell 1')
ss304_1_cell.fill = ss304
ss304_1_cell.region = -ss304_outer_radius & +ss304_1_min & -ss304_1_max
control_universe.add_cell(ss304_1_cell)
# Create ss304 2 cell
ss304_2_cell = openmc.Cell(name='Stainless Steel 304 Cell 2')
ss304_2_cell.fill = ss304
ss304_2_cell.region = -ss304_outer_radius & +ss304_2_min & -ss304_2_max
control_universe.add_cell(ss304_2_cell)
# Create ss304 3 cell
ss304_3_cell = openmc.Cell(name='Stainless Steel 304 Cell 3')
ss304_3_cell.fill = ss304
ss304_3_cell.region = -ss304_outer_radius & +ss304_3_min & -ss304_3_max
control_universe.add_cell(ss304_3_cell)
# Create ss304 4 cell
ss304_4_cell = openmc.Cell(name='Stainless Steel 304 Cell 4')
ss304_4_cell.fill = ss304
ss304_4_cell.region = -ss304_outer_radius & +ss304_4_min & -ss304_4_max
control_universe.add_cell(ss304_4_cell)
# Create ss304 5 cell
ss304_5_cell = openmc.Cell(name='Stainless Steel 304 Cell 5')
ss304_5_cell.fill = ss304
ss304_5_cell.region = -ss304_outer_radius & +ss304_5_min & -ss304_5_max
control_universe.add_cell(ss304_5_cell)
# Create b4c 1 cell
b4c_1_cell = openmc.Cell(name='B4C cell')
b4c_1_cell.fill = b4c
b4c_1_cell.region = -b4c_outer_radius & +b4c_1_min & -b4c_1_max
control_universe.add_cell(b4c_1_cell)
# Create a clad Cell
clad_cell = openmc.Cell(name='Stainless Steel 304 Cladding')
clad_cell.fill = ss304
clad_cell.region = -clad_outer_radius & +ss304_outer_radius & +clad_min & -clad_max #Miriam: cladding is only the exterior coat.
control_universe.add_cell(clad_cell)
# Create water universe to surround the lattice
all_water_cell = openmc.Cell(fill=water)
outer_universe = openmc.Universe(cells=(all_water_cell,))
# Create surfaces that will divide rings in the circular lattice
ring_radii = np.array([0.0, 8.0, 16.0, 24.0, 32.0, 40.0])
radial_surf = [openmc.ZCylinder(r=r) for r in
(ring_radii[:-1] + ring_radii[1:])/2]
water_cells = []
for i in range(ring_radii.size):
# Create annular region
if i == 0:
water_region = -radial_surf[i]
elif i == ring_radii.size - 1:
water_region = +radial_surf[i-1]
else:
water_region = +radial_surf[i-1] & -radial_surf[i]
water_cells.append(openmc.Cell(fill=water, region=water_region))
# Plot the rings to visualize the circular lattice, without rods
plot_args = {'width': (2*24.1, 2*24.1)}
bundle_universe = openmc.Universe(cells=water_cells)
bundle_universe.plot(**plot_args)
# Arrange the pins in the circular lattice
num_pins = [1, 6, 12, 18, 24, 30]
angles = [0, 0, 0, 0, 0, 0]
controlRods = {'numPins' :[num_pins[1], num_pins[3], num_pins[5]],
'howLeftFrom3oclock':[4 , 2 , 0]}
transientRods = {'numPins' :[num_pins[3], num_pins[2]],
'howLeftFrom3oclock':[1 , 0]}
waterRods = {'numPins' :[num_pins[5], num_pins[2]],
'howLeftFrom3oclock':[1 , 8]}
def ControlRod(controlRods,n,j):
for irod in range(len(controlRods['numPins'])):
if n == controlRods['numPins'][irod] and \
j-1 == controlRods['howLeftFrom3oclock'][irod]:
return True
return False
def TransientRod(transientRods,n,j):
for irod in range(len(transientRods['numPins'])):
if n == transientRods['numPins'][irod] and \
j-1 == transientRods['howLeftFrom3oclock'][irod]:
return True
return False
def WaterRod(waterRods,n,j):
for irod in range(len(waterRods['numPins'])):
if n == waterRods['numPins'][irod] and \
j-1 == waterRods['howLeftFrom3oclock'][irod]:
return True
return False
for i, (r, n, a) in enumerate(zip(ring_radii, num_pins, angles)):
for j in range(n):
# Determine location of center of pin
theta = (a + j/n*360.) * np.pi/180.
x = r*np.cos(theta)
y = r*np.sin(theta)
pin_boundary = openmc.ZCylinder(x0=x, y0=y, r=clad_outer_radius.r)
water_cells[i].region &= +pin_boundary
# Create each fuel pin -- note that we explicitly assign an ID so
# that we can identify the pin later when looking at tallies
if ControlRod(controlRods,n,j):
print('Adding in a control rod...')
pin = openmc.Cell(fill=control_universe, region=-pin_boundary)
elif TransientRod(transientRods,n,j):
print('Adding in a transient rod...')
pin = openmc.Cell(fill=transient_universe, region=-pin_boundary)
elif WaterRod(waterRods,n,j):
print('Adding in a water rod...')
pin = openmc.Cell(fill=outer_universe, region=-pin_boundary)
else:
pin = openmc.Cell(fill=fuel_universe, region=-pin_boundary)
pin.translation = (x, y, 0)
pin.id = (i + 1)*100 + j
bundle_universe.add_cell(pin)
# Plot the rings to visualize the filled circular lattice
bundle_universe.plot(width=(100, 100), origin=[0,0,-40],
basis='xy', color_by='material',
colors={water:'blue',uzrh:'orange',
zirconium:'green',graphite:'gray',
b4c:'yellow'})
# Plotting fuel rod
fuel_universe.plot(width=(20, 150), origin=[0,0,-40], basis='yz', color_by='material', colors={ss304:'fuchsia'})
# Plotting transient rod
transient_universe.plot(width=(20, 150), origin=[0,0,-40], basis='yz', color_by='material', colors={water:'fuchsia'})
# Plotting control rod
control_universe.plot(width=(20, 150), origin=[0,0,-40], basis='yz', color_by='material', colors={ss304:'fuchsia'})
# Geometry definitions for the reactor
reactor_wall = openmc.ZCylinder(r=50.0, boundary_type='vacuum')
reactor_top = openmc.ZPlane(z0=0.0, boundary_type='vacuum')
reactor_bottom = openmc.ZPlane(z0=-114.3, boundary_type='vacuum')
reactor = openmc.Cell()
reactor.region = -reactor_wall & -reactor_top & +reactor_bottom
reactor.fill = bundle_universe
reactor_universe = openmc.Universe(cells=[reactor])
reactor_universe.plot(width=(100, 100), origin=[0,0,-40],
basis='yz', color_by='material',
colors={water:'blue',uzrh:'orange',
zirconium:'green',graphite:'gray',
b4c:'yellow'})
geometry = openmc.Geometry(reactor_universe)
geometry.export_to_xml()
# OpenMC simulation parameters
batches = 100
inactive = 10
particles = 5000
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
bounds = [-28.527375, -28.527375, -28.527375, 28.527375, 28.527375, 28.527375]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.Source(space=uniform_dist)
settings_file.export_to_xml()
openmc.run()

521
stress-test/validation/neutron_physics.py Normal file
View file

@ -0,0 +1,521 @@
#!/usr/bin/env python3
import argparse
import os
from pathlib import Path
import re
import shutil
import subprocess
import h5py
from matplotlib import pyplot as plt
import numpy as np
import openmc
from openmc.data import K_BOLTZMANN
from .utils import zaid, szax, create_library, read_results
def main():
parser = argparse.ArgumentParser()
parser.add_argument('nuclide', type=str,
help='Name of the nuclide, e.g. "U235"')
parser.add_argument('-d', '--density', type=float, default=1.,
help='Density of the material in g/cm^3')
parser.add_argument('-e', '--energy', type=float, default=1e6,
help='Energy of the source in eV')
parser.add_argument('-p', '--particles', type=int, default=100000,
help='Number of source particles')
parser.add_argument('-c', '--code', choices=['mcnp', 'serpent'],
default='mcnp',
help='Code to validate OpenMC against.')
parser.add_argument('-s', '--suffix', type=str, default='70c',
help='MCNP cross section suffix')
parser.add_argument('-x', '--xsdir', type=str, help='XSDIR directory '
'file. If specified, it will be used to locate the '
'ACE table corresponding to the given nuclide and '
'suffix, and an HDF5 library that can be used by '
'OpenMC will be created from the data.')
parser.add_argument('-t', '--thermal', type=str, help='ZAID of the '
'thermal scattering data, e.g. "grph.10t". If '
'specified, thermal scattering data will be assigned '
'to the material.')
parser.add_argument('-o', '--output-name', type=str,
help='Name used for output.')
args = parser.parse_args()
model = NeutronPhysicsModel(
args.nuclide, args.density, args.energy, args.particles, args.code,
args.suffix, args.xsdir, args.thermal, args.output_name
)
model.run()
class NeutronPhysicsModel:
"""Monoenergetic, isotropic point source in an infinite geometry.
Parameters
----------
nuclide : str
Name of the nuclide
density : float
Density of the material in g/cm^3.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Cross section suffix
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given nuclide and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
thermal : str
ZAID of the thermal scattering data. If specified, thermal scattering
data will be assigned to the material.
name : str
Name used for output.
Attributes
----------
nuclide : str
Name of the nuclide
density : float
Density of the material in g/cm^3.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Cross section suffix for MCNP
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given nuclide and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
thermal : str
ZAID of the thermal scattering data. If specified, thermal scattering
data will be assigned to the material.
name : str
Name used for output.
temperature : float
Temperature (Kelvin) of the cross section data
bins : int
Number of bins in the energy grid
batches : int
Number of batches to simulate
min_energy : float
Lower limit of energy grid (eV)
openmc_dir : pathlib.Path
Working directory for OpenMC
other_dir : pathlib.Path
Working directory for MCNP or Serpent
table_names : list of str
Names of the ACE tables used in the model
"""
def __init__(self, nuclide, density, energy, particles, code, suffix,
xsdir=None, thermal=None, name=None):
self._temperature = None
self._bins = 500
self._batches = 100
self._min_energy = 1.e-5
self._openmc_dir = None
self._other_dir = None
self.nuclide = nuclide
self.density = density
self.energy = energy
self.particles = particles
self.code = code
self.suffix = suffix
self.xsdir = xsdir
self.thermal = thermal
self.name = name
@property
def energy(self):
return self._energy
@property
def particles(self):
return self._particles
@property
def code(self):
return self._code
@property
def suffix(self):
return self._suffix
@property
def xsdir(self):
return self._xsdir
@property
def openmc_dir(self):
if self._openmc_dir is None:
self._openmc_dir = Path('openmc')
os.makedirs(self._openmc_dir, exist_ok=True)
return self._openmc_dir
@property
def other_dir(self):
if self._other_dir is None:
self._other_dir = Path(self.code)
os.makedirs(self._other_dir, exist_ok=True)
return self._other_dir
@property
def table_names(self):
table_names = [zaid(self.nuclide, self.suffix)]
if self.thermal is not None:
table_names.append(self.thermal)
return table_names
@energy.setter
def energy(self, energy):
if energy <= self._min_energy:
msg = (f'Energy {energy} eV must be above the minimum energy '
f'{self._min_energy} eV.')
raise ValueError(msg)
self._energy = energy
@particles.setter
def particles(self, particles):
if particles % self._batches != 0:
msg = (f'Number of particles {particles} must be divisible by '
f'the number of batches {self._batches}.')
raise ValueError(msg)
self._particles = particles
@code.setter
def code(self, code):
if code not in ('mcnp', 'serpent'):
msg = (f'Unsupported code {code}: code must be either "mcnp" or '
'"serpent".')
raise ValueError(msg)
executable = 'mcnp6' if code == 'mcnp' else 'sss2'
if not shutil.which(executable, os.X_OK):
msg = f'Unable to locate executable {executable} in path.'
raise ValueError(msg)
self._code = code
@suffix.setter
def suffix(self, suffix):
match = '(7[0-4]c)|(8[0-6]c)|(71[0-6]nc)|[0][3,6,9]c|[1][2,5,8]c'
if not re.match(match, suffix):
msg = f'Unsupported cross section suffix {suffix}.'
raise ValueError(msg)
self._suffix = suffix
@xsdir.setter
def xsdir(self, xsdir):
if xsdir is not None:
xsdir = Path(xsdir)
if not xsdir.is_file():
msg = f'Could not locate the XSDIR file {xsdir}.'
raise ValueError(msg)
self._xsdir = xsdir
def _make_openmc_input(self):
"""Generate the OpenMC input XML
"""
# Define material
mat = openmc.Material()
mat.add_nuclide(self.nuclide, 1.0)
if self.thermal is not None:
name, suffix = self.thermal.split('.')
thermal_name = openmc.data.thermal.get_thermal_name(name)
mat.add_s_alpha_beta(thermal_name)
mat.set_density('g/cm3', self.density)
materials = openmc.Materials([mat])
if self.xsdir is not None:
xs_path = (self.openmc_dir / 'cross_sections.xml').resolve()
materials.cross_sections = str(xs_path)
materials.export_to_xml(self.openmc_dir / 'materials.xml')
# Set up geometry
x1 = openmc.XPlane(x0=-1.e9, boundary_type='reflective')
x2 = openmc.XPlane(x0=+1.e9, boundary_type='reflective')
y1 = openmc.YPlane(y0=-1.e9, boundary_type='reflective')
y2 = openmc.YPlane(y0=+1.e9, boundary_type='reflective')
z1 = openmc.ZPlane(z0=-1.e9, boundary_type='reflective')
z2 = openmc.ZPlane(z0=+1.e9, boundary_type='reflective')
cell = openmc.Cell(fill=materials)
cell.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2
geometry = openmc.Geometry([cell])
geometry.export_to_xml(self.openmc_dir / 'geometry.xml')
# Define source
source = openmc.Source()
source.space = openmc.stats.Point((0,0,0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([self.energy], [1.])
# Settings
settings = openmc.Settings()
if self._temperature is not None:
settings.temperature = {'default': self._temperature}
settings.source = source
settings.particles = self.particles // self._batches
settings.run_mode = 'fixed source'
settings.batches = self._batches
settings.create_fission_neutrons = False
settings.export_to_xml(self.openmc_dir / 'settings.xml')
# Define tallies
energy_bins = np.logspace(np.log10(self._min_energy),
np.log10(1.0001*self.energy), self._bins+1)
energy_filter = openmc.EnergyFilter(energy_bins)
tally = openmc.Tally(name='tally')
tally.filters = [energy_filter]
tally.scores = ['flux']
tallies = openmc.Tallies([tally])
tallies.export_to_xml(self.openmc_dir / 'tallies.xml')
def _make_mcnp_input(self):
"""Generate the MCNP input file
"""
# Create the problem description
lines = ['Point source in infinite geometry']
# Create the cell cards: material 1 inside sphere, void outside
lines.append('c --- Cell cards ---')
if self._temperature is not None:
kT = self._temperature * K_BOLTZMANN * 1e-6
lines.append(f'1 1 -{self.density} -1 imp:n=1 tmp={kT}')
else:
lines.append(f'1 1 -{self.density} -1 imp:n=1')
lines.append('2 0 1 imp:n=0')
lines.append('')
# Create the surface cards: box centered on origin with 2e9 cm sides`
# and reflective boundary conditions
lines.append('c --- Surface cards ---')
lines.append('*1 rpp -1.e9 1e9 -1.e9 1.e9 -1.e9 1.e9')
lines.append('')
# Create the data cards
lines.append('c --- Data cards ---')
# Materials
if re.match('(71[0-6]nc)', self.suffix):
name = szax(self.nuclide, self.suffix)
else:
name = zaid(self.nuclide, self.suffix)
lines.append(f'm1 {name} 1.0')
if self.thermal is not None:
lines.append(f'mt1 {self.thermal}')
lines.append('nonu 2')
# Physics: neutron transport
lines.append('mode n')
# Source definition: isotropic point source at center of sphere
energy = self.energy * 1e-6
lines.append(f'sdef cel=1 erg={energy}')
# Tallies: neutron flux over cell
lines.append('f4:n 1')
min_energy = self._min_energy * 1e-6
lines.append(f'e4 {min_energy} {self._bins-1}ilog {1.0001*energy}')
# Problem termination: number of particles to transport
lines.append(f'nps {self.particles}')
# Write the problem
with open(self.other_dir / 'inp', 'w') as f:
f.write('\n'.join(lines))
def _make_serpent_input(self):
"""Generate the Serpent input file
"""
# Create the problem description
lines = ['% Point source in infinite geometry']
lines.append('')
# Set the cross section library directory
if self.xsdir is not None:
xsdata = (self.other_dir / 'xsdata').resolve()
lines.append(f'set acelib "{xsdata}"')
lines.append('')
# Create the cell cards: material 1 inside sphere, void outside
lines.append('% --- Cell cards ---')
lines.append('cell 1 0 m1 -1')
lines.append('cell 2 0 outside 1')
lines.append('')
# Create the surface cards: box centered on origin with 2e9 cm sides`
# and reflective boundary conditions
lines.append('% --- Surface cards ---')
lines.append('surf 1 cube 0.0 0.0 0.0 1.e9')
# Reflective boundary conditions
lines.append('set bc 2')
lines.append('')
# Create the material cards
lines.append('% --- Material cards ---')
name = zaid(self.nuclide, self.suffix)
if self.thermal is not None:
Z, A, m = openmc.data.zam(self.nuclide)
lines.append(f'mat m1 -{self.density} moder t1 {1000*Z + A}')
else:
lines.append(f'mat m1 -{self.density}')
lines.append(f'{name} 1.0')
# Add thermal scattering library associated with the nuclide
if self.thermal is not None:
lines.append(f'therm t1 {self.thermal}')
lines.append('')
# External source mode with isotropic point source at center of sphere
lines.append('% --- Set external source mode ---')
lines.append(f'set nps {self.particles} {self._batches}')
energy = self.energy * 1e-6
lines.append(f'src 1 n se {energy} sp 0.0 0.0 0.0')
lines.append('')
# Detector definition: flux energy spectrum
lines.append('% --- Detector definition ---')
lines.append('det 1 de 1 dc 1')
# Energy grid definition: equal lethargy spacing
min_energy = self._min_energy * 1e-6
lines.append(f'ene 1 3 {self._bins} {min_energy} {1.0001*energy}')
lines.append('')
# Treat fission as capture
lines.append('set nphys 0')
# Turn on unresolved resonance probability treatment
lines.append('set ures 1')
# Write the problem
with open(self.other_dir / 'input', 'w') as f:
f.write('\n'.join(lines))
def _plot(self):
"""Extract and plot the results
"""
# Read results
path = self.openmc_dir / f'statepoint.{self._batches}.h5'
x1, y1, _ = read_results('openmc', path)
if self.code == 'serpent':
path = self.other_dir / 'input_det0.m'
else:
path = self.other_dir / 'outp'
x2, y2, sd = read_results(self.code, path)
# Convert energies to eV
x1 *= 1e6
x2 *= 1e6
# Normalize the spectra
y1 /= np.diff(np.insert(x1, 0, self._min_energy))*sum(y1)
y2 /= np.diff(np.insert(x2, 0, self._min_energy))*sum(y2)
# Compute the relative error
err = np.zeros_like(y2)
idx = np.where(y2 > 0)
err[idx] = (y1[idx] - y2[idx])/y2[idx]
# Set up the figure
fig = plt.figure(1, facecolor='w', figsize=(8,8))
ax1 = fig.add_subplot(111)
# Create a second y-axis that shares the same x-axis, keeping the first
# axis in front
ax2 = ax1.twinx()
ax1.set_zorder(ax2.get_zorder() + 1)
ax1.patch.set_visible(False)
# Plot the spectra
label = 'Serpent' if self.code == 'serpent' else 'MCNP'
ax1.loglog(x2, y2, 'r', linewidth=1, label=label)
ax1.loglog(x1, y1, 'b', linewidth=1, label='OpenMC', linestyle='--')
# Plot the relative error and uncertainties
ax2.semilogx(x2, err, color=(0.2, 0.8, 0.0), linewidth=1)
ax2.semilogx(x2, 2*sd, color='k', linestyle='--', linewidth=1)
ax2.semilogx(x2, -2*sd, color='k', linestyle='--', linewidth=1)
# Set grid and tick marks
ax1.tick_params(axis='both', which='both', direction='in', length=10)
ax1.grid(b=False, axis='both', which='both')
ax2.tick_params(axis='y', which='both', right=False)
ax2.grid(b=True, which='both', axis='both', alpha=0.5, linestyle='--')
# Set axes labels and limits
ax1.set_xlim([self._min_energy, self.energy])
ax1.set_xlabel('Energy (eV)', size=12)
ax1.set_ylabel('Spectrum', size=12)
ax1.legend()
ax2.set_ylabel("Relative error", size=12)
title = f'{self.nuclide}'
if self.thermal is not None:
name, suffix = self.thermal.split('.')
thermal_name = openmc.data.thermal.get_thermal_name(name)
title += f' + {thermal_name}'
title += f', {self.energy:.1e} eV Source'
plt.title(title)
# Save plot
os.makedirs('plots', exist_ok=True)
if self.name is not None:
name = self.name
else:
name = f'{self.nuclide}'
if self.thermal is not None:
name += f'-{thermal_name}'
name += f'-{self.energy:.1e}eV'
if self._temperature is not None:
name += f'-{self._temperature:.1f}K'
plt.savefig(Path('plots') / f'{name}.png', bbox_inches='tight')
plt.close()
def run(self):
"""Generate inputs, run problem, and plot results.
"""
# Create HDF5 cross section library and Serpent XSDATA file
if self.xsdir is not None:
path = self.other_dir if self.code == 'serpent' else None
create_library(self.xsdir, self.table_names, self.openmc_dir, path)
# Get the temperature of the cross section data
f = h5py.File(self.openmc_dir / (self.nuclide + '.h5'), 'r')
temperature = list(f[self.nuclide]['kTs'].values())[0][()]
self._temperature = temperature / K_BOLTZMANN
# Generate input files
self._make_openmc_input()
if self.code == 'serpent':
self._make_serpent_input()
args = ['sss2', 'input']
else:
self._make_mcnp_input()
args = ['mcnp6']
if self.xsdir is not None:
args.append(f'XSDIR={self.xsdir}')
# Remove old MCNP output files
for f in ('outp', 'runtpe'):
try:
os.remove(self.other_dir / f)
except OSError:
pass
# Run code and capture and print output
p = subprocess.Popen(args, cwd=self.code, stdout=subprocess.PIPE,
stderr=subprocess.STDOUT, universal_newlines=True)
while True:
line = p.stdout.readline()
if not line and p.poll() is not None:
break
print(line, end='')
openmc.run(cwd='openmc')
self._plot()

566
stress-test/validation/photon_physics.py Normal file
View file

@ -0,0 +1,566 @@
#!/usr/bin/env python3
import argparse
import os
from pathlib import Path
import re
import shutil
import subprocess
from matplotlib import pyplot as plt
import numpy as np
import openmc
from openmc.data import ATOMIC_NUMBER, NEUTRON_MASS, K_BOLTZMANN
from .utils import create_library, read_results
def main():
parser = argparse.ArgumentParser()
parser.add_argument('element', type=str,
help='Name of the element, e.g. "U"')
parser.add_argument('-d', '--density', type=float, default=1.,
help='Density of the material in g/cm^3')
parser.add_argument('-e', '--energy', type=float, default=1e6,
help='Energy of the source in eV')
parser.add_argument('-p', '--particles', type=int, default=1000000,
help='Number of source particles')
parser.add_argument('-t', '--electron-treatment', choices=('ttb', 'led'),
default='ttb', help='Whether to use local energy'
'deposition or thick-target bremsstrahlung treatment '
'for electrons and positrons.')
parser.add_argument('-c', '--code', choices=['mcnp', 'serpent'],
default='mcnp',
help='Code to validate OpenMC against.')
parser.add_argument('-s', '--suffix', default='12p',
help='Photon cross section suffix')
parser.add_argument('-x', '--xsdir', type=str, help='XSDIR directory '
'file. If specified, it will be used to locate the '
'ACE table corresponding to the given nuclide and '
'suffix, and an HDF5 library that can be used by '
'OpenMC will be created from the data.')
parser.add_argument('-g', '--serpent_pdata', type=str, help='Directory '
'containing the additional data files needed for '
'photon physics in Serpent.')
parser.add_argument('-o', '--output-name', type=str,
help='Name used for output.')
args = parser.parse_args()
model = PhotonPhysicsModel(
args.element, args.density, [(args.element, 1.)], args.energy,
args.particles, args.electron_treatment, args.code, args.suffix,
args.xsdir, args.serpent_pdata, args.output_name
)
model.run()
class PhotonPhysicsModel:
"""Monoenergetic, isotropic point source in an infinite geometry.
Parameters
----------
material : str
Name of the material.
density : float
Density of the material in g/cm^3.
elements : list of tuple
List in which each item is a 2-tuple consisting of an element string and
the atom fraction.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
electron_treatment : {'led' or 'ttb'}
Whether to deposit electron energy locally ('led') or create secondary
bremsstrahlung photons ('ttb').
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Photon cross section suffix
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given element and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
serpent_pdata : str
Directory containing the additional data files needed for photon
physics in Serpent.
name : str
Name used for output.
Attributes
----------
material : str
Name of the material.
density : float
Density of the material in g/cm^3.
elements : list of tuple
List in which each item is a 2-tuple consisting of an element string and
the atom fraction.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
electron_treatment : {'led' or 'ttb'}
Whether to deposit electron energy locally ('led') or create secondary
bremsstrahlung photons ('ttb').
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Photon cross section suffix
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given element and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
serpent_pdata : str
Directory containing the additional data files needed for photon
physics in Serpent.
name : str
Name used for output.
bins : int
Number of bins in the energy grid
batches : int
Number of batches to simulate
cutoff_energy: float
Photon cutoff energy (eV)
openmc_dir : pathlib.Path
Working directory for OpenMC
other_dir : pathlib.Path
Working directory for MCNP or Serpent
table_names : list of str
Names of the ACE tables used in the model
"""
def __init__(self, material, density, elements, energy, particles,
electron_treatment, code, suffix, xsdir=None,
serpent_pdata=None, name=None):
self._bins = 500
self._batches = 100
self._cutoff_energy = 1.e3
self._openmc_dir = None
self._other_dir = None
self.material = material
self.density = density
self.elements = elements
self.energy = energy
self.particles = particles
self.electron_treatment = electron_treatment
self.code = code
self.suffix = suffix
self.xsdir = xsdir
self.serpent_pdata = serpent_pdata
self.name = name
@property
def energy(self):
return self._energy
@property
def particles(self):
return self._particles
@property
def code(self):
return self._code
@property
def suffix(self):
return self._suffix
@property
def xsdir(self):
return self._xsdir
@property
def serpent_pdata(self):
return self._serpent_pdata
@property
def openmc_dir(self):
if self._openmc_dir is None:
self._openmc_dir = Path('openmc')
os.makedirs(self._openmc_dir, exist_ok=True)
return self._openmc_dir
@property
def other_dir(self):
if self._other_dir is None:
self._other_dir = Path(self.code)
os.makedirs(self._other_dir, exist_ok=True)
return self._other_dir
@property
def table_names(self):
table_names = []
for element, _ in self.elements:
Z = ATOMIC_NUMBER[element]
table_names.append(f'{1000*Z}.{self.suffix}')
return table_names
@energy.setter
def energy(self, energy):
if energy <= self._cutoff_energy:
msg = (f'Energy {energy} eV must be above the cutoff energy '
f'{self._cutoff_energy} eV.')
raise ValueError(msg)
self._energy = energy
@particles.setter
def particles(self, particles):
if particles % self._batches != 0:
msg = (f'Number of particles {particles} must be divisible by '
f'the number of batches {self._batches}.')
raise ValueError(msg)
self._particles = particles
@code.setter
def code(self, code):
if code not in ('mcnp', 'serpent'):
msg = (f'Unsupported code {code}: code must be either "mcnp" or '
'"serpent".')
raise ValueError(msg)
executable = 'mcnp6' if code == 'mcnp' else 'sss2'
if not shutil.which(executable, os.X_OK):
msg = f'Unable to locate executable {executable} in path.'
raise ValueError(msg)
self._code = code
@suffix.setter
def suffix(self, suffix):
if not re.match('12p', suffix):
msg = f'Unsupported cross section suffix {suffix}.'
raise ValueError(msg)
self._suffix = suffix
@xsdir.setter
def xsdir(self, xsdir):
if xsdir is not None:
xsdir = Path(xsdir)
if not xsdir.is_file():
msg = f'Could not locate the XSDIR file {xsdir}.'
raise ValueError(msg)
self._xsdir = xsdir
@serpent_pdata.setter
def serpent_pdata(self, serpent_pdata):
if self.code == 'serpent':
if serpent_pdata is None:
msg = ('Serpent photon data path is required to run a '
'calculation with Serpent.')
raise ValueError(msg)
serpent_pdata = Path(serpent_pdata).resolve()
if not serpent_pdata.is_dir():
msg = (f'Could not locate the Serpent photon data directory '
f'{serpent_pdata}.')
raise ValueError(msg)
self._serpent_pdata = serpent_pdata
def _make_openmc_input(self):
"""Generate the OpenMC input XML
"""
# Define material
mat = openmc.Material()
for element, fraction in self.elements:
mat.add_element(element, fraction)
mat.set_density('g/cm3', self.density)
materials = openmc.Materials([mat])
if self.xsdir is not None:
xs_path = (self.openmc_dir / 'cross_sections.xml').resolve()
materials.cross_sections = str(xs_path)
materials.export_to_xml(self.openmc_dir / 'materials.xml')
# Set up geometry
x1 = openmc.XPlane(x0=-1.e9, boundary_type='reflective')
x2 = openmc.XPlane(x0=+1.e9, boundary_type='reflective')
y1 = openmc.YPlane(y0=-1.e9, boundary_type='reflective')
y2 = openmc.YPlane(y0=+1.e9, boundary_type='reflective')
z1 = openmc.ZPlane(z0=-1.e9, boundary_type='reflective')
z2 = openmc.ZPlane(z0=+1.e9, boundary_type='reflective')
cell = openmc.Cell(fill=materials)
cell.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2
geometry = openmc.Geometry([cell])
geometry.export_to_xml(self.openmc_dir / 'geometry.xml')
# Define source
source = openmc.Source()
source.space = openmc.stats.Point((0,0,0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([self.energy], [1.])
source.particle = 'photon'
# Settings
settings = openmc.Settings()
settings.source = source
settings.particles = self.particles // self._batches
settings.run_mode = 'fixed source'
settings.batches = self._batches
settings.photon_transport = True
settings.electron_treatment = self.electron_treatment
settings.cutoff = {'energy_photon' : self._cutoff_energy}
settings.export_to_xml(self.openmc_dir / 'settings.xml')
# Define tallies
energy_bins = np.logspace(np.log10(self._cutoff_energy),
np.log10(1.0001*self.energy), self._bins+1)
energy_filter = openmc.EnergyFilter(energy_bins)
particle_filter = openmc.ParticleFilter('photon')
tally = openmc.Tally(name='tally')
tally.filters = [energy_filter, particle_filter]
tally.scores = ['flux']
tallies = openmc.Tallies([tally])
tallies.export_to_xml(self.openmc_dir / 'tallies.xml')
def _make_mcnp_input(self):
"""Generate the MCNP input file
"""
# Create the problem description
lines = ['Point source in infinite geometry']
# Create the cell cards: material 1 inside sphere, void outside
lines.append('c --- Cell cards ---')
lines.append(f'1 1 -{self.density} -1 imp:p=1')
lines.append('2 0 1 imp:p=0')
lines.append('')
# Create the surface cards: box centered on origin with 2e9 cm sides`
# and reflective boundary conditions
lines.append('c --- Surface cards ---')
lines.append('*1 rpp -1.e9 1e9 -1.e9 1.e9 -1.e9 1.e9')
lines.append('')
# Create the data cards
lines.append('c --- Data cards ---')
# Materials
material_card = 'm1'
for element, fraction in self.elements:
Z = openmc.data.ATOMIC_NUMBER[element]
material_card += f' {Z}000.{self.suffix} -{fraction}'
lines.append(material_card)
# Energy in MeV
energy = self.energy * 1e-6
cutoff_energy = self._cutoff_energy * 1e-6
# Physics: photon transport, 1 keV photon cutoff energy
if self.electron_treatment == 'led':
flag = 1
else:
flag = 'j'
lines.append('mode p')
lines.append(f'phys:p j {flag} j j j')
lines.append(f'cut:p j {cutoff_energy}')
# Source definition: isotropic point source at center of sphere
lines.append(f'sdef cel=1 erg={energy}')
# Tallies: photon flux over cell
lines.append('f4:p 1')
lines.append(f'e4 {cutoff_energy} {self._bins-1}ilog {1.0001*energy}')
# Problem termination: number of particles to transport
lines.append(f'nps {self.particles}')
# Write the problem
with open(self.other_dir / 'inp', 'w') as f:
f.write('\n'.join(lines))
def _make_serpent_input(self):
"""Generate the Serpent input file
"""
# Create the problem description
lines = ['% Point source in infinite geometry']
lines.append('')
# Set the cross section library directory
if self.xsdir is not None:
xsdata = (self.other_dir / 'xsdata').resolve()
lines.append(f'set acelib "{xsdata}"')
# Set the photon data directory
lines.append(f'set pdatadir "{self.serpent_pdata}"')
lines.append('')
# Create the cell cards: material 1 inside sphere, void outside
lines.append('% --- Cell cards ---')
lines.append('cell 1 0 m1 -1')
lines.append('cell 2 0 outside 1')
lines.append('')
# Create the surface cards: box centered on origin with 2e9 cm sides`
# and reflective boundary conditions
lines.append('% --- Surface cards ---')
lines.append('surf 1 cube 0.0 0.0 0.0 1.e9')
# Reflective boundary conditions
lines.append('set bc 2')
lines.append('')
# Create the material cards
lines.append('% --- Material cards ---')
lines.append(f'mat m1 -{self.density}')
# Add element data
for element, fraction in self.elements:
Z = ATOMIC_NUMBER[element]
name = f'{1000*Z}.{self.suffix}'
lines.append(f'{name} {fraction}')
# Turn on unresolved resonance probability treatment
lines.append('set ures 1')
# Set electron treatment
if self.electron_treatment == 'led':
lines.append('set ttb 0')
else:
lines.append('set ttb 1')
# Energy in MeV
energy = self.energy * 1e-6
cutoff_energy = self._cutoff_energy * 1e-6
# Set cutoff energy
lines.append(f'set ecut 0 {cutoff_energy}')
lines.append('')
# External source mode with isotropic point source at center of sphere
lines.append('% --- Set external source mode ---')
lines.append(f'set nps {self.particles} {self._batches}')
lines.append(f'src 1 g se {energy} sp 0.0 0.0 0.0')
lines.append('')
# Detector definition: flux energy spectrum
lines.append('% --- Detector definition ---')
lines.append('det 1 de 1 dc 1')
# Energy grid definition: equal lethargy spacing
lines.append(f'ene 1 3 {self._bins} {cutoff_energy} {1.0001*energy}')
lines.append('')
# Write the problem
with open(self.other_dir / 'input', 'w') as f:
f.write('\n'.join(lines))
def _plot(self):
"""Extract and plot the results
"""
# Read results
path = self.openmc_dir / f'statepoint.{self._batches}.h5'
x1, y1, _ = read_results('openmc', path)
if self.code == 'serpent':
path = self.other_dir / 'input_det0.m'
else:
path = self.other_dir / 'outp'
x2, y2, sd = read_results(self.code, path)
# Normalize the spectra
cutoff_energy = self._cutoff_energy * 1e-6
y1 /= np.diff(np.insert(x1, 0, cutoff_energy))*sum(y1)
y2 /= np.diff(np.insert(x2, 0, cutoff_energy))*sum(y2)
# Compute the relative error
err = np.zeros_like(y2)
idx = np.where(y2 > 0)
err[idx] = (y1[idx] - y2[idx])/y2[idx]
# Set up the figure
fig = plt.figure(1, facecolor='w', figsize=(8,8))
ax1 = fig.add_subplot(111)
# Create a second y-axis that shares the same x-axis, keeping the first
# axis in front
ax2 = ax1.twinx()
ax1.set_zorder(ax2.get_zorder() + 1)
ax1.patch.set_visible(False)
# Plot the spectra
label = 'Serpent' if self.code == 'serpent' else 'MCNP'
ax1.loglog(x2, y2, 'r', linewidth=1, label=label)
ax1.loglog(x1, y1, 'b', linewidth=1, label='OpenMC', linestyle='--')
# Plot the relative error and uncertainties
ax2.semilogx(x2, err, color=(0.2, 0.8, 0.0), linewidth=1)
ax2.semilogx(x2, 2*sd, color='k', linestyle='--', linewidth=1)
ax2.semilogx(x2, -2*sd, color='k', linestyle='--', linewidth=1)
# Set grid and tick marks
ax1.tick_params(axis='both', which='both', direction='in', length=10)
ax1.grid(b=False, axis='both', which='both')
ax2.tick_params(axis='y', which='both', right=False)
ax2.grid(b=True, which='both', axis='both', alpha=0.5, linestyle='--')
# Energy in MeV
energy = self.energy * 1e-6
# Set axes labels and limits
ax1.set_xlim([cutoff_energy, energy])
ax1.set_xlabel('Energy (MeV)', size=12)
ax1.set_ylabel('Spectrum', size=12)
ax1.legend()
ax2.set_ylabel("Relative error", size=12)
title = f'{self.material}, {energy:.1e} MeV Source'
plt.title(title)
# Save plot
os.makedirs('plots', exist_ok=True)
if self.name is not None:
name = self.name
else:
name = f'{self.material}-{energy:.1e}MeV'
plt.savefig(Path('plots') / f'{name}.png', bbox_inches='tight')
plt.close()
def run(self):
"""Generate inputs, run problem, and plot results.
"""
# Create the HDF5 library
if self.xsdir is not None:
path = self.other_dir if self.code == 'serpent' else None
create_library(self.xsdir, self.table_names, self.openmc_dir, path)
# TODO: Currently the neutron libraries are still read in to OpenMC
# even when doing pure photon transport, so we need to locate them and
# register them with the library.
path = os.getenv('OPENMC_CROSS_SECTIONS')
lib = openmc.data.DataLibrary.from_xml(path)
path = self.openmc_dir / 'cross_sections.xml'
data_lib = openmc.data.DataLibrary.from_xml(path)
for element, fraction in self.elements:
element = openmc.Element(element)
for nuclide, _, _ in element.expand(fraction, 'ao'):
h5_file = lib.get_by_material(nuclide)['path']
data_lib.register_file(h5_file)
data_lib.export_to_xml(path)
# Generate input files
self._make_openmc_input()
if self.code == 'serpent':
self._make_serpent_input()
args = ['sss2', 'input']
else:
self._make_mcnp_input()
args = ['mcnp6']
if self.xsdir is not None:
args.append(f'XSDIR={self.xsdir}')
# Remove old MCNP output files
for f in ('outp', 'runtpe'):
try:
os.remove(self.other_dir / f)
except OSError:
pass
# Run code and capture and print output
p = subprocess.Popen(
args, cwd=self.other_dir, stdout=subprocess.PIPE,
stderr=subprocess.STDOUT, universal_newlines=True
)
while True:
line = p.stdout.readline()
if not line and p.poll() is not None:
break
print(line, end='')
openmc.run(cwd=self.openmc_dir)
self._plot()

639
stress-test/validation/photon_production.py Normal file
View file

@ -0,0 +1,639 @@
#!/usr/bin/env python3
import argparse
import os
from pathlib import Path
import re
import shutil
import subprocess
import h5py
from matplotlib import pyplot as plt
import numpy as np
import openmc
from openmc.data import K_BOLTZMANN, NEUTRON_MASS
from .utils import zaid, szax, create_library, read_results
def main():
parser = argparse.ArgumentParser()
parser.add_argument('nuclide', type=str,
help='Name of the nuclide, e.g. "U235"')
parser.add_argument('-d', '--density', type=float, default=1.,
help='Density of the material in g/cm^3')
parser.add_argument('-e', '--energy', type=float, default=1e6,
help='Energy of the source in eV')
parser.add_argument('-p', '--particles', type=int, default=1000000,
help='Number of source particles')
parser.add_argument('-t', '--electron-treatment', choices=('ttb', 'led'),
default='ttb', help='Whether to use local energy'
'deposition or thick-target bremsstrahlung treatment '
'for electrons and positrons.')
parser.add_argument('-c', '--code', choices=['mcnp', 'serpent'],
default='mcnp',
help='Code to validate OpenMC against.')
parser.add_argument('-s', '--suffix', default='70c',
help='Neutron cross section suffix')
parser.add_argument('-k', '--photon-suffix', default='12p',
help='Photon cross section suffix')
parser.add_argument('-x', '--xsdir', type=str, help='XSDIR directory '
'file. If specified, it will be used to locate the '
'ACE table corresponding to the given nuclide and '
'suffix, and an HDF5 library that can be used by '
'OpenMC will be created from the data.')
parser.add_argument('-g', '--serpent_pdata', type=str, help='Directory '
'containing the additional data files needed for '
'photon physics in Serpent.')
parser.add_argument('-o', '--output-name', type=str,
help='Name used for output.')
args = parser.parse_args()
model = PhotonProductionModel(
args.nuclide, args.density, [(args.nuclide, 1.)], args.energy,
args.particles, args.electron_treatment, args.code, args.suffix,
args.photon_suffix, args.xsdir, args.serpent_pdata, args.output_name
)
model.run()
class PhotonProductionModel:
"""Monoenergetic, monodirectional neutron source directed down a thin,
infinitely long cylinder ('Broomstick' problem).
Parameters
----------
material : str
Name of the material.
density : float
Density of the material in g/cm^3.
nuclides : list of tuple
List in which each item is a 2-tuple consisting of a nuclide string and
the atom fraction.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
electron_treatment : {'led' or 'ttb'}
Whether to deposit electron energy locally ('led') or create secondary
bremsstrahlung photons ('ttb').
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Neutron cross section suffix
photon_suffix : str
Photon cross section suffix
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given nuclide and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
serpent_pdata : str
Directory containing the additional data files needed for photon
physics in Serpent.
name : str
Name used for output.
Attributes
----------
material : str
Name of the material.
density : float
Density of the material in g/cm^3.
nuclides : list of tuple
List in which each item is a 2-tuple consisting of a nuclide string and
the atom fraction.
energy : float
Energy of the source (eV)
particles : int
Number of source particles.
electron_treatment : {'led' or 'ttb'}
Whether to deposit electron energy locally ('led') or create secondary
bremsstrahlung photons ('ttb').
code : {'mcnp', 'serpent'}
Code to validate against
suffix : str
Neutron cross section suffix
photon_suffix : str
Photon cross section suffix
xsdir : str
XSDIR directory file. If specified, it will be used to locate the ACE
table corresponding to the given nuclide and suffix, and an HDF5
library that can be used by OpenMC will be created from the data.
serpent_pdata : str
Directory containing the additional data files needed for photon
physics in Serpent.
name : str
Name used for output.
temperature : float
Temperature (Kelvin) of the cross section data
bins : int
Number of bins in the energy grid
batches : int
Number of batches to simulate
max_energy : float
Upper limit of energy grid (eV)
cutoff_energy: float
Photon cutoff energy (eV)
openmc_dir : pathlib.Path
Working directory for OpenMC
other_dir : pathlib.Path
Working directory for MCNP or Serpent
table_names : list of str
Names of the ACE tables used in the model
"""
def __init__(self, material, density, nuclides, energy, particles,
electron_treatment, code, suffix, photon_suffix, xsdir=None,
serpent_pdata=None, name=None):
self._temperature = None
self._bins = 500
self._batches = 100
self._cutoff_energy = 1.e3
self._openmc_dir = None
self._other_dir = None
self.material = material
self.density = density
self.nuclides = nuclides
self.energy = energy
self.particles = particles
self.electron_treatment = electron_treatment
self.code = code
self.suffix = suffix
self.photon_suffix = photon_suffix
self.xsdir = xsdir
self.serpent_pdata = serpent_pdata
self.name = name
@property
def particles(self):
return self._particles
@property
def code(self):
return self._code
@property
def suffix(self):
return self._suffix
@property
def photon_suffix(self):
return self._photon_suffix
@property
def xsdir(self):
return self._xsdir
@property
def serpent_pdata(self):
return self._serpent_pdata
@property
def max_energy(self):
if self.energy < 1.e6:
return 1.e7
else:
return self.energy * 10
@property
def openmc_dir(self):
if self._openmc_dir is None:
self._openmc_dir = Path('openmc')
os.makedirs(self._openmc_dir, exist_ok=True)
return self._openmc_dir
@property
def other_dir(self):
if self._other_dir is None:
self._other_dir = Path(self.code)
os.makedirs(self._other_dir, exist_ok=True)
return self._other_dir
@property
def table_names(self):
table_names = []
for nuclide, _ in self.nuclides:
table_names.append(zaid(nuclide, self.suffix))
Z, A, m = openmc.data.zam(nuclide)
photon_table = f'{1000*Z}.{self.photon_suffix}'
if photon_table not in table_names:
table_names.append(photon_table)
return table_names
@particles.setter
def particles(self, particles):
if particles % self._batches != 0:
msg = (f'Number of particles {particles} must be divisible by '
f'the number of batches {self._batches}.')
raise ValueError(msg)
self._particles = particles
@code.setter
def code(self, code):
if code not in ('mcnp', 'serpent'):
msg = (f'Unsupported code {code}: code must be either "mcnp" or '
'"serpent".')
raise ValueError(msg)
executable = 'mcnp6' if code == 'mcnp' else 'sss2'
if not shutil.which(executable, os.X_OK):
msg = f'Unable to locate executable {executable} in path.'
raise ValueError(msg)
self._code = code
@suffix.setter
def suffix(self, suffix):
match = '(7[0-4]c)|(8[0-6]c)|(71[0-6]nc)|[0][3,6,9]c|[1][2,5,8]c'
if not re.match(match, suffix):
msg = f'Unsupported cross section suffix {suffix}.'
raise ValueError(msg)
self._suffix = suffix
@photon_suffix.setter
def photon_suffix(self, photon_suffix):
if not re.match('12p', photon_suffix):
msg = f'Unsupported photon cross section suffix {photon_suffix}.'
raise ValueError(msg)
self._photon_suffix = photon_suffix
@xsdir.setter
def xsdir(self, xsdir):
if xsdir is not None:
xsdir = Path(xsdir)
if not xsdir.is_file():
msg = f'Could not locate the XSDIR file {xsdir}.'
raise ValueError(msg)
self._xsdir = xsdir
@serpent_pdata.setter
def serpent_pdata(self, serpent_pdata):
if self.code == 'serpent':
if serpent_pdata is None:
msg = ('Serpent photon data path is required to run a '
'calculation with Serpent.')
raise ValueError(msg)
serpent_pdata = Path(serpent_pdata).resolve()
if not serpent_pdata.is_dir():
msg = (f'Could not locate the Serpent photon data directory '
f'{serpent_pdata}.')
raise ValueError(msg)
self._serpent_pdata = serpent_pdata
def _make_openmc_input(self):
"""Generate the OpenMC input XML
"""
# Define material
mat = openmc.Material()
for nuclide, fraction in self.nuclides:
mat.add_nuclide(nuclide, fraction)
mat.set_density('g/cm3', self.density)
materials = openmc.Materials([mat])
if self.xsdir is not None:
xs_path = (self.openmc_dir / 'cross_sections.xml').resolve()
materials.cross_sections = str(xs_path)
materials.export_to_xml(self.openmc_dir / 'materials.xml')
# Instantiate surfaces
cyl = openmc.XCylinder(boundary_type='vacuum', r=1.e-6)
px1 = openmc.XPlane(boundary_type='vacuum', x0=-1.)
px2 = openmc.XPlane(boundary_type='transmission', x0=1.)
px3 = openmc.XPlane(boundary_type='vacuum', x0=1.e9)
# Instantiate cells
inner_cyl_left = openmc.Cell()
inner_cyl_right = openmc.Cell()
outer_cyl = openmc.Cell()
# Set cells regions and materials
inner_cyl_left.region = -cyl & +px1 & -px2
inner_cyl_right.region = -cyl & +px2 & -px3
outer_cyl.region = ~(-cyl & +px1 & -px3)
inner_cyl_right.fill = mat
# Create root universe and export to XML
geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
geometry.export_to_xml(self.openmc_dir / 'geometry.xml')
# Define source
source = openmc.Source()
source.space = openmc.stats.Point((0,0,0))
source.angle = openmc.stats.Monodirectional()
source.energy = openmc.stats.Discrete([self.energy], [1.])
source.particle = 'neutron'
# Settings
settings = openmc.Settings()
if self._temperature is not None:
settings.temperature = {'default': self._temperature}
settings.source = source
settings.particles = self.particles // self._batches
settings.run_mode = 'fixed source'
settings.batches = self._batches
settings.photon_transport = True
settings.electron_treatment = self.electron_treatment
settings.cutoff = {'energy_photon' : self._cutoff_energy}
settings.export_to_xml(self.openmc_dir / 'settings.xml')
# Define filters
surface_filter = openmc.SurfaceFilter(cyl)
particle_filter = openmc.ParticleFilter('photon')
energy_bins = np.logspace(np.log10(self._cutoff_energy),
np.log10(self.max_energy), self._bins+1)
energy_filter = openmc.EnergyFilter(energy_bins)
# Create tallies and export to XML
tally = openmc.Tally(name='tally')
tally.filters = [surface_filter, energy_filter, particle_filter]
tally.scores = ['current']
tallies = openmc.Tallies([tally])
tallies.export_to_xml(self.openmc_dir / 'tallies.xml')
def _make_mcnp_input(self):
"""Generate the MCNP input file
"""
# Create the problem description
lines = ['Broomstick problem']
# Create the cell cards: material 1 inside cylinder, void outside
lines.append('c --- Cell cards ---')
if self._temperature is not None:
kT = self._temperature * openmc.data.K_BOLTZMANN * 1e-6
lines.append(f'1 1 -{self.density} -4 6 -7 imp:n,p=1 tmp={kT}')
else:
lines.append(f'1 1 -{self.density} -4 6 -7 imp:n,p=1')
lines.append('2 0 -4 5 -6 imp:n,p=1')
lines.append('3 0 #(-4 5 -7) imp:n,p=0')
lines.append('')
# Create the surface cards: cylinder with radius 1e-6 cm along x-axis
lines.append('c --- Surface cards ---')
lines.append('4 cx 1.0e-6')
lines.append('5 px -1.0')
lines.append('6 px 1.0')
lines.append('7 px 1.0e9')
lines.append('')
# Create the data cards
lines.append('c --- Data cards ---')
# Materials
material_card = 'm1'
for nuclide, fraction in self.nuclides:
if re.match('(71[0-6]nc)', self.suffix):
name = szax(nuclide, self.suffix)
else:
name = zaid(nuclide, self.suffix)
material_card += f' {name} -{fraction} plib={self.photon_suffix}'
lines.append(material_card)
# Energy in MeV
energy = self.energy * 1e-6
max_energy = self.max_energy * 1e-6
cutoff_energy = self._cutoff_energy * 1e-6
# Physics: neutron and neutron-induced photon, 1 keV photon cutoff energy
if self.electron_treatment == 'led':
flag = 1
else:
flag = 'j'
lines.append('mode n p')
lines.append(f'phys:p j {flag} j j j')
lines.append(f'cut:p j {cutoff_energy}')
# Source definition: point source at origin monodirectional along
# positive x-axis
lines.append(f'sdef cel=2 erg={energy} vec=1 0 0 dir=1 par=1')
# Tallies: Photon current over surface
lines.append('f1:p 4')
lines.append(f'e1 {cutoff_energy} {self._bins-1}ilog {max_energy}')
# Problem termination: number of particles to transport
lines.append(f'nps {self.particles}')
# Write the problem
with open(self.other_dir / 'inp', 'w') as f:
f.write('\n'.join(lines))
def _make_serpent_input(self):
"""Generate the Serpent input file
"""
# Create the problem description
lines = ['% Broomstick problem']
lines.append('')
# Set the cross section library directory
if self.xsdir is not None:
xsdata = (self.other_dir / 'xsdata').resolve()
lines.append(f'set acelib "{xsdata}"')
# Set the photon data directory
lines.append(f'set pdatadir "{self.serpent_pdata}"')
lines.append('')
# Create the cell cards: material 1 inside cylinder, void outside
lines.append('% --- Cell cards ---')
lines.append('cell 1 0 m1 -1 3 -4')
lines.append('cell 2 0 void -1 2 -3')
lines.append('cell 3 0 outside 1')
lines.append('cell 4 0 outside -2')
lines.append('cell 5 0 outside 4')
lines.append('')
# Create the surface cards: cylinder with radius 1e-6 cm along x-axis
lines.append('% --- Surface cards ---')
lines.append('surf 1 cylx 0.0 0.0 1.0e-6')
lines.append('surf 2 px -1.0')
lines.append('surf 3 px 1.0')
lines.append('surf 4 px 1.0e9')
lines.append('')
# Create the material cards
lines.append('% --- Material cards ---')
lines.append(f'mat m1 -{self.density}')
elements = {}
for nuclide, fraction in self.nuclides:
# Add nuclide data
name = zaid(nuclide, self.suffix)
lines.append(f'{name} {fraction}')
# Sum element fractions
Z, A, m = openmc.data.zam(nuclide)
name = f'{1000*Z}.{self.photon_suffix}'
if name not in elements:
elements[name] = fraction
else:
elements[name] += fraction
# Add element data
for name, fraction in elements.items():
lines.append(f'{name} {fraction}')
lines.append('')
# Turn on unresolved resonance probability treatment
lines.append('set ures 1')
# Set electron treatment
if self.electron_treatment == 'led':
lines.append('set ttb 0')
else:
lines.append('set ttb 1')
# Turn on Doppler broadening of Compton scattered photons (on by
# default)
lines.append('set cdop 1')
# Coupled neutron-gamma calculations (0 is off, 1 is analog, 2 is
# implicit)
lines.append('set ngamma 1')
# Energy in MeV
energy = self.energy * 1e-6
max_energy = self.max_energy * 1e-6
cutoff_energy = self._cutoff_energy * 1e-6
# Set cutoff energy
lines.append(f'set ecut 0 {cutoff_energy}')
lines.append('')
# External source mode with isotropic point source at center of sphere
lines.append('% --- Set external source mode ---')
lines.append(f'set nps {self.particles} {self._batches}')
lines.append(f'src 1 n se {energy} sp 0.0 0.0 0.0 sd 1.0 0.0 0.0')
lines.append('')
# Detector definition: photon current over surface
lines.append('% --- Detector definition ---')
lines.append('det 1 p de 1 ds 1 1')
# Energy grid definition: equal lethargy spacing
lines.append(f'ene 1 3 {self._bins} {cutoff_energy} {max_energy}')
lines.append('')
# Write the problem
with open(self.other_dir / 'input', 'w') as f:
f.write('\n'.join(lines))
def _plot(self):
"""Extract and plot the results
"""
# Read results
path = self.openmc_dir / f'statepoint.{self._batches}.h5'
x1, y1, _ = read_results('openmc', path)
if self.code == 'serpent':
path = self.other_dir / 'input_det0.m'
else:
path = self.other_dir / 'outp'
x2, y2, sd = read_results(self.code, path)
# Normalize the spectra
cutoff_energy = self._cutoff_energy * 1e-6
y1 /= np.diff(np.insert(x1, 0, cutoff_energy))*sum(y1)
y2 /= np.diff(np.insert(x2, 0, cutoff_energy))*sum(y2)
# Compute the relative error
err = np.zeros_like(y2)
idx = np.where(y2 > 0)
err[idx] = (y1[idx] - y2[idx])/y2[idx]
# Set up the figure
fig = plt.figure(1, facecolor='w', figsize=(8,8))
ax1 = fig.add_subplot(111)
# Create a second y-axis that shares the same x-axis, keeping the first
# axis in front
ax2 = ax1.twinx()
ax1.set_zorder(ax2.get_zorder() + 1)
ax1.patch.set_visible(False)
# Plot the spectra
label = 'Serpent' if self.code == 'serpent' else 'MCNP'
ax1.loglog(x2, y2, 'r', linewidth=1, label=label)
ax1.loglog(x1, y1, 'b', linewidth=1, label='OpenMC', linestyle='--')
# Plot the relative error and uncertainties
ax2.semilogx(x2, err, color=(0.2, 0.8, 0.0), linewidth=1)
ax2.semilogx(x2, 2*sd, color='k', linestyle='--', linewidth=1)
ax2.semilogx(x2, -2*sd, color='k', linestyle='--', linewidth=1)
# Set grid and tick marks
ax1.tick_params(axis='both', which='both', direction='in', length=10)
ax1.grid(b=False, axis='both', which='both')
ax2.tick_params(axis='y', which='both', right=False)
ax2.grid(b=True, which='both', axis='both', alpha=0.5, linestyle='--')
# Energy in MeV
energy = self.energy * 1e-6
max_energy = self.max_energy * 1e-6
# Set axes labels and limits
ax1.set_xlim([cutoff_energy, max_energy])
ax1.set_xlabel('Energy (MeV)', size=12)
ax1.set_ylabel('Particle Current', size=12)
ax1.legend()
ax2.set_ylabel("Relative error", size=12)
title = f'{self.material}, {energy:.1e} MeV Source'
plt.title(title)
# Save plot
os.makedirs('plots', exist_ok=True)
if self.name is not None:
name = self.name
else:
name = f'{self.material}-{energy:.1e}MeV'
if self._temperature is not None:
name += f'-{self._temperature:.1f}K'
plt.savefig(Path('plots') / f'{name}.png', bbox_inches='tight')
plt.close()
def run(self):
"""Generate inputs, run problem, and plot results.
"""
# Create HDF5 cross section library and Serpent XSDATA file
if self.xsdir is not None:
path = self.other_dir if self.code == 'serpent' else None
create_library(self.xsdir, self.table_names, self.openmc_dir, path)
# Get the temperature of the cross section data
nuclide = self.nuclides[0][0]
f = h5py.File(self.openmc_dir / (nuclide + '.h5'), 'r')
temperature = list(f[nuclide]['kTs'].values())[0][()]
self._temperature = temperature / K_BOLTZMANN
# Generate input files
self._make_openmc_input()
if self.code == 'serpent':
self._make_serpent_input()
args = ['sss2', 'input']
else:
self._make_mcnp_input()
args = ['mcnp6']
if self.xsdir is not None:
args.append(f'XSDIR={self.xsdir}')
# Remove old MCNP output files
for f in ('outp', 'runtpe'):
try:
os.remove(self.other_dir / f)
except OSError:
pass
# Run code and capture and print output
p = subprocess.Popen(args, cwd=self.code, stdout=subprocess.PIPE,
stderr=subprocess.STDOUT, universal_newlines=True)
while True:
line = p.stdout.readline()
if not line and p.poll() is not None:
break
print(line, end='')
openmc.run(cwd='openmc')
self._plot()

487
stress-test/validation/utils.py Normal file
View file

@ -0,0 +1,487 @@
from pathlib import Path
import re
import numpy as np
import openmc.data
from openmc.data import K_BOLTZMANN, NEUTRON_MASS
class XSDIR:
"""XSDIR directory file
Parameters
----------
filename : str
Path of the XSDIR file to load.
Attributes
----------
filename : str
Path of the XSDIR file.
datapath : str
Directory where the data libraries are stored.
atomic_weight_ratio : dict of int to double
Dictionary whose keys are ZAIDs and values are atomic weight ratios.
directory : dict of str to XSDIRTable
Dictionary whose keys are table names and values the entries in an
XSDIR cross section table description.
"""
def __init__(self, filename):
self.filename = filename
self.datapath = None
self.atomic_weight_ratio = {}
self.directory = {}
self._read()
def _read(self):
"""Read the XSDIR directory file.
"""
with open(self.filename) as f:
# First section: read the datapath if it is specified
line = f.readline()
tokens = re.split('\s|=', line)
if tokens[0].lower() == 'datapath':
self.datapath = tokens[1]
line = f.readline()
while line.strip().lower() != 'atomic weight ratios':
line = f.readline()
# Second section: read the ZAID/atomic weight ratio pairs
line = f.readline()
while line.strip().lower() != 'directory':
tokens = line.split()
if len(tokens) > 1:
items = {int(tokens[i]): float(tokens[i+1])
for i in range(0, len(tokens), 2)}
self.atomic_weight_ratio.update(items)
line = f.readline()
# Third section: read the available data tables
line = f.readline()
while line:
# Handle continuation lines
while line[-2] == '+':
line += f.readline()
line = line.replace('+\n', '')
# Store the entry if we need this table
tokens = line.split()
self.directory[tokens[0]] = XSDIRTable(line)
line = f.readline()
def export_to_xsdata(self, path='xsdata', table_names=None):
"""Create a Serpent XSDATA directory file.
Parameters
----------
path : str
Path to file to write. Defaults to 'xsdata'.
table_names : None, str, or iterable, optional
Tables from the XSDIR file to write to the XSDATA file. If None,
all of the entries are written. If str or iterable, only the
entries matching the table names are written.
"""
if table_names is None:
table_names = self.directory.keys()
else:
table_names = set(table_names)
# Classes of data included in the XSDATA file (continuous-energy
# neutron, neutron dosimetry, thermal scattering, and continuous-energy
# photoatomic)
data_classes = {'c': 1, 'y': 2, 't': 3, 'p': 5}
lines = []
for name in table_names:
table = self.directory.get(name)
if table is None:
msg = f'Could not find table {name} in {self.filename}.'
raise ValueError(msg)
# Check file format
if table.file_type != 'ascii':
msg = f'Unsupported file type {table.file_type} for {name}.'
raise ValueError(msg)
if self.datapath is None:
# Set the access route as the datapath if it is specified;
# otherwise, set the parent directory of XSDIR as the datapath
if table.access_route is not None:
datapath = Path(table.access_route)
else:
datapath = Path(self.filename).parent
else:
datapath = Path(self.datapath)
# Get the full path to the ace library
ace_path = datapath / table.file_name
if not ace_path.is_file():
raise ValueError(f'Could not find ACE file {ace_path}.')
zaid, suffix = name.split('.')
# Skip this table if it is not one of the data classes included in
# XSDATA
if suffix[-1] not in data_classes:
continue
# Get information about material and type of cross section data
data_class = data_classes[suffix[-1]]
if data_class == 3:
ZA = 0
m = 0
else:
zaid = int(zaid)
_, element, Z, A, m = openmc.data.get_metadata(zaid, 'nndc')
ZA = 1000*Z + A
alias = f'{element}-'
if A == 0:
alias += 'nat.'
elif m == 0:
alias += f'{A}.'
else:
alias += f'{A}m.'
alias += suffix
# Calculate the atomic weight
if zaid in self.atomic_weight_ratio:
atomic_weight = self.atomic_weight_ratio[zaid] * NEUTRON_MASS
else:
atomic_weight = table.atomic_weight_ratio * NEUTRON_MASS
# Calculate the temperature in Kelvin
temperature = table.temperature / K_BOLTZMANN * 1e6
# Entry in the XSDATA file
lines.append(f'{name} {name} {data_class} {ZA} {m} '
f'{atomic_weight:.8f} {temperature:.1f} 0 {ace_path}')
# Also write an entry with the alias if this is not a thermal
# scattering table
if data_class != 3:
lines.append(f'{alias} {name} {data_class} {ZA} {m} '
f'{atomic_weight:.8f} {temperature:.1f} 0 {ace_path}')
# Write the XSDATA file
with open(path, 'w') as f:
f.write('\n'.join(lines))
def get_tables(self, table_names):
"""Read ACE cross section tables from an XSDIR directory file.
Parameters
----------
table_names : str or iterable
Names of the ACE tables to load
Returns
-------
list of openmc.data.ace.Table
ACE cross section tables
"""
if isinstance(table_names, str):
table_names = [table_names]
else:
table_names = set(table_names)
tables = []
for name in table_names:
table = self.directory.get(name)
if table is None:
msg = f'Could not find table {name} in {self.filename}.'
raise ValueError(msg)
if self.datapath is None:
# Set the access route as the datapath if it is specified;
# otherwise, set the parent directory of XSDIR as the datapath
if table.access_route is not None:
datapath = Path(table.access_route)
else:
datapath = Path(self.filename).parent
else:
datapath = Path(self.datapath)
# Get the full path to the ace library
ace_path = datapath / table.file_name
if not ace_path.is_file():
raise ValueError(f'Could not find ACE file {ace_path}.')
zaid, suffix = name.split('.')
if re.match('(8[0-6]c)|(71[0-6]nc)', suffix):
nuclide, _, _, _, _ = openmc.data.get_metadata(int(zaid))
name = szax(nuclide, suffix)
# Get the ACE table
print(f'Converting table {name} from library {ace_path}...')
tables.append(openmc.data.ace.get_table(ace_path, name))
return tables
class XSDIRTable:
"""XSDIR description of a cross section table
Parameters
----------
line : str
Cross section table description from an XSDIR directory file.
Attributes
----------
name : str
ZAID of the table.
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
file_name : str
Name of the library that contains the table.
access_route : str
Path to the library.
file_type : {'ascii', 'binary'}
File format.
address : int
For type 1 files the address is the line number in the file where the
table starts. For type 2 files it is the record number of the first
record of the table.
table_length : int
Length (total number of words) of the table.
record_length : int
For type 1 files the record length is unused. For type 2 files it is a
multiple of the number of entries per record.
entries_per_record : int
For type 1 files this is unused. For type 2 files it is the number of
entries per record.
temperature : float
Temperature in MeV at which a neutron table is processed. This is used
only for neutron data.
ptables : bool
If true, it indicates a continuous-energy neutron nuclide has
unresolved resonance range probability tables.
"""
def __init__(self, line):
entries = line.split()
num_entries = len(entries)
self.name = entries[0]
self.atomic_weight_ratio = float(entries[1])
self.file_name = entries[2]
if entries[3] != '0':
self.access_route = entries[3]
else:
self.access_route = None
if entries[4] == '1':
self.file_type = 'ascii'
else:
self.file_type = 'binary'
self.address = int(entries[5])
self.table_length = int(entries[6])
if num_entries > 7:
self.record_length = int(entries[7])
else:
self.record_length = 0
if num_entries > 8:
self.entries_per_record = int(entries[8])
else:
self.entries_per_record = 0
if num_entries > 9:
self.temperature = float(entries[9])
else:
self.temperature = 0.0
if num_entries > 10:
self.ptables = entries[10].lower() == 'ptable'
else:
self.ptables = False
def zaid(nuclide, suffix):
"""Return ZAID for a given nuclide and cross section suffix.
Parameters
----------
nuclide : str
Name of the nuclide
suffix : str
Cross section suffix for MCNP
Returns
-------
str
ZA identifier
"""
Z, A, m = openmc.data.zam(nuclide)
# Serpent metastable convention
if re.match('[0][3,6,9]c|[1][2,5,8]c', suffix):
# Increase mass number above 300
if m > 0:
while A < 300:
A += 100
# MCNP metastable convention
else:
# Correct the ground state and first excited state of Am242, which
# are the reverse of the convention
if A == 242 and m == 0:
m = 1
elif A == 242 and m == 1:
m = 0
if m > 0:
A += 300 + 100*m
if re.match('(71[0-6]nc)', suffix):
suffix = f'8{suffix[2]}c'
return f'{1000*Z + A}.{suffix}'
def szax(nuclide, suffix):
"""Return SZAX for a given nuclide and cross section suffix.
Parameters
----------
nuclide : str
Name of the nuclide
suffix : str
Cross section suffix for MCNP
Returns
-------
str
SZA identifier
"""
Z, A, m = openmc.data.zam(nuclide)
# Correct the ground state and first excited state of Am242, which are
# the reverse of the convention
if A == 242 and m == 0:
m = 1
elif A == 242 and m == 1:
m = 0
if re.match('(7[0-4]c)|(8[0-6]c)', suffix):
suffix = f'71{suffix[1]}nc'
return f'{1000000*m + 1000*Z + A}.{suffix}'
def create_library(xsdir, table_names, hdf5_dir, xsdata_dir=None):
"""Convert the ACE data from the MCNP or Serpent distribution into an
HDF5 library that can be used by OpenMC and create and XSDATA directory
file for use with Serpent.
Parameters
----------
xsdir : str
Path of the XSDIR directory file
table_names : str or iterable
Names of the ACE tables to convert
hdf5_dir : str
Directory to write the HDF5 library to
xsdata_dir : str
If specified, an XSDATA directory file containing entries for each of
the table names provided will be written to this directory.
"""
# Create data library
data_lib = openmc.data.DataLibrary()
# Load the XSDIR directory file
xsdir = XSDIR(xsdir)
# Get the ACE cross section tables
tables = xsdir.get_tables(table_names)
for table in tables:
zaid, suffix = table.name.split('.')
# Convert cross section data
if suffix[-1] == 'c':
match = '(7[0-4]c)|(8[0-6]c)|(71[0-6]nc)'
scheme = 'mcnp' if re.match(match, suffix) else 'nndc'
data = openmc.data.IncidentNeutron.from_ace(table, scheme)
elif suffix[-1] == 'p':
data = openmc.data.IncidentPhoton.from_ace(table)
elif suffix[-1] == 't':
data = openmc.data.ThermalScattering.from_ace(table)
else:
msg = ('Unknown data class: cannot convert cross section data '
f'from table {table.name}')
raise ValueError(msg)
# Export HDF5 files and register with library
h5_file = Path(hdf5_dir) / f'{data.name}.h5'
data.export_to_hdf5(h5_file, 'w')
data_lib.register_file(h5_file)
# Write cross_sections.xml
data_lib.export_to_xml(Path(hdf5_dir) / 'cross_sections.xml')
# Write the Serpent XSDATA file
if xsdata_dir is not None:
xsdir.export_to_xsdata(Path(xsdata_dir) / 'xsdata', table_names)
def read_results(code, filename):
"""Read the energy, mean, and standard deviation from the output
Parameters
----------
code : {'openmc', 'mcnp', 'serpent'}
Code which produced the output file
filename : str
Path to the output file
Returns
-------
energy : numpy.ndarray
Energy bin values [MeV]
mean : numpy.ndarray
Sample mean of the tally
std_dev : numpy.ndarray
Sample standard deviation of the tally
"""
if code == 'openmc':
with openmc.StatePoint(filename) as sp:
t = sp.get_tally(name='tally')
energy = t.find_filter(openmc.EnergyFilter).bins[:,1]*1e-6
mean = t.mean[:,0,0]
std_dev = t.std_dev[:,0,0]
elif code == 'mcnp':
with open(filename, 'r') as f:
text = f.read()
p = text.find('1tally')
p = text.find('energy', p) + 10
q = text.find('total', p)
t = np.fromiter(text[p:q].split(), float)
t.shape = (len(t) // 3, 3)
energy = t[1:,0]
mean = t[1:,1]
std_dev = t[1:,2]
elif code == 'serpent':
with open(filename, 'r') as f:
text = re.split('\[|\]', f.read())
t = np.fromiter(text[1].split(), float)
t = t.reshape(len(t) // 12, 12)
e = np.fromiter(text[3].split(), float)
e = e.reshape(len(e) // 3, 3)
energy = e[:,1]
mean = t[:,10]
std_dev = t[:,11]
return energy, mean, std_dev