mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Implement LinLinEnergyFilter in Fortran
This commit is contained in:
Sterling Harper
parent
20a628e537
commit
429667e8b8
4 changed files with 216 additions and 13 deletions
|
|
@ -1436,3 +1436,79 @@ class DelayedGroupFilter(IntegralFilter):
|
|||
filter's bins.
|
||||
|
||||
"""
|
||||
|
||||
|
||||
class LinLinEnergyFilter(Filter):
|
||||
"""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energy : Iterable of Real
|
||||
A grid of energy values in eV.
|
||||
y : iterable of Real
|
||||
A grid of interpolant values in eV.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
energy : Iterable of Real
|
||||
A grid of energy values in eV.
|
||||
y : iterable of Real
|
||||
A grid of interpolant values in eV.
|
||||
num_bins : Integral
|
||||
The number of filter bins (always 1 for this filter)
|
||||
stride : Integral
|
||||
The number of filter, nuclide and score bins within each of this
|
||||
filter's bins.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, energy, y):
|
||||
self.energy = energy
|
||||
self.y = y
|
||||
self._stride = None
|
||||
|
||||
@property
|
||||
def energy(self):
|
||||
return self._energy
|
||||
|
||||
@property
|
||||
def y(self):
|
||||
return self._y
|
||||
|
||||
@energy.setter
|
||||
def energy(self, energy):
|
||||
# Make sure the energy grid is a numpy array.
|
||||
energy = np.array(energy)
|
||||
|
||||
# If the grid is 0D numpy array, promote to 1D.
|
||||
if energy.shape == ():
|
||||
energy.shape = (1,)
|
||||
|
||||
# Make sure the values are Real and positive.
|
||||
cv.check_type('filter energy grid', energy, Iterable, Real)
|
||||
for E in energy:
|
||||
cv.check_greater_than('filter energy grid', E, 0, equality=True)
|
||||
|
||||
self._energy = energy
|
||||
|
||||
@y.setter
|
||||
def y(self, y):
|
||||
# Make sure the values are in a numpy array.
|
||||
y = np.array(y)
|
||||
|
||||
# If the array is 0D, promote to 1D.
|
||||
if y.shape == ():
|
||||
y.shape = (1,)
|
||||
|
||||
# Make sure the values are Real.
|
||||
cv.check_type('filter interpolant values', y, Iterable, Real)
|
||||
|
||||
self._y = y
|
||||
|
||||
def to_xml(self):
|
||||
"""Return XML Element representing the Filter."""
|
||||
element = ET.Element('filter')
|
||||
element.set('type', self.short_name.lower())
|
||||
element.set('energy', ' '.join(str(e) for e in self.energy))
|
||||
element.set('y', ' '.join(str(y) for y in self.y))
|
||||
return element
|
||||
|
|
|
|||
|
|
@ -342,7 +342,7 @@ module constants
|
|||
integer, parameter :: NO_BIN_FOUND = -1
|
||||
|
||||
! Tally filter and map types
|
||||
integer, parameter :: N_FILTER_TYPES = 13
|
||||
integer, parameter :: N_FILTER_TYPES = 14
|
||||
integer, parameter :: &
|
||||
FILTER_UNIVERSE = 1, &
|
||||
FILTER_MATERIAL = 2, &
|
||||
|
|
@ -356,7 +356,8 @@ module constants
|
|||
FILTER_MU = 10, &
|
||||
FILTER_POLAR = 11, &
|
||||
FILTER_AZIMUTHAL = 12, &
|
||||
FILTER_DELAYEDGROUP = 13
|
||||
FILTER_DELAYEDGROUP = 13, &
|
||||
FILTER_LINLINENERGY = 14
|
||||
|
||||
! Mesh types
|
||||
integer, parameter :: &
|
||||
|
|
|
|||
|
|
@ -2973,18 +2973,21 @@ contains
|
|||
temp_str = to_lower(temp_str)
|
||||
|
||||
! Determine number of bins
|
||||
if (check_for_node(node_filt, "bins")) then
|
||||
if (temp_str == 'energy' .or. temp_str == 'energyout' .or. &
|
||||
temp_str == 'mu' .or. temp_str == 'polar' .or. &
|
||||
temp_str == 'azimuthal') then
|
||||
n_words = get_arraysize_double(node_filt, "bins")
|
||||
else
|
||||
n_words = get_arraysize_integer(node_filt, "bins")
|
||||
select case(temp_str)
|
||||
case ("energy", "energyout", "mu", "polar", "azimuthal")
|
||||
if (.not. check_for_node(node_filt, "bins")) then
|
||||
call fatal_error("Bins not set in filter on tally " &
|
||||
// trim(to_str(t % id)))
|
||||
end if
|
||||
else
|
||||
call fatal_error("Bins not set in filter on tally " &
|
||||
// trim(to_str(t % id)))
|
||||
end if
|
||||
n_words = get_arraysize_double(node_filt, "bins")
|
||||
case ("mesh", "universe", "material", "cell", "distribcell", &
|
||||
"cellborn", "surface", "delayedgroup")
|
||||
if (.not. check_for_node(node_filt, "bins")) then
|
||||
call fatal_error("Bins not set in filter on tally " &
|
||||
// trim(to_str(t % id)))
|
||||
end if
|
||||
n_words = get_arraysize_integer(node_filt, "bins")
|
||||
end select
|
||||
|
||||
! Determine type of filter
|
||||
select case (temp_str)
|
||||
|
|
@ -3282,6 +3285,39 @@ contains
|
|||
! Set the filter index in the tally find_filter array
|
||||
t % find_filter(FILTER_AZIMUTHAL) = j
|
||||
|
||||
case ('linlinenergy')
|
||||
! Allocate and declare the filter type.
|
||||
allocate(LinLinEnergyFilter::t % filters(j) % obj)
|
||||
select type (filt => t % filters(j) % obj)
|
||||
type is (LinLinEnergyFilter)
|
||||
filt % n_bins = 1
|
||||
! Make sure this is continuous-energy mode.
|
||||
if (.not. run_CE) then
|
||||
call fatal_error("LinLinEnergy filters are only supported for &
|
||||
&continuous-energy transport calculations")
|
||||
end if
|
||||
|
||||
! Allocate and store energy grid.
|
||||
if (.not. check_for_node(node_filt, "energy")) then
|
||||
call fatal_error("Energy grid not specified for LinLinEnergy &
|
||||
&filter on tally " // trim(to_str(t % id)))
|
||||
end if
|
||||
n_words = get_arraysize_double(node_filt, "energy")
|
||||
allocate(filt % energy(n_words))
|
||||
call get_node_array(node_filt, "energy", filt % energy)
|
||||
|
||||
! Allocate and store interpolant values.
|
||||
if (.not. check_for_node(node_filt, "y")) then
|
||||
call fatal_error("y values not specified for LinLinEnergy &
|
||||
&filter on tally " // trim(to_str(t % id)))
|
||||
end if
|
||||
n_words = get_arraysize_double(node_filt, "y")
|
||||
allocate(filt % y(n_words))
|
||||
call get_node_array(node_filt, "y", filt % y)
|
||||
end select
|
||||
! Set the filter index in the tally find_filter array
|
||||
t % find_filter(FILTER_LINLINENERGY) = j
|
||||
|
||||
case default
|
||||
! Specified tally filter is invalid, raise error
|
||||
call fatal_error("Unknown filter type '" &
|
||||
|
|
|
|||
|
|
@ -191,6 +191,19 @@ module tally_filter
|
|||
procedure :: text_label => text_label_azimuthal
|
||||
end type AzimuthalFilter
|
||||
|
||||
!===============================================================================
|
||||
! LinLinEnergyFilter
|
||||
!===============================================================================
|
||||
type, extends(TallyFilter) :: LinLinEnergyFilter
|
||||
real(8), allocatable :: energy(:)
|
||||
real(8), allocatable :: y(:)
|
||||
|
||||
contains
|
||||
procedure :: get_next_bin => get_next_bin_linlinenergy
|
||||
procedure :: to_statepoint => to_statepoint_linlinenergy
|
||||
procedure :: text_label => text_label_linlinenergy
|
||||
end type LinLinEnergyFilter
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -1209,6 +1222,83 @@ contains
|
|||
// trim(to_str(E1)) // ")"
|
||||
end function text_label_azimuthal
|
||||
|
||||
!===============================================================================
|
||||
! LinLinEnergyFilter methods
|
||||
!===============================================================================
|
||||
subroutine get_next_bin_linlinenergy(this, p, estimator, current_bin, &
|
||||
next_bin, weight)
|
||||
class(LinLinEnergyFilter), intent(in) :: this
|
||||
type(Particle), intent(in) :: p
|
||||
integer, intent(in) :: estimator
|
||||
integer, value, intent(in) :: current_bin
|
||||
integer, intent(out) :: next_bin
|
||||
real(8), intent(out) :: weight
|
||||
|
||||
integer :: n, indx
|
||||
real(8) :: E, f
|
||||
|
||||
select type(this)
|
||||
type is (LinLinEnergyFilter)
|
||||
if (current_bin == NO_BIN_FOUND) then
|
||||
n = size(this % energy)
|
||||
|
||||
! Make sure the correct energy is used.
|
||||
if (estimator == ESTIMATOR_TRACKLENGTH) then
|
||||
E = p % E
|
||||
else
|
||||
E = p % last_E
|
||||
end if
|
||||
|
||||
! Check if energy of the particle is within energy bins.
|
||||
if (E < this % energy(1) .or. E > this % energy(n)) then
|
||||
next_bin = NO_BIN_FOUND
|
||||
weight = ONE
|
||||
|
||||
else
|
||||
! Search to find incoming energy bin.
|
||||
indx = binary_search(this % energy, n, E)
|
||||
|
||||
! Compute an interpolation factor between nearest bins.
|
||||
f = (E - this % energy(indx)) &
|
||||
/ (this % energy(indx+1) - this % energy(indx))
|
||||
|
||||
! Interpolate on the lin-lin grid.
|
||||
next_bin = 1
|
||||
weight = (ONE - f) * this % y(indx) + f * this % y(indx+1)
|
||||
end if
|
||||
|
||||
else
|
||||
next_bin = NO_BIN_FOUND
|
||||
weight = ONE
|
||||
end if
|
||||
end select
|
||||
end subroutine get_next_bin_linlinenergy
|
||||
|
||||
subroutine to_statepoint_linlinenergy(this, filter_group)
|
||||
class(LinLinEnergyFilter), intent(in) :: this
|
||||
integer(HID_T), intent(in) :: filter_group
|
||||
|
||||
select type(this)
|
||||
type is (LinLinEnergyFilter)
|
||||
call write_dataset(filter_group, "type", "linlinenergy")
|
||||
call write_dataset(filter_group, "energy", this % energy)
|
||||
call write_dataset(filter_group, "y", this % y)
|
||||
end select
|
||||
end subroutine to_statepoint_linlinenergy
|
||||
|
||||
function text_label_linlinenergy(this, bin) result(label)
|
||||
class(LinLinEnergyFilter), intent(in) :: this
|
||||
integer, intent(in) :: bin
|
||||
character(MAX_LINE_LEN) :: label
|
||||
|
||||
select type(this)
|
||||
type is (LinLinEnergyFilter)
|
||||
label = "Lin-Lin Energy Interpolation [" &
|
||||
// trim(to_str(this % energy(1))) // ", ..., " &
|
||||
// trim(to_str(this % energy(size(this % energy)))) // "]"
|
||||
end select
|
||||
end function text_label_linlinenergy
|
||||
|
||||
!===============================================================================
|
||||
! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and
|
||||
! a target offset to build a string which is the path from the base universe to
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue