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

Merge remote-tracking branch 'upstream/develop' into fission_q

Browse source
This commit is contained in:
Sterling Harper 2016-08-03 10:18:18 -05:00
commit 3b8b55ab8d
54 changed files with 4286 additions and 2744 deletions
Download patch file Download diff file Expand all files Collapse all files

13
docs/source/pythonapi/examples/mgxs-part-i.ipynb
View file

@ -383,6 +383,9 @@
"* `ScatterMatrixXS`\n",
"* `NuScatterMatrixXS`\n",
"* `Chi`\n",
"* `ChiPrompt`\n",
"* `InverseVelocity`\n",
"* `PromptNuFissionXS`\n",
"\n",
"These classes provide us with an interface to generate the tally inputs as well as perform post-processing of OpenMC's tally data to compute the respective multi-group cross sections. In this case, let's create the multi-group total, absorption and scattering cross sections with our 2-group structure."
]
@ -1164,21 +1167,21 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 2",
"language": "python",
"name": "python3"
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,

15
docs/source/pythonapi/examples/mgxs-part-iii.ipynb
View file

@ -549,6 +549,9 @@
"* `ScatterMatrixXS` (`\"scatter matrix\"`)\n",
"* `NuScatterMatrixXS` (`\"nu-scatter matrix\"`)\n",
"* `Chi` (`\"chi\"`)\n",
"* `ChiPrompt` (`\"chi prompt\"`)\n",
"* `InverseVelocity` (`\"inverse-velocity\"`)\n",
"* `PromptNuFissionXS` (`\"prompt-nu-fission\"`)\n",
"\n",
"In this case, let's create the multi-group cross sections needed to run an OpenMOC simulation to verify the accuracy of our cross sections. In particular, we will define `\"transport\"`, `\"nu-fission\"`, `'\"fission\"`, `\"nu-scatter matrix\"` and `\"chi\"` cross sections for our `Library`.\n",
"\n",
@ -571,7 +574,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports `\"material,\"` `\"cell,\"` and `\"universe\"` domain types. We will use a `\"cell\"` domain type here to compute cross sections in each of the cells in the fuel assembly geometry.\n",
"Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports `\"material\"`, `\"cell\"`, `\"universe\"`, and `\"mesh\"` domain types. We will use a `\"cell\"` domain type here to compute cross sections in each of the cells in the fuel assembly geometry.\n",
"\n",
"**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell or universe) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property. In our case, we wish to compute multi-group cross sections in each and every cell since they will be needed in our downstream OpenMOC calculation on the identical combinatorial geometry mesh."
]
@ -1596,21 +1599,21 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 2",
"language": "python",
"name": "python3"
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,

14
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

@ -519,9 +519,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports \"material,\" \"cell,\" and \"universe\" domain types. In this simple example, we wish to compute multi-group cross sections only for each material andtherefore will use a \"material\" domain type.\n",
"Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports \"material\" \"cell\", \"universe\", and \"mesh\" domain types. In this simple example, we wish to compute multi-group cross sections only for each material and therefore will use a \"material\" domain type.\n",
"\n",
"**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell or universe) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property."
"**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell, universe, or mesh) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property."
]
},
{
@ -1437,21 +1437,21 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 2",
"language": "python",
"name": "python3"
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,

506
docs/source/pythonapi/examples/nuclear-data.ipynb
View file

File diff suppressed because one or more lines are too long

3
docs/source/pythonapi/index.rst
View file

@ -269,13 +269,16 @@ Multi-group Cross Sections
openmc.mgxs.AbsorptionXS
openmc.mgxs.CaptureXS
openmc.mgxs.Chi
openmc.mgxs.ChiPrompt
openmc.mgxs.FissionXS
openmc.mgxs.InverseVelocity
openmc.mgxs.KappaFissionXS
openmc.mgxs.MultiplicityMatrixXS
openmc.mgxs.NuFissionXS
openmc.mgxs.NuFissionMatrixXS
openmc.mgxs.NuScatterXS
openmc.mgxs.NuScatterMatrixXS
openmc.mgxs.PromptNuFissionXS
openmc.mgxs.ScatterXS
openmc.mgxs.ScatterMatrixXS
openmc.mgxs.TotalXS

4
docs/source/usersguide/input.rst
View file

@ -1764,6 +1764,10 @@ The ``<tally>`` element accepts the following sub-elements:
| |fission. This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
|prompt-nu-fission |Total production of prompt neutrons due to |
| |fission. This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
|nu-scatter, |These scores are similar in functionality to their |

2
openmc/__init__.py
View file

@ -9,8 +9,8 @@ from openmc.plots import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
from openmc.mgxs_library import *
from openmc.mesh import *
from openmc.mgxs_library import *
from openmc.filter import *
from openmc.trigger import *
from openmc.tallies import *

2
openmc/checkvalue.py
View file

@ -212,7 +212,7 @@ def check_less_than(name, value, maximum, equality=False):
raise ValueError(msg)
def check_greater_than(name, value, minimum, equality=False):
"""Ensure that an object's value is less than a given value.
"""Ensure that an object's value is greater than a given value.
Parameters
----------

9
openmc/data/data.py
View file

@ -122,10 +122,11 @@ ATOMIC_SYMBOL = {1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N',
114: 'Fl', 116: 'Lv'}
ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)', 5: '(n,misc)', 11: '(n,2nd)',
16: '(n,2n)', 17: '(n,3n)', 18: '(n,fission)', 19: '(n,f)',
20: '(n,nf)', 21: '(n,2nf)', 22: '(n,na)', 23: '(n,n3a)',
24: '(n,2na)', 25: '(n,3na)', 28: '(n,np)', 29: '(n,n2a)',
REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
5: '(n,misc)', 11: '(n,2nd)', 16: '(n,2n)', 17: '(n,3n)',
18: '(n,fission)', 19: '(n,f)', 20: '(n,nf)', 21: '(n,2nf)',
22: '(n,na)', 23: '(n,n3a)', 24: '(n,2na)', 25: '(n,3na)',
27: '(n,absorption)', 28: '(n,np)', 29: '(n,n2a)',
30: '(n,2n2a)', 32: '(n,nd)', 33: '(n,nt)', 34: '(n,nHe-3)',
35: '(n,nd2a)', 36: '(n,nt2a)', 37: '(n,4n)', 38: '(n,3nf)',
41: '(n,2np)', 42: '(n,3np)', 44: '(n,n2p)', 45: '(n,npa)',

48
openmc/data/function.py
View file

@ -80,50 +80,46 @@ class Tabulated1D(object):
# Get indices for interpolation
idx = np.searchsorted(self.x, x, side='right') - 1
# Find lowest valid index
i_low = np.searchsorted(idx, 0)
# Loop over interpolation regions
for k in range(len(self.breakpoints)):
# Determine which x values are within this interpolation range
i_high = np.searchsorted(idx, self.breakpoints[k] - 1)
# Get indices for the begining and ending of this region
i_begin = self.breakpoints[k-1] - 1 if k > 0 else 0
i_end = self.breakpoints[k] - 1
# Get x values and bounding (x,y) pairs
xk = x[i_low:i_high]
xi = self.x[idx[i_low:i_high]]
xi1 = self.x[idx[i_low:i_high] + 1]
yi = self.y[idx[i_low:i_high]]
yi1 = self.y[idx[i_low:i_high] + 1]
# Figure out which idx values lie within this region
contained = (idx >= i_begin) & (idx < i_end)
xk = x[contained] # x values in this region
xi = self.x[idx[contained]] # low edge of corresponding bins
xi1 = self.x[idx[contained] + 1] # high edge of corresponding bins
yi = self.y[idx[contained]]
yi1 = self.y[idx[contained] + 1]
if self.interpolation[k] == 1:
# Histogram
y[i_low:i_high] = yi
y[contined] = yi
elif self.interpolation[k] == 2:
# Linear-linear
y[i_low:i_high] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
y[contained] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
elif self.interpolation[k] == 3:
# Linear-log
y[i_low:i_high] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
y[contained] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
elif self.interpolation[k] == 4:
# Log-linear
y[i_low:i_high] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
y[contained] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
elif self.interpolation[k] == 5:
# Log-log
y[i_low:i_high] = yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)*np.log(yi1/yi))
y[contained] = (yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)
*np.log(yi1/yi)))
i_low = i_high
# In some cases, the first/last point of x may be less than the first
# value of self.x due only to precision, so we check if they're close
# and set them equal if so. Otherwise, the interpolated value might be
# out of range (and thus zero)
if np.isclose(x[0], self.x[0], 1e-8):
y[0] = self.y[0]
if np.isclose(x[-1], self.x[-1], 1e-8):
y[-1] = self.y[-1]
# In some cases, x values might be outside the tabulated region due only
# to precision, so we check if they're close and set them equal if so.
y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]
return y if iterable else y[0]

11
openmc/data/neutron.py
View file

@ -346,6 +346,17 @@ class IncidentNeutron(object):
tgroup = group['total_nu']
rx.derived_products.append(Product.from_hdf5(tgroup))
# Build summed reactions. Start from the highest MT number because high
# MTs never depend on lower MTs.
for mt_sum in sorted(SUM_RULES, reverse=True):
if mt_sum not in data:
xs_components = [data[mt].xs for mt in SUM_RULES[mt_sum]
if mt in data]
if len(xs_components) > 0:
rxn = Reaction(mt_sum)
rxn.xs = Sum(xs_components)
data.summed_reactions[mt_sum] = rxn
# Read unresolved resonance probability tables
if 'urr' in group:
urr_group = group['urr']

2
openmc/data/reaction.py
View file

@ -400,7 +400,7 @@ class Reaction(object):
# Read cross section
if 'xs' in group:
xs = group['xs'].value
rx.xs = Tabulated1D(energy, xs)
rx.xs = Tabulated1D(energy[rx.threshold_idx:], xs)
# Determine number of products
n_product = 0

7
openmc/filter.py
View file

@ -162,6 +162,11 @@ class Filter(object):
if not isinstance(bins, Iterable):
bins = [bins]
# If the bin is 0D numpy array, promote to 1D
elif isinstance(bins, np.ndarray):
if bins.shape == ():
bins.shape = (1,)
# If the bins are in a collection, convert it to a list
else:
bins = list(bins)
@ -563,7 +568,7 @@ class Filter(object):
# Initialize dictionary to build Pandas Multi-index column
filter_dict = {}
# Append Mesh ID as outermost index of mult-index
# Append Mesh ID as outermost index of multi-index
mesh_key = 'mesh {0}'.format(self.mesh.id)
# Find mesh dimensions - use 3D indices for simplicity

16
openmc/mgxs/library.py
View file

