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

Merge pull request #39 from nelsonag/mg_update

Browse source 
Update to MG benchmark suite
This commit is contained in:
Paul Romano 2019-11-11 10:35:26 -06:00 • committed by GitHub
commit 88ef77e56f
4 changed files with 92 additions and 54 deletions
Download patch file Download diff file Expand all files Collapse all files

37
mg-criticality/make_library.py
View file

@ -16,6 +16,27 @@ GROUP_STRUCT = {1: openmc.mgxs.EnergyGroups(group_edges=[0.0, 20.0e6]),
###############################################################################
def check_it(scatter_matrix):
import scipy.special as ss
groups = scatter_matrix.shape[0]
orders = scatter_matrix.shape[2]
Nmu = 50
mu = np.linspace(-1, 1, Nmu)
f = np.zeros(shape=(groups, groups, Nmu))
neg = False
for gin in range(groups):
for gout in range(groups):
data = np.array(scatter_matrix[gin, gout, :], copy=True)
if np.sum(data) > 0:
data /= data[0]
for l in range(orders):
f[gin, gout, :] += ((float(l) + 0.5) *
ss.eval_legendre(l, mu) * data[l])
if np.any(f[gin, gout, :] < 0):
neg = True
return neg
def plot_it(scatter_matrix):
import matplotlib.pyplot as plt
import scipy.special as ss
@ -44,6 +65,9 @@ def plot_it(scatter_matrix):
def set_it(name, groups, order, fission, nu, absorption, scatt, total, chi):
neg = check_it(scatt)
if neg:
print(name + ' is negative!')
xsd = openmc.XSdata(name, groups)
xsd.order = order
if fission is not None:
@ -305,7 +329,6 @@ def create_2g():
scatter = np.rollaxis(scatter, 0, 3)
total = [0.2208, 0.3360]
chi = [0.575, 0.425]
Pu = set_it('Pu', groups, 0, fiss, nu, absorption, scatter, total, chi)
mg_cross_sections_file.add_xsdata(Pu)
@ -462,7 +485,7 @@ def create_2g():
capture = [0.0087078, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.0045552],
[[[0.31980, 0.004555],
[0.000000, 0.42410]]])
scatter = np.rollaxis(scatter, 0, 3)
total = [0.33588, 0.54628]
@ -491,8 +514,8 @@ def create_2g():
###########################################################################
# 5.2.2 UD2O-1
nu = [2.50, 2.50]
fiss = np.array([0.0028172, 0.097])
capture = [0.0087078, 0.02518]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
@ -526,7 +549,7 @@ def create_3g():
[0.000, 0.000, 2.000]]])
scatter = np.rollaxis(scatter, 0, 3)
total = [0.240, 0.975, 3.10]
chi = [0.96, 0.05, 0.0]
chi = [0.96, 0.04, 0.0]
URR = set_it('URR', groups, 0, fiss, nu, absorption, scatter, total, chi)
mg_cross_sections_file.add_xsdata(URR)
@ -552,10 +575,10 @@ def create_6g():
[0.000, 0.000, 2.000, 0.00, 0.00, 0.00],
[0.000, 0.000, 0.000, 2.00, 0.00, 0.00],
[0.000, 0.000, 0.000, 0.275, 0.60, 0.00],
[0.000, 0.000, 2.000, 0.033, 0.171, 0.024]]])
[0.000, 0.000, 0.000, 0.033, 0.171, 0.024]]])
scatter = np.rollaxis(scatter, 0, 3)
total = [0.240, 0.975, 3.10, 3.10, 0.975, 0.240]
chi = [0.48, 0.02, 0.0, 0.0, 0.02, 0.048]
chi = [0.48, 0.02, 0.0, 0.0, 0.02, 0.48]
URR = set_it('URR', groups, 0, fiss, nu, absorption, scatter, total, chi)
mg_cross_sections_file.add_xsdata(URR)

31
mg-criticality/make_model.py
View file

