Skip to content

Commit

Permalink
Updating the directionMap to degrees, and updating the docs
Browse files Browse the repository at this point in the history
  • Loading branch information
V. Raffuzzi(vr339) committed May 5, 2024
1 parent 98ec4df commit 6989dfb
Show file tree
Hide file tree
Showing 7 changed files with 143 additions and 95 deletions.
26 changes: 14 additions & 12 deletions Tallies/TallyMaps/Maps1D/Tests/directionMap_test.f90
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,10 @@ subroutine setUp(this)
call tempDict % kill()

! Build map with default plane
call tempDict % init(1)
call tempDict % store('N', 8)
call tempDict % init(3)
call tempDict % store('N', 9)
call tempDict % store('min', 90)
call tempDict % store('max', 180)

call this % map2 % init(tempDict)
call tempDict % kill()
Expand Down Expand Up @@ -70,8 +72,8 @@ end subroutine tearDown
subroutine testMap1(this)
class(test_directionMap), intent(inout) :: this
real(defReal),dimension(4),parameter :: x = [0.44_defReal, 15.8_defReal, 83.2_defReal, 999.1_defReal]
real(defReal),dimension(4),parameter :: phi = [-0.1_defReal, 0.84_defReal, 0.33_defReal, 0.64_defReal]*PI
integer(shortInt),dimension(4),parameter :: RES_IDX = [2, 4, 3, 4]
real(defReal),dimension(4),parameter :: phi = [100.1_defReal, 20.84_defReal, 333.9_defReal, 264.8_defReal]*PI/180.0_defReal
integer(shortInt),dimension(4),parameter :: RES_IDX = [4, 3, 2, 1]
integer(shortInt),dimension(4) :: idx
type(particleState),dimension(4) :: states

Expand All @@ -92,19 +94,19 @@ end subroutine testMap1
@Test
subroutine testMap2(this)
class(test_directionMap), intent(inout) :: this
real(defReal),dimension(4),parameter :: x = [0.9999_defReal, 0.87_defReal, -0.3_defReal, 0.18_defReal]
real(defReal),dimension(4),parameter :: y = [0.45_defReal, 0.88_defReal, 0.51_defReal, -0.92_defReal]
real(defReal),dimension(4),parameter :: z = [1.0_defReal, 39.8_defReal, 0.05_defReal, -12.2_defReal]
integer(shortInt),dimension(4),parameter :: RES_IDX = [5, 6, 7, 3]
real(defReal),dimension(4),parameter :: z = [0.44_defReal, 15.8_defReal, 83.2_defReal, 999.1_defReal]
real(defReal),dimension(4),parameter :: phi = [170.1_defReal, 90.84_defReal, 133.9_defReal, 264.8_defReal]*PI/180.0_defReal
integer(shortInt),dimension(4),parameter :: RES_IDX = [9, 1, 5, 0]
integer(shortInt),dimension(4) :: idx
type(particleState),dimension(4) :: states

! Initialise states
states(:) % dir(1) = x
states(:) % dir(2) = y
states(:) % dir(1) = cos(phi)
states(:) % dir(2) = sin(phi)
states(:) % dir(3) = z

idx = this % map2 % map(states)

@assertEqual(RES_IDX, idx)

end subroutine testMap2
Expand All @@ -120,7 +122,7 @@ subroutine testBinNumber(this)
@assertEqual(4, this % map1 % bins(0),'All bins')
@assertEqual(4, this % map1 % bins(1),'1st dimension')
@assertEqual(0, this % map2 % bins(2),'2nd dimension')
@assertEqual(8, this % map2 % bins(1),'1st dimension')
@assertEqual(9, this % map2 % bins(1),'1st dimension')

! Get dimensionality
@assertEqual(1, this % map1 % dimensions())
Expand All @@ -130,7 +132,7 @@ end subroutine testBinNumber

!!
!! Test correctness of print subroutine
!! Does not checks that values are correct, but that calls sequence is without errors
!! Does not check that values are correct, but that call sequence is without errors
!!
@Test
subroutine testPrint(this)
Expand Down
4 changes: 2 additions & 2 deletions Tallies/TallyMaps/Maps1D/Tests/radialMap_test.f90
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ subroutine setUp(this)
! Build cylindrical map with different origin & unstruct bins
call tempDict % init(4)
call tempDict % store('axis','z')
call tempDict % store('origin',[ONE, ONE])
call tempDict % store('origin',[ONE, ONE, ZERO])
call tempDict % store('grid','unstruct')
call tempDict % store('bins', [1.5_defReal, 2.3_defReal, 3.8_defReal, 8.0_defReal])

Expand Down Expand Up @@ -292,7 +292,7 @@ end subroutine testBinNumber

!!
!! Test correctness of print subroutine
!! Does not checks that values are correct, but that calls sequence is without errors
!! Does not check that values are correct, but that call sequence is without errors
!!
@Test
subroutine testPrint(this)
Expand Down
56 changes: 42 additions & 14 deletions Tallies/TallyMaps/Maps1D/directionMap_class.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,20 +13,29 @@ module directionMap_class
private

!!
!! Maps the particle direction in linear bins in the range 0-360 degrees
!! Maps the particle direction in linear bins in the range -180 -> 180 degrees (default)
!! or boundaries defined by the user
!!
!! The angle is calculated using as a reference the positive direction of the
!! axis x in the 'xy' and 'xz' cases, or y in the 'yz' case
!!
!! NOTE: the map is built in radians for convenience when using the function atan2, but
!! the user input and the results are in degrees for easy of interpretation
!!
!! Interface:
!! tallyMap1D interface
!!
!! NOTE:
!! Behaviour of points exactly at the boundary of bins is undefined.
!! particle can end-up in either of the two
!! Particle can end-up in either of the two
!!
!! Sample Dictionary Input:
!! directionMap {
!! type directionMap;
!! #plane xz;# // Optional. Default xy
!! #plane xz;# // Optional. Default xy
!! N 10;
!! #min 60;# // Optional. Default -180
!! #max 120;# // Optional. Default 180
!! }
!!
type, public, extends (tallyMap1D) :: directionMap
Expand Down Expand Up @@ -55,9 +64,10 @@ module directionMap_class
!! See tallyMap for specification.
!!
subroutine init(self, dict)
class(directionMap), intent(inout) :: self
class(dictionary), intent(in) :: dict
character(nameLen) :: type
class(directionMap), intent(inout) :: self
class(dictionary), intent(in) :: dict
character(nameLen) :: type
real(defReal) :: min, max
character(100), parameter :: Here = 'init (directionMap_class.f90)'

! Check orientation of the cylinder
Expand All @@ -81,13 +91,28 @@ subroutine init(self, dict)
self % DIM2 = Y_AXIS