@ -54,9 +54,9 @@ class Library(object):
If true, computes cross sections for each nuclide in each domain
mgxs_types : Iterable of str
The types of cross sections in the library (e.g., ['total', 'scatter'])
domain_type : {'material', 'cell', 'distribcell', 'universe'}
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
domains : Iterable of openmc.Material, openmc.Cell or openmc.Universe
domains : Iterable of openmc.Material, openmc.Cell, openmc.Universe or openmc.Mesh
The spatial domain(s) for which MGXS in the Library are computed
correction : {'P0', None}
Apply the P0 correction to scattering matrices if set to 'P0'
@ -183,6 +183,8 @@ class Library(object):
return self.openmc_geometry.get_all_material_cells()
elif self.domain_type == 'universe':
return self.openmc_geometry.get_all_universes()
elif self.domain_type == 'mesh':
raise ValueError('Unable to get domains for Mesh domain type')
else:
raise ValueError('Unable to get domains without a domain type')
else:
@ -273,6 +275,12 @@ class Library(object):
elif self.domain_type == 'universe':
cv.check_iterable_type('domain', domains, openmc.Universe)
all_domains = self.openmc_geometry.get_all_universes()
elif self.domain_type == 'mesh':
cv.check_iterable_type('domain', domains, openmc.Mesh)
# The mesh and geometry are independent, so set all_domains
# to the input domains
all_domains = domains
else:
raise ValueError('Unable to set domains with domain '
'type "{}"'.format(self.domain_type))
@ -452,7 +460,7 @@ class Library(object):
----------
domain : Material or Cell or Universe or Integral
The material, cell, or universe object of interest (or its ID)
mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi'}
mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission'}
The type of multi-group cross section object to return
Returns
@ -474,6 +482,8 @@ class Library(object):
cv.check_type('domain', domain, (openmc.Cell, Integral))
elif self.domain_type == 'universe':
cv.check_type('domain', domain, (openmc.Universe, Integral))
elif self.domain_type == 'mesh':
cv.check_type('domain', domain, (openmc.Mesh, Integral))
# Check that requested domain is included in library
if isinstance(domain, Integral):

762
openmc/mgxs/mgxs.py
View file

File diff suppressed because it is too large Load diff

2
openmc/tallies.py
View file

@ -1293,7 +1293,7 @@ class Tally(object):
# Create list of 2- or 3-tuples tuples for mesh cell bins
if self_filter.type == 'mesh':
dimension = self_filter.mesh.dimension
xyz = map(lambda x: np.arange(1, x+1), dimension)
xyz = [range(1, x+1) for x in dimension]
bins = list(itertools.product(*xyz))
# Create list of 2-tuples for energy boundary bins

59
src/cmfd_data.F90
View file

@ -6,6 +6,7 @@ module cmfd_data
!==============================================================================
use constants
use tally_filter, only: MeshFilter
implicit none
private
@ -94,7 +95,10 @@ contains
! Associate tallies and mesh
t => cmfd_tallies(1)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
select type(filt => t % filters(t % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
m => meshes(i_mesh)
! Set mesh widths
@ -109,7 +113,10 @@ contains
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
select type(filt => t % filters(t % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
m => meshes(i_mesh)
i_filter_mesh = t % find_filter(FILTER_MESH)
@ -138,7 +145,7 @@ contains
TALLY: if (ital == 1) then
! Reset all bins to 1
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! Set ijk as mesh indices
ijk = (/ i, j, k /)
@ -152,7 +159,8 @@ contains
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1
! Get flux
flux = t % results(1,score_index) % sum
@ -181,7 +189,7 @@ contains
INGROUP: do g = 1, ng
! Reset all bins to 1
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! Set ijk as mesh indices
ijk = (/ i, j, k /)
@ -198,7 +206,8 @@ contains
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1
! Get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
@ -220,7 +229,7 @@ contains
else if (ital == 3) then
! Initialize and filter for energy
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
if (i_filter_ein > 0) then
matching_bins(i_filter_ein) = ng - h + 1
end if
@ -229,60 +238,72 @@ contains
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
! Right surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
! Back surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
! Front surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
! Bottom surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
! Top surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
end if TALLY

97
src/cmfd_input.F90
View file

@ -267,14 +267,16 @@ contains
use mesh_header, only: RegularMesh
use string
use tally, only: setup_active_cmfdtallies
use tally_header, only: TallyObject, TallyFilter
use tally_header, only: TallyObject
use tally_filter_header
use tally_filter
use tally_initialize, only: add_tallies
use xml_interface
type(Node), pointer :: doc ! pointer to XML doc info
character(MAX_LINE_LEN) :: temp_str ! temp string
integer :: i ! loop counter
integer :: i, j ! loop counter
integer :: n ! size of arrays in mesh specification
integer :: ng ! number of energy groups (default 1)
integer :: n_filters ! number of filters
@ -283,7 +285,7 @@ contains
real(8) :: rarray3(3) ! temp double array
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(TallyFilter) :: filters(N_FILTER_TYPES) ! temporary filters
type(TallyFilterContainer) :: filters(N_FILTER_TYPES) ! temporary filters
type(Node), pointer :: node_mesh
! Set global variables if they are 0 (this can happen if there is no tally
@ -410,21 +412,25 @@ contains
! Set up mesh filter
n_filters = 1
filters(n_filters) % type = FILTER_MESH
filters(n_filters) % n_bins = product(m % dimension)
allocate(filters(n_filters) % int_bins(1))
filters(n_filters) % int_bins(1) = n_user_meshes + 1
allocate(MeshFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (MeshFilter)
filt % n_bins = product(m % dimension)
filt % mesh = n_user_meshes + 1
end select
t % find_filter(FILTER_MESH) = n_filters
! Read and set incoming energy mesh filter
if (check_for_node(node_mesh, "energy")) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYIN
ng = get_arraysize_double(node_mesh, "energy")
filters(n_filters) % n_bins = ng - 1
allocate(filters(n_filters) % real_bins(ng))
call get_node_array(node_mesh, "energy", &
filters(n_filters) % real_bins)
allocate(EnergyFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (EnergyFilter)
ng = get_arraysize_double(node_mesh, "energy")
filt % n_bins = ng - 1
allocate(filt % bins(ng))
call get_node_array(node_mesh, "energy", filt % bins)
end select
t % find_filter(FILTER_ENERGYIN) = n_filters
end if
@ -448,9 +454,10 @@ contains
t % type = TALLY_VOLUME
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate scoring bins
allocate(t % score_bins(3))
@ -481,23 +488,22 @@ contains
! read and set outgoing energy mesh filter
if (check_for_node(node_mesh, "energy")) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYOUT
ng = get_arraysize_double(node_mesh, "energy")
filters(n_filters) % n_bins = ng - 1
allocate(filters(n_filters) % real_bins(ng))
call get_node_array(node_mesh, "energy", &
filters(n_filters) % real_bins)
allocate(EnergyoutFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (EnergyoutFilter)
ng = get_arraysize_double(node_mesh, "energy")
filt % n_bins = ng - 1
allocate(filt % bins(ng))
call get_node_array(node_mesh, "energy", filt % bins)
end select
t % find_filter(FILTER_ENERGYOUT) = n_filters
end if
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! deallocate filters bins array
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(n_filters) % real_bins)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate macro reactions
allocate(t % score_bins(2))
@ -522,25 +528,25 @@ contains
! Add extra filter for surface
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_SURFACE
filters(n_filters) % n_bins = 2 * m % n_dimension
allocate(filters(n_filters) % int_bins(2 * m % n_dimension))
if (m % n_dimension == 2) then
filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT /)
elseif (m % n_dimension == 3) then
filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP /)
end if
allocate(SurfaceFilter :: filters(n_filters) % obj)
select type(filt => filters(n_filters) % obj)
type is(SurfaceFilter)
filt % n_bins = 2 * m % n_dimension
allocate(filt % surfaces(2 * m % n_dimension))
if (m % n_dimension == 2) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT, &
IN_TOP, OUT_TOP /)
end if
end select
t % find_filter(FILTER_SURFACE) = n_filters
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! Deallocate filters bins array
deallocate(filters(n_filters) % int_bins)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate macro reactions
allocate(t % score_bins(1))
@ -558,15 +564,10 @@ contains
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
i_filter_mesh = t % find_filter(FILTER_MESH)
t % filters(i_filter_mesh) % n_bins = product(m % dimension + 1)
t % filters(i_filter_mesh) % obj % n_bins = product(m % dimension + 1)
end if
! Deallocate filter bins
deallocate(filters(1) % int_bins)
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(2) % real_bins)
end do
! Put cmfd tallies into active tally array and turn tallies on

9
src/constants.F90
View file

@ -289,7 +289,7 @@ module constants
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 22
integer, parameter :: N_SCORE_TYPES = 23
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -310,9 +310,10 @@ module constants
SCORE_NU_SCATTER_YN = -17, & ! angular flux-weighted nu-scattering moment (0:N)
SCORE_EVENTS = -18, & ! number of events
SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate
SCORE_INVERSE_VELOCITY = -20, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -21, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -22 ! recoverable fission Q-value
SCORE_PROMPT_NU_FISSION = -20, & ! prompt neutron production rate
SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -23 ! recoverable fission Q-value
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10

2
src/endf.F90
View file

@ -42,6 +42,8 @@ contains
string = "nu-fission"
case (SCORE_DELAYED_NU_FISSION)
string = "delayed-nu-fission"
case (SCORE_PROMPT_NU_FISSION)
string = "prompt-nu-fission"
case (SCORE_KAPPA_FISSION)
string = "kappa-fission"
case (SCORE_CURRENT)

4
src/geometry.F90
View file

@ -499,7 +499,7 @@ contains
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
@ -563,7 +563,7 @@ contains
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then

11
src/global.F90
View file

@ -14,7 +14,7 @@ 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, TallyResult
use trigger_header, only: KTrigger
use timer_header, only: Timer
@ -137,6 +137,7 @@ module global
type(RegularMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
integer, allocatable :: matching_bins(:)
real(8), allocatable :: filter_weights(:)
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
@ -176,9 +177,6 @@ module global
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
!$omp& global_tally_tracklength, global_tally_leakage)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
integer :: n_meshes = 0 ! # of structured meshes
integer :: n_user_meshes = 0 ! # of structured user meshes
integer :: n_tallies = 0 ! # of tallies
@ -443,7 +441,8 @@ module global
type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
!$omp& trace, thread_id, current_work, matching_bins)
!$omp& trace, thread_id, current_work, matching_bins, &
!$omp& filter_weights)
contains
@ -509,7 +508,7 @@ contains
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) deallocate(tallies)
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(tally_maps)) deallocate(tally_maps)
if (allocated(filter_weights)) deallocate(filter_weights)
! Deallocate fission and source bank and entropy
!$omp parallel

111
src/initialize.F90
View file

