diff options
Diffstat (limited to 'Dragon/src/SNGMRE.f')
| -rw-r--r-- | Dragon/src/SNGMRE.f | 229 |
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 |