mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Refactor volume calculation to allow material/universe volumes as well
This commit is contained in:
Paul Romano
parent
99adb712ab
commit
a5763ae599
10 changed files with 298 additions and 140 deletions
|
|
@ -369,10 +369,13 @@ class Cell(object):
|
|||
Results from a stochastic volume calculation
|
||||
|
||||
"""
|
||||
for cell_id in volume_calc.results:
|
||||
if cell_id == self.id:
|
||||
self._volume_information = volume_calc.results[cell_id]
|
||||
break
|
||||
if volume_calc.domain_type == 'cell':
|
||||
for cell_id in volume_calc.results:
|
||||
if cell_id == self.id:
|
||||
self._volume_information = volume_calc.results[cell_id]
|
||||
break
|
||||
else:
|
||||
raise ValueError('No volume information found for this cell.')
|
||||
else:
|
||||
raise ValueError('No volume information found for this cell.')
|
||||
|
||||
|
|
|
|||
|
|
@ -59,9 +59,10 @@ class Geometry(object):
|
|||
Results from a stochastic volume calculation
|
||||
|
||||
"""
|
||||
for cell in self.get_all_cells():
|
||||
if cell.id in volume_calc.results:
|
||||
cell.add_volume_information(volume_calc)
|
||||
if volume_calc.domain_type == 'cell':
|
||||
for cell in self.get_all_cells():
|
||||
if cell.id in volume_calc.results:
|
||||
cell.add_volume_information(volume_calc)
|
||||
|
||||
def export_to_xml(self):
|
||||
"""Create a geometry.xml file that can be used for a simulation.
|
||||
|
|
|
|||
101
openmc/volume.py
101
openmc/volume.py
|
|
@ -5,7 +5,7 @@ from xml.etree import ElementTree as ET
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from openmc import Cell, Union
|
||||
import openmc
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
|
||||
|
|
@ -14,8 +14,8 @@ class VolumeCalculation(object):
|
|||
|
||||
Parameters
|
||||
----------
|
||||
cells : Iterable of Cell
|
||||
Cells to find volumes of
|
||||
domains : Iterable of openmc.Cell, openmc.Material, or openmc.Universe
|
||||
Domains to find volumes of
|
||||
samples : int
|
||||
Number of samples used to generate volume estimates
|
||||
lower_left : Iterable of float
|
||||
|
|
@ -29,8 +29,10 @@ class VolumeCalculation(object):
|
|||
|
||||
Attributes
|
||||
----------
|
||||
cell_ids : Iterable of int
|
||||
IDs of cells to find volumes of
|
||||
ids : Iterable of int
|
||||
IDs of domains to find volumes of
|
||||
domain_type : {'cell', 'material', 'universe'}
|
||||
Type of each domain
|
||||
samples : int
|
||||
Number of samples used to generate volume estimates
|
||||
lower_left : Iterable of float
|
||||
|
|
@ -38,23 +40,31 @@ class VolumeCalculation(object):
|
|||
upper_right : Iterable of float
|
||||
Upper-right coordinates of bounding box used to sample points
|
||||
results : dict
|
||||
Dictionary whose keys are unique IDs of cells and values are
|
||||
Dictionary whose keys are unique IDs of domains and values are
|
||||
dictionaries with calculated volumes and total number of atoms for each
|
||||
nuclide present in the cell.
|
||||
nuclide present in the domain.
|
||||
volumes : dict
|
||||
Dictionary whose keys are unique IDs of cells and values are the
|
||||
Dictionary whose keys are unique IDs of domains and values are the
|
||||
estimated volumes
|
||||
atoms_dataframe : pandas.DataFrame
|
||||
DataFrame showing the estimated number of atoms for each nuclide present
|
||||
in each cell specified.
|
||||
in each domain specified.
|
||||
|
||||
"""
|
||||
def __init__(self, cells, samples, lower_left=None,
|
||||
def __init__(self, domains, samples, lower_left=None,
|
||||
upper_right=None):
|
||||
self._results = None
|
||||
|
||||
cv.check_type('cells', cells, Iterable, Cell)
|
||||
self.cell_ids = [c.id for c in cells]
|
||||
cv.check_type('domains', domains, Iterable,
|
||||
(openmc.Cell, openmc.Material, openmc.Universe))
|
||||
if isinstance(domains[0], openmc.Cell):
|
||||
self._domain_type = 'cell'
|
||||
elif isinstance(domains[0], openmc.Material):
|
||||
self._domain_type = 'material'
|
||||
elif isinstance(domains[0], openmc.Universe):
|
||||
self._domain_type = 'universe'
|
||||
self.ids = [d.id for d in domains]
|
||||
|
||||
self.samples = samples
|
||||
|
||||
if lower_left is not None:
|
||||
|
|
@ -64,17 +74,21 @@ class VolumeCalculation(object):
|
|||
'should be specified')
|
||||
self.upper_right = upper_right
|
||||
else:
|
||||
ll, ur = Union(*[c.region for c in cells]).bounding_box
|
||||
if np.any(np.isinf(ll)) or np.any(np.isinf(ur)):
|
||||
if self.domain_type == 'cell':
|
||||
ll, ur = openmc.Union(*[c.region for c in domains]).bounding_box
|
||||
if np.any(np.isinf(ll)) or np.any(np.isinf(ur)):
|
||||
raise ValueError('Could not automatically determine bounding '
|
||||
'box for stochastic volume calculation.')
|
||||
else:
|
||||
self.lower_left = ll
|
||||
self.upper_right = ur
|
||||
else:
|
||||
raise ValueError('Could not automatically determine bounding box '
|
||||
'for stochastic volume calculation.')
|
||||
else:
|
||||
self.lower_left = ll
|
||||
self.upper_right = ur
|
||||
|
||||
@property
|
||||
def cell_ids(self):
|
||||
return self._cell_ids
|
||||
def ids(self):
|
||||
return self._ids
|
||||
|
||||
@property
|
||||
def samples(self):
|
||||
|
|
@ -92,6 +106,10 @@ class VolumeCalculation(object):
|
|||
def results(self):
|
||||
return self._results
|
||||
|
||||
@property
|
||||
def domain_type(self):
|
||||
return self._domain_type
|
||||
|
||||
@property
|
||||
def volumes(self):
|
||||
return {uid: results['volume'] for uid, results in self.results.items()}
|
||||
|
|
@ -99,17 +117,18 @@ class VolumeCalculation(object):
|
|||
@property
|
||||
def atoms_dataframe(self):
|
||||
items = []
|
||||
columns = ['Cell', 'Nuclide', 'Atoms', 'Uncertainty']
|
||||
for cell_id, results in self.results.items():
|
||||
columns = [self.domain_type.capitalize(), 'Nuclide', 'Atoms',
|
||||
'Uncertainty']
|
||||
for uid, results in self.results.items():
|
||||
for name, atoms in results['atoms']:
|
||||
items.append((cell_id, name, atoms[0], atoms[1]))
|
||||
items.append((uid, name, atoms[0], atoms[1]))
|
||||
|
||||
return pd.DataFrame.from_records(items, columns=columns)
|
||||
|
||||
@cell_ids.setter
|
||||
def cell_ids(self, cell_ids):
|
||||
cv.check_type('cell IDs', cell_ids, Iterable, Real)
|
||||
self._cell_ids = cell_ids
|
||||
@ids.setter
|
||||
def ids(self, ids):
|
||||
cv.check_type('domain IDs', ids, Iterable, Real)
|
||||
self._ids = ids
|
||||
|
||||
@samples.setter
|
||||
def samples(self, samples):
|
||||
|
|
@ -154,16 +173,17 @@ class VolumeCalculation(object):
|
|||
import h5py
|
||||
|
||||
with h5py.File(filename, 'r') as f:
|
||||
domain_type = f.attrs['domain_type'].decode()
|
||||
samples = f.attrs['samples']
|
||||
lower_left = f.attrs['lower_left']
|
||||
upper_right = f.attrs['upper_right']
|
||||
|
||||
results = {}
|
||||
cell_ids = []
|
||||
ids = []
|
||||
for obj_name in f:
|
||||
if obj_name.startswith('cell_'):
|
||||
cell_id = int(obj_name[5:])
|
||||
cell_ids.append(cell_id)
|
||||
if obj_name.startswith('domain_'):
|
||||
domain_id = int(obj_name[7:])
|
||||
ids.append(domain_id)
|
||||
group = f[obj_name]
|
||||
volume = tuple(group['volume'].value)
|
||||
nucnames = group['nuclides'].value
|
||||
|
|
@ -172,14 +192,19 @@ class VolumeCalculation(object):
|
|||
atom_list = []
|
||||
for name_i, atoms_i in zip(nucnames, atoms):
|
||||
atom_list.append((name_i.decode(), tuple(atoms_i)))
|
||||
results[cell_id] = {'volume': volume, 'atoms': atom_list}
|
||||
results[domain_id] = {'volume': volume, 'atoms': atom_list}
|
||||
|
||||
# Instantiate some throw-away cells that are used by the constructor to
|
||||
# assign IDs
|
||||
cells = [Cell(uid) for uid in cell_ids]
|
||||
# Instantiate some throw-away domains that are used by the constructor
|
||||
# to assign IDs
|
||||
if domain_type == 'cell':
|
||||
domains = [openmc.Cell(uid) for uid in ids]
|
||||
elif domain_type == 'material':
|
||||
domains = [openmc.Material(uid) for uid in ids]
|
||||
elif domain_type == 'universe':
|
||||
domains = [openmc.Universe(uid) for uid in ids]
|
||||
|
||||
# Instantiate the class and assign results
|
||||
vol = cls(cells, samples, lower_left, upper_right)
|
||||
vol = cls(domains, samples, lower_left, upper_right)
|
||||
vol.results = results
|
||||
return vol
|
||||
|
||||
|
|
@ -193,8 +218,10 @@ class VolumeCalculation(object):
|
|||
|
||||
"""
|
||||
element = ET.Element("volume_calc")
|
||||
cell_elem = ET.SubElement(element, "cells")
|
||||
cell_elem.text = ' '.join(str(uid) for uid in self.cell_ids)
|
||||
dt_elem = ET.SubElement(element, "domain_type")
|
||||
dt_elem.text = self.domain_type
|
||||
id_elem = ET.SubElement(element, "domain_ids")
|
||||
id_elem.text = ' '.join(str(uid) for uid in self.ids)
|
||||
samples_elem = ET.SubElement(element, "samples")
|
||||
samples_elem.text = str(self.samples)
|
||||
ll_elem = ET.SubElement(element, "lower_left")
|
||||
|
|
|
|||
|
|
@ -143,8 +143,10 @@ element settings {
|
|||
element verbosity { xsd:positiveInteger }? &
|
||||
|
||||
element volume_calc {
|
||||
(element cells { list { xsd:positiveInteger+ } } |
|
||||
attribute cells { list { xsd:positiveInteger+ } }) &
|
||||
(element domain_type { xsd:string } |
|
||||
attribute domain_type { xsd:string }) &
|
||||
(element domain_ids { list { xsd:integer+ } } |
|
||||
attribute domain_ids { list { xsd:integer+ } }) &
|
||||
(element samples { xsd:positiveInteger } |
|
||||
attribute samples { xsd:positiveInteger }) &
|
||||
(element lower_left { list { xsd:double+ } } |
|
||||
|
|
|
|||
|
|
@ -629,17 +629,25 @@
|
|||
<element name="volume_calc">
|
||||
<interleave>
|
||||
<choice>
|
||||
<element name="cells">
|
||||
<element name="domain_type">
|
||||
<data type="string"/>
|
||||
</element>
|
||||
<attribute name="domain_type">
|
||||
<data type="string"/>
|
||||
</attribute>
|
||||
</choice>
|
||||
<choice>
|
||||
<element name="domain_ids">
|
||||
<list>
|
||||
<oneOrMore>
|
||||
<data type="positiveInteger"/>
|
||||
<data type="integer"/>
|
||||
</oneOrMore>
|
||||
</list>
|
||||
</element>
|
||||
<attribute name="cells">
|
||||
<attribute name="domain_ids">
|
||||
<list>
|
||||
<oneOrMore>
|
||||
<data type="positiveInteger"/>
|
||||
<data type="integer"/>
|
||||
</oneOrMore>
|
||||
</list>
|
||||
</attribute>
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@ module volume_calc
|
|||
use geometry, only: find_cell
|
||||
use global
|
||||
use hdf5_interface, only: file_create, file_close, write_attribute, &
|
||||
create_group, close_group, write_dataset
|
||||
create_group, close_group, write_dataset, write_attribute_string
|
||||
use output, only: write_message, header
|
||||
use message_passing
|
||||
use particle_header, only: Particle
|
||||
|
|
@ -33,7 +33,8 @@ contains
|
|||
subroutine run_volume_calculations()
|
||||
integer :: i, j
|
||||
integer :: n
|
||||
real(8), allocatable :: volume(:,:) ! volume mean/stdev in each cell
|
||||
real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain
|
||||
character(10) :: domain_type
|
||||
character(MAX_FILE_LEN) :: filename ! filename for HDF5 file
|
||||
type(Timer) :: time_volume ! timer for volume calculation
|
||||
type(VectorInt), allocatable :: nuclide_vec(:) ! indices in nuclides array
|
||||
|
|
@ -46,7 +47,7 @@ contains
|
|||
end if
|
||||
|
||||
do i = 1, size(volume_calcs)
|
||||
n = size(volume_calcs(i) % cell_id)
|
||||
n = size(volume_calcs(i) % domain_id)
|
||||
allocate(nuclide_vec(n))
|
||||
allocate(atoms_vec(n), uncertainty_vec(n))
|
||||
allocate(volume(2,n))
|
||||
|
|
@ -60,11 +61,20 @@ contains
|
|||
uncertainty_vec)
|
||||
|
||||
if (master) then
|
||||
! Display cell volumes
|
||||
do j = 1, size(volume_calcs(i) % cell_id)
|
||||
call write_message(" Cell " // trim(to_str(volume_calcs(i) % &
|
||||
cell_id(j))) // ": " // trim(to_str(volume(1,j))) // " +/- " // &
|
||||
trim(to_str(volume(2,j))) // " cm^3")
|
||||
select case (volume_calcs(i) % domain_type)
|
||||
case (FILTER_CELL)
|
||||
domain_type = ' Cell'
|
||||
case (FILTER_MATERIAL)
|
||||
domain_type = ' Material'
|
||||
case (FILTER_UNIVERSE)
|
||||
domain_type = ' Universe'
|
||||
end select
|
||||
|
||||
! Display domain volumes
|
||||
do j = 1, size(volume_calcs(i) % domain_id)
|
||||
call write_message(trim(domain_type) // " " // trim(to_str(&
|
||||
volume_calcs(i) % domain_id(j))) // ": " // trim(to_str(&
|
||||
volume(1,j))) // " +/- " // trim(to_str(volume(2,j))) // " cm^3")
|
||||
end do
|
||||
call write_message("")
|
||||
|
||||
|
|
@ -85,13 +95,13 @@ contains
|
|||
end subroutine run_volume_calculations
|
||||
|
||||
!===============================================================================
|
||||
! GET_VOLUME stochastically determines the volume of a set of cells along with
|
||||
! the average number densities of nuclides within the cell
|
||||
! GET_VOLUME stochastically determines the volume of a set of domains along with
|
||||
! the average number densities of nuclides within the domain
|
||||
!===============================================================================
|
||||
|
||||
subroutine get_volume(this, volume, nuclide_vec, atoms_vec, uncertainty_vec)
|
||||
type(VolumeCalculation), intent(in) :: this
|
||||
real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each cell
|
||||
real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each domain
|
||||
type(VectorInt), intent(out) :: nuclide_vec(:) ! indices in nuclides array
|
||||
type(VectorReal), intent(out) :: atoms_vec(:) ! total # of atoms of each nuclide
|
||||
type(VectorReal), intent(out) :: uncertainty_vec(:) ! uncertainty of total # of atoms
|
||||
|
|
@ -99,22 +109,22 @@ contains
|
|||
! Variables that are private to each thread
|
||||
integer(8) :: i
|
||||
integer :: j, k
|
||||
integer :: i_cell ! index in cell_id array
|
||||
integer :: i_domain ! index in domain_id array
|
||||
integer :: i_material ! index in materials array
|
||||
integer :: level ! local coordinate level
|
||||
logical :: found_cell
|
||||
type(VectorInt) :: indices(size(this % cell_id)) ! List of material indices
|
||||
type(VectorInt) :: hits(size(this % cell_id)) ! Number of hits for each material
|
||||
type(VectorInt) :: indices(size(this % domain_id)) ! List of material indices
|
||||
type(VectorInt) :: hits(size(this % domain_id)) ! Number of hits for each material
|
||||
type(Particle) :: p
|
||||
|
||||
! Shared variables
|
||||
integer :: i_start, i_end ! Starting/ending sample for each process
|
||||
type(VectorInt) :: master_indices(size(this % cell_id))
|
||||
type(VectorInt) :: master_hits(size(this % cell_id))
|
||||
type(VectorInt) :: master_indices(size(this % domain_id))
|
||||
type(VectorInt) :: master_hits(size(this % domain_id))
|
||||
|
||||
! Variables used outside of parallel region
|
||||
integer :: i_nuclide ! index in nuclides array
|
||||
integer :: total_hits ! total hits for a single cell (summed over materials)
|
||||
integer :: total_hits ! total hits for a single domain (summed over materials)
|
||||
integer :: min_samples ! minimum number of samples per process
|
||||
integer :: remainder ! leftover samples from uneven divide
|
||||
#ifdef MPI
|
||||
|
|
@ -140,15 +150,15 @@ contains
|
|||
|
||||
call p % initialize()
|
||||
|
||||
!$omp parallel private(i, j, k, i_cell, i_material, level, found_cell, &
|
||||
!$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, &
|
||||
!$omp& indices, hits) firstprivate(p)
|
||||
|
||||
! Reset vectors -- this is really to get around a gfortran 4.6 bug. Ideally,
|
||||
! indices and hits would just be firstprivate but 4.6 complains because they
|
||||
! have allocatable components...
|
||||
do i_cell = 1, size(this % cell_id)
|
||||
call indices(i_cell) % clear()
|
||||
call hits(i_cell) % clear()
|
||||
do i_domain = 1, size(this % domain_id)
|
||||
call indices(i_domain) % clear()
|
||||
call hits(i_domain) % clear()
|
||||
end do
|
||||
|
||||
call prn_set_stream(STREAM_VOLUME)
|
||||
|
|
@ -174,32 +184,89 @@ contains
|
|||
call find_cell(p, found_cell)
|
||||
if (.not. found_cell) cycle
|
||||
|
||||
! Determine if point is within desired cell
|
||||
LEVEL_LOOP: do level = 1, p % n_coord
|
||||
CELL_CHECK_LOOP: do i_cell = 1, size(this % cell_id)
|
||||
if (cells(p % coord(level) % cell) % id == this % cell_id(i_cell)) then
|
||||
|
||||
! Determine what material this is
|
||||
i_material = p % material
|
||||
if (this % domain_type == FILTER_MATERIAL) then
|
||||
! ======================================================================
|
||||
! MATERIAL VOLUME
|
||||
|
||||
i_material = p % material
|
||||
MATERIAL_LOOP: do i_domain = 1, size(this % domain_id)
|
||||
if (i_material == materials(i_domain) % id) then
|
||||
! Check if we've already had a hit in this material and if so,
|
||||
! simply add one
|
||||
do j = 1, indices(i_cell) % size()
|
||||
if (indices(i_cell) % data(j) == i_material) then
|
||||
hits(i_cell) % data(j) = hits(i_cell) % data(j) + 1
|
||||
cycle CELL_CHECK_LOOP
|
||||
do j = 1, indices(i_domain) % size()
|
||||
if (indices(i_domain) % data(j) == i_material) then
|
||||
hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1
|
||||
cycle MATERIAL_LOOP
|
||||
end if
|
||||
end do
|
||||
|
||||
! If we make it here, that means we haven't yet had a hit in this
|
||||
! material. Add an entry to both the indices list and the hits list
|
||||
call indices(i_cell) % push_back(i_material)
|
||||
call hits(i_cell) % push_back(1)
|
||||
call indices(i_domain) % push_back(i_material)
|
||||
call hits(i_domain) % push_back(1)
|
||||
end if
|
||||
end do CELL_CHECK_LOOP
|
||||
end do LEVEL_LOOP
|
||||
end do MATERIAL_LOOP
|
||||
|
||||
elseif (this % domain_type == FILTER_CELL) THEN
|
||||
! ======================================================================
|
||||
! CELL VOLUME
|
||||
|
||||
do level = 1, p % n_coord
|
||||
CELL_LOOP: do i_domain = 1, size(this % domain_id)
|
||||
if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then
|
||||
|
||||
! Determine what material this is
|
||||
i_material = p % material
|
||||
|
||||
! Check if we've already had a hit in this material and if so,
|
||||
! simply add one
|
||||
do j = 1, indices(i_domain) % size()
|
||||
if (indices(i_domain) % data(j) == i_material) then
|
||||
hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1
|
||||
cycle CELL_LOOP
|
||||
end if
|
||||
end do
|
||||
|
||||
! If we make it here, that means we haven't yet had a hit in this
|
||||
! material. Add an entry to both the indices list and the hits list
|
||||
call indices(i_domain) % push_back(i_material)
|
||||
call hits(i_domain) % push_back(1)
|
||||
end if
|
||||
end do CELL_LOOP
|
||||
end do
|
||||
|
||||
elseif (this % domain_type == FILTER_UNIVERSE) then
|
||||
! ======================================================================
|
||||
! UNIVERSE VOLUME
|
||||
|
||||
do level = 1, p % n_coord
|
||||
UNIVERSE_LOOP: do i_domain = 1, size(this % domain_id)
|
||||
if (universes(p % coord(level) % universe) % id == &
|
||||
this % domain_id(i_domain)) then
|
||||
|
||||
! Determine what material this is
|
||||
i_material = p % material
|
||||
|
||||
! Check if we've already had a hit in this material and if so,
|
||||
! simply add one
|
||||
do j = 1, indices(i_domain) % size()
|
||||
if (indices(i_domain) % data(j) == i_material) then
|
||||
hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1
|
||||
cycle UNIVERSE_LOOP
|
||||
end if
|
||||
end do
|
||||
|
||||
! If we make it here, that means we haven't yet had a hit in this
|
||||
! material. Add an entry to both the indices list and the hits list
|
||||
call indices(i_domain) % push_back(i_material)
|
||||
call hits(i_domain) % push_back(1)
|
||||
end if
|
||||
end do UNIVERSE_LOOP
|
||||
end do
|
||||
|
||||
end if
|
||||
end do SAMPLE_LOOP
|
||||
!$omp end do
|
||||
!$omp end do
|
||||
|
||||
! ==========================================================================
|
||||
! REDUCE HITS ONTO MASTER THREAD
|
||||
|
|
@ -212,14 +279,14 @@ contains
|
|||
!$omp do ordered schedule(static)
|
||||
THREAD_LOOP: do i = 1, omp_get_num_threads()
|
||||
!$omp ordered
|
||||
do i_cell = 1, size(this % cell_id)
|
||||
INDEX_LOOP: do j = 1, indices(i_cell) % size()
|
||||
do i_domain = 1, size(this % domain_id)
|
||||
INDEX_LOOP: do j = 1, indices(i_domain) % size()
|
||||
! Check if this material has been added to the master list and if so,
|
||||
! accumulate the number of hits
|
||||
do k = 1, master_indices(i_cell) % size()
|
||||
if (indices(i_cell) % data(j) == master_indices(i_cell) % data(k)) then
|
||||
master_hits(i_cell) % data(k) = &
|
||||
master_hits(i_cell) % data(k) + hits(i_cell) % data(j)
|
||||
do k = 1, master_indices(i_domain) % size()
|
||||
if (indices(i_domain) % data(j) == master_indices(i_domain) % data(k)) then
|
||||
master_hits(i_domain) % data(k) = &
|
||||
master_hits(i_domain) % data(k) + hits(i_domain) % data(j)
|
||||
cycle INDEX_LOOP
|
||||
end if
|
||||
end do
|
||||
|
|
@ -227,8 +294,8 @@ contains
|
|||
! If we made it here, this means the material hasn't yet been added to
|
||||
! the master list, so add an entry to both the master indices and master
|
||||
! hits lists
|
||||
call master_indices(i_cell) % push_back(indices(i_cell) % data(j))
|
||||
call master_hits(i_cell) % push_back(hits(i_cell) % data(j))
|
||||
call master_indices(i_domain) % push_back(indices(i_domain) % data(j))
|
||||
call master_hits(i_domain) % push_back(hits(i_domain) % data(j))
|
||||
end do INDEX_LOOP
|
||||
end do
|
||||
!$omp end ordered
|
||||
|
|
@ -247,7 +314,7 @@ contains
|
|||
|
||||
volume_sample = product(this % upper_right - this % lower_left)
|
||||
|
||||
do i_cell = 1, size(this % cell_id)
|
||||
do i_domain = 1, size(this % domain_id)
|
||||
atoms(:, :) = ZERO
|
||||
total_hits = 0
|
||||
|
||||
|
|
@ -261,9 +328,9 @@ contains
|
|||
call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, &
|
||||
MPI_STATUS_IGNORE, mpi_err)
|
||||
do k = 0, n - 1
|
||||
do m = 1, master_indices(i_cell) % size()
|
||||
if (data(2*k + 1) == master_indices(i_cell) % data(m)) then
|
||||
master_hits(i_cell) % data(m) = master_hits(i_cell) % data(m) + &
|
||||
do m = 1, master_indices(i_domain) % size()
|
||||
if (data(2*k + 1) == master_indices(i_domain) % data(m)) then
|
||||
master_hits(i_domain) % data(m) = master_hits(i_domain) % data(m) + &
|
||||
data(2*k + 2)
|
||||
end if
|
||||
end do
|
||||
|
|
@ -272,12 +339,12 @@ contains
|
|||
end do
|
||||
#endif
|
||||
|
||||
do j = 1, master_indices(i_cell) % size()
|
||||
total_hits = total_hits + master_hits(i_cell) % data(j)
|
||||
f = real(master_hits(i_cell) % data(j), 8) / this % samples
|
||||
do j = 1, master_indices(i_domain) % size()
|
||||
total_hits = total_hits + master_hits(i_domain) % data(j)
|
||||
f = real(master_hits(i_domain) % data(j), 8) / this % samples
|
||||
var_f = f*(ONE - f) / this % samples
|
||||
|
||||
i_material = master_indices(i_cell) % data(j)
|
||||
i_material = master_indices(i_domain) % data(j)
|
||||
if (i_material == MATERIAL_VOID) cycle
|
||||
|
||||
associate (mat => materials(i_material))
|
||||
|
|
@ -293,9 +360,9 @@ contains
|
|||
end do
|
||||
|
||||
! Determine volume
|
||||
volume(1, i_cell) = real(total_hits, 8) / this % samples * volume_sample
|
||||
volume(2, i_cell) = sqrt(volume(1, i_cell) * (volume_sample - &
|
||||
volume(1, i_cell)) / this % samples)
|
||||
volume(1, i_domain) = real(total_hits, 8) / this % samples * volume_sample
|
||||
volume(2, i_domain) = sqrt(volume(1, i_domain) * (volume_sample - &
|
||||
volume(1, i_domain)) / this % samples)
|
||||
|
||||
! Determine total number of atoms. At this point, we have values in
|
||||
! atoms/b-cm. To get to atoms we multiple by 10^24 V.
|
||||
|
|
@ -307,19 +374,19 @@ contains
|
|||
! Convert full arrays to vectors
|
||||
do j = 1, size(nuclides)
|
||||
if (atoms(1, j) > ZERO) then
|
||||
call nuclide_vec(i_cell) % push_back(j)
|
||||
call atoms_vec(i_cell) % push_back(atoms(1, j))
|
||||
call uncertainty_vec(i_cell) % push_back(atoms(2, j))
|
||||
call nuclide_vec(i_domain) % push_back(j)
|
||||
call atoms_vec(i_domain) % push_back(atoms(1, j))
|
||||
call uncertainty_vec(i_domain) % push_back(atoms(2, j))
|
||||
end if
|
||||
end do
|
||||
|
||||
else
|
||||
#ifdef MPI
|
||||
n = master_indices(i_cell) % size()
|
||||
n = master_indices(i_domain) % size()
|
||||
allocate(data(2*n))
|
||||
do k = 0, n - 1
|
||||
data(2*k + 1) = master_indices(i_cell) % data(k + 1)
|
||||
data(2*k + 2) = master_hits(i_cell) % data(k + 1)
|
||||
data(2*k + 1) = master_indices(i_domain) % data(k + 1)
|
||||
data(2*k + 2) = master_hits(i_domain) % data(k + 1)
|
||||
end do
|
||||
|
||||
call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err)
|
||||
|
|
@ -340,7 +407,7 @@ contains
|
|||
uncertainty_vec)
|
||||
type(VolumeCalculation), intent(in) :: this
|
||||
character(*), intent(in) :: filename ! filename for HDF5 file
|
||||
real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each cell
|
||||
real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each domain
|
||||
type(VectorInt), intent(in) :: nuclide_vec(:) ! indices in nuclides array
|
||||
type(VectorReal), intent(in) :: atoms_vec(:) ! total # of atoms of each nuclide
|
||||
type(VectorReal), intent(in) :: uncertainty_vec(:) ! uncertainty of total # of atoms
|
||||
|
|
@ -357,15 +424,23 @@ contains
|
|||
file_id = file_create(filename)
|
||||
|
||||
! Write basic metadata
|
||||
select case (this % domain_type)
|
||||
case (FILTER_CELL)
|
||||
call write_attribute_string(file_id, ".", "domain_type", "cell")
|
||||
case (FILTER_MATERIAL)
|
||||
call write_attribute_string(file_id, ".", "domain_type", "material")
|
||||
case (FILTER_UNIVERSE)
|
||||
call write_attribute_string(file_id, ".", "domain_type", "universe")
|
||||
end select
|
||||
call write_attribute(file_id, "samples", this % samples)
|
||||
call write_attribute(file_id, "lower_left", this % lower_left)
|
||||
call write_attribute(file_id, "upper_right", this % upper_right)
|
||||
|
||||
do i = 1, size(this % cell_id)
|
||||
group_id = create_group(file_id, "cell_" // trim(to_str(&
|
||||
this % cell_id(i))))
|
||||
do i = 1, size(this % domain_id)
|
||||
group_id = create_group(file_id, "domain_" // trim(to_str(&
|
||||
this % domain_id(i))))
|
||||
|
||||
! Write volume for cell
|
||||
! Write volume for domain
|
||||
call write_dataset(group_id, "volume", volume(:, i))
|
||||
|
||||
! Create array of nuclide names from the vector
|
||||
|
|
|
|||
|
|
@ -1,12 +1,14 @@
|
|||
module volume_header
|
||||
|
||||
use error, only: fatal_error
|
||||
use constants, only: FILTER_CELL, FILTER_MATERIAL, FILTER_UNIVERSE
|
||||
use error, only: fatal_error
|
||||
use xml_interface
|
||||
|
||||
implicit none
|
||||
|
||||
type VolumeCalculation
|
||||
integer, allocatable :: cell_id(:)
|
||||
integer :: domain_type
|
||||
integer, allocatable :: domain_id(:)
|
||||
real(8) :: lower_left(3)
|
||||
real(8) :: upper_right(3)
|
||||
integer :: samples
|
||||
|
|
@ -20,16 +22,31 @@ contains
|
|||
class(VolumeCalculation), intent(out) :: this
|
||||
type(Node), pointer :: node_vol
|
||||
|
||||
integer :: num_cells
|
||||
integer :: num_domains
|
||||
character(10) :: temp_str
|
||||
|
||||
! Check domain type
|
||||
call get_node_value(node_vol, "domain_type", temp_str)
|
||||
select case (temp_str)
|
||||
case ('cell')
|
||||
this % domain_type = FILTER_CELL
|
||||
case ('material')
|
||||
this % domain_type = FILTER_MATERIAL
|
||||
case ('universe')
|
||||
this % domain_type = FILTER_UNIVERSE
|
||||
case default
|
||||
call fatal_error("Unrecognized domain type for stochastic volume &
|
||||
&calculation: " // trim(temp_str))
|
||||
end select
|
||||
|
||||
! Read cell IDs
|
||||
if (check_for_node(node_vol, "cells")) then
|
||||
num_cells = get_arraysize_integer(node_vol, "cells")
|
||||
if (check_for_node(node_vol, "domain_ids")) then
|
||||
num_domains = get_arraysize_integer(node_vol, "domain_ids")
|
||||
else
|
||||
call fatal_error("Must specify at least one cell for a volume calculation")
|
||||
end if
|
||||
allocate(this % cell_id(num_cells))
|
||||
call get_node_array(node_vol, "cells", this % cell_id)
|
||||
allocate(this % domain_id(num_domains))
|
||||
call get_node_array(node_vol, "domain_ids", this % domain_id)
|
||||
|
||||
! Read lower-left and upper-right bounding coordinates
|
||||
call get_node_array(node_vol, "lower_left", this % lower_left)
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
b2de6ac20ca2ca38b00d61adae707d05519f6130eea25ab050220bac005cb6d23fa5812bf876307e92e90a8806f0b371500c5611ab23f80b7c073f118eb31649
|
||||
102569289552d021b6803f404a0c17a9c17a40578fdba43a6ba08b77a731e0368fffa6a8a7abd48555167cb9997c6dba9ec5044c8593b12056957b7e3ec44ed0
|
||||
|
|
@ -1,7 +1,8 @@
|
|||
k-combined: 4.165451e-02 3.582531e-04
|
||||
Cell 1: 31.4693 +/- 0.0721 cm^3
|
||||
Cell 2: 2.0933 +/- 0.0310 cm^3
|
||||
Cell 3: 2.0486 +/- 0.0307 cm^3
|
||||
Volume calculation 0
|
||||
Domain 1: 31.4693 +/- 0.0721 cm^3
|
||||
Domain 2: 2.0933 +/- 0.0310 cm^3
|
||||
Domain 3: 2.0486 +/- 0.0307 cm^3
|
||||
Cell Nuclide Atoms Uncertainty
|
||||
0 1 U235.71c 3.481769e+23 7.979991e+20
|
||||
1 1 Mo99.71c 3.481769e+22 7.979991e+19
|
||||
|
|
@ -11,3 +12,20 @@ Cell 3: 2.0486 +/- 0.0307 cm^3
|
|||
5 3 H1.71c 1.369920e+23 2.051689e+21
|
||||
6 3 O16.71c 6.849599e+22 1.025844e+21
|
||||
7 3 B10.71c 6.849599e+18 1.025844e+17
|
||||
Volume calculation 1
|
||||
Domain 1: 4.1419 +/- 0.0426 cm^3
|
||||
Domain 2: 31.4693 +/- 0.0721 cm^3
|
||||
Material Nuclide Atoms Uncertainty
|
||||
0 1 H1.71c 2.769690e+23 2.850068e+21
|
||||
1 1 O16.71c 1.384845e+23 1.425034e+21
|
||||
2 1 B10.71c 1.384845e+19 1.425034e+17
|
||||
3 2 U235.71c 3.481769e+23 7.979991e+20
|
||||
4 2 Mo99.71c 3.481769e+22 7.979991e+19
|
||||
Volume calculation 2
|
||||
Domain 0: 35.6112 +/- 0.0664 cm^3
|
||||
Universe Nuclide Atoms Uncertainty
|
||||
0 0 H1.71c 2.769690e+23 2.850068e+21
|
||||
1 0 O16.71c 1.384845e+23 1.425034e+21
|
||||
2 0 B10.71c 1.384845e+19 1.425034e+17
|
||||
3 0 U235.71c 3.481769e+23 7.979991e+20
|
||||
4 0 Mo99.71c 3.481769e+22 7.979991e+19
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
import glob
|
||||
import sys
|
||||
sys.path.insert(0, os.pardir)
|
||||
from testing_harness import PyAPITestHarness
|
||||
|
|
@ -43,9 +44,12 @@ class VolumeTest(PyAPITestHarness):
|
|||
geometry.export_to_xml()
|
||||
|
||||
# Set up stochastic volume calculation
|
||||
vol_calc = openmc.VolumeCalculation(
|
||||
[inside_cyl, top_hemisphere, bottom_hemisphere],
|
||||
100000)
|
||||
ll, ur = openmc.Union(*[c.region for c in root.cells.values()]).bounding_box
|
||||
vol_calcs = [
|
||||
openmc.VolumeCalculation(list(root.cells.values()), 100000),
|
||||
openmc.VolumeCalculation([water, fuel], 100000, ll, ur),
|
||||
openmc.VolumeCalculation([root], 100000, ll, ur)
|
||||
]
|
||||
|
||||
# Define settings
|
||||
settings = openmc.Settings()
|
||||
|
|
@ -54,7 +58,7 @@ class VolumeTest(PyAPITestHarness):
|
|||
settings.inactive = 0
|
||||
settings.source = openmc.Source(space=openmc.stats.Box(
|
||||
[-1., -1., -5.], [1., 1., 5.]))
|
||||
settings.volume_calculations = vol_calc
|
||||
settings.volume_calculations = vol_calcs
|
||||
settings.export_to_xml()
|
||||
|
||||
def _get_results(self):
|
||||
|
|
@ -65,15 +69,18 @@ class VolumeTest(PyAPITestHarness):
|
|||
# Write out k-combined.
|
||||
outstr = 'k-combined: {:12.6e} {:12.6e}\n'.format(*sp.k_combined)
|
||||
|
||||
# Read volume calculation results
|
||||
vol = openmc.VolumeCalculation.from_hdf5(
|
||||
os.path.join(os.getcwd(), 'volume_1.h5'))
|
||||
for i, filename in enumerate(sorted(glob.glob(os.path.join(
|
||||
os.getcwd(), 'volume_*.h5')))):
|
||||
outstr += 'Volume calculation {}\n'.format(i)
|
||||
|
||||
# Write cell volumes and total # of atoms for each nuclide
|
||||
for cell_id, results in sorted(vol.results.items()):
|
||||
outstr += 'Cell {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format(
|
||||
cell_id, results['volume'])
|
||||
outstr += str(vol.atoms_dataframe) + '\n'
|
||||
# Read volume calculation results
|
||||
vol = openmc.VolumeCalculation.from_hdf5(filename)
|
||||
|
||||
# Write cell volumes and total # of atoms for each nuclide
|
||||
for uid, results in sorted(vol.results.items()):
|
||||
outstr += 'Domain {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format(
|
||||
uid, results['volume'])
|
||||
outstr += str(vol.atoms_dataframe) + '\n'
|
||||
|
||||
return outstr
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue