From 4dade286301afe495b560ef75fc4ff13a1f14fde Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 11 Mar 2018 14:56:37 -0500 Subject: [PATCH 1/9] Start adding support for reading photon data from ACE files --- openmc/data/__init__.py | 4 +- openmc/data/ace.py | 75 +++++++++++++++- openmc/data/neutron.py | 71 +-------------- openmc/data/photon.py | 162 +++++++++++++++++++++++++++++++++-- scripts/openmc-update-inputs | 3 +- src/photon_header.F90 | 14 +-- 6 files changed, 241 insertions(+), 88 deletions(-) diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 6dd6a9218..aaf98c574 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -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 * diff --git a/openmc/data/ace.py b/openmc/data/ace.py index 385408bd4..cbbc4e06d 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -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)]) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 99847be44..3a76cfe2a 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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] diff --git a/openmc/data/photon.py b/openmc/data/photon.py index be04ec911..4fe529900 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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 @@ -96,6 +97,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. @@ -272,7 +274,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 +371,68 @@ 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 + + return data + @classmethod def from_endf(cls, photoatomic, relaxation=None): """Generate incident photon data from an ENDF evaluation @@ -548,15 +612,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 +777,89 @@ 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 = np.exp(ace.xss[idx : idx+n]) + rx.xs = Tabulated1D(energy, xs, [n], [5]) + + # Get form factors for incoherent/coherent scattering + if mt == 502: + idx = ace.jxs[3] + if ace.nxs[6] > 0: + n = (ace.jxs[4] - ace.jxs[3]) // 2 + 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 ace.nxs[6] > 0: + 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 +875,7 @@ class PhotonReaction(EqualityMixin): Photon reaction data """ - rx = PhotonReaction(mt) + rx = cls(mt) # Read photon cross section if (23, mt) in ev.section: diff --git a/scripts/openmc-update-inputs b/scripts/openmc-update-inputs index ef7cdf47b..0d25b5b97 100755 --- a/scripts/openmc-update-inputs +++ b/scripts/openmc-update-inputs @@ -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' diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 27fe62f18..e50a871d2 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -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') From 0bbdd412c9c741e48c8c5a516e917c48a140ffb6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Mar 2018 16:31:39 -0500 Subject: [PATCH 2/9] Support atomic relaxation data and subshell xs from ACE --- openmc/data/photon.py | 130 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 113 insertions(+), 17 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 4fe529900..8644f8e7c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -24,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', @@ -196,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 @@ -227,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 = {} @@ -243,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] @@ -252,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)) @@ -263,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 @@ -431,6 +490,40 @@ class IncidentPhoton(EqualityMixin): 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 @@ -819,14 +912,17 @@ class PhotonReaction(EqualityMixin): 'data for MT={}.'.format(mt)) # Store cross section - xs = np.exp(ace.xss[idx : idx+n]) + 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 ace.nxs[6] > 0: - n = (ace.jxs[4] - ace.jxs[3]) // 2 + if new_format: + n = (ace.jxs[4] - ace.jxs[3]) // 3 x = ace.xss[idx : idx+n] idx += n else: @@ -843,7 +939,7 @@ class PhotonReaction(EqualityMixin): elif mt == 504: idx = ace.jxs[2] - if ace.nxs[6] > 0: + if new_format: n = (ace.jxs[3] - ace.jxs[2]) // 2 x = ace.xss[idx : idx+n] idx += n From f86ff8bf1ab474eca4990af08c3ac751bad69da7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 07:32:11 -0500 Subject: [PATCH 3/9] Set root universe in find_cell if not set --- src/geometry.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 10474848f..6096563d8 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 From b670fe4421cffa10e248459a0f195a581930f693 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 11:53:33 -0500 Subject: [PATCH 4/9] Make sure URR RNG stream advances after energy change --- src/physics.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 0f27a5361..1bc7e4729 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -128,12 +128,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 +135,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 !=============================================================================== From 0358901e46ada5e99aa296ed7c6659a0ab90f3eb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 11:58:59 -0500 Subject: [PATCH 5/9] Add back post-collision energy cutoff check --- src/physics.F90 | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index 1bc7e4729..982584844 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 @@ -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. From a2c4b60f6fed6e08e40b28293e7e669e139e4c87 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 12:46:33 -0500 Subject: [PATCH 6/9] Fix spacing in constants.F90 --- src/constants.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 93495ab9c..d2de91dae 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 From 093a22384fdf233728948db413619829ad6b6f98 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 13:35:38 -0500 Subject: [PATCH 7/9] Make sure photon data gets cleared at end of simulation --- src/api.F90 | 6 +++++- src/photon_header.F90 | 17 +++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/api.F90 b/src/api.F90 index d79a32d94..8f33e85e3 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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() diff --git a/src/photon_header.F90 b/src/photon_header.F90 index e50a871d2..15203e4ff 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -434,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 From e8fc084c2f23fa41362bc37a19f7b06836a3b3db Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 16:31:13 -0500 Subject: [PATCH 8/9] Document photon data format --- docs/source/io_formats/nuclear_data.rst | 100 ++++++++++++++++++++++-- 1 file changed, 94 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 2e553a4ed..68019a875 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -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 **//** @@ -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 **//kTs/** @@ -31,7 +33,7 @@ temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. :Datasets: - - **K** (*double*) -- kT values (in eV) for each temperature + - **K** (*double*) -- kT values in [eV] for each temperature TTT (in Kelvin) **//reactions/reaction_/** @@ -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 + +**//** + +:Attributes: - **Z** (*int*) -- Atomic number + +:Datasets: + - **energy** (*double[]*) -- Energies in [eV] at which cross sections + are tabulated + +**//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] + +**//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 + +**//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)` + +**//incoherent/** + +:Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b] + - **scattering_factor** (:ref:`tabulated <1d_tabulated>`) -- + +**//pair_production_electron/** + +:Datasets: - **xs** (*double[]*) -- Pair production (electron field) cross section in [b] + +**//pair_production_nuclear/** + +:Datasets: - **xs** (*double[]*) -- Pair production (nuclear field) cross section in [b] + +**//photoelectric/** + +:Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b] + +**//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] + +**//subshells/** + +:Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2' + +**//subshells//** + +: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 ------------------------------- From 4b002ef714731b445b0b16b0e136ac0126144a09 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2018 20:31:50 -0500 Subject: [PATCH 9/9] Update inputs in scripts --- scripts/openmc-get-photon-data | 1 - scripts/openmc-make-compton | 4 +--- scripts/openmc-make-stopping-powers | 8 +++----- 3 files changed, 4 insertions(+), 9 deletions(-) diff --git a/scripts/openmc-get-photon-data b/scripts/openmc-get-photon-data index b7025db74..182fdf589 100755 --- a/scripts/openmc-get-photon-data +++ b/scripts/openmc-get-photon-data @@ -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 diff --git a/scripts/openmc-make-compton b/scripts/openmc-make-compton index 1f920f366..2591dce94 100755 --- a/scripts/openmc-make-compton +++ b/scripts/openmc-make-compton @@ -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 diff --git a/scripts/openmc-make-stopping-powers b/scripts/openmc-make-stopping-powers index 5a9c255a1..8ee9602ef 100755 --- a/scripts/openmc-make-stopping-powers +++ b/scripts/openmc-make-stopping-powers @@ -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