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

Address @paulromano comments in #1023

Browse source
This commit is contained in:
amandalund 2018-07-05 23:25:42 -05:00
parent 7dcd6239eb
commit 4ebeb737ab
11 changed files with 95 additions and 70 deletions
Download patch file Download diff file Expand all files Collapse all files

5
openmc/data/endf.py
View file

@ -24,8 +24,9 @@ from openmc.stats.univariate import Uniform, Tabular, Legendre
_LIBRARY = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF',
4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL',
31: 'INDL/V', 32: 'INDL/A', 33: 'FENDL', 34: 'IRDF',
35: 'BROND', 36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'}
17: 'TENDL', 18: 'ROSFOND', 21: 'SG-21', 31: 'INDL/V',
32: 'INDL/A', 33: 'FENDL', 34: 'IRDF', 35: 'BROND',
36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'}
_SUBLIBRARY = {
0: 'Photo-nuclear data',

35
openmc/data/photon.py
View file

@ -105,6 +105,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.
@ -328,7 +329,7 @@ class AtomicRelaxation(EqualityMixin):
class IncidentPhoton(EqualityMixin):
"""Photon interaction data.
r"""Photon interaction data.
This class stores photo-atomic, photo-nuclear, atomic relaxation,
Compton profile, stopping power, and bremsstrahlung data assembled from
@ -350,9 +351,9 @@ class IncidentPhoton(EqualityMixin):
Atomic relaxation data
bremsstrahlung : dict
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
(incident electron kinetic energy values in eV), 'photon_energy'
(incident electron kinetic energy values in [eV]), 'photon_energy'
(ratio of the energy of the emitted photon to the incident electron
kinetic energy), and 'dcs' (cross sectin values in mb). The cross
kinetic energy), and 'dcs' (cross sectin values in [b]). The cross
sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`,
where :math:`E_k` is the energy of the emitted photon.
compton_profiles : dict
@ -366,10 +367,10 @@ class IncidentPhoton(EqualityMixin):
Contains the cross sections for each photon reaction. The keys are MT
values and the values are instances of :class:`PhotonReaction`.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in eV), 'I' (mean
Dictionary of stopping power data with keys 'energy' (in [eV]), 'I' (mean
excitation energy), 's_collision' (collision stopping power in
eV cm:sup:`2`/g), and 's_radiative' (radiative stopping power in
eV cm:sup:`2`/g)
[eV cm\ :sup:`2`/g]), and 's_radiative' (radiative stopping power in
[eV cm\ :sup:`2`/g])
summed_reactions : collections.OrderedDict
Contains summed cross sections. The keys are MT values and the values
are instances of :class:`PhotonReaction`.
@ -587,13 +588,13 @@ class IncidentPhoton(EqualityMixin):
_STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV
for i in range(1, 99):
group = f['{:03}'.format(i)]
_STOPPING_POWERS[i] = {'I': group.attrs['I'],
's_collision': group['s_collision'].value,
's_radiative': group['s_radiative'].value}
# Units are in MeV cm^2/g; convert to eV cm^2/g
_STOPPING_POWERS[i]['s_collision'] *= EV_PER_MEV
_STOPPING_POWERS[i]['s_radiative'] *= EV_PER_MEV
_STOPPING_POWERS[i] = {
'I': group.attrs['I'],
's_collision': group['s_collision'].value*EV_PER_MEV,
's_radiative': group['s_radiative'].value*EV_PER_MEV
}
# Add stopping power data
if Z < 99:
@ -616,8 +617,9 @@ class IncidentPhoton(EqualityMixin):
# Index in data
p = 39
# Get log of incident electron kinetic energy values, used for cubic
# spline interpolation in log energy. Units are in MeV, so convert to eV.
# Get log of incident electron kinetic energy values, used for
# cubic spline interpolation in log energy. Units are in MeV, so
# convert to eV.
logx = np.log(np.fromiter(brem[p:p+n], float, n)*EV_PER_MEV)
p += n
@ -628,9 +630,10 @@ class IncidentPhoton(EqualityMixin):
for i in range(1, 101):
dcs = np.empty([len(log_energy), k])
# Get the scaled cross section values for each electron energy and
# reduced photon energy for this Z
y = np.reshape(np.fromiter(brem[p:p+n*k], float, n*k), (n, k))
# Get the scaled cross section values for each electron energy
# and reduced photon energy for this Z. Units are in mb, so
# convert to b.
y = np.reshape(np.fromiter(brem[p:p+n*k], float, n*k), (n, k))*1.0e-3
p += k*n
for j in range(k):

5
scripts/openmc-get-nndc-data
View file

@ -54,10 +54,7 @@ for f in files:
req = urlopen(url)
# Get file size from header
if sys.version_info[0] < 3:
file_size = int(req.info().getheaders('Content-Length')[0])
else:
file_size = req.length
file_size = req.length
downloaded = 0
# Check if file already downloaded

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

@ -12,8 +12,7 @@ import shutil
import zipfile
import argparse
from io import BytesIO
import requests
from urllib.request import urlopen
import openmc.data
@ -32,6 +31,7 @@ args = parser.parse_args()
base_url = 'http://www.nndc.bnl.gov/endf/b7.1/zips/'
files = ['ENDF-B-VII.1-photoat.zip', 'ENDF-B-VII.1-atomic_relax.zip']
block_size = 16384
# ==============================================================================
# DOWNLOAD FILES FROM NNDC SITE
@ -39,13 +39,44 @@ files = ['ENDF-B-VII.1-photoat.zip', 'ENDF-B-VII.1-atomic_relax.zip']
if not os.path.exists('photon_hdf5'):
os.mkdir('photon_hdf5')
for f in files:
# Establish connection to URL
print('Downloading {}...'.format(f))
url = base_url + f
r = requests.get(url, stream=True)
zipfile.ZipFile(BytesIO(r.content)).extractall()
req = urlopen(url)
# Get file size from header
file_size = req.length
downloaded = 0
# Check if file already downloaded
if os.path.exists(f):
if os.path.getsize(f) == file_size:
print('Skipping ' + f)
continue
else:
overwrite = input('Overwrite {}? ([y]/n) '.format(f))
if overwrite.lower().startswith('n'):
continue
# Copy file to disk
print('Downloading {}... '.format(f), end='')
with open(f, 'wb') as fh:
while True:
chunk = req.read(block_size)
if not chunk: break
fh.write(chunk)
downloaded += len(chunk)
status = '{0:10} [{1:3.2f}%]'.format(
downloaded, downloaded * 100. / file_size)
print(status + chr(8)*len(status), end='')
print('')
# ==============================================================================
# EXTRACT FILES
for f in files:
print('Extracting {0}...'.format(f))
zipfile.ZipFile(f).extractall()
# ==============================================================================
# GENERATE HDF5 DATA LIBRARY

5
scripts/openmc-make-compton
View file

@ -20,10 +20,7 @@ block_size = 16384
req = urlopen(base_url + filename)
# Get file size from header
if sys.version_info[0] < 3:
file_size = int(req.info().getheaders('Content-Length')[0])
else:
file_size = req.length
file_size = req.length
downloaded = 0
# Check if file already downloaded

11
src/material_header.F90
View file

@ -735,11 +735,7 @@ contains
mat => materials(i_material)
! Determine whether we are generating electron or positron data
if (particle == POSITRON) then
positron_ = .true.
else
positron_ = .false.
end if
positron_ = (particle == POSITRON)
! Get the size of the energy grids
n_k = size(ttb_k_grid)
@ -853,7 +849,7 @@ contains
beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON)
! Compute the integrand of the PDF
f(j) = (1.0e-3_8 * x) / (beta**2 * stopping_power(j) * w)
f(j) = x / (beta**2 * stopping_power(j) * w)
end do
! Number of points to integrate
@ -913,9 +909,6 @@ contains
this % yield = -500.0_8
end where
deallocate(stopping_power_collision, stopping_power_radiative, &
stopping_power, dcs, f, z)
end subroutine bremsstrahlung_init
end module material_header

13
src/math_functions.cpp
View file

@ -691,7 +691,7 @@ void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
}
void spline_c(int n, double x[], double y[], double z[])
void spline_c(int n, const double x[], const double y[], double z[])
{
double c_new[n-1];
@ -718,7 +718,8 @@ void spline_c(int n, double x[], double y[], double z[])
}
double spline_interpolate_c(int n, double x[], double y[], double z[], double xint)
double spline_interpolate_c(int n, const double x[], const double y[],
const double z[], double xint)
{
// Find the lower bounding index in x of xint
int i = n - 1;
@ -734,14 +735,12 @@ double spline_interpolate_c(int n, double x[], double y[], double z[], double xi
double c = z[i]/2.0;
double d = (z[i+1] - z[i])/(h*6.0);
double yint = y[i] + b*r + c*r*r + d*r*r*r;
return yint;
return y[i] + b*r + c*r*r + d*r*r*r;
}
double spline_integrate_c(int n, double x[], double y[], double z[], double xa,
double xb)
double spline_integrate_c(int n, const double x[], const double y[],
const double z[], double xa, double xb)
{
// Find the lower bounding index in x of the lower limit of integration.
int ia = n - 1;

22
src/math_functions.h
View file

@ -157,15 +157,15 @@ extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
//! used in any subsequent calls to spline_interpolate or spline_integrate for
//! the same set of x and y values.
//!
//! @param n Number of points
//! @param x Values of the independent variable, which must be strictly
//! @param n Number of points
//! @param x Values of the independent variable, which must be strictly
//! increasing.
//! @param y Values of the dependent variable.
//! @param z The second derivative of the interpolating function at each value
//! of x.
//! @param y Values of the dependent variable.
//! @param[out] z The second derivative of the interpolating function at each
//! value of x.
//==============================================================================
extern "C" void spline_c(int n, double x[], double y[], double z[]);
extern "C" void spline_c(int n, const double x[], const double y[], double z[]);
//==============================================================================
//! Determine the cubic spline interpolated y-value for a given x-value.
@ -180,8 +180,8 @@ extern "C" void spline_c(int n, double x[], double y[], double z[]);
//! @result Interpolated value
//==============================================================================
extern "C" double spline_interpolate_c(int n, double x[], double y[], double z[],
double xint);
extern "C" double spline_interpolate_c(int n, const double x[], const double y[],
const double z[], double xint);
//==============================================================================
//! Evaluate the definite integral of the interpolating cubic spline between
@ -198,8 +198,8 @@ extern "C" double spline_interpolate_c(int n, double x[], double y[], double z[]
//! @result Integral
//==============================================================================
extern "C" double spline_integrate_c(int n, double x[], double y[], double z[],
double xa, double xb);
extern "C" double spline_integrate_c(int n, const double x[], const double y[],
const double z[], double xa, double xb);
} // namespace openmc
#endif // MATH_FUNCTIONS_H
#endif // MATH_FUNCTIONS_H

4
tests/regression_tests/photon_source/inputs_true.dat
View file

@ -35,7 +35,11 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="1" type="particle">
<bins>2</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>

20
tests/regression_tests/photon_source/test.py
View file

@ -19,9 +19,7 @@ class SourceTestHarness(PyAPITestHarness):
inside_sphere = openmc.Cell()
inside_sphere.region = -sphere
inside_sphere.fill = mat
root = openmc.Universe()
root.add_cell(inside_sphere)
geometry = openmc.Geometry(root)
geometry = openmc.Geometry([inside_sphere])
geometry.export_to_xml()
source = openmc.Source()
@ -40,20 +38,22 @@ class SourceTestHarness(PyAPITestHarness):
settings.source = source
settings.export_to_xml()
particle_filter = openmc.ParticleFilter('photon')
tally = openmc.Tally()
tally.filters = [particle_filter]
tally.scores = ['flux']
tallies = openmc.Tallies([tally])
tallies.export_to_xml()
def _get_results(self):
sp = openmc.StatePoint(self._sp_name)
outstr = ''
t = sp.get_tally()
outstr += 'tally {}:\n'.format(t.id)
outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0])
outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0])
with openmc.StatePoint(self._sp_name) as sp:
outstr = ''
t = sp.get_tally()
outstr += 'tally {}:\n'.format(t.id)
outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0])
outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0])
return outstr
return outstr
def test_source():

2
tools/ci/travis-before-script.sh
View file

@ -13,7 +13,7 @@ fi
# Download ENDF/B-VII.1 distribution
ENDF=$HOME/endf-b-vii.1/
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
wget https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz -O - | tar -C $HOME -xvJ
wget https://anl.box.com/shared/static/yw7xe3k9gbps0e6muyf0cg134tq0punw.xz -O - | tar -C $HOME -xvJ
fi
# Download multipole library