*DECK MCGBIC SUBROUTINE MCGBIC(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS, 1 IPRINT,IPTRK,IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT, 2 NZON,MATALB,M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,SOUR, 3 IAAC,ISCR,LC,LFORW,PACA,ITST,MAXI,QFR,PHIIN,CAZ0, 4 CAZ1,CAZ2,CPO,ZMU,WZMU,V,EPS,EPSI,REPSI,SIGAL,LPS, 5 NG,NGEFF,NGIND,NCONV,LNCONV,NLIN,NFUNL,KEYFLX,KEYCUR, 6 STIS,NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM, 7 IDIR,NBATCH) * *----------------------------------------------------------------------- * *Purpose: * Solve the linear system obtained by the characteristics formalism * with BiCGSTAB iterative approach. * *Copyright: * Copyright (C) 2002 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 * *Parameters: * SUBFFI flux integration subroutine with isotropic source. * SUBFFA flux integration subroutine with anisotropic source. * SUBDLC 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. * K total number of volumes for which specific values * of the neutron flux and reactions rates are required. * N total number of unknowns in vectors SUNKNO and FUNKNO. * 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 plan. * NREG number of regions (volumes). * NSOUT number of outer surfaces. * SOUR scratch. * 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. * V volumes. * EPS precision reached after min(MAXI,ITST) iterations. * EPSI 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. * SIGAL total cross-section and albedo array. * LPS used in scr acceleration. * NG 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 FUNKNO vector. * KEYCUR position of current elements in FUNKNO 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 NONE *---- * SUBROUTINE ARGUMENTS *---- TYPE(C_PTR) KPSYS(NGEFF),IPTRK,IPMACR INTEGER NGEFF,IPRINT,IFTRAK,NDIM,K,N,NLONG,NG,NZON(NLONG),M, 1 NANI,NMU,NMAX,IAAC,LC,PACA,ITST(NGEFF),MAXI,NANGL,NREG,NSOUT, 2 ISCR,LPS,NGIND(NGEFF),LNCONV,NLIN,NFUNL,KEYFLX(NREG,NLIN,NFUNL), 3 KEYCUR(NLONG-NREG),STIS,NPJJM,MATALB(-NSOUT:NREG),N2REG,N2SOU, 4 NZP,IDIR,NBATCH REAL QFR(N,NGEFF),PHIIN(N,NGEFF),CPO(NMU),ZMU(NMU),WZMU(NMU), 1 V(NLONG),EPS(NGEFF),EPSI,REPSI(MAXI,NGEFF),SIGAL(-6:M,NGEFF), 2 DELU,FACSYM DOUBLE PRECISION CAZ0(NANGL),CAZ1(NANGL),CAZ2(NANGL), 1 PHIOUT(N,NGEFF),SOUR(N,NGEFF) LOGICAL LFORW,CYCLIC,NCONV(NGEFF),REBFLG,LPRISM EXTERNAL SUBFFI,SUBFFA,SUBDLC,SUBSCH *--- * LOCAL VARIABLES *--- INTEGER J,II,ITER,MAXINT REAL R,BI,WI,RT1,EPSINT,REPSMAX REAL SDOT DOUBLE PRECISION DDOT LOGICAL RHSFLG INTRINSIC SQRT,ABS *---- * ALLOCATABLE ARRAYS *---- REAL, ALLOCATABLE, DIMENSION(:,:) :: RHS,PI,RI,SI,ROT,API,AUX *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(RHS(N,NGEFF),PI(N,NGEFF),RI(N,NGEFF),SI(N,NGEFF), 1 ROT(N,NGEFF),API(N,NGEFF),AUX(2,NGEFF)) *--- IF(MAXI.LT.4) CALL XABORT('MCGBIC: MAXI MUST BE >= 4.') MAXINT=MAXI-1 EPSINT=0.01*EPSI RHSFLG=.TRUE. *--- ITER=1 * a first iteration CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT,IPTRK, 1 IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB,M,NANI,NMU, 2 NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC,ISCR,LC,LFORW, 3 PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR,QFR,PHIIN(1,1), 4 CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV,.FALSE.,STIS, 5 NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH) DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N RI(J,II)=REAL(PHIOUT(J,II))-PHIIN(J,II) PHIIN(J,II)=REAL(PHIOUT(J,II)) ENDDO R=SDOT(N,RI(1,II),1,RI(1,II),1) REPSI(ITER,II)=REAL(SQRT(R/DDOT(N,PHIOUT(1,II),1, 1 PHIOUT(1,II),1))) IF (REPSI(ITER,II).LE.EPSI) THEN NCONV(II)=.FALSE. ITST(II)=ITER EPS(II)=REPSI(ITER,II) LNCONV=LNCONV-1 ENDIF IF (LNCONV.EQ.0) RETURN ENDIF ENDDO *--- 10 ITER=ITER+1 * compute initial residual vector CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT,IPTRK, 1 IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB,M,NANI,NMU, 2 NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC,ISCR,LC,LFORW, 3 PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR,QFR,PHIIN(1,1), 4 CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV,.FALSE.,STIS, 5 NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH) DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N RI(J,II)=REAL(PHIOUT(J,II))-PHIIN(J,II) ENDDO R=SDOT(N,RI(1,II),1,RI(1,II),1) REPSI(ITER,II)=REAL(SQRT(R/DDOT(N,PHIOUT(1,II),1, 1 PHIOUT(1,II),1))) IF (REPSI(ITER,II).LE.EPSI) THEN NCONV(II)=.FALSE. ITST(II)=ITER EPS(II)=REPSI(ITER,II) LNCONV=LNCONV-1 DO J=1,N PHIIN(J,II)=REAL(PHIOUT(J,II)) ENDDO ENDIF IF (LNCONV.EQ.0) RETURN DO J=1,N PI(J,II)=RI(J,II) ROT(J,II)=RI(J,II) ENDDO * RT2=R !!SDOT(N,RI,1,ROT,1) AUX(1,II)=R !!SDOT(N,RI(1,II),1,ROT(1,II),1) ENDIF ENDDO *--- IF (RHSFLG) THEN * evaluate RHS of the linear system IF (IPRINT.GT.3) THEN WRITE(6,100) IAAC,ISCR ENDIF RHS(:N,:NGEFF)=0.0 CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT, 1 IPTRK,IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB, 2 M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC, 2 ISCR,LC,LFORW,PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR, 3 QFR,RHS(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS, 4 NCONV,.FALSE.,STIS,NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP, 5 DELU,FACSYM,IDIR,NBATCH) DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N RHS(J,II)=REAL(PHIOUT(J,II)) ENDDO ENDIF ENDDO RHSFLG=.FALSE. ENDIF * CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT,IPTRK, 1 IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB,M,NANI,NMU, 2 NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC,ISCR,LC,LFORW, 3 PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR,QFR,PI(1,1),CAZ0, 4 CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV,.FALSE.,STIS,NPJJM, 5 REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH) ITER=ITER+1 DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N API(J,II)=PI(J,II)-REAL(PHIOUT(J,II))+RHS(J,II) ENDDO REPSI(ITER,II)=REPSI((ITER-1),II) ENDIF ENDDO * DO WHILE (ITER.LT.(MAXI-1)) * BiCGSTAB iterations ITER=ITER+1 * DO II=1,NGEFF IF (NCONV(II)) THEN AUX(2,II)=AUX(1,II)/SDOT(N,API(1,II),1,ROT(1,II),1) DO J=1,N SI(J,II)=RI(J,II)-AUX(2,II)*API(J,II) ENDDO ENDIF ENDDO * CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT, 1 IPTRK,IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB, 2 M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC, 3 ISCR,LC,LFORW,PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR, 4 QFR,SI(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV, 5 .FALSE.,STIS,NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU, 6 FACSYM,IDIR,NBATCH) REPSMAX=0.0 DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N RI(J,II)=SI(J,II)-REAL(PHIOUT(J,II))+RHS(J,II) ENDDO WI=SDOT(N,RI(1,II),1,SI(1,II),1)/ 1 SDOT(N,RI(1,II),1,RI(1,II),1) DO J=1,N PHIIN(J,II)=PHIIN(J,II)+AUX(2,II)*PI(J,II)+WI*SI(J,II) RI(J,II)=SI(J,II)-WI*RI(J,II) ENDDO R=SDOT(N,RI(1,II),1,RI(1,II),1) REPSI(ITER,II)=SQRT(R/SDOT(N,PHIIN(1,II),1, 1 PHIIN(1,II),1)) REPSMAX=MAX(REPSMAX,REPSI(ITER,II)) IF (REPSI(ITER,II).LE.EPSI) THEN NCONV(II)=.FALSE. ITST(II)=ITER EPS(II)=REPSI(ITER,II) LNCONV=LNCONV-1 ENDIF IF (LNCONV.EQ.0) GO TO 20 RT1=AUX(1,II) AUX(1,II)=SDOT(N,RI(1,II),1,ROT(1,II),1) BI=AUX(1,II)/RT1*AUX(2,II)/WI DO J=1,N PI(J,II)=RI(J,II)+BI*(PI(J,II)-WI*API(J,II)) ENDDO ENDIF ENDDO * CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT, 1 IPTRK,IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB, 2 M,NANI,NMU,NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,IAAC, 3 ISCR,LC,LFORW,PACA,EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR, 4 QFR,PI(1,1),CAZ0,CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV, 5 .FALSE.,STIS,NPJJM,REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU, 6 FACSYM,IDIR,NBATCH) ITER=ITER+1 DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N API(J,II)=PI(J,II)-REAL(PHIOUT(J,II))+RHS(J,II) ENDDO REPSI(ITER,II)=REPSI((ITER-1),II) ENDIF ENDDO ENDDO * 20 CONTINUE * determine final residual norm ITER=ITER+1 DO II=1,NGEFF IF (NCONV(II)) THEN ITST(II)=ITER ELSE IF (ITST(II).NE.1) THEN NCONV(II)=.TRUE. ITST(II)=ITST(II)+1 ENDIF ENDIF ENDDO CALL MCGFL1(SUBFFI,SUBFFA,SUBDLC,SUBSCH,CYCLIC,KPSYS,IPRINT,IPTRK, 1 IFTRAK,IPMACR,NDIM,K,N,NLONG,PHIOUT,NZON,MATALB,M,NANI,NMU, 2 NMAX,NANGL,NREG,NSOUT,NG,NGEFF,NGIND,SOUR,0,0,LC,LFORW,PACA, 3 EPSINT,MAXINT,NLIN,NFUNL,KEYFLX,KEYCUR,QFR,PHIIN(1,1),CAZ0, 4 CAZ1,CAZ2,CPO,ZMU,WZMU,V,SIGAL,LPS,NCONV,.FALSE.,STIS,NPJJM, 5 REBFLG,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,IDIR,NBATCH) DO II=1,NGEFF IF (NCONV(II)) THEN DO J=1,N RI(J,II)=REAL(PHIOUT(J,II))-PHIIN(J,II) ENDDO R=SDOT(N,RI(1,II),1,RI(1,II),1) REPSI(ITST(II),II)=SQRT(R/SDOT(N,PHIIN(1,II),1, 1 PHIIN(1,II),1)) DO J=1,N PHIIN(J,II)=REAL(PHIOUT(J,II)) ENDDO EPS(II)=REPSI(ITST(II),II) ENDIF ENDDO LNCONV=0 IF (ITER.LT.MAXI) THEN DO II=1,NGEFF IF (EPS(II).GT.EPSI) THEN IF ((IAAC.GT.0).OR.(ISCR.GT.0)) THEN IAAC=0 ISCR=0 RHSFLG=.TRUE. ENDIF IF (IPRINT.GT.2) WRITE(6,200) ITER NCONV(II)=.TRUE. LNCONV=LNCONV+1 ELSE NCONV(II)=.FALSE. ENDIF ENDDO IF (LNCONV.GT.0) GO TO 10 ENDIF *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(AUX,API,ROT,SI,RI,PI,RHS) RETURN * 100 FORMAT(31H RHS CALCULATED WITH AAC-SCR : ,I1,1H-,I1) 200 FORMAT(37H WARNING : BAD PREVISION, RESTART AT ,I4,10H ITERATION) END