Split distribution_header into distribution_univariate and

distribution_multivariate and extend capabilities for sampling on unit sphere.

src/distribution_multivariate.F90

@@ -0,0 +1,93 @@
+module distribution_multivariate
+
+  use constants, only: ONE, TWO, PI
+  use distribution_univariate, only: Distribution
+  use math, only: rotate_angle
+  use random_lcg, only: prn
+
+!===============================================================================
+! UNITSPHEREDISTRIBUTION type defines a probability density function for points
+! on the unit sphere. Extensions of this type are used to sample angular
+! distributions for starting soures
+!===============================================================================
+
+  type, abstract :: UnitSphereDistribution
+    real(8) :: reference_uvw(3)
+  contains
+    procedure(iSample), deferred :: sample
+  end type UnitSphereDistribution
+
+  abstract interface
+    function iSample(this) result(uvw)
+      import UnitSphereDistribution
+      class(UnitSphereDistribution), intent(in) :: this
+      real(8) :: uvw(3)
+    end function iSample
+  end interface
+
+!===============================================================================
+! Derived classes of UnitSphereDistribution
+!===============================================================================
+
+  ! Explicit distribution of polar and azimuthal angles
+  type, extends(UnitSphereDistribution) :: PolarAzimuthal
+    class(Distribution), allocatable :: mu
+    class(Distribution), allocatable :: phi
+  contains
+    procedure :: sample => polar_azimuthal_sample
+  end type PolarAzimuthal
+
+  ! Uniform distribution on the unit sphere
+  type, extends(UnitSphereDistribution) :: Isotropic
+  contains
+    procedure :: sample => isotropic_sample
+  end type Isotropic
+
+  ! Monodirectional distribution
+  type, extends(UnitSphereDistribution) :: Monodirectional
+  contains
+    procedure :: sample => monodirectional_sample
+  end type Monodirectional
+
+contains
+
+  function polar_azimuthal_sample(this) result(uvw)
+    class(PolarAzimuthal), intent(in) :: this
+    real(8) :: uvw(3)
+
+    real(8) :: mu     ! cosine of polar angle
+    real(8) :: phi    ! azimuthal angle
+
+    ! Sample cosine of polar angle
+    mu = this%mu%sample()
+    if (mu == ONE) then
+      uvw(:) = this%reference_uvw
+    else
+      ! Sample azimuthal angle
+      phi = this%phi%sample()
+      uvw(:) = rotate_angle(this%reference_uvw, mu, phi)
+    end if
+  end function polar_azimuthal_sample
+
+  function isotropic_sample(this) result(uvw)
+    class(Isotropic), intent(in) :: this
+    real(8) :: uvw(3)
+
+    real(8) :: phi
+    real(8) :: mu
+
+    phi = TWO*PI*prn()
+    mu = TWO*prn() - ONE
+    uvw(1) = mu
+    uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
+    uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
+  end function isotropic_sample
+
+  function monodirectional_sample(this) result(uvw)
+    class(Monodirectional), intent(in) :: this
+    real(8) :: uvw(3)
+
+    uvw(:) = this%reference_uvw
+  end function monodirectional_sample
+
+end module distribution_multivariate

src/distribution_univariate.F90

@@ -1,8 +1,8 @@
-module distribution_header
+module distribution_univariate
 
-  use constants, only: ZERO, HALF, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
+  use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR
   use error, only: fatal_error
-  use math, only: rotate_angle, maxwell_spectrum, watt_spectrum
+  use math, only: maxwell_spectrum, watt_spectrum
   use random_lcg, only: prn
 
 !===============================================================================
@@ -14,6 +14,10 @@ module distribution_header
     procedure(iSample), deferred :: sample
   end type Distribution
 
+  type DistributionContainer
+    class(Distribution), allocatable :: obj
+  end type DistributionContainer
+
   abstract interface
     function iSample(this) result(x)
       import Distribution
@@ -61,23 +65,12 @@ module distribution_header
     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
+    real(8), allocatable :: 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)
@@ -199,28 +192,4 @@ contains
     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) :: phi
-    real(8) :: mu
-
-    select type (polar_cos => this%mu)
-    type is (Uniform)
-      phi = TWO*PI*prn()
-      mu = TWO*prn() - ONE
-      uvw(1) = mu
-      uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
-      uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
-    class default
-      mu = polar_cos%sample()
-      if (mu == ONE) then
-        uvw(:) = this%reference_uvw
-      else
-        uvw(:) = rotate_angle(this%reference_uvw, mu)
-      end if
-    end select
-  end function angle_sample
-
-end module distribution_header
+end module distribution_univariate

src/global.F90

@@ -452,8 +452,8 @@ contains
     ! Deallocate external source
     if (allocated(external_source % params_space)) &
          deallocate(external_source % params_space)
