Initial commit from Polytechnique Montreal

1 files changed, 342 insertions, 0 deletions

diff --git a/Dragon/src/MCGFL1.f b/Dragon/src/MCGFL1.f
new file mode 100644
index 0000000..efaacc1
--- /dev/null
+++ b/Dragon/src/MCGFL1.f
@@ -0,0 +1,342 @@
+*DECK MCGFL1
+      SUBROUTINE MCGFL1(SUBFFI,SUBFFA,SUBLDC,SUBSCH,CYCLIC,KPSYS,IPRINT,
+     1           IPTRK,IFTRAK,IPMACR,NDIM,K,KPN,NLONG,PHIOUT,NZON,
+     2           MATALB,M,NANI,NMU,N2MAX,NANGL,NREG,NSOUT,NG,NGEFF,
+     3           NGIND,S,IAAC,ISCR,LC,LFORW,PACA,EPSACC,MAXACC,NLIN,
+     4           NFUNL,KEYFLX,KEYCUR,QFR,PHIIN,CAZ0,CAZ1,CAZ2,CPO,ZMU,
+     5           WZMU,V,SIGAL,LPS,NCONV,LAST,STIS,NPJJM,REBFLG,LPRISM,
+     6           N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one inner iteration of the characteristics method.
+*
+*Copyright:
+* Copyright (C) 2002 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): R. Le Tellier
+*
+*Parameters: input/output
+* SUBFFI  flux integration subroutine with isotropic source.
+* SUBFFA  flux integration subroutine with anisotropic source.
+* SUBLDC  flux integration subroutine with linear-discontinuous source.
+* SUBSCH  track coefficients calculation subroutine.
+* CYCLIC  cyclic tracking flag.
+* KPSYS   pointer array for each group properties.
+* IPRINT  print parameter (equal to zero for no print).
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IFTRAK  tracking file unit number.
+* IPMACR  pointer to the macrolib LCM object.
+* NDIM    number of dimensions for the geometry.
+* K       total number of volumes for which specific values
+*         of the neutron flux and reactions rates are required.
+* KPN     total number of unknowns per group in vector PHIIN.
+* NLONG   number of spatial unknowns.
+* PHIOUT  output flux vector.
+* NZON    mixture-albedo index array in MCCG format.
+* MATALB  albedo-mixture index array in MOCC format.
+* M       number of material mixtures.
+* NANI    scattering anisotropy (=1 for isotropic scattering).
+* NMU     order of the polar quadrature in 2D / 1 in 3D.
+* N2MAX   maximum number of elements in a track.
+* NANGL   number of tracking angles in the plane.
+* NREG    number of volumes.
+* NSOUT   number of outer surfaces.
+* NG      number of groups.
+* NGEFF   number of groups to process.
+* NGIND   index of the groups to process.
+* S       scratch.
+* IAAC    no acceleration / CDD acceleration of inner iterations (0/1).
+* ISCR    no acceleration / SCR acceleration of inner iterations (0/1).
+* LC      dimension of profiled matrices MCU and CQ.
+* LFORW   flag set to .false. to transpose the coefficient matrix.
+* PACA    type of preconditioner to solve the ACA corrective system.
+* EPSACC  stopping criterion for BICGSTAB in ACA resolution.
+* MAXACC  maximum number of iterations allowed for BICGSTAB in ACA
+*         resolution.
+* NLIN    number of polynomial components in flux spatial expansion.
+* NFUNL   number of moments of the flux (in 2D: NFUNL=NANI*(NANI+1)/2).
+* KEYFLX  position of flux elements in PHIIN vector.
+* KEYCUR  position of current elements in PHIIN vector.
+* QFR     input source vector.
+* PHIIN   input flux vector.
+* CAZ0    cosines of the tracking polar angles in 3D.
+* CAZ1    first cosines of the different tracking azimuthal angles.
+* CAZ2    second cosines of the different tracking azimuthal angles.
+* CPO     cosines of the different tracking polar angles in 2D.
+* ZMU     polar quadrature set in 2D.
+* WZMU    polar quadrature set in 2D.
+* V       volumes.
+* SIGAL   total cross-section and albedo array.
+* LPS     used in scr acceleration.
+* NCONV   logical array of convergence status for each group (.TRUE.
+*         not converged).
+* LAST    flag for SCR and ACA rebalancing.
+* STIS    source term isolation' option for flux integration.
+* NPJJM   number of pjj modes to store for STIS option.
+* REBFLG  ACA or SCR rebalancing flag.
+* LPRISM  3D prismatic extended tracking flag.
+* N2REG   number of regions in the 2D tracking if LPRISM.
+* N2SOU   number of external surfaces in the 2D tracking if LPRISM.
+* NZP     number of z-plans if LPRISM.
+* DELU    input track spacing for 3D track reconstruction if LPRISM.
+* FACSYM  tracking symmetry factor for maximum track length if LPRISM.
+* IDIR    direction of fundamental current for TIBERE with MoC 
+*         (=0,1,2,3). 
+* NBATCH  number of tracks processed in each OpenMP core (default: =1).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*---
+* SUBROUTINES ARGUMENTS
+*---
+      TYPE(C_PTR) KPSYS(NGEFF),IPTRK,IPMACR
+      INTEGER NGEFF,IPRINT,IFTRAK,NDIM,K,KPN,NLONG,NG,NGIND(NGEFF),
+     1 NZON(NLONG),M,NANI,NMU,N2MAX,NANGL,IAAC,LC,PACA,NREG,NSOUT,
+     2 ISCR,LPS,NLIN,NFUNL,KEYFLX(NREG,NLIN,NFUNL),KEYCUR(NLONG-NREG),
+     3 MAXACC,STIS,NPJJM,MATALB(-NSOUT:NREG),N2REG,N2SOU,NZP,IDIR,NBATCH
+      REAL QFR(KPN,NGEFF),PHIIN(KPN,NGEFF),CPO(NMU),ZMU(NMU),WZMU(NMU),
+     1 V(NLONG),SIGAL(-6:M,NGEFF),EPSACC,DELU,FACSYM
+      DOUBLE PRECISION CAZ0(NANGL),CAZ1(NANGL),CAZ2(NANGL),
+     1 PHIOUT(KPN,NGEFF),S(KPN,NGEFF)
+      LOGICAL LFORW,CYCLIC,NCONV(NGEFF),LAST,REBFLG,LPRISM
+      EXTERNAL SUBFFI,SUBFFA,SUBLDC,SUBSCH
+*---
+* LOCAL VARIABLES
+*---
+      TYPE(C_PTR) JPSYS
+      REAL T1,T2,T3
+      INTEGER NCODE(6),SSYM
+      INTEGER, DIMENSION(1), TARGET :: IDUMMY
+      CHARACTER TEXT4*4
+      LOGICAL MACFLG,COMBFLG
+      INTEGER ICREB,ITYLCM
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      TYPE(C_PTR) IBC_PTR,XSSC_PTR,PJJM_PTR
+      TYPE(C_PTR) INDREG_PTR,Z_PTR,VNORF_PTR,CMU_PTR,CMUI_PTR,SMU_PTR,
+     1 SMUI_PTR,TMU_PTR,TMUI_PTR
+      INTEGER, ALLOCATABLE, DIMENSION(:) ::KEYANI
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ISGNR
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: XSIXYZ
+      INTEGER, POINTER, DIMENSION(:) :: IBC,INDREG,PJJM
+      REAL, POINTER, DIMENSION(:) :: XSSC,Z
+      DOUBLE PRECISION, POINTER, DIMENSION(:) :: VNORF,CMU,CMUI,SMU,
+     1 SMUI,TMU,TMUI
+*
+      PHIOUT(:KPN,:NGEFF)=0.0D0
+*
+      REWIND IFTRAK
+      READ(IFTRAK) TEXT4,NCOMNT,NBTR,IFMT
+      DO ICOM=1,NCOMNT
+         READ(IFTRAK)
+      ENDDO
+      READ(IFTRAK) (NITMA, II=1,7),MXSUB,NITMA
+      DO ICOM=1,6
+         READ(IFTRAK)
+      ENDDO
+*
+      IF(LPRISM) THEN
+         CALL LCMGET(IPTRK,'NCODE',NCODE)
+         IF(NCODE(6).EQ.30) THEN
+            IF(NCODE(5).EQ.30) THEN
+*           Z- and Z+ surfaces symmetry
+               SSYM=2
+            ELSE
+*           Z+ symmetry
+               SSYM=1
+            ENDIF
+         ELSE
+            SSYM=0
+         ENDIF
+         NDIMB=3
+         NMAX=(INT(FACSYM)+1)*N2MAX*(NZP+2)
+      ELSE
+         NDIMB=NDIM
+         NMAX=N2MAX
+      ENDIF
+*----
+*  SOURCE ELEMENTS CALCULATION
+*----
+      IF(NLONG-NREG.GT.0) THEN
+        CALL LCMGPD(IPTRK,'BC-REFL+TRAN',IBC_PTR)
+        CALL C_F_POINTER(IBC_PTR,IBC,(/ NLONG-NREG /))
+      ELSE
+        IBC => IDUMMY
+      ENDIF
+      DO II=1,NGEFF
+         IF(NCONV(II)) THEN
+           JPSYS=KPSYS(II)
+           CALL LCMGPD(JPSYS,'DRAGON-S0XSC',XSSC_PTR)
+           CALL C_F_POINTER(XSSC_PTR,XSSC,(/ (M+1)*NANI /))
+           CALL MCGFCS(NLONG,NDIMB,NZON,QFR(1,II),PHIIN(1,II),M,NANI,
+     1        NLIN,NFUNL,XSSC,S(1,II),KPN,NREG,IPRINT,KEYFLX,KEYCUR,
+     2        IBC,SIGAL(-6,II),STIS)
+         ENDIF
+      ENDDO
+*---
+* GENERATE ALL SIGNS FOR SPHERICAL HARMONICS
+*---
+      NANIX=NANI
+      IF(NLIN.EQ.3) NANIX=MAX(NANI,3)
+      IF(NDIM.EQ.1) THEN
+         NFUNLX=NANIX
+         NMOD=2
+      ELSE IF((.NOT.LPRISM).AND.(NDIM.EQ.2)) THEN
+         NFUNLX=NANIX*(NANIX+1)/2
+         NMOD=4
+      ELSE ! NDIM.EQ.3
+         NFUNLX=NANIX*NANIX
+         NMOD=8
+      ENDIF
+      ALLOCATE(ISGNR(NMOD,NFUNLX),KEYANI(NFUNLX))
+      CALL MOCIK3(NANIX-1,NFUNLX,NMOD,ISGNR,KEYANI)
+*----
+*  FLUX INTEGRATION UPON THE TRACKING FILE
+*----
+      CALL KDRCPU(T1)
+      IF(CYCLIC) THEN
+*     --------------------------------
+*     Method of Cyclic Characteristics
+*     --------------------------------
+         CALL MOCFCF(SUBFFI,SUBFFA,SUBLDC,SUBSCH,IFTRAK,NBTR,MXSUB,
+     1        N2MAX,NDIM,KPN,NREG,NSOUT,M,6,NGEFF,NANGL,NMU,NANI,
+     2        NFUNL,NMOD,NANIX,NLIN,NFUNLX,KEYFLX,MATALB,NCONV,SIGAL,
+     3        CAZ1,CAZ2,CPO,ZMU,WZMU,S,ISGNR,IDIR,NBATCH,PHIOUT)
+      ELSE
+*     ------------------------------------
+*     Method of Non-Cyclic Characteristics
+*     ------------------------------------
+         IF(LPRISM) THEN
+*        3D PRISMATIC GEOMETRY CONSTRUCTED FROM A 2D TRACKING
+            CALL LCMSIX(IPTRK,'PROJECTION',1)
+            CALL LCMGPD(IPTRK,'IND2T3',INDREG_PTR)
+            CALL LCMGPD(IPTRK,'ZCOORD',Z_PTR)
+            CALL LCMGPD(IPTRK,'VNORF',VNORF_PTR)
+            CALL LCMGPD(IPTRK,'CMU',CMU_PTR)
+            CALL LCMGPD(IPTRK,'CMUI',CMUI_PTR)
+            CALL LCMGPD(IPTRK,'SMU',SMU_PTR)
+            CALL LCMGPD(IPTRK,'SMUI',SMUI_PTR)
+            CALL LCMGPD(IPTRK,'TMU',TMU_PTR)
+            CALL LCMGPD(IPTRK,'TMUI',TMUI_PTR)
+            CALL LCMSIX(IPTRK,'PROJECTION',2)
+*
+            CALL C_F_POINTER(INDREG_PTR,INDREG,
+     1                       (/ (N2REG+N2SOU+1)*(NZP+2) /))
+            CALL C_F_POINTER(Z_PTR,Z,(/ NZP+1 /))
+            CALL C_F_POINTER(VNORF_PTR,VNORF,(/ NREG*NANGL*NMU*2 /))
+            CALL C_F_POINTER(CMU_PTR,CMU,(/ NMU /))
+            CALL C_F_POINTER(CMUI_PTR,CMUI,(/ NMU /))
+            CALL C_F_POINTER(SMU_PTR,SMU,(/ NMU /))
+            CALL C_F_POINTER(SMUI_PTR,SMUI,(/ NMU /))
+            CALL C_F_POINTER(TMU_PTR,TMU,(/ NMU /))
+            CALL C_F_POINTER(TMUI_PTR,TMUI,(/ NMU /))
+            CALL MCGPTF(SUBFFI,SUBFFA,SUBSCH,IFTRAK,NBTR,N2MAX,KPN,
+     1           K,NREG,M,NGEFF,NANGL,NMU,NANI,NFUNL,NMOD,KEYFLX,
+     2           KEYCUR,NZON,NCONV,CAZ1,CAZ2,CPO,WZMU,PHIOUT,S,SIGAL,
+     3           ISGNR,NMAX,NZP,N2REG,N2SOU,DELU,INDREG,Z,VNORF,CMU,
+     4           CMUI,SMU,SMUI,TMU,TMUI,SSYM,IDIR)
+         ELSE
+*        REGULAR 2D OR 3D GEOMETRY
+            ALLOCATE(XSIXYZ(NSOUT,3))
+            XSIXYZ(:NSOUT,:3)=0.0D0
+            CALL LCMLEN(IPTRK,'XSI$MCCG',ICREB,ITYLCM)
+            IF(ICREB.EQ.3*NSOUT) CALL LCMGET(IPTRK,'XSI$MCCG',XSIXYZ)
+            CALL MCGFCF(SUBFFI,SUBFFA,SUBLDC,SUBSCH,IFTRAK,NBTR,N2MAX,
+     1           NDIM,KPN,K,NREG,M,NGEFF,NANGL,NMU,NANI,NFUNL,NMOD,
+     2           NANIX,NLIN,NFUNLX,KEYFLX,KEYCUR,NZON,NCONV,CAZ0,CAZ1,
+     3           CAZ2,CPO,ZMU,WZMU,S,SIGAL,ISGNR,IDIR,NSOUT,NBATCH,
+     4           XSIXYZ(1,IDIR),PHIOUT)
+            DEALLOCATE(XSIXYZ)
+         ENDIF
+      ENDIF
+*
+      IF(STIS.EQ.1) THEN
+         CALL LCMGPD(IPTRK,'PJJIND$MCCG',PJJM_PTR)
+         CALL C_F_POINTER(PJJM_PTR,PJJM,(/ NPJJM*2 /))
+         CALL MCGFST(NGEFF,KPSYS,NCONV,KPN,NLONG,NREG,NANI,NFUNL,NPJJM,
+     1        KEYFLX,KEYCUR,PJJM,NZON,V,S,PHIOUT,IDIR)
+      ELSEIF(STIS.EQ.-1) THEN
+         DO II=1,NGEFF
+         IF(NCONV(II)) THEN
+            CALL MCGFMC(KPN,NLONG,NREG,M,NANI,NFUNL,NZON,KEYFLX,KEYCUR,
+     1           PHIOUT(1,II),V,S(1,II),SIGAL(0,II),KEYANI)
+         ENDIF
+         ENDDO
+      ELSE
+         DO II=1,NGEFF
+         IF(NCONV(II)) THEN
+            DO I=1,NLONG
+            IF(V(I).GT.0.) THEN 
+               IF(NZON(I).LT.0) THEN
+                  IND=KEYCUR(I-NREG)
+                  PHIOUT(IND,II)=PHIOUT(IND,II)/V(I)
+               ELSE
+                  DO IL=1,NFUNL
+                     DO IU=1,NLIN
+                        IND=KEYFLX(I,IU,IL)
+                        PHIOUT(IND,II)=PHIOUT(IND,II)/V(I)
+                     ENDDO
+                  ENDDO
+               ENDIF
+            ENDIF
+            ENDDO
+         ENDIF
+         ENDDO
+      ENDIF
+      CALL KDRCPU(T2)
+*----
+*  PRECONDITIONING TECHNIQUES
+*----
+      MACFLG=((LAST).AND.(C_ASSOCIATED(IPMACR)).AND.((IAAC.GT.1).OR.
+     1       (ISCR.GT.1)))
+      REBFLG=(REBFLG.AND.MACFLG)
+      COMBFLG=((IAAC.GT.0).AND.(ISCR.GT.0))
+      IF(IAAC.GT.0) THEN
+*     ---------------------------------
+*     Algebraic Collapsing Acceleration
+*     ---------------------------------
+         IF(REBFLG) MAXACC=IAAC
+         CALL MCGFCA(IPTRK,KPSYS,IPMACR,IPRINT,NLONG,NG,NGEFF,KPN,NREG,
+     1        NANI,NFUNL,M,LC,LFORW,PACA,KEYFLX,KEYCUR,NZON,NGIND,NCONV,
+     2        MAXACC,EPSACC,MACFLG,REBFLG,PHIOUT,PHIIN,COMBFLG,NPJJM,
+     3        KEYANI,IDIR)
+         CALL KDRCPU(T3)
+*     ---------------------------------
+      ENDIF
+      NLON2=NLONG
+      IF(COMBFLG) NLON2=NREG
+      IF(ISCR.GT.0) THEN
+*     ---------------------------------
+*     Self-Collision Rebalancing Method
+*     ---------------------------------
+         IF(REBFLG) THEN
+*        rebalancing            
+            MAXSCR=ISCR
+         ELSE
+            MAXSCR=1
+         ENDIF
+         CALL MCGSCR(IPTRK,KPSYS,IPMACR,IPRINT,NLON2,NG,NGEFF,KPN,K,
+     1        NREG,NANI,NFUNL,M,LPS,KEYFLX,KEYCUR,NZON,NGIND,NCONV,
+     2        MAXSCR,EPSACC,MACFLG,PHIOUT,PHIIN,V,NPJJM,KEYANI,IDIR)
+         CALL KDRCPU(T3)
+*     ---------------------------------
+      ENDIF
+      DEALLOCATE(KEYANI,ISGNR)
+      IF(IPRINT.GT.1) THEN
+         WRITE(6,100) ' FLUX INTEGRATION       ',(T2-T1)
+         IF((IAAC.GT.0).OR.(ISCR.NE.0)) THEN
+            WRITE(6,100) ' ACCELERATION           ',(T3-T2)
+         ENDIF
+      ENDIF
+      RETURN
+*
+ 100  FORMAT(10X,'MCGFL1: -->>TIME SPENT IN ',A24,':',F13.3)
+      END