crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Make tally derivatives seperate objects

Browse source 
Fix the OpenMP error by making these seperate objects threadprivate
This commit is contained in:
Sterling Harper 2016-01-29 16:54:07 -05:00
parent 5634993bf9
commit 55ab085e74
15 changed files with 614 additions and 383 deletions
Download patch file Download diff file Expand all files Collapse all files

1
openmc/__init__.py
View file

@ -9,6 +9,7 @@ from openmc.universe import *
from openmc.mesh import *
from openmc.filter import *
from openmc.trigger import *
from openmc.tally_derivative import *
from openmc.tallies import *
from openmc.cmfd import *
from openmc.executor import *

62
openmc/statepoint.py
View file

@ -83,6 +83,9 @@ class StatePoint(object):
Dictionary whose keys are tally IDs and whose values are Tally objects
tallies_present : bool
Indicate whether user-defined tallies are present
tally_derivatives : dict
Dictionary whose keys are tally derivative IDs and whose values are
TallyDerivative objects
version: tuple of Integral
Version of OpenMC
summary : None or openmc.summary.Summary
@ -115,6 +118,7 @@ class StatePoint(object):
self._summary = False
self._global_tallies = None
self._sparse = False
self._derivs_read = False
def close(self):
self._f.close()
@ -360,27 +364,10 @@ class StatePoint(object):
tally.num_realizations = n_realizations
# Read derivative information.
derivatives_present = self._f[
'{0}{1}/derivative present'.format(base, tally_key)].value
assert derivatives_present in (0, 1)
if derivatives_present == 1:
tally.diff_variable = self._f[
'{0}{1}/derivative/dependent variable'.format(
base, tally_key)].value.decode()
if tally.diff_variable == 'density':
tally.diff_material = self._f[
'{0}{1}/derivative/material'.format(
base, tally_key)].value
elif tally.diff_variable == 'nuclide_density':
tally.diff_material = self._f[
'{0}{1}/derivative/material'.format(
base, tally_key)].value
tally.diff_nuclide = self._f[
'{0}{1}/derivative/nuclide'.format(
base, tally_key)].value.decode()
else:
raise RuntimeError('Unrecognized tally differential'
'variable')
if 'derivative' in self._f['{0}{1}'.format(base, tally_key)]:
deriv_id = self._f['{0}{1}/derivative'.format(
base, tally_key)].value
tally.derivative = self.tally_derivatives[deriv_id]
# Read the number of Filters
n_filters = \
@ -462,6 +449,39 @@ class StatePoint(object):
def tallies_present(self):
return self._f['tallies/tallies_present'].value
@property
def tally_derivatives(self):
if not self._derivs_read:
# Initialize dictionaries for the Meshes
# Keys - Derivative IDs
# Values - TallyDerivative objects
self._derivs = {}
# Populate the dictionary if any derivatives are present.
if 'derivatives' in self._f['tallies']:
# Read the derivative ids.
base = 'tallies/derivatives'
deriv_ids = [int(k.split(' ')[1]) for k in self._f[base].keys()]
# Create each derivative object and add it to the dictionary.
for d_id in deriv_ids:
base = 'tallies/derivatives/derivative {:d}'.format(d_id)
deriv = openmc.TallyDerivative(derivative_id=d_id)
deriv.variable = self._f[base + '/dependent variable'].value
if deriv.variable == 'density':
deriv.material = self._f[base + '/material'].value
elif deriv.variable == 'nuclide_density':
deriv.material = self._f[base + '/material'].value
deriv.nuclide = self._f[base + '/nuclide'].value.decode()
else:
raise RuntimeError('Unrecognized tally differential '
'variable')
self._derivs[d_id] = deriv
self._derivs_read = True
return self._derivs
@property
def version(self):
return (self._f['version_major'].value,

117
openmc/tallies.py
View file

@ -11,7 +11,7 @@ import sys
import numpy as np
from openmc import Mesh, Filter, Trigger, Nuclide
from openmc import Mesh, Filter, Trigger, Nuclide, TallyDerivative
from openmc.cross import CrossScore, CrossNuclide, CrossFilter
from openmc.aggregate import AggregateScore, AggregateNuclide, AggregateFilter
from openmc.filter import _FILTER_TYPES
@ -95,6 +95,8 @@ class Tally(object):
sparse : bool
Whether or not the tally uses SciPy's LIL sparse matrix format for
compressed data storage
derivative : openmc.tally_derivative.TallyDerivative
A material perturbation derivative to apply to all scores in the tally.
"""
@ -107,9 +109,7 @@ class Tally(object):
self._scores = []
self._estimator = None
self._triggers = []
self._diff_variable = None
self._diff_material = None
self._diff_nuclide = None
self._derivative = None
self._num_realizations = 0
self._with_summary = False
@ -134,9 +134,7 @@ class Tally(object):
clone.id = self.id
clone.name = self.name
clone.estimator = self.estimator
clone._diff_variable = self.diff_variable
clone._diff_material = self.diff_material
clone._diff_nuclide = self.diff_nuclide
clone._derivative = self.derivative
clone.num_realizations = self.num_realizations
clone._sum = copy.deepcopy(self._sum, memo)
clone._sum_sq = copy.deepcopy(self._sum_sq, memo)
@ -194,11 +192,7 @@ class Tally(object):
return False
# Check derivatives
if self.diff_variable != other.diff_variable:
return False
if self.diff_material != other.diff_material:
return False
if self.diff_nuclide != other.diff_nuclide:
if self.derivative != other.derivative:
return False
# Check all scores
@ -225,20 +219,9 @@ class Tally(object):
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name)
if self.diff_variable is not None:
string += '{0: <16}{1}{2}\n'.format('\tDerivative', '=\t',
self.diff_variable)
if self.diff_variable == 'density':
string += '{0: <16}{1}{2}\n'.format('\tDiff_material', '=\t',
self.diff_material)
elif self.diff_variable == 'nuclide_density':
string += '{0: <16}{1}{2}\n'.format('\tDiff_material', '=\t',
self.diff_material)
string += '{0: <16}{1}{2}\n'.format('\tDiff_nuclide', '=\t',
self.diff_nuclide)
else:
raise RuntimeError("Encountered unrecognized differential "
"variable in a tally.")
if self.derivative is not None:
string += '{0: <16}id =\t{1}\n'.format('\tDerivative', '=\t',
str(self.derivative.id))
string += '{0: <16}{1}\n'.format('\tFilters', '=\t')
@ -443,16 +426,8 @@ class Tally(object):
return self._derived
@property
def diff_variable(self):
return self._diff_variable
@property
def diff_material(self):
return self._diff_material
@property
def diff_nuclide(self):
return self._diff_nuclide
def derivative(self):
return self._derivative
@property
def sparse(self):
@ -501,26 +476,11 @@ class Tally(object):
else:
self._name = ''
@diff_variable.setter
def diff_variable(self, var):
if var is not None:
cv.check_type('differential variable', var, basestring)
if var not in ('density', 'nuclide_density'):
raise ValueError("A tally differential variable must be either "
"'density' or 'nuclide_density'")
self._diff_variable = var
@diff_material.setter
def diff_material(self, mat):
if mat is not None:
cv.check_type('differential material', mat, Integral)
self._diff_material = mat
@diff_nuclide.setter
def diff_nuclide(self, nuc):
if nuc is not None:
cv.check_type('differential nuclide', nuc, basestring)
self._diff_nuclide = nuc
@derivative.setter
def derivative(self, deriv):
if deriv is not None:
cv.check_type('tally derivative', deriv, TallyDerivative)
self._derivative = deriv
def add_filter(self, new_filter):
"""Add a filter to the tally
@ -884,21 +844,6 @@ class Tally(object):
subelement = ET.SubElement(element, "nuclides")
subelement.text = nuclides.rstrip(' ')
# Optional derivative
if self.diff_variable is not None:
if self.diff_variable == 'density':
subelement = ET.SubElement(element, "derivative")
subelement.set("variable", self.diff_variable)
subelement.set("material", str(self.diff_material))
elif self.diff_variable == 'nuclide_density':
subelement = ET.SubElement(element, "derivative")
subelement.set("variable", self.diff_variable)
subelement.set("material", str(self.diff_material))
subelement.set("nuclide", self.diff_nuclide)
else:
raise RuntimeError("Encountered unrecognized differential "
"variable in a tally.")
# Scores
if len(self.scores) == 0:
msg = 'Unable to get XML for Tally ID="{0}" since it does not ' \
@ -922,6 +867,11 @@ class Tally(object):
for trigger in self.triggers:
trigger.get_trigger_xml(element)
# Optional derivatives
if self.derivative is not None:
subelement = ET.SubElement(element, "derivative")
subelement.text = str(self.derivative.id)
return element
def find_filter(self, filter_type):
@ -1421,12 +1371,12 @@ class Tally(object):
df['score'] = np.tile(self.scores, int(tile_factor))
# Include columns for derivatives if user requested it
if derivative and (self.diff_variable is not None):
df['d_variable'] = self.diff_variable
if self.diff_material is not None:
df['d_material'] = self.diff_material
if self.diff_nuclide is not None:
df['d_nuclide'] = self.diff_nuclide
if derivative and (self.derivative is not None):
df['d_variable'] = self.derivative.variable
if self.derivative.material is not None:
df['d_material'] = self.derivative.material
if self.derivative.nuclide is not None:
df['d_nuclide'] = self.derivative.nuclide
# Append columns with mean, std. dev. for each tally bin
df['mean'] = self.mean.ravel()
@ -3199,6 +3149,18 @@ class TalliesFile(object):
xml_element = mesh.get_mesh_xml()
self._tallies_file.append(xml_element)
def _create_derivative_subelements(self):
# Get a list of all derivatives referenced in a tally.
derivs = []
for tally in self._tallies:
deriv = tally.derivative
if deriv is not None and deriv not in derivs:
derivs.append(deriv)
# Add the derivatives to the XML tree.
for d in derivs:
self._tallies_file.append(d.get_derivative_xml())
def export_to_xml(self):
"""Create a tallies.xml file that can be used for a simulation.
@ -3209,6 +3171,7 @@ class TalliesFile(object):
self._create_mesh_subelements()
self._create_tally_subelements()
self._create_derivative_subelements()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._tallies_file)

182
openmc/tally_derivative.py Normal file
View file

@ -0,0 +1,182 @@
from __future__ import division
#from collections import Iterable, defaultdict
#import copy
#import os
#import pickle
#import itertools
from numbers import Integral
from xml.etree import ElementTree as ET
#import sys
#
#import numpy as np
#
#from openmc import Mesh, Filter, Trigger, Nuclide
#from openmc.cross import CrossScore, CrossNuclide, CrossFilter
#from openmc.aggregate import AggregateScore, AggregateNuclide, AggregateFilter
#from openmc.filter import _FILTER_TYPES
import openmc.checkvalue as cv
#from openmc.clean_xml import *
#
#
#if sys.version_info[0] >= 3:
# basestring = str
# "Static" variable for auto-generated TallyDerivative IDs
AUTO_TALLY_DERIV_ID = 10000
def reset_auto_tally_deriv_id():
global AUTO_TALLY_ID
AUTO_TALLY_DERIV_ID = 10000
class TallyDerivative(object):
"""A material perturbation derivative to apply to a tally.
Parameters
----------
derivative_id : Integral, optional
Unique identifier for the tally derivative. If none is specified, an
identifier will automatically be assigned
Attributes
----------
id : Integral
Unique identifier for the tally derivative
variable : str
Either 'density' or 'nuclide_density'
material : Integral
The perutrubed material
nuclide : str
The perturbed nuclide. Only needed for 'nuclide_density' derivatives.
Ex: 'Xe-135'
"""
def __init__(self, derivative_id=None):
# Initialize Tally class attributes
self.id = derivative_id
self._variable = None
self._material = None
self._nuclide = None
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a
# copy
if existing is None:
clone = type(self).__new__(type(self))
clone.id = self.id
clone.variable = self.variable
clone.material = self.material
clone.nuclide = self.nuclide
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __eq__(self, other):
if not isinstance(other, TallyDerivative):
return False
return (self.variable == other.variable and
self.material == other.material and
self.nuclide == other.nuclide)
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Tally Derivative\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
string += '{0: <16}{1}{2}\n'.format('\tVariable', '=\t',
self.variable)
if self.variable == 'density':
string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t',
self.material)
elif self.variable == 'nuclide_density':
string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t',
self.material)
string += '{0: <16}{1}{2}\n'.format('\tNuclide', '=\t',
self.nuclide)
else:
raise RuntimeError("Encountered unrecognized differential "
"variable in a tally.")
return string
@property
def id(self):
return self._id
@property
def variable(self):
return self._variable
@property
def material(self):
return self._material
@property
def nuclide(self):
return self._nuclide
@id.setter
def id(self, deriv_id):
if deriv_id is None:
global AUTO_TALLY_DERIV_ID
self._id = AUTO_TALLY_DERIV_ID
AUTO_TALLY_DERIV_ID += 1
else:
cv.check_type('tally derivative ID', deriv_id, Integral)
cv.check_greater_than('tally derivative ID', deriv_id, 0,
equality=True)
self._id = deriv_id
@variable.setter
def variable(self, var):
if var is not None:
cv.check_type('derivative variable', var, basestring)
if var not in ('density', 'nuclide_density'):
raise ValueError("A tally differential variable must be either "
"'density' or 'nuclide_density'")
self._variable = var
@material.setter
def material(self, mat):
if mat is not None:
cv.check_type('derivative material', mat, Integral)
self._material = mat
@nuclide.setter
def nuclide(self, nuc):
if nuc is not None:
cv.check_type('derivative nuclide', nuc, basestring)
self._nuclide = nuc
def get_derivative_xml(self):
"""Return XML representation of the tally derivative
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing derivative data
"""
element = ET.Element("derivative")
element.set("id", str(self.id))
element.set("variable", self.variable)
element.set("material", str(self.material))
if self.variable == 'nuclide_density':
element.set("nuclide", self.nuclide)
return element

7
src/global.F90
View file

@ -13,7 +13,8 @@ module global
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: SourceDistribution
use tally_header, only: TallyObject, TallyMap, TallyResult
use tally_header, only: TallyObject, TallyMap, TallyResult, &
TallyDerivative
use trigger_header, only: KTrigger
use timer_header, only: Timer
@ -145,6 +146,10 @@ module global
integer :: n_tallies = 0 ! # of tallies
integer :: n_user_tallies = 0 ! # of user tallies
! Tally derivatives
type(TallyDerivative), allocatable :: tally_derivs(:)
!$omp threadprivate(tally_derivs)
! Normalization for statistics
integer :: n_realizations = 0 ! # of independent realizations
real(8) :: total_weight ! total starting particle weight in realization

150
src/input_xml.F90
View file

@ -2334,12 +2334,13 @@ contains
type(Node), pointer :: node_mesh => null()
type(Node), pointer :: node_tal => null()
type(Node), pointer :: node_filt => null()
type(Node), pointer :: node_trigger=>null()
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_deriv => null()
type(NodeList), pointer :: node_mesh_list => null()
type(NodeList), pointer :: node_tal_list => null()
type(NodeList), pointer :: node_filt_list => null()
type(NodeList), pointer :: node_trigger_list => null()
type(NodeList), pointer :: node_deriv_list => null()
type(ElemKeyValueCI), pointer :: scores
type(ElemKeyValueCI), pointer :: next
@ -2525,6 +2526,81 @@ contains
! We only need the mesh info for plotting
if (run_mode == MODE_PLOTTING) return
! ==========================================================================
! READ DATA FOR DERIVATIVES
! Get pointer list to XML <derivative>.
call get_node_list(doc, "derivative", node_deriv_list)
! Allocate TallyDerivative array.
allocate(tally_derivs(get_list_size(node_deriv_list)))
! Read derivative attributes.
do i = 1, get_list_size(node_deriv_list)
associate(deriv => tally_derivs(i))
! Get pointer to derivative node.
call get_list_item(node_deriv_list, i, node_deriv)
! Copy the derivative id.
if (check_for_node(node_deriv, "id")) then
call get_node_value(node_deriv, "id", deriv % id)
else
call fatal_error("Must specify an ID for <derivative> elements in the&
& tally XML file")
end if
! Make sure the id is > 0.
if (deriv % id <= 0) call fatal_error("<derivative> IDs must be an &
&integer greater than zero")
! Make sure this id has not already been used.
do j = 1, i-1
if (tally_derivs(j) % id == deriv % id) call fatal_error("Two or more&
& <derivative>'s use the same unique ID: " &
// trim(to_str(deriv % id)))
end do
! Read the dependent variable name.
temp_str = ""
call get_node_value(node_deriv, "variable", temp_str)
temp_str = to_lower(temp_str)
select case(temp_str)
case("density")
deriv % variable = DIFF_DENSITY
call get_node_value(node_deriv, "material", deriv % diff_material)
case("nuclide_density")
deriv % variable = DIFF_NUCLIDE_DENSITY
call get_node_value(node_deriv, "material", deriv % diff_material)
call get_node_value(node_deriv, "nuclide", word)
word = trim(to_lower(word))
pair_list => nuclide_dict % keys()
do while (associated(pair_list))
if (starts_with(pair_list % key, word)) then
word = pair_list % key(1:150)
exit
end if
! Advance to next
pair_list => pair_list % next
end do
! Check if no nuclide was found
if (.not. associated(pair_list)) then
call fatal_error("Could not find the nuclide " &
&// trim(word) // " specified in derivative " &
&// trim(to_str(deriv % id)) // " in any material.")
end if
deallocate(pair_list)
deriv % diff_nuclide = nuclide_dict % get_key(word)
end select
end associate
end do
! ==========================================================================
! READ TALLY DATA
@ -2982,52 +3058,6 @@ contains
end do
end if
! =======================================================================
! READ DATA FOR DERIVATIVES
if (check_for_node(node_tal, "derivative")) then
call get_node_ptr(node_tal, "derivative", node_deriv)
allocate(t % deriv)
temp_str = ""
call get_node_value(node_deriv, "variable", temp_str)
temp_str = to_lower(temp_str)
select case(temp_str)
case("density")
t % deriv % variable = DIFF_DENSITY
call get_node_value(node_deriv, "material", t % deriv % diff_material)
case("nuclide_density")
t % deriv % variable = DIFF_NUCLIDE_DENSITY
call get_node_value(node_deriv, "material", t % deriv % diff_material)
call get_node_value(node_deriv, "nuclide", word)
word = trim(to_lower(word))
pair_list => nuclide_dict % keys()
do while (associated(pair_list))
if (starts_with(pair_list % key, word)) then
word = pair_list % key(1:150)
exit
end if
! Advance to next
pair_list => pair_list % next
end do
! Check if no nuclide was found
if (.not. associated(pair_list)) then
call fatal_error("Could not find the nuclide " &
&// trim(word) // " specified in tally " &
&// trim(to_str(t % id)) // " in any material.")
end if
deallocate(pair_list)
t % deriv % diff_nuclide = nuclide_dict % get_key(word)
end select
end if
! =======================================================================
! READ DATA FOR SCORES
@ -3479,10 +3509,28 @@ contains
&// trim(to_str(t % id)) // ".")
end if
! If a derivative is present, we can only use analog or collision
! estimators.
if (allocated(t % deriv) .and. t % estimator == ESTIMATOR_TRACKLENGTH) &
t % estimator = ESTIMATOR_COLLISION
! Check for a tally derivative.
if (check_for_node(node_tal, "derivative")) then
! Temporarily store the derivative id.
call get_node_value(node_tal, "derivative", t % deriv)
! Find the derivative with the given id, and store it's index.
do j = 1, size(tally_derivs)
if (tally_derivs(j) % id == t % deriv) then
t % deriv = j
! Only analog or collision estimators are supported for differential
! tallies.
if (t % estimator == ESTIMATOR_TRACKLENGTH) &
t % estimator = ESTIMATOR_COLLISION
exit
end if
if (j == size(tally_derivs)) then
call fatal_error("Could not find derivative " &
// trim(to_str(t % deriv)) // " specified on tally " &
// trim(to_str(t % id)))
end if
end do
end if
! If settings.xml trigger is turned on, create tally triggers
if (trigger_on) then

30
src/output.F90
View file

@ -1017,21 +1017,21 @@ contains
endif
! Write derivative information.
if (allocated(t % deriv)) then
select case (t % deriv % variable)
case (DIFF_DENSITY)
write(unit=unit_tally, fmt="(' Density derivative Material ',A)") &
to_str(t % deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
i_listing = nuclides(t % deriv % diff_nuclide) % listing
write(unit=unit_tally, fmt="(' Nuclide density derivative Material '&
&,A,' Nuclide ',A)") trim(to_str(t % deriv % diff_material)), &
trim(xs_listings(i_listing) % alias)
case default
call fatal_error("Differential tally dependent variable for tally " &
// trim(to_str(t % id)) // " not defined in output.F90.")
end select
end if
!if (allocated(t % deriv)) then
! select case (t % deriv % variable)
! case (DIFF_DENSITY)
! write(unit=unit_tally, fmt="(' Density derivative Material ',A)") &
! to_str(t % deriv % diff_material)
! case (DIFF_NUCLIDE_DENSITY)
! i_listing = nuclides(t % deriv % diff_nuclide) % listing
! write(unit=unit_tally, fmt="(' Nuclide density derivative Material '&
! &,A,' Nuclide ',A)") trim(to_str(t % deriv % diff_material)), &
! trim(xs_listings(i_listing) % alias)
! case default
! call fatal_error("Differential tally dependent variable for tally " &
! // trim(to_str(t % id)) // " not defined in output.F90.")
! end select
!end if
! Handle surface current tallies separately
if (t % type == TALLY_SURFACE_CURRENT) then

2
src/simulation.F90
View file

@ -80,7 +80,7 @@ contains
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work

58
src/state_point.F90
View file

@ -52,7 +52,7 @@ contains
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group, tally_group
integer(HID_T) :: meshes_group, mesh_group
integer(HID_T) :: filter_group, deriv_group
integer(HID_T) :: filter_group, derivs_group, deriv_group
character(20), allocatable :: str_array(:)
character(MAX_FILE_LEN) :: filename
type(RegularMesh), pointer :: meshp
@ -194,6 +194,35 @@ contains
call close_group(meshes_group)
! Write information for derivatives.
if (size(tally_derivs) > 0) then
derivs_group = create_group(tallies_group, "derivatives")
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
deriv_group = create_group(derivs_group, "derivative " &
// trim(to_str(deriv % id)))
select case (deriv % variable)
case (DIFF_DENSITY)
call write_dataset(deriv_group, "dependent variable", "density")
call write_dataset(deriv_group, "material", deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
call write_dataset(deriv_group, "dependent variable", &
"nuclide_density")
call write_dataset(deriv_group, "material", deriv % diff_material)
i_list = nuclides(deriv % diff_nuclide) % listing
call write_dataset(deriv_group, "nuclide", &
xs_listings(i_list) % alias)
case default
call fatal_error("Dependent variable for derivative " &
// trim(to_str(deriv % id)) // " not defined in &
&state_point.F90.")
end select
call close_group(deriv_group)
end associate
end do
call close_group(derivs_group)
end if
! Write number of tallies
call write_dataset(tallies_group, "n_tallies", n_tallies)
@ -308,30 +337,9 @@ contains
deallocate(str_array)
! Write derivative information.
if (allocated(tally % deriv)) then
call write_dataset(tally_group, "derivative present", 1)
deriv_group = create_group(tally_group, "derivative")
select case (tally % deriv % variable)
case (DIFF_DENSITY)
call write_dataset(deriv_group, "dependent variable", "density")
call write_dataset(deriv_group, "material", &
tally % deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
call write_dataset(deriv_group, "dependent variable", &
"nuclide_density")
call write_dataset(deriv_group, "material", &
tally % deriv % diff_material)
i_list = nuclides(tally % deriv % diff_nuclide) % listing
call write_dataset(deriv_group, "nuclide", &
xs_listings(i_list) % alias)
case default
call fatal_error("Differential tally dependent variable for &
&tally " // trim(to_str(tally % id)) // " not defined in &
&state_point.F90.")
end select
call close_group(deriv_group)
else
call write_dataset(tally_group, "derivative present", 0)
if (tally % deriv /= NONE) then
call write_dataset(tally_group, "derivative", &
tally_derivs(tally % deriv) % id)
end if
! Write scores.

324
src/tally.F90
View file

@ -746,7 +746,7 @@ contains
!#########################################################################
! Add derivative information on score for differential tallies.
if (allocated(t % deriv)) then
if (t % deriv /= NONE) then
call apply_derivative_to_score(p, t, i_nuclide, atom_density, &
score_bin, score)
end if
@ -1049,7 +1049,7 @@ contains
! Add derivative information for differenetial tallies. Note that the
! i_nuclide and atom_density arguments do not matter since this is an
! analog estimator.
if (allocated(t % deriv)) then
if (t % deriv /= NONE) then
call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score)
end if
@ -2254,166 +2254,163 @@ contains
integer :: l ! loop index for nuclides in material
logical :: scoring_diff_nuclide
select case (t % deriv % variable)
associate(deriv => tally_derivs(t % deriv))
select case (deriv % variable)
case (DIFF_DENSITY)
select case (t % estimator)
case (DIFF_DENSITY)
select case (t % estimator)
case (ESTIMATOR_ANALOG)
if (materials(p % material) % id == t % deriv % diff_material) then
score = score * (t % deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * t % deriv % flux_deriv
end if
case (ESTIMATOR_COLLISION)
select case (score_bin)
case (SCORE_FLUX)
score = score * t % deriv % flux_deriv
case (SCORE_TOTAL)
if (materials(p % material) % id == t % deriv % diff_material &
.and. material_xs % total /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE &
case (ESTIMATOR_ANALOG)
if (materials(p % material) % id == deriv % diff_material) then
score = score * (deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * t % deriv % flux_deriv
score = score * deriv % flux_deriv
end if
case (SCORE_ABSORPTION)
if (materials(p % material) % id == t % deriv % diff_material &
.and. material_xs % absorption /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * t % deriv % flux_deriv
end if
case (ESTIMATOR_COLLISION)
case (SCORE_FISSION)
if (materials(p % material) % id == t % deriv % diff_material &
.and. material_xs % fission /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * t % deriv % flux_deriv
end if
select case (score_bin)
case (SCORE_NU_FISSION)
if (materials(p % material) % id == t % deriv % diff_material &
.and. material_xs % nu_fission /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * t % deriv % flux_deriv
end if
case (SCORE_FLUX)
score = score * deriv % flux_deriv
case (SCORE_TOTAL)
if (materials(p % material) % id == deriv % diff_material &
.and. material_xs % total /= ZERO) then
score = score * (deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * deriv % flux_deriv
end if
case (SCORE_ABSORPTION)
if (materials(p % material) % id == deriv % diff_material &
.and. material_xs % absorption /= ZERO) then
score = score * (deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * deriv % flux_deriv
end if
case (SCORE_FISSION)
if (materials(p % material) % id == deriv % diff_material &
.and. material_xs % fission /= ZERO) then
score = score * (deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * deriv % flux_deriv
end if
case (SCORE_NU_FISSION)
if (materials(p % material) % id == deriv % diff_material &
.and. material_xs % nu_fission /= ZERO) then
score = score * (deriv % flux_deriv + ONE &
/ materials(p % material) % density_gpcc)
else
score = score * deriv % flux_deriv
end if
case default
call fatal_error('Tally derivative not defined for a score on &
&tally ' // trim(to_str(t % id)))
end select
case default
call fatal_error('Tally derivative not defined for a score on &
&tally ' // trim(to_str(t % id)))
call fatal_error("Differential tallies are only implemented for &
&analog and collision estimators.")
end select
case default
call fatal_error("Differential tallies are only implemented for &
&analog and collision estimators.")
case (DIFF_NUCLIDE_DENSITY)
select case (t % estimator)
end select
case (DIFF_NUCLIDE_DENSITY)
select case (t % estimator)
case (ESTIMATOR_ANALOG)
if (materials(p % material) % id == t % deriv % diff_material &
.and. p % event_nuclide == t % deriv % diff_nuclide) then
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
if (mat % nuclide(l) == t % deriv % diff_nuclide) exit
end do
score = score * (t % deriv % flux_deriv &
+ ONE / mat % atom_density(l))
end associate
else
score = score * t % deriv % flux_deriv
end if
case (ESTIMATOR_COLLISION)
scoring_diff_nuclide = &
(materials(p % material) % id == t % deriv % diff_material) &
.and. (i_nuclide == t % deriv % diff_nuclide)
select case (score_bin)
case (SCORE_FLUX)
score = score * t % deriv % flux_deriv
case (SCORE_TOTAL)
if (i_nuclide == -1 .and. &
materials(p % material) % id == t % deriv % diff_material) then
score = score * (t % deriv % flux_deriv &
+ micro_xs(t % deriv % diff_nuclide) % total &
/ material_xs % total)
else if (scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % total /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE / atom_density)
case (ESTIMATOR_ANALOG)
if (materials(p % material) % id == deriv % diff_material &
.and. p % event_nuclide == deriv % diff_nuclide) then
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
if (mat % nuclide(l) == deriv % diff_nuclide) exit
end do
score = score * (deriv % flux_deriv &
+ ONE / mat % atom_density(l))
end associate
else
score = score * t % deriv % flux_deriv
score = score * deriv % flux_deriv
end if
case (SCORE_ABSORPTION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == t % deriv % diff_material) then
score = score * (t % deriv % flux_deriv &
+ micro_xs(t % deriv % diff_nuclide) % absorption &
/ material_xs % absorption )
else if (scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % absorption /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE / atom_density)
else
score = score * t % deriv % flux_deriv
end if
case (ESTIMATOR_COLLISION)
scoring_diff_nuclide = &
(materials(p % material) % id == deriv % diff_material) &
.and. (i_nuclide == deriv % diff_nuclide)
case (SCORE_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == t % deriv % diff_material) then
score = score * (t % deriv % flux_deriv &
+ micro_xs(t % deriv % diff_nuclide) % fission &
/ material_xs % fission)
else if (scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % fission /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE / atom_density)
else
score = score * t % deriv % flux_deriv
end if
select case (score_bin)
case (SCORE_NU_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == t % deriv % diff_material) then
score = score * (t % deriv % flux_deriv &
+ micro_xs(t % deriv % diff_nuclide) % nu_fission &
/ material_xs % nu_fission)
else if (scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO) then
score = score * (t % deriv % flux_deriv + ONE / atom_density)
else
score = score * t % deriv % flux_deriv
end if
case (SCORE_FLUX)
score = score * deriv % flux_deriv
case (SCORE_TOTAL)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material) then
score = score * (deriv % flux_deriv &
+ micro_xs(deriv % diff_nuclide) % total &
/ material_xs % total)
else if (scoring_diff_nuclide .and. &
micro_xs(deriv % diff_nuclide) % total /= ZERO) then
score = score * (deriv % flux_deriv + ONE / atom_density)
else
score = score * deriv % flux_deriv
end if
case (SCORE_ABSORPTION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material) then
score = score * (deriv % flux_deriv &
+ micro_xs(deriv % diff_nuclide) % absorption &
/ material_xs % absorption )
else if (scoring_diff_nuclide .and. &
micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then
score = score * (deriv % flux_deriv + ONE / atom_density)
else
score = score * deriv % flux_deriv
end if
case (SCORE_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material) then
score = score * (deriv % flux_deriv &
+ micro_xs(deriv % diff_nuclide) % fission &
/ material_xs % fission)
else if (scoring_diff_nuclide .and. &
micro_xs(deriv % diff_nuclide) % fission /= ZERO) then
score = score * (deriv % flux_deriv + ONE / atom_density)
else
score = score * deriv % flux_deriv
end if
case (SCORE_NU_FISSION)
if (i_nuclide == -1 .and. &
materials(p % material) % id == deriv % diff_material) then
score = score * (deriv % flux_deriv &
+ micro_xs(deriv % diff_nuclide) % nu_fission &
/ material_xs % nu_fission)
else if (scoring_diff_nuclide .and. &
micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then
score = score * (deriv % flux_deriv + ONE / atom_density)
else
score = score * deriv % flux_deriv
end if
case default
call fatal_error('Tally derivative not defined for a score on &
&tally ' // trim(to_str(t % id)))
end select
case default
call fatal_error('Tally derivative not defined for a score on &
&tally ' // trim(to_str(t % id)))
call fatal_error("Differential tallies are only implemented for &
&analog and collision estimators.")
end select
case default
call fatal_error("Differential tallies are only implemented for &
&analog and collision estimators.")
end select
end select
end associate
end subroutine apply_derivative_to_score
!===============================================================================
@ -2427,33 +2424,32 @@ contains
integer :: i
! A void material cannot be perturbed so it will not affect flux derivatives
if (p % material == MATERIAL_VOID) return
do i = 1, active_tallies % size()
associate (t => tallies(active_tallies % get_item(i)))
if (.not. allocated(t % deriv)) cycle ! Ignore non-differential tallies
select case (t % deriv % variable)
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
select case (deriv % variable)
case (DIFF_DENSITY)
associate (mat => materials(p % material))
if (mat % id == t % deriv % diff_material) then
if (mat % id == deriv % diff_material) then
! phi = e^(-Sigma_tot * dist)
! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
t % deriv % flux_deriv = t % deriv % flux_deriv &
deriv % flux_deriv = deriv % flux_deriv &
- distance * material_xs % total / mat % density_gpcc
end if
end associate
case (DIFF_NUCLIDE_DENSITY)
associate (mat => materials(p % material))
if (mat % id == t % deriv % diff_material) then
if (mat % id == deriv % diff_material) then
! phi = e^(-Sigma_tot * dist)
! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
t % deriv % flux_deriv = t % deriv % flux_deriv &
- distance * micro_xs(t % deriv % diff_nuclide) % total
deriv % flux_deriv = deriv % flux_deriv &
- distance * micro_xs(deriv % diff_nuclide) % total
end if
end associate
end select
@ -2471,41 +2467,41 @@ contains
integer :: i, j
! A void material cannot be perturbed so it will not affect flux derivatives
if (p % material == MATERIAL_VOID) return
do i = 1, active_tallies % size()
associate (t => tallies(active_tallies % get_item(i)))
if (.not. allocated(t % deriv)) cycle ! Ignore non-differential tallies
select case (t % deriv % variable)
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
select case (deriv % variable)
case (DIFF_DENSITY)
associate (mat => materials(p % material))
if (mat % id == t % deriv % diff_material) then
if (mat % id == deriv % diff_material) then
! phi = Sigma_MT
! (1 / phi) * (d_phi / d_rho) = (d_Sigma_MT / d_rho) / Sigma_MT
! (1 / phi) * (d_phi / d_rho) = 1 / rho
t % deriv % flux_deriv = t % deriv % flux_deriv &
deriv % flux_deriv = deriv % flux_deriv &
+ ONE / mat % density_gpcc
end if
end associate
case (DIFF_NUCLIDE_DENSITY)
associate (mat => materials(p % material))
if (mat % id == t % deriv % diff_material &
.and. p % event_nuclide == t % deriv % diff_nuclide) then
if (mat % id == deriv % diff_material &
.and. p % event_nuclide == deriv % diff_nuclide) then
! Find the index in this material for the diff_nuclide.
do j = 1, mat % n_nuclides
if (mat % nuclide(j) == t % deriv % diff_nuclide) exit
if (mat % nuclide(j) == deriv % diff_nuclide) exit
end do
! Make sure we found the nuclide.
if (mat % nuclide(j) /= t % deriv % diff_nuclide) then
if (mat % nuclide(j) /= deriv % diff_nuclide) then
call fatal_error("Couldn't find the right nuclide.")
end if
! phi = Sigma_MT
! (1 / phi) * (d_phi / d_N) = (d_Sigma_MT / d_N) / Sigma_MT
! (1 / phi) * (d_phi / d_N) = sigma_MT / Sigma_MT
! (1 / phi) * (d_phi / d_N) = 1 / N
t % deriv % flux_deriv = t % deriv % flux_deriv &
deriv % flux_deriv = deriv % flux_deriv &
+ ONE / mat % atom_density(j)
end if
end associate
@ -2520,10 +2516,8 @@ contains
subroutine zero_flux_derivs()
integer :: i
do i = 1, n_tallies
if (allocated(tallies(i) % deriv)) then
tallies(i) % deriv % flux_deriv = ZERO
end if
do i = 1, size(tally_derivs)
tally_derivs(i) % flux_deriv = ZERO
end do
end subroutine zero_flux_derivs

5
src/tally_header.F90
View file

@ -66,6 +66,7 @@ module tally_header
!===============================================================================
type TallyDerivative
integer :: id
real(8) :: flux_deriv
integer :: variable
integer :: diff_material
@ -138,8 +139,8 @@ module tally_header
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Derivative for differentially tallies
type(TallyDerivative), allocatable :: deriv
! Index for the TallyDerivative for differential tallies.
integer :: deriv = NONE
end type TallyObject
end module tally_header

6
src/tracking.F90
View file

@ -63,7 +63,7 @@ contains
endif
! Every particle starts with no accumulated flux derivative.
call zero_flux_derivs()
if (active_tallies % size() > 0) call zero_flux_derivs()
EVENT_LOOP: do
! If the cell hasn't been determined based on the particle's location,
@ -120,7 +120,7 @@ contains
end if
! Score flux derivative accumulators for differential tallies.
call score_track_derivative(p, distance)
if (active_tallies % size() > 0) call score_track_derivative(p, distance)
if (d_collision > d_boundary) then
! ====================================================================
@ -197,7 +197,7 @@ contains
end do
! Score flux derivative accumulators for differential tallies.
call score_collision_derivative(p)
if (active_tallies % size() > 0) call score_collision_derivative(p)
end if
! If particle has too many events, display warning and kill it

2
tests/test_diff_tally/inputs_true.dat
View file

@ -1 +1 @@
cb132b2211831bf333bc9be0b10d092ad3bbe9a188743cf4c430ae97037987763b268997757a0602c164d8c2deb3eb9843590d860ba5a4858092a9ad0b8adeba
e398b5d359cab6a448e1273beaf2d5c9125d4ad06546ddfd327111b101ac223e11b6f98d3b490b8544972969f37456c24a78016ce68773abba27828ba43e3173

49
tests/test_diff_tally/test_diff_tally.py
View file

@ -10,8 +10,8 @@ except:
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
from openmc import Filter, Mesh, Tally, TalliesFile, Summary, StatePoint
#from openmc.statepoint import StatePoint
from openmc import Filter, Mesh, Tally, TalliesFile, Summary, StatePoint, \
TallyDerivative
from openmc.source import Source
from openmc.stats import Box
@ -33,26 +33,35 @@ class DiffTallyTestHarness(PyAPITestHarness):
filt_mats = Filter(type='material', bins=(1, 3))
filt_eout = Filter(type='energyout', bins=(0.0, 1.0, 20.0))
# We want density derivatives for both water and fuel to get coverage
# for both fissile and non-fissile materials.
d1 = TallyDerivative(derivative_id=1)
d1.variable = 'density'
d1.material = 3
d2 = TallyDerivative(derivative_id=2)
d2.variable = 'density'
d2.material = 1
# O-16 is a good nuclide to test against because it is present in both
# water and fuel. Some routines need to recognize that they have the
# perturbed nuclide but not the perturbed material.
d3 = TallyDerivative(derivative_id=3)
d3.variable = 'nuclide_density'
d3.material = 1
d3.nuclide = 'O-16'
# A fissile nuclide, just for good measure.
d4 = TallyDerivative(derivative_id=4)
d4.variable = 'nuclide_density'
d4.material = 1
d4.nuclide = 'U-235'
def add_derivs(tally_list):
assert len(tally_list) == 4
# We want density derivatives for both water and fuel to get
# coverage for both fissile and non-fissile materials.
tally_list[0].diff_variable = 'density'
tally_list[0].diff_material = 3
tally_list[1].diff_variable = 'density'
tally_list[1].diff_material = 1
# O-16 is a good nuclide to test against because it is present
# in both water and fuel. Some routines need to recognize that they
# have the perturbed nuclide but not the perturbed material.
tally_list[2].diff_variable = 'nuclide_density'
tally_list[2].diff_material = 1
tally_list[2].diff_nuclide = 'O-16'
# A fissile nuclide, just for good measure.
tally_list[3].diff_variable = 'nuclide_density'
tally_list[3].diff_material = 1
tally_list[3].diff_nuclide = 'U-235'
tally_list[0].derivative = d1
tally_list[1].derivative = d2
tally_list[2].derivative = d3
tally_list[3].derivative = d4
# Cover the flux score.
tallies = [Tally() for i in range(4)]

2
tests/testing_harness.py
View file

@ -26,7 +26,7 @@ class TestHarness(object):
self.parser = OptionParser()
self.parser.add_option('--exe', dest='exe', default='openmc')
self.parser.add_option('--mpi_exec', dest='mpi_exec', default=None)
self.parser.add_option('--mpi_np', dest='mpi_np', type=int, default=3)
self.parser.add_option('--mpi_np', dest='mpi_np', type=int, default=1)
self.parser.add_option('--update', dest='update', action='store_true',
default=False)
self._opts = None