diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 30e9ed07b9..0f9a3be18a 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -987,6 +987,15 @@ Each ```` element can have the following attributes or sub-elements: *Default*: A region filling all space. + :temperature: + The temperature of the cell in Kelvin. If windowed-multipole data is + avalable, this temperature will be used to Doppler broaden some cross + sections in the resolved resonance region. A list of temperatures can be + specified for the "distributed temperature" feature. This will give each + unique instance of the cell its own temperature. + + *Default*: The temperature of the coldest nuclide in the cell's material(s) + :rotation: If the cell is filled with a universe, this element specifies the angles in degrees about the x, y, and z axes that the filled universe should be diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index 83602e5065..3be0cfdfe9 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -98,6 +98,10 @@ The current revision of the summary file format is 1. material. The data is an array if the cell uses distributed materials, otherwise it is a scalar. +**/geometry/cells/cell /temperature** (*double[]*) + + Temperature of the cell in Kelvin. + **/geometry/cells/cell /offset** (*int[]*) Offsets used for distribcell tally filter. This dataset is present only if diff --git a/openmc/summary.py b/openmc/summary.py index eb14d3bb81..d8f8978306 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -264,6 +264,9 @@ class Summary(object): self._f['geometry/cells'][key]['rotation'][...] rotation = np.asarray(rotation, dtype=np.int) cell.rotation = rotation + elif fill_type == 'normal': + cell.temperature = \ + self._f['geometry/cells'][key]['temperature'][...] # Store Cell fill information for after Universe/Lattice creation self._cell_fills[index] = (fill_type, fill) diff --git a/openmc/universe.py b/openmc/universe.py index 237ddce1f7..027d43fc77 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -51,9 +51,13 @@ class Cell(object): name : str Name of the cell fill : Material or Universe or Lattice or 'void' or iterable of Material - Indicates what the region of space is filled with + Indicates what the region of space is filled with. Multiple materials + can be given to give each distributed cell instance a unique material. region : openmc.region.Region Region of space that is assigned to the cell. + temperature : float or iterable of float + Temperature of the cell in Kelvin. Multiple temperatures can be given + to give each distributed cell instance a unique temperature. rotation : ndarray If the cell is filled with a universe, this array specifies the angles in degrees about the x, y, and z axes that the filled universe should be @@ -75,6 +79,7 @@ class Cell(object): self._fill = None self._type = None self._region = None + self._temperature = None self._rotation = None self._translation = None self._offsets = None @@ -91,6 +96,8 @@ class Cell(object): return False elif self.region != other.region: return False + elif self.temperature != other.temperature: + return False elif self.rotation != other.rotation: return False elif self.translation != other.translation: @@ -126,6 +133,10 @@ class Cell(object): string += '{0: <16}{1}{2}\n'.format('\tRegion', '=\t', self._region) + if self.fill_type == 'material': + string += '\t{0: <15}=\t{1}\n'.format('Temperature', + self.temperature) + string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t', self._rotation) string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t', @@ -150,7 +161,7 @@ class Cell(object): @property def fill_type(self): - if isinstance(self.fill, openmc.Material): + if isinstance(self.fill, (openmc.Material, Iterable)): return 'material' elif isinstance(self.fill, openmc.Universe): return 'universe' @@ -163,6 +174,10 @@ class Cell(object): def region(self): return self._region + @property + def temperature(self): + return self._temperature + @property def rotation(self): return self._rotation @@ -251,6 +266,17 @@ class Cell(object): cv.check_type('cell region', region, Region) self._region = region + @temperature.setter + def temperature(self, temperature): + cv.check_type('cell temperature', temperature, (Iterable, Real)) + if isinstance(temperature, Iterable): + cv.check_type('cell temperature', temperature, Iterable, Real) + for T in temperature: + cv.check_greater_than('cell temperature', T, 0.0, True) + else: + cv.check_greater_than('cell temperature', temperature, 0.0, True) + self._temperature = temperature + @distribcell_index.setter def distribcell_index(self, ind): cv.check_type('distribcell index', ind, Integral) @@ -445,6 +471,13 @@ class Cell(object): # Call the recursive function from the top node create_surface_elements(self.region, xml_element) + if self.temperature is not None: + if isinstance(self.temperature, Iterable): + element.set("temperature", ' '.join( + [str(t) for t in self.temperature])) + else: + element.set("temperature", str(self.temperature)) + if self.translation is not None: element.set("translation", ' '.join(map(str, self.translation))) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b62d9b389c..033e80b174 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -105,10 +105,13 @@ contains i_nuclide = mat % nuclide(i) ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(i_nuclide) % last_E .or. mat % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT) + if (p % E /= micro_xs(i_nuclide) % last_E & + .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & + i_grid, p % sqrtkT) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & + i_grid, p % sqrtkT) end if ! ======================================================================== @@ -145,7 +148,8 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union, sqrtkT) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, & + i_log_union, sqrtkT) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy @@ -590,14 +594,20 @@ contains ! sections in the resolved resonance regions !=============================================================================== - subroutine multipole_eval(multipole, Emev, sqrtkT, sigT, sigA, sigF) - type(MultipoleArray), intent(in) :: multipole ! The windowed multipole object to process. - real(8), intent(in) :: Emev ! The energy at which to evaluate the cross section in MeV - real(8), intent(in) :: sqrtkT ! The temperature in the form sqrt(kT (in eV)), at which to evaluate the cross section. - real(8), intent(out) :: sigT ! Total cross section - real(8), intent(out) :: sigA ! Absorption cross section - real(8), intent(out) :: sigF ! Fission cross section - complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic form + subroutine multipole_eval(multipole, Emev, sqrtkT_, sigT, sigA, sigF) + type(MultipoleArray), intent(in) :: multipole ! The windowed multipole + ! object to process. + real(8), intent(in) :: Emev ! The energy at which to + ! evaluate the cross section + ! in MeV + real(8), intent(in) :: sqrtkT_ ! The temperature in the form + ! sqrt(kT (in MeV)), at which + ! to evaluate the XS. + real(8), intent(out) :: sigT ! Total cross section + real(8), intent(out) :: sigA ! Absorption cross section + real(8), intent(out) :: sigF ! Fission cross section + complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic + ! form complex(8) :: c_temp ! complex temporary variable complex(8) :: w_val ! The faddeeva function evaluated at Z complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole) @@ -607,6 +617,7 @@ contains real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E real(8) :: temp ! real temporary value real(8) :: E ! energy, eV + real(8) :: sqrtkT ! sqrt(kT (in eV)) integer :: iP ! index of pole integer :: iC ! index of curvefit integer :: iW ! index of window @@ -617,6 +628,7 @@ contains ! Convert to eV E = Emev * 1.0e6_8 + sqrtkT = sqrtkT_ * 1.0e3_8 sqrtE = sqrt(E) invE = ONE/E diff --git a/src/geometry.F90 b/src/geometry.F90 index 8a38f982b5..5f37f1047b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -279,6 +279,36 @@ contains p % material = c % material(offset + 1) end if + ! Set the particle temperature + if (size(c % sqrtkT) == 1) then + ! Only one temperature for this cell; assign that one to the particle. + p % sqrtkT = c % sqrtkT(1) + else + ! Distributed instances of this cell have different temperatures. + ! Determine which instance this is and assign the matching temp. + 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 % sqrtkT = c % sqrtkT(offset + 1) + end if + elseif (c % type == CELL_FILL) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 519b15bff4..0756e23dc4 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -139,6 +139,9 @@ module geometry_header ! only) integer :: distribcell_index ! Index corresponding to this cell in ! distribcell arrays + real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann * + ! temperature in MeV. Multiple for + ! distribcell ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index de681fb949..610bfe538b 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -120,6 +120,9 @@ contains ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() + ! Set undefined cell temperatures to match the material data. + call lookup_material_temperatures() + ! Construct unionized or log energy grid for cross-sections select case (grid_method) case (GRID_NUCLIDE) @@ -967,10 +970,11 @@ contains end do end do - ! We also need distribcell if any distributed materials are present. + ! We also need distribcell if any distributed materials or distributed + ! temperatues are present. if (.not. distribcell_active) then do i = 1, n_cells - if (size(cells(i) % material) > 1) then + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then distribcell_active = .true. exit end if @@ -995,8 +999,8 @@ contains end do end do - ! Make sure the number of materials matches the number of cell instances for - ! distributed materials. + ! Make sure the number of materials and temperatures matches the number of + ! cell instances. do i = 1, n_cells associate (c => cells(i)) if (size(c % material) > 1) then @@ -1008,6 +1012,15 @@ contains &equal one or the number of instances.") end if end if + if (size(c % sqrtkT) > 1) then + if (size(c % sqrtkT) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % sqrtkT))) & + // " temperatures but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of temperatures must & + &equal one or the number of instances.") + end if + end if end associate end do @@ -1049,9 +1062,9 @@ contains end do end do - ! List all cells with multiple (distributed) materials. + ! List all cells with multiple (distributed) materials or temperatures. do i = 1, n_cells - if (size(cells(i) % material) > 1) then + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then call cell_list % add(i) end if end do @@ -1120,4 +1133,57 @@ contains end subroutine allocate_offsets +!=============================================================================== +! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to +! find their temperatures from material data. +!=============================================================================== + + subroutine lookup_material_temperatures() + integer :: i, j, k + real(8) :: min_temp + logical :: warning_given + + warning_given = .false. + do i = 1, n_cells + ! Ignore non-normal cells and cells with defined temperature. + if (cells(i) % type /= CELL_NORMAL) cycle + if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle + + ! Set the number of temperatures equal to the number of materials. + deallocate(cells(i) % sqrtkT) + allocate(cells(i) % sqrtkT(size(cells(i) % material))) + + ! Check each of the cell materials for temperature data. + do j = 1, size(cells(i) % material) + ! Arbitrarily set void regions to 0K. + if (cells(i) % material(j) == MATERIAL_VOID) then + cells(i) % sqrtkT(j) = ZERO + cycle + end if + + associate (mat => materials(cells(i) % material(j))) + ! Find the temperature of the coldest nuclide. + min_temp = nuclides(mat % nuclide(1)) % kT + do k = 2, mat % n_nuclides + ! Warn the user if the nuclides don't have identical temperatues. + if (nuclides(mat % nuclide(k)) % kT /= min_temp & + .and. .not. warning_given) then + call warning("OpenMC cannot & + &identify the temperature of at least one cell. For the & + &purposes of multipole cross section evaluations, all cells & + &with unknown temperature will be set to the coldest & + &temperature found in the nuclear data for that cell's & + &material") + warning_given = .true. + end if + min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT) + end do + + ! Set the temperature for this cell instance. + cells(i) % sqrtkT(j) = sqrt(min_temp) + end associate + end do + end do + end subroutine lookup_material_temperatures + end module initialize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bbf4da4f9e..34cd0f1b95 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1297,6 +1297,40 @@ contains call get_node_array(node_cell, "translation", c % translation) end if + ! Read cell temperatures. If the temperature is not specified, set it to + ! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with + ! the temperature from the material data. + if (check_for_node(node_cell, "temperature")) then + n = get_arraysize_double(node_cell, "temperature") + if (n > 0) then + ! Make sure this is a "normal" cell. + if (c % material(1) == NONE) call fatal_error("Cell " & + // trim(to_str(c % id)) // " was specified with a temperature & + &but no material. Temperature specification is only valid for & + &cells filled with a material.") + + ! Copy in temperatures + allocate(c % sqrtkT(n)) + call get_node_array(node_cell, "temperature", c % sqrtkT) + + ! Make sure all temperatues are positive + do j = 1, size(c % sqrtkT) + if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " & + // trim(to_str(c % id)) // " was specified with a negative & + &temperature. All cell temperatures must be non-negative.") + end do + + ! Convert to sqrt(kT) + c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:)) + else + allocate(c % sqrtkT(1)) + c % sqrtkT(1) = ERROR_REAL + end if + else + allocate(c % sqrtkT(1)) + c % sqrtkT = ERROR_REAL + end if + ! Add cell to dictionary call cell_dict % add_key(c % id, i) @@ -1854,7 +1888,6 @@ contains real(8) :: temp_dble ! temporary double prec. real logical :: file_exists ! does materials.xml exist? logical :: sum_density ! density is taken to be sum of nuclide densities - logical :: temp_known ! Has the temperature yet been defined? character(12) :: name ! name of isotope, e.g. 92235.03c character(12) :: alias ! alias of nuclide, e.g. U-235.03c character(MAX_WORD_LEN) :: units ! units on density @@ -1938,17 +1971,6 @@ contains cycle end if - ! ======================================================================= - ! READ AND PARSE TAG - if (check_for_node(node_mat, "temperature")) then - call get_node_ptr(node_mat, "temperature", node_temp) - call get_node_value(node_temp, "value", temp_dble) - mat % sqrtkT = sqrt(temp_dble * K_BOLTZMANN * 1.0D6) - temp_known = .true. - else - temp_known = .false. - end if - ! ======================================================================= ! READ AND PARSE TAG @@ -2055,16 +2077,6 @@ contains call get_node_value(node_nuc, "xs", name) name = trim(temp_str) // "." // trim(name) - ! If needed, look up temperature - if (.not. temp_known) then - ! Find xs_listing and set the name/alias according to the listing - index_list = xs_listing_dict % get_key(to_lower(name)) - if(xs_listings(index_list) % kT /= 0.0_8) then - mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6) - temp_known = .true. - end if - end if - ! save name and density to list call list_names % append(name) @@ -2153,16 +2165,6 @@ contains temp_str = "data" end if - ! If still needed, look up temperature - if (.not. temp_known) then - ! Find xs_listing and set kT - index_list = xs_listing_dict % get_key(to_lower(list_names % tail % data)) - if(xs_listings(index_list) % kT /= 0.0_8) then - mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6) - temp_known = .true. - end if - end if - ! Set ace or iso-in-lab scattering for each nuclide in element do k = 1, n_nuc_ele if (adjustl(to_lower(temp_str)) == "iso-in-lab") then @@ -2177,12 +2179,6 @@ contains end do NATURAL_ELEMENTS - ! If still undefined, set the temperature to zero - if (.not. temp_known) then - mat % sqrtkT = 0.0_8 - temp_known = .true. - end if - ! ======================================================================== ! COPY NUCLIDES TO ARRAYS IN MATERIAL diff --git a/src/material_header.F90 b/src/material_header.F90 index 397fa0cd96..91c4cdfb82 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,7 +13,6 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - real(8) :: sqrtkT ! sqrt(kT), kT in eV ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 0b9b251ee5..f1d7375d7a 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -2,7 +2,7 @@ module particle_header use bank_header, only: Bank use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, & - MAX_DELAYED_GROUPS + MAX_DELAYED_GROUPS, ERROR_REAL use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE @@ -79,6 +79,9 @@ module particle_header integer :: material ! index for current material integer :: last_material ! index for last material + ! Temperature of the current cell + real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV + ! Statistical data integer :: n_collision ! # of collisions @@ -124,6 +127,7 @@ contains this % absorb_wgt = ZERO this % n_bank = 0 this % wgt_bank = ZERO + this % sqrtkT = ERROR_REAL this % n_collision = 0 this % fission = .false. this % delayed_group = 0 diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 8d25789f5a..1bca009a61 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -9,6 +9,8 @@ element geometry { (element material { ( xsd:int | "void" )+ } | attribute material { ( xsd:int | "void" )+ }) ) & + (element temperature { list { xsd:double+ } } | + attribute temperature { list { xsd:double+ } } )? & (element region { xsd:string } | attribute region { xsd:string })? & (element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? & (element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })? diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index d40401b281..ba8eab73d3 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -64,6 +64,24 @@ + + + + + + + + + + + + + + + + + + diff --git a/src/summary.F90 b/src/summary.F90 index c06d2549c2..2d8001e845 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -108,6 +108,7 @@ contains integer :: i, j, k, m integer, allocatable :: lattice_universes(:,:,:) integer, allocatable :: cell_materials(:) + real(8), allocatable :: cell_temperatures(:) integer(HID_T) :: geom_group integer(HID_T) :: cells_group, cell_group integer(HID_T) :: surfaces_group, surface_group @@ -151,6 +152,7 @@ contains select case (c%type) case (CELL_NORMAL) call write_dataset(cell_group, "fill_type", "normal") + if (size(c % material) == 1) then if (c % material(1) == MATERIAL_VOID) then call write_dataset(cell_group, "material", MATERIAL_VOID) @@ -171,6 +173,12 @@ contains deallocate(cell_materials) end if + allocate(cell_temperatures(size(c % sqrtkT))) + cell_temperatures(:) = c % sqrtkT(:) + cell_temperatures(:) = cell_temperatures(:)**2 / K_BOLTZMANN + call write_dataset(cell_group, "temperature", cell_temperatures) + deallocate(cell_temperatures) + case (CELL_FILL) call write_dataset(cell_group, "fill_type", "universe") call write_dataset(cell_group, "fill", universes(c%fill)%id) diff --git a/tests/test_distribmat/results_true.dat b/tests/test_distribmat/results_true.dat index 70464fbc6c..73a48907cb 100644 --- a/tests/test_distribmat/results_true.dat +++ b/tests/test_distribmat/results_true.dat @@ -5,6 +5,7 @@ Cell Name = Material = [2, 3, void, 2] Region = -10000 + Temperature = [ 293.60594237 293.60594237 0. 293.60594237] Rotation = None Translation = None Offset = None diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat new file mode 100644 index 0000000000..a429c9b8d0 --- /dev/null +++ b/tests/test_multipole/inputs_true.dat @@ -0,0 +1 @@ +5c1cec635da5c4c869bdf58f62924a4cb1648e4ddecf0ce71b823cf45178767ba7f2e089b76d36e8616b1b21ffa43b32ab1b17d20bb74120f900b9e3e9ab9bcc \ No newline at end of file diff --git a/tests/test_multipole/results_true.dat b/tests/test_multipole/results_true.dat new file mode 100644 index 0000000000..19dd4d26b7 --- /dev/null +++ b/tests/test_multipole/results_true.dat @@ -0,0 +1,12 @@ +k-combined: +1.445285E+00 9.521660E-03 +Cell + ID = 11 + Name = + Material = 2 + Region = -10000 + Temperature = [ 500. 0. 700. 800.] + Rotation = None + Translation = None + Offset = None + Distribcell index= 1 diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py new file mode 100644 index 0000000000..a5437e433a --- /dev/null +++ b/tests/test_multipole/test_multipole.py @@ -0,0 +1,133 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness, PyAPITestHarness +import openmc +from openmc.stats import Box +from openmc.source import Source + + +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) + + mats_file = openmc.MaterialsFile() + mats_file.default_xs = '71c' + mats_file.add_materials([moderator, dense_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 + c11.temperature = [500, 0, 700, 800] + 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.source = 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() + + 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): + 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.*') + harness.main()