From 14b8ef6a1f6bd248a22e6c4cf806d27fcd812fca Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 5 Aug 2016 13:48:12 -0500 Subject: [PATCH] Make PyAPI Polynomial class, remove Constant1D --- openmc/data/function.py | 91 ++++++++++++++++++++++++++++++++++++++++- openmc/data/product.py | 28 +++---------- openmc/data/reaction.py | 6 +-- src/endf_header.F90 | 30 -------------- src/nuclide_header.F90 | 8 ++-- src/product_header.F90 | 8 ++-- src/tally.F90 | 44 +++++--------------- 7 files changed, 115 insertions(+), 100 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index bea6f5e9a..e827e1283 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,3 +1,4 @@ +from abc import ABCMeta, abstractmethod from collections import Iterable, Callable from numbers import Real, Integral @@ -9,7 +10,53 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} -class Tabulated1D(object): +class Function1D(object): + """A function of one independent variable with HDF5 support.""" + + __meta__class = ABCMeta + + def __init__(self): pass + + @abstractmethod + def __call__(self): pass + + @abstractmethod + def to_hdf5(self, group, name='xy'): + """Write function to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + pass + + @classmethod + def from_hdf5(cls, dataset): + """Generate function from an HDF5 dataset + + Parameters + ---------- + dataset : h5py.Dataset + Dataset to read from + + Returns + ------- + openmc.data.Function1D + Function read from dataset + + """ + for subclass in cls.__subclasses__(): + if dataset.attrs['type'].decode() == subclass.__name__: + return subclass.from_hdf5(dataset) + raise ValueError("Unrecognized Function1D class: '" + + dataset.attrs['type'].decode() + "'") + + +class Tabulated1D(Function1D): """A one-dimensional tabulated function. This class mirrors the TAB1 type from the ENDF-6 format. A tabulated @@ -239,7 +286,7 @@ class Tabulated1D(object): """ dataset = group.create_dataset(name, data=np.vstack( [self.x, self.y])) - dataset.attrs['type'] = np.string_('tab1') + dataset.attrs['type'] = np.string_(type(self).__name__) dataset.attrs['breakpoints'] = self.breakpoints dataset.attrs['interpolation'] = self.interpolation @@ -258,6 +305,10 @@ class Tabulated1D(object): Function read from dataset """ + if dataset.attrs['type'].decode() != cls.__name__: + raise ValueError("Expected an HDF5 attribute 'type' equal to '" + + cls.__name__ + "'") + x = dataset.value[0, :] y = dataset.value[1, :] breakpoints = dataset.attrs['breakpoints'] @@ -304,6 +355,42 @@ class Tabulated1D(object): return Tabulated1D(x, y, breakpoints, interpolation) +class Polynomial(np.polynomial.Polynomial, Function1D): + def to_hdf5(self, group, name='xy'): + """Write polynomial function to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + dataset = group.create_dataset(name, data=self.coef) + dataset.attrs['type'] = np.string_(type(self).__name__) + + @classmethod + def from_hdf5(cls, dataset): + """Generate function from an HDF5 dataset + + Parameters + ---------- + dataset : h5py.Dataset + Dataset to read from + + Returns + ------- + openmc.data.Function1D + Function read from dataset + + """ + if dataset.attrs['type'].decode() != cls.__name__: + raise ValueError("Expected an HDF5 attribute 'type' equal to '" + + cls.__name__ + "'") + return cls(dataset.value) + + class Sum(object): """Sum of multiple functions. diff --git a/openmc/data/product.py b/openmc/data/product.py index dd276daac..116905f7a 100644 --- a/openmc/data/product.py +++ b/openmc/data/product.py @@ -3,10 +3,9 @@ from numbers import Real import sys import numpy as np -from numpy.polynomial.polynomial import Polynomial import openmc.checkvalue as cv -from .function import Tabulated1D +from .function import Tabulated1D, Polynomial, Function1D from .angle_energy import AngleEnergy if sys.version_info[0] >= 3: @@ -36,7 +35,7 @@ class Product(object): yield represents particles from prompt and delayed sources. particle : str What particle the reaction product is. - yield_ : float or openmc.data.Tabulated1D or numpy.polynomial.Polynomial + yield_ : float or openmc.data.Tabulated1D or openmc.data.Polynomial Yield of secondary particle in the reaction. """ @@ -47,7 +46,7 @@ class Product(object): self.emission_mode = 'prompt' self.distribution = [] self.applicability = [] - self.yield_ = 1 + self.yield_ = Polynomial((1,)) # 0-order polynomial i.e. a constant def __repr__(self): if isinstance(self.yield_, Real): @@ -120,7 +119,7 @@ class Product(object): @yield_.setter def yield_(self, yield_): cv.check_type('product yield', yield_, - (Real, Tabulated1D, Polynomial)) + (Tabulated1D, Polynomial)) self._yield = yield_ def to_hdf5(self, group): @@ -138,16 +137,7 @@ class Product(object): group.attrs['decay_rate'] = self.decay_rate # Write yield - if isinstance(self.yield_, Tabulated1D): - self.yield_.to_hdf5(group, 'yield') - dset = group['yield'] - dset.attrs['type'] = np.string_('tabulated') - elif isinstance(self.yield_, Polynomial): - dset = group.create_dataset('yield', data=self.yield_.coef) - dset.attrs['type'] = np.string_('polynomial') - else: - dset = group.create_dataset('yield', data=float(self.yield_)) - dset.attrs['type'] = np.string_('constant') + self.yield_.to_hdf5(group, 'yield') # Write applicability/distribution group.attrs['n_distribution'] = len(self.distribution) @@ -180,13 +170,7 @@ class Product(object): p.decay_rate = group.attrs['decay_rate'] # Read yield - yield_type = group['yield'].attrs['type'].decode() - if yield_type == 'constant': - p.yield_ = group['yield'].value - elif yield_type == 'polynomial': - p.yield_ = Polynomial(group['yield'].value) - elif yield_type == 'tabulated': - p.yield_ = Tabulated1D.from_hdf5(group['yield']) + p.yield_ = Function1D.from_hdf5(group['yield']) # Read applicability/distribution n_distribution = group.attrs['n_distribution'] diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index ad707d276..d35599579 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -5,13 +5,12 @@ from numbers import Real from warnings import warn import numpy as np -from numpy.polynomial import Polynomial import openmc.checkvalue as cv from openmc.stats import Uniform from .angle_distribution import AngleDistribution from .angle_energy import AngleEnergy -from .function import Tabulated1D +from .function import Tabulated1D, Polynomial from .data import REACTION_NAME from .product import Product from .uncorrelated import UncorrelatedAngleEnergy @@ -465,7 +464,8 @@ class Reaction(object): idx = ace.jxs[11] + abs(ty) - 101 yield_ = Tabulated1D.from_ace(ace, idx) else: - yield_ = abs(ty) + # 0-order polynomial i.e. a constant + yield_ = Polynomial((abs(ty),)) neutron = Product('neutron') neutron.yield_ = yield_ diff --git a/src/endf_header.F90 b/src/endf_header.F90 index c4b4ef3ad..8d8aefaa3 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -30,17 +30,6 @@ module endf_header end subroutine function1d_from_hdf5_ end interface -!=============================================================================== -! CONSTANT1D represents a constant one-dimensional function -!=============================================================================== - - type, extends(Function1D) :: Constant1D - real(8) :: y - contains - procedure :: from_hdf5 => constant1d_from_hdf5 - procedure :: evaluate => constant1d_evaluate - end type Constant1D - !=============================================================================== ! POLYNOMIAL represents a one-dimensional function expressed as a polynomial !=============================================================================== @@ -72,25 +61,6 @@ module endf_header contains -!=============================================================================== -! Constant1D implementation -!=============================================================================== - - subroutine constant1d_from_hdf5(this, dset_id) - class(Constant1D), intent(inout) :: this - integer(HID_T), intent(in) :: dset_id - - call read_dataset(this % y, dset_id) - end subroutine constant1d_from_hdf5 - - pure function constant1d_evaluate(this, x) result(y) - class(Constant1D), intent(in) :: this - real(8), intent(in) :: x - real(8) :: y - - y = this % y - end function constant1d_evaluate - !=============================================================================== ! Polynomial implementation !=============================================================================== diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8ff83482e..514958a07 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -10,7 +10,7 @@ module nuclide_header use constants use dict_header, only: DictIntInt use endf, only: reaction_name, is_fission, is_disappearance - use endf_header, only: Function1D, Constant1D, Polynomial, Tabulated1D + use endf_header, only: Function1D, Polynomial, Tabulated1D use error, only: fatal_error, warning use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, read_dataset, close_dataset, get_shape @@ -283,11 +283,9 @@ module nuclide_header total_nu = open_dataset(nu_group, 'yield') call read_attribute(temp, total_nu, 'type') select case (temp) - case ('constant') - allocate(Constant1D :: this % total_nu) - case ('tabulated') + case ('Tabulated1D') allocate(Tabulated1D :: this % total_nu) - case ('polynomial') + case ('Polynomial') allocate(Polynomial :: this % total_nu) end select call this % total_nu % from_hdf5(total_nu) diff --git a/src/product_header.F90 b/src/product_header.F90 index f20adf0d4..a69929473 100644 --- a/src/product_header.F90 +++ b/src/product_header.F90 @@ -5,7 +5,7 @@ module product_header use angleenergy_header, only: AngleEnergyContainer use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, & EMISSION_TOTAL, NEUTRON, PHOTON - use endf_header, only: Tabulated1D, Function1D, Constant1D, Polynomial + use endf_header, only: Tabulated1D, Function1D, Polynomial use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, close_dataset, read_dataset use random_lcg, only: prn @@ -109,11 +109,9 @@ contains yield = open_dataset(group_id, 'yield') call read_attribute(temp, yield, 'type') select case (temp) - case ('constant') - allocate(Constant1D :: this % yield) - case ('tabulated') + case ('Tabulated1D') allocate(Tabulated1D :: this % yield) - case ('polynomial') + case ('Polynomial') allocate(Polynomial :: this % yield) end select call this % yield % from_hdf5(yield) diff --git a/src/tally.F90 b/src/tally.F90 index ec46a6194..0a910a969 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,7 +1,6 @@ module tally use constants - use endf_header, only: Constant1D use error, only: fatal_error use geometry_header use global @@ -247,24 +246,17 @@ contains ! reaction with neutrons in the exit channel if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. & (p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then - ! Don't waste time on very common reactions we know have multiplicities - ! of one. + ! Don't waste time on very common reactions we know have + ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p%event_nuclide)%reaction_index% & + m = nuclides(p % event_nuclide) % reaction_index % & get_key(p % event_MT) ! Get yield and apply to score - associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + associate (rxn => nuclides(p % event_nuclide) % reactions(m)) + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if @@ -289,16 +281,9 @@ contains get_key(p % event_MT) ! Get yield and apply to score - associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + associate (rxn => nuclides(p % event_nuclide) % reactions(m)) + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if @@ -324,15 +309,8 @@ contains ! Get yield and apply to score associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if