case default
call fatalError(Here, 'Keyword orientation must be x, y or z. It is: '//type)
call fatalError(Here, 'Keyword orientation must be xy, yz or xz. It is: '//type)

end select

! Build grid
! Get grid details
call dict % get(self % N, 'N')
call self % bounds % init(-PI, PI, self % N, 'lin')
call dict % getOrDefault(min, 'min', -180.0_defReal)
call dict % getOrDefault(max, 'max', 180.0_defReal)

! Check boundaries
if (min < -180.0_defReal .or. min > 180.0_defReal) then
call fatalError(Here, 'Minimum angle must be between -180 and 180 degrees. It is: '//type)
end if

if (max < -180.0_defReal .or. max > 180.0_defReal) then
call fatalError(Here, 'Maximum angle must be between -180 and 180 degrees. It is: '//type)
end if

! Build map in radians
min = min*PI/180.0_defReal
max = max*PI/180.0_defReal
call self % bounds % init(min, max, self % N, 'lin')

end subroutine init

Expand Down Expand Up @@ -133,11 +158,14 @@ elemental function map(self, state) result(idx)
integer(shortInt) :: idx
real(defReal) :: x, y, theta

! Map the angle
x = state % dir(self % DIM1)
y = state % dir(self % DIM2)

! Returns angle in radians in the range -PI to PI
theta = atan2(y,x)
! Search along the azimuthal dimension and return 0 if index is out-of-bounds

! Search in the grid return 0 if index is out-of-bounds
idx = self % bounds % search(theta)

! Should never happen
Expand Down Expand Up @@ -165,11 +193,11 @@ subroutine print(self,out)
call out % startArray(name, [2,self % N])

do i = 1, self % N
! Print lower bin boundary
call out % addValue(self % bounds % bin(i))
! Print lower bin boundary in degrees
call out % addValue(self % bounds % bin(i)/PI*180.0_defReal)

! Print upper bin boundar
call out % addValue(self % bounds % bin(i + 1))
! Print upper bin boundary in degrees
call out % addValue(self % bounds % bin(i + 1)/PI*180.0_defReal)
end do

call out % endArray()
Expand Down
24 changes: 13 additions & 11 deletions Tallies/TallyMaps/Maps1D/radialMap_class.f90
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ module radialMap_class
!! Sample Dictionary Input:
!! radialMap {
!! type radialMap;
!! axis x; // Optional. Default is spherical map
!! axis x; // Optional. Default is 'xyz' (spherical map)
!! #origin (1.0 0.0 0.0);# // Optional. Default (0.0 0.0 0.0)
!! grid lin;
!! #min 2.0;# // Optional. Default 0.0
Expand Down Expand Up @@ -61,7 +61,7 @@ module radialMap_class
subroutine init(self, dict)
class(radialMap), intent(inout) :: self
class(dictionary), intent(in) :: dict
real(defReal), dimension(:), allocatable :: grid
real(defReal), dimension(:), allocatable :: temp, grid
character(nameLen) :: type
real(defReal) :: min, max, vol, exp
integer(shortInt) :: i
Expand Down Expand Up @@ -91,17 +91,19 @@ subroutine init(self, dict)
spherical = .true.

case default
call fatalError(Here, 'Keyword orientation must be x, y, z or nothing for spherical. It is: '//type)
call fatalError(Here, 'Keyword orientation must be x, y, z or xyz for spherical. It is: '//type)

end select

! Check & load origin
if (spherical) then
call dict % getOrDefault(self % origin, 'origin', [ZERO, ZERO, ZERO])
else
call dict % getOrDefault(self % origin, 'origin', [ZERO, ZERO])
call dict % getOrDefault(temp, 'origin', [ZERO, ZERO, ZERO])

if (size(temp) /= 3) then
call fatalError(Here, 'Expected 3 values for origin. Got: ' // numToChar(size(temp)))
end if

self % origin = temp(self % axis)

! Load radial grid information
if (.not. dict % isPresent('grid')) call fatalError(Here, 'Keyword grid must be present')
call dict % get(type, 'grid')
Expand Down Expand Up @@ -208,10 +210,10 @@ end function getAxisName
!! See tallyMap for specification.
!!
elemental function map(self, state) result(idx)
class(radialMap), intent(in) :: self
class(particleState), intent(in) :: state
integer(shortInt) :: idx
real(defReal) :: r
class(radialMap), intent(in) :: self
class(particleState), intent(in) :: state
integer(shortInt) :: idx
real(defReal) :: r

! Calculate the distance from the origin
r = norm2(state % r(self % axis) - self % origin)
Expand Down
4 changes: 2 additions & 2 deletions Tallies/TallyMaps/Tests/cylindricalMap_test.f90
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ subroutine setUp(this)

! Build map with different origin & unstruct bins
call tempDict % init(3)
call tempDict % store('origin',[ONE, ONE])
call tempDict % store('origin',[ONE, ONE, TWO])
call tempDict % store('rGrid','unstruct')
call tempDict % store('bins', [1.5_defReal, 2.3_defReal, 3.8_defReal, 8.0_defReal])

Expand Down Expand Up @@ -179,7 +179,7 @@ end subroutine testBinNumber

!!
!! Test correctness of print subroutine
!! Does not checks that values are correct, but that calls sequence is without errors
!! Does not check that values are correct, but that call sequence is without errors
!!
@Test
subroutine testPrint(this)
Expand Down
20 changes: 11 additions & 9 deletions Tallies/TallyMaps/cylindricalMap_class.f90
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ module cylindricalMap_class
!! cylindricalMap {
!! type cylindricalMap;
!! #orientation x;# // Optional. Defeult z
!! #origin (1.0 0.0);# // Optional. Default (0 0). Order is xy, xz or yz.
!! #origin (1.0 0.0 0.0);# // Optional. Default (0 0 0)
!! rGrid lin;
!! #Rmin 2.0;# // Optional. Default 0.0
!! Rmax 10.0;
Expand Down Expand Up @@ -66,6 +66,7 @@ module cylindricalMap_class
procedure :: map
procedure :: print
procedure :: kill

end type cylindricalMap

contains
Expand All @@ -84,14 +85,6 @@ subroutine init(self, dict)
integer(shortInt) :: i
character(100), parameter :: Here = 'init (cylindricalMap_class.f90)'

! Check & load origin
call dict % getOrDefault(temp, 'origin', [ZERO, ZERO])

if (size(temp) /= 2) then
call fatalError(Here, 'Expected 2 values for origin. Got: ' // numToChar(size(temp)))
end if
self % origin = temp

! Check orientation of the cylinder
if (dict % isPresent('orientation')) then
call dict % get(type, 'orientation')
Expand Down Expand Up @@ -120,6 +113,15 @@ subroutine init(self, dict)

end select

! Check & load origin
call dict % getOrDefault(temp, 'origin', [ZERO, ZERO, ZERO])

if (size(temp) /= 3) then
call fatalError(Here, 'Expected 3 values for origin. Got: ' // numToChar(size(temp)))
end if

self % origin = temp([self % DIM1, self % DIM2])

! Load radial grid information
if (.not. dict % isPresent('rGrid')) call fatalError(Here, 'Keyword rGrid must be present')
call dict % get(type, 'rGrid')
Expand Down
Loading

0 comments on commit 6989dfb

Please sign in to comment.