From 9f3590cd39e018b7412cc60ea37d057e414e4982 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 20 Dec 2015 18:40:07 -0500 Subject: [PATCH 01/17] Make cell % material an allocatable array --- src/geometry.F90 | 2 +- src/geometry_header.F90 | 14 ++++++++------ src/initialize.F90 | 6 +++--- src/input_xml.F90 | 13 +++++++------ src/plot.F90 | 6 +++--- src/summary.F90 | 4 ++-- 6 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 9a084a77cb..19289f9a45 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -248,7 +248,7 @@ contains ! set material p % last_material = p % material - p % material = c % material + p % material = c % material(1) elseif (c % type == CELL_FILL) then CELL_TYPE ! ====================================================================== diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 6ca13740ba..30d4a119b6 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -126,15 +126,17 @@ module geometry_header integer :: fill ! universe # filling this cell integer :: instances ! number of instances of this cell in ! the geom - integer :: material ! Material within cell (0 for - ! universe) - integer, allocatable :: offset (:) ! Distribcell offset for tally + integer, allocatable :: material(:) ! Material within cell. Multiple + ! materials for distribcell + ! instances. 0 signifies a universe + integer, allocatable :: offset(:) ! Distribcell offset for tally ! counter integer, allocatable :: region(:) ! Definition of spatial region as - ! Boolean expression of half-spaces + ! Boolean expression of half-spaces integer, allocatable :: rpn(:) ! Reverse Polish notation for region - ! expression - logical :: simple ! Is the region simple (intersections only) + ! expression + logical :: simple ! Is the region simple (intersections + ! only) ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index 9d63eb4e93..802f890aa9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -614,13 +614,13 @@ contains ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - id = c%material + id = c%material(1) if (id == MATERIAL_VOID) then c%type = CELL_NORMAL elseif (id /= 0) then if (material_dict%has_key(id)) then c%type = CELL_NORMAL - c%material = material_dict%get_key(id) + c%material(1) = material_dict%get_key(id) else call fatal_error("Could not find material " // trim(to_str(id)) & &// " specified on cell " // trim(to_str(c%id))) @@ -1137,7 +1137,7 @@ contains ! Allocate offset table for fill cells do i = 1, n_cells - if (cells(i)%material == NONE) then + if (cells(i) % type /= CELL_NORMAL) then allocate(cells(i)%offset(n_maps)) end if end do diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5449b169dd..016b811ac4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1085,36 +1085,37 @@ contains end if ! Read material + allocate(c % material(1)) word = '' if (check_for_node(node_cell, "material")) & call get_node_value(node_cell, "material", word) select case(to_lower(word)) case ('void') - c % material = MATERIAL_VOID + c % material(1) = MATERIAL_VOID case ('') ! This case is called if no material was specified - c % material = NONE + c % material(1) = NONE case default - c % material = int(str_to_int(word), 4) + c % material(1) = int(str_to_int(word), 4) ! Check for error - if (c % material == ERROR_INT) then + if (c % material(1) == ERROR_INT) then call fatal_error("Invalid material specified on cell " & &// to_str(c % id)) end if end select ! Check to make sure that either material or fill was specified - if (c % material == NONE .and. c % fill == NONE) then + if (c % material(1) == NONE .and. c % fill == NONE) then call fatal_error("Neither material nor fill was specified for cell " & &// trim(to_str(c % id))) end if ! Check to make sure that both material and fill haven't been ! specified simultaneously - if (c % material /= NONE .and. c % fill /= NONE) then + if (c % material(1) /= NONE .and. c % fill /= NONE) then call fatal_error("Cannot specify material and fill simultaneously") end if diff --git a/src/plot.F90 b/src/plot.F90 index a5497bc203..99d7e3ae94 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -82,7 +82,7 @@ contains if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material c => cells(p % coord(j) % cell) - if (c % material == MATERIAL_VOID) then + if (c % material(1) == MATERIAL_VOID) then ! By default, color void cells white rgb = 255 id = -1 @@ -91,8 +91,8 @@ contains rgb = pl % not_found % rgb id = -1 else - rgb = pl % colors(c % material) % rgb - id = materials(c % material) % id + rgb = pl % colors(c % material(1)) % rgb + id = materials(c % material(1)) % id end if else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell diff --git a/src/summary.F90 b/src/summary.F90 index f2c261ec75..e9c809cf0c 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -150,10 +150,10 @@ contains select case (c%type) case (CELL_NORMAL) call write_dataset(cell_group, "fill_type", "normal") - if (c%material == MATERIAL_VOID) then + if (c%material(1) == MATERIAL_VOID) then call write_dataset(cell_group, "material", -1) else - call write_dataset(cell_group, "material", materials(c%material)%id) + call write_dataset(cell_group, "material", materials(c%material(1))%id) end if case (CELL_FILL) From fef4a50c5e817a7aa7907b3afc724e9c0f3b6e47 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 28 Dec 2015 21:52:34 -0500 Subject: [PATCH 02/17] Build distribcell maps for cells with distribmats --- src/geometry.F90 | 28 +++++++++++++++++++++-- src/initialize.F90 | 57 +++++++++++++++++++++++++++------------------- src/input_xml.F90 | 51 +++++++++++++++++++++++++---------------- 3 files changed, 92 insertions(+), 44 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 19289f9a45..f386dcf008 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -194,6 +194,7 @@ contains integer, optional :: search_cells(:) integer :: i ! index over cells integer :: j ! coordinate level index + integer :: k ! distribcell instance integer :: i_xyz(3) ! indices in lattice integer :: n ! number of cells to search integer :: index_cell ! index in cells array @@ -246,9 +247,32 @@ contains ! ====================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH - ! set material + ! Set the particle material p % last_material = p % material - p % material = c % material(1) +! if (size(c % material) == 1) then + ! Only one material for this cell; assign that one to the particle. + p % material = c % material(1) +! else +! ! Distributed instances of this cell have different materials. +! ! Determine which instance this is and assign the matching material. +! offset = 0 +! do k = 1, p % n_coord +! if (cells(p % coord(k) % cell) % type == CELL_FILL) then +! offset = offset + cells(p % coord(k) % cell) % & +! offset(t % filters(i) % offset) +! elseif(cells(p % coord(k) % cell) % type == CELL_LATTICE) then +! if (lattices(p % coord(k + 1) % lattice) % obj & +! % are_valid_indices([& +! p % coord(k + 1) % lattice_x, & +! p % coord(k + 1) % lattice_y, & +! p % coord(k + 1) % lattice_z])) then +! offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & +! offset(t % filters(i) % offset, & +! p % coord(k + 1) % lattice_x, & +! p % coord(k + 1) % lattice_y, & +! p % coord(k + 1) % lattice_z) +! end if +! end if elseif (c % type == CELL_FILL) then CELL_TYPE ! ====================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index 802f890aa9..f28959618f 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -614,31 +614,33 @@ contains ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - id = c%material(1) - if (id == MATERIAL_VOID) then - c%type = CELL_NORMAL - elseif (id /= 0) then - if (material_dict%has_key(id)) then - c%type = CELL_NORMAL - c%material(1) = material_dict%get_key(id) - else - call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c%id))) - end if - else - id = c%fill - if (universe_dict%has_key(id)) then - c%type = CELL_FILL - c%fill = universe_dict%get_key(id) - elseif (lattice_dict%has_key(id)) then - lid = lattice_dict%get_key(id) - c%type = CELL_LATTICE - c%fill = lid + if (c % material(1) == NONE) then + id = c % fill + if (universe_dict % has_key(id)) then + c % type = CELL_FILL + c % fill = universe_dict % get_key(id) + elseif (lattice_dict % has_key(id)) then + lid = lattice_dict % get_key(id) + c % type = CELL_LATTICE + c % fill = lid else call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "& - &// trim(to_str(c%id)) // " is neither a universe nor a & + &// trim(to_str(c % id)) // " is neither a universe nor a & &lattice.") end if + else + do j = 1, size(c % material) + id = c % material(j) + if (id == MATERIAL_VOID) then + c % type = CELL_NORMAL + else if (material_dict % has_key(id)) then + c % type = CELL_NORMAL + c % material(j) = material_dict % get_key(id) + else + call fatal_error("Could not find material " // trim(to_str(id)) & + &// " specified on cell " // trim(to_str(c % id))) + end if + end do end if end do @@ -1042,11 +1044,12 @@ contains class(Lattice), pointer :: lat ! pointer to lattice type(TallyObject), pointer :: t ! pointer to tally type(TallyFilter), pointer :: filter ! pointer to filter + type(Cell), pointer :: c ! pointer to cell ! Begin gathering list of cells in distribcell tallies n_maps = 0 - ! Populate list of distribcells to track + ! List all cells referenced in distribcell filters. do i = 1, n_tallies t => tallies(i) @@ -1062,6 +1065,14 @@ contains end do end do + ! List all cells with multiple (distributed) materials. + do i = 1, n_cells + c => cells(i) + if (size(c % material) > 1) then + if (.not. cell_list % contains(i)) call cell_list % add(i) + end if + end do + ! Compute the number of unique universes containing these distribcells ! to determine the number of offset tables to allocate do i = 1, n_universes @@ -1076,7 +1087,7 @@ contains ! Allocate the list of offset tables for each unique universe allocate(univ_list(n_maps)) - ! Allocate list to accumulate target distribccell counts in each universe + ! Allocate list to accumulate target distribcell counts in each universe allocate(counts(n_universes, n_maps)) ! Allocate list to track if target distribcells are found in each universe diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 016b811ac4..24b36f3fb8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -982,8 +982,7 @@ contains subroutine read_geometry_xml() integer :: i, j, k, m, i_x, i_a, input_index - integer :: n - integer :: n_x, n_y, n_z, n_rings, n_rlats, n_hlats + integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats integer :: universe_num integer :: n_cells_in_univ integer :: coeffs_reqd @@ -994,6 +993,7 @@ contains logical :: boundary_exists character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word + character(MAX_WORD_LEN), allocatable :: sarray(:) character(1000) :: region_spec type(Cell), pointer :: c class(Surface), pointer :: s @@ -1085,27 +1085,40 @@ contains end if ! Read material - allocate(c % material(1)) - word = '' - if (check_for_node(node_cell, "material")) & - call get_node_value(node_cell, "material", word) - select case(to_lower(word)) - case ('void') - c % material(1) = MATERIAL_VOID + if (check_for_node(node_cell, "material")) then + n_mats = get_arraysize_string(node_cell, "material") - case ('') - ! This case is called if no material was specified - c % material(1) = NONE + if (n_mats > 0) then + allocate(sarray(n_mats)) + call get_node_array(node_cell, "material", sarray) - case default - c % material(1) = int(str_to_int(word), 4) + allocate(c % material(n_mats)) + do j = 1, n_mats + select case(trim(to_lower(sarray(j)))) + case ('void') + c % material(j) = MATERIAL_VOID + case default + c % material(j) = int(str_to_int(sarray(j)), 4) - ! Check for error - if (c % material(1) == ERROR_INT) then - call fatal_error("Invalid material specified on cell " & - &// to_str(c % id)) + ! Check for error + if (c % material(j) == ERROR_INT) then + call fatal_error("Invalid material specified on cell " & + &// to_str(c % id)) + end if + end select + end do + + deallocate(sarray) + + else + allocate(c % material(1)) + c % material(1) = NONE end if - end select + + else + allocate(c % material(1)) + c % material(1) = NONE + end if ! Check to make sure that either material or fill was specified if (c % material(1) == NONE .and. c % fill == NONE) then From 20d6cc471a46ccb442bd2839c58abc5997b8c19c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 18:20:24 -0500 Subject: [PATCH 03/17] Implement distribmat in find_cell --- src/geometry.F90 | 53 ++++++++++++++++++++++++++---------------------- src/plot.F90 | 14 ++++++------- 2 files changed, 36 insertions(+), 31 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index f386dcf008..571ed6d407 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -193,8 +193,9 @@ contains logical, intent(inout) :: found integer, optional :: search_cells(:) integer :: i ! index over cells - integer :: j ! coordinate level index - integer :: k ! distribcell instance + integer :: j, k ! coordinate level index + integer :: offset ! instance # of a distributed cell + integer :: distribcell_index integer :: i_xyz(3) ! indices in lattice integer :: n ! number of cells to search integer :: index_cell ! index in cells array @@ -249,30 +250,34 @@ contains ! Set the particle material p % last_material = p % material -! if (size(c % material) == 1) then + if (size(c % material) == 1) then ! Only one material for this cell; assign that one to the particle. p % material = c % material(1) -! else -! ! Distributed instances of this cell have different materials. -! ! Determine which instance this is and assign the matching material. -! offset = 0 -! do k = 1, p % n_coord -! if (cells(p % coord(k) % cell) % type == CELL_FILL) then -! offset = offset + cells(p % coord(k) % cell) % & -! offset(t % filters(i) % offset) -! elseif(cells(p % coord(k) % cell) % type == CELL_LATTICE) then -! if (lattices(p % coord(k + 1) % lattice) % obj & -! % are_valid_indices([& -! p % coord(k + 1) % lattice_x, & -! p % coord(k + 1) % lattice_y, & -! p % coord(k + 1) % lattice_z])) then -! offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & -! offset(t % filters(i) % offset, & -! p % coord(k + 1) % lattice_x, & -! p % coord(k + 1) % lattice_y, & -! p % coord(k + 1) % lattice_z) -! end if -! end if + else + ! Distributed instances of this cell have different materials. + ! Determine which instance this is and assign the matching material. + distribcell_index = c % distribcell_index + offset = 0 + do k = 1, p % n_coord + if (cells(p % coord(k) % cell) % type == CELL_FILL) then + offset = offset + cells(p % coord(k) % cell) % & + offset(distribcell_index) + elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then + if (lattices(p % coord(k + 1) % lattice) % obj & + % are_valid_indices([& + p % coord(k + 1) % lattice_x, & + p % coord(k + 1) % lattice_y, & + p % coord(k + 1) % lattice_z])) then + offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & + offset(distribcell_index, & + p % coord(k + 1) % lattice_x, & + p % coord(k + 1) % lattice_y, & + p % coord(k + 1) % lattice_z) + end if + end if + end do + p % material = c % material(offset + 1) + end if elseif (c % type == CELL_FILL) then CELL_TYPE ! ====================================================================== diff --git a/src/plot.F90 b/src/plot.F90 index 99d7e3ae94..75d0bc9442 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -82,17 +82,17 @@ contains if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material c => cells(p % coord(j) % cell) - if (c % material(1) == MATERIAL_VOID) then - ! By default, color void cells white - rgb = 255 - id = -1 - else if (c % type == CELL_FILL) then + if (c % type == CELL_FILL) then ! If we stopped on a middle universe level, treat as if not found rgb = pl % not_found % rgb id = -1 + else if (p % material == MATERIAL_VOID) then + ! By default, color void cells white + rgb = 255 + id = -1 else - rgb = pl % colors(c % material(1)) % rgb - id = materials(c % material(1)) % id + rgb = pl % colors(p % material) % rgb + id = materials(p % material) % id end if else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell From 60732fd7dfeacbc5f92e2b769a19caf10658d5d8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 19:36:53 -0500 Subject: [PATCH 04/17] Enable distributed materials in PyAPI --- openmc/geometry.py | 7 +++++-- openmc/universe.py | 12 +++++++++++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/openmc/geometry.py b/openmc/geometry.py index fc8e19f071..9788671f5c 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -1,4 +1,4 @@ -from collections import OrderedDict +from collections import Iterable, OrderedDict from xml.etree import ElementTree as ET import openmc @@ -140,7 +140,10 @@ class Geometry(object): materials = set() for cell in material_cells: - materials.add(cell._fill) + if isinstance(cell.fill, Iterable): + for m in cell.fill: materials.add(m) + else: + materials.add(cell.fill) materials = list(materials) materials.sort(key=lambda x: x.id) diff --git a/openmc/universe.py b/openmc/universe.py index 8417f27c71..2fa535265c 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -50,7 +50,7 @@ class Cell(object): Unique identifier for the cell name : str Name of the cell - fill : Material or Universe or Lattice or 'void' + fill : Materials or Universe or Lattice or 'void' Indicates what the region of space is filled with region : openmc.region.Region Region of space that is assigned to the cell. @@ -112,6 +112,9 @@ class Cell(object): if isinstance(self._fill, openmc.Material): string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', self._fill._id) + elif isinstance(self._fill, Iterable): + string += ('{0: <16}{1}'.format('\tMaterial', '=\t') + + '[' + ', '.join([str(m.id) for m in self.fill]) + ']\n') elif isinstance(self._fill, (Universe, Lattice)): string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill._id) @@ -205,6 +208,10 @@ class Cell(object): elif isinstance(fill, openmc.Material): self._type = 'normal' + elif isinstance(fill, Iterable): + cv.check_type('cell.fill', fill, Iterable, openmc.Material) + self._type = 'normal' + elif isinstance(fill, Universe): self._type = 'fill' @@ -394,6 +401,9 @@ class Cell(object): if isinstance(self._fill, openmc.Material): element.set("material", str(self._fill._id)) + elif isinstance(self._fill, Iterable): + element.set("material", ' '.join([str(m.id) for m in self.fill])) + elif isinstance(self._fill, (Universe, Lattice)): element.set("fill", str(self._fill._id)) self._fill.create_xml_subelement(xml_element) From 12844d0dd1b4cf0c6ed279025eaae3f297cb2b8f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 20:42:40 -0500 Subject: [PATCH 05/17] Add distributed material test --- tests/test_distribmat/inputs_true.dat | 1 + tests/test_distribmat/results_true.dat | 2 + tests/test_distribmat/test_distribmat.py | 122 +++++++++++++++++++++++ 3 files changed, 125 insertions(+) create mode 100644 tests/test_distribmat/inputs_true.dat create mode 100644 tests/test_distribmat/results_true.dat create mode 100644 tests/test_distribmat/test_distribmat.py diff --git a/tests/test_distribmat/inputs_true.dat b/tests/test_distribmat/inputs_true.dat new file mode 100644 index 0000000000..72d296c9b1 --- /dev/null +++ b/tests/test_distribmat/inputs_true.dat @@ -0,0 +1 @@ +dcb9a612432305763ad10aaaab3ed095716cbb73f6b357ec180d300ccf365f9781aac7b1f3228c8756ceccf6aa230d0e03cf18ef2c0bc21fba1f082c3d5d6896 \ No newline at end of file diff --git a/tests/test_distribmat/results_true.dat b/tests/test_distribmat/results_true.dat new file mode 100644 index 0000000000..9481371d85 --- /dev/null +++ b/tests/test_distribmat/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.433669E+00 7.069157E-03 diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py new file mode 100644 index 0000000000..b6bffd90d9 --- /dev/null +++ b/tests/test_distribmat/test_distribmat.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness, PyAPITestHarness +import openmc + + +class DistribmatTestHarness(PyAPITestHarness): + def _build_inputs(self): + #################### + # Materials + #################### + + moderator = openmc.Material(material_id=1) + moderator.set_density('g/cc', 1.0) + moderator.add_nuclide('H-1', 2.0) + moderator.add_nuclide('O-16', 1.0) + + dense_fuel = openmc.Material(material_id=2) + dense_fuel.set_density('g/cc', 4.5) + dense_fuel.add_nuclide('U-235', 1.0) + + light_fuel = openmc.Material(material_id=3) + light_fuel.set_density('g/cc', 2.0) + light_fuel.add_nuclide('U-235', 1.0) + + mats_file = openmc.MaterialsFile() + mats_file.default_xs = '71c' + mats_file.add_materials([moderator, dense_fuel, light_fuel]) + mats_file.export_to_xml() + + + #################### + # Geometry + #################### + + c1 = openmc.Cell(cell_id=1) + c1.fill = moderator + mod_univ = openmc.Universe(universe_id=1) + mod_univ.add_cell(c1) + + r0 = openmc.ZCylinder(R=0.3) + c11 = openmc.Cell(cell_id=11) + c11.region = -r0 + c11.fill = [dense_fuel, light_fuel] + [dense_fuel]*2 + c12 = openmc.Cell(cell_id=12) + c12.region = +r0 + c12.fill = moderator + fuel_univ = openmc.Universe(universe_id=11) + fuel_univ.add_cells((c11, c12)) + + lat = openmc.RectLattice(lattice_id=101) + lat.dimension = [2, 2] + lat.lower_left = [-2.0, -2.0] + lat.pitch = [2.0, 2.0] + lat.universes = [[fuel_univ]*2]*2 + lat.outer = mod_univ + + x0 = openmc.XPlane(x0=-3.0) + x1 = openmc.XPlane(x0=3.0) + y0 = openmc.YPlane(y0=-3.0) + y1 = openmc.YPlane(y0=3.0) + for s in [x0, x1, y0, y1]: + s.boundary_type = 'reflective' + c101 = openmc.Cell(cell_id=101) + c101.region = +x0 & -x1 & +y0 & -y1 + c101.fill = lat + root_univ = openmc.Universe(universe_id=0) + root_univ.add_cell(c101) + + geometry = openmc.Geometry() + geometry.root_universe = root_univ + geo_file = openmc.GeometryFile() + geo_file.geometry = geometry + geo_file.export_to_xml() + + + #################### + # Settings + #################### + + sets_file = openmc.SettingsFile() + sets_file.batches = 5 + sets_file.inactive = 0 + sets_file.particles = 1000 + sets_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) + #sets_file.output = {'summary': True} + sets_file.export_to_xml() + + +# #################### +# # Plots +# #################### +# +# plots_file = openmc.PlotsFile() +# +# plot = openmc.Plot(plot_id=1) +# plot.basis = 'xy' +# plot.color = 'cell' +# plot.filename = 'cellplot' +# plot.origin = (0, 0, 0) +# plot.width = (7, 7) +# plot.pixels = (400, 400) +# plots_file.add_plot(plot) +# +# plot = openmc.Plot(plot_id=2) +# plot.basis = 'xy' +# plot.color = 'mat' +# plot.filename = 'matplot' +# plot.origin = (0, 0, 0) +# plot.width = (7, 7) +# plot.pixels = (400, 400) +# plots_file.add_plot(plot) +# +# plots_file.export_to_xml() + + +if __name__ == '__main__': + harness = DistribmatTestHarness('statepoint.5.*') + harness.main() From d5d486536bc7be21b40b782714b97c330d0819f2 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 23:02:29 -0500 Subject: [PATCH 06/17] Update summary.h5 files for distribmat --- openmc/summary.py | 11 ++++++++--- src/initialize.F90 | 2 +- src/summary.F90 | 13 ++++++++++--- tests/test_distribmat/inputs_true.dat | 2 +- tests/test_distribmat/results_true.dat | 9 +++++++++ tests/test_distribmat/test_distribmat.py | 11 ++++++++++- 6 files changed, 39 insertions(+), 9 deletions(-) diff --git a/openmc/summary.py b/openmc/summary.py index 46fec01d50..ad1dd67cae 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -1,3 +1,4 @@ +from collections import Iterable import numpy as np import re @@ -477,10 +478,14 @@ class Summary(object): # Retrieve the object corresponding to the fill type and ID if fill_type == 'normal': - if fill_id > 0: - fill = self.get_material_by_id(fill_id) + if isinstance(fill_id, Iterable): + fill = [self.get_material_by_id(mat) if mat > 0 else 'void' + for mat in fill_id] else: - fill = 'void' + if fill_id > 0: + fill = self.get_material_by_id(fill_id) + else: + fill = 'void' elif fill_type == 'universe': fill = self.get_universe_by_id(fill_id) else: diff --git a/src/initialize.F90 b/src/initialize.F90 index 7457be90d8..2c44755d20 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1038,7 +1038,7 @@ contains integer, intent(out), allocatable :: counts(:,:) ! Target count logical, intent(out), allocatable :: found(:,:) ! Target found - integer :: i, j, k, l, m ! Loop counters + integer :: i, j, k ! Loop counters type(SetInt) :: cell_list ! distribells to track type(Universe), pointer :: univ ! pointer to universe class(Lattice), pointer :: lat ! pointer to lattice diff --git a/src/summary.F90 b/src/summary.F90 index 3fcc3e5c91..deccaddf59 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -107,6 +107,7 @@ contains integer :: i, j, k, m integer, allocatable :: lattice_universes(:,:,:) + integer, allocatable :: cell_materials(:) integer(HID_T) :: geom_group integer(HID_T) :: cells_group, cell_group integer(HID_T) :: surfaces_group, surface_group @@ -150,10 +151,16 @@ contains select case (c%type) case (CELL_NORMAL) call write_dataset(cell_group, "fill_type", "normal") - if (c%material(1) == MATERIAL_VOID) then - call write_dataset(cell_group, "material", -1) + if (size(c % material) == 1) then + call write_dataset(cell_group, "material", & + materials(c % material(1)) % id) else - call write_dataset(cell_group, "material", materials(c%material(1))%id) + allocate(cell_materials(size(c % material))) + do j = 1, size(c % material) + cell_materials(j) = materials(c % material(j)) % id + end do + call write_dataset(cell_group, "material", cell_materials) + deallocate(cell_materials) end if case (CELL_FILL) diff --git a/tests/test_distribmat/inputs_true.dat b/tests/test_distribmat/inputs_true.dat index 72d296c9b1..48044e2872 100644 --- a/tests/test_distribmat/inputs_true.dat +++ b/tests/test_distribmat/inputs_true.dat @@ -1 +1 @@ -dcb9a612432305763ad10aaaab3ed095716cbb73f6b357ec180d300ccf365f9781aac7b1f3228c8756ceccf6aa230d0e03cf18ef2c0bc21fba1f082c3d5d6896 \ No newline at end of file +24a3537446feafcf79fc84d268b43130e9521ad7cccfcf733ba35ca50bf3c644e0b8e91bfdef7bdb4073783197cfdd730d1e459fc622aa3ce6dcf3143c805a1f \ No newline at end of file diff --git a/tests/test_distribmat/results_true.dat b/tests/test_distribmat/results_true.dat index 9481371d85..5e93e4e473 100644 --- a/tests/test_distribmat/results_true.dat +++ b/tests/test_distribmat/results_true.dat @@ -1,2 +1,11 @@ k-combined: 1.433669E+00 7.069157E-03 +Cell + ID = 11 + Name = + Material = [2, 3, 2, 2] + Region = -10000 + Rotation = None + Translation = None + Offset = None + Distribcell index= 1 diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py index b6bffd90d9..61e40e64c6 100644 --- a/tests/test_distribmat/test_distribmat.py +++ b/tests/test_distribmat/test_distribmat.py @@ -86,7 +86,7 @@ class DistribmatTestHarness(PyAPITestHarness): sets_file.inactive = 0 sets_file.particles = 1000 sets_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) - #sets_file.output = {'summary': True} + sets_file.output = {'summary': True} sets_file.export_to_xml() @@ -116,6 +116,15 @@ class DistribmatTestHarness(PyAPITestHarness): # # plots_file.export_to_xml() + def _get_results(self): + outstr = super(DistribmatTestHarness, self)._get_results() + su = openmc.Summary('summary.h5') + outstr += str(su.get_cell_by_id(11)) + return outstr + +# def _cleanup(self): +# return None + if __name__ == '__main__': harness = DistribmatTestHarness('statepoint.5.*') From e2a8b96cfdf0a4bb218c9a7eecc32490cb679074 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 23:34:28 -0500 Subject: [PATCH 07/17] Update docs for distribmat --- docs/source/usersguide/input.rst | 4 +++- docs/source/usersguide/output/summary.rst | 8 +++++--- src/relaxng/geometry.rnc | 4 ++-- src/relaxng/geometry.rng | 20 ++++++++++++-------- 4 files changed, 22 insertions(+), 14 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index fa43ca1d2c..c4538b889f 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -891,7 +891,9 @@ Each ```` element can have the following attributes or sub-elements: :material: The ``id`` of the material that this cell contains. If the cell should - contain no material, this can also be set to "void". + contain no material, this can also be set to "void". A list of materials + can be specified for the "distributed material" feature. This will give each + unique instance of the cell its own material. .. note:: If a material is specified, no fill should be given. diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index 66a5cad154..d100fab3bf 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -91,10 +91,12 @@ The current revision of the summary file format is 1. Type of fill for the cell. Can be 'normal', 'universe', or 'lattice'. -**/geometry/cells/cell /material** (*int*) +**/geometry/cells/cell /material** (*int* or *int[]*) - Unique ID of the material assigned to the cell. This dataset is present only - if fill_type is set to 'normal'. + Unique ID of the material(s) assigned to the cell. This dataset is present + only if fill_type is set to 'normal'. The data is an array if the cell uses + distributed materials, otherwise it is a scalar. The value '-1' signifies + void material. **/geometry/cells/cell /offset** (*int[]*) diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 2a8d07b8c0..8d25789f5a 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -6,8 +6,8 @@ element geometry { (element universe { xsd:int } | attribute universe { xsd:int })? & ( (element fill { xsd:int } | attribute fill { xsd:int }) | - (element material { ( xsd:int | "void" ) } | - attribute material { ( xsd:int | "void" ) }) + (element material { ( xsd:int | "void" )+ } | + attribute material { ( xsd:int | "void" )+ }) ) & (element region { xsd:string } | attribute region { xsd:string })? & (element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? & diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index fcb310d0b3..d40401b281 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -47,16 +47,20 @@ - - - void - + + + + void + + - - - void - + + + + void + + From 95057442914b2d36d87809dd23e55c724cd44fc6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 9 Jan 2016 23:44:26 -0500 Subject: [PATCH 08/17] Fix void distribmats --- openmc/universe.py | 13 +++++++++---- src/summary.F90 | 14 +++++++++++--- tests/test_distribmat/inputs_true.dat | 2 +- tests/test_distribmat/results_true.dat | 4 ++-- tests/test_distribmat/test_distribmat.py | 5 +---- 5 files changed, 24 insertions(+), 14 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 2fa535265c..048b6a0e9f 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -113,8 +113,11 @@ class Cell(object): string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', self._fill._id) elif isinstance(self._fill, Iterable): - string += ('{0: <16}{1}'.format('\tMaterial', '=\t') - + '[' + ', '.join([str(m.id) for m in self.fill]) + ']\n') + string += '{0: <16}{1}'.format('\tMaterial', '=\t') + string += '[' + string += ', '.join(['void' if m == 'void' else str(m.id) + for m in self.fill]) + string += ']\n' elif isinstance(self._fill, (Universe, Lattice)): string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill._id) @@ -209,7 +212,8 @@ class Cell(object): self._type = 'normal' elif isinstance(fill, Iterable): - cv.check_type('cell.fill', fill, Iterable, openmc.Material) + cv.check_type('cell.fill', fill, Iterable, + (openmc.Material, basestring)) self._type = 'normal' elif isinstance(fill, Universe): @@ -402,7 +406,8 @@ class Cell(object): element.set("material", str(self._fill._id)) elif isinstance(self._fill, Iterable): - element.set("material", ' '.join([str(m.id) for m in self.fill])) + element.set("material", ' '.join([m if m == 'void' else str(m.id) + for m in self.fill])) elif isinstance(self._fill, (Universe, Lattice)): element.set("fill", str(self._fill._id)) diff --git a/src/summary.F90 b/src/summary.F90 index deccaddf59..c06d2549c2 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -152,12 +152,20 @@ contains case (CELL_NORMAL) call write_dataset(cell_group, "fill_type", "normal") if (size(c % material) == 1) then - call write_dataset(cell_group, "material", & - materials(c % material(1)) % id) + if (c % material(1) == MATERIAL_VOID) then + call write_dataset(cell_group, "material", MATERIAL_VOID) + else + call write_dataset(cell_group, "material", & + materials(c % material(1)) % id) + end if else allocate(cell_materials(size(c % material))) do j = 1, size(c % material) - cell_materials(j) = materials(c % material(j)) % id + if (c % material(j) == MATERIAL_VOID) then + cell_materials(j) = MATERIAL_VOID + else + cell_materials(j) = materials(c % material(j)) % id + end if end do call write_dataset(cell_group, "material", cell_materials) deallocate(cell_materials) diff --git a/tests/test_distribmat/inputs_true.dat b/tests/test_distribmat/inputs_true.dat index 48044e2872..0fb8cc8bf8 100644 --- a/tests/test_distribmat/inputs_true.dat +++ b/tests/test_distribmat/inputs_true.dat @@ -1 +1 @@ -24a3537446feafcf79fc84d268b43130e9521ad7cccfcf733ba35ca50bf3c644e0b8e91bfdef7bdb4073783197cfdd730d1e459fc622aa3ce6dcf3143c805a1f \ No newline at end of file +bdc2acd3a4f5078c61d7cb83ff30e07b76c89b21d8131c1b0ce86a43156da7bd0b9ebb5d6acc0b98a7d8128667f99403ef12a6d76ed2ec6cb187810222e4f6b3 \ No newline at end of file diff --git a/tests/test_distribmat/results_true.dat b/tests/test_distribmat/results_true.dat index 5e93e4e473..70464fbc6c 100644 --- a/tests/test_distribmat/results_true.dat +++ b/tests/test_distribmat/results_true.dat @@ -1,9 +1,9 @@ k-combined: -1.433669E+00 7.069157E-03 +1.309285E+00 1.263629E-02 Cell ID = 11 Name = - Material = [2, 3, 2, 2] + Material = [2, 3, void, 2] Region = -10000 Rotation = None Translation = None diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py index 61e40e64c6..f0d09e963a 100644 --- a/tests/test_distribmat/test_distribmat.py +++ b/tests/test_distribmat/test_distribmat.py @@ -44,7 +44,7 @@ class DistribmatTestHarness(PyAPITestHarness): r0 = openmc.ZCylinder(R=0.3) c11 = openmc.Cell(cell_id=11) c11.region = -r0 - c11.fill = [dense_fuel, light_fuel] + [dense_fuel]*2 + c11.fill = [dense_fuel, light_fuel, 'void', dense_fuel] c12 = openmc.Cell(cell_id=12) c12.region = +r0 c12.fill = moderator @@ -122,9 +122,6 @@ class DistribmatTestHarness(PyAPITestHarness): outstr += str(su.get_cell_by_id(11)) return outstr -# def _cleanup(self): -# return None - if __name__ == '__main__': harness = DistribmatTestHarness('statepoint.5.*') From d0a0dda0b50cd7b13506e4e9b6663ea3a3c07de8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 10 Jan 2016 00:26:14 -0500 Subject: [PATCH 09/17] Check number of distributed materials --- src/initialize.F90 | 74 +++++++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 33 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 2c44755d20..5637585ae4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -946,7 +946,6 @@ contains subroutine prepare_distribcell() integer :: i, j ! Tally, filter loop counters - integer :: n_filt ! Number of filters originally in tally logical :: count_all ! Count all cells type(TallyObject), pointer :: t ! Current tally type(Universe), pointer :: univ ! Pointer to universe @@ -955,58 +954,67 @@ contains integer, allocatable :: counts(:,:) ! Target count logical, allocatable :: found(:,:) ! Target found + ! Do we need to count cell instances? count_all = .false. - ! Loop over tallies + ! We need to count instances if any distribcell filters are present. do i = 1, n_tallies - - ! Get pointer to tally - t => tallies(i) - - n_filt = t%n_filters - - ! Loop over the filters to determine how many additional filters - ! need to be added to this tally - do j = 1, t%n_filters - - ! Determine type of filter - if (t%filters(j)%type == FILTER_DISTRIBCELL) then + do j = 1, tallies(i) % n_filters + if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then count_all = .true. - if (size(t%filters(j)%int_bins) > 1) then + 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 do - end do + ! We also need to count instnaces if any distributed materials are present. + if (.not. count_all) then + do i = 1, n_cells + if (size(cells(i) % material) > 1) then + count_all = .true. + exit + end if + end do + end if + + ! Count the number of instances of each cell. if (count_all) then + call count_instance(universes(BASE_UNIVERSE)) + end if - univ => universes(BASE_UNIVERSE) - - ! sum the number of occurrences of all cells - call count_instance(univ) - - ! Loop over tallies + ! Set the number of bins in all distribcell filters. + if (count_all) then do i = 1, n_tallies - - ! Get pointer to tally t => tallies(i) - - ! Initialize the filters do j = 1, t%n_filters - - ! Set the number of bins to the number of instances of the cell - if (t%filters(j)%type == FILTER_DISTRIBCELL) then - c => cells(t%filters(j)%int_bins(1)) - t%filters(j)%n_bins = c%instances + if (t % filters(j) % type == FILTER_DISTRIBCELL) then + ! Set the number of bins to the number of instances of the cell + c => cells(t % filters(j) % int_bins(1)) + t % filters(j) % n_bins = c % instances end if - end do end do + end if + ! Make sure the number of materials matches the number of cell instances for + ! distributed materials + if (count_all) then + do i = 1, n_cells + associate (c => cells(i)) + if (size(c % material) > 1) then + if (size(c % material) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % material))) & + // " materials but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of materials must & + &equal one or the number of instances.") + end if + end if + end associate + end do end if ! Allocate offset maps at each level in the geometry From a5dc14a33bbf9953b496c7108ef8088e6870d549 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 11 Jan 2016 13:26:23 -0500 Subject: [PATCH 10/17] Replace 'count_all' with 'distribcell_active' --- docs/source/usersguide/output/summary.rst | 6 +- openmc/universe.py | 2 +- src/initialize.F90 | 73 +++++++++++------------ 3 files changed, 39 insertions(+), 42 deletions(-) diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index d100fab3bf..83602e5065 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -94,9 +94,9 @@ The current revision of the summary file format is 1. **/geometry/cells/cell /material** (*int* or *int[]*) Unique ID of the material(s) assigned to the cell. This dataset is present - only if fill_type is set to 'normal'. The data is an array if the cell uses - distributed materials, otherwise it is a scalar. The value '-1' signifies - void material. + only if fill_type is set to 'normal'. The value '-1' signifies void + material. The data is an array if the cell uses distributed materials, + otherwise it is a scalar. **/geometry/cells/cell /offset** (*int[]*) diff --git a/openmc/universe.py b/openmc/universe.py index 048b6a0e9f..ccd0cb27c3 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -50,7 +50,7 @@ class Cell(object): Unique identifier for the cell name : str Name of the cell - fill : Materials or Universe or Lattice or 'void' + fill : Material or Universe or Lattice or 'void' or iterable of Material Indicates what the region of space is filled with region : openmc.region.Region Region of space that is assigned to the cell. diff --git a/src/initialize.F90 b/src/initialize.F90 index 5637585ae4..cfbee78915 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -945,8 +945,8 @@ contains subroutine prepare_distribcell() - integer :: i, j ! Tally, filter loop counters - logical :: count_all ! Count all cells + integer :: i, j ! Tally, filter loop counters + logical :: distribcell_active ! Does simulation use distribcell? type(TallyObject), pointer :: t ! Current tally type(Universe), pointer :: univ ! Pointer to universe type(Cell), pointer :: c ! Pointer to cell @@ -954,14 +954,14 @@ contains integer, allocatable :: counts(:,:) ! Target count logical, allocatable :: found(:,:) ! Target found - ! Do we need to count cell instances? - count_all = .false. + ! Assume distribcell is not needed until proven otherwise. + distribcell_active = .false. - ! We need to count instances if any distribcell filters are present. + ! 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 - count_all = .true. + 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.") @@ -970,52 +970,49 @@ contains end do end do - ! We also need to count instnaces if any distributed materials are present. - if (.not. count_all) then + ! We also need distribcell if any distributed materials are present. + if (.not. distribcell_active) then do i = 1, n_cells if (size(cells(i) % material) > 1) then - count_all = .true. + distribcell_active = .true. exit end if end do end if + ! If distribcell isn't used in this simulation then no more work left to do. + if (.not. distribcell_active) return + ! Count the number of instances of each cell. - if (count_all) then - call count_instance(universes(BASE_UNIVERSE)) - end if + call count_instance(universes(BASE_UNIVERSE)) ! Set the number of bins in all distribcell filters. - if (count_all) then - do i = 1, n_tallies - t => tallies(i) - do j = 1, t%n_filters - if (t % filters(j) % type == FILTER_DISTRIBCELL) then - ! Set the number of bins to the number of instances of the cell - c => cells(t % filters(j) % int_bins(1)) - t % filters(j) % n_bins = c % instances - end if - end do + do i = 1, n_tallies + t => tallies(i) + do j = 1, t%n_filters + if (t % filters(j) % type == FILTER_DISTRIBCELL) then + ! Set the number of bins to the number of instances of the cell + c => cells(t % filters(j) % int_bins(1)) + t % filters(j) % n_bins = c % instances + end if end do - end if + end do ! Make sure the number of materials matches the number of cell instances for - ! distributed materials - if (count_all) then - do i = 1, n_cells - associate (c => cells(i)) - if (size(c % material) > 1) then - if (size(c % material) /= c % instances) then - call fatal_error("Cell " // trim(to_str(c % id)) // " was & - &specified with " // trim(to_str(size(c % material))) & - // " materials but has " // trim(to_str(c % instances)) & - // " distributed instances. The number of materials must & - &equal one or the number of instances.") - end if + ! distributed materials. + do i = 1, n_cells + associate (c => cells(i)) + if (size(c % material) > 1) then + if (size(c % material) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % material))) & + // " materials but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of materials must & + &equal one or the number of instances.") end if - end associate - end do - end if + end if + end associate + end do ! Allocate offset maps at each level in the geometry call allocate_offsets(univ_list, counts, found) From 2ea4deb2882893f8a573811b1c4efc23cb22c9a5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 11 Jan 2016 13:34:52 -0500 Subject: [PATCH 11/17] More prepare_distribcell clean-up --- src/initialize.F90 | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index cfbee78915..c5bfe22e7c 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -947,9 +947,6 @@ contains integer :: i, j ! Tally, filter loop counters logical :: distribcell_active ! Does simulation use distribcell? - type(TallyObject), pointer :: t ! Current tally - type(Universe), pointer :: univ ! Pointer to universe - type(Cell), pointer :: c ! Pointer to cell integer, allocatable :: univ_list(:) ! Target offsets integer, allocatable :: counts(:,:) ! Target count logical, allocatable :: found(:,:) ! Target found @@ -988,13 +985,13 @@ contains ! Set the number of bins in all distribcell filters. do i = 1, n_tallies - t => tallies(i) - do j = 1, t%n_filters - if (t % filters(j) % type == FILTER_DISTRIBCELL) then - ! Set the number of bins to the number of instances of the cell - c => cells(t % filters(j) % int_bins(1)) - t % filters(j) % n_bins = c % instances - end if + 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 end do end do @@ -1020,16 +1017,10 @@ contains ! Calculate offsets for each target distribcell do i = 1, n_maps do j = 1, n_universes - univ => universes(j) - call calc_offsets(univ_list(i), i, univ, counts, found) + call calc_offsets(univ_list(i), i, universes(j), counts, found) end do end do - ! Deallocate temporary target variable arrays - deallocate(counts) - deallocate(found) - deallocate(univ_list) - end subroutine prepare_distribcell !=============================================================================== From f56f613593e7f0dee3ab80b609501580ea9ed7b4 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 11 Jan 2016 13:57:18 -0500 Subject: [PATCH 12/17] Clean-up allocate offsets --- src/initialize.F90 | 70 ++++++++++++++++++---------------------------- 1 file changed, 27 insertions(+), 43 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index c5bfe22e7c..b7e3417856 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1034,47 +1034,33 @@ contains integer, intent(out), allocatable :: counts(:,:) ! Target count logical, intent(out), allocatable :: found(:,:) ! Target found - integer :: i, j, k ! Loop counters - type(SetInt) :: cell_list ! distribells to track - type(Universe), pointer :: univ ! pointer to universe - class(Lattice), pointer :: lat ! pointer to lattice - type(TallyObject), pointer :: t ! pointer to tally - type(TallyFilter), pointer :: filter ! pointer to filter - type(Cell), pointer :: c ! pointer to cell + integer :: i, j, k ! Loop counters + type(SetInt) :: cell_list ! distribells to track ! Begin gathering list of cells in distribcell tallies n_maps = 0 ! List all cells referenced in distribcell filters. do i = 1, n_tallies - t => tallies(i) - - do j = 1, t % n_filters - filter => t % filters(j) - - if (filter % type == FILTER_DISTRIBCELL) then - if (.not. cell_list % contains(filter % int_bins(1))) then - call cell_list % add(filter % int_bins(1)) - end if + 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 - end do end do ! List all cells with multiple (distributed) materials. do i = 1, n_cells - c => cells(i) - if (size(c % material) > 1) then - if (.not. cell_list % contains(i)) call cell_list % add(i) + if (size(cells(i) % material) > 1) then + call cell_list % add(i) end if end do ! Compute the number of unique universes containing these distribcells ! to determine the number of offset tables to allocate do i = 1, n_universes - univ => universes(i) - do j = 1, univ % n_cells - if (cell_list % contains(univ % cells(j))) then + do j = 1, universes(i) % n_cells + if (cell_list % contains(universes(i) % cells(j))) then n_maps = n_maps + 1 end if end do @@ -1085,22 +1071,21 @@ contains ! Allocate list to accumulate target distribcell counts in each universe allocate(counts(n_universes, n_maps)) + counts(:,:) = 0 ! Allocate list to track if target distribcells are found in each universe allocate(found(n_universes, n_maps)) - - counts(:,:) = 0 found(:,:) = .false. - k = 1 + ! Search through universes for distributed cells and assign each one a ! unique distribcell array index. + k = 1 do i = 1, n_universes - univ => universes(i) - do j = 1, univ % n_cells - if (cell_list % contains(univ % cells(j))) then - cells(univ % cells(j)) % distribcell_index = k - univ_list(k) = univ % id + do j = 1, universes(i) % n_cells + if (cell_list % contains(universes(i) % cells(j))) then + cells(universes(i) % cells(j)) % distribcell_index = k + univ_list(k) = universes(i) % id k = k + 1 end if end do @@ -1108,20 +1093,19 @@ contains ! Allocate the offset tables for lattices do i = 1, n_lattices - lat => lattices(i) % obj + associate(lat => lattices(i) % obj) + select type(lat) - select type(lat) - - type is (RectLattice) - allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & - lat % n_cells(3))) - type is (HexLattice) - allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & - 2 * lat % n_rings - 1, lat % n_axial)) - end select - - lat % offset(:, :, :, :) = 0 + type is (RectLattice) + allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & + lat % n_cells(3))) + type is (HexLattice) + allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & + 2 * lat % n_rings - 1, lat % n_axial)) + end select + lat % offset(:, :, :, :) = 0 + end associate end do ! Allocate offset table for fill cells From a7ac2db068de2d81a5831c4817470170ff7e4efe Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 12 Jan 2016 13:31:18 -0500 Subject: [PATCH 13/17] Typo fix --- src/initialize.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index b7e3417856..a50fe09a01 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -959,7 +959,7 @@ contains do j = 1, tallies(i) % n_filters if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then distribcell_active = .true. - if (size(tallies(i) % filters(j) %int_bins) > 1) then + 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 From c37639d90dacaaa89ca33385ea9ed95e966bbc77 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 12 Jan 2016 16:04:50 -0500 Subject: [PATCH 14/17] Make max lost particles relative to total # part.s --- src/constants.F90 | 3 +++ src/geometry.F90 | 8 +++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index 8eabdaf604..10b41e5602 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -139,6 +139,9 @@ module constants ! Maximum number of lost particles integer, parameter :: MAX_LOST_PARTICLES = 10 + ! Maximum number of lost particles, relative to the total number of particles + real(8), parameter :: REL_MAX_LOST_PARTICLES = 1e-5_8 + ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS diff --git a/src/geometry.F90 b/src/geometry.F90 index 9a084a77cb..667dc02bfa 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -947,6 +947,8 @@ contains type(Particle), intent(inout) :: p character(*) :: message + integer(8) :: tot_n_particles + ! Print warning and write lost particle file call warning(message) call write_particle_restart(p) @@ -956,9 +958,13 @@ contains !$omp atomic n_lost_particles = n_lost_particles + 1 + ! Count the total number of simulated particles + tot_n_particles = n_batches * gen_per_batch * n_particles + ! Abort the simulation if the maximum number of lost particles has been ! reached - if (n_lost_particles == MAX_LOST_PARTICLES) then + if (n_lost_particles >= MAX_LOST_PARTICLES .and. & + n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then call fatal_error("Maximum number of lost particles has been reached.") end if From 445c96980445703f3c8561c499eb788a092c51f2 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 13 Jan 2016 11:36:05 -0500 Subject: [PATCH 15/17] Small fixes for #548 --- src/initialize.F90 | 7 ++- src/input_xml.F90 | 4 +- tests/test_distribmat/test_distribmat.py | 56 +++++++++++++----------- 3 files changed, 38 insertions(+), 29 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index a50fe09a01..de681fb949 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -625,7 +625,7 @@ contains c % fill = lid else call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "& - &// trim(to_str(c % id)) // " is neither a universe nor a & + // trim(to_str(c % id)) // " is neither a universe nor a & &lattice.") end if else @@ -638,7 +638,7 @@ contains c % material(j) = material_dict % get_key(id) else call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c % id))) + // " specified on cell " // trim(to_str(c % id))) end if end do end if @@ -1115,6 +1115,9 @@ contains end if end do + ! Free up memory + call cell_list % clear() + end subroutine allocate_offsets end module initialize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 41ee5aab2a..1ca8cc8799 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1104,7 +1104,7 @@ contains ! Check for error if (c % material(j) == ERROR_INT) then call fatal_error("Invalid material specified on cell " & - &// to_str(c % id)) + // to_str(c % id)) end if end select end do @@ -1124,7 +1124,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material(1) == NONE .and. c % fill == NONE) then call fatal_error("Neither material nor fill was specified for cell " & - &// trim(to_str(c % id))) + // trim(to_str(c % id))) end if ! Check to make sure that both material and fill haven't been diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py index f0d09e963a..b1252a8e7b 100644 --- a/tests/test_distribmat/test_distribmat.py +++ b/tests/test_distribmat/test_distribmat.py @@ -90,31 +90,31 @@ class DistribmatTestHarness(PyAPITestHarness): sets_file.export_to_xml() -# #################### -# # Plots -# #################### -# -# plots_file = openmc.PlotsFile() -# -# plot = openmc.Plot(plot_id=1) -# plot.basis = 'xy' -# plot.color = 'cell' -# plot.filename = 'cellplot' -# plot.origin = (0, 0, 0) -# plot.width = (7, 7) -# plot.pixels = (400, 400) -# plots_file.add_plot(plot) -# -# plot = openmc.Plot(plot_id=2) -# plot.basis = 'xy' -# plot.color = 'mat' -# plot.filename = 'matplot' -# plot.origin = (0, 0, 0) -# plot.width = (7, 7) -# plot.pixels = (400, 400) -# plots_file.add_plot(plot) -# -# plots_file.export_to_xml() + #################### + # Plots + #################### + + plots_file = openmc.PlotsFile() + + plot = openmc.Plot(plot_id=1) + plot.basis = 'xy' + plot.color = 'cell' + plot.filename = 'cellplot' + plot.origin = (0, 0, 0) + plot.width = (7, 7) + plot.pixels = (400, 400) + plots_file.add_plot(plot) + + plot = openmc.Plot(plot_id=2) + plot.basis = 'xy' + plot.color = 'mat' + plot.filename = 'matplot' + plot.origin = (0, 0, 0) + plot.width = (7, 7) + plot.pixels = (400, 400) + plots_file.add_plot(plot) + + plots_file.export_to_xml() def _get_results(self): outstr = super(DistribmatTestHarness, self)._get_results() @@ -122,6 +122,12 @@ class DistribmatTestHarness(PyAPITestHarness): outstr += str(su.get_cell_by_id(11)) return outstr + def _cleanup(self): + f = os.path.join(os.getcwd(), 'plots.xml') + if os.path.exists(f): + os.remove(f) + super(DistribmatTestHarness, self)._cleanup() + if __name__ == '__main__': harness = DistribmatTestHarness('statepoint.5.*') From d5400cb1dafd51dd09fcd0369d30d3730ead4b60 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 13 Jan 2016 11:45:44 -0500 Subject: [PATCH 16/17] Small fixes for #555 --- src/constants.F90 | 2 +- src/geometry.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 10b41e5602..e5585fcb92 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -140,7 +140,7 @@ module constants integer, parameter :: MAX_LOST_PARTICLES = 10 ! Maximum number of lost particles, relative to the total number of particles - real(8), parameter :: REL_MAX_LOST_PARTICLES = 1e-5_8 + real(8), parameter :: REL_MAX_LOST_PARTICLES = 1e-6_8 ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS diff --git a/src/geometry.F90 b/src/geometry.F90 index 667dc02bfa..ffa07c2978 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -958,8 +958,8 @@ contains !$omp atomic n_lost_particles = n_lost_particles + 1 - ! Count the total number of simulated particles - tot_n_particles = n_batches * gen_per_batch * n_particles + ! Count the total number of simulated particles (on this processor) + tot_n_particles = n_batches * gen_per_batch * work ! Abort the simulation if the maximum number of lost particles has been ! reached From c5eb4b79cc1a5ba28317532148ef71e954ca6ebc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 13 Jan 2016 12:24:25 -0500 Subject: [PATCH 17/17] Update distribmat test --- tests/test_distribmat/inputs_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_distribmat/inputs_true.dat b/tests/test_distribmat/inputs_true.dat index 0fb8cc8bf8..9d0a70c471 100644 --- a/tests/test_distribmat/inputs_true.dat +++ b/tests/test_distribmat/inputs_true.dat @@ -1 +1 @@ -bdc2acd3a4f5078c61d7cb83ff30e07b76c89b21d8131c1b0ce86a43156da7bd0b9ebb5d6acc0b98a7d8128667f99403ef12a6d76ed2ec6cb187810222e4f6b3 \ No newline at end of file +83a9f1412b1ddfdd8b9fa8c7e8b44be8137dc6aa785d44c3eea9f0928242475aea4ba73ba4eef01afdf63bd1bff18116e40b798eeafab50c856186de4c68d5a6 \ No newline at end of file