Initial commit from Polytechnique Montreal

1 files changed, 316 insertions, 0 deletions

diff --git a/Dragon/src/MCGMRE.f b/Dragon/src/MCGMRE.f
new file mode 100644
index 0000000..5b5e6b0
--- /dev/null
+++ b/Dragon/src/MCGMRE.f
@@ -0,0 +1,316 @@
+*DECK MCGMRE
+      SUBROUTINE MCGMRE(SUBFFI,SUBFFA,SUBLDC,SUBSCH,CYCLIC,KPSYS,
+     1           IPRINT,IPTRK,IFTRAK,IPMACR,NDIM,KV,NUN,NLONG,PHIOUT,
+     2           NZON,MATALB,M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,SOUR,IAAC,
+     3           ISCR,LC,LFORW,PACA,ITST,MAXI,QFR,PHIIN,CAZ0,CAZ1,CAZ2,
+     4           CPO,ZMU,WZMU,VOL,EPS,ERRTOL,REPSI,NSTART,SIGAL,LPS,
+     5           NGROUP,NGEFF,NGIND,NCONV,LNCONV,NLIN,NFUNL,KEYFLX,
+     6           KEYCUR,STIS,NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,
+     7           FACSYM,IDIR,NBATCH)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solve the linear system obtained by the characteristics formalism
+* with GMRES iterative approach.
+*
+*Copyright:
+* Copyright (C) 2015 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 and A. Hebert
+*
+*Parameters:
+* 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.
+* KV      total number of volumes for which specific values
+*         of the neutron flux and reactions rates are required.
+* NUN     total number of unknowns in vectors QFR and 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 set.
+* NMAX    maximum number of elements in a track.
+* NANGL   number of tracking angles in the plane.
+* NREG    number of regions (volumes).
+* NSOUT   number of outer surfaces.
+* SOUR    undefined.
+* 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.
+* ITST    output: number of inner iterations.
+* MAXI    maximum number of inner iterations allowed.
+* 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.
+* WZMU    polar quadrature set.
+* VOL     volumes.
+* EPS     precision reached after min(MAXI,ITST) iterations.
+* ERRTOL  tolerance for stopping criterion. process is stopped
+*         as soon as ||Phi(n+1)-Phi(n)||/||Phi(n)|| <= EPS
+*         with ||.|| the euclidean norm.
+* REPSI   array containing precision of each iteration.
+* NSTART  undefined.
+* SIGAL   total cross-section and albedo array.
+* LPS     used in scr acceleration.
+* NGROUP  number of groups.
+* NGEFF   number of groups to process.
+* NGIND   index of the groups to process.
+* NCONV   array of convergence flag for each group.
+* LNCONV  number of unconverged groups.
+* 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 PHIOUT vector.
+* KEYCUR  position of current elements in PHIOUT vector.
+* 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
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+*----
+*  Subroutine arguments
+*----
+      TYPE(C_PTR) KPSYS(NGEFF),IPTRK,IPMACR
+      INTEGER NGEFF,IPRINT,IFTRAK,NDIM,KV,NUN,NLONG,NGROUP,NZON(NLONG),
+     1 M,NANI,NMU,NMAX,IAAC,LC,PACA,ITST(NGEFF),NSTART,MAXI,NANGL,NREG,
+     2 NSOUT,ISCR,LPS,NGIND(NGEFF),LNCONV,NLIN,NFUNL,KEYFLX(NREG,NLIN,
+     3 NFUNL),KEYCUR(NLONG-NREG),STIS,NPJJM,MATALB(-NSOUT:NREG),N2REG,
+     4 N2SOU,NZP,IDIR,NBATCH
+      REAL QFR(NUN,NGEFF),PHIIN(NUN,NGEFF),CPO(NMU),ZMU(NMU),
+     1 WZMU(NMU),VOL(NLONG),EPS(NGEFF),ERRTOL,REPSI(MAXI,NGEFF),
+     2 SIGAL(-6:M,NGEFF),DELU,FACSYM
+      DOUBLE PRECISION PHIOUT(NUN,NGEFF),CAZ0(NANGL),CAZ1(NANGL),
+     1 CAZ2(NANGL),SOUR(NUN,NGEFF)
+      LOGICAL LFORW,CYCLIC,NCONV(NGEFF),REBFLG,LPRISM
+      EXTERNAL SUBFFI,SUBFFA,SUBLDC,SUBSCH
+*---
+*   Local variables
+*---
+      REAL EPSINTO,FAC
+      LOGICAL RHSFLG
+      PARAMETER (FAC=100.0,MAXINT=200)
+*---
+*   Allocatable arrays
+*---
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: KMAX
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: RHS, GAR
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DENOM, RHO
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: R, G, C, S, FLOUT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: V, H
+*---
+*   Scratch storage allocation
+*---
+      ALLOCATE(DENOM(NGEFF),RHO(NGEFF),KMAX(NGEFF),FLOUT(NUN,NGEFF))
+      ALLOCATE(RHS(NUN,NGEFF),GAR(NUN,NGEFF),V(NUN,NGEFF,NSTART+1),
+     1 G(NSTART+1,NGEFF),H(NSTART+1,NSTART,NGEFF),C(NSTART+1,NGEFF),
+     2 S(NSTART+1,NGEFF),R(NUN,NGEFF))
+*
+      IF(MAXI.LT.3) CALL XABORT('MCGMRE: MAXI MUST BE >= 3.')
+      RHSFLG=.TRUE.
+      EPSINTO=ERRTOL/FAC
+      MAXIT=MAXI-1
+      NCONV(:)=.FALSE.
+      RHO(:)=0.0D0
+      LNCONV=0
+      DO II=1,NGEFF
+        DENOM(II)=SQRT(DOT_PRODUCT(QFR(:,II),QFR(:,II)))
+        RHO(II)=1.0D20
+        NCONV(II)=(DENOM(II) /= 0.0D0)
+        IF(NCONV(II)) LNCONV=LNCONV+1
+        ITST(II)=0
+        EPS(II)=0.0
+      ENDDO
+*---
+*   Global GMRES(M) iteration
+*---
+      ITER=0
+      DO WHILE((LNCONV /= 0) .AND. (ITER < MAXIT))
+        ITER=ITER+1
+        CALL MCGFL1(SUBFFI,SUBFFA,SUBLDC,SUBSCH,CYCLIC,KPSYS,IPRINT,
+     1       IPTRK,IFTRAK,IPMACR,NDIM,KV,NUN,NLONG,PHIOUT(1,1),NZON,
+     2       MATALB,M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NGROUP,NGEFF,NGIND,
+     3       SOUR,IAAC,ISCR,LC,LFORW,PACA,EPSINTO,MAXINT,NLIN,NFUNL,
+     4       KEYFLX,KEYCUR,QFR,PHIIN(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,
+     5       VOL,SIGAL,LPS,NCONV,.FALSE.,STIS,NPJJM,REBFLG,LPRISM,N2REG,
+     6       N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH)
+        ERROR=0.0D0
+        DO II=1,NGEFF
+          REPSI(ITER,II)=0.0
+          IF(.NOT.NCONV(II)) CYCLE
+          R(:,II)=PHIOUT(:,II)-PHIIN(:,II)
+          RHO(II)=SQRT(DOT_PRODUCT(R(:,II),R(:,II)))
+          REPSI(ITER,II)=REAL(RHO(II)/DENOM(II))
+          IF(IPRINT.GT.4) WRITE(6,200) ITER,II,REPSI(ITER,II)
+          EPS(II)=REPSI(ITER,II)
+          ITST(II)=ITER
+          ERROR=MAX(ERROR,RHO(II)/DENOM(II))
+          IF(RHO(II) < ERRTOL*DENOM(II)) THEN
+            NCONV(II)=.FALSE.
+            LNCONV=LNCONV-1
+          ENDIF
+        ENDDO
+*
+*       Test do termination on entry
+        IF(LNCONV == 0) EXIT
+*
+        H(:,:,:)=0.0D0
+        V(:,:,:)=0.0D0
+        C(:,:)=0.0D0
+        S(:,:)=0.0D0
+        G(:,:)=0.0D0
+        KMAX(:)=0
+        DO II=1,NGEFF
+          IF(.NOT.NCONV(II)) CYCLE
+          G(1,II)=RHO(II)
+          V(:,II,1)=R(:,II)/RHO(II)
+        ENDDO
+*---
+*   Evaluate RHS of the linear system
+*---
+        IF(RHSFLG) THEN
+          IF(IPRINT > 3) WRITE(6,100) IAAC,ISCR
+          RHS(:,:)=0.0
+          CALL MCGFL1(SUBFFI,SUBFFA,SUBLDC,SUBSCH,CYCLIC,KPSYS,IPRINT,
+     1         IPTRK,IFTRAK,IPMACR,NDIM,KV,NUN,NLONG,FLOUT(1,1),NZON,
+     2         MATALB,M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NGROUP,NGEFF,
+     3         NGIND,SOUR,IAAC,ISCR,LC,LFORW,PACA,EPSINTO,MAXINT,NLIN,
+     4         NFUNL,KEYFLX,KEYCUR,QFR,RHS(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,
+     5         WZMU,VOL,SIGAL,LPS,NCONV,.FALSE.,STIS,NPJJM,REBFLG,
+     6         LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH)
+          DO II=1,NGEFF
+            IF(NCONV(II)) RHS(:,II)=REAL(FLOUT(:,II))
+          ENDDO
+          RHSFLG=.FALSE.
+        ENDIF
+*---
+*   GMRES(1) iteration
+*---
+        K=0
+        DO WHILE((LNCONV /= 0) .AND. (K < NSTART) .AND. (ITER < MAXIT))
+          K=K+1
+          ITER=ITER+1
+          GAR(:,:)=REAL(V(:,:,K))
+          CALL MCGFL1(SUBFFI,SUBFFA,SUBLDC,SUBSCH,CYCLIC,KPSYS,IPRINT,
+     1         IPTRK,IFTRAK,IPMACR,NDIM,KV,NUN,NLONG,FLOUT(1,1),NZON,
+     2         MATALB,M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NGROUP,NGEFF,
+     3         NGIND,SOUR,IAAC,ISCR,LC,LFORW,PACA,EPSINTO,MAXINT,NLIN,
+     4         NFUNL,KEYFLX,KEYCUR,QFR,GAR(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,
+     5         WZMU,VOL,SIGAL,LPS,NCONV,.FALSE.,STIS,NPJJM,REBFLG,
+     6         LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH)
+          V(:,:,K+1)=V(:,:,K)-FLOUT(:,:)+RHS(:,:)
+          ERROR=0.0D0
+          DO II=1,NGEFF
+            REPSI(ITER,II)=0.0
+            IF(.NOT.NCONV(II)) CYCLE
+            KMAX(II)=K
+*
+*           Modified Gram-Schmidt
+            DO J=1,K
+              H(J,K,II)=DOT_PRODUCT(V(:,II,J),V(:,II,K+1))
+              V(:,II,K+1)=V(:,II,K+1)-H(J,K,II)*V(:,II,J)
+            ENDDO
+            H(K+1,K,II)=SQRT(DOT_PRODUCT(V(:,II,K+1),V(:,II,K+1)))
+*
+*           Reorthogonalize
+            DO J=1,K
+              HR=DOT_PRODUCT(V(:,II,J),V(:,II,K+1))
+              H(J,K,II)=H(J,K,II)+HR
+              V(:,II,K+1)=V(:,II,K+1)-HR*V(:,II,J)
+            ENDDO
+            H(K+1,K,II)=SQRT(DOT_PRODUCT(V(:,II,K+1),V(:,II,K+1)))
+*
+            ! Watch out do happy breakdown 
+            IF(H(K+1,K,II) /= 0.0D0) V(:,II,K+1)=V(:,II,K+1)/H(K+1,K,II)
+*
+*           Form and store the information for the new Givens rotation
+            DO I=1,K-1
+              W1=C(I,II)*H(I,K,II)-S(I,II)*H(I+1,K,II)
+              W2=S(I,II)*H(I,K,II)+C(I,II)*H(I+1,K,II)
+              H(I,K,II)=W1
+              H(I+1,K,II)=W2
+            ENDDO
+            ZNU=SQRT(H(K,K,II)**2+H(K+1,K,II)**2)
+            IF(ZNU /= 0.0D0) THEN
+              C(K,II)=H(K,K,II)/ZNU
+              S(K,II)=-H(K+1,K,II)/ZNU
+              H(K,K,II)=C(K,II)*H(K,K,II)-S(K,II)*H(K+1,K,II)
+              H(K+1,K,II)=0.0D0
+              W1=C(K,II)*G(K,II)-S(K,II)*G(K+1,II)
+              W2=S(K,II)*G(K,II)+C(K,II)*G(K+1,II)
+              G(K,II)=W1
+              G(K+1,II)=W2
+            ENDIF
+*
+*           Update the residual norm
+            RHO(II)=ABS(G(K+1,II))
+            REPSI(ITER,II)=REAL(RHO(II)/DENOM(II))
+            IF(IPRINT.GT.4) WRITE(6,200) ITER,II,REPSI(ITER,II)
+            EPS(II)=REPSI(ITER,II)
+            ITST(II)=ITER
+            IF(RHO(II) < ERRTOL*DENOM(II)) THEN
+              NCONV(II)=.FALSE.
+              LNCONV=LNCONV-1
+            ENDIF
+            ERROR=MAX(ERROR,RHO(II)/DENOM(II))
+          ENDDO
+        ENDDO
+*---
+*   At this point either K > NSTART or RHOGRP < ERRTOL. It's time to
+*   compute PHIIN and cycle.
+*---
+        DO II=1,NGEFF
+          K=KMAX(II)
+          IF(K == 0) CYCLE
+          G(K,II)=G(K,II)/H(K,K,II)
+          DO L=K-1,1,-1
+            W1=G(L,II)-DOT_PRODUCT(H(L,L+1:K,II),G(L+1:K,II))
+            G(L,II)=W1/H(L,L,II)
+          ENDDO
+          DO J=1,K
+            PHIIN(:,II)=PHIIN(:,II)+REAL(G(J,II)*V(:,II,J))
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  Scratch storage deallocation
+*----
+      DEALLOCATE(R, S, C, H, G, V, GAR, RHS)
+      DEALLOCATE(FLOUT, KMAX, RHO, DENOM)
+      RETURN
+*
+  100 FORMAT(' MCGMRE: RHS CALCULATED WITH AAC-SCR : ',I2,1H-,I1)
+  200 FORMAT(' MCGMRE: EPS Iteration ',I5,' Group ',I5,2X,1P,E16.7)
+      END