diff options
Diffstat (limited to 'Dragon/src/SNTSFH.f')
| -rw-r--r-- | Dragon/src/SNTSFH.f | 253 |
1 files changed, 253 insertions, 0 deletions
diff --git a/Dragon/src/SNTSFH.f b/Dragon/src/SNTSFH.f new file mode 100644 index 0000000..630f38a --- /dev/null +++ b/Dragon/src/SNTSFH.f @@ -0,0 +1,253 @@ +*DECK SNTSFH + SUBROUTINE SNTSFH (IMPX,ITYPE,NHEX,LZ,MCELL,ISPLH,MAT,LOZSWP, + > COORDMAP) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Output arrays for lozenge sweep order (direction-dependent) and +* coordinate map, both needed for resolution of the discrete ordinates +* transport equation in hexagonal geometry. +* +*Copyright: +* Copyright (C) 2025 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. A. Calloo +* +*Parameters: input +* IMPX print parameter. +* ITYPE geometry type (8:hexagonal 2D, 9:hexagonal 3D). +* NHEX number of hexagons (for 3D, in one plane only). +* LZ number of mesh elements in z-axis (including split). +* MCELL number of macrocells to use along z-axis. +* ISPLH mesh-splitting in 3*ISPLH**2 lozenges per hexagon. +* MAT mixture index assigned to each element. +* +*Parameters: local +* NRINGS number of hexagonal rings in the domain, assuming the centre +* hexagon counts as 1 ring. +* +*Parameters: output +* LOZSWP lozenge sweep order depending on direction. +* COORDMAP coordinate map: mapping hexagon from DRAGON geometry indices +* to the axial coordinate system, using p, r, s axes. The s +* axis is redundant, which means that using p and r axes +* effectively maps the hexagon geometry to a 2D map. Refer to +* the redblobgames blog for more information. +* +*Comments: +* The lozenge under consideration is given by the position within the +* matrix. See user manual and/or data manual and/or thesis +* _____ +* / / \ +* / B / \ +* ,----(---- A )----. +* / \ C \ / \ +* / \___\_/ \ +* \ 4 / \ 2 / +* \ / \ / +* )----( 1 )----( +* / \ / \ +* / \_____/ \ +* \ 5 / \ 7 / +* \ / \ / +* `----( 6 )----' +* \ / +* \_____/ +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IMPX,ITYPE,NHEX,LZ,MCELL,ISPLH,MAT(ISPLH**2,3,NHEX,LZ), + > COORDMAP(3,NHEX) +*---- +* LOCAL VARIABLES +*---- + INTEGER, DIMENSION(3,6) :: LOZSWP,MAPCODE + INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: TMPMAT + INTEGER, ALLOCATABLE, DIMENSION(:) :: TASKSPERWAVE +*---- +* LOZENGE SWEEP ORDERING WITHIN HEXAGONS - DIRECTION DEPENDENT +*---- + LOZSWP = RESHAPE((/ 3, 2, 1, 3, 1, 2, 1, 3, 2, 1, 2, 3, 2, 1, + > 3, 2, 3, 1 /), SHAPE(LOZSWP)) +*---- +* CREATE COORDIDATE MAP FROM DRAGON INDEX TO AXIAL COORDINATES +*---- + NRINGS=INT((SQRT( REAL((4*NHEX-1)/3) )+1.)/2.) + IF(NRINGS.EQ.1) CALL XABORT('NOT IMPLEMENTED FOR SINGLE HEX YET.') + IF(NHEX.NE.1+3*NRINGS*(NRINGS-1)) CALL XABORT('SNTSFH: INVALID ' + 1 //'VALUE OF NHEX(1).') +* + MAPCODE = RESHAPE((/ -1, 0, 1, -1, 1, 0, 0, 1, -1, 1, 0, -1, 1, + > -1, 0, 0, -1, 1 /), SHAPE(MAPCODE)) +* + ! It should be noted that the algorithm below effectively + ! reverses the y-axis. However, this should be of no consequence + ! whatsoever as it would akin to the user defining the domain + ! somewhat differently in the geometry. Calculations and results + ! should be unaffected. + IHEX_DOM=1 + DO IRING=1,NRINGS + ! Initialise first 'ring', i.e., centre hexagon + IF(IRING.EQ.1) THEN + ITMP1 = NRINGS + ITMP2 = NRINGS + ITMP3 = -2*(NRINGS) + COORDMAP(1,IHEX_DOM)=ITMP1 + COORDMAP(2,IHEX_DOM)=ITMP2 + COORDMAP(3,IHEX_DOM)=ITMP3 + + IHEX_DOM = IHEX_DOM+1 + ! Ignore rest of this loop and move on to next ring + CYCLE + ENDIF + + ! Find coordinates for hexagon when moving from ring n-1 to n + ITMP1 = ITMP1+1 + ITMP2 = ITMP2-1 + ITMP3 = ITMP3+0 + COORDMAP(1,IHEX_DOM)=ITMP1 + COORDMAP(2,IHEX_DOM)=ITMP2 + COORDMAP(3,IHEX_DOM)=ITMP3 + IHEX_DOM = IHEX_DOM+1 + + ! 'Sweep' each of the 3 axes of the hexagonal plane and their + ! negative directions + DO IND=1,6 + ! Step through each hexagon per each axis + DO IHEX=1,IRING-1 + ITMP1 = ITMP1+MAPCODE(1,IND) + ITMP2 = ITMP2+MAPCODE(2,IND) + ITMP3 = ITMP3+MAPCODE(3,IND) + ! Store each of the coordinates except the last hexagon + ! in the last direction. This is because we already + ! computed that hexagon when moving from ring n-1 to n + IF((IND.EQ.6).AND.(IHEX.EQ.IRING-1))THEN + CONTINUE + ELSE + COORDMAP(1,IHEX_DOM)=ITMP1 + COORDMAP(2,IHEX_DOM)=ITMP2 + COORDMAP(3,IHEX_DOM)=ITMP3 + IHEX_DOM = IHEX_DOM+1 + ENDIF + ENDDO ! ihex + ENDDO ! ind + ENDDO ! iring + +*---- +* COMPUTE NUMBER OF CONCURRENT HEXAGONS PER WAVEFRONT FOR PRINTING +* PURPOSES ONLY +*---- + IF(MCELL > 0)THEN + + ! Build material array in axial coordinates + NCOLS=2*NRINGS -1 + ALLOCATE(TMPMAT(ISPLH**2,3,NCOLS,NCOLS,LZ)) + TMPMAT(:,:,:,:,:)=-1 + DO IZ=1,LZ + DO IHEX_XY=1,NHEX + TMPMAT(:,:,COORDMAP(1,IHEX_XY),COORDMAP(2,IHEX_XY),IZ) = + > MAT(:,:,IHEX_XY,IZ) + ENDDO + ENDDO + + ! Build TasksPerWave array + IF(ITYPE==8)THEN + ! 2D Hexagonal + NWAVES=NCOLS+NCOLS-1 + ALLOCATE(TASKSPERWAVE(NWAVES)) + TASKSPERWAVE(:)=0 + + DO IWAVE=1,NWAVES + ICOUNT = 0 + DO J=MAX(1,IWAVE-NCOLS+1),MIN(NCOLS,IWAVE) + I=IWAVE-J+1 + I=NCOLS+1-I + IF((I.GT.NCOLS).OR.(I.LT.1)) CYCLE + IF((J.GT.NCOLS).OR.(J.LT.1)) CYCLE + ! If within corners of Cartesian axial coordinate map + ! (where there are no hexagons), skip loop + IF(TMPMAT(1,1,I,J,1).EQ.-1) CYCLE + ICOUNT = ICOUNT + 1 + ENDDO + TASKSPERWAVE(IWAVE) = ICOUNT + ENDDO + ELSEIF(ITYPE==9)THEN + ! 3D Hexagonal + MCELLZ = MCELL + NWAVES=NCOLS+NCOLS+MCELLZ-2 + ALLOCATE(TASKSPERWAVE(NWAVES)) + TASKSPERWAVE(:)=0 + + DO IWAVE=1,NWAVES + ICOUNT = 0 + J_STT=MAX(1,IWAVE-NCOLS-MCELLZ+2) + J_END=MIN(NCOLS,IWAVE) + DO J_MC=J_STT,J_END + J=J_MC + I_STT=MAX(1,IWAVE-J_MC-MCELLZ+2) + I_END=MIN(NCOLS,IWAVE-J_MC+1) + DO I_MC=I_STT,I_END + I=I_MC + I=NCOLS+1-I + ! If within corners of Cartesian axial coordinate map + ! (where there are no hexagons), skip loop + IF(TMPMAT(1,1,I,J,1).EQ.-1) CYCLE + K_MC=IWAVE-I_MC-J_MC+2 + ICOUNT = ICOUNT + 1 + ENDDO + ENDDO + TASKSPERWAVE(IWAVE) = ICOUNT + ENDDO + ENDIF + + DEALLOCATE(TMPMAT) + + ENDIF + +*---- +* PRINT A FEW GEOMETRY CHARACTERISTICS +*---- + IF(IMPX.GT.2) THEN + WRITE(*, 100) + WRITE(*, 101) NCOLS + WRITE(*, 102) NRINGS + WRITE(*, 103) + DO I=1,6 + WRITE(*,104) I, LOZSWP(:,I) + ENDDO + IF(MCELL > 0)THEN + WRITE(*, 105) NWAVES + WRITE(*, 106) + DO I = 1, NWAVES + WRITE(*, 107) TASKSPERWAVE(I) + END DO + DEALLOCATE(TASKSPERWAVE) + ENDIF + ENDIF + IF(IMPX.GT.4) THEN + WRITE(*, 109) + WRITE(*, 110) + DO I = 1, NHEX + WRITE(*, 111) I, COORDMAP(:, I) + END DO + ENDIF + + RETURN + 100 FORMAT (' ') + 101 FORMAT ('NCOLS =', I4) + 102 FORMAT ('NRINGS =', I4) + 103 FORMAT ('LOZENGE SWEEP ORDER') + 104 FORMAT ('IND_XY:', I4, ' LOZ. ORDER:', 3I4) + 105 FORMAT ('NWAVES =', I4) + 106 FORMAT ('TASKS PER WAVE') + 107 FORMAT (I4) + 109 FORMAT (' ') + 110 FORMAT ('COORDINATE MAP IS GIVEN BELOW:') + 111 FORMAT ('DRAGON IND:', I4, ' AXIAL COORD:', 3I4) + END |