diff --git a/openmc/cell.py b/openmc/cell.py index c4a0952e9..fe8b4a52c 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -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.') diff --git a/openmc/geometry.py b/openmc/geometry.py index c827f4b2a..b2a9b5e4a 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -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. diff --git a/openmc/volume.py b/openmc/volume.py index 167971bdb..e9b1dac1d 100644 --- a/openmc/volume.py +++ b/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") diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 78a32171b..e23cbdd6c 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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+ } } | diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 64fe42239..38acea598 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -629,17 +629,25 @@ - + + + + + + + + + - + - + - + diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index eac96f9ba..03e241a66 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -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 diff --git a/src/volume_header.F90 b/src/volume_header.F90 index 3a34f1231..c70345f15 100644 --- a/src/volume_header.F90 +++ b/src/volume_header.F90 @@ -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) diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index 2a3f25e6d..e146703f3 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -1 +1 @@ -b2de6ac20ca2ca38b00d61adae707d05519f6130eea25ab050220bac005cb6d23fa5812bf876307e92e90a8806f0b371500c5611ab23f80b7c073f118eb31649 \ No newline at end of file +102569289552d021b6803f404a0c17a9c17a40578fdba43a6ba08b77a731e0368fffa6a8a7abd48555167cb9997c6dba9ec5044c8593b12056957b7e3ec44ed0 \ No newline at end of file diff --git a/tests/test_volume_calc/results_true.dat b/tests/test_volume_calc/results_true.dat index e7a39a931..da6dfd2af 100644 --- a/tests/test_volume_calc/results_true.dat +++ b/tests/test_volume_calc/results_true.dat @@ -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 diff --git a/tests/test_volume_calc/test_volume_calc.py b/tests/test_volume_calc/test_volume_calc.py index a6d779a25..e2c3eba79 100644 --- a/tests/test_volume_calc/test_volume_calc.py +++ b/tests/test_volume_calc/test_volume_calc.py @@ -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