@ -22,8 +22,9 @@ module initialize
use state_point, only: load_state_point
use string, only: to_str, starts_with, ends_with, str_to_int
use summary, only: write_summary
use tally_header, only: TallyObject, TallyResult, TallyFilter
use tally_header, only: TallyObject, TallyResult
use tally_initialize,only: configure_tallies
use tally_filter
use tally, only: init_tally_routines
#ifdef MPI
@ -711,76 +712,15 @@ contains
! =======================================================================
! ADJUST INDICES FOR EACH TALLY FILTER
FILTER_LOOP: do j = 1, t%n_filters
select case (t%filters(j)%type)
case (FILTER_DISTRIBCELL)
do k = 1, size(t%filters(j)%int_bins)
id = t%filters(j)%int_bins(k)
if (cell_dict%has_key(id)) then
t%filters(j)%int_bins(k) = cell_dict%get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_CELL, FILTER_CELLBORN)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (cell_dict%has_key(id)) then
t%filters(j)%int_bins(k) = cell_dict%get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_SURFACE)
FILTER_LOOP: do j = 1, size(t % filters)
select type(filt => t % filters(j) % obj)
type is (SurfaceFilter)
! Check if this is a surface filter only for surface currents
if (any(t%score_bins == SCORE_CURRENT)) cycle FILTER_LOOP
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (surface_dict%has_key(id)) then
t%filters(j)%int_bins(k) = surface_dict%get_key(id)
else
call fatal_error("Could not find surface " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_UNIVERSE)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (universe_dict%has_key(id)) then
t%filters(j)%int_bins(k) = universe_dict%get_key(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_MATERIAL)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (material_dict%has_key(id)) then
t%filters(j)%int_bins(k) = material_dict%get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_MESH)
! The mesh filter already has been set to the index in meshes rather
! than the user-specified id, so it doesn't need to be changed.
if (.not. any(t % score_bins == SCORE_CURRENT)) &
call filt % initialize()
class default
call filt % initialize()
end select
end do FILTER_LOOP
@ -894,14 +834,11 @@ contains
! We need distribcell if any tallies have distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
distribcell_active = .true.
if (size(tallies(i) % filters(j) % int_bins) > 1) then
call fatal_error("A distribcell filter was specified with &
&multiple bins. This feature is not supported.")
end if
end if
end select
end do
end do
@ -924,13 +861,12 @@ contains
! Set the number of bins in all distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
associate (filt => tallies(i) % filters(j))
if (filt % type == FILTER_DISTRIBCELL) then
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % int_bins(1)) % instances
end if
end associate
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % cell) % instances
end select
end do
end do
@ -990,10 +926,11 @@ contains
! List all cells referenced in distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
call cell_list % add(tallies(i) % filters(j) % int_bins(1))
end if
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
call cell_list % add(filt % cell)
end select
end do
end do

552
src/input_xml.F90
View file

