crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #2 from paulromano/photon-from-ace

Browse source 
Add ability to read photon/relaxation data from ACE files
This commit is contained in:
Amanda Lund 2018-03-16 13:11:51 -07:00 • committed by GitHub
commit 32790b9df2
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 493 additions and 128 deletions
Download patch file Download diff file Expand all files Collapse all files

100
docs/source/io_formats/nuclear_data.rst
View file

@ -1,8 +1,8 @@
.. _io_nuclear_data:
========================
Nuclear Data File Format
========================
=========================
Nuclear Data File Formats
=========================
---------------------
Incident Neutron Data
@ -10,7 +10,7 @@ Incident Neutron Data
**/**
:Attributes:
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
- **version** (*int[2]*) -- Major and minor version of the data
**/<nuclide name>/**
@ -22,7 +22,9 @@ Incident Neutron Data
- **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses
- **n_reaction** (*int*) -- Number of reactions
:Datasets: - **energy** (*double[]*) -- Energy points at which cross sections are tabulated
:Datasets:
- **energy** (*double[]*) -- Energies in [eV] at which cross sections
are tabulated
**/<nuclide name>/kTs/**
@ -31,7 +33,7 @@ temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Datasets:
- **<TTT>K** (*double*) -- kT values (in eV) for each temperature
- **<TTT>K** (*double*) -- kT values in [eV] for each temperature
TTT (in Kelvin)
**/<nuclide name>/reactions/reaction_<mt>/**
@ -113,6 +115,92 @@ temperature-dependent data set. For example, the data set corresponding to
:ref:`tabulated <1d_tabulated>`) -- The recoverable fission Q-value
(Q_prompt + delayed neutrons + delayed photons + betas)
--------------------
Incident Photon Data
--------------------
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
- **version** (*int[2]*) -- Major and minor version of the data
**/<element>/**
:Attributes: - **Z** (*int*) -- Atomic number
:Datasets:
- **energy** (*double[]*) -- Energies in [eV] at which cross sections
are tabulated
**/<element>/bremsstrahlung/**
:Datasets: - **electron_energy** (*double[]*) -- Incident electron energy in [eV]
- **photon_energy** (*double[]*) -- Outgoing photon energy as
fraction of incident electron energy
- **dcs** (*double[][]*) -- Bremsstrahlung differential cross section
at each incident energy in [mb/eV]
**/<element>/coherent/**
:Datasets: - **xs** (*double[]*) -- Coherent scattering cross section in [b]
- **integrated_scattering_factor** (:ref:`tabulated <1d_tabulated>`)
-- Integrated coherent scattering form factor
- **anomalous_real** (:ref:`tabulated <1d_tabulated>`) -- Real part
of the anomalous scattering factor
- **anomalous_imag** (:ref:`tabulated <1d_tabulated>`) -- Imaginary
part of the anomalous scattering factor
**/<element>/compton_profiles/**
:Datasets: - **binding_energy** (*double[]*) -- Binding energy for each subshell in [eV]
- **num_electrons** (*double[]*) -- Number of electrons in each subshell
- **pz** (*double[]*) -- Projection of the electron momentum on the
scattering vector in units of :math:`me^2 / \hbar` where :math:`m`
is the electron rest mass and :math:`e` is the electron charge
- **J** (*double[][]*) -- Compton profile for each subshell in units
of :math:`\hbar / (me^2)`
**/<element>/incoherent/**
:Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b]
- **scattering_factor** (:ref:`tabulated <1d_tabulated>`) --
**/<element>/pair_production_electron/**
:Datasets: - **xs** (*double[]*) -- Pair production (electron field) cross section in [b]
**/<element>/pair_production_nuclear/**
:Datasets: - **xs** (*double[]*) -- Pair production (nuclear field) cross section in [b]
**/<element>/photoelectric/**
:Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b]
**/<element>/stopping_powers/**
:Datasets: - **density_effect** (*double[]*) -- Density effect parameter
- **energy** (*double[]*) -- Energies in [eV]
- **s_collision** (*double[]*) -- Collisiong stopping power in [eV-cm\ :sup:`2`\ /g]
- **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g]
**/<element>/subshells/**
:Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2'
**/<element>/subshells/<designator>/**
:Attributes: - **binding_energy** (*double*) -- Binding energy of the subshell in [eV]
- **num_electrons** (*double*) -- Number of electrons in the subshell
:Datasets: - **transitions** (*double[][]*) -- Atomic relaxation data
- **xs** (*double[]*) -- Photoionization cross section for subshell
in [b] tabulated against the main energy grid
:Attributes:
- **threshold_idx** (*int*) -- Index on the energy
grid that the reaction threshold
-------------------------------
Thermal Neutron Scattering Data
-------------------------------

4
openmc/data/__init__.py
View file

@ -12,10 +12,10 @@ from .neutron import *
from .photon import *
from .decay import *
from .reaction import *
from .ace import *
from . import ace
from .angle_distribution import *
from .function import *
from .endf import *
from . import endf
from .energy_distribution import *
from .product import *
from .angle_energy import *

75
openmc/data/ace.py
View file

@ -16,13 +16,84 @@ generates ACE-format cross sections.
"""
from os import SEEK_CUR
from pathlib import PurePath
import struct
import sys
import numpy as np
from openmc.mixin import EqualityMixin
from openmc.data.endf import ENDF_FLOAT_RE
import openmc.checkvalue as cv
from .data import ATOMIC_SYMBOL
from .endf import ENDF_FLOAT_RE
def get_metadata(zaid, metastable_scheme='nndc'):
"""Return basic identifying data for a nuclide with a given ZAID.
Parameters
----------
zaid : int
ZAID (1000*Z + A) obtained from a library
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
Returns
-------
name : str
Name of the table
element : str
The atomic symbol of the isotope in the table; e.g., Zr.
Z : int
Number of protons in the nucleus
mass_number : int
Number of nucleons in the nucleus
metastable : int
Metastable state of the nucleus. A value of zero indicates ground state.
"""
cv.check_type('zaid', zaid, int)
cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp'])
Z = zaid // 1000
mass_number = zaid % 1000
if metastable_scheme == 'mcnp':
if zaid > 1000000:
# New SZA format
Z = Z % 1000
if zaid == 1095242:
metastable = 0
else:
metastable = zaid // 1000000
else:
if zaid == 95242:
metastable = 1
elif zaid == 95642:
metastable = 0
else:
metastable = 1 if mass_number > 300 else 0
elif metastable_scheme == 'nndc':
metastable = 1 if mass_number > 300 else 0
while mass_number > 3 * Z:
mass_number -= 100
# Determine name
element = ATOMIC_SYMBOL[Z]
name = '{}{}'.format(element, mass_number)
if metastable > 0:
name += '_m{}'.format(metastable)
return (name, element, Z, mass_number, metastable)
def ascii_to_binary(ascii_file, binary_file):
"""Convert an ACE file in ASCII format (type 1) to binary format (type 2).
@ -160,7 +231,7 @@ class Library(EqualityMixin):
# Determine whether file is ASCII or binary
try:
fh = open(filename, 'rb')
fh = open(str(filename), 'rb')
# Grab 10 lines of the library
sb = b''.join([fh.readline() for i in range(10)])

71
openmc/data/neutron.py
View file

@ -14,7 +14,7 @@ import numpy as np
import h5py
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Library, Table, get_table
from .ace import Library, Table, get_table, get_metadata
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record
from .fission_energy import FissionEnergyRelease
@ -33,73 +33,6 @@ from openmc.mixin import EqualityMixin
_RESONANCE_ENERGY_GRID = np.logspace(-3, 3, 61)
def _get_metadata(zaid, metastable_scheme='nndc'):
"""Return basic identifying data for a nuclide with a given ZAID.
Parameters
----------
zaid : int
ZAID (1000*Z + A) obtained from a library
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
Returns
-------
name : str
Name of the table
element : str
The atomic symbol of the isotope in the table; e.g., Zr.
Z : int
Number of protons in the nucleus
mass_number : int
Number of nucleons in the nucleus
metastable : int
Metastable state of the nucleus. A value of zero indicates ground state.
"""
cv.check_type('zaid', zaid, int)
cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp'])
Z = zaid // 1000
mass_number = zaid % 1000
if metastable_scheme == 'mcnp':
if zaid > 1000000:
# New SZA format
Z = Z % 1000
if zaid == 1095242:
metastable = 0
else:
metastable = zaid // 1000000
else:
if zaid == 95242:
metastable = 1
elif zaid == 95642:
metastable = 0
else:
metastable = 1 if mass_number > 300 else 0
elif metastable_scheme == 'nndc':
metastable = 1 if mass_number > 300 else 0
while mass_number > 3 * Z:
mass_number -= 100
# Determine name
element = ATOMIC_SYMBOL[Z]
name = '{}{}'.format(element, mass_number)
if metastable > 0:
name += '_m{}'.format(metastable)
return (name, element, Z, mass_number, metastable)
class IncidentNeutron(EqualityMixin):
"""Continuous-energy neutron interaction data.
@ -676,7 +609,7 @@ class IncidentNeutron(EqualityMixin):
# If mass number hasn't been specified, make an educated guess
zaid, xs = ace.name.split('.')
name, element, Z, mass_number, metastable = \
_get_metadata(int(zaid), metastable_scheme)
get_metadata(int(zaid), metastable_scheme)
# Assign temperature to the running list
kTs = [ace.temperature*EV_PER_MEV]

284
openmc/data/photon.py
View file

@ -12,6 +12,7 @@ from scipy.interpolate import CubicSpline
from openmc.mixin import EqualityMixin
import openmc.checkvalue as cv
from . import HDF5_VERSION
from .ace import Table, get_metadata, get_table
from .data import ATOMIC_SYMBOL, EV_PER_MEV
from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
from .function import Tabulated1D
@ -23,6 +24,15 @@ _SUBSHELLS = ['K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11',
'Q1', 'Q2', 'Q3']
# Helper function to map designator to subshell string or None
def _subshell(i):
if i == 0:
return None
else:
return _SUBSHELLS[i - 1]
_REACTION_NAME = {
501: 'Total photon interaction',
502: 'Photon coherent scattering',
@ -96,6 +106,7 @@ _STOPPING_POWERS = {}
# for each element are in a 2D array with shape (n, k) stored on the key 'Z'.
_BREMSSTRAHLUNG = {}
class AtomicRelaxation(EqualityMixin):
"""Atomic relaxation data.
@ -194,6 +205,65 @@ class AtomicRelaxation(EqualityMixin):
cv.check_type('transitions', df, pd.DataFrame)
self._transitions = transitions
@classmethod
def from_ace(cls, ace):
"""Generate atomic relaxation data from an ACE file
Parameters
----------
ace : openmc.data.ace.Table
ACE table to read from
Returns
-------
openmc.data.AtomicRelaxation
Atomic relaxation data
"""
# Create data dictionaries
binding_energy = {}
num_electrons = {}
transitions = {}
# Get shell designators
n = ace.nxs[7]
idx = ace.jxs[11]
shells = [_subshell(int(i)) for i in ace.xss[idx : idx+n]]
# Get number of electrons for each shell
idx = ace.jxs[12]
for shell, num in zip(shells, ace.xss[idx : idx+n]):
num_electrons[shell] = num
# Get binding energy for each shell
idx = ace.jxs[13]
for shell, e in zip(shells, ace.xss[idx : idx+n]):
binding_energy[shell] = e*EV_PER_MEV
# Get transition table
columns = ['secondary', 'tertiary', 'energy (eV)', 'probability']
idx = ace.jxs[18]
for i, subi in enumerate(shells):
n_transitions = int(ace.xss[ace.jxs[15] + i])
if n_transitions > 0:
records = []
for j in range(n_transitions):
subj = _subshell(int(ace.xss[idx]))
subk = _subshell(int(ace.xss[idx + 1]))
etr = ace.xss[idx + 2]*EV_PER_MEV
if j == 0:
ftr = ace.xss[idx + 3]
else:
ftr = ace.xss[idx + 3] - ace.xss[idx - 1]
records.append((subj, subk, etr, ftr))
idx += 4
# Create dataframe for transitions
transitions[subi] = pd.DataFrame.from_records(
records, columns=columns)
return cls(binding_energy, num_electrons, transitions)
@classmethod
def from_endf(cls, ev_or_filename):
"""Generate atomic relaxation data from an ENDF evaluation
@ -225,13 +295,6 @@ class AtomicRelaxation(EqualityMixin):
params = get_head_record(file_obj)
n_subshells = params[4]
# Helper function to map designator to subshell string or None
def subshell(i):
if i == 0:
return None
else:
return _SUBSHELLS[i - 1]
# Create data dictionaries
binding_energy = {}
num_electrons = {}
@ -241,7 +304,7 @@ class AtomicRelaxation(EqualityMixin):
# Read data for each subshell
for i in range(n_subshells):
params, list_items = get_list_record(file_obj)
subi = subshell(int(params[0]))
subi = _subshell(int(params[0]))
n_transitions = int(params[5])
binding_energy[subi] = list_items[0]
num_electrons[subi] = list_items[1]
@ -250,8 +313,8 @@ class AtomicRelaxation(EqualityMixin):
# Read transition data
records = []
for j in range(n_transitions):
subj = subshell(int(list_items[6*(j+1)]))
subk = subshell(int(list_items[6*(j+1) + 1]))
subj = _subshell(int(list_items[6*(j+1)]))
subk = _subshell(int(list_items[6*(j+1) + 1]))
etr = list_items[6*(j+1) + 2]
ftr = list_items[6*(j+1) + 3]
records.append((subj, subk, etr, ftr))
@ -261,9 +324,7 @@ class AtomicRelaxation(EqualityMixin):
records, columns=columns)
# Return instance of class
data = cls(binding_energy, num_electrons, transitions)
return data
return cls(binding_energy, num_electrons, transitions)
def to_hdf5(self, group):
raise NotImplementedError
@ -272,7 +333,7 @@ class AtomicRelaxation(EqualityMixin):
class IncidentPhoton(EqualityMixin):
"""Photon interaction data.
This class stores photo-atomic, photo-nuclear, atomic relaxation,
This class stores photo-atomic, photo-nuclear, atomic relaxation,
Compton profile, stopping power, and bremsstrahlung data assembled from
different sources. To create an instance, the factory method
:meth:`IncidentPhoton.from_endf` can be used. To add atomic relaxation or
@ -369,6 +430,102 @@ class IncidentPhoton(EqualityMixin):
AtomicRelaxation)
self._atomic_relaxation = atomic_relaxation
@classmethod
def from_ace(cls, ace_or_filename):
"""Generate incident photon data from an ACE table
Parameters
----------
ace_or_filename : str or openmc.data.ace.Table
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
Returns
-------
openmc.data.IncidentPhoton
Photon interaction data
"""
# First obtain the data for the first provided ACE table/file
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
ace = get_table(ace_or_filename)
# Get atomic number based on name of ACE table
zaid = ace.name.split('.')[0]
Z = get_metadata(int(zaid))[2]
# Read each reaction
data = cls(Z)
for mt in (502, 504, 515, 522):
data.reactions[mt] = PhotonReaction.from_ace(ace, mt)
# Compton profiles
n_shell = ace.nxs[5]
if n_shell != 0:
# Get number of electrons in each shell
idx = ace.jxs[6]
data.compton_profiles['num_electrons'] = ace.xss[idx : idx+n_shell]
# Get binding energy for each shell
idx = ace.jxs[7]
data.compton_profiles['binding_energy'] = ace.xss[idx : idx+n_shell]
# Create Compton profile for each electron shell
profiles = []
for k in range(n_shell):
# Get number of momentum values and interpolation scheme
loca = int(ace.xss[ace.jxs[9] + k])
jj = int(ace.xss[ace.jxs[10] + loca - 1])
m = int(ace.xss[ace.jxs[10] + loca])
# Read momentum and PDF
idx = ace.jxs[10] + loca + 1
pz = ace.xss[idx : idx+m]
pdf = ace.xss[idx+m : idx+2*m]
# Create proflie function
J_k = Tabulated1D(pz, pdf, [m], [jj])
profiles.append(J_k)
data.compton_profiles['J'] = profiles
# Subshell photoelectric xs and atomic relaxation data
if ace.nxs[7] > 0:
data.atomic_relaxation = AtomicRelaxation.from_ace(ace)
# Get subshell designators
n_subshells = ace.nxs[7]
idx = ace.jxs[11]
designators = [int(i) for i in ace.xss[idx : idx+n_subshells]]
# Get energy grid for subshell photoionization
n_energy = ace.nxs[3]
idx = ace.jxs[1]
energy = np.exp(ace.xss[idx : idx+n_energy])*EV_PER_MEV
# Get cross section for each subshell
idx = ace.jxs[16]
for d in designators:
# Create photon reaction
mt = 533 + d
rx = PhotonReaction(mt)
data.reactions[mt] = rx
# Store cross section
xs = ace.xss[idx : idx+n_energy].copy()
nonzero = (xs != 0.0)
xs[nonzero] = np.exp(xs[nonzero])
rx.xs = Tabulated1D(energy, xs, [n_energy], [5])
idx += n_energy
# Copy binding energy
shell = _subshell(d)
e = data.atomic_relaxation.binding_energy[shell]
rx.subshell_binding_energy = e
return data
@classmethod
def from_endf(cls, photoatomic, relaxation=None):
"""Generate incident photon data from an ENDF evaluation
@ -548,15 +705,15 @@ class IncidentPhoton(EqualityMixin):
if rx.scattering_factor is not None:
rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor')
# Write pair production cross section
# Write electron-field pair production cross section
if 515 in self:
pair_group = group.create_group('pair_production')
pair_group = group.create_group('pair_production_electron')
pair_group.create_dataset('xs', data=self[515].xs(union_grid))
# Write triplet production cross section
# Write nuclear-field pair production cross section
if 517 in self:
triplet_group = group.create_group('triplet_production')
triplet_group.create_dataset('xs', data=self[517].xs(union_grid))
pair_group = group.create_group('pair_production_nuclear')
pair_group.create_dataset('xs', data=self[517].xs(union_grid))
# Write photoelectric cross section
photoelec_group = group.create_group('photoelectric')
@ -713,7 +870,92 @@ class PhotonReaction(EqualityMixin):
self._xs = xs
@classmethod
def from_endf(self, ev, mt):
def from_ace(cls, ace, mt):
"""Generate photon reaction from an ACE table
Parameters
----------
ace : openmc.data.ace.Table
ACE table to read from
mt : int
The MT value of the reaction to get data for
Returns
-------
openmc.data.PhotonReaction
Photon reaction data
"""
# Create instance
rx = cls(mt)
# Get energy grid (stored as logarithms)
n = ace.nxs[3]
idx = ace.jxs[1]
energy = np.exp(ace.xss[idx : idx+n])*EV_PER_MEV
# Get index for appropriate reaction
if mt == 502:
# Coherent scattering
idx = ace.jxs[1] + 2*n
elif mt == 504:
# Incoherent scattering
idx = ace.jxs[1] + n
elif mt == 515:
# Pair production
idx = ace.jxs[1] + 4*n
elif mt == 522:
# Photoelectric
idx = ace.jxs[1] + 3*n
else:
raise ValueError('ACE photoatomic cross sections do not have '
'data for MT={}.'.format(mt))
# Store cross section
xs = ace.xss[idx : idx+n].copy()
nonzero = (xs != 0.0)
xs[nonzero] = np.exp(xs[nonzero])
rx.xs = Tabulated1D(energy, xs, [n], [5])
# Get form factors for incoherent/coherent scattering
new_format = (ace.nxs[6] > 0)
if mt == 502:
idx = ace.jxs[3]
if new_format:
n = (ace.jxs[4] - ace.jxs[3]) // 3
x = ace.xss[idx : idx+n]
idx += n
else:
x = np.array([
0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12,
0.15, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55,
0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4,
3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6,
5.8, 6.0])
n = x.size
ff = ace.xss[idx+n : idx+2*n]
rx.scattering_factor = Tabulated1D(x, ff)
elif mt == 504:
idx = ace.jxs[2]
if new_format:
n = (ace.jxs[3] - ace.jxs[2]) // 2
x = ace.xss[idx : idx+n]
idx += n
else:
x = np.array([
0.0, 0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6,
0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 8.0
])
n = x.size
ff = ace.xss[idx : idx+n]
rx.scattering_factor = Tabulated1D(x, ff)
return rx
@classmethod
def from_endf(cls, ev, mt):
"""Generate photon reaction from an ENDF evaluation
Parameters
@ -729,7 +971,7 @@ class PhotonReaction(EqualityMixin):
Photon reaction data
"""
rx = PhotonReaction(mt)
rx = cls(mt)
# Read photon cross section
if (23, mt) in ev.section:

1
scripts/openmc-get-photon-data
View file

@ -6,7 +6,6 @@ relaxation data and convert it to an HDF5 library for use with OpenMC.
This data is used for photon transport in OpenMC.
"""
from __future__ import print_function
import os
import sys
import shutil

4
scripts/openmc-make-compton
View file

@ -1,12 +1,10 @@
#!/usr/bin/env python
from __future__ import print_function, division
import os
import sys
import tarfile
from urllib.request import urlopen
from six.moves import input
from six.moves.urllib.request import urlopen
import numpy as np
import h5py

8
scripts/openmc-make-stopping-powers
View file

@ -1,13 +1,11 @@
#!/usr/bin/env python
from __future__ import print_function
from six.moves.urllib.parse import urlencode
from six.moves.urllib.request import urlopen
from urllib.parse import urlencode
from urllib.request import urlopen
from lxml import html
import numpy as np
import h5py
from openmc.data import ATOMIC_SYMBOL

3
scripts/openmc-update-inputs
View file

@ -248,8 +248,7 @@ def update_materials(root):
# If a nuclide name is in the ZAID notation (e.g., a number),
# convert it to the proper nuclide name.
if nucname.strip().isnumeric():
nucname = \
openmc.data.neutron._get_metadata(int(nucname))[0]
nucname = openmc.data.ace.get_metadata(int(nucname))[0]
nucname = nucname.replace('Nat', '0')
if nucname.endswith('m'):
nucname = nucname[:-1] + '_m1'

6
src/api.F90
View file

@ -116,7 +116,8 @@ contains
check_overlaps = .false.
confidence_intervals = .false.
create_fission_neutrons = .true.
energy_cutoff = ZERO
electron_treatment = ELECTRON_LED
energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO]
energy_max_neutron = INFINITY
energy_min_neutron = ZERO
entropy_on = .false.
@ -135,6 +136,7 @@ contains
output_summary = .true.
output_tallies = .true.
particle_restart_run = .false.
photon_transport = .false.
pred_batches = .false.
reduce_tallies = .true.
res_scat_on = .false.
@ -305,6 +307,7 @@ contains
use cmfd_header
use mgxs_header
use photon_header
use plot_header
use sab_header
use settings
@ -321,6 +324,7 @@ contains
call free_memory_volume()
call free_memory_simulation()
call free_memory_nuclide()
call free_memory_photon()
call free_memory_settings()
call free_memory_mgxs()
call free_memory_sab()

4
src/constants.F90
View file

@ -461,8 +461,8 @@ module constants
MODE_PARTICLE = 4, & ! Particle restart mode
MODE_VOLUME = 5 ! Volume calculation mode
! Electron treatments
integer, parameter :: &
! Electron treatments
integer, parameter :: &
ELECTRON_LED = 1, & ! Local Energy Deposition
ELECTRON_TTB = 2 ! Thick Target Bremsstrahlung

7
src/geometry.F90
View file

@ -206,9 +206,12 @@ contains
end do
j = p % n_coord
! Determine universe (if not set, use root universe
! Determine universe (if not yet set, use root universe)
i_universe = p % coord(j) % universe
if (i_universe == NONE) i_universe = root_universe
if (i_universe == NONE) then
p % coord(j) % universe = root_universe
i_universe = root_universe
end if
! set size of list to search
if (present(search_cells)) then

31
src/photon_header.F90
View file

@ -179,14 +179,18 @@ contains
call close_group(rgroup)
! Read pair production
rgroup = open_group(group_id, 'pair_production')
call read_dataset(this % pair_production_nuclear, rgroup, 'xs')
rgroup = open_group(group_id, 'pair_production_electron')
call read_dataset(this % pair_production_electron, rgroup, 'xs')
call close_group(rgroup)
! Read pair production
rgroup = open_group(group_id, 'triplet_production')
call read_dataset(this % pair_production_electron, rgroup, 'xs')
call close_group(rgroup)
if (object_exists(group_id, 'pair_production_nuclear')) then
rgroup = open_group(group_id, 'pair_production_nuclear')
call read_dataset(this % pair_production_nuclear, rgroup, 'xs')
call close_group(rgroup)
else
this % pair_production_nuclear(:) = ZERO
end if
! Read photoelectric
rgroup = open_group(group_id, 'photoelectric')
@ -430,4 +434,21 @@ contains
end subroutine photon_calculate_xs
!===============================================================================
! FREE_MEMORY_PHOTON deallocates/resets global variables in this module
!===============================================================================
subroutine free_memory_photon()
! Deallocate photon cross section data
if (allocated(elements)) deallocate(elements)
if (allocated(compton_profile_pz)) deallocate(compton_profile_pz)
n_elements = 0
call element_dict % clear()
! Clear TTB-related arrays
if (allocated(ttb_e_grid)) deallocate(ttb_e_grid)
if (allocated(ttb_k_grid)) deallocate(ttb_k_grid)
if (allocated(ttb)) deallocate(ttb)
end subroutine free_memory_photon
end module photon_header

23
src/physics.F90
View file

@ -51,6 +51,13 @@ contains
call sample_positron_reaction(p)
end if
! Kill particle if energy falls below cutoff
if (p % E < energy_cutoff(p % type)) then
p % alive = .false.
p % wgt = ZERO
p % last_wgt = ZERO
end if
! Display information about collision
if (verbosity >= 10 .or. trace) then
if (p % type == NEUTRON) then
@ -128,12 +135,6 @@ contains
! exiting neutron
call scatter(p, i_nuclide, i_nuc_mat)
! Play russian roulette if survival biasing is turned on
if (survival_biasing) then
call russian_roulette(p)
if (.not. p % alive) return
end if
! Advance URR seed stream 'N' times after energy changes
if (p % E /= p % last_E) then
call prn_set_stream(STREAM_URR_PTABLE)
@ -141,6 +142,12 @@ contains
call prn_set_stream(STREAM_TRACKING)
end if
! Play russian roulette if survival biasing is turned on
if (survival_biasing) then
call russian_roulette(p)
if (.not. p % alive) return
end if
end subroutine sample_neutron_reaction
!===============================================================================
@ -170,7 +177,9 @@ contains
real(8) :: uvw(3) ! new direction
real(8) :: rel_vel ! relative velocity of electron
! Kill photon if below energy cutoff
! Kill photon if below energy cutoff -- an extra check is made here because
! photons with energy below the cutoff may have been produced by neutrons
! reactions or atomic relaxation
if (p % E < energy_cutoff(PHOTON)) then
p % E = ZERO
p % alive = .false.