Initial commit from Polytechnique Montreal

1 files changed, 229 insertions, 0 deletions

diff --git a/Dragon/src/SNGMRE.f b/Dragon/src/SNGMRE.f
new file mode 100644
index 0000000..10a439f
--- /dev/null
+++ b/Dragon/src/SNGMRE.f
@@ -0,0 +1,229 @@
+*DECK SNGMRE
+      SUBROUTINE SNGMRE (KPSYS,NGIND,IPTRK,IMPX,NGEFF,NREG,NBMIX,NUN,
+     1 NSTART,MAXIT,EPSINR,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,NBS,KPSOU1,
+     2 KPSOU2,FLUXC)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solve N-group transport equation for fluxes using the discrete
+* ordinates (SN) method with a GMRES(m) acceleration.
+*
+*Copyright:
+* Copyright (C) 2007 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. Hebert
+*
+*Parameters: input
+* KPSYS   pointer to the assembly matrices. KPSYS is an array of
+*         directories.
+* NGIND   energy group indices assign to the NGEFF set.
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IMPX    print flag (equal to zero for no print).
+* NGEFF   number of energy groups processed in parallel.
+* NREG    total number of regions for which specific values of the
+*         neutron flux and reactions rates are required.
+* NBMIX   number of mixtures.
+* NUN     total number of unknowns in vectors SUNKNO and FUNKNO.
+* NSTART  restarts the GMRES method every NSTART iterations.
+* MAXIT   maximum number of inner iterations.
+* EPSINR  convergence criterion.
+* MAT     index-number of the mixture type assigned to each volume.
+* VOL     volumes.
+* KEYFLX  position of averaged flux elements in FUNKNO vector.
+* SUNKNO  input source vector.
+*
+*Parameters: input/output
+* FUNKNO  unknown vector.
+* FLUXC   flux at the cutoff energy.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER     NGEFF,NGIND(NGEFF),IMPX,NREG,NBMIX,NUN,NSTART,MAXIT,
+     1            MAT(NREG),KEYFLX(NREG),NBS(NGEFF)
+      TYPE(C_PTR) KPSYS(NGEFF),IPTRK,KPSOU1(NGEFF),KPSOU2(NGEFF)
+      REAL        EPSINR,VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF)
+      REAL,OPTIONAL :: FLUXC(NREG)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER  (IUNOUT=6)
+      REAL        SDOT
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: KMAX
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: RR,QQ
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: EPS1,RHO
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: G,C,S,X
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: V,H
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: INCONV
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(G(NSTART+1,NGEFF),C(NSTART+1,NGEFF),S(NSTART+1,NGEFF),
+     1 X(NUN,NGEFF),V(NUN,NSTART+1,NGEFF),H(NSTART+1,NSTART,NGEFF),
+     1 EPS1(NGEFF),INCONV(NGEFF),KMAX(NGEFF),RHO(NGEFF))
+      ALLOCATE(RR(NUN,NGEFF),QQ(NUN,NGEFF))
+*----
+*  GLOBAL GMRES ITERATION.
+*----
+      DO II=1,NGEFF
+        EPS1(II)=EPSINR*SQRT(SDOT(NUN,SUNKNO(1,II),1,SUNKNO(1,II),1))
+        RHO(II)=1.0E10
+        INCONV(II)=(EPS1(II).NE.0.0)
+      ENDDO
+      LNCONV=NGEFF
+      ITER=0
+      NITER=1
+      DO WHILE((LNCONV.GT.0).AND.(ITER.LT.MAXIT))
+        RR(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF)
+        CALL SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG,
+     1  NBMIX,NUN,MAT,VOL,KEYFLX,RR,SUNKNO(1,1),NITER,NBS,KPSOU1,
+     2  KPSOU2,FLUXC)
+        NITER=NITER+1
+        DO II=1,NGEFF
+          IF(.NOT.INCONV(II)) CYCLE
+          X(:NUN,II)=RR(:NUN,II)-FUNKNO(:NUN,II)
+          RHO(II)=SQRT(DOT_PRODUCT(X(:NUN,II),X(:NUN,II)))
+          IF(RHO(II).LT.EPS1(II)) THEN
+            LNCONV=LNCONV-1
+            INCONV(II)=.FALSE.
+          ENDIF
+        ENDDO
+*----
+*  TEST FOR TERMINATION ON ENTRY
+*----
+        IF(LNCONV.EQ.0) GO TO 100
+*
+        G(:NSTART+1,:NGEFF)=0.0D0
+        C(:NSTART+1,:NGEFF)=0.0D0
+        S(:NSTART+1,:NGEFF)=0.0D0
+        V(:NUN,:NSTART+1,:NGEFF)=0.0D0
+        H(:NSTART+1,:NSTART,:NGEFF)=0.0D0
+        KMAX(:NGEFF)=0
+        DO II=1,NGEFF
+          IF(.NOT.INCONV(II)) CYCLE
+          G(1,II)=RHO(II)
+          DO I=1,NUN
+            V(I,1,II)=X(I,II)/RHO(II)
+            X(I,II)=FUNKNO(I,II)
+          ENDDO
+        ENDDO
+*----
+*  GMRES(1) ITERATION
+*----
+        K=0
+        DO WHILE((LNCONV.GT.0).AND.(K.LT.NSTART).AND.(ITER.LT.MAXIT))
+          K=K+1
+          ITER=ITER+1
+          IF(IMPX.GT.2) WRITE(IUNOUT,300) ITER,MAXVAL(RHO(:NGEFF)),
+     1    LNCONV
+          RR(:NUN,:NGEFF)=0.0
+          QQ(:NUN,:NGEFF)=0.0
+          DO II=1,NGEFF
+            IF(.NOT.INCONV(II)) CYCLE
+            RR(:NUN,II)=REAL(V(:NUN,K,II))
+          ENDDO
+          CALL SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG,
+     1    NBMIX,NUN,MAT,VOL,KEYFLX,RR,QQ(1,1),NITER,NBS,KPSOU1,
+     2    KPSOU2,FLUXC)
+          NITER=NITER+1
+          DO II=1,NGEFF
+            IF(.NOT.INCONV(II)) CYCLE
+            V(:NUN,K+1,II)=V(:NUN,K,II)-RR(:NUN,II)
+            KMAX(II)=K
+*----
+*  MODIFIED GRAM-SCHMIDT
+*----
+            DO J=1,K
+              HR=DOT_PRODUCT(V(:NUN,J,II),V(:NUN,K+1,II))
+              H(J,K,II)=HR
+              V(:NUN,K+1,II)=V(:NUN,K+1,II)-HR*V(:NUN,J,II)
+            ENDDO
+            H(K+1,K,II)=SQRT(DOT_PRODUCT(V(:NUN,K+1,II),V(:NUN,K+1,II)))
+*----
+*  REORTHOGONALIZE
+*----
+            DO J=1,K
+              HR=DOT_PRODUCT(V(:NUN,J,II),V(:NUN,K+1,II))
+              H(J,K,II)=H(J,K,II)+HR
+              V(:NUN,K+1,II)=V(:NUN,K+1,II)-HR*V(:NUN,J,II)
+            ENDDO
+            H(K+1,K,II)=SQRT(DOT_PRODUCT(V(:NUN,K+1,II),V(:NUN,K+1,II)))
+*----
+*  WATCH OUT FOR HAPPY BREAKDOWN 
+*----
+            IF(H(K+1,K,II).NE.0.0) THEN
+              V(:NUN,K+1,II)=V(:NUN,K+1,II)/H(K+1,K,II)
+            ENDIF
+*----
+*  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.NE.0.0) 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))
+            IF(RHO(II).LE.EPS1(II)) THEN
+              INCONV(II)=.FALSE.
+              LNCONV=LNCONV-1
+            ENDIF
+          ENDDO
+        ENDDO
+*----
+*  AT THIS POINT EITHER K > NSTART OR RHO < EPS1.
+*  IT'S TIME TO COMPUTE X AND CYCLE.
+*----
+        DO II=1,NGEFF
+          K=KMAX(II)
+          IF(K.EQ.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
+            X(:,II)=X(:,II)+G(J,II)*V(:,J,II)
+          ENDDO
+          FUNKNO(:,II)=REAL(X(:,II))
+        ENDDO
+      ENDDO
+*
+      IF(IMPX.GT.2) WRITE(IUNOUT,'(32H SNGMRE: NUMBER OF ONE-SPEED ITE,
+     1 8HRATIONS=,I5,1H.)') ITER
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+  100 DEALLOCATE(QQ,RR)
+      DEALLOCATE(RHO,KMAX,INCONV,EPS1,H,V,X,S,C,G)
+      RETURN
+*
+  300 FORMAT(28H SNGMRE: ONE-SPEED ITERATION,I4,10H  L2 NORM=,1P,E11.4,
+     1 23H  NON-CONVERGED GROUPS=,I5)
+      END