mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Make PyAPI Polynomial class, remove Constant1D
This commit is contained in:
Sterling Harper
parent
ab432ac151
commit
14b8ef6a1f
7 changed files with 115 additions and 100 deletions
|
|
@ -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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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']
|
||||
|
|
|
|||
|
|
@ -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_
|
||||
|
|
|
|||
|
|
@ -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
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue