diff --git a/Tallies/TallyMaps/Maps1D/Tests/directionMap_test.f90 b/Tallies/TallyMaps/Maps1D/Tests/directionMap_test.f90 index ff5c3e508..f09c43e9d 100644 --- a/Tallies/TallyMaps/Maps1D/Tests/directionMap_test.f90 +++ b/Tallies/TallyMaps/Maps1D/Tests/directionMap_test.f90 @@ -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() @@ -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 @@ -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 @@ -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()) @@ -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) diff --git a/Tallies/TallyMaps/Maps1D/Tests/radialMap_test.f90 b/Tallies/TallyMaps/Maps1D/Tests/radialMap_test.f90 index eea303f19..dedd16637 100644 --- a/Tallies/TallyMaps/Maps1D/Tests/radialMap_test.f90 +++ b/Tallies/TallyMaps/Maps1D/Tests/radialMap_test.f90 @@ -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]) @@ -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) diff --git a/Tallies/TallyMaps/Maps1D/directionMap_class.f90 b/Tallies/TallyMaps/Maps1D/directionMap_class.f90 index c6d345b37..a1ce82ed9 100644 --- a/Tallies/TallyMaps/Maps1D/directionMap_class.f90 +++ b/Tallies/TallyMaps/Maps1D/directionMap_class.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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() diff --git a/Tallies/TallyMaps/Maps1D/radialMap_class.f90 b/Tallies/TallyMaps/Maps1D/radialMap_class.f90 index 710b6fd50..52aa5f122 100644 --- a/Tallies/TallyMaps/Maps1D/radialMap_class.f90 +++ b/Tallies/TallyMaps/Maps1D/radialMap_class.f90 @@ -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 @@ -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 @@ -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') @@ -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) diff --git a/Tallies/TallyMaps/Tests/cylindricalMap_test.f90 b/Tallies/TallyMaps/Tests/cylindricalMap_test.f90 index 71aafb283..d73ec4dd4 100644 --- a/Tallies/TallyMaps/Tests/cylindricalMap_test.f90 +++ b/Tallies/TallyMaps/Tests/cylindricalMap_test.f90 @@ -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]) @@ -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) diff --git a/Tallies/TallyMaps/cylindricalMap_class.f90 b/Tallies/TallyMaps/cylindricalMap_class.f90 index f55d1e33f..7d326b64d 100644 --- a/Tallies/TallyMaps/cylindricalMap_class.f90 +++ b/Tallies/TallyMaps/cylindricalMap_class.f90 @@ -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; @@ -66,6 +66,7 @@ module cylindricalMap_class procedure :: map procedure :: print procedure :: kill + end type cylindricalMap contains @@ -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') @@ -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') diff --git a/docs/User Manual.rst b/docs/User Manual.rst index ac7575ed0..fba5d5a6f 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -1123,6 +1123,27 @@ Example: :: map { type cellMap; cells (1 5 3 2 4 100); undefBin T; } +* collNumMap (1D map), filters the particles tallied over number of collisions they underwent + + - collNumbers: list of collision numbers (integers) to be used as map bins + +Examples: :: + + map1 { type collNumMap; collNumbers ( 0 1 2 3 4 5 10 20); } + +* directionMap (1D map), angular map for the particle's direction with a linear grid + + - axis (*optional*, default = ``xy``): ``xy``, ``yz``, ``xz`` define the plane + of the direction to map + - N: number of bins + - min (*optional*, default = -180): grid lower limit [degrees] + - max (*optional*, default = 180): grid upper limit [degrees] + +Examples: :: + + map1 { type directionMap; axis xz; min 0.0; max 90.0; N 6; } + map2 { type directionMap; N 36; } + * energyMap (1D map), defines an energy group structure - grid: ``log`` for logarithmically spaced bins or ``lin`` for linearly spaced bins @@ -1177,23 +1198,33 @@ Example: :: map { type materialMap; materials (fuel water cladding reflector fuelGd); undefBin T; } -* multiMap, ensemble of multiple 1D maps +* radialMap, spherical or cylindrical radial map - - maps: list of the names of the maps that will compose the ``multiMap``. This - is followed by dictionaries that define the requested maps + - axis (*optional*, default = ``xyz``): ``x``, ``y``, ``z``, is the normal of + the cylindrical plane, or ``xyz`` to indicate spherical coordinates + - origin (*optional*, default = (0.0 0.0 0.0)): (x y z) vector with the origin + of the radial map. If the map is cylindrical, only the two coordinates perpendicular + to the cylinder's normal matter + - grid: ``lin`` for linearly spaced bins or ``equivolume`` -Example: :: + + min (*optional*, default = 0.0): minimum radius [cm] + + max: maximum radius [cm] + + N: number of radial bins - map { type multiMap; maps (map1 map2 map10); - map1 { <1D map definition> } - map2 { <1D map definition> } - map10 { <1D map definition> } - } + - grid: ``unstruct`` for unstructured grids, to be manually defined + + + bins: array with the explicit definition of the spherical bin boundaries + to be used + +Examples: :: + + map1 { type radialMap; axis xyz; origin (2.0 1.0 0.0); grid lin; min 3.0; max 10.0; N 14; } + map2 { type radialMap; axis z; grid equivolume; max 20.0; N 10; } + map3 { type radialMap; grid unstruct; bins (1.0 2.0 2.5 3.0 5.0); } * spaceMap (1D map), geometric cartesian map - axis: ``x``, ``y`` or ``z`` - - grid: ``lin`` for linearly spaced bins + min: bottom coordinate [cm] @@ -1209,27 +1240,23 @@ Examples: :: map1 { type spaceMap; axis x; grid lin; min -50.0; max 50.0; N 100; } map2 { type spaceMap; axis z; grid unstruct; bins (0.0 0.2 0.3 0.5 0.7 0.8 1.0); } -* sphericalMap, geometric spherical map - - - origin (*optional*, default = (0.0 0.0 .0.)): (x y z) vector with the origin - of the spherical map +* weightMap (1D map), divides weight into number of discrete bins - - grid: ``lin`` for linearly spaced bins or ``equivolume`` for spherical shells + - grid: ``log`` for logarithmically spaced bins or ``lin`` for linearly spaced bins - + Rmin (*optional*, default = 0.0): minimum radius [cm] - + Rmax: maximum radius [cm] - + N: number of radial bins + + min: bottom weight + + max: top weight + + N: number of bins - grid: ``unstruct`` for unstructured grids, to be manually defined - + bins: array with the explicit definition of the spherical bin boundaries - to be used + + bins: array with the explicit definition of the weight bin boundaries to be used Examples: :: - map1 { type sphericalMap; origin (2.0 1.0 0.0); grid lin; Rmin 3.0; Rmax 10.0; N 14; } - map2 { type sphericalMap; grid equivolume; Rmax 20.0; N 10; } - map3 { type sphericalMap; grid unstruct; bins (1.0 2.0 2.5 3.0 5.0); } + map1 { type weightMap; grid log; min 1.0e-3; max 100.0; N 100; } + map2 { type weightMap; grid lin; min 0.1; max 2.0; N 20; } + map3 { type weightMap; bins (0.0 0.2 0.4 0.6 0.8 1.0 2.0 5.0 10.0); } * cylindricalMap, geometric cylindrical map; other than the radial discretisation, one could add axial and azimuthal discretisation @@ -1261,31 +1288,18 @@ Example: :: map1 { type cylindricalMap; orientation y; origin (7.0 0.0); rGrid lin; Rmax 5.0; rN 10; } map2 { type cylindricalMap; rGrid unstruct; bins (2.0 3.0 4.5 5.0); axGrid lin; axMin 0.0; axMax 6.0 axN 24; azimuthalN 8; } -* collNumMap (1D map), filters the particles tallied over number of collisions they underwent - - - collNumbers: list of collision numbers (integers) to be used as map bins - -Examples: :: - - map1 { type collNumMap; collNumbers ( 0 1 2 3 4 5 10 20); } - -* weightMap (1D map), divides weight into number of discrete bins - - - grid: ``log`` for logarithmically spaced bins or ``lin`` for linearly spaced bins - - + min: bottom weight - + max: top weight - + N: number of bins - - - grid: ``unstruct`` for unstructured grids, to be manually defined +* multiMap, ensemble of multiple 1D maps - + bins: array with the explicit definition of the weight bin boundaries to be used + - maps: list of the names of the maps that will compose the ``multiMap``. This + is followed by dictionaries that define the requested maps -Examples: :: +Example: :: - map1 { type weightMap; grid log; min 1.0e-3; max 100.0; N 100; } - map2 { type weightMap; grid lin; min 0.1; max 2.0; N 20; } - map3 { type weightMap; bins (0.0 0.2 0.4 0.6 0.8 1.0 2.0 5.0 10.0); } + map { type multiMap; maps (map1 map2 map10); + map1 { <1D map definition> } + map2 { <1D map definition> } + map10 { <1D map definition> } + } Tally Filters #############