From 05285cd188576659fb2326469cd4184c0423d210 Mon Sep 17 00:00:00 2001 From: Adam Parler Date: Sat, 16 Dec 2023 19:41:24 -0800 Subject: [PATCH] Added new stress test, 'validation' --- TRIGA/model_creation.py | 509 ++++++++++++++++ stress-test/validation/neutron_physics.py | 521 ++++++++++++++++ stress-test/validation/photon_physics.py | 566 +++++++++++++++++ stress-test/validation/photon_production.py | 639 ++++++++++++++++++++ stress-test/validation/utils.py | 487 +++++++++++++++ 5 files changed, 2722 insertions(+) create mode 100644 TRIGA/model_creation.py create mode 100644 stress-test/validation/neutron_physics.py create mode 100644 stress-test/validation/photon_physics.py create mode 100644 stress-test/validation/photon_production.py create mode 100644 stress-test/validation/utils.py diff --git a/TRIGA/model_creation.py b/TRIGA/model_creation.py new file mode 100644 index 0000000..ed5eef7 --- /dev/null +++ b/TRIGA/model_creation.py @@ -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() \ No newline at end of file diff --git a/stress-test/validation/neutron_physics.py b/stress-test/validation/neutron_physics.py new file mode 100644 index 0000000..925812d --- /dev/null +++ b/stress-test/validation/neutron_physics.py @@ -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() diff --git a/stress-test/validation/photon_physics.py b/stress-test/validation/photon_physics.py new file mode 100644 index 0000000..943b7c1 --- /dev/null +++ b/stress-test/validation/photon_physics.py @@ -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() diff --git a/stress-test/validation/photon_production.py b/stress-test/validation/photon_production.py new file mode 100644 index 0000000..1faa51f --- /dev/null +++ b/stress-test/validation/photon_production.py @@ -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() diff --git a/stress-test/validation/utils.py b/stress-test/validation/utils.py new file mode 100644 index 0000000..703651f --- /dev/null +++ b/stress-test/validation/utils.py @@ -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