From 2abc9f1e2940cc188e80f0491b3b99bbd95147b5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Nov 2015 16:41:44 -0600 Subject: [PATCH] Support for tabular angular and energy distributions. Implemented new polymorphic Distribution type which is used by AngleDistribution, both of which appear in distribution_header. --- docs/source/usersguide/input.rst | 60 +++++++-- openmc/settings.py | 78 ++++++++--- src/distribution_header.F90 | 216 +++++++++++++++++++++++++++++++ src/error.F90 | 15 +-- src/global.F90 | 11 +- src/input_xml.F90 | 214 ++++++++++++++++++++---------- src/math.F90 | 49 +++++++ src/physics.F90 | 51 +------- src/random_lcg.F90 | 5 +- src/relaxng/settings.rnc | 9 +- src/relaxng/settings.rng | 36 +----- src/source.F90 | 67 +++------- src/source_header.F90 | 8 +- src/state_point.F90 | 1 + 14 files changed, 563 insertions(+), 257 deletions(-) create mode 100644 src/distribution_header.F90 diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index fa43ca1d2c..ef9be8123f 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -440,7 +440,6 @@ attributes/sub-elements: has the following attributes: :type: - The type of spatial distribution. Valid options are "box", "fission", and "point". A "box" spatial distribution has coordinates sampled uniformly in a parallelepiped. A "fission" spatial distribution samples locations from @@ -468,14 +467,23 @@ attributes/sub-elements: has the following attributes: :type: - The type of angular distribution. Valid options are "isotropic" and - "monodirectional". The angle of the particle emitted from a source site is - isotropic if the "isotropic" option is given. The angle of the particle - emitted from a source site is the direction specified in the - attribute if "monodirectional" option is given. + The type of angular distribution. Valid options are "isotropic", + "monodirectional", and "tabular". The angle of the particle emitted from a + source site is isotropic if the "isotropic" option is given. The angle of + the particle emitted from a source site is the direction specified in the + attribute if "monodirectional" option is given. The "tabular" + option produces directions with polar angles sampled from a tabulated + distribution. *Default*: isotropic + :interpolation: + For a "tabular" angular distribution, ``interpolation`` can be set to + "histogram" or "linear-linear" thereby specifying how tabular points are + to be interpolated. + + *Default*: histogram + :parameters: For an "isotropic" angular distribution, ``parameters`` should not be specified. @@ -484,6 +492,22 @@ attributes/sub-elements: given as three real numbers which specify the angular cosines with respect to each axis. + For a "tabular" angular distribution, ``parameters`` provides the + :math:`(\mu,p)` pairs defining the tabular distribution. All :math:`\mu` + points are given first followed by corresponding :math:`p` points. The + following example gives a histogram distribution with even probability of + selecting a polar angle in the range [-1,-0.5] and [0.5,1] (Note that the + last :math:`p` point is inconsequential): + + .. code-block:: xml + + + + -1.0 -0.5 0.5 1.0 + 1.0 0.0 1.0 0.0 + + + *Default*: None :energy: @@ -491,15 +515,23 @@ attributes/sub-elements: has the following attributes: :type: - The type of energy distribution. Valid options are "monoenergetic", - "watt", and "maxwell". The "monoenergetic" option produces source sites at - a single energy. The "watt" option produces source sites whose energy is - sampled from a Watt fission spectrum. The "maxwell" option produce source - sites whose energy is sampled from a Maxwell fission spectrum. + "watt", "maxwell", and "tabular". The "monoenergetic" option produces + source sites at a single energy. The "watt" option produces source sites + whose energy is sampled from a Watt fission spectrum. The "maxwell" option + produce source sites whose energy is sampled from a Maxwell fission + spectrum. The "tabular" option produces source sites whose energy is + sampled from a tabulated distribution. *Default*: watt + :interpolation: + For a "tabular" angular distribution, ``interpolation`` can be set to + "histogram" or "linear-linear" thereby specifying how tabular points are + to be interpolated. + + *Default*: histogram + :parameters: For a "monoenergetic" energy distribution, ``parameters`` should be given as the energy in MeV of the source sites. @@ -512,6 +544,10 @@ attributes/sub-elements: real number :math:`a` that parameterizes the distribution :math:`p(E) dE = c E e^{-E/a} dE`. + For a "tabular" energy distribution, ``parameters`` provides the + :math:`(E,p)` pairs defining the tabular distribution. All :math:`E` + points are given first followed by corresponding :math:`p` points. + *Default*: 0.988 2.249 :write_initial: @@ -519,7 +555,7 @@ attributes/sub-elements: the beginning of the first batch. The output file is named "initial_source.binary(h5)" - *Default*: false + *Default*: false ```` Element ------------------------- diff --git a/openmc/settings.py b/openmc/settings.py index e8ec6463af..2cb6a14460 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -139,8 +139,10 @@ class SettingsFile(object): self._source_space_type = None self._source_space_params = None self._source_angle_type = None + self._source_angle_interpolation = None self._source_angle_params = None self._source_energy_type = None + self._source_energy_interpolation = None self._source_energy_params = None self._confidence_intervals = None @@ -243,6 +245,10 @@ class SettingsFile(object): def source_angle_type(self): return self._source_angle_type + @property + def source_angle_interpolation(self): + return self._source_angle_interpolation + @property def source_angle_params(self): return self._source_angle_params @@ -251,6 +257,10 @@ class SettingsFile(object): def source_energy_type(self): return self._source_energy_type + @property + def source_energy_interpolation(self): + return self._source_energy_interpolation + @property def source_energy_params(self): return self._source_energy_params @@ -511,17 +521,19 @@ class SettingsFile(object): self._source_space_type = stype self._source_space_params = params - def set_source_angle(self, stype, params=[]): + def set_source_angle(self, stype, params=[], interp='histogram'): """Defined the angular distribution of the external/starting source. Parameters ---------- stype : str - The type of angular distribution. Valid options are "isotropic" and - "monodirectional". The angle of the particle emitted from a source - site is isotropic if the "isotropic" option is given. The angle of - the particle emitted from a source site is the direction specified - in ``params`` if the "monodirectional" option is given. + The type of angular distribution. Valid options are "isotropic", + "monodirectional", and "tabular". The angle of the particle emitted + from a source site is isotropic if the "isotropic" option is + given. The angle of the particle emitted from a source site is the + direction specified in ``params`` if the "monodirectional" option is + given. The "tabular" option produces directions with polar angles + sampled from a tabulated distribution. params : Iterable of float For an "isotropic" angular distribution, ``params`` should not be specified. @@ -530,11 +542,20 @@ class SettingsFile(object): be given as three floats which specify the angular cosines with respect to each axis. + For a "tabular" angular distribution, ``parameters`` provides the + :math:`(\mu,p)` pairs defining the tabular distribution. All + :math:`\mu` points are given first followed by corresponding + :math:`p` points. + interp : { 'histogram', 'linear-linear' } + For a "tabular" angular distribution, ``interpolation`` can be set + to "histogram" or "linear-linear" thereby specifying how tabular + points are to be interpolated. + """ check_type('source angle type', stype, basestring) check_value('source angle type', stype, - ['isotropic', 'monodirectional']) + ['isotropic', 'monodirectional', 'tabular']) check_type('source angle parameters', params, Iterable, Real) if stype == 'isotropic' and params is not None: msg = 'Unable to set source angle parameters since they are not ' \ @@ -543,22 +564,30 @@ class SettingsFile(object): elif stype == 'monodirectional': check_length('source angle parameters for a monodirectional ' 'source', params, 3) + elif stype == 'tabular': + check_type('source angle interpolation', interp, basestring) + check_value('source angle interpolation', interp, + ['histogram', 'linear-linear']) + self._source_angle_interpolation = interp self._source_angle_type = stype self._source_angle_params = params - def set_source_energy(self, stype, params=[]): + def set_source_energy(self, stype, params=[], interp='histogram'): """Defined the energy distribution of the external/starting source. Parameters ---------- stype : str + The type of energy distribution. Valid options are "monoenergetic", - "watt", and "maxwell". The "monoenergetic" option produces source - sites at a single energy. The "watt" option produces source sites - whose energy is sampled from a Watt fission spectrum. The "maxwell" - option produce source sites whose energy is sampled from a Maxwell - fission spectrum. + "watt", "maxwell", and "tabular". The "monoenergetic" option + produces source sites at a single energy. The "watt" option produces + source sites whose energy is sampled from a Watt fission + spectrum. The "maxwell" option produce source sites whose energy is + sampled from a Maxwell fission spectrum. The "tabular" option + produces source sites whose energy is sampled from a tabulated + distribution. params : Iterable of float For a "monoenergetic" energy distribution, ``params`` should be given as the energy in MeV of the source sites. @@ -571,11 +600,19 @@ class SettingsFile(object): one real number :math:`a` that parameterizes the distribution :math:`p(E) dE = c E e^{-E/a} dE`. + For a "tabular" energy distribution, ``parameters`` provides the + :math:`(E,p)` pairs defining the tabular distribution. All :math:`E` + points are given first followed by corresponding :math:`p` points. + interp : { 'histogram', 'linear-linear' } + For a "tabular" energy distribution, ``interpolation`` can be set + to "histogram" or "linear-linear" thereby specifying how tabular + points are to be interpolated. + """ check_type('source energy type', stype, basestring) check_value('source energy type', stype, - ['monoenergetic', 'watt', 'maxwell']) + ['monoenergetic', 'watt', 'maxwell', 'tabular']) check_type('source energy parameters', params, Iterable, Real) if stype in ['monoenergetic', 'maxwell']: check_length('source energy parameters for a monoenergetic ' @@ -583,6 +620,11 @@ class SettingsFile(object): elif stype == 'watt': check_length('source energy parameters for a Watt source', params, 2) + elif stype == 'tabular': + check_type('source energy interpolation', interp, basestring) + check_value('source energy interpolation', interp, + ['histogram', 'linear-linear']) + self._source_energy_interpolation = interp self._source_energy_type = stype self._source_energy_params = params @@ -949,6 +991,9 @@ class SettingsFile(object): element = ET.SubElement(self._source_subelement, "angle") element.set("type", self._source_angle_type) + if self.source_angle_interpolation is not None: + element.set("interpolation", self.source_angle_interpolation) + subelement = ET.SubElement(element, "parameters") subelement.text = ' '.join(map(str, self._source_angle_params)) @@ -961,7 +1006,10 @@ class SettingsFile(object): element = ET.SubElement(self._source_subelement, "energy") element.set("type", self._source_energy_type) - subelement = ET.SubElement(element, "parameters") + if self.source_energy_interpolation is not None: + element.set("interpolation", self.source_energy_interpolation) + + subelement = ET.SubElement(element, "parameters") subelement.text = ' '.join(map(str, self._source_energy_params)) def _create_output_subelement(self): diff --git a/src/distribution_header.F90 b/src/distribution_header.F90 new file mode 100644 index 0000000000..9b10f636b3 --- /dev/null +++ b/src/distribution_header.F90 @@ -0,0 +1,216 @@ +module distribution_header + + use constants, only: ZERO, HALF, ONE, HISTOGRAM, LINEAR_LINEAR + use error, only: fatal_error + use math, only: rotate_angle, maxwell_spectrum, watt_spectrum + use random_lcg, only: prn + +!=============================================================================== +! DISTRIBUTION type defines a probability density function +!=============================================================================== + + type, abstract :: Distribution + contains + procedure(iSample), deferred :: sample + end type Distribution + + abstract interface + function iSample(this) result(x) + import Distribution + class(Distribution), intent(in) :: this + real(8) :: x + end function iSample + end interface + +!=============================================================================== +! Derived classes of Distribution +!=============================================================================== + + ! delta function at a single point + type, extends(Distribution) :: Delta + real(8) :: x0 + contains + procedure :: sample => delta_sample + end type Delta + + ! Uniform distribution over the interval [a,b] + type, extends(Distribution) :: Uniform + real(8) :: a + real(8) :: b + contains + procedure :: sample => uniform_sample + end type Uniform + + ! Maxwellian distribution of form c*E*exp(-E/a) + type, extends(Distribution) :: Maxwell + real(8) :: theta + contains + procedure :: sample => maxwell_sample + end type Maxwell + + ! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E)) + type, extends(Distribution) :: Watt + real(8) :: a + real(8) :: b + contains + procedure :: sample => watt_sample + end type Watt + + ! Histogram or linear-linear interpolated tabular distribution + type, extends(Distribution) :: Tabular + integer :: interpolation + real(8), allocatable :: x(:) ! tabulated independent variable + real(8), allocatable :: p(:) ! tabulated probability density + real(8), allocatable, private :: c(:) ! cumulative distribution at tabulated values + contains + procedure :: sample => tabular_sample + procedure :: initialize => tabular_initialize + end type Tabular + +!=============================================================================== +! AngleDistribution +!=============================================================================== + + type :: AngleDistribution + real(8) :: reference_uvw(3) + class(Distribution), allocatable :: mu + contains + procedure :: sample => angle_sample + end type AngleDistribution + +contains + + function delta_sample(this) result(x) + class(Delta), intent(in) :: this + real(8) :: x + + x = this%x0 + end function delta_sample + + function uniform_sample(this) result(x) + class(Uniform), intent(in) :: this + real(8) :: x + + x = this%a + prn()*(this%b - this%a) + end function uniform_sample + + function maxwell_sample(this) result(x) + class(Maxwell), intent(in) :: this + real(8) :: x + + x = maxwell_spectrum(this%theta) + end function maxwell_sample + + function watt_sample(this) result(x) + class(Watt), intent(in) :: this + real(8) :: x + + x = watt_spectrum(this%a, this%b) + end function watt_sample + + function tabular_sample(this) result(x) + class(Tabular), intent(in) :: this + real(8) :: x + + integer :: i + real(8) :: c ! sampled cumulative frequency + real(8) :: m ! slope of PDF + real(8) :: x_i, x_i1 ! i-th and (i+1)th x values + real(8) :: c_i, c_i1 ! i-th and (i+1)th cumulative distribution values + real(8) :: p_i, p_i1 ! i-th and (i+1)th probability density values + + ! Sample value of CDF + c = prn() + + ! Find first CDF bin which is above the sampled value + c_i = this%c(1) + do i = 1, size(this%c) - 1 + c_i1 = this%c(i + 1) + if (c <= c_i1) exit + c_i = c_i1 + end do + + ! Determine bounding PDF values + x_i = this%x(i) + p_i = this%p(i) + + if (this%interpolation == HISTOGRAM) then + ! Histogram interpolation + if (p_i > ZERO) then + x = x_i + (c - c_i)/p_i + else + x = x_i + end if + else + ! Linear-linear interpolation + x_i1 = this%x(i + 1) + p_i1 = this%p(i + 1) + + m = (p_i1 - p_i)/(x_i1 - x_i) + if (m == ZERO) then + x = x_i + (c - c_i)/p_i + else + x = x_i + (sqrt(max(ZERO, p_i*p_i + 2*m*(c - c_i))) - p_i)/m + end if + end if + end function tabular_sample + + subroutine tabular_initialize(this, x, p, interp) + class(Tabular), intent(inout) :: this + real(8), intent(in) :: x(:) + real(8), intent(in) :: p(:) + integer, intent(in) :: interp + + integer :: n + + ! Check interpolation parameter + + if (interp /= HISTOGRAM .and. interp /= LINEAR_LINEAR) then + call fatal_error('Only histogram and linear-linear interpolation for tabular & + &distribution is supported.') + end if + + ! Check length of x, p arrays + if (size(x) /= size(p)) then + call fatal_error('Tabulated probabilities not of same length as & + &independent variable.') + end if + + ! Copy probability density function and interpolation parameter + n = size(x) + allocate(this%x(n), this%p(n), this%c(n)) + this%interpolation = interp + this%x(:) = x(:) + this%p(:) = p(:) + + ! Calculate cumulative distribution function + this%c(1) = ZERO + do i = 2, n + if (this%interpolation == HISTOGRAM) then + this%c(i) = this%c(i-1) + this%p(i-1)*(this%x(i) - this%x(i-1)) + elseif (this%interpolation == LINEAR_LINEAR) then + this%c(i) = this%c(i-1) + HALF*(this%p(i-1) + this%p(i)) * & + (this%x(i) - this%x(i-1)) + end if + end do + + ! Normalize density and distribution functions + this%p(:) = this%p(:)/this%c(n) + this%c(:) = this%c(:)/this%c(n) + end subroutine tabular_initialize + + function angle_sample(this) result(uvw) + class(AngleDistribution), intent(in) :: this + real(8) :: uvw(3) + + real(8) :: mu + + mu = this%mu%sample() + if (mu == ONE) then + uvw(:) = this%reference_uvw + else + uvw(:) = rotate_angle(this%reference_uvw, mu) + end if + end function angle_sample + +end module distribution_header diff --git a/src/error.F90 b/src/error.F90 index ecd7eb4530..9e5ca84792 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -3,8 +3,6 @@ module error use, intrinsic :: ISO_FORTRAN_ENV use constants - use global - #ifdef MPI use message_passing #endif @@ -87,6 +85,9 @@ contains integer :: line_wrap ! length of line integer :: length ! length of message integer :: indent ! length of indentation +#ifdef MPI + integer :: mpi_err +#endif ! set default error code @@ -136,16 +137,6 @@ contains end if end do - ! Write information on current batch, generation, and particle - if (current_batch > 0) then - write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch - write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen - write(ERROR_UNIT,*) - end if - - ! Release memory from all allocatable arrays - call free_memory() - #ifdef MPI ! Abort MPI call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err) diff --git a/src/global.F90 b/src/global.F90 index 88ace73b6d..b030aab545 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -289,9 +289,6 @@ module global character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart character(MAX_FILE_LEN) :: path_output = '' ! Path to output directory - ! Random number seed - integer(8) :: seed = 1_8 - ! The verbosity controls how much information will be printed to the ! screen and in logs integer :: verbosity = 7 @@ -455,10 +452,10 @@ contains ! Deallocate external source if (allocated(external_source % params_space)) & deallocate(external_source % params_space) - if (allocated(external_source % params_angle)) & - deallocate(external_source % params_angle) - if (allocated(external_source % params_energy)) & - deallocate(external_source % params_energy) + if (allocated(external_source % angle % mu)) & + deallocate(external_source % angle % mu) + if (allocated(external_source % energy)) & + deallocate(external_source % energy) ! Deallocate k and entropy if (allocated(k_generation)) deallocate(k_generation) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 07d739b889..eb48e1d6f6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3,6 +3,7 @@ module input_xml use cmfd_input, only: configure_cmfd use constants use dict_header, only: DictIntInt, ElemKeyValueCI + use distribution_header use energy_grid, only: grid_method, n_log_bins use error, only: fatal_error, warning use geometry_header, only: Cell, Lattice, RectLattice, HexLattice @@ -11,7 +12,7 @@ module input_xml use mesh_header, only: RegularMesh use output, only: write_message use plot_header - use random_lcg, only: prn + use random_lcg, only: prn, seed use surface_header use stl_vector, only: VectorInt use string, only: to_lower, to_str, str_to_int, str_to_real, & @@ -59,6 +60,7 @@ contains integer :: temp_int_array3(3) integer, allocatable :: temp_int_array(:) integer(8) :: temp_long + real(8), allocatable :: temp_real(:) integer :: n_tracks logical :: file_exists character(MAX_FILE_LEN) :: env_variable @@ -423,24 +425,6 @@ contains ! Get pointer to angular distribution call get_node_ptr(node_source, "angle", node_dist) - ! Check for type of angular distribution - type = '' - if (check_for_node(node_dist, "type")) & - call get_node_value(node_dist, "type", type) - select case (to_lower(type)) - case ('isotropic') - external_source % type_angle = SRC_ANGLE_ISOTROPIC - coeffs_reqd = 0 - case ('monodirectional') - external_source % type_angle = SRC_ANGLE_MONO - coeffs_reqd = 3 - case ('tabular') - external_source % type_angle = SRC_ANGLE_TABULAR - case default - call fatal_error("Invalid angular distribution for external source: "& - &// trim(type)) - end select - ! Determine number of parameters specified if (check_for_node(node_dist, "parameters")) then n = get_arraysize_double(node_dist, "parameters") @@ -448,21 +432,76 @@ contains n = 0 end if + ! Check for type of angular distribution + type = '' + if (check_for_node(node_dist, "type")) & + call get_node_value(node_dist, "type", type) + select case (to_lower(type)) + case ('isotropic') + allocate(Uniform :: external_source%angle%mu) + + case ('monodirectional') + allocate(Delta :: external_source%angle%mu) + if (n /= 3) then + call fatal_error('Monodirectional angular distribution must have & + &three parameters specified.') + end if + + case ('tabular') + allocate(Tabular :: external_source%angle%mu) + + case default + call fatal_error("Invalid angular distribution for external source: "& + // trim(type)) + end select + + ! Set reference unit vector to be positive z-direction + external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE] + ! Read parameters for angle distribution - if (n < coeffs_reqd) then - call fatal_error("Not enough parameters specified for angle & - &distribution of external source.") - elseif (n > coeffs_reqd) then - call fatal_error("Too many parameters specified for angle & - &distribution of external source.") - elseif (n > 0) then - allocate(external_source % params_angle(n)) + select type (mu => external_source%angle%mu) + type is (Uniform) + mu%a = -ONE + mu%b = ONE + + type is (Delta) + mu%x0 = ONE call get_node_array(node_dist, "parameters", & - external_source % params_angle) - end if + external_source%angle%reference_uvw) + + type is (Tabular) + ! Read interpolation + if (check_for_node(node_source, "interpolation")) then + call get_node_value(node_source, "interpolation", temp_str) + select case(to_lower(temp_str)) + case ('histogram') + temp_int = HISTOGRAM + case ('linear-linear') + temp_int = LINEAR_LINEAR + case default + call fatal_error("Unknown interpolation type for source & + &angular distribution: " // trim(temp_str)) + end select + else + temp_int = HISTOGRAM + end if + + ! Read and initialize tabular distribution + allocate(temp_real(n)) + call get_node_array(node_dist, "parameters", temp_real) + call mu%initialize(temp_real(1:n), temp_real(n+1:2*n), temp_int) + deallocate(temp_real) + end select + else ! Set default angular distribution isotropic - external_source % type_angle = SRC_ANGLE_ISOTROPIC + allocate(Uniform :: external_source%angle%mu) + select type(mu => external_source%angle%mu) + type is (Uniform) + mu%a = -ONE + mu%b = ONE + end select + external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE] end if ! Determine external source energy distribution @@ -471,27 +510,6 @@ contains ! Get pointer to energy distribution call get_node_ptr(node_source, "energy", node_dist) - ! Check for type of energy distribution - type = '' - if (check_for_node(node_dist, "type")) & - call get_node_value(node_dist, "type", type) - select case (to_lower(type)) - case ('monoenergetic') - external_source % type_energy = SRC_ENERGY_MONO - coeffs_reqd = 1 - case ('maxwell') - external_source % type_energy = SRC_ENERGY_MAXWELL - coeffs_reqd = 1 - case ('watt') - external_source % type_energy = SRC_ENERGY_WATT - coeffs_reqd = 2 - case ('tabular') - external_source % type_energy = SRC_ENERGY_TABULAR - case default - call fatal_error("Invalid energy distribution for external source: " & - &// trim(type)) - end select - ! Determine number of parameters specified if (check_for_node(node_dist, "parameters")) then n = get_arraysize_double(node_dist, "parameters") @@ -499,23 +517,87 @@ contains n = 0 end if + ! Check for type of energy distribution + type = '' + if (check_for_node(node_dist, "type")) & + call get_node_value(node_dist, "type", type) + select case (to_lower(type)) + case ('monoenergetic') + allocate(Delta :: external_source%energy) + if (n /= 1) then + call fatal_error('Monoenergetic energy distribution must have one & + ¶meter specified.') + end if + + case ('maxwell') + allocate(Maxwell :: external_source%energy) + if (n /= 1) then + call fatal_error('Maxwell energy distribution must have one & + ¶meter specified.') + end if + + case ('watt') + allocate(Watt :: external_source%energy) + if (n /= 2) then + call fatal_error('Watt energy distribution must have two & + ¶meter specified.') + end if + + case ('tabular') + allocate(Tabular :: external_source%energy) + + case default + call fatal_error("Invalid energy distribution for external source: " & + // trim(type)) + end select + ! Read parameters for energy distribution - if (n < coeffs_reqd) then - call fatal_error("Not enough parameters specified for energy & - &distribution of external source.") - elseif (n > coeffs_reqd) then - call fatal_error("Too many parameters specified for energy & - &distribution of external source.") - elseif (n > 0) then - allocate(external_source % params_energy(n)) - call get_node_array(node_dist, "parameters", & - external_source % params_energy) - end if + select type(energy => external_source%energy) + type is (Delta) + call get_node_value(node_dist, "parameters", energy%x0) + + type is (Maxwell) + call get_node_value(node_dist, "parameters", energy%theta) + + type is (Watt) + allocate(temp_real(2)) + call get_node_array(node_dist, "parameters", temp_real) + energy%a = temp_real(1) + energy%b = temp_real(2) + deallocate(temp_real) + + type is (Tabular) + ! Read interpolation + if (check_for_node(node_dist, "interpolation")) then + call get_node_value(node_dist, "interpolation", temp_str) + select case(to_lower(temp_str)) + case ('histogram') + temp_int = HISTOGRAM + case ('linear-linear') + temp_int = LINEAR_LINEAR + case default + call fatal_error("Unknown interpolation type for source & + &angular distribution: " // trim(temp_str)) + end select + else + temp_int = HISTOGRAM + end if + + ! Read and initialize tabular distribution + allocate(temp_real(n)) + call get_node_array(node_dist, "parameters", temp_real) + call energy%initialize(temp_real(1:n), temp_real(n+1:2*n), temp_int) + deallocate(temp_real) + end select + else - ! Set default energy distribution to Watt fission spectrum - external_source % type_energy = SRC_ENERGY_WATT - allocate(external_source % params_energy(2)) - external_source % params_energy = (/ 0.988_8, 2.249_8 /) + ! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 + allocate(Watt :: external_source%energy) + select type(energy => external_source%energy) + type is (Watt) + energy%a = 0.988_8 + energy%b = 2.249_8 + end select end if end if diff --git a/src/math.F90 b/src/math.F90 index 15aa672e15..6b964684a2 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -557,6 +557,55 @@ contains end function calc_rn +!=============================================================================== +! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is +! mu and through an azimuthal angle sampled uniformly. Note that this is done +! with direct sampling rather than rejection as is done in MCNP and SERPENT. +!=============================================================================== + + function rotate_angle(uvw0, mu) result(uvw) + real(8), intent(in) :: uvw0(3) ! directional cosine + real(8), intent(in) :: mu ! cosine of angle in lab or CM + real(8) :: uvw(3) ! rotated directional cosine + + real(8) :: phi ! azimuthal angle + real(8) :: sinphi ! sine of azimuthal angle + real(8) :: cosphi ! cosine of azimuthal angle + real(8) :: a ! sqrt(1 - mu^2) + real(8) :: b ! sqrt(1 - w^2) + real(8) :: u0 ! original cosine in x direction + real(8) :: v0 ! original cosine in y direction + real(8) :: w0 ! original cosine in z direction + + ! Copy original directional cosines + u0 = uvw0(1) + v0 = uvw0(2) + w0 = uvw0(3) + + ! Sample azimuthal angle in [0,2pi) + phi = TWO * PI * prn() + + ! Precompute factors to save flops + sinphi = sin(phi) + cosphi = cos(phi) + a = sqrt(max(ZERO, ONE - mu*mu)) + b = sqrt(max(ZERO, ONE - w0*w0)) + + ! Need to treat special case where sqrt(1 - w**2) is close to zero by + ! expanding about the v component rather than the w component + if (b > 1e-10) then + uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b + uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b + uvw(3) = mu*w0 - a*b*cosphi + else + b = sqrt(ONE - v0*v0) + uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b + uvw(2) = mu*v0 - a*b*cosphi + uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + end if + + end function rotate_angle + !=============================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based ! on a direct sampling scheme. The probability distribution function for a diff --git a/src/physics.F90 b/src/physics.F90 index a01a3a30cb..217d232182 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -9,7 +9,7 @@ module physics use global use interpolation, only: interpolate_tab1 use material_header, only: Material - use math, only: maxwell_spectrum, watt_spectrum + use math, only: rotate_angle, maxwell_spectrum, watt_spectrum use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: Particle @@ -1517,55 +1517,6 @@ contains end function sample_angle -!=============================================================================== -! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is -! mu and through an azimuthal angle sampled uniformly. Note that this is done -! with direct sampling rather than rejection as is done in MCNP and SERPENT. -!=============================================================================== - - function rotate_angle(uvw0, mu) result(uvw) - real(8), intent(in) :: uvw0(3) ! directional cosine - real(8), intent(in) :: mu ! cosine of angle in lab or CM - real(8) :: uvw(3) ! rotated directional cosine - - real(8) :: phi ! azimuthal angle - real(8) :: sinphi ! sine of azimuthal angle - real(8) :: cosphi ! cosine of azimuthal angle - real(8) :: a ! sqrt(1 - mu^2) - real(8) :: b ! sqrt(1 - w^2) - real(8) :: u0 ! original cosine in x direction - real(8) :: v0 ! original cosine in y direction - real(8) :: w0 ! original cosine in z direction - - ! Copy original directional cosines - u0 = uvw0(1) - v0 = uvw0(2) - w0 = uvw0(3) - - ! Sample azimuthal angle in [0,2pi) - phi = TWO * PI * prn() - - ! Precompute factors to save flops - sinphi = sin(phi) - cosphi = cos(phi) - a = sqrt(max(ZERO, ONE - mu*mu)) - b = sqrt(max(ZERO, ONE - w0*w0)) - - ! Need to treat special case where sqrt(1 - w**2) is close to zero by - ! expanding about the v component rather than the w component - if (b > 1e-10) then - uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b - uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b - uvw(3) = mu*w0 - a*b*cosphi - else - b = sqrt(ONE - v0*v0) - uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b - uvw(2) = mu*v0 - a*b*cosphi - uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b - end if - - end function rotate_angle - !=============================================================================== ! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary ! neutron from a collision or for a prompt/delayed fission neutron diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 1ebe651d3c..8f50477c5b 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -7,6 +7,9 @@ module random_lcg private save + ! Random number seed + integer(8), public :: seed = 1_8 + integer(8) :: prn_seed0 ! original seed integer(8) :: prn_seed(N_STREAMS) ! current seed integer(8) :: prn_mult ! multiplication factor, g @@ -56,8 +59,6 @@ contains subroutine initialize_prng() - use global, only: seed - integer :: i prn_seed0 = seed diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index a95b1f297e..e29b4e22bf 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -69,27 +69,24 @@ element settings { element space { (element type { xsd:string { maxLength = "16" } } | attribute type { xsd:string { maxLength = "16" } }) & - (element length { xsd:int } | attribute length { xsd:int })? & (element interpolation { xsd:string { maxLength = "10" } } | - attribute interplation { xsd:string { maxLength = "10" } })? & + attribute interpolation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? }? & element angle { (element type { xsd:string { maxLength = "16" } } | attribute type { xsd:string { maxLength = "16" } }) & - (element length { xsd:int } | attribute length { xsd:int })? & (element interpolation { xsd:string { maxLength = "10" } } | - attribute interplation { xsd:string { maxLength = "10" } })? & + attribute interpolation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? }? & element energy { (element type { xsd:string { maxLength = "16" } } | attribute type { xsd:string { maxLength = "16" } }) & - (element length { xsd:int } | attribute length { xsd:int })? & (element interpolation { xsd:string { maxLength = "10" } } | - attribute interplation { xsd:string { maxLength = "10" } })? & + attribute interpolation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? }? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 8dd1a6febc..3536f3bd35 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -289,16 +289,6 @@ - - - - - - - - - - @@ -306,7 +296,7 @@ 10 - + 10 @@ -349,16 +339,6 @@ - - - - - - - - - - @@ -366,7 +346,7 @@ 10 - + 10 @@ -409,16 +389,6 @@ - - - - - - - - - - @@ -426,7 +396,7 @@ 10 - + 10 diff --git a/src/source.F90 b/src/source.F90 index a16eb245d0..c0181e844e 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -2,6 +2,7 @@ module source use bank_header, only: Bank use constants + use distribution_header, only: Delta use error, only: fatal_error use geometry, only: find_cell use geometry_header, only: BASE_UNIVERSE @@ -100,12 +101,8 @@ contains integer :: i ! dummy loop index real(8) :: r(3) ! sampled coordinates - real(8) :: phi ! azimuthal angle - real(8) :: mu ! cosine of polar angle real(8) :: p_min(3) ! minimum coordinates of source real(8) :: p_max(3) ! maximum coordinates of source - real(8) :: a ! Arbitrary parameter 'a' - real(8) :: b ! Arbitrary parameter 'b' logical :: found ! Does the source particle exist within geometry? type(Particle) :: p ! Temporary particle for using find_cell integer, save :: num_resamples = 0 ! Number of resamples encountered @@ -188,58 +185,28 @@ contains end select ! Sample angle - select case (external_source%type_angle) - case (SRC_ANGLE_ISOTROPIC) - ! Sample isotropic distribution - phi = TWO*PI*prn() - mu = TWO*prn() - ONE - site%uvw(1) = mu - site%uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - site%uvw(3) = sqrt(ONE - mu*mu) * sin(phi) + site%uvw(:) = external_source%angle%sample() - case (SRC_ANGLE_MONO) - ! Monodirectional source - site%uvw = external_source%params_angle - - case default - call fatal_error("No angle distribution specified for external source!") - end select - - ! Sample energy distribution - select case (external_source%type_energy) - case (SRC_ENERGY_MONO) - ! Monoenergtic source - site%E = external_source%params_energy(1) - if (site%E >= energy_max_neutron) then + ! Check for monoenergetic source above maximum neutron energy + select type (energy => external_source%energy) + type is (Delta) + if (energy%x0 >= energy_max_neutron) then call fatal_error("Source energy above range of energies of at least & &one cross section table") end if - - case (SRC_ENERGY_MAXWELL) - a = external_source%params_energy(1) - do - ! Sample Maxwellian fission spectrum - site%E = maxwell_spectrum(a) - - ! resample if energy is greater than maximum neutron energy - if (site%E < energy_max_neutron) exit - end do - - case (SRC_ENERGY_WATT) - a = external_source%params_energy(1) - b = external_source%params_energy(2) - do - ! Sample Watt fission spectrum - site%E = watt_spectrum(a, b) - - ! resample if energy is greater than maximum neutron energy - if (site%E < energy_max_neutron) exit - end do - - case default - call fatal_error("No energy distribution specified for external source!") end select + do + ! Sample energy spectrum + site%E = external_source%energy%sample() + + ! resample if energy is greater than maximum neutron energy + if (site%E < energy_max_neutron) exit + end do + + ! Set delayed group + site%delayed_group = 0 + ! Set the random number generator back to the tracking stream. call prn_set_stream(STREAM_TRACKING) diff --git a/src/source_header.F90 b/src/source_header.F90 index 11c74d110c..767b1bb654 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -1,5 +1,7 @@ module source_header + use distribution_header, only: Distribution, AngleDistribution + implicit none !=============================================================================== @@ -9,11 +11,9 @@ module source_header type ExtSource integer :: type_space ! spacial distribution, e.g. 'box' or 'point' - integer :: type_angle ! angle distribution, e.g. 'isotropic' - integer :: type_energy ! energy distribution, e.g. 'Watt' real(8), allocatable :: params_space(:) ! parameters for spatial distribution - real(8), allocatable :: params_angle(:) ! parameters for angle distribution - real(8), allocatable :: params_energy(:) ! parameters for energy distribution + type(AngleDistribution) :: angle ! angle distribution + class(Distribution), allocatable :: energy ! energy distribution end type ExtSource end module source_header diff --git a/src/state_point.F90 b/src/state_point.F90 index 046c857fd2..e89cdf12a8 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -22,6 +22,7 @@ module state_point use tally_header, only: TallyObject use mesh_header, only: RegularMesh use dict_header, only: ElemKeyValueII, ElemKeyValueCI + use random_lcg, only: seed #ifdef MPI use message_passing