-    if (allocated(external_source % angle % mu)) &
-         deallocate(external_source % angle % mu)
+    if (allocated(external_source % angle)) &
+         deallocate(external_source % angle)
     if (allocated(external_source % energy)) &
          deallocate(external_source % energy)
 

src/input_xml.F90

@@ -3,7 +3,8 @@ module input_xml
   use cmfd_input,       only: configure_cmfd
   use constants
   use dict_header,      only: DictIntInt, ElemKeyValueCI
-  use distribution_header
+  use distribution_multivariate
+  use distribution_univariate
   use energy_grid,      only: grid_method, n_log_bins
   use error,            only: fatal_error, warning
   use geometry_header,  only: Cell, Lattice, RectLattice, HexLattice
@@ -438,17 +439,29 @@ contains
              call get_node_value(node_dist, "type", type)
         select case (to_lower(type))
         case ('isotropic')
-          allocate(Uniform :: external_source%angle%mu)
+          allocate(Isotropic :: external_source%angle)
 
         case ('monodirectional')
-          allocate(Delta :: external_source%angle%mu)
+          allocate(Monodirectional :: external_source%angle)
           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)
+          allocate(PolarAzimuthal :: external_source%angle)
+          select type (angle => external_source%angle)
+          type is (PolarAzimuthal)
+            allocate(Tabular :: angle%mu)
+
+            ! For now, azimuthal is uniform
+            allocate(Uniform :: angle%phi)
+            select type (phi => angle%phi)
+            type is (Uniform)
+              phi%a = ZERO
+              phi%b = TWO*PI
+            end select
+          end select
 
         case default
           call fatal_error("Invalid angular distribution for external source: "&
@@ -459,48 +472,41 @@ contains
         external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
 
         ! Read parameters for angle distribution
-        select type (mu => external_source%angle%mu)
-        type is (Uniform)
-          mu%a = -ONE
-          mu%b = ONE
-
-        type is (Delta)
-          mu%x0 = ONE
+        select type (angle => external_source%angle)
+        type is (Monodirectional)
           call get_node_array(node_dist, "parameters", &
                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')
+        type is (PolarAzimuthal)
+          select type (mu => angle%mu)
+          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
-            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
+            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)
+            ! 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
         end select
 
       else
         ! Set default angular distribution 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
+        allocate(Isotropic :: external_source%angle)
         external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
       end if
 

src/math.F90

@@ -563,12 +563,13 @@ contains
 ! with direct sampling rather than rejection as is done in MCNP and SERPENT.
 !===============================================================================
 
-  function rotate_angle(uvw0, mu) result(uvw)
+  function rotate_angle(uvw0, mu, phi) result(uvw)
     real(8), intent(in) :: uvw0(3) ! directional cosine
     real(8), intent(in) :: mu      ! cosine of angle in lab or CM
+    real(8), optional   :: phi     ! azimuthal angle
     real(8)             :: uvw(3)  ! rotated directional cosine
 
-    real(8) :: phi    ! azimuthal angle
+    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)
@@ -582,12 +583,16 @@ contains
     v0 = uvw0(2)
     w0 = uvw0(3)
 
-    ! Sample azimuthal angle in [0,2pi)
-    phi = TWO * PI * prn()
+    ! Sample azimuthal angle in [0,2pi) if none provided
+    if (present(phi)) then
+      phi_ = phi
+    else
+      phi_ = TWO * PI * prn()
+    end if
 
     ! Precompute factors to save flops
-    sinphi = sin(phi)
-    cosphi = cos(phi)
+    sinphi = sin(phi_)
+    cosphi = cos(phi_)
     a = sqrt(max(ZERO, ONE - mu*mu))
     b = sqrt(max(ZERO, ONE - w0*w0))
 

src/source.F90

@@ -2,7 +2,7 @@ module source
 
   use bank_header,      only: Bank
   use constants
-  use distribution_header, only: Delta
+  use distribution_univariate, only: Delta
   use error,            only: fatal_error
   use geometry,         only: find_cell
   use geometry_header,  only: BASE_UNIVERSE

src/source_header.F90

@@ -1,6 +1,7 @@
 module source_header
 
-  use distribution_header, only: Distribution, AngleDistribution
+  use distribution_univariate, only: Distribution
+  use distribution_multivariate, only: UnitSphereDistribution
 
   implicit none
 
@@ -12,7 +13,7 @@ module source_header
   type ExtSource
     integer :: type_space              ! spacial distribution, e.g. 'box' or 'point'
     real(8), allocatable :: params_space(:)  ! parameters for spatial distribution
-    type(AngleDistribution) :: angle           ! angle distribution
+    class(UnitSphereDistribution), allocatable :: angle  ! angle distribution
     class(Distribution), allocatable :: energy ! energy distribution
   end type ExtSource