@ -24,7 +24,8 @@ module input_xml
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with, tokenize, split_string
use tally_header, only: TallyObject, TallyFilter
use tally_header, only: TallyObject
use tally_filter
use tally_initialize, only: add_tallies
use xml_interface
@ -2667,7 +2668,7 @@ contains
type(ElemKeyValueCI), pointer :: pair_list
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(TallyFilter), allocatable :: filters(:) ! temporary filters
type(TallyFilterContainer), allocatable :: filters(:) ! temporary filters
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mesh => null()
type(Node), pointer :: node_tal => null()
@ -2907,115 +2908,122 @@ contains
call get_node_list(node_tal, "filter", node_filt_list)
n_filters = get_list_size(node_filt_list)
if (n_filters /= 0) then
! Allocate filters array
allocate(t % filters(n_filters))
! Allocate filters array
t % n_filters = n_filters
allocate(t % filters(n_filters))
READ_FILTERS: do j = 1, n_filters
! Get pointer to filter xml node
call get_list_item(node_filt_list, j, node_filt)
READ_FILTERS: do j = 1, n_filters
! Get pointer to filter xml node
call get_list_item(node_filt_list, j, node_filt)
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
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")
end if
! 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
call fatal_error("Bins not set in filter on tally " &
// trim(to_str(t % id)))
n_words = get_arraysize_integer(node_filt, "bins")
end if
else
call fatal_error("Bins not set in filter on tally " &
// trim(to_str(t % id)))
end if
! Determine type of filter
select case (temp_str)
! Determine type of filter
select case (temp_str)
case ('distribcell')
! Set type of filter
t % filters(j) % type = FILTER_DISTRIBCELL
! Going to add new filters to this tally if n_words > 1
case ('distribcell')
! Allocate and declare the filter type
allocate(DistribcellFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (DistribcellFilter)
if (n_words /= 1) call fatal_error("Only one cell can be &
&specified per distribcell filter.")
! Store bins
call get_node_value(node_filt, "bins", filt % cell)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_DISTRIBCELL) = j
case ('cell')
! Allocate and declare the filter type
allocate(CellFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (CellFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('cell')
! Set type of filter
t % filters(j) % type = FILTER_CELL
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_CELL) = j
case ('cellborn')
! Allocate and declare the filter type
allocate(CellbornFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (CellbornFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('cellborn')
! Set type of filter
t % filters(j) % type = FILTER_CELLBORN
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_CELLBORN) = j
case ('material')
! Allocate and declare the filter type
allocate(MaterialFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MaterialFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('material')
! Set type of filter
t % filters(j) % type = FILTER_MATERIAL
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % materials(n_words))
call get_node_array(node_filt, "bins", filt % materials)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MATERIAL) = j
case ('universe')
! Allocate and declare the filter type
allocate(UniverseFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (UniverseFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('universe')
! Set type of filter
t % filters(j) % type = FILTER_UNIVERSE
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % universes(n_words))
call get_node_array(node_filt, "bins", filt % universes)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_UNIVERSE) = j
case ('surface')
call fatal_error("Surface filter is not yet supported!")
! Allocate and declare the filter type
allocate(SurfaceFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (SurfaceFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
filt % n_bins = n_words
allocate(filt % surfaces(n_words))
call get_node_array(node_filt, "bins", filt % surfaces)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_SURFACE) = j
case ('surface')
call fatal_error("Surface filter is not yet supported!")
! Set type of filter
t % filters(j) % type = FILTER_SURFACE
! Set number of bins
t % filters(j) % n_bins = n_words
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('mesh')
! Set type of filter
t % filters(j) % type = FILTER_MESH
! Check to make sure multiple meshes weren't given
if (n_words /= 1) then
call fatal_error("Can only have one mesh filter specified.")
end if
case ('mesh')
! Allocate and declare the filter type
allocate(MeshFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MeshFilter)
if (n_words /= 1) call fatal_error("Only one mesh can be &
&specified per mesh filter.")
! Determine id of mesh
call get_node_value(node_filt, "bins", id)
@ -3032,214 +3040,220 @@ contains
! Determine number of bins -- this is assuming that the tally is
! a volume tally and not a surface current tally. If it is a
! surface current tally, the number of bins will get reset later
t % filters(j) % n_bins = product(m % dimension)
filt % n_bins = product(m % dimension)
! Allocate and store index of mesh
allocate(t % filters(j) % int_bins(1))
t % filters(j) % int_bins(1) = i_mesh
case ('energy')
! Set type of filter
t % filters(j) % type = FILTER_ENERGYIN
! Set number of bins
t % filters(j) % n_bins = n_words - 1
! Store the index of the mesh
filt % mesh = i_mesh
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MESH) = j
case ('energy')
! Allocate and declare the filter type
allocate(EnergyFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (EnergyFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! We can save tallying time if we know that the tally bins
! match the energy group structure. In that case, the matching bin
! We can save tallying time if we know that the tally bins match
! the energy group structure. In that case, the matching bin
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(t % filters(j) % real_bins == &
energy_bins(energy_groups + 1:1:-1))) &
t % energy_matches_groups = .true.
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_ENERGYIN) = j
case ('energyout')
! Set type of filter
t % filters(j) % type = FILTER_ENERGYOUT
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('energyout')
! Allocate and declare the filter type
allocate(EnergyoutFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (EnergyoutFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! We can save tallying time if we know that the tally bins
! match the energy group structure. In that case, the matching bin
! We can save tallying time if we know that the tally bins match
! the energy group structure. In that case, the matching bin
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(t % filters(j) % real_bins == &
energy_bins(energy_groups + 1:1:-1))) &
t % energyout_matches_groups = .true.
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_ENERGYOUT) = j
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('delayedgroup')
! Check to see if running in MG mode, because if so, the current
! system isnt set up yet to support delayed group data and thus
! these tallies
if (.not. run_CE) then
call fatal_error("delayedgroup filter on tally " &
// trim(to_str(t % id)) // " not yet supported&
& for multi-group mode.")
end if
! Set type of filter
t % filters(j) % type = FILTER_DELAYEDGROUP
! Set number of bins
t % filters(j) % n_bins = n_words
case ('delayedgroup')
! Check to see if running in MG mode, because if so, the current
! system isnt set up yet to support delayed group data and thus
! these tallies
if (.not. run_CE) then
call fatal_error("delayedgroup filter on tally " &
// trim(to_str(t % id)) // " not yet supported&
& for multi-group mode.")
end if
! Allocate and declare the filter type
allocate(DelayedGroupFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (DelayedGroupFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
filt % n_bins = n_words
allocate(filt % groups(n_words))
call get_node_array(node_filt, "bins", filt % groups)
! Check bins to make sure all are between 1 and MAX_DELAYED_GROUPS
! Check that bins are all are between 1 and MAX_DELAYED_GROUPS
do d = 1, n_words
if (t % filters(j) % int_bins(d) < 1 .or. &
t % filters(j) % int_bins(d) > MAX_DELAYED_GROUPS) then
if (filt % groups(d) < 1 .or. &
filt % groups(d) > MAX_DELAYED_GROUPS) then
call fatal_error("Encountered delayedgroup bin with index " &
// trim(to_str(t % filters(j) % int_bins(d))) // " that is&
& outside the range of 1 to MAX_DELAYED_GROUPS ( " &
// trim(to_str(filt % groups(d))) // " that is outside &
&the range of 1 to MAX_DELAYED_GROUPS ( " &
// trim(to_str(MAX_DELAYED_GROUPS)) // ")")
end if
end do
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_DELAYEDGROUP) = j
case ('mu')
! Set type of filter
t % filters(j) % type = FILTER_MU
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('mu')
! Allocate and declare the filter type
allocate(MuFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MuFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you subivide [-1,1] evenly with the input being the number of bins
! Allow a user to input a lone number which will mean that you
! subdivide [-1,1] evenly with the input being the number of bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = TWO / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = -ONE + (iangle - 1) * dangle
filt % bins(iangle) = -ONE + (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = ONE
filt % bins(Nangle + 1) = ONE
else
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MU) = j
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('polar')
! Set type of filter
t % filters(j) % type = FILTER_POLAR
! Set number of bins
t % filters(j) % n_bins = n_words - 1
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('polar')
! Allocate and declare the filter type
allocate(PolarFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (PolarFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you subivide [0,pi] evenly with the input being the number of bins
! Allow a user to input a lone number which will mean that you
! subdivide [0,pi] evenly with the input being the number of bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = PI / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = (iangle - 1) * dangle
filt % bins(iangle) = (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = PI
filt % bins(Nangle + 1) = PI
else
call fatal_error("Number of bins for polar filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_POLAR) = j
case ('azimuthal')
! Set type of filter
t % filters(j) % type = FILTER_AZIMUTHAL
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('azimuthal')
! Allocate and declare the filter type
allocate(AzimuthalFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (AzimuthalFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you sub-divide [-pi,pi) evenly with the input being the number of
! Allow a user to input a lone number which will mean that you
! subdivide [-pi,pi) evenly with the input being the number of
! bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = TWO * PI / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = -PI + (iangle - 1) * dangle
filt % bins(iangle) = -PI + (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = PI
filt % bins(Nangle + 1) = PI
else
call fatal_error("Number of bins for azimuthal filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
case default
! Specified tally filter is invalid, raise error
call fatal_error("Unknown filter type '" &
// trim(temp_str) // "' on tally " &
// trim(to_str(t % id)) // ".")
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_AZIMUTHAL) = j
! Set find_filter, e.g. if filter(3) has type FILTER_CELL, then
! find_filter(FILTER_CELL) would be set to 3.
case default
! Specified tally filter is invalid, raise error
call fatal_error("Unknown filter type '" &
// trim(temp_str) // "' on tally " &
// trim(to_str(t % id)) // ".")
t % find_filter(t % filters(j) % type) = j
end select
end do READ_FILTERS
end do READ_FILTERS
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
call fatal_error("Cannot specify both cell and surface filters for &
&tally " // trim(to_str(t % id)))
end if
else
! No filters were specified
t % n_filters = 0
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
call fatal_error("Cannot specify both cell and surface filters for &
&tally " // trim(to_str(t % id)))
end if
! =======================================================================
@ -3275,8 +3289,9 @@ contains
if (trim(sarray(j)) == 'total') then
! Check if a delayedgroup filter is present for this tally
do l = 1, t % n_filters
if (t % filters(l) % type == FILTER_DELAYEDGROUP) then
do l = 1, size(t % filters)
select type(filt => t % filters(l) % obj)
type is (DelayedGroupFilter)
call warning("A delayedgroup filter was used on a total &
&nuclide tally. Cross section libraries are not &
&guaranteed to have the same delayed group structure &
@ -3285,7 +3300,7 @@ contains
&all isotopes while the JEFF 3.1.1 library has the same &
&delayed group structure across all isotopes. Use with &
&caution!")
end if
end select
end do
t % nuclide_bins(j) = -1
@ -3334,8 +3349,9 @@ contains
t % n_nuclide_bins = 1
! Check if a delayedgroup filter is present for this tally
do l = 1, t % n_filters
if (t % filters(l) % type == FILTER_DELAYEDGROUP) then
do l = 1, size(t % filters)
select type(filt => t % filters(l) % obj)
type is (DelayedGroupFilter)
call warning("A delayedgroup filter was used on a total nuclide &
&tally. Cross section libraries are not guaranteed to have the&
& same delayed group structure across all isotopes. In &
@ -3343,7 +3359,7 @@ contains
&group structure across all isotopes while the JEFF 3.1.1 &
&library has the same delayed group structure across all &
&isotopes. Use with caution!")
end if
end select
end do
end if
@ -3619,6 +3635,12 @@ contains
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('prompt-nu-fission')
t % score_bins(j) = SCORE_PROMPT_NU_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
! Disallow for MG mode since data not present
if (.not. run_CE) then
@ -3657,36 +3679,42 @@ contains
&filter.")
end if
! Get pointer to mesh
i_mesh = t % filters(k) % int_bins(1)
m => meshes(i_mesh)
! Declare the type of the mesh filter
select type(filt => t % filters(k) % obj)
type is (MeshFilter)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
t % filters(k) % n_bins = product(m % dimension + 1)
! Get pointer to mesh
i_mesh = filt % mesh
m => meshes(i_mesh)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
filt % n_bins = product(m % dimension + 1)
end select
! Copy filters to temporary array
allocate(filters(t % n_filters + 1))
filters(1:t % n_filters) = t % filters
allocate(filters(size(t % filters) + 1))
filters(1:size(t % filters)) = t % filters
! Move allocation back -- filters becomes deallocated during
! this call
call move_alloc(FROM=filters, TO=t%filters)
! Add surface filter
t % n_filters = t % n_filters + 1
t % filters(t % n_filters) % type = FILTER_SURFACE
t % filters(t % n_filters) % n_bins = 2 * m % n_dimension
allocate(t % filters(t % n_filters) % int_bins(&
2 * m % n_dimension))
if (m % n_dimension == 2) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /)
end if
t % find_filter(FILTER_SURFACE) = t % n_filters
n_filters = size(t % filters)
allocate(SurfaceFilter :: t % filters(n_filters) % obj)
select type (filt => t % filters(size(t % filters)) % obj)
type is (SurfaceFilter)
filt % n_bins = 2 * m % n_dimension
allocate(filt % surfaces(2 * m % n_dimension))
if (m % n_dimension == 2) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT,&
IN_TOP, OUT_TOP /)
end if
end select
t % find_filter(FILTER_SURFACE) = size(t % filters)
case ('events')
t % score_bins(j) = SCORE_EVENTS
@ -4396,9 +4424,11 @@ contains
&meshlines on plot " // trim(to_str(pl % id)))
end if
i_mesh = cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
select type(filt => cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
pl % meshlines_mesh => meshes(i_mesh)
case ('entropy')

112
src/mesh.F90
View file

@ -290,34 +290,42 @@ contains
! Check if line intersects left surface -- calculate the intersection point
! y
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! x
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection
! point y
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! x
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_2d
@ -359,56 +367,68 @@ contains
! Check if line intersects left surface -- calculate the intersection point
! (y,z)
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects bottom surface -- calculate the intersection
! point (x,y)
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection point
! (y,z)
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects top surface -- calculate the intersection point
! (x,y)
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_3d

423
src/output.F90
View file

@ -17,6 +17,7 @@ module output
use sab_header, only: SAlphaBeta
use string, only: to_upper, to_str
use tally_header, only: TallyObject
use tally_filter
implicit none
@ -736,7 +737,6 @@ contains
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: l ! loop index for user scores
integer :: type ! type of tally filter
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
@ -747,7 +747,6 @@ contains
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
character(MAX_FILE_LEN) :: filename ! name of output file
character(16) :: filter_name(N_FILTER_TYPES) ! names of tally filters
character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function
character(36) :: score_name ! names of scoring function
! to be applied at write-time
@ -756,21 +755,6 @@ contains
! Skip if there are no tallies
if (n_tallies == 0) return
! Initialize names for tally filter types
filter_name(FILTER_UNIVERSE) = "Universe"
filter_name(FILTER_MATERIAL) = "Material"
filter_name(FILTER_DISTRIBCELL) = "Distributed Cell"
filter_name(FILTER_CELL) = "Cell"
filter_name(FILTER_CELLBORN) = "Birth Cell"
filter_name(FILTER_SURFACE) = "Surface"
filter_name(FILTER_MESH) = "Mesh"
filter_name(FILTER_ENERGYIN) = "Incoming Energy"
filter_name(FILTER_ENERGYOUT) = "Outgoing Energy"
filter_name(FILTER_MU) = "Change-in-Angle"
filter_name(FILTER_POLAR) = "Polar Angle"
filter_name(FILTER_AZIMUTHAL) = "Azimuthal Angle"
filter_name(FILTER_DELAYEDGROUP) = "Delayed Group"
! Initialize names for scores
score_names(abs(SCORE_FLUX)) = "Flux"
score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate"
@ -790,6 +774,7 @@ contains
score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
score_names(abs(SCORE_PROMPT_NU_FISSION)) = "Prompt-Nu-Fission Rate"
score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity"
score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power"
score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power"
@ -842,14 +827,14 @@ contains
! to be used for a given tally.
! Initialize bins, filter level, and indentation
matching_bins(1:t%n_filters) = 0
matching_bins(1:size(t % filters)) = 0
j = 1
indent = 0
print_bin: do
find_bin: do
! Check for no filters
if (t % n_filters == 0) exit find_bin
if (size(t % filters) == 0) exit find_bin
! Increment bin combination
matching_bins(j) = matching_bins(j) + 1
@ -857,7 +842,7 @@ contains
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (matching_bins(j) > t % filters(j) % n_bins) then
if (matching_bins(j) > t % filters(j) % obj % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
@ -870,12 +855,11 @@ contains
else
! Check if this is last filter
if (j == t % n_filters) exit find_bin
if (j == size(t % filters)) exit find_bin
! Print current filter information
type = t % filters(j) % type
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(t % filters(j) % obj % text_label(matching_bins(j)))
indent = indent + 2
j = j + 1
end if
@ -883,25 +867,25 @@ contains
end do find_bin
! Print filter information
if (t % n_filters > 0) then
type = t % filters(j) % type
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
if (size(t % filters) > 0) then
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(t % filters(j) % obj % text_label(matching_bins(j)))
end if
! Determine scoring index for this bin combination -- note that unlike
! in the score_tally subroutine, we have to use max(bins,1) since all
! bins below the lowest filter level will be zeros
if (t % n_filters > 0) then
filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1
if (size(t % filters) > 0) then
filter_index = sum((max(matching_bins(1:size(t % filters)),1) - 1) &
* t % stride) + 1
else
filter_index = 1
end if
! Write results for this filter bin combination
score_index = 0
if (t % n_filters > 0) indent = indent + 2
if (size(t % filters) > 0) indent = indent + 2
do n = 1, t % n_nuclide_bins
! Write label for nuclide
i_nuclide = t % nuclide_bins(n)
@ -973,7 +957,7 @@ contains
end do
indent = indent - 2
if (t % n_filters == 0) exit print_bin
if (size(t % filters) == 0) exit print_bin
end do print_bin
@ -1010,16 +994,19 @@ contains
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
select type(filt => t % filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! determine how many energy in bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = t % filters(i_filter_ein) % n_bins
n = t % filters(i_filter_ein) % obj % n_bins
else
print_ebin = .false.
n = 1
@ -1042,22 +1029,25 @@ contains
matching_bins(i_filter_ein) = l
! Write incoming energy bin
write(UNIT=unit_tally, FMT='(3X,A,1X,A)') &
"Incoming Energy", trim(get_label(t, i_filter_ein))
write(UNIT=unit_tally, FMT='(3X,A)') &
trim(t % filters(i_filter_ein) % obj % text_label( &
matching_bins(i_filter_ein)))
end if
! Left Surface
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % results(1,filter_index) % sum), &
@ -1067,14 +1057,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % results(1,filter_index) % sum), &
@ -1084,14 +1076,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % results(1,filter_index) % sum), &
@ -1101,14 +1095,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % results(1,filter_index) % sum), &
@ -1118,14 +1114,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % results(1,filter_index) % sum), &
@ -1135,14 +1133,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % results(1,filter_index) % sum), &
@ -1155,329 +1155,4 @@ contains
end subroutine write_surface_current
!===============================================================================
! GET_LABEL returns a label for a cell/surface/etc given a tally, filter type,
! and corresponding bin
!===============================================================================
function get_label(t, i_filter) result(label)
type(TallyObject), intent(in) :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(MAX_LINE_LEN) :: label ! user-specified identifier
integer :: i ! index in cells/surfaces/etc array
integer :: bin
integer :: offset
integer, allocatable :: ijk(:) ! indices in mesh
real(8) :: E0 ! lower bound for energy bin
real(8) :: E1 ! upper bound for energy bin
type(RegularMesh), pointer :: m
type(Universe), pointer :: univ
bin = matching_bins(i_filter)
select case(t % filters(i_filter) % type)
case (FILTER_UNIVERSE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(universes(i) % id)
case (FILTER_MATERIAL)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(materials(i) % id)
case (FILTER_CELL, FILTER_CELLBORN)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(cells(i) % id)
case (FILTER_DISTRIBCELL)
label = ''
univ => universes(BASE_UNIVERSE)
offset = 0
call find_offset(t % filters(i_filter) % int_bins(1), &
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(surfaces(i)%obj%id)
case (FILTER_MESH)
m => meshes(t % filters(i_filter) % int_bins(1))
allocate(ijk(m % n_dimension))
call bin_to_mesh_indices(m, bin, ijk)
if (m % n_dimension == 2) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
case (FILTER_ENERGYIN, FILTER_ENERGYOUT, FILTER_MU, FILTER_POLAR, &
FILTER_AZIMUTHAL)
E0 = t % filters(i_filter) % real_bins(bin)
E1 = t % filters(i_filter) % real_bins(bin + 1)
label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")"
case (FILTER_DELAYEDGROUP)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(i)
end select
end function get_label
!===============================================================================
! FIND_OFFSET 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 the target cell
! with the given offset
!===============================================================================
recursive subroutine find_offset(goal, univ, final, offset, path)
integer, intent(in) :: goal ! The target cell index
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: map ! Index in maps vector
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
integer :: n_x, n_y, n_z ! Lattice cell array dimensions
integer :: n ! Number of cells to search
integer :: cell_index ! Index in cells array
integer :: lat_offset ! Offset from lattice
integer :: temp_offset ! Looped sum of offsets
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
type(Cell), pointer :: c ! Pointer to current cell
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice
! Get the distribcell index for this cell
map = cells(goal) % distribcell_index
n = univ % n_cells
! Write to the geometry stack
if (univ%id == 0) then
path = trim(path) // to_str(univ%id)
else
path = trim(path) // "->" // to_str(univ%id)
end if
! Look through all cells in this universe
do i = 1, n
! If the cell matches the goal and the offset matches final, write to the
! geometry stack
if (univ % cells(i) == goal .and. offset == final) then
c => cells(univ % cells(i))
path = trim(path) // "->" // to_str(c % id)
return
end if
end do
! Find the fill cell or lattice cell that we need to enter
do i = 1, n
later_cell = .false.
cell_index = univ % cells(i)
c => cells(cell_index)
this_cell = .false.
! If we got here, we still think the target is in this universe
! or further down, but it's not this exact cell.
! Compare offset to next cell to see if we should enter this cell
if (i /= n) then
do j = i+1, n
cell_index = univ % cells(j)
c => cells(cell_index)
! Skip normal cells which do not have offsets
if (c % type == CELL_NORMAL) then
cycle
end if
! Break loop once we've found the next cell with an offset
exit
end do
! Ensure we didn't just end the loop by iteration
if (c % type /= CELL_NORMAL) then
! There are more cells in this universe that it could be in
later_cell = .true.
! Two cases, lattice or fill cell
if (c % type == CELL_FILL) then
temp_offset = c % offset(map)
! Get the offset of the first lattice location
else
lat => lattices(c % fill) % obj
temp_offset = lat % offset(map, 1, 1, 1)
end if
! If the final offset is in the range of offset - temp_offset+offset
! then the goal is in this cell
if (final < temp_offset + offset) then
this_cell = .true.
end if
end if
end if
if (n == 1 .and. c % type /= CELL_NORMAL) then
this_cell = .true.
end if
if (.not. later_cell) then
this_cell = .true.
end if
! Get pointer to THIS cell because target must be in this cell
if (this_cell) then
cell_index = univ % cells(i)
c => cells(cell_index)
path = trim(path) // "->" // to_str(c%id)
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == CELL_FILL) then
! Enter this cell to update the current offset
offset = c % offset(map) + offset
next_univ => universes(c % fill)
call find_offset(goal, next_univ, final, offset, path)
return
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
! ==================================================================
! RECTANGULAR LATTICES
type is (RectLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
n_z = lat % n_cells(3)
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do k = 1, n_x
do l = 1, n_y
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
path = trim(path) // "(" // trim(to_str(k)) // &
"," // trim(to_str(l)) // "," // &
trim(to_str(m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
old_l = l
old_k = k
cycle
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
return
end if
end do
end do
end do
! ==================================================================
! HEXAGONAL LATTICES
type is (HexLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
n_z = lat % n_axial
n_y = 2 * lat % n_rings - 1
n_x = 2 * lat % n_rings - 1
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do m = 1, n_z
do l = 1, n_y
do k = 1, n_x
! This array position is never used
if (k + l < lat % n_rings + 1) then
cycle
! This array position is never used
else if (k + l > 3*lat % n_rings - 1) then
cycle
end if
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == lat % n_rings .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
path = trim(path) // "(" // &
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
old_l = l
old_k = k
cycle
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(goal, next_univ, final, offset, path)
return
end if
end do
end do
end do
end select
end if
end if
end do
end subroutine find_offset
end module output

34
src/particle_header.F90
View file

@ -59,10 +59,13 @@ module particle_header
logical :: alive ! is particle alive?
! Pre-collision physical data
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_uvw(3) ! previous direction coordinates
real(8) :: last_wgt ! pre-collision particle weight
real(8) :: absorb_wgt ! weight absorbed for survival biasing
real(8) :: last_xyz_current(3) ! coordinates of the last collision or
! reflective/periodic surface crossing
! for current tallies
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_uvw(3) ! previous direction coordinates
real(8) :: last_wgt ! pre-collision particle weight
real(8) :: absorb_wgt ! weight absorbed for survival biasing
! What event last took place
logical :: fission ! did the particle cause implicit fission
@ -193,20 +196,21 @@ contains
call this % initialize()
! copy attributes from source bank site
this % wgt = src % wgt
this % last_wgt = src % wgt
this % coord(1) % xyz = src % xyz
this % coord(1) % uvw = src % uvw
this % last_xyz = src % xyz
this % last_uvw = src % uvw
this % wgt = src % wgt
this % last_wgt = src % wgt
this % coord(1) % xyz = src % xyz
this % coord(1) % uvw = src % uvw
this % last_xyz_current = src % xyz
this % last_xyz = src % xyz
this % last_uvw = src % uvw
if (run_CE) then
this % E = src % E
this % E = src % E
else
this % g = int(src % E)
this % last_g = int(src % E)
this % E = energy_bin_avg(this % g)
this % g = int(src % E)
this % last_g = int(src % E)
this % E = energy_bin_avg(this % g)
end if
this % last_E = this % E
this % last_E = this % E
end subroutine initialize_from_source

11
src/particle_restart.F90
View file

@ -107,11 +107,12 @@ contains
end if
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
p % last_wgt = p % wgt
p % last_xyz_current = p % coord(1)%xyz
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
! Close hdf5 file
call file_close(file_id)

50
src/state_point.F90
View file

@ -237,57 +237,13 @@ contains
end select
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
call write_dataset(tally_group, "n_filters", tally % n_filters)
call write_dataset(tally_group, "n_filters", size(tally % filters))
! Write filter information
FILTER_LOOP: do j = 1, tally % n_filters
FILTER_LOOP: do j = 1, size(tally % filters)
filter_group = create_group(tally_group, "filter " // &
trim(to_str(j)))
! Write name of type
select case (tally % filters(j) % type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type", "energyout")
case(FILTER_MU)
call write_dataset(filter_group, "type", "mu")
case(FILTER_POLAR)
call write_dataset(filter_group, "type", "polar")
case(FILTER_AZIMUTHAL)
call write_dataset(filter_group, "type", "azimuthal")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
case(FILTER_DELAYEDGROUP)
call write_dataset(filter_group, "type", "delayedgroup")
end select
call write_dataset(filter_group, "n_bins", &
tally % filters(j) % n_bins)
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT .or. &
tally % filters(j) % type == FILTER_MU .or. &
tally % filters(j) % type == FILTER_POLAR .or. &
tally % filters(j) % type == FILTER_AZIMUTHAL) then
call write_dataset(filter_group, "bins", &
tally % filters(j) % real_bins)
else
call write_dataset(filter_group, "bins", &
tally % filters(j) % int_bins)
end if
call tally % filters(j) % obj % to_statepoint(filter_group)
call close_group(filter_group)
end do FILTER_LOOP

76
src/summary.F90
View file

@ -3,7 +3,7 @@ module summary
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
&HexLattice, BASE_UNIVERSE
&HexLattice
use global
use hdf5_interface
use material_header, only: Material
@ -13,7 +13,6 @@ module summary
use surface_header
use string, only: to_str
use tally_header, only: TallyObject
use output, only: find_offset
use hdf5
@ -584,10 +583,6 @@ contains
type(RegularMesh), pointer :: m
type(TallyObject), pointer :: t
integer :: offset ! distibcell offset
character(MAX_LINE_LEN), allocatable :: paths(:) ! distribcell paths array
character(MAX_LINE_LEN) :: path ! distribcell path
tallies_group = create_group(file_id, "tallies")
! Write total number of meshes
@ -628,74 +623,11 @@ contains
call write_dataset(tally_group, "name", t%name)
! Write number of filters
call write_dataset(tally_group, "n_filters", t%n_filters)
call write_dataset(tally_group, "n_filters", size(t % filters))
FILTER_LOOP: do j = 1, t % n_filters
FILTER_LOOP: do j = 1, size(t % filters)
filter_group = create_group(tally_group, "filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_dataset(filter_group, "n_bins", t % filters(j) % n_bins)
! Write filter bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j)% type == FILTER_ENERGYOUT .or. &
t % filters(j) % type == FILTER_MU .or. &
t % filters(j) % type == FILTER_POLAR .or. &
t % filters(j) % type == FILTER_AZIMUTHAL) then
call write_dataset(filter_group, "bins", t % filters(j) % real_bins)
else
call write_dataset(filter_group, "bins", t % filters(j) % int_bins)
end if
! Write paths to reach each distribcell instance
if (t % filters(j) % type == FILTER_DISTRIBCELL) then
! Allocate array of strings for each distribcell path
allocate(paths(t % filters(j) % n_bins))
! Store path for each distribcell instance
do k = 1, t % filters(j) % n_bins
path = ''
offset = 1
call find_offset(t % filters(j) % int_bins(1), &
universes(BASE_UNIVERSE), k, offset, path)
paths(k) = path
end do
! Write array of distribcell paths to summary file
call write_dataset(filter_group, "paths", paths)
deallocate(paths)
end if
! Write name of type
select case (t%filters(j)%type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type", "energyout")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
case(FILTER_MU)
call write_dataset(filter_group, "type", "mu")
case(FILTER_POLAR)
call write_dataset(filter_group, "type", "polar")
case(FILTER_AZIMUTHAL)
call write_dataset(filter_group, "type", "azimuthal")
case(FILTER_DELAYEDGROUP)
call write_dataset(filter_group, "type", "delayedgroup")
end select
call t % filters(j) % obj % to_summary(filter_group)
call close_group(filter_group)
end do FILTER_LOOP

1942
src/tally.F90
View file

File diff suppressed because it is too large Load diff

1491
src/tally_filter.F90 Normal file
View file

File diff suppressed because it is too large Load diff

105
src/tally_filter_header.F90 Normal file
View file

@ -0,0 +1,105 @@
module tally_filter_header
use constants, only: MAX_LINE_LEN
use particle_header, only: Particle
use hdf5
implicit none
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell
! should score to the tally.
!===============================================================================
type, abstract :: TallyFilter
integer :: n_bins = 0
contains
procedure(get_next_bin_), deferred :: get_next_bin
procedure(to_statepoint_), deferred :: to_statepoint
procedure :: to_summary => filter_to_summary
procedure(text_label_), deferred :: text_label
procedure :: initialize => filter_initialize
end type TallyFilter
abstract interface
!===============================================================================
! GET_NEXT_BIN gives the index for the next valid filter bin and a weight that
! will be applied to the flux.
!
! In principle, a filter can have multiple valid bins. If current_bin =
! NO_BIN_FOUND, then this method should give the first valid bin. Providing the
! first valid bin should then give the second valid bin, and so on. When there
! are no valid bins left, the next_bin should be NO_VALID_BIN.
subroutine get_next_bin_(this, p, estimator, current_bin, next_bin, weight)
import TallyFilter
import Particle
class(TallyFilter), 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
end subroutine get_next_bin_
!===============================================================================
! TO_STATEPOINT writes all the information needed to reconstruct the filter to
! the given filter_group.
subroutine to_statepoint_(this, filter_group)
import TallyFilter
import HID_T
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
end subroutine to_statepoint_
!===============================================================================
! TEXT_LABEL returns a string describing the given filter bin. For example, an
! energy filter might return the string "Incoming Energy [0.625E-6, 20.0)".
! This is used to write the tallies.out file.
function text_label_(this, bin) result(label)
import TallyFilter
import MAX_LINE_LEN
class(TallyFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
end function text_label_
end interface
!===============================================================================
! TALLYFILTERCONTAINER contains an allocatable TallyFilter object for arrays of
! TallyFilters
!===============================================================================
type TallyFilterContainer
class(TallyFilter), allocatable :: obj
end type TallyFilterContainer
contains
!===============================================================================
! TO_SUMMARY writes all the information needed to reconstruct the filter to the
! given filter_group. If this procedure is not overridden by the derived class,
! then it will call to_statepoint by default.
subroutine filter_to_summary(this, filter_group)
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call this % to_statepoint(filter_group)
end subroutine filter_to_summary
!===============================================================================
! INITIALIZE sets up any internal data, as necessary. If this procedure is not
! overriden by the derived class, then it will do nothing by default.
subroutine filter_initialize(this)
class(TallyFilter), intent(inout) :: this
end subroutine filter_initialize
end module tally_filter_header

59
src/tally_header.F90
View file

@ -1,41 +1,13 @@
module tally_header
use constants, only: NONE, N_FILTER_TYPES
use trigger_header, only: TriggerObject
use constants, only: NONE, N_FILTER_TYPES
use tally_filter_header, only: TallyFilterContainer
use trigger_header, only: TriggerObject
use, intrinsic :: ISO_C_BINDING
implicit none
!===============================================================================
! TALLYMAPELEMENT gives an index to a tally which is to be scored and the
! corresponding bin for the filter variable
!===============================================================================
type TallyMapElement
integer :: index_tally
integer :: index_bin
end type TallyMapElement
!===============================================================================
! TALLYMAPITEM contains a list of tally/bin combinations for each mappable
! filter bin specified.
!===============================================================================
type TallyMapItem
type(TallyMapElement), allocatable :: elements(:)
end type TallyMapItem
!===============================================================================
! TALLYMAP contains a list of pairs of indices to tallies and the corresponding
! bin for a given filter. There is one TallyMap for each mappable filter
! type. The items array is as long as the corresponding array for that filter,
! e.g. for tally_maps(FILTER_CELL), items is n_cells long.
!===============================================================================
type TallyMap
type(TallyMapItem), allocatable :: items(:)
end type TallyMap
!===============================================================================
! TALLYRESULT provides accumulation of results in a particular tally bin
!===============================================================================
@ -46,19 +18,6 @@ module tally_header
real(C_DOUBLE) :: sum_sq = 0.
end type TallyResult
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell
! should score to the tally.
!===============================================================================
type TallyFilter
integer :: type = NONE
integer :: n_bins = 0
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
end type TallyFilter
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyFilters and the results are stored in a
@ -73,11 +32,7 @@ module tally_header
integer :: type ! volume, surface current
integer :: estimator ! collision, track-length
real(8) :: volume ! volume of region
! Information about what filters should be used
integer :: n_filters ! Number of filters
type(TallyFilter), allocatable :: filters(:) ! Filter data (type/bins)
type(TallyFilterContainer), allocatable :: filters(:)
! The stride attribute is used for determining the index in the results
! array for a matching_bin combination. Since multiple dimensions are
@ -124,10 +79,6 @@ module tally_header
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Multi-Group Specific Information To Enable Rapid Tallying
logical :: energy_matches_groups = .false.
logical :: energyout_matches_groups = .false.
end type TallyObject
end module tally_header

103
src/tally_initialize.F90
View file

@ -2,7 +2,7 @@ module tally_initialize
use constants
use global
use tally_header, only: TallyObject, TallyMapElement, TallyMapItem
use tally_header, only: TallyObject
implicit none
private
@ -23,7 +23,6 @@ contains
allocate(global_tallies(N_GLOBAL_TALLIES))
call setup_tally_arrays()
call setup_tally_maps()
end subroutine configure_tallies
@ -45,17 +44,17 @@ contains
t => tallies(i)
! Allocate stride and matching_bins arrays
allocate(t % stride(t % n_filters))
max_n_filters = max(max_n_filters, t % n_filters)
allocate(t % stride(size(t % filters)))
max_n_filters = max(max_n_filters, size(t % filters))
! The filters are traversed in opposite order so that the last filter has
! the shortest stride in memory and the first filter has the largest
! stride
n = 1
STRIDE: do j = t % n_filters, 1, -1
STRIDE: do j = size(t % filters), 1, -1
t % stride(j) = n
n = n * t % filters(j) % n_bins
n = n * t % filters(j) % obj % n_bins
end do STRIDE
! Set total number of filter and scoring bins
@ -70,101 +69,11 @@ contains
! Allocate array for matching filter bins
!$omp parallel
allocate(matching_bins(max_n_filters))
allocate(filter_weights(max_n_filters))
!$omp end parallel
end subroutine setup_tally_arrays
!===============================================================================
! SETUP_TALLY_MAPS creates a map that allows a quick determination of which
! tallies and bins need to be scored to when a particle makes a collision. This
! subroutine also sets the stride attribute for each tally as well as allocating
! storage for the results array.
!===============================================================================
subroutine setup_tally_maps()
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: k ! loop index for bins
integer :: bin ! filter bin entries
integer :: type ! type of tally filter
type(TallyObject), pointer :: t
! allocate tally map array -- note that we don't need a tally map for the
! energy_in and energy_out filters
allocate(tally_maps(N_FILTER_TYPES - 3))
! allocate list of items for each different filter type
allocate(tally_maps(FILTER_UNIVERSE) % items(n_universes))
allocate(tally_maps(FILTER_MATERIAL) % items(n_materials))
allocate(tally_maps(FILTER_CELL) % items(n_cells))
allocate(tally_maps(FILTER_CELLBORN) % items(n_cells))
allocate(tally_maps(FILTER_SURFACE) % items(n_surfaces))
TALLY_LOOP: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! No need to set up tally maps for surface current tallies
if (t % type == TALLY_SURFACE_CURRENT) cycle
FILTER_LOOP: do j = 1, t % n_filters
! Determine type of filter
type = t % filters(j) % type
if (type == FILTER_CELL .or. type == FILTER_SURFACE .or. &
type == FILTER_MATERIAL .or. type == FILTER_UNIVERSE .or. &
type == FILTER_CELLBORN) then
! Add map elements
BIN_LOOP: do k = 1, t % filters(j) % n_bins
bin = t % filters(j) % int_bins(k)
call add_map_element(tally_maps(type) % items(bin), i, k)
end do BIN_LOOP
end if
end do FILTER_LOOP
end do TALLY_LOOP
end subroutine setup_tally_maps
!===============================================================================
! ADD_MAP_ELEMENT adds a pair of tally and bin indices to the list for a given
! cell/surface/etc.
!===============================================================================
subroutine add_map_element(item, index_tally, index_bin)
type(TallyMapItem), intent(inout) :: item
integer, intent(in) :: index_tally ! index in tallies array
integer, intent(in) :: index_bin ! index in bins array
integer :: n ! size of elements array
type(TallyMapElement), allocatable :: temp(:)
if (.not. allocated(item % elements)) then
allocate(item % elements(1))
item % elements(1) % index_tally = index_tally
item % elements(1) % index_bin = index_bin
else
! determine size of elements array
n = size(item % elements)
! allocate temporary storage and copy elements
allocate(temp(n+1))
temp(1:n) = item % elements
! move allocation back to main array
call move_alloc(FROM=temp, TO=item%elements)
! set new element
item % elements(n+1) % index_tally = index_tally
item % elements(n+1) % index_bin = index_bin
end if
end subroutine add_map_element
!===============================================================================
! ADD_TALLIES extends the tallies array with a new group of tallies and assigns
! pointers to each group. This is called once for user tallies, once for CMFD

5
src/tracking.F90
View file

@ -191,7 +191,7 @@ contains
p % fission = .false.
! Save coordinates for tallying purposes
p % last_xyz = p % coord(1) % xyz
p % last_xyz_current = p % coord(1) % xyz
! Set last material to none since cross sections will need to be
! re-evaluated
@ -211,6 +211,9 @@ contains
end do
end if
! Save coordinates for tallying purposes
p % last_xyz = p % coord(1) % xyz
! If particle has too many events, display warning and kill it
n_event = n_event + 1
if (n_event == MAX_EVENTS) then

38
src/trigger.F90
View file

@ -12,6 +12,7 @@ module trigger
use mesh_header, only: RegularMesh
use trigger_header, only: TriggerObject
use tally, only: TallyObject
use tally_filter, only: MeshFilter
implicit none
@ -165,7 +166,7 @@ contains
else
! Initialize bins, filter level
matching_bins(1:t % n_filters) = 0
matching_bins(1:size(t % filters)) = 0
FILTER_LOOP: do filter_index = 1, t % total_filter_bins
@ -265,7 +266,7 @@ contains
end if
end if
end do NUCLIDE_LOOP
if (t % n_filters == 0) exit FILTER_LOOP
if (size(t % filters) == 0) exit FILTER_LOOP
end do FILTER_LOOP
end if
end do TRIGGER_LOOP
@ -300,16 +301,19 @@ contains
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
select type(filt => t % filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:t % n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! determine how many energyin bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = t % filters(i_filter_ein) % n_bins
n = t % filters(i_filter_ein) % obj % n_bins
else
print_ebin = .false.
n = 1
@ -329,7 +333,7 @@ contains
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -341,7 +345,7 @@ contains
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -356,7 +360,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -368,7 +372,7 @@ contains
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -383,7 +387,7 @@ contains
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -396,7 +400,7 @@ contains
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -411,7 +415,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -423,7 +427,7 @@ contains
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -438,7 +442,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -450,7 +454,7 @@ contains
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -465,7 +469,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -477,7 +481,7 @@ contains
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev

2
tests/test_mgxs_library_condense/inputs_true.dat
View file

@ -1 +1 @@
855919f7a333acff6423527b82656d6a472ea8416002fb475c2d21c676e2f75f5c040d209a3e9a3a6118cf944f2c1bfb8cf2318273b890861d790b1927791a35
e2cdca7ea5b3532050af5b12fac26d7ef212d2696bb1b73cdd00929b2243c40d100ad02438c7b090555b49815d0de6c48cf1b4ebf437a48bc80c2d2b4bad292e

18
tests/test_mgxs_library_condense/results_true.dat
View file

@ -34,6 +34,12 @@
0 10000 1 1 total 0.085835 0.005592
material group out nuclide mean std. dev.
0 10000 1 total 1.0 0.046071
material group out nuclide mean std. dev.
0 10000 1 total 1.0 0.051471
material group in nuclide mean std. dev.
0 10000 1 total 4.996730e-07 3.650635e-08
material group in nuclide mean std. dev.
0 10000 1 total 0.090004 0.006367
material group in nuclide mean std. dev.
0 10001 1 total 0.311594 0.013793
material group in nuclide mean std. dev.
@ -70,6 +76,12 @@
0 10001 1 1 total 0.0 0.0
material group out nuclide mean std. dev.
0 10001 1 total 0.0 0.0
material group out nuclide mean std. dev.
0 10001 1 total 0.0 0.0
material group in nuclide mean std. dev.
0 10001 1 total 5.454760e-07 4.949800e-08
material group in nuclide mean std. dev.
0 10001 1 total 0.0 0.0
material group in nuclide mean std. dev.
0 10002 1 total 0.904999 0.043964
material group in nuclide mean std. dev.
@ -106,3 +118,9 @@
0 10002 1 1 total 0.0 0.0
material group out nuclide mean std. dev.
0 10002 1 total 0.0 0.0
material group out nuclide mean std. dev.
0 10002 1 total 0.0 0.0
material group in nuclide mean std. dev.
0 10002 1 total 5.773006e-07 5.322132e-08
material group in nuclide mean std. dev.
0 10002 1 total 0.0 0.0

2
tests/test_mgxs_library_distribcell/inputs_true.dat
View file

@ -1 +1 @@
3a3b7f75b326c94a8e5c7efe3046b2fdb887e9f75ecf6eb27587f9450c77cf8fd6acc4198c15bffb4e7ceead6d7b4327c19536bbf9cc35dfaae3f4ce4c26cc1a
2d948f3b12293294eaeca231a3df9d51195379e8bb38dd3e68d3bc512a7d08ed52a1109054ca381684ec127268710f6d6e9210ac8154c9b379608e996627624a

6
tests/test_mgxs_library_distribcell/results_true.dat
View file

@ -34,3 +34,9 @@
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.0 0.0
avg(distribcell) group out nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.0 0.0
avg(distribcell) group out nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.0 0.0
avg(distribcell) group in nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.000001 6.946255e-07
avg(distribcell) group in nuclide mean std. dev.
0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.0 0.0

2
tests/test_mgxs_library_hdf5/inputs_true.dat
View file

@ -1 +1 @@
855919f7a333acff6423527b82656d6a472ea8416002fb475c2d21c676e2f75f5c040d209a3e9a3a6118cf944f2c1bfb8cf2318273b890861d790b1927791a35
e2cdca7ea5b3532050af5b12fac26d7ef212d2696bb1b73cdd00929b2243c40d100ad02438c7b090555b49815d0de6c48cf1b4ebf437a48bc80c2d2b4bad292e

27
tests/test_mgxs_library_hdf5/results_true.dat
View file

@ -63,6 +63,15 @@ domain=10000 type=nu-fission matrix
domain=10000 type=chi
[ 1. 0.]
[ 0.04607052 0. ]
domain=10000 type=chi-prompt
[ 1. 0.]
[ 0.05147146 0. ]
domain=10000 type=inverse-velocity
[ 5.70932329e-08 2.85573950e-06]
[ 4.68792969e-09 2.44216503e-07]
domain=10000 type=prompt-nu-fission
[ 0.01923922 0.46671903]
[ 0.00130951 0.04141087]
domain=10001 type=total
[ 0.31373767 0.3008214 ]
[ 0.0155819 0.02805245]
@ -128,6 +137,15 @@ domain=10001 type=nu-fission matrix
domain=10001 type=chi
[ 0. 0.]
[ 0. 0.]
domain=10001 type=chi-prompt
[ 0. 0.]
[ 0. 0.]
domain=10001 type=inverse-velocity
[ 5.99598048e-08 2.98549021e-06]
[ 4.55309296e-09 3.41701554e-07]
domain=10001 type=prompt-nu-fission
[ 0. 0.]
[ 0. 0.]
domain=10002 type=total
[ 0.66457226 2.05238401]
[ 0.03121475 0.22434291]
@ -193,3 +211,12 @@ domain=10002 type=nu-fission matrix
domain=10002 type=chi
[ 0. 0.]
[ 0. 0.]
domain=10002 type=chi-prompt
[ 0. 0.]
[ 0. 0.]
domain=10002 type=inverse-velocity
[ 6.02207831e-08 3.04495537e-06]
[ 3.78043696e-09 3.60007673e-07]
domain=10002 type=prompt-nu-fission
[ 0. 0.]
[ 0. 0.]

1
tests/test_mgxs_library_mesh/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
a4cd030bea212e45fdb159e75a7fb3d1947e9bf3d0384ac5d37a72298d67dcfdd1b9eb5c6af8ac6e5983bd5b47de9c17a2ea472b467b7222a4909ee070bf1ca3

132
tests/test_mgxs_library_mesh/results_true.dat Normal file
View file

@ -0,0 +1,132 @@
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.640786 0.044177
1 1 2 1 1 total 0.660597 0.128423
2 2 1 1 1 total 0.615276 0.104046
3 2 2 1 1 total 0.646999 0.186709
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.36665 0.048814
1 1 2 1 1 total 0.40784 0.096486
2 2 1 1 1 total 0.36356 0.074111
3 2 2 1 1 total 0.41456 0.160443
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.366650 0.048814
1 1 2 1 1 total 0.407840 0.096486
2 2 1 1 1 total 0.363560 0.074111
3 2 2 1 1 total 0.414593 0.160436
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.025749 0.002863
1 1 2 1 1 total 0.028400 0.005275
2 2 1 1 1 total 0.022988 0.004099
3 2 2 1 1 total 0.027589 0.010350
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.015861 0.002876
1 1 2 1 1 total 0.017280 0.004371
2 2 1 1 1 total 0.014403 0.003542
3 2 2 1 1 total 0.018061 0.010110
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.009888 0.001077
1 1 2 1 1 total 0.011121 0.002456
2 2 1 1 1 total 0.008585 0.001552
3 2 2 1 1 total 0.009527 0.003659
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.026065 0.002907
1 1 2 1 1 total 0.029084 0.006430
2 2 1 1 1 total 0.022596 0.004062
3 2 2 1 1 total 0.025066 0.009687
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 1.938476 0.211550
1 1 2 1 1 total 2.177360 0.480780
2 2 1 1 1 total 1.682799 0.303764
3 2 2 1 1 total 1.864890 0.715661
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.615037 0.041754
1 1 2 1 1 total 0.632196 0.123878
2 2 1 1 1 total 0.592288 0.100439
3 2 2 1 1 total 0.619410 0.177190
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.584014 0.054315
1 1 2 1 1 total 0.622514 0.111323
2 2 1 1 1 total 0.587256 0.084833
3 2 2 1 1 total 0.613792 0.168612
mesh 1 group in group out nuclide moment mean std. dev.
x y z
0 1 1 1 1 1 total P0 0.584014 0.054315
1 1 1 1 1 1 total P1 0.243427 0.025488
2 1 1 1 1 1 total P2 0.089236 0.007357
3 1 1 1 1 1 total P3 0.008994 0.005768
4 1 2 1 1 1 total P0 0.622514 0.111323
5 1 2 1 1 1 total P1 0.239376 0.042594
6 1 2 1 1 1 total P2 0.088386 0.017200
7 1 2 1 1 1 total P3 -0.001243 0.005639
8 2 1 1 1 1 total P0 0.587256 0.084833
9 2 1 1 1 1 total P1 0.245120 0.041033
10 2 1 1 1 1 total P2 0.086784 0.016255
11 2 1 1 1 1 total P3 0.008660 0.004755
12 2 2 1 1 1 total P0 0.612950 0.167940
13 2 2 1 1 1 total P1 0.226176 0.061882
14 2 2 1 1 1 total P2 0.086593 0.026126
15 2 2 1 1 1 total P3 0.009672 0.011995
mesh 1 group in group out nuclide moment mean std. dev.
x y z
0 1 1 1 1 1 total P0 0.584014 0.054315
1 1 1 1 1 1 total P1 0.243427 0.025488
2 1 1 1 1 1 total P2 0.089236 0.007357
3 1 1 1 1 1 total P3 0.008994 0.005768
4 1 2 1 1 1 total P0 0.622514 0.111323
5 1 2 1 1 1 total P1 0.239376 0.042594
6 1 2 1 1 1 total P2 0.088386 0.017200
7 1 2 1 1 1 total P3 -0.001243 0.005639
8 2 1 1 1 1 total P0 0.587256 0.084833
9 2 1 1 1 1 total P1 0.245120 0.041033
10 2 1 1 1 1 total P2 0.086784 0.016255
11 2 1 1 1 1 total P3 0.008660 0.004755
12 2 2 1 1 1 total P0 0.613792 0.168612
13 2 2 1 1 1 total P1 0.226142 0.061856
14 2 2 1 1 1 total P2 0.086174 0.025979
15 2 2 1 1 1 total P3 0.009721 0.012027
mesh 1 group in group out nuclide mean std. dev.
x y z
0 1 1 1 1 1 total 1.000000 0.088094
1 1 2 1 1 1 total 1.000000 0.160891
2 2 1 1 1 1 total 1.000000 0.126864
3 2 2 1 1 1 total 1.001374 0.305883
mesh 1 group in group out nuclide mean std. dev.
x y z
0 1 1 1 1 1 total 0.027395 0.004680
1 1 2 1 1 1 total 0.022914 0.006025
2 2 1 1 1 1 total 0.019384 0.002846
3 2 2 1 1 1 total 0.029629 0.006292
mesh 1 group out nuclide mean std. dev.
x y z
0 1 1 1 1 total 1.0 0.220956
1 1 2 1 1 total 1.0 0.316565
2 2 1 1 1 total 1.0 0.132140
3 2 2 1 1 total 1.0 0.181577
mesh 1 group out nuclide mean std. dev.
x y z
0 1 1 1 1 total 1.0 0.222246
1 1 2 1 1 total 1.0 0.316565
2 2 1 1 1 total 1.0 0.132140
3 2 2 1 1 total 1.0 0.181577
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 3.610522e-07 3.169931e-08
1 1 2 1 1 total 3.942353e-07 8.459167e-08
2 2 1 1 1 total 3.097784e-07 5.252025e-08
3 2 2 1 1 total 3.799163e-07 1.806470e-07
mesh 1 group in nuclide mean std. dev.
x y z
0 1 1 1 1 total 0.025920 0.002893
1 1 2 1 1 total 0.028922 0.006394
2 2 1 1 1 total 0.022467 0.004039
3 2 2 1 1 total 0.024923 0.009632

81
tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py Normal file
View file

@ -0,0 +1,81 @@
#!/usr/bin/env python
import os
import sys
import glob
import hashlib
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import openmc
import openmc.mgxs
class MGXSTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Generate inputs using parent class routine
super(MGXSTestHarness, self)._build_inputs()
# Initialize a one-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 20.])
# Initialize MGXS Library for a few cross section types
# for one material-filled cell in the geometry
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)
self.mgxs_lib.by_nuclide = False
# Test all MGXS types
self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES
self.mgxs_lib.energy_groups = energy_groups
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'mesh'
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [2, 2]
mesh.lower_left = [-100., -100.]
mesh.width = [100., 100.]
self.mgxs_lib.domains = [mesh]
self.mgxs_lib.build_library()
# Initialize a tallies file
self._input_set.tallies = openmc.Tallies()
self.mgxs_lib.add_to_tallies_file(self._input_set.tallies, merge=False)
self._input_set.tallies.export_to_xml()
def _get_results(self, hash_output=False):
"""Digest info in the statepoint and return as a string."""
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = openmc.StatePoint(statepoint)
# Load the MGXS library from the statepoint
self.mgxs_lib.load_from_statepoint(sp)
# Build a string from Pandas Dataframe for each 1-group MGXS
outstr = ''
for domain in self.mgxs_lib.domains:
for mgxs_type in self.mgxs_lib.mgxs_types:
mgxs = self.mgxs_lib.get_mgxs(domain, mgxs_type)
df = mgxs.get_pandas_dataframe()
outstr += df.to_string() + '\n'
# Hash the results if necessary
if hash_output:
sha512 = hashlib.sha512()
sha512.update(outstr.encode('utf-8'))
outstr = sha512.hexdigest()
return outstr
def _cleanup(self):
super(MGXSTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = MGXSTestHarness('statepoint.10.*', True)
harness.main()

2
tests/test_mgxs_library_no_nuclides/inputs_true.dat
View file

@ -1 +1 @@
855919f7a333acff6423527b82656d6a472ea8416002fb475c2d21c676e2f75f5c040d209a3e9a3a6118cf944f2c1bfb8cf2318273b890861d790b1927791a35
e2cdca7ea5b3532050af5b12fac26d7ef212d2696bb1b73cdd00929b2243c40d100ad02438c7b090555b49815d0de6c48cf1b4ebf437a48bc80c2d2b4bad292e

27
tests/test_mgxs_library_no_nuclides/results_true.dat
View file

@ -75,6 +75,15 @@
material group out nuclide mean std. dev.
1 10000 1 total 1.0 0.046071
0 10000 2 total 0.0 0.000000
material group out nuclide mean std. dev.
1 10000 1 total 1.0 0.051471
0 10000 2 total 0.0 0.000000
material group in nuclide mean std. dev.
1 10000 1 total 5.709323e-08 4.687930e-09
0 10000 2 total 2.855740e-06 2.442165e-07
material group in nuclide mean std. dev.
1 10000 1 total 0.019239 0.001310
0 10000 2 total 0.466719 0.041411
material group in nuclide mean std. dev.
1 10001 1 total 0.313738 0.015582
0 10001 2 total 0.300821 0.028052
@ -152,6 +161,15 @@
material group out nuclide mean std. dev.
1 10001 1 total 0.0 0.0
0 10001 2 total 0.0 0.0
material group out nuclide mean std. dev.
1 10001 1 total 0.0 0.0
0 10001 2 total 0.0 0.0
material group in nuclide mean std. dev.
1 10001 1 total 5.995980e-08 4.553093e-09
0 10001 2 total 2.985490e-06 3.417016e-07
material group in nuclide mean std. dev.
1 10001 1 total 0.0 0.0
0 10001 2 total 0.0 0.0
material group in nuclide mean std. dev.
1 10002 1 total 0.664572 0.031215
0 10002 2 total 2.052384 0.224343
@ -229,3 +247,12 @@
material group out nuclide mean std. dev.
1 10002 1 total 0.0 0.0
0 10002 2 total 0.0 0.0
material group out nuclide mean std. dev.
1 10002 1 total 0.0 0.0
0 10002 2 total 0.0 0.0
material group in nuclide mean std. dev.
1 10002 1 total 6.022078e-08 3.780437e-09
0 10002 2 total 3.044955e-06 3.600077e-07
material group in nuclide mean std. dev.
1 10002 1 total 0.0 0.0
0 10002 2 total 0.0 0.0

2
tests/test_mgxs_library_nuclides/inputs_true.dat
View file

@ -1 +1 @@
739796983940a1bad601998cf9ea2f90453a994477c7f675c2fd404d1864fe04a7fbfb5a15c5fe7cf9bad016b78432ba0910baba6f9cc026143761e9f526b62a
e4a5f03ab6167e96462c4ef537533fe33b98d7878ae00824c5619356bda8d548b3c71af01ba8c88d5a9b46dd1471d331e6f678a164af922200f2ee3642be6340

2
tests/test_mgxs_library_nuclides/results_true.dat
View file

@ -1 +1 @@
d56c6bae6bf3cd8950d3f50f089458c1c6c807be780fc97570532c6af6eb7e3057968d3345bd3c363f01315271129ef7b2ca028b05767353e601dc37b035f8a7
e494320a213b5704a2ac915a2ba504857be91961ceb6735b6ad05d81eb31c44c9584d5bd9d40baececf1dcb5b030e6ecec63cfbd20639baf69bcb596c5c46591

2
tests/test_score_current/results_true.dat
View file

@ -1 +1 @@
a9310752363eb059ff40f16ac9716b41ccab6ec6607d29f498069318745e485d18d784264304cc2586865bd58cef7587203cc22a1d485c58ddd63c14c0defdb9
bafab1921a12146abb2bb29603b52b9cc28a5a950a7a6bb1e3f012c05891c310fad643760d4f148b04d0fef3d1f3e141d146e3a278d81cc6fc8187c37717c5e7

2
tests/test_tallies/inputs_true.dat
View file

@ -1 +1 @@
a392a7a8f27fd2f959b06a6809df29b3482e215da175b19846c248af44dc2ee7e2b05269b802dd9d6ed434506b05092f349436a0048411adab6b296dba0bd683
930af242a043f2676a000dbc5a2db6b148edcb31ed8c87dbaa35a8efb37a3be8cff30cdf4dc03f9c5c7eb4021f7e4c3327e64681cdd8fd8722c95c69db850227

2
tests/test_tallies/results_true.dat
View file

@ -1 +1 @@
a0c7d6ca246ecd7dd5fed06373af142390971401c4e97744f29e55810ab9c231c97c4d8947cdf0b3d2df0ae829a9ddf768e5b2d889bbea34f2b6db0e567db884
a51db2a4efc681805f85968e04411dc33beee0532c202f5179b9a82880ab60a75e53fa9141c81045ea1d2842372f2d8da900326f09382ea61dd80a3c9b43bba1

10
tests/test_tallies/test_tallies.py
View file

@ -113,21 +113,23 @@ class TalliesTestHarness(PyAPITestHarness):
polar_tally4.estimator = 'tracklength'
universe_tally = Tally()
universe_tally.filters = [Filter(type='universe', bins=(1, 2, 3, 4))]
universe_tally.filters = [
Filter(type='universe', bins=(1, 2, 3, 4, 6, 8))]
universe_tally.scores = ['total']
cell_filter = Filter(type='cell', bins=(10, 21, 22, 23))
cell_filter = Filter(type='cell', bins=(10, 21, 22, 23, 60))
score_tallies = [Tally(), Tally(), Tally()]
for t in score_tallies:
t.filters = [cell_filter]
t.scores = ['absorption', 'delayed-nu-fission', 'events', 'fission',
'inverse-velocity', 'kappa-fission', '(n,2n)', '(n,n1)',
'(n,gamma)', 'nu-fission', 'scatter', 'elastic', 'total']
'(n,gamma)', 'nu-fission', 'scatter', 'elastic', 'total',
'prompt-nu-fission']
score_tallies[0].estimator = 'tracklength'
score_tallies[1].estimator = 'analog'
score_tallies[2].estimator = 'collision'
cell_filter2 = Filter(type='cell', bins=(21, 22, 23, 27, 28, 29))
cell_filter2 = Filter(type='cell', bins=(21, 22, 23, 27, 28, 29, 60))
flux_tallies = [Tally() for i in range(4)]
for t in flux_tallies:
t.filters = [cell_filter2]