diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index a4d0f7d24..b9043a944 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -66,6 +66,8 @@ Other Methods A good survey of other energy grid techniques, including unionized energy grids, can be found in a paper by Leppanen_. +.. _windowed_multipole: + Windowed Multipole Representation --------------------------------- @@ -141,6 +143,40 @@ but not always the case. Future library versions may eliminate this issue. The data format used by OpenMC to represent windowed multipole data is specified in :ref:`io_data_wmp`. +.. _temperature_treatment: + +Temperature Treatment +--------------------- + +At the beginning of a simulation, OpenMC collects a list of all temperatures +that are present in a model. It then uses this list to determine what cross +sections to load. The data that is loaded depends on what temperature method has +been selected. There are three methods available: + +:Nearest: Cross sections are loaded only if they are within a specified + tolerance of the actual temperatures in the model. + +:Interpolation: Cross sections are loaded at temperatures that bound the actual + temperatures in the model. During transport, cross sections for + each material are calculated using statistical linear-linear + interpolation between bounding temperature. Suppose cross + sections are available at temperatures :math:`T_1, T_2, ..., + T_n` and a material is assigned a temperature :math:`T` where + :math:`T_i < T < T_{i+1}`. Statistical interpolation is applied + as follows: a uniformly-distributed random number of the unit + interval, :math:`\xi`, is sampled. If :math:`\xi < (T - + T_i)/(T_{i+1} - T_i)`, then cross sections at temperature + :math:`T_{i+1}` are used. Otherwise, cross sections at + :math:`T_i` are used. This procedure is applied for pointwise + cross sections in the resolved resonance range, unresolved + resonance probability tables, and :math:`S(\alpha,\beta)` + thermal scattering tables. + +:Multipole: Resolved resonance cross sections are calculated on-the-fly using + techniques/data described in :ref:`windowed_multipole`. Cross + section data is loaded for a single temperature and is used in the + unresolved resonance and fast energy ranges. + ---------------- Multi-Group Data ---------------- diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index c689883d0..c2e7270b6 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -265,6 +265,8 @@ Again, we need to check whether the denominator is zero. If so, this means that the particle's direction of flight is parallel to the plane and it will therefore never hit the plane. +.. _cylinder_distance: + Cylinder Parallel to an Axis ---------------------------- @@ -366,7 +368,74 @@ will then be either both positive or both negative. If they are both positive, the smaller (closer) one will be the solution with a negative sign on the square root of the discriminant. -.. TODO: Need to add derivation for x-cone, y-cone, and z-cone. +Cone Parallel to an Axis +------------------------ + +The equation for a cone parallel to, for example, the x-axis is :math:`(y - +y_0)^2 + (z - z_0)^2 = R^2(x - x_0)^2`. Thus, we need to solve :math:`(y + dv - +y_0)^2 + (z + dw - z_0)^2 = R^2(x + du - x_0)^2`. Let us define :math:`\bar{x} = +x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. We then have + +.. math:: + :label: dist-xcone-1 + + (\bar{y} + dv)^2 + (\bar{z} + dw)^2 = R^2(\bar{x} + du)^2 + +Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain + +.. math:: + :label: dist-xcylinder-2 + + (v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d + + (\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0 + +Defining the terms + +.. math:: + :label: dist-quadric-terms + + a = v^2 + w^2 - R^2u^2 + + k = \bar{y}v + \bar{z}w - R^2\bar{x}u + + c = \bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2 + +we then have the simple quadratic equation :math:`ad^2 + 2kd + c = 0` which can +be solved as described in :ref:`cylinder_distance`. + +General Quadric +--------------- + +The equation for a general quadric surface is :math:`Ax^2 + By^2 + Cz^2 + Dxy + +Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, we need to solve the equation + +.. math:: + :label: dist-quadric-1 + + A(x+du)^2 + B(y+dv)^2 + C(z+dw)^2 + D(x+du)(y+dv) + E(y+dv)(z+dw) + \\ + F(x+du)(z+dw) + G(x+du) + H(y+dv) + J(z+dw) + K = 0 + +Expanding equation :eq:`dist-quadric-1` and rearranging terms, we obtain + +.. math:: + :label: dist-quadric-2 + + d^2(uv + vw + uw) + 2d(Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + \\ + F(wx + uz))/2) + (x(Ax + Dy) + y(By + Ez) + z(Cz + Fx)) = 0 + +Defining the terms + +.. math:: + :label: dist-quadric-terms + + a = uv + vw + uw + + k = Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + F(wx + uz))/2 + + c = x(Ax + Dy) + y(By + Ez) + z(Cz + Fx) + +we then have the simple quadratic equation :math:`ad^2 + 2kd + c = 0` which can +be solved as described in :ref:`cylinder_distance`. .. _find-cell: @@ -810,6 +879,18 @@ form of the solution: w' = w + \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w)}{R^2 (1 + R^2) \bar{z}} +General Quadric +--------------- + +A general quadric surface has the form :math:`f(x,y,z) = Ax^2 + By^2 + Cz^2 + +Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is + +.. math:: + :label: reflection-quadric-grad + + \nabla f = \left ( \begin{array}{c} 2Ax + Dy + Fz + G \\ 2By + Dx + Ez + H + \\ 2Cz + Ey + Fx + J \end{array} \right ). + .. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry .. _surfaces: http://en.wikipedia.org/wiki/Surface diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 66300dc52..c5b29b194 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -53,6 +53,16 @@ Benchmarking Coupling and Multi-physics -------------------------- +- Matthew Ellis, Benoit Forget, Kord Smith, and Derek Gaston, "Continuous + Temperature Representation in Coupled OpenMC/MOOSE Simulations," *Proc. PHYSOR + 2016*, Sun Valley, Idaho, May 1-5, 2016. + +- Antonios G. Mylonakis, Melpomeni Varvayanni, and Nicolas Catsaros, + "Investigating a Matrix-free, Newton-based, Neutron-Monte + Carlo/Thermal-Hydraulic Coupling Scheme", *Proc. Int. Conf. Nuclear Energy for + New Europe*, Portoroz, Slovenia, Sep .14-17 + (2015). ``_ + - Matt Ellis, Benoit Forget, Kord Smith, and Derek Gaston, "Preliminary coupling of the Monte Carlo code OpenMC and the Multiphysics Object-Oriented Simulation Environment (MOOSE) for analyzing Doppler feedback in Monte Carlo @@ -80,8 +90,17 @@ Geometry Miscellaneous ------------- +- Yunzhao Li, Qingming He, Liangzhi Cao, Hongchun Wu, and Tiejun Zu, "Resonance + Elastic Scattering and Interference Effects Treatments in Subgroup Method," + *Nucl. Eng. Tech.*, **48**, 339-350 + (2016). ``_ + - William Boyd, Sterling Harper, and Paul K. Romano, "Equipping OpenMC for the - big data era," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016. + big data era," *Proc. PHYSOR*, Sun Valley, Idaho, May 1-5, 2016. + +- Michal Kostal, Vojtech Rypar, Jan Milcak, Vlastimil Juricek, Evzen Losa, + Benoit Forget, and Sterling Harper, *Ann. Nucl. Energy*, **87**, 601-611 + (2016). ``_ - Qicang Shen, William Boyd, Benoit Forget, and Kord Smith, "Tally precision triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**, @@ -95,6 +114,11 @@ Miscellaneous Multi-group Cross Section Generation ------------------------------------ +- Zhaoyuan Liu, Kord Smith, and Benoit Forget, "A Cumulative Migration Method + for Computing Rigorous Transport Cross Sections and Diffusion Coefficients for + LWR Lattices with Monte Carlo," *Proc. PHYSOR*, Sun Valley, Idaho, May + 1-5, 2016. + - Adam G. Nelson and William R. Martin, "Improved Monte Carlo tallying of multi-group scattering moments using the NDPP code," *Trans. Am. Nucl. Soc.*, **113**, 645-648 (2015) @@ -108,18 +132,43 @@ Multi-group Cross Section Generation Computational Methods Applied to Nuclear Science and Engineering*, Sun Valley, Idaho, May 5--9 (2013). ------------- -Nuclear Data ------------- + +------------------ +Doppler Broadening +------------------ - Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole - for cross section Doppler broadening," *J. Comput. Phys.*, In Press + for cross section Doppler broadening," *J. Comput. Phys.*, **307**, 715-727 (2016). ``_ +- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown, + "On-the-fly Doppler Broadening of Unresolved Resonance Region Cross Sections + via Probability Band Interpolation," *Proc. PHYSOR*, Sun Valley, Idaho, May + 1-5, 2016. + - Colin Josey, Benoit Forget, and Kord Smith, "Windowed multipole sensitivity to target accuracy of the optimization procedure," *J. Nucl. Sci. Technol.*, **52**, 987-992 (2015). ``_ +- Paul K. Romano and Timothy H. Trumbull, "Comparison of algorithms for Doppler + broadening pointwise tabulated cross sections," *Ann. Nucl. Energy*, **75**, + 358--364 (2015). ``_ + +- Tuomas Viitanen, Jaakko Leppanen, and Benoit Forget, "Target motion sampling + temperature treatment technique with track-length esimators in OpenMC -- + Preliminary results," *Proc. PHYSOR*, Kyoto, Japan, Sep. 28--Oct. 3 (2014). + +- Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte + Carlo simulations using the multipole representation," *Ann. Nucl. Energy*, + **64**, 78--85 (2014). ``_ + +------------ +Nuclear Data +------------ + +- Paul K. Romano and Sterling M. Harper, "Nuclear data processing capabilities + in OpenMC", *Proc. Nuclear Data*, Sep. 11-16, 2016. + - Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith, "Optimizations of the energy grid search algorithm in continuous-energy Monte Carlo particle transport codes", *Comput. Phys. Commun.*, **196**, 134-142 @@ -139,29 +188,17 @@ Nuclear Data performance analysis for varying cross section parameter regimes," *Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). -- Paul K. Romano and Timothy H. Trumbull, "Comparison of algorithms for Doppler - broadening pointwise tabulated cross sections," *Ann. Nucl. Energy*, **75**, - 358--364 (2015). ``_ - -- Tuomas Viitanen, Jaakko Leppanen, and Benoit Forget, "Target motion sampling - temperature treatment technique with track-length esimators in OpenMC -- - Preliminary results," *Proc. PHYSOR*, Kyoto, Japan, Sep. 28--Oct. 3 (2014). - - Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated sampling of the free gas resonance elastic scattering kernel," *Ann. Nucl. Energy*, **69**, 116--124 (2014). ``_ -- Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte - Carlo simulations using the multipole representation," *Ann. Nucl. Energy*, - **64**, 78--85 (2014). ``_ - ----------- Parallelism ----------- - Paul K. Romano, John R. Tramm, and Andrew R. Siegel, "Efficacy of hardware threading for Monte Carlo particle transport calculations on multi- and - many-core systems," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016. + many-core systems," *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016. - David Ozog, Allen D. Malony, and Andrew R. Siegel, "A performance analysis of SIMD algorithms for Monte Carlo simulations of nuclear reactor cores," @@ -228,3 +265,11 @@ Parallelism - Paul K. Romano and Benoit Forget, "Parallel Fission Bank Algorithms in Monte Carlo Criticality Calculations," *Nucl. Sci. Eng.*, **170**, 125--135 (2012). ``_ + +--------- +Depletion +--------- + +- Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion + chain simplification based of significance analysis," *Proc. PHYSOR*, Sun + Valley, Idaho, May 1-5, 2016. diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 960bee4ca..9b2d78607 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -725,13 +725,16 @@ a material default temperature. ```` Element -------------------------------- -The ```` element has an accepted value of "nearest" or -"interpolation". A value of "nearest" indicates that for each cell, the nearest -temperature at which cross sections are given is to be applied, within a given -tolerance (see :ref:`temperature_tolerance`). A value of "multipole" indicates -that the windowed multipole method should be used to evaluate -temperature-dependent cross sections in the resolved resonance range (a -:ref:`windowed multipole library ` must also be available). +The ```` element has an accepted value of "nearest", +"interpolation", or "multipole". A value of "nearest" indicates that for each +cell, the nearest temperature at which cross sections are given is to be +applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of +"interpolation" indicates that cross sections are to be linear-linear +interpolated between temperatures at which nuclear data are present (see +:ref:`temperature_treatment`). A value of "multipole" indicates that the +windowed multipole method should be used to evaluate temperature-dependent cross +sections in the resolved resonance range (a :ref:`windowed multipole library +` must also be available). *Default*: "nearest" diff --git a/openmc/geometry.py b/openmc/geometry.py index b2a9b5e4a..6d9330dc6 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -64,7 +64,7 @@ class Geometry(object): if cell.id in volume_calc.results: cell.add_volume_information(volume_calc) - def export_to_xml(self): + def export_to_xml(self, path='geometry.xml'): """Create a geometry.xml file that can be used for a simulation. """ @@ -82,8 +82,7 @@ class Geometry(object): # Write the XML Tree to the geometry.xml file tree = ET.ElementTree(geometry_file) - tree.write("geometry.xml", xml_declaration=True, encoding='utf-8', - method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") def find(self, point): """Find cells/universes/lattices which contain a given point diff --git a/openmc/material.py b/openmc/material.py index f56af485a..2861d5be6 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -41,19 +41,16 @@ class Material(object): name : str, optional Name of the material. If not specified, the name will be the empty string. - temperature : str, optional - The temperature identifier applied to this material. The units are - in Kelvin and the temperature rounded to the nearest integer. - For example, a tempreature of 293.6K would be provided as '294K' + temperature : float, optional + Temperature of the material in Kelvin. If not specified, the material + inherits the default temperature applied to the model. Attributes ---------- id : int Unique identifier for the material - temperature : str - The temperature identifier applied to this material. The units are - in Kelvin and the temperature rounded to the nearest integer. - For example, a tempreature of 293.6K would be provided as '294K' + temperature : float + Temperature of the material in Kelvin. density : float Density of the material (units defined separately) density_units : str @@ -217,12 +214,9 @@ class Material(object): @temperature.setter def temperature(self, temperature): - if temperature is not None: - cv.check_type('Temperature for Material ID="{0}"'.format(self._id), - temperature, basestring) - self._temperature = temperature - else: - self._temperature = '' + cv.check_type('Temperature for Material ID="{0}"'.format(self._id), + temperature, (Real, type(None))) + self._temperature = temperature def set_density(self, units, density=None): """Set the density of the material @@ -631,9 +625,9 @@ class Material(object): element.set("name", str(self._name)) # Create temperature XML subelement - if len(self.temperature) > 0: + if self.temperature is not None: subelement = ET.SubElement(element, "temperature") - subelement.text = self.temperature + subelement.text = str(self.temperature) # Create density XML subelement subelement = ET.SubElement(element, "density") @@ -817,7 +811,7 @@ class Materials(cv.CheckedList): xml_element = material.get_material_xml() self._materials_file.append(xml_element) - def export_to_xml(self): + def export_to_xml(self, path='materials.xml'): """Create a materials.xml file that can be used for a simulation. """ @@ -833,5 +827,4 @@ class Materials(cv.CheckedList): # Write the XML Tree to the materials.xml file tree = ET.ElementTree(self._materials_file) - tree.write("materials.xml", xml_declaration=True, - encoding='utf-8', method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") diff --git a/openmc/plots.py b/openmc/plots.py index cc5c0d44b..f5c1f2b2a 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -616,7 +616,7 @@ class Plots(cv.CheckedList): self._plots_file.append(xml_element) - def export_to_xml(self): + def export_to_xml(self, path='plots.xml'): """Create a plots.xml file that can be used by OpenMC. """ @@ -631,5 +631,4 @@ class Plots(cv.CheckedList): # Write the XML Tree to the plots.xml file tree = ET.ElementTree(self._plots_file) - tree.write("plots.xml", xml_declaration=True, - encoding='utf-8', method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") diff --git a/openmc/settings.py b/openmc/settings.py index c7d6debb9..15898a2b8 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1118,7 +1118,7 @@ class Settings(object): for r in self.resonance_scattering: elem.append(r.to_xml_element()) - def export_to_xml(self): + def export_to_xml(self, path='settings.xml'): """Create a settings.xml file that can be used for a simulation. """ @@ -1162,8 +1162,7 @@ class Settings(object): # Write the XML Tree to the settings.xml file tree = ET.ElementTree(self._settings_file) - tree.write("settings.xml", xml_declaration=True, - encoding='utf-8', method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") class ResonanceScattering(object): diff --git a/src/cross_section.F90 b/src/cross_section.F90 index e2cc31d5e..4bb323c51 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -162,15 +162,33 @@ contains ! temperature. Note that there is no tolerance here, so this ! temperature could be very far off! kT = sqrtkT**2 + i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) end if else - ! If not using multipole data, do a linear search on temperature kT = sqrtkT**2 - do i_temp = 1, size(nuclides(i_nuclide) % kTs) - if (abs(nuclides(i_nuclide) % kTs(i_temp) - kT) < & - K_BOLTZMANN*temperature_tolerance) exit - end do + + select case (temperature_method) + case (TEMPERATURE_NEAREST) + ! If using nearest temperature, do linear search on temperature + do i_temp = 1, size(nuc % kTs) + if (abs(nuc % kTs(i_temp) - kT) < K_BOLTZMANN * & + temperature_tolerance) exit + end do + case (TEMPERATURE_INTERPOLATION) + ! Find temperatures that bound the actual temperature + do i_temp = 1, size(nuc % kTs) - 1 + if (nuc % kTs(i_temp) <= kT .and. kT < nuc % kTs(i_temp + 1)) exit + end do + + ! Randomly sample between temperature i and i+1 + f = (kT - nuc % kTs(i_temp)) / & + (nuc % kTs(i_temp + 1) - nuc % kTs(i_temp)) + if (f > prn()) i_temp = i_temp + 1 + case (TEMPERATURE_MULTIPOLE) + i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) + end select + end if ! Evaluate multipole or interpolate @@ -317,10 +335,25 @@ contains ! Determine temperature for S(a,b) table kT = sqrtkT**2 - do i_temp = 1, size(sab_tables(i_sab) % kTs) - if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < & - K_BOLTZMANN*temperature_tolerance) exit - end do + if (temperature_method == TEMPERATURE_NEAREST) then + ! If using nearest temperature, do linear search on temperature + do i_temp = 1, size(sab_tables(i_sab) % kTs) + if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < & + K_BOLTZMANN*temperature_tolerance) exit + end do + else + ! Find temperatures that bound the actual temperature + do i_temp = 1, size(sab_tables(i_sab) % kTs) - 1 + if (sab_tables(i_sab) % kTs(i_temp) <= kT .and. & + kT < sab_tables(i_sab) % kTs(i_temp + 1)) exit + end do + + ! Randomly sample between temperature i and i+1 + f = (kT - sab_tables(i_sab) % kTs(i_temp)) / & + (sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp)) + if (f > prn()) i_temp = i_temp + 1 + end if + ! Get pointer to S(a,b) table associate (sab => sab_tables(i_sab) % data(i_temp)) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 04db716ba..ca7d4d2c0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4922,7 +4922,7 @@ contains call names % push_back('Ga0') call densities % push_back(density) else - call names % push_back('Ha69') + call names % push_back('Ga69') call densities % push_back(density * 0.60108_8) call names % push_back('Ga71') call densities % push_back(density * 0.39892_8) @@ -5814,7 +5814,7 @@ contains file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), & - temperature_tolerance) + temperature_method, temperature_tolerance) call close_group(group_id) call file_close(file_id) @@ -5967,6 +5967,7 @@ contains integer :: i, j integer :: i_library + integer :: method integer(HID_T) :: file_id integer(HID_T) :: group_id real(8) :: xs_cdf_sum @@ -5994,8 +5995,9 @@ contains ! Read nuclide data from HDF5 file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) + method = TEMPERATURE_NEAREST call resonant_nuc % from_hdf5(group_id, temperature, & - TEMPERATURE_NEAREST, 1000.0_8) + method, 1000.0_8) call close_group(group_id) call file_close(file_id) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 1fbd691ee..4235f4234 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -190,7 +190,7 @@ module nuclide_header class(Nuclide), intent(inout) :: this integer(HID_T), intent(in) :: group_id type(VectorReal), intent(in) :: temperature ! list of desired temperatures - integer, intent(in) :: method + integer, intent(inout) :: method real(8), intent(in) :: tolerance integer :: i @@ -241,11 +241,22 @@ module nuclide_header call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) temps_available(i) = temps_available(i) / K_BOLTZMANN end do + call sort(temps_available) + ! If only one temperature is available, revert to nearest temperature + if (size(temps_available) == 1 .and. & + method == TEMPERATURE_INTERPOLATION) then + call warning("Cross sections for " // trim(this % name) // " are only & + &available at one temperature. Reverting to nearest temperature & + &method.") + method = TEMPERATURE_NEAREST + end if + + ! Determine actual temperatures to read select case (method) case (TEMPERATURE_NEAREST) - ! Determine actual temperatures to read - TEMP_LOOP: do i = 1, temperature % size() + ! Find nearest temperatures + do i = 1, temperature % size() temp_desired = temperature % data(i) i_closest = minloc(abs(temps_available - temp_desired), dim=1) temp_actual = temps_available(i_closest) @@ -265,11 +276,31 @@ module nuclide_header &for " // trim(this % name) // " at or near " // & trim(to_str(nint(temp_desired))) // " K.") end if - end do TEMP_LOOP + end do case (TEMPERATURE_INTERPOLATION) - ! TODO: Get bounding temperatures - call fatal_error("Temperature interpolation not yet implemented") + ! If temperature interpolation or multipole is selected, get a list of + ! bounding temperatures for each actual temperature present in the model + TEMP_LOOP: do i = 1, temperature % size() + temp_desired = temperature % data(i) + + do j = 1, size(temps_available) - 1 + if (temps_available(j) <= temp_desired .and. & + temp_desired < temps_available(j + 1)) then + if (find(temps_to_read, nint(temps_available(j))) == -1) then + call temps_to_read % push_back(nint(temps_available(j))) + end if + if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then + call temps_to_read % push_back(nint(temps_available(j + 1))) + end if + cycle TEMP_LOOP + end if + end do + + call fatal_error("Nuclear data library does not contain cross sections & + &for " // trim(this % name) // " at temperatures that bound " // & + trim(to_str(nint(temp_desired))) // " K.") + end do TEMP_LOOP case (TEMPERATURE_MULTIPOLE) ! Add first available temperature diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 147643065..21a2ab032 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -80,10 +80,11 @@ module sab_header contains - subroutine salphabeta_from_hdf5(this, group_id, temperature, tolerance) + subroutine salphabeta_from_hdf5(this, group_id, temperature, method, tolerance) class(SAlphaBeta), intent(inout) :: this integer(HID_T), intent(in) :: group_id type(VectorReal), intent(in) :: temperature ! list of temperatures + integer, intent(in) :: method real(8), intent(in) :: tolerance integer :: i, j @@ -142,25 +143,55 @@ contains call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) temps_available(i) = temps_available(i) / K_BOLTZMANN end do + call sort(temps_available) - ! Determine actual temperatures to read - TEMP_LOOP: do i = 1, temperature % size() - temp_desired = temperature % data(i) - i_closest = minloc(abs(temps_available - temp_desired), dim=1) - temp_actual = temps_available(i_closest) - if (abs(temp_actual - temp_desired) < tolerance) then - if (find(temps_to_read, nint(temp_actual)) == -1) then - call temps_to_read % push_back(nint(temp_actual)) + select case (method) + case (TEMPERATURE_NEAREST) + ! Determine actual temperatures to read + do i = 1, temperature % size() + temp_desired = temperature % data(i) + i_closest = minloc(abs(temps_available - temp_desired), dim=1) + temp_actual = temps_available(i_closest) + if (abs(temp_actual - temp_desired) < tolerance) then + if (find(temps_to_read, nint(temp_actual)) == -1) then + call temps_to_read % push_back(nint(temp_actual)) + end if + else + call fatal_error("Nuclear data library does not contain cross sections & + &for " // trim(this % name) // " at or near " // & + trim(to_str(nint(temp_desired))) // " K.") end if - else - call fatal_error("Nuclear data library does not contain cross sections & - &for " // trim(this % name) // " at or near " // & - trim(to_str(nint(temp_desired))) // " K.") - end if - end do TEMP_LOOP + end do - ! TODO: If using interpolation, add a block to add bounding temperatures for - ! each + case (TEMPERATURE_INTERPOLATION) + ! If temperature interpolation or multipole is selected, get a list of + ! bounding temperatures for each actual temperature present in the model + TEMP_LOOP: do i = 1, temperature % size() + temp_desired = temperature % data(i) + + do j = 1, size(temps_available) - 1 + if (temps_available(j) <= temp_desired .and. & + temp_desired < temps_available(j + 1)) then + if (find(temps_to_read, nint(temps_available(j))) == -1) then + call temps_to_read % push_back(nint(temps_available(j))) + end if + if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then + call temps_to_read % push_back(nint(temps_available(j + 1))) + end if + cycle TEMP_LOOP + end if + end do + + call fatal_error("Nuclear data library does not contain cross sections & + &for " // trim(this % name) // " at temperatures that bound " // & + trim(to_str(nint(temp_desired))) // " K.") + end do TEMP_LOOP + + case (TEMPERATURE_MULTIPOLE) + ! Add first available temperature + call temps_to_read % push_back(nint(temps_available(1))) + + end select ! Sort temperatures to read call sort(temps_to_read) @@ -301,6 +332,12 @@ contains end do end associate end do + + ! Clear data on correlated angle-energy object + deallocate(correlated_dist % breakpoints) + deallocate(correlated_dist % interpolation) + deallocate(correlated_dist % energy) + deallocate(correlated_dist % distribution) end if call close_group(inelastic_group) diff --git a/tests/test_source/test_source.py b/tests/test_source/test_source.py index ce9012bc2..26844ea4e 100644 --- a/tests/test_source/test_source.py +++ b/tests/test_source/test_source.py @@ -13,7 +13,7 @@ import openmc class SourceTestHarness(PyAPITestHarness): def _build_inputs(self): - mat1 = openmc.Material(material_id=1, temperature='294') + mat1 = openmc.Material(material_id=1, temperature=294) mat1.set_density('g/cm3', 4.5) mat1.add_nuclide(openmc.Nuclide('U235'), 1.0) materials = openmc.Materials([mat1])