Initial commit from Polytechnique Montreal

1 files changed, 169 insertions, 0 deletions

diff --git a/Dragon/src/SNEST.f b/Dragon/src/SNEST.f
new file mode 100644
index 0000000..0f3f52a
--- /dev/null
+++ b/Dragon/src/SNEST.f
@@ -0,0 +1,169 @@
+*DECK SNEST
+      SUBROUTINE SNEST (IPTRK,IMPX,NREG,NUN,MAT,IG,KEYFLX,KEYSPN,
+     1   FUNKNO) 
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Rearrange SPn flux in the Sn order so that SPn can be used to 
+* initialise Sn calculation. Use SPn flux to obtain rough estimate
+* of boundary fluxes.
+*
+*Copyright:
+* Copyright (C) 2020 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
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IMPX    print level.
+* NREG    total number of regions for which specific values of the
+*         neutron flux and reactions rates are required.
+* NUN     total number of unknowns in vectors SUNKNO and FUNKNO.
+* MAT     index-number of the mixture type assigned to each volume.
+* IG      group number.
+* KEYFLX  position of averaged flux elements in FUNKNO vector.
+* KEYSPN  position of SPn unknowns in FUNKNO vector.
+*
+*Parameters: input/output
+* FUNKNO  SPn (in) / SN (out) unknown vector.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK
+      INTEGER     IMPX,NREG,NUN,MAT(NREG),IG,KEYFLX(NREG),KEYSPN(NREG) 
+      REAL        FUNKNO(NUN)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (IUNOUT=6,NSTATE=40)
+      INTEGER ISTRK(NSTATE),NLOZH,SPLTL,SBMSH,REM,ISPLH
+      REAL ZCODE(6)
+*
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: TMPKEY,ORIKEY
+      REAL, ALLOCATABLE, DIMENSION(:) :: FUNSPN
+*
+*----
+*  RECOVER TRACKING INFORMATION
+*----
+      CALL LCMGET(IPTRK,'ZCODE',ZCODE)
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTRK)
+      ITYPE=ISTRK(6)
+      NSCT=ISTRK(7)
+      IELEM=ISTRK(8)
+      NDIM=ISTRK(9)
+      LL4=ISTRK(11)
+      LX=ISTRK(12)
+      LY=ISTRK(13)
+      LZ=ISTRK(14)
+      ISPLH=1
+      NHEX=1
+      IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9)) THEN
+         ISPLH=ISTRK(26)
+         NHEX =LX/(3*ISPLH**2)
+      ENDIF
+*----
+*  PARAMETER VALIDATION
+*---- 
+      IF((KEYSPN(NREG)).GT.LL4)THEN
+         CALL XABORT('SNEST: MORE SPN UNKNOWNS THAN SN UNKNOWNS. '
+     1   //'CANNOT GUARANTEE INTEGRITY OF IMPORTED SOLUTION. '
+     2   //'CONSIDER INCREASING SPATIAL ORDER FOR SN OR DECREASING '
+     3   //'SPATIAL ORDER FOR SPN/DIFF.')
+      ENDIF
+      IF((ITYPE.NE.2).AND.(ITYPE.NE.5).AND.(ITYPE.NE.7).AND.
+     1   (ITYPE.NE.8).AND.(ITYPE.NE.9))CALL XABORT('SNEST: TYPE '
+     2   //'OF DISCRETIZATION NOT IMPLEMENTED.')
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(FUNSPN(NUN))
+*----
+*  PRINT IMPORTED FLUX
+*----
+      IF(IMPX.GT.5) THEN
+         WRITE(IUNOUT,'(//33H I M P O R T E D    F L U X E S (,I5,
+     1   3H ):)') IG
+         WRITE(IUNOUT,'(1P,4(5X,E15.7))') (FUNKNO(I),I=1,NUN)
+      ENDIF
+*
+*----
+*  REBUILD KEYFLX FOR 2D HEXAGONAL CASE
+*----
+      ! NLOZH - num. of loz. per hexagon
+      ! SBMSH - num. of submeshes per lozenge (integer)
+      ! SPLTL - split of the lozenge (ISPLH)
+      IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9))THEN
+         ALLOCATE(TMPKEY(NREG),ORIKEY(NREG))
+         TMPKEY(:) = 0
+         ORIKEY(1:NREG) = KEYFLX(1:NREG)
+         IND = 0
+         JND = 0
+         NLOZH  = 3*ISPLH**2
+         SBMSH  = NLOZH/3
+         SPLTL  = ISPLH
+         DO IZ=1,LZ
+         DO IH=1,NHEX
+            DO IM=1,SBMSH
+               REM=MOD(IM-1,SPLTL)
+               IF((REM.EQ.0).AND.(SBMSH.NE.1))THEN
+                  JND = (IH-1)*NLOZH + SBMSH - (IM/SPLTL)
+                  JND = JND + (IZ-1)*LX
+               ELSEIF((REM.NE.0).AND.(SBMSH.NE.1))THEN
+                  JND = JND - (SBMSH*3) - SPLTL
+               ENDIF
+               DO ILZ=1,3
+                  IND = (IZ-1)*LX +(IH-1)*NLOZH +(IM-1)*3 +(ILZ-1) +1
+                  IF(SBMSH.EQ.1) JND = IND
+                  TMPKEY(IND) = KEYFLX(JND)
+                  JND = JND + SBMSH
+               ENDDO
+            ENDDO
+         ENDDO
+         ENDDO
+         KEYFLX(:) = TMPKEY(:)
+         DEALLOCATE(TMPKEY)
+      ENDIF
+*
+*----
+*  REARRANGE SPN FLUX IN SN ORDER, FOR P0 ISOTROPIC FLUX ONLY
+*----
+      FUNSPN(1:NUN) = FUNKNO(1:NUN)
+      FUNKNO(1:NUN) = 0.0
+
+      DO 100 IR=1,NREG
+         IF(MAT(IR).LE.0) GO TO 100
+         INDSN=KEYFLX(IR)
+         INDPN=KEYSPN(IR)
+         FUNKNO(INDSN)=FUNSPN(INDPN)
+  100 CONTINUE
+*
+*----
+*  RECUPERATE ORIGINAL KEYFLX FOR HEXAGONAL CASES
+*----
+      IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9))THEN
+         KEYFLX(1:NREG) = ORIKEY(1:NREG)
+         DEALLOCATE(ORIKEY)
+      ENDIF
+*----
+*  PRINT REARRANGED FLUX
+*----
+      IF(IMPX.GT.3) THEN
+         WRITE(IUNOUT,'(//37H R E A R R A N G E D    F L U X E S (,I5,
+     1   3H ):)') IG
+         WRITE(IUNOUT,'(1P,4(5X,E15.7))') (FUNKNO(I),I=1,NUN)
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(FUNSPN)
+      RETURN
+      END