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