diff options
Diffstat (limited to 'Dragon/src/MOCFCF.f')
| -rw-r--r-- | Dragon/src/MOCFCF.f | 412 |
1 files changed, 412 insertions, 0 deletions
diff --git a/Dragon/src/MOCFCF.f b/Dragon/src/MOCFCF.f new file mode 100644 index 0000000..fdf9515 --- /dev/null +++ b/Dragon/src/MOCFCF.f @@ -0,0 +1,412 @@ +*DECK MOCFCF + SUBROUTINE MOCFCF(SUBFFI,SUBFFA,SUBLDC,SUBSCH,IFTRAK,NBTR,MXSUB, + 1 MXSEG,NDIM,KPN,NREG,NSOUT,NMAT,NALB,NGEFF,NPHI, + 2 NGSS,NLF,NFUNL,NMOD,NLFX,NLIN,NFUNLX,KEYFLX, + 3 MATALB,NCONV,SIGANG,CAZ1,CAZ2,XGSS,YGSS,WGSS, + 4 SOUR,ISGNR,IDIR,NBATCH,PHIOUT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Flux integration upon the cyclic tracking. +* +*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. Roy and R. Le Tellier +* +*Parameters: input +* SUBFFI Isotropic flux integration subroutine. +* SUBFFA Anisotropic flux integration subroutine. +* SUBLDC flux integration subroutine with linear-discontinuous source. +* SUBSCH Track coefficients calculation subroutine. +* IFTRAK tracking file unit number. +* NGEFF number of groups to process. +* NMAT number of mixtures. +* NLF number of Legendre orders for the flux. +* NREG number of regions. +* KPN number of unknowns per energy group including spherical +* harmonic terms and fundamental currents. +* NGSS number of polar angles. +* NSOUT number of surfaces. +* MXSUB maximun number of subtracks in a track. +* MXSEG maximun number of segments in a track. +* NBTR number of tracks. +* NPHI number of angles in the plane. +* NFUNL number of moments of the flux (in 2D: NFUNL=NLF*(NLF+1)/2). +* NMOD first dimension of ISGNR. +* NLFX scattering anisotropy used to compute spherical harmonics. +* NLIN linear discontinuous flag (=1 SC/DD0; =3 LDC/DD1). +* NFUNLX number of spherical harmonics components. +* ISGNR array of spherical harmonics signs. +* KEYFLX position of flux elements in PHIIN vector. +* MATALB mixture and albedo indices. +* WGSS polar weights. +* XGSS polar angle cosines. +* YGSS polar angle sines. +* NALB number of albedos. +* SIGANG arrays of total cross-sections and albedos. +* CAZ1 first cosines of the different tracking azimuthal angles. +* CAZ2 second cosines of the different tracking azimuthal angles. +* SOUR total source vector components. +* NCONV logical array of convergence status for each group (.TRUE. for +* not converged). +* IDIR direction of fundamental current for TIBERE with MoC +* (=0,1,2,3). +* NBATCH number of tracks processed in each OpenMP core (default: =1). +* +*Parameters: output +* PHIOUT vector containing the zonal flux moments. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IFTRAK,NGEFF,NMAT,NLF,NREG,NDIM,KPN,NGSS,NSOUT,MXSUB, + 1 MXSEG,NBTR,NPHI,NFUNL,NMOD,NLFX,NLIN,NFUNLX,ISGNR(NMOD,NFUNLX), + 2 KEYFLX(NREG,NLIN,NFUNL),MATALB(-NSOUT:NREG),NALB,IDIR,NBATCH + REAL WGSS(NGSS),XGSS(NGSS),YGSS(NGSS), + 1 SIGANG(-NALB:NMAT,NGEFF) + DOUBLE PRECISION CAZ1(NPHI),CAZ2(NPHI),SOUR(KPN,NGEFF), + 1 PHIOUT(KPN,NGEFF) + LOGICAL NCONV(NGEFF) + EXTERNAL SUBFFI,SUBFFA,SUBLDC,SUBSCH +*---- +* LOCAL VARIABLES +*---- + INTEGER MXE + PARAMETER (MXE=64) + INTEGER ILINE,IANG,JANG,ISUB,I,IE,II,NOMI,NZI,IND,JF,INDX, + 1 INDY,IREG,I0,IL1,IBATCH,NDFUNLX + DOUBLE PRECISION DWEIG(MXE),Q0,Q1,Q0X,Q1X,Q0Y,Q1Y + LOGICAL LNEW +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: NSUB,NSEG + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NRSEG,KANGL + REAL, ALLOCATABLE, DIMENSION(:,:) :: RHARM + REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: TRHAR + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DSIG,WEIGHT,EXPT, + 1 EXP2,FLUX + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: SEGLEN,COEFI, + 1 OMEGAX,OMEGAY,OMG2,FLM,FLP,CYM,CYP,DFLM,DFLP,CYM2,CYP2,FLUV, + 2 DFLUV + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: PHIV,DPHIV +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(NSUB(NBATCH),NSEG(NBATCH),WEIGHT(NBATCH), + 1 KANGL(MXSUB,NBATCH),NRSEG(MXSEG,NBATCH),SEGLEN(MXSEG,NBATCH)) + ALLOCATE(FLM(NGSS,MXSEG),FLP(NGSS,MXSEG),CYM(NGSS,MXSEG), + 1 CYP(NGSS,MXSEG),EXPT(NGSS*MXSEG)) + ALLOCATE(OMEGAX(NPHI,NGSS),OMEGAY(NPHI,NGSS),OMG2(NGSS,3)) +*--- +* Compute flux and currents for this tracking line +*--- + PHIOUT(:KPN,:NGEFF)=0.0D0 + IF((NLF.EQ.1).AND.(NLIN.EQ.1)) THEN +*---- +* ISOTROPIC SCATTERING +*---- + ALLOCATE(FLUX(KPN),EXP2(2*NGSS*MXSEG)) + DO IE=1,NGSS + DO IANG=1,NPHI + OMEGAX(IANG,IE)=CAZ1(IANG)/YGSS(IE) + OMEGAY(IANG,IE)=CAZ2(IANG)/YGSS(IE) + ENDDO + ENDDO + DO IBATCH=1,(NBTR-1)/NBATCH+1 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + READ(IFTRAK) NSUB(IL1),NSEG(IL1),WEIGHT(IL1), + 1 (KANGL(I,IL1),I=1,NSUB(IL1)),(NRSEG(I,IL1),I=1,NSEG(IL1)), + 2 (SEGLEN(I,IL1),I=1,NSEG(IL1)) + IF(NSUB(IL1).GT.MXSUB) CALL XABORT('MOCFCF: MXSUB OVERFLOW.') + ENDDO +*$OMP PARALLEL DO +*$OMP1 PRIVATE(IL1,IE,DWEIG,ISUB,LNEW,NOMI,NZI,IND,FLM,FLP,OMG2) +*$OMP2 PRIVATE(FLUX,I0,ILINE,EXPT,EXP2,CYM,CYP) + DO II=1,NGEFF + IF(NCONV(II)) THEN + FLUX(:KPN)=0.0D0 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + OMG2(:NGSS,:3)=0.0D0 + DO IE=1,NGSS + DWEIG(IE)=WEIGHT(IL1)*WGSS(IE) + ENDDO + FLM(:NGSS,:MXSEG)=0.D0 + FLP(:NGSS,:MXSEG)=0.D0 + ISUB=0 + LNEW=.TRUE. + DO I0=1,NSEG(IL1) + NOMI=NRSEG(I0,IL1) + NZI=MATALB(NOMI) + IF(NZI.LE.0) THEN + LNEW=.TRUE. + ELSE + IF(LNEW) THEN + ISUB=ISUB+1 + LNEW=.FALSE. + ENDIF + IND=KEYFLX(NOMI,1,1) + IF(IDIR.EQ.0) THEN + DO IE=1,NGSS + FLM(IE,I0)=DWEIG(IE)*SOUR(IND,II) + FLP(IE,I0)=FLM(IE,I0) + ENDDO + ELSEIF(IDIR.EQ.1) THEN + DO IE=1,NGSS + OMG2(IE,1)=3.0D0*OMEGAX(KANGL(ISUB,IL1),IE)**2 + FLM(IE,I0)=DWEIG(IE)*SOUR(IND,II)*OMG2(IE,1) + FLP(IE,I0)=FLM(IE,I0) + ENDDO + ELSEIF(IDIR.EQ.2) THEN + DO IE=1,NGSS + OMG2(IE,2)=3.0D0*OMEGAY(KANGL(ISUB,IL1),IE)**2 + FLM(IE,I0)=DWEIG(IE)*SOUR(IND,II)*OMG2(IE,2) + FLP(IE,I0)=FLM(IE,I0) + ENDDO + ELSEIF(IDIR.EQ.3) THEN + DO IE=1,NGSS + OMG2(IE,3)=3.0D0*(1.0-1.0/YGSS(IE)**2) + FLM(IE,I0)=DWEIG(IE)*SOUR(IND,II)*OMG2(IE,3) + FLP(IE,I0)=FLM(IE,I0) + ENDDO + ENDIF + ENDIF + ENDDO +* MOCFFIR: 'Source Term Isolation' Strategy turned on +* MOCFFIS: 'Source Term Isolation' Strategy turned off +* MOCFFIT: 'MOCC/MCI' Iterative Strategy + CALL SUBFFI(SUBSCH,NREG,NSOUT,KPN,NMAT,NSEG(IL1), + 1 SEGLEN(1,IL1),NRSEG(1,IL1),NGSS,MATALB, + 2 SIGANG(-NALB,II),KEYFLX,YGSS,FLUX,EXPT, + 3 EXP2,FLM,FLP,CYM,CYP,IDIR,OMG2) + ENDDO ! ILINE + PHIOUT(:KPN,II)=PHIOUT(:KPN,II)+FLUX(:KPN) + ENDIF + ENDDO ! II +*$OMP END PARALLEL DO + ENDDO ! IBATCH + DEALLOCATE(OMG2,OMEGAY,OMEGAX) + DEALLOCATE(EXP2,FLUX) + ELSE IF(NLIN.EQ.1) THEN +*---- +* ANISOTROPIC SCATTERING +*---- + ALLOCATE(FLUX(KPN),EXP2(2*NGSS*MXSEG)) + ALLOCATE(TRHAR(NGSS,NFUNL,NPHI,2),RHARM(NGSS,NFUNL)) + DO IANG=1,NPHI + CALL MOCCHR(2,NLF-1,NFUNL,NGSS,XGSS,CAZ1(IANG),CAZ2(IANG), + 1 RHARM) + DO 15 JF=1,NFUNL + DO 10 IE=1,NGSS + TRHAR(IE,JF,IANG,1)=ISGNR(1,JF)*RHARM(IE,JF) + TRHAR(IE,JF,IANG,2)=ISGNR(NMOD,JF)*RHARM(IE,JF) + 10 CONTINUE + 15 CONTINUE + ENDDO + DEALLOCATE(RHARM) + DO IBATCH=1,(NBTR-1)/NBATCH+1 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + READ(IFTRAK) NSUB(IL1),NSEG(IL1),WEIGHT(IL1), + 1 (KANGL(I,IL1),I=1,NSUB(IL1)),(NRSEG(I,IL1),I=1,NSEG(IL1)), + 1 (SEGLEN(I,IL1),I=1,NSEG(IL1)) + IF(NSUB(IL1).GT.MXSUB) CALL XABORT('MOCFCF: MXSUB OVERFLOW.') + ENDDO +*$OMP PARALLEL DO +*$OMP1 PRIVATE(IL1,IE,DWEIG,ISUB,LNEW,NOMI,NZI,Q0,Q1,FLM,FLP) +*$OMP2 PRIVATE(FLUX,I0,II,EXPT,EXP2,CYM,CYP) + DO II=1,NGEFF + IF(NCONV(II)) THEN + FLUX(:KPN)=0.0D0 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + DO IE=1,NGSS + DWEIG(IE)=WEIGHT(IL1)*WGSS(IE) + ENDDO + FLM(:NGSS,:MXSEG)=0.0D0 + FLP(:NGSS,:MXSEG)=0.0D0 + ISUB=0 + LNEW=.TRUE. + DO I0=1,NSEG(IL1) + NOMI=NRSEG(I0,IL1) + NZI=MATALB(NOMI) + IF(NZI.LE.0) THEN + LNEW=.TRUE. + ELSE + IF(LNEW) THEN + ISUB=ISUB+1 + LNEW=.FALSE. + ENDIF + DO IE=1,NGSS + Q0=0.D0 + Q1=0.D0 + DO JF=1,NFUNL + IND=KEYFLX(NOMI,1,JF) + Q0=Q0+SOUR(IND,II)*TRHAR(IE,JF,KANGL(ISUB,IL1),2) + Q1=Q1+SOUR(IND,II)*TRHAR(IE,JF,KANGL(ISUB,IL1),1) + ENDDO + FLM(IE,I0)=DWEIG(IE)*Q0 + FLP(IE,I0)=DWEIG(IE)*Q1 + ENDDO + ENDIF + ENDDO + IF(ISUB.NE.NSUB(IL1)) CALL XABORT('MOCFCF: NSUB ERROR.') +* MOCFFAR: 'Source Term Isolation' Strategy turned on +* MOCFFAS: 'Source Term Isolation' Strategy turned off +* MOCFFAT: 'MOCC/MCI' Iterative Strategy + CALL SUBFFA(SUBSCH,NREG,NSOUT,KPN,NMAT,NSEG(IL1), + 1 SEGLEN(1,IL1),NRSEG(1,IL1),NGSS,NFUNL,MATALB, + 2 SIGANG(-NALB,II),KEYFLX,YGSS,FLUX,EXPT,EXP2,FLM, + 3 FLP,CYM,CYP,NPHI,NSUB(IL1),KANGL(1,IL1),TRHAR) + ENDDO ! ILINE + PHIOUT(:KPN,II)=PHIOUT(:KPN,II)+FLUX(:KPN) + ENDIF + ENDDO ! II +*$OMP END PARALLEL DO + ENDDO ! IBATCH + DEALLOCATE(EXP2,FLUX,TRHAR) + ELSE IF(NLIN.EQ.3) THEN +*---- +* LINEAR DISCONTINUOUS SOURCE APPROXIMATION +*---- + NDFUNLX=NDIM*NFUNLX + ALLOCATE(PHIV(NFUNLX,NREG,NGEFF),DPHIV(NDFUNLX,NREG,NGEFF)) + ALLOCATE(FLUV(NFUNLX,NREG),DFLUV(NDFUNLX,NREG)) + ALLOCATE(EXP2(5*NGSS*MXSEG)) + ALLOCATE(TRHAR(NGSS,NFUNLX,NPHI,2),RHARM(NGSS,NFUNLX)) + DO IANG=1,NPHI + CALL MOCCHR(2,NLFX-1,NFUNLX,NGSS,XGSS,CAZ1(IANG),CAZ2(IANG), + 1 RHARM) + DO 25 JF=1,NFUNLX + DO 20 IE=1,NGSS + TRHAR(IE,JF,IANG,1)=ISGNR(1,JF)*RHARM(IE,JF) + TRHAR(IE,JF,IANG,2)=ISGNR(NMOD,JF)*RHARM(IE,JF) + 20 CONTINUE + 25 CONTINUE + ENDDO + DEALLOCATE(RHARM) + DO II=1,NGEFF + IF(NCONV(II)) THEN + PHIV(:NFUNLX,:NREG,II)=0.0D0 + DPHIV(:NDFUNLX,:NREG,II)=0.0D0 + ENDIF + ENDDO + ALLOCATE(DFLM(NGSS,MXSEG),DFLP(NGSS,MXSEG),DSIG(MXSEG), + 1 CYM2(NGSS,MXSEG),CYP2(NGSS,MXSEG)) + DO IBATCH=1,(NBTR-1)/NBATCH+1 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + READ(IFTRAK) NSUB(IL1),NSEG(IL1),WEIGHT(IL1), + 1 (KANGL(I,IL1),I=1,NSUB(IL1)),(NRSEG(I,IL1),I=1,NSEG(IL1)), + 1 (SEGLEN(I,IL1),I=1,NSEG(IL1)) + IF(NSUB(IL1).GT.MXSUB) CALL XABORT('MOCFCF: MXSUB OVERFLOW.') + ENDDO +*$OMP PARALLEL DO +*$OMP1 PRIVATE(IL1,IE,DWEIG,ISUB,JANG,LNEW,NOMI,NZI,Q0,Q1,Q0X,Q1X,Q0Y) +*$OMP2 PRIVATE(Q1Y,IND,INDX,INDY,FLM,FLP,DFLM,DFLP,FLUV,DFLUV,I0,ILINE) +*$OMP3 PRIVATE(DSIG,EXPT,EXP2,CYM,CYP,CYM2,CYP2) + DO II=1,NGEFF + IF(NCONV(II)) THEN + FLUV(:NFUNLX,:NREG)=0.0D0 + DFLUV(:NDFUNLX,:NREG)=0.0D0 + DO ILINE=(IBATCH-1)*NBATCH+1,MIN(IBATCH*NBATCH,NBTR) + IL1=ILINE-(IBATCH-1)*NBATCH + DO IE=1,NGSS + DWEIG(IE)=WEIGHT(IL1)*WGSS(IE) + ENDDO + FLM(:NGSS,:MXSEG)=0.0D0 + FLP(:NGSS,:MXSEG)=0.0D0 + DFLM(:NGSS,:MXSEG)=0.0D0 + DFLP(:NGSS,:MXSEG)=0.0D0 + ISUB=0 + JANG=0 + LNEW=.TRUE. + DO I0=1,NSEG(IL1) + NOMI=NRSEG(I0,IL1) + NZI=MATALB(NOMI) + IF(NZI.LE.0) THEN + LNEW=.TRUE. + ELSE + IF(LNEW) THEN + ISUB=ISUB+1 + JANG=KANGL(ISUB,IL1) + LNEW=.FALSE. + ENDIF + DO IE=1,NGSS + Q0=0.D0 + Q1=0.D0 + Q0X=0.0D0 + Q1X=0.0D0 + Q0Y=0.0D0 + Q1Y=0.0D0 + DO JF=1,NFUNL + IND=KEYFLX(NOMI,1,JF) + INDX=KEYFLX(NOMI,2,JF) + INDY=KEYFLX(NOMI,3,JF) + Q0=Q0+SOUR(IND,II)*TRHAR(IE,JF,JANG,2) + Q1=Q1+SOUR(IND,II)*TRHAR(IE,JF,JANG,1) + Q0X=Q0X+SOUR(INDX,II)*TRHAR(IE,JF,JANG,2) + Q1X=Q1X+SOUR(INDX,II)*TRHAR(IE,JF,JANG,1) + Q0Y=Q0Y+SOUR(INDY,II)*TRHAR(IE,JF,JANG,2) + Q1Y=Q1Y+SOUR(INDY,II)*TRHAR(IE,JF,JANG,1) + ENDDO + FLM(IE,I0)=Q0 + FLP(IE,I0)=Q1 + DFLM(IE,I0)=-Q0X*CAZ1(JANG)-Q0Y*CAZ2(JANG) + DFLP(IE,I0)=Q1X*CAZ1(JANG)+Q1Y*CAZ2(JANG) + ENDDO + ENDIF + ENDDO + IF(ISUB.NE.NSUB(IL1)) CALL XABORT('MOCFCF: NSUB ERROR.') +* MOCFFAL: 'Source Term Isolation' Strategy turned off + CALL SUBLDC(SUBSCH,NREG,NSOUT,NMAT,NSEG(IL1), + 1 SEGLEN(1,IL1),NRSEG(1,IL1),NGSS,NFUNLX,MATALB, + 2 DWEIG,SIGANG(-NALB,II),YGSS,FLM,FLP,DFLM,DFLP, + 3 NPHI,NSUB(IL1),KANGL(1,IL1),TRHAR,FLUV, + 4 DFLUV,DSIG,EXPT,EXP2,CYM,CYP,CYM2,CYP2) + ENDDO ! ILINE + PHIV(:NFUNLX,:NREG,II)=PHIV(:NFUNLX,:NREG,II)+ + 1 FLUV(:NFUNLX,:NREG) + DPHIV(:NDFUNLX,:NREG,II)=DPHIV(:NDFUNLX,:NREG,II)+ + 1 DFLUV(:NDFUNLX,:NREG) + ENDIF + ENDDO ! II +*$OMP END PARALLEL DO + ENDDO ! IBATCH + DEALLOCATE(CYP2,CYM2,DSIG,DFLP,DFLM) + DEALLOCATE(EXP2,TRHAR) + ALLOCATE(COEFI(2*NFUNLX,2*NFUNLX)) + CALL MCGCOEF(NFUNLX,NGSS,YGSS,WGSS,NPHI,CAZ1,CAZ2,COEFI) + DO II=1,NGEFF + IF(NCONV(II)) THEN + DO IREG=1,NREG + DPHIV(:,IREG,II)=MATMUL(COEFI,DPHIV(:,IREG,II)) + DO JF=1,NFUNL + PHIOUT(KEYFLX(IREG,1,JF),II)=PHIV(JF,IREG,II) + PHIOUT(KEYFLX(IREG,2,JF),II)=DPHIV(JF,IREG,II) + PHIOUT(KEYFLX(IREG,3,JF),II)=DPHIV(NFUNLX+JF,IREG,II) + ENDDO + ENDDO + ENDIF + ENDDO + DEALLOCATE(COEFI,DFLUV,FLUV,DPHIV,PHIV) + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(EXPT,CYP,CYM,FLP,FLM) + DEALLOCATE(SEGLEN,NRSEG,KANGL,WEIGHT,NSEG,NSUB) +* + RETURN + END |