diff --git a/.travis.yml b/.travis.yml index e922eba6ae..e80ab4a9c1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -43,7 +43,7 @@ install: true before_script: - if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget https://anl.box.com/shared/static/1ow880up90hgwynvh91twwfteameuwre.xz -O - | tar -C $HOME -xvJ; + wget https://anl.box.com/shared/static/oy85y1i3gboho82ifdtjlu1q4ozhd7el.xz -O - | tar -C $HOME -xvJ; fi - export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml diff --git a/LICENSE b/LICENSE index fe18b53f9f..660806c886 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011-2016 Massachusetts Institute of Technology +Copyright (c) 2011-2017 Massachusetts Institute of Technology Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/source/conf.py b/docs/source/conf.py index 4e21ec20b6..d7fc23c43e 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -66,7 +66,7 @@ master_doc = 'index' # General information about the project. project = u'OpenMC' -copyright = u'2011-2016, Massachusetts Institute of Technology' +copyright = u'2011-2017, Massachusetts Institute of Technology' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index c8644e45cc..936f0c8388 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -153,6 +153,9 @@ Avoid extraneous whitespace in the following situations: Yes: if (variable == 2) then No: if ( variable==2 ) then +The structure component designator ``%`` should be surrounded by one space on +each side. + Do not leave trailing whitespace at the end of a line. ------ diff --git a/docs/source/license.rst b/docs/source/license.rst index c51902d6ff..1e5f88bdc5 100644 --- a/docs/source/license.rst +++ b/docs/source/license.rst @@ -4,7 +4,7 @@ License Agreement ================= -Copyright © 2011-2016 Massachusetts Institute of Technology +Copyright © 2011-2017 Massachusetts Institute of Technology Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 7f18ab28a2..1cf4fb97a0 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -56,7 +56,7 @@ Indicates the default path to a directory containing windowed multipole data if the user has not specified the tag in .I settings.xml\fP. .SH LICENSE -Copyright \(co 2011-2016 Massachusetts Institute of Technology. +Copyright \(co 2011-2017 Massachusetts Institute of Technology. .PP Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 53c6469d90..61eb30ba4c 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -1,6 +1,7 @@ from collections import Iterable from difflib import get_close_matches from numbers import Real +import re from warnings import warn import numpy as np @@ -9,7 +10,7 @@ import h5py import openmc.checkvalue as cv from openmc.mixin import EqualityMixin from . import HDF5_VERSION, HDF5_VERSION_MAJOR -from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV +from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE from .ace import Table, get_table from .angle_energy import AngleEnergy from .function import Tabulated1D @@ -17,51 +18,80 @@ from .correlated import CorrelatedAngleEnergy from openmc.stats import Discrete, Tabular -_THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27', - 'be': 'c_Be', - 'beo': 'c_BeO', - 'bebeo': 'c_Be_in_BeO', 'be-o': 'c_Be_in_BeO', 'be/o': 'c_Be_in_BeO', - 'benz': 'c_Benzine', - 'cah': 'c_Ca_in_CaH2', - 'dd2o': 'c_D_in_D2O', 'hwtr': 'c_D_in_D2O', 'hw': 'c_D_in_D2O', - 'fe': 'c_Fe56', 'fe56': 'c_Fe56', - 'graph': 'c_Graphite', 'grph': 'c_Graphite', 'gr': 'c_Graphite', - 'hca': 'c_H_in_CaH2', - 'hch2': 'c_H_in_CH2', 'poly': 'c_H_in_CH2', 'pol': 'c_H_in_CH2', - 'hh2o': 'c_H_in_H2O', 'lwtr': 'c_H_in_H2O', 'lw': 'c_H_in_H2O', - 'hzrh': 'c_H_in_ZrH', 'h-zr': 'c_H_in_ZrH', 'h/zr': 'c_H_in_ZrH', 'hzr': 'c_H_in_ZrH', - 'lch4': 'c_liquid_CH4', 'lmeth': 'c_liquid_CH4', - 'mg': 'c_Mg24', - 'obeo': 'c_O_in_BeO', 'o-be': 'c_O_in_BeO', 'o/be': 'c_O_in_BeO', - 'orthod': 'c_ortho_D', 'dortho': 'c_ortho_D', - 'orthoh': 'c_ortho_H', 'hortho': 'c_ortho_H', - 'ouo2': 'c_O_in_UO2', 'o2-u': 'c_O_in_UO2', 'o2/u': 'c_O_in_UO2', - 'sio2': 'c_SiO2', - 'parad': 'c_para_D', 'dpara': 'c_para_D', - 'parah': 'c_para_H', 'hpara': 'c_para_H', - 'sch4': 'c_solid_CH4', 'smeth': 'c_solid_CH4', - 'uuo2': 'c_U_in_UO2', 'u-o2': 'c_U_in_UO2', 'u/o2': 'c_U_in_UO2', - 'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH', 'zr/h': 'c_Zr_in_ZrH'} +_THERMAL_NAMES = { + 'c_Al27': ('al', 'al27'), + 'c_Be': ('be', 'be-metal'), + 'c_BeO': ('beo',), + 'c_Be_in_BeO': ('bebeo', 'be-o', 'be/o'), + 'c_C6H6': ('benz', 'c6h6'), + 'c_C_in_SiC': ('csic',), + 'c_Ca_in_CaH2': ('cah',), + 'c_D_in_D2O': ('dd2o', 'hwtr', 'hw'), + 'c_Fe56': ('fe', 'fe56'), + 'c_Graphite': ('graph', 'grph', 'gr'), + 'c_H_in_CaH2': ('hcah2',), + 'c_H_in_CH2': ('hch2', 'poly', 'pol'), + 'c_H_in_CH4_liquid': ('lch4', 'lmeth'), + 'c_H_in_CH4_solid': ('sch4', 'smeth'), + 'c_H_in_H2O': ('hh2o', 'lwtr', 'lw'), + 'c_H_in_H2O_solid': ('hice',), + 'c_H_in_C5O2H8': ('lucite', 'c5o2h8'), + 'c_H_in_YH2': ('hyh2',), + 'c_H_in_ZrH': ('hzrh', 'h-zr', 'h/zr', 'hzr'), + 'c_Mg24': ('mg', 'mg24'), + 'c_O_in_BeO': ('obeo', 'o-be', 'o/be'), + 'c_O_in_D2O': ('od2o',), + 'c_O_in_H2O_ice': ('oice',), + 'c_O_in_UO2': ('ouo2', 'o2-u', 'o2/u'), + 'c_ortho_D': ('orthod', 'dortho'), + 'c_ortho_H': ('orthoh', 'hortho'), + 'c_Si_in_SiC': ('sisic',), + 'c_SiO2_alpha': ('sio2', 'sio2a'), + 'c_SiO2_beta': ('sio2b'), + 'c_para_D': ('parad', 'dpara'), + 'c_para_H': ('parah', 'hpara'), + 'c_U_in_UO2': ('uuo2', 'u-o2', 'u/o2'), + 'c_Y_in_YH2': ('yyh2',), + 'c_Zr_in_ZrH': ('zrzrh', 'zr-h', 'zr/h') +} def get_thermal_name(name): - """Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'""" + """Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O' - if name in _THERMAL_NAMES.values(): + Parameters + ---------- + name : str + Name of an ACE thermal scattering table + + Returns + ------- + str + GND-format thermal scattering name + + """ + if name in _THERMAL_NAMES: return name - elif name.lower() in _THERMAL_NAMES: - return _THERMAL_NAMES[name.lower()] else: - # Make an educated guess?? This actually works well for - # JEFF-3.2 which stupidly uses names like lw00.32t, - # lw01.32t, etc. for different temperatures - matches = get_close_matches( - name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) - if len(matches) > 0: - return _THERMAL_NAMES[matches[0]] + for proper_name, names in _THERMAL_NAMES.items(): + if name.lower() in names: + return proper_name else: - # OK, we give up. Just use the ACE name. - return 'c_' + name + # Make an educated guess?? This actually works well for + # JEFF-3.2 which stupidly uses names like lw00.32t, + # lw01.32t, etc. for different temperatures + for proper_name, names in _THERMAL_NAMES.items(): + matches = get_close_matches( + name.lower(), names, cutoff=0.5) + if len(matches) > 0: + warn('Thermal scattering material "{}" is not recognized. ' + 'Assigning a name of {}.'.format(name, proper_name)) + return proper_name + else: + # OK, we give up. Just use the ACE name. + warn('Thermal scattering material "{0}" is not recognized. ' + 'Assigning a name of c_{0}.'.format(name)) + return 'c_' + name class CoherentElastic(EqualityMixin): @@ -425,22 +455,7 @@ class ThermalScattering(EqualityMixin): # Get new name that is GND-consistent ace_name, xs = ace.name.split('.') - if name is None: - if ace_name.lower() in _THERMAL_NAMES: - name = _THERMAL_NAMES[ace_name.lower()] - else: - # Make an educated guess?? This actually works well for JEFF-3.2 - # which stupidly uses names like lw00.32t, lw01.32t, etc. for - # different temperatures - matches = get_close_matches( - ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) - if len(matches) > 0: - name = _THERMAL_NAMES[matches[0]] - else: - # OK, we give up. Just use the ACE name. - name = 'c_' + ace.name - warn('Thermal scattering material "{}" is not recognized. ' - 'Assigning a name of {}.'.format(ace.name, name)) + name = get_thermal_name(ace_name) # Assign temperature to the running list kTs = [ace.temperature*EV_PER_MEV] @@ -523,30 +538,47 @@ class ThermalScattering(EqualityMixin): # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] + n_mu = ace.nxs[6] + 1 if idx != 0: n_energy = int(ace.xss[idx]) energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] if ace.nxs[5] == 4: + # Coherent elastic table.elastic_xs[temperatures[0]] = CoherentElastic( energy, P*EV_PER_MEV) + + # Coherent elastic shouldn't have angular distributions listed + assert n_mu == 0 else: + # Incoherent elastic table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P) - # Angular distribution - n_mu = ace.nxs[6] - if n_mu != -1: + # Angular distribution + assert n_mu > 0 idx = ace.jxs[6] table.elastic_mu_out[temperatures[0]] = \ ace.xss[idx:idx + n_energy * n_mu] table.elastic_mu_out[temperatures[0]].shape = \ (n_energy, n_mu) - # Get relevant nuclides + # Get relevant nuclides -- NJOY only allows one to specify three + # nuclides that the S(a,b) table applies to. Thus, for all elements + # other than H and Fe, we automatically add all the naturally-occurring + # isotopes. for zaid, awr in ace.pairs: if zaid > 0: Z, A = divmod(zaid, 1000) - table.nuclides.append(ATOMIC_SYMBOL[Z] + str(A)) + element = ATOMIC_SYMBOL[Z] + if element in ['H', 'Fe']: + table.nuclides.append(element + str(A)) + else: + if element + '0' not in table.nuclides: + table.nuclides.append(element + '0') + for isotope in sorted(NATURAL_ABUNDANCE): + if re.match(r'{}\d+'.format(element), isotope): + if isotope not in table.nuclides: + table.nuclides.append(isotope) return table diff --git a/scripts/openmc-get-nndc-data b/scripts/openmc-get-nndc-data index fe9d3ff5bd..92beaf9888 100755 --- a/scripts/openmc-get-nndc-data +++ b/scripts/openmc-get-nndc-data @@ -114,15 +114,19 @@ for filename in glob.glob('nndc/293.6K/ENDF-B-VII.1-neutron-293.6K/*'): shutil.move(filename, 'nndc/293.6K/') # ============================================================================== -# EDIT GRAPHITE ZAID (6012 to 6000) +# FIX ZAID ASSIGNMENTS FOR VARIOUS S(A,B) TABLES -print('Changing graphite ZAID from 6012 to 6000') -graphite = os.path.join('nndc', 'tsl', 'graphite.acer') -with open(graphite) as fh: - text = fh.read() -text = text.replace('6012', '6000', 1) -with open(graphite, 'w') as fh: - fh.write(text) +def fix_zaid(table, old, new): + filename = os.path.join('nndc', 'tsl', table) + with open(filename, 'r') as fh: + text = fh.read() + text = text.replace(old, new, 1) + with open(filename, 'w') as fh: + fh.write(text) + +print('Fixing ZAIDs for S(a,b) tables') +fix_zaid('bebeo.acer', '8016', ' 0') +fix_zaid('obeo.acer', '4009', ' 0') # ============================================================================== # PROMPT USER TO DELETE .TAR.GZ FILES diff --git a/scripts/openmc-update-inputs b/scripts/openmc-update-inputs index 613e4a2241..586e7b6dc7 100755 --- a/scripts/openmc-update-inputs +++ b/scripts/openmc-update-inputs @@ -29,7 +29,7 @@ from shutil import move import xml.etree.ElementTree as ET import openmc.data -from openmc.data.thermal import _THERMAL_NAMES + description = "Update OpenMC's input XML files to the latest format." epilog = """\ diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 39552a7169..b471ae54f6 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -80,7 +80,7 @@ contains j = j + 1 ! Don't check for S(a,b) tables if there are no more left - if (j > mat % n_sab) check_sab = .false. + if (j > size(mat % i_sab_tables)) check_sab = .false. end if end if diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4b80923544..8f25d6301f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1,5 +1,6 @@ module input_xml + use algorithm, only: find use cmfd_input, only: configure_cmfd use constants use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI @@ -2502,14 +2503,13 @@ contains ! Set number of S(a,b) tables mat % n_sab = n_sab - ! Allocate names and indices for nuclides and tables + ! Allocate names and indices for nuclides and tables -- for now we + ! allocate these as the number of S(a,b) tables listed. Since a single + ! table might apply to multiple nuclides, they are resized later if a + ! table is indeed applied to multiple nuclides. allocate(mat % sab_names(n_sab)) - allocate(mat % i_sab_nuclides(n_sab)) allocate(mat % i_sab_tables(n_sab)) - ! Initialize i_sab_nuclides - mat % i_sab_nuclides = NONE - do j = 1, n_sab ! Get pointer to S(a,b) table call get_list_item(node_sab_list, j, node_sab) @@ -4990,6 +4990,8 @@ contains integer :: m ! position for sorting integer :: temp_nuclide ! temporary value for sorting integer :: temp_table ! temporary value for sorting + type(VectorInt) :: i_sab_tables + type(VectorInt) :: i_sab_nuclides do i = 1, size(materials) ! Skip materials with no S(a,b) tables @@ -5004,20 +5006,32 @@ contains associate (sab => sab_tables(mat % i_sab_tables(k))) FIND_NUCLIDE: do j = 1, size(mat % nuclide) if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then - mat % i_sab_nuclides(k) = j - exit FIND_NUCLIDE + call i_sab_tables % push_back(k) + call i_sab_nuclides % push_back(j) end if end do FIND_NUCLIDE end associate ! Check to make sure S(a,b) table matched a nuclide - if (mat % i_sab_nuclides(k) == NONE) then + if (find(i_sab_tables, k) == -1) then call fatal_error("S(a,b) table " // trim(mat % & sab_names(k)) // " did not match any nuclide on material " & // trim(to_str(mat % id))) end if end do ASSIGN_SAB + ! Update i_sab_tables and i_sab_nuclides + deallocate(mat % i_sab_tables) + m = i_sab_tables % size() + allocate(mat % i_sab_tables(m)) + allocate(mat % i_sab_nuclides(m)) + mat % i_sab_tables(:) = i_sab_tables % data(1:m) + mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m) + + ! Clear entries in vectors for next material + call i_sab_tables % clear() + call i_sab_nuclides % clear() + ! If there are multiple S(a,b) tables, we need to make sure that the ! entries in i_sab_nuclides are sorted or else they won't be applied ! correctly in the cross_section module. The algorithm here is a simple diff --git a/src/output.F90 b/src/output.F90 index f5766187ee..a0b4a99b54 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -67,7 +67,7 @@ contains write(UNIT=OUTPUT_UNIT, FMT=*) & ' | The OpenMC Monte Carlo Code' write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Copyright | 2011-2016 Massachusetts Institute of Technology' + ' Copyright | 2011-2017 Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT=*) & ' License | http://openmc.readthedocs.io/en/latest/license.html' write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I1,".",I1)') & diff --git a/src/physics.F90 b/src/physics.F90 index 813aadf414..d644b5a364 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -546,6 +546,7 @@ contains real(8) :: c_j, c_j1 ! cumulative probability real(8) :: frac ! interpolation factor on outgoing energy real(8) :: r1 ! RNG for outgoing energy + real(8) :: mu_left, mu_right ! adjacent mu values i_temp = micro_xs(i_nuclide) % index_temp_sab @@ -737,23 +738,36 @@ contains E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) end if - ! Find angular distribution for closest outgoing energy bin - if (r1 - c_j < c_j1 - r1) then - j = j - else - j = j + 1 - end if - ! Sample outgoing cosine bin k = 1 + int(prn() * sab % n_inelastic_mu) - ! Will use mu from the randomly chosen incoming and closest outgoing - ! energy bins - mu = sab % inelastic_data(l) % mu(k, j) + ! Rather than use the sampled discrete mu directly, it is smeared over + ! a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the + ! discrete mu value itself. + associate (mu_l => sab % inelastic_data(l) % mu) + f = (r1 - c_j)/(c_j1 - c_j) + + ! Determine (k-1)th mu value + if (k == 1) then + mu_left = -ONE + else + mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j)) + end if + + ! Determine kth mu value + mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j)) + + ! Determine (k+1)th mu value + if (k == sab % n_inelastic_mu) then + mu_right = ONE - mu + else + mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu + end if + end associate + + ! Smear angle + mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF) - else - call fatal_error("Invalid secondary energy mode on S(a,b) table " & - // trim(sab_tables(i_sab) % name)) end if ! (inelastic secondary energy treatment) end if ! (elastic or inelastic) end associate diff --git a/tests/test_salphabeta/results_true.dat b/tests/test_salphabeta/results_true.dat index fb691f1686..7b889c859f 100644 --- a/tests/test_salphabeta/results_true.dat +++ b/tests/test_salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.331430E-01 3.074913E-03 +8.465449E-01 5.743852E-03