From f747bb4fbd253ba06c3d8527ee8aaff56d9ceade Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 5 Jun 2016 15:09:32 -0400 Subject: [PATCH] Exposed tabular_legendre to the openmc.Library interface to the MGXS Data library. Also ever-so slightly sped up the tabular scattering runtime --- openmc/mgxs/library.py | 68 ++++++++++++++++++++++++++++++++++++---- openmc/mgxs_library.py | 2 +- src/scattdata_header.F90 | 21 +++++++------ 3 files changed, 75 insertions(+), 16 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 62dde28ab..76bdd885a 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -753,7 +753,8 @@ class Library(object): return pickle.load(open(full_filename, 'rb')) def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', - xs_id='1m', order=None): + xs_id='1m', order=None, tabular_legendre=None, + tabular_points=33): """Generates an openmc.XSdata object describing a multi-group cross section data set for eventual combination in to an openmc.MGXSLibrary object (i.e., the library). @@ -773,9 +774,23 @@ class Library(object): nuclide this will be set to 'macro' regardless. xs_ids : str Cross section set identifier. Defaults to '1m'. - order : Scattering order for this data entry. Default is None, + order : int + Scattering order for this data entry. Default is None, which will set the XSdata object to use the order of the Library. + tabular_legendre : {None, bool} + Flag to denote whether or not the Legendre expansion of the + scattering angular distribution is to be converted to a tabular + representation by OpenMC. A value of `True` means that it is to be + converted while a value of 'False' means that it will not be. + Defaults to `None` which leaves the default behavior of OpenMC in + place (the distribution is not converted to a tabular + representation). + tabular_points : {int} + This parameter is not used unless the `tabular_legendre` is set to + `True`. In this case, this parameter sets the number of + equally-spaced points in the domain of [-1,1] to be used in + building the tabular distribution. Default is `33`. Returns ------- @@ -804,6 +819,10 @@ class Library(object): if order is not None: cv.check_greater_than('order', order, 0, equality=True) cv.check_less_than('order', order, 10, equality=True) + cv.check_type('tabular_legendre', tabular_legendre, + (type(None), bool)) + if tabular_points is not None: + cv.check_greater_than('tabular_points', tabular_points, 1) # Make sure statepoint has been loaded if self._sp_filename is None: @@ -830,6 +849,11 @@ class Library(object): # the provided order or the Library's order. xsdata.order = min(order, self.legendre_order) + # Set the tabular_legendre option if needed + if tabular_legendre is not None: + xsdata.tabular_legendre = {'enable': tabular_legendre, + 'num_points': tabular_points} + if nuclide is not 'total': xsdata.zaid = self._nuclides[nuclide][0] xsdata.awr = self._nuclides[nuclide][1] @@ -906,7 +930,8 @@ class Library(object): return xsdata def create_mg_library(self, xs_type='macro', xsdata_names=None, - xs_ids=None): + xs_ids=None, tabular_legendre=None, + tabular_points=33): """Creates an openmc.MGXSLibrary object to contain the MGXS data for the Multi-Group mode of OpenMC. @@ -923,6 +948,19 @@ class Library(object): Cross section set identifier (i.e., '71c') for all data sets (if only str) or for each individual one (if iterable of str). Defaults to '1m'. + tabular_legendre : {None, bool} + Flag to denote whether or not the Legendre expansion of the + scattering angular distribution is to be converted to a tabular + representation by OpenMC. A value of `True` means that it is to be + converted while a value of 'False' means that it will not be. + Defaults to `None` which leaves the default behavior of OpenMC in + place (the distribution is not converted to a tabular + representation). + tabular_points : {int} + This parameter is not used unless the `tabular_legendre` is set to + `True`. In this case, this parameter sets the number of + equally-spaced points in the domain of [-1,1] to be used in + building the tabular distribution. Default is `33`. Returns ------- @@ -983,13 +1021,16 @@ class Library(object): xsdata_name += '_' + nuclide xsdata = self.get_xsdata(domain, xsdata_name, nuclide=nuclide, - xs_type=xs_type, xs_id=xs_ids[i]) + xs_type=xs_type, xs_id=xs_ids[i], + tabular_legendre=tabular_legendre, + tabular_points=tabular_points) mgxs_file.add_xsdata(xsdata) return mgxs_file - def create_mg_mode(self, xsdata_names=None, xs_ids=None): + def create_mg_mode(self, xsdata_names=None, xs_ids=None, + tabular_legendre=None, tabular_points=33): """Creates an openmc.MGXSLibrary object to contain the MGXS data for the Multi-Group mode of OpenMC as well as the associated openmc.Materials and openmc.Geometry objects. The created Geometry is the same as that @@ -1007,6 +1048,19 @@ class Library(object): Cross section set identifier (i.e., '71c') for all data sets (if only str) or for each individual one (if iterable of str). Defaults to '1m'. + tabular_legendre : {None, bool} + Flag to denote whether or not the Legendre expansion of the + scattering angular distribution is to be converted to a tabular + representation by OpenMC. A value of `True` means that it is to be + converted while a value of 'False' means that it will not be. + Defaults to `None` which leaves the default behavior of OpenMC in + place (the distribution is not converted to a tabular + representation). + tabular_points : {int} + This parameter is not used unless the `tabular_legendre` is set to + `True`. In this case, this parameter sets the number of + equally-spaced points in the domain of [-1,1] to be used in + building the tabular distribution. Default is `33`. Returns ------- @@ -1071,7 +1125,9 @@ class Library(object): # Create XSdata and Macroscopic for this domain xsdata = self.get_xsdata(domain, xsdata_name, nuclide='total', - xs_type=xs_type, xs_id=xs_ids[i]) + xs_type=xs_type, xs_id=xs_ids[i], + tabular_legendre=tabular_legendre, + tabular_points=tabular_points) mgxs_file.add_xsdata(xsdata) macroscopic = openmc.Macroscopic(name=xsdata_name, xs=xs_ids[i]) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 35b48d873..f34f3b71c 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -445,7 +445,7 @@ class XSdata(object): raise ValueError(msg) if 'num_points' in tabular_legendre: num_points = tabular_legendre['num_points'] - check_value('num_points', num_points, Integral) + check_type('num_points', num_points, Integral) check_greater_than('num_points', num_points, 0) else: if not enable: diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 12c11e2e4..06e052eec 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -329,7 +329,7 @@ contains real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use integer :: imu, gin, gout, groups, order - real(8) :: norm + real(8) :: norm, m, mu0, mu1, p0 real(8), allocatable :: energy(:, :) real(8), allocatable :: matrix(:, :, :) @@ -410,13 +410,17 @@ contains this % fmu(gin) % data(:, gout) / norm end if - ! Now create CDF from fmu with trapezoidal rule + ! Now create CDF from fmu with the analytical integral this % dist(gin) % data(1, gout) = ZERO do imu = 2, order - this % dist(gin) % data(imu, gout) = & - this % dist(gin) % data(imu - 1, gout) + & - HALF * this % dmu * (this % fmu(gin) % data(imu - 1, gout) + & - this % fmu(gin) % data(imu, gout)) + p0 = this % fmu(gin) % data(imu - 1, gout) + mu0 = this % mu(imu - 1) + mu1 = this % mu(imu) + m = (this % fmu(gin) % data(imu, gout) - p0) / (mu1 - mu0) + this % dist(gin) % data(imu, gout) = HALF * m * mu1 * mu1 + & + (p0 - m * mu0) * mu1 + & + (HALF * m * mu0 * mu0 - p0 * mu0) + end do ! Ensure we normalize to 1 still norm = this % dist(gin) % data(order, gout) @@ -628,11 +632,10 @@ contains p1 = this % fmu(gin) % data(k + 1, gout) mu1 = this % mu(k + 1) - frac = (p1 - p0) / (mu1 - mu0) - - if (frac == ZERO) then + if (p0 == p1) then mu = mu0 + (xi - c_k) / p0 else + frac = (p1 - p0) / (mu1 - mu0) mu = mu0 + & (sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac end if