Initial commit from Polytechnique Montreal

1 files changed, 222 insertions, 0 deletions

diff --git a/Trivac/src/GPTMRA.f b/Trivac/src/GPTMRA.f
new file mode 100755
index 0000000..e2dda52
--- /dev/null
+++ b/Trivac/src/GPTMRA.f
@@ -0,0 +1,222 @@
+*DECK GPTMRA
+      SUBROUTINE GPTMRA(IPTRK,IPSYS,IPFLUP,LADJ,LL4,ITY,NUN,NGRP,ICL1,
+     1 ICL2,IMPX,NADI,MAXINR,NSTART,MAXX0,EPS2,EPSINR,EVAL,EVECT,ADECT,
+     2 EASS,SOUR,TKT,TKB,ZNORM,ITER)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of a multigroup fixed source eigenvalue problem for the
+* calculation of a gpt solution in Trivac. Use the preconditioned power
+* method with GMRES(m) acceleration.
+*
+*Copyright:
+* Copyright (C) 2019 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
+* IPTRK   L_TRACK pointer to the tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPFLUP  L_FLUX pointer to the gpt solution
+* LADJ    flag set to .TRUE. for adjoint solution acceleration.
+* LL4     order of the system matrices.
+* ITY     type of solution (2: classical Trivac; 3: Thomas-Raviart).
+* NUN     number of unknowns in each energy group.
+* NGRP    number of energy groups.
+* ICL1    number of free up-scattering iterations in one cycle of the
+*         inverse power method.
+* ICL2    number of accelerated up-scattering iterations in one cycle.
+* IMPX    print parameter (set to 0 for no printing).
+* NADI    initial number of inner ADI iterations per outer iteration.
+* MAXINR  maximum number of thermal iterations.
+* NSTART  restarts the GMRES method every NSTART iterations.
+* MAXX0   maximum number of outer iterations
+* EPS2    outer iteration convergence criterion
+* EPSINR  thermal iteration convergence criterion
+* EVAL    eigenvalue.
+* EVECT   unknown vector for the non perturbed direct flux
+* ADECT   unknown vector for the non perturbed adjoint flux
+* SOUR    fixed source
+*
+*Parameters: input/output
+* EASS    solution of the fixed source eigenvalue problem
+* TKT     CPU time spent to compute the solution of linear systems.
+* TKB     CPU time spent to compute the bilinear products.
+* ZNORM   Hotelling deflation accuracy.
+* ITER    number of iterations.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPSYS,IPFLUP
+      LOGICAL LADJ
+      INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,NADI,MAXINR,NSTART,MAXX0,
+     1 ITER
+      REAL EPS2,EPSINR,EVECT(NUN,NGRP),ADECT(NUN,NGRP),EASS(NUN*NGRP),
+     1 SOUR(NUN,NGRP),TKT,TKB,SDOT
+      DOUBLE PRECISION EVAL,ZNORM
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER  (IUNOUT=6)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: RR,QQ,VV,GAR1
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: V,H
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: G,C,S,X
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(V(NUN*NGRP,NSTART+1),G(NSTART+1),H(NSTART+1,NSTART+1),
+     1 C(NSTART+1),S(NSTART+1),X(NUN*NGRP),GAR1(NUN*NGRP))
+*----
+*  GLOBAL GMRES ITERATION.
+*----
+      ALLOCATE(RR(NUN*NGRP),QQ(NUN*NGRP),VV(NUN*NGRP))
+      
+      EPS1=EPS2*SQRT(SDOT(NUN*NGRP,SOUR,1,SOUR,1))
+      RHO=1.0E10
+      ITER=0
+      NITER=1
+      NNADI=NADI
+      DO WHILE((RHO.GT.EPS1).AND.(ITER.LT.MAXX0))
+        CALL GPTGRA(IPTRK,IPSYS,IPFLUP,LADJ,.TRUE.,LL4,ITY,NUN,NGRP,
+     1  ICL1,ICL2,IMPX,NNADI,MAXINR,EPSINR,EVAL,EVECT,ADECT,EASS(1),
+     2  SOUR(1,1),GAR1,JTER0,TKT,TKB,ZNORM,RR)
+        NITER=NITER+1
+        DO I=1,NUN*NGRP
+          X(I)=RR(I)
+        ENDDO
+        RHO=SQRT(DDOT(NUN*NGRP,X(1),1,X(1),1))
+*----
+*  TEST FOR TERMINATION ON ENTRY
+*----
+        IF(RHO.LT.EPS1) THEN
+           DEALLOCATE(VV,QQ,RR)
+           GO TO 100
+        ENDIF
+*
+        V(:NUN*NGRP,:NSTART+1)=0.0D0
+        G(:NSTART+1)=0.0D0
+        H(:NSTART+1,:NSTART+1)=0.0D0
+        C(:NSTART+1)=0.0D0
+        S(:NSTART+1)=0.0D0
+        G(1)=RHO
+        DO I=1,NUN*NGRP
+          V(I,1)=X(I)/RHO
+        ENDDO
+*----
+*  GMRES(1) ITERATION
+*----
+        K=0
+        DO WHILE((RHO.GT.EPS1).AND.(K.LT.NSTART).AND.(ITER.LT.MAXX0))
+          K=K+1
+          ITER=ITER+1
+          IF(IMPX.GT.1) WRITE(IUNOUT,300) ITER,RHO,JTER0
+          DO I=1,NUN*NGRP
+            VV(I)=REAL(V(I,K))
+            QQ(I)=0.0
+          ENDDO
+          CALL GPTGRA(IPTRK,IPSYS,IPFLUP,LADJ,.TRUE.,LL4,ITY,NUN,NGRP,
+     1    ICL1,ICL2,IMPX,NNADI,MAXINR,EPSINR,EVAL,EVECT,ADECT,VV(1),
+     2    QQ(1),GAR1,JTER,TKT,TKB,ZNORM,RR)
+          IF(JTER.NE.JTER0) CALL XABORT('GPTMRA: INCONSISTENT PRECONDIT'
+     1    //'IONING IN GMRES.')
+          NITER=NITER+1
+          DO I=1,NUN*NGRP
+            V(I,K+1)=-RR(I)
+          ENDDO
+*----
+*  MODIFIED GRAM-SCHMIDT
+*----
+          DO J=1,K
+            HR=DDOT(NUN*NGRP,V(1,J),1,V(1,K+1),1)
+            H(J,K)=HR
+            DO I=1,NUN*NGRP
+              V(I,K+1)=V(I,K+1)-HR*V(I,J)
+            ENDDO
+          ENDDO
+          H(K+1,K)=SQRT(DDOT(NUN*NGRP,V(1,K+1),1,V(1,K+1),1))
+*----
+*  REORTHOGONALIZE
+*----
+          DO J=1,K
+            HR=DDOT(NUN*NGRP,V(1,J),1,V(1,K+1),1)
+            H(J,K)=H(J,K)+HR
+            DO I=1,NUN*NGRP
+              V(I,K+1)=V(I,K+1)-HR*V(I,J)
+            ENDDO
+          ENDDO
+          H(K+1,K)=SQRT(DDOT(NUN*NGRP,V(1,K+1),1,V(1,K+1),1))
+*----
+*  WATCH OUT FOR HAPPY BREAKDOWN 
+*----
+          IF(H(K+1,K).NE.0.0) THEN
+            DO I=1,NUN*NGRP
+              V(I,K+1)=V(I,K+1)/H(K+1,K)
+            ENDDO
+          ENDIF
+*----
+*  FORM AND STORE THE INFORMATION FOR THE NEW GIVENS ROTATION
+*----
+          DO I=1,K-1
+            W1=C(I)*H(I,K)-S(I)*H(I+1,K)
+            W2=S(I)*H(I,K)+C(I)*H(I+1,K)
+            H(I,K)=W1
+            H(I+1,K)=W2
+          ENDDO
+          ZNU=SQRT(H(K,K)**2+H(K+1,K)**2)
+          IF(ZNU.NE.0.0) THEN
+            C(K)=H(K,K)/ZNU
+            S(K)=-H(K+1,K)/ZNU
+            H(K,K)=C(K)*H(K,K)-S(K)*H(K+1,K)
+            H(K+1,K)=0.0D0
+            W1=C(K)*G(K)-S(K)*G(K+1)
+            W2=S(K)*G(K)+C(K)*G(K+1)
+            G(K)=W1
+            G(K+1)=W2
+          ENDIF
+*----
+*  UPDATE THE RESIDUAL NORM
+*----
+          RHO=ABS(G(K+1))
+        ENDDO
+*----
+*  AT THIS POINT EITHER K > NSTART OR RHO < EPS1.
+*  IT'S TIME TO COMPUTE X AND CYCLE.
+*----
+        DO J=1,K
+          H(J,K+1)=G(J)
+        ENDDO
+        CALL ALSBD(K,1,H,IER,NSTART+1)
+        IF(IER.NE.0) CALL XABORT('GPTMRA: SINGULAR MATRIX.')
+        DO I=1,NUN*NGRP
+          EASS(I)=EASS(I)+REAL(DDOT(K,V(I,1),NUN*NGRP,H(1,K+1),1))
+        ENDDO
+        IF(K.EQ.NSTART) THEN
+          NNADI=NNADI+1
+          IF(IMPX.NE.0) WRITE (6,310) NNADI
+        ENDIF
+      ENDDO
+      DEALLOCATE(VV,QQ,RR)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+  100 DEALLOCATE(GAR1,X,S,C,H,G,V)
+      RETURN
+*
+  300 FORMAT(24H GPTMRA: OUTER ITERATION,I4,10H  L2 NORM=,1P,E11.4,
+     1 28H (NB. OF THERMAL ITERATIONS=,I4,1H))
+  310 FORMAT(/53H GPTMRA: INCREASING THE NUMBER OF INNER ITERATIONS TO,
+     1 I3,36H ADI ITERATIONS PER OUTER ITERATION./)
+      END