@ -1,6 +1,5 @@
#!/usr/bin/env python
import numpy as np
from subprocess import CalledProcessError
import openmc
import openmc.checkvalue as cv
@ -55,7 +54,7 @@ class Case(object):
mats.append(openmc.Material(name=self.mat_names[i]))
mats[-1].set_density('macro', 1.0)
mats[-1].add_macroscopic(macros[-1])
materials_file.add_material(mats[-1])
materials_file.append(mats[-1])
materials_file.cross_sections = GROUP_FILES[self.groups]
@ -101,6 +100,8 @@ class Case(object):
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.output = {'summary': False}
return settings_file
def make_tallies(self, r=None):
@ -223,10 +224,10 @@ class Case(object):
surfs = []
for r, rad in enumerate(self.rad):
if r == len(self.rad) - 1:
surfs.append(openmc.ZCylinder(R=rad,
surfs.append(openmc.ZCylinder(r=rad,
boundary_type='vacuum'))
else:
surfs.append(openmc.ZCylinder(R=rad))
surfs.append(openmc.ZCylinder(r=rad))
# Instantiate Cells
cells = []
@ -242,9 +243,9 @@ class Case(object):
surfs = []
for r, rad in enumerate(self.rad):
if r == len(self.rad) - 1:
surfs.append(openmc.Sphere(R=rad, boundary_type='vacuum'))
surfs.append(openmc.Sphere(r=rad, boundary_type='vacuum'))
else:
surfs.append(openmc.Sphere(R=rad))
surfs.append(openmc.Sphere(r=rad))
# Instantiate Cells
cells = []
@ -301,9 +302,15 @@ class Case(object):
geometry.export_to_xml()
def execute(self, quiet=True):
returncode = openmc.run(output=(not quiet))
spfile = 'statepoint.' + str(self.batches) + '.h5'
sp = openmc.StatePoint(spfile)
self.keff = sp.k_combined[:]
success = True
try:
openmc.run(output=(not quiet))
except CalledProcessError:
success = False
return returncode
if success:
spfile = 'statepoint.' + str(self.batches) + '.h5'
sp = openmc.StatePoint(spfile, autolink=False)
self.keff = sp.k_combined
return success

28
mg-criticality/repo.py
View file

@ -1,12 +1,10 @@
#!/usr/bin/env python
from make_model import Case
# GLOBAL DATA
batches = 5000
inactive = 500
particles = 10000
tab_leg = {'enable': True, 'num_points': 33}
# GLOBAL OPTIONS
batches = 1000
inactive = 50
particles = 2000
tab_leg = {'enable': True}
def build_cases():
@ -574,22 +572,6 @@ def build_cases():
params))
names[name] = len(cases) - 1
# # CASE -11
# case = -11
# name = 'Ua-1-1-IN'
# mat_names = ['Ua1']
# groups = 1
# order = 1
# geom = 'IN'
# rad = [1.]
# ref_k = 2.25
# mesh_dim = [20, 1, 1]
# params = {'mesh_dim': mesh_dim, 'tab_leg': tab_leg, 'batches': batches,
# 'inactive': inactive, 'particles': particles}
# cases.append(Case(case, name, mat_names, groups, order, geom, rad, ref_k,
# params))
# names[name] = len(cases) - 1
# CASE 36
case = 36
name = 'Ua-1-1-CY'

50
mg-criticality/run_cases.py
View file

@ -1,20 +1,18 @@
#!/usr/bin/env python
from __future__ import print_function
from optparse import OptionParser
import matplotlib.pyplot as plt
import numpy as np
import csv
import repo
import make_library
def print_case(case):
bias = np.abs(1.0E5 * (case.keff[0] - case.ref_k))
print('\tCalculated keff = {0:1.6f}'.format(case.keff[0]))
bias = 1.0E5 * (case.keff - case.ref_k)
print('\tCalculated keff = {0:1.6f}'.format(case.keff))
print('\tReference keff = {0:1.6f}'.format(case.ref_k))
print('\tBias [pcm] = {0:1.1f}'.format(bias))
return bias
# Command line parsing
@ -33,6 +31,10 @@ if options.list_cases:
exit()
# Run specific case, if requested
case_nums = []
case_names = []
biases = []
codes = []
if options.case_name:
if options.case_name in repo.names:
case_num = repo.names[options.case_name]
@ -44,10 +46,15 @@ if options.case_name:
print('Running Case\n')
case.make_model()
code = case.execute(False)
if code != 0:
print(case.name + ' Failed Execution!')
if code:
bias = print_case(case)
else:
print_case(case)
print(case.name + ' Failed Execution!')
bias = 0.
biases.append(bias)
case_names.append(case.name)
case_nums.append(case_num)
codes.append(code)
else:
print('Invalid Case Name: ' + options.case)
@ -62,7 +69,26 @@ else:
print(case.number, case.name)
case.make_model()
code = case.execute(True)
if code != 0:
print(case.name + ' Failed Execution!')
if code:
bias = print_case(case)
else:
print_case(case)
print(case.name + ' Failed Execution!')
bias = 0.
biases.append(bias)
case_names.append(case.name)
codes.append(code)
case_nums.append(case.number)
# Write to a CSV
with open("results.csv", mode="w") as csv_file:
writer = csv.writer(csv_file, delimiter=",", quotechar='"',
quoting=csv.QUOTE_NONNUMERIC)
# write the header
writer.writerow(['Case ID', 'Case Name', 'Eigenvalue Bias [pcm]',
'Eigenvalue Bias Std. Dev. [pcm]', 'Successful Execution'])
# Now run each case
for i in range(len(case_nums)):
writer.writerow([case_nums[i], case_names[i], biases[i].n,
biases[i].s, codes[i]])