diff options
Diffstat (limited to 'Trivac/src/NSSEIG.f')
| -rwxr-xr-x | Trivac/src/NSSEIG.f | 572 |
1 files changed, 572 insertions, 0 deletions
diff --git a/Trivac/src/NSSEIG.f b/Trivac/src/NSSEIG.f new file mode 100755 index 0000000..b82fba8 --- /dev/null +++ b/Trivac/src/NSSEIG.f @@ -0,0 +1,572 @@ +*DECK NSSEIG + SUBROUTINE NSSEIG(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL, + > VOL,MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,CHI,SIGF,SCAT,A11X,A11Y, + > A11Z,EPSTHR,MAXTHR,NADI,EPSOUT,MAXOUT,ICL1,ICL2,ITER,EVECT, + > FKEFF,IMPX) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solution of a multigroup eigenvalue system for the calculation of the +* direct neutron flux in Trivac. Use the preconditioned power method +* with a two-parameter SVAT acceleration technique. CMFD solution. +* +*Copyright: +* Copyright (C) 2023 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 +* NMAX first dimension of array A11X. +* NMAY first dimension of array A11Y. +* NMAZ first dimension of array A11Z. +* LL4F number of unknowns per energy group. +* NDIM number of dimensions (1, 2 or 3). +* NEL number of nodes. +* NMIX number of mixtures in the nodal calculation. +* NG number of energy groups. +* MAT material mixtures. +* IDL position of averaged fluxes in unknown vector. +* VOL node volumes. +* MUX X-oriented compressed storage mode indices. +* MUY Y-oriented compressed storage mode indices. +* MUZ Z-oriented compressed storage mode indices. +* IMAX X-oriented position of each first non-zero column element. +* IMAY Y-oriented position of each first non-zero column element. +* IMAZ Z-oriented position of each first non-zero column element. +* IPY Y-oriented permutation matrices. +* IPZ Z-oriented permutation matrices. +* CHI fission spectra. +* SIGF nu times fission cross section. +* SCAT scattering cross section. +* A11X X-oriented sparse coefficient matrix. +* A11Y Y-oriented sparse coefficient matrix. +* A11Z Z-oriented sparse coefficient matrix. +* EPSTHR thermal iteration epsilon. +* MAXTHR maximum number of thermal iterations. +* NADI number of inner ADI iterations. +* EPSOUT convergence epsilon for the power method. +* MAXOUT maximum number of iterations for the power method. +* ICL1 number of free iretations in one cycle of the up-scattering +* iterations. +* ICL2 number of accelerated up-scattering iterations in one cycle. +* EVECT initial estimate of fundamental eigenvalue. +* IMPX print parameter. +* FKEFF initial estimate of fundamental eigenvalue. +* +*Parameters: output +* ITER number of iterations. +* EVECT corresponding eigenvector. +* FKEFF fundamental eigenvalue. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER, INTENT(IN) :: NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG, + > MAT(NEL),IDL(NEL),MUX(LL4F),MUY(LL4F),MUZ(LL4F),IMAX(LL4F), + > IMAY(LL4F),IMAZ(LL4F),IPY(LL4F),IPZ(LL4F),MAXTHR,NADI,MAXOUT, + > ICL1,ICL2,IMPX + REAL, INTENT(IN) :: VOL(NEL),CHI(NMIX,NG),SIGF(NMIX,NG), + > SCAT(NMIX,NG,NG),A11X(NMAX,NG),A11Y(NMAY,NG),A11Z(NMAZ,NG) + INTEGER, INTENT(OUT) :: ITER + REAL, INTENT(IN) :: EPSTHR,EPSOUT + REAL, INTENT(INOUT) :: EVECT(LL4F,NG),FKEFF +*---- +* LOCAL VARIABLES +*---- + REAL, PARAMETER :: EPS1=1.0E-5 + REAL(KIND=8), PARAMETER :: ALP_TAB(24) = (/ 0.2, 0.4, 0.6, + 1 0.8, 1.0, 1.2, 1.5, 2.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, + 2 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0 /) + REAL(KIND=8), PARAMETER :: BET_TAB(11) = (/ -1.0, -0.8, -0.6, + 1 -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 /) + REAL(KIND=8) :: AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,BEBE,BEBG,BEBH, + 1 BGBG,BGBH,BHBH,AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH, + 2 X,DXDA,DXDB,Y,DYDA,DYDB,Z,DZDA,DZDB,F,D2F(2,3),EVAL,ALP,BET, + 3 FMIN,VVV + LOGICAL LOGTES + CHARACTER(LEN=3) :: TEXT3 +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, ALLOCATABLE, DIMENSION(:) :: S2,F1,GARM1,GARM2 + REAL, ALLOCATABLE, DIMENSION(:,:) :: S,GRAD1,GRAD2,GAR1,GAR2, + 1 GAR3,GAF1,GAF2,GAF3 + REAL, ALLOCATABLE, DIMENSION(:,:) :: IA11X,IA11Y,IA11Z + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WORK +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(S2(LL4F),S(LL4F,NG),GRAD1(LL4F,NG),GRAD2(LL4F,NG), + > GAR1(LL4F,NG),GAR2(LL4F,NG),GAR3(LL4F,NG),GAF1(LL4F,NG), + > GAF2(LL4F,NG),GAF3(LL4F,NG),WORK(LL4F,NG,3),IA11X(NMAX,NG), + > IA11Y(NMAY,NG),IA11Z(NMAZ,NG)) +*---- +* LU MATRIX FACTORIZATION +*---- + IA11X(:NMAX,:NG)=A11X(:NMAX,:NG) + DO IG=1,NG + CALL ALLUF(LL4F,IA11X(1,IG),MUX,IMAX) + ENDDO + IF(NDIM.GT.1) THEN + IA11Y(:NMAY,:NG)=A11Y(:NMAY,:NG) + DO IG=1,NG + CALL ALLUF(LL4F,IA11Y(1,IG),MUY,IMAY) + ENDDO + ENDIF + IF(NDIM.EQ.3) THEN + IA11Z(:NMAZ,:NG)=A11Z(:NMAZ,:NG) + DO IG=1,NG + CALL ALLUF(LL4F,IA11Z(1,IG),MUZ,IMAZ) + ENDDO + ENDIF +*---- +* POWER METHOD +*---- + NCTOT=ICL1+ICL2 + IF(ICL2.EQ.0) THEN + NCPTM=NCTOT+1 + ELSE + NCPTM=ICL1 + ENDIF + EVAL=1.0D0/FKEFF + VVV=0.0D0 + ISTART=1 + NNADI=NADI + TEST=0.0 + IF(IMPX.GE.2) WRITE (6,600) NADI + ITER=0 + DO + ITER=ITER+1 + IF(ITER > MAXOUT) CALL XABORT('NSSEIG: OUTER ITER. FAILURE.') +*---- +* EIGENVALUE EVALUATION +*---- + CALL NSSMPA(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL, + > VOL,MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,SCAT,A11X,A11Y,A11Z, + > EVECT,WORK(1,1,1)) + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,EVECT, + > WORK(1,1,2)) + AEBE=0.0D0 + BEBE=0.0D0 + DO IG=1,NG + DO I=1,LL4F + AEBE=AEBE+WORK(I,IG,1)*WORK(I,IG,2) + BEBE=BEBE+WORK(I,IG,2)**2 + ENDDO + ENDDO + EVAL=AEBE/BEBE + S(:LL4F,:NG)=REAL(EVAL)*WORK(:LL4F,:NG,2)-WORK(:LL4F,:NG,1) +*---- +* PERFORM THERMAL (UP-SCATTERING) ITERATIONS +*---- + WORK(:LL4F,:NG,:3)=0.0D0 + IGDEB=1 + TEXT3='NO ' + JTER=1 + ALLOCATE(F1(LL4F),GARM1(LL4F),GARM2(LL4F)) + DO + WORK(:LL4F,:NG,1)=WORK(:LL4F,:NG,2) + WORK(:LL4F,:NG,2)=WORK(:LL4F,:NG,3) + WORK(:LL4F,:NG,3)=0.0D0 + GRAD1(:LL4F,:NG)=0.0D0 + DO IG=IGDEB,NG + S2(:LL4F)=S(:LL4F,IG) + DO JG=1,NG + IF(JG.EQ.IG) CYCLE + DO IEL=1,NEL + IBM=MAT(IEL) + IF(IBM.LE.0) CYCLE + IND=IDL(IEL) + IF(IND.EQ.0) CYCLE + S2(IND)=S2(IND)+VOL(IEL)*SCAT(IBM,IG,JG)*GRAD1(IND,JG) + ENDDO + ENDDO +* + WORK(:LL4F,IG,3)=0.0 + DO IADI=1,NNADI + IF(IADI.EQ.1) THEN + F1(:LL4F)=S2(:LL4F) + ELSE +* scalar multiplication for a x-oriented matrix. + CALL ALLUM(LL4F,A11X(1,IG),WORK(1,IG,3),F1(1),MUX, + 1 IMAX,1) + IF(NDIM.GE.2) THEN +* scalar multiplication for a y-oriented matrix. + GARM1(IPY(:LL4F))=WORK(:LL4F,IG,3) + GARM2(IPY(:LL4F))=F1(:LL4F) + CALL ALLUM(LL4F,A11Y(1,IG),GARM1(1),GARM2(1),MUY, + 1 IMAY,2) + F1(:LL4F)=GARM2(IPY(:LL4F)) + ENDIF + IF(NDIM.EQ.3) THEN +* scalar multiplication for a z-oriented matrix. + GARM1(IPZ(:LL4F))=WORK(:LL4F,IG,3) + GARM2(IPZ(:LL4F))=F1(:LL4F) + CALL ALLUM(LL4F,A11Z(1,IG),GARM1(1),GARM2(1),MUZ, + 1 IMAZ,2) + F1(:LL4F)=GARM2(IPZ(:LL4F)) + ENDIF + F1(:LL4F)=S2(:LL4F)-F1(:LL4F) + ENDIF +* scalar solution for a x-oriented linear system. + CALL ALLUS(LL4F,MUX,IMAX,IA11X(1,IG),F1) + IF(NDIM.GE.2) THEN +* scalar solution for a y-oriented linear system. + DO I=1,LL4F + II=IPY(I) + GARM1(II)=F1(I)*A11Y(MUY(II),IG) + ENDDO + CALL ALLUS(LL4F,MUY,IMAY,IA11Y(1,IG),GARM1) + F1(:LL4F)=GARM1(IPY(:LL4F)) + ENDIF + IF(NDIM.EQ.3) THEN +* scalar solution for a z-oriented linear system. + DO I=1,LL4F + II=IPZ(I) + GARM1(II)=F1(I)*A11Z(MUZ(II),IG) + ENDDO + CALL ALLUS(LL4F,MUZ,IMAZ,IA11Z(1,IG),GARM1) + F1(:LL4F)=GARM1(IPZ(:LL4F)) + ENDIF + WORK(:LL4F,IG,3)=WORK(:LL4F,IG,3)+F1(:LL4F) + GRAD1(:LL4F,IG)=WORK(:LL4F,IG,3) + ENDDO + + ENDDO + IF(MAXTHR.EQ.0) EXIT + IF(MOD(JTER-1,NCTOT).GE.NCPTM) THEN + CALL NSS2AC(NG,LL4F,IGDEB,WORK,ZMU) + ELSE + ZMU=1.0D0 + ENDIF + IGDEBO=IGDEB + DO IG=IGDEBO,NG + GINN=0.0D0 + FINN=0.0D0 + DO I=1,LL4F + GINN=MAX(GINN,ABS(WORK(I,IG,2)-WORK(I,IG,3))) + FINN=MAX(FINN,ABS(WORK(I,IG,3))) + ENDDO + GINN=GINN/FINN + IF((GINN.LT.EPSTHR).AND.(IGDEB.EQ.IG)) IGDEB=IGDEB+1 + ENDDO + IF(GINN.LT.EPSTHR) TEXT3='YES' + IF(IMPX.GT.2) WRITE(6,610) JTER,GINN,EPSTHR,IGDEB,ZMU,TEXT3 + IF((GINN.LT.EPSTHR).OR.(JTER.EQ.MAXTHR)) EXIT + JTER=JTER+1 + ENDDO + DEALLOCATE(GARM2,GARM1,F1) +*---- +* DISPLACEMENT EVALUATION +*---- + F=0.0D0 + DELS=ABS(REAL((EVAL-VVV)/EVAL)) + VVV=EVAL +*---- +* EVALUATION OF THE TWO ACCELERATION PARAMETERS ALP AND BET +*---- + ALP=1.0D0 + BET=0.0D0 + N=0 + AEAE=0.0D0 + AEAG=0.0D0 + AEAH=0.0D0 + AGAG=0.0D0 + AGAH=0.0D0 + AHAH=0.0D0 + BEBG=0.0D0 + BEBH=0.0D0 + BGBG=0.0D0 + BGBH=0.0D0 + BHBH=0.0D0 + AEBG=0.0D0 + AEBH=0.0D0 + AGBE=0.0D0 + AGBG=0.0D0 + AGBH=0.0D0 + AHBE=0.0D0 + AHBG=0.0D0 + AHBH=0.0D0 + CALL NSSMPA(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL,VOL, + > MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,SCAT,A11X,A11Y,A11Z,EVECT, + > GAR1) + CALL NSSMPA(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL,VOL, + > MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,SCAT,A11X,A11Y,A11Z,GRAD1, + > GAR2) + IF(1+MOD(ITER-ISTART,ICL1+ICL2).GT.ICL1) THEN + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,EVECT,GAF1) + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,GRAD1,GAF2) + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,GRAD2,GAF3) + DO IG=1,NG + DO I=1,LL4F +* COMPUTE (A ,A ) + AEAE=AEAE+GAR1(I,IG)**2 + AEAG=AEAG+GAR1(I,IG)*GAR2(I,IG) + AEAH=AEAH+GAR1(I,IG)*GAR3(I,IG) + AGAG=AGAG+GAR2(I,IG)**2 + AGAH=AGAH+GAR2(I,IG)*GAR3(I,IG) + AHAH=AHAH+GAR3(I,IG)**2 +* COMPUTE (B ,B ) + BEBG=BEBG+GAF1(I,IG)*GAF2(I,IG) + BEBH=BEBH+GAF1(I,IG)*GAF3(I,IG) + BGBG=BGBG+GAF2(I,IG)**2 + BGBH=BGBH+GAF2(I,IG)*GAF3(I,IG) + BHBH=BHBH+GAF3(I,IG)**2 +* COMPUTE (A ,B ) + AEBG=AEBG+GAR1(I,IG)*GAF2(I,IG) + AEBH=AEBH+GAR1(I,IG)*GAF3(I,IG) + AGBE=AGBE+GAR2(I,IG)*GAF1(I,IG) + AGBG=AGBG+GAR2(I,IG)*GAF2(I,IG) + AGBH=AGBH+GAR2(I,IG)*GAF3(I,IG) + AHBE=AHBE+GAR3(I,IG)*GAF1(I,IG) + AHBG=AHBG+GAR3(I,IG)*GAF2(I,IG) + AHBH=AHBH+GAR3(I,IG)*GAF3(I,IG) + ENDDO + ENDDO +* + 210 N=N+1 + IF(N.GT.10) GO TO 215 +* COMPUTE X(ITER+1) + X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH + > +ALP*BET*BGBH) + DXDA=2.0D0*(BEBG+ALP*BGBG+BET*BGBH) + DXDB=2.0D0*(BEBH+ALP*BGBH+BET*BHBH) +* COMPUTE Y(ITER+1) + Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH + > +ALP*BET*AGAH) + DYDA=2.0D0*(AEAG+ALP*AGAG+BET*AGAH) + DYDB=2.0D0*(AEAH+ALP*AGAH+BET*AHAH) +* COMPUTE Z(ITER+1) + Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE) + > +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG) + DZDA=AEBG+AGBE+2.0D0*ALP*AGBG+BET*(AGBH+AHBG) + DZDB=AEBH+AHBE+ALP*(AGBH+AHBG)+2.0D0*BET*AHBH +* COMPUTE F(ITER+1) + F=X*Y-Z*Z + D2F(1,1)=2.0D0*(BGBG*Y+DXDA*DYDA+X*AGAG-DZDA**2-2.0D0*Z*AGBG) + D2F(1,2)=2.0D0*BGBH*Y+DXDA*DYDB+DXDB*DYDA+2.0D0*X*AGAH + > -2.0D0*DZDA*DZDB-2.0D0*Z*(AGBH+AHBG) + D2F(2,2)=2.0D0*(BHBH*Y+DXDB*DYDB+X*AHAH-DZDB**2-2.0D0*Z*AHBH) + D2F(2,1)=D2F(1,2) + D2F(1,3)=DXDA*Y+X*DYDA-2.0D0*Z*DZDA + D2F(2,3)=DXDB*Y+X*DYDB-2.0D0*Z*DZDB +* SOLUTION OF A LINEAR SYSTEM. + CALL ALSBD(2,1,D2F,IER,2) + IF(IER.NE.0) GO TO 215 + ALP=ALP-D2F(1,3) + BET=BET-D2F(2,3) + IF(ALP.GT.100.0D0) GO TO 215 + IF((ABS(D2F(1,3)).LE.1.0D-4).AND.(ABS(D2F(2,3)).LE.1.0D-4)) + > GO TO 220 + GO TO 210 +* +* alternative algorithm in case of Newton-Raphton failure + 215 IF(IMPX.GT.0) WRITE(6,'(/30H NSSEIG: FAILURE OF THE NEWTON, + > 55H-RAPHTON ALGORIHTHM FOR COMPUTING THE OVERRELAXATION PA, + > 9HRAMETERS.)') + IAMIN=999 + IBMIN=999 + FMIN=HUGE(FMIN) + DO IA=1,SIZE(ALP_TAB) + ALP=ALP_TAB(IA) + DO IB=1,SIZE(BET_TAB) + BET=BET_TAB(IB) +* COMPUTE X + X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH + > +ALP*BET*BGBH) +* COMPUTE Y + Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH + > +ALP*BET*AGAH) +* COMPUTE Z + Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE) + > +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG) +* COMPUTE F + F=X*Y-Z*Z + IF(F.LT.FMIN) THEN + IAMIN=IA + IBMIN=IB + FMIN=F + ENDIF + ENDDO + ENDDO + ALP=ALP_TAB(IAMIN) + BET=BET_TAB(IBMIN) + 220 BET=BET/ALP + IF((ALP.LT.1.0D0).AND.(ALP.GT.0.0D0)) THEN + ALP=1.0D0 + BET=0.0D0 + ELSE IF(ALP.LE.0.0D0) THEN + ISTART=ITER+1 + ALP=1.0D0 + BET=0.0D0 + ENDIF + DO IG=1,NG + DO I=1,LL4F + GRAD1(I,IG)=REAL(ALP)*(GRAD1(I,IG)+REAL(BET)*GRAD2(I,IG)) + GAR2(I,IG)=REAL(ALP)*(GAR2(I,IG)+REAL(BET)*GAR3(I,IG)) + ENDDO + ENDDO + ENDIF +* + LOGTES=(ITER.LT.ICL1).OR.(MOD(ITER-ISTART,ICL1+ICL2).EQ.ICL1-1) + DELT=0.0D0 + IF(LOGTES.AND.(DELS.LE.EPS1)) THEN + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,EVECT,GAF1) + CALL NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,GRAD1,GAF2) + DO IG=1,NG + DELN=0.0D0 + DELD=0.0D0 + DO I=1,LL4F + EVECT(I,IG)=EVECT(I,IG)+GRAD1(I,IG) + GAR1(I,IG)=GAR1(I,IG)+GAR2(I,IG) + GRAD2(I,IG)=GRAD1(I,IG) + GAR3(I,IG)=GAR2(I,IG) + DELN=MAX(DELN,ABS(GAF2(I,IG))) + DELD=MAX(DELD,ABS(GAF1(I,IG))) + ENDDO + IF(DELD.NE.0.0D0) DELT=MAX(DELT,DELN/DELD) + ENDDO + IF(IMPX.GE.2) WRITE (6,620) ITER,AEAE,AEAG,AEAH,AGAG,AGAH, + > AHAH,BEBE,ALP,BET,EVAL,F,DELS,DELT,N,BEBG,BEBH,BGBG,BGBH, + > BHBH,AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH + IF(DELT.LE.EPSOUT) EXIT + ELSE + DO IG=1,NG + DO I=1,LL4F + EVECT(I,IG)=EVECT(I,IG)+GRAD1(I,IG) + GAR1(I,IG)=GAR1(I,IG)+GAR2(I,IG) + GRAD2(I,IG)=GRAD1(I,IG) + GAR3(I,IG)=GAR2(I,IG) + ENDDO + ENDDO + IF(IMPX.GE.2) WRITE (6,620) ITER,AEAE,AEAG,AEAH,AGAG,AGAH, + > AHAH,BEBE,ALP,BET,EVAL,F,DELS,DELT,N,BEBG,BEBH,BGBG,BGBH, + > BHBH,AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH + ENDIF +* + IF(ITER.EQ.1) TEST=DELS + IF((ITER.GT.5).AND.(DELS.GT.TEST)) CALL XABORT('NSSEIG: CONVER' + > //'GENCE FAILURE.') + IF(ITER.GE.MAXOUT) THEN + WRITE (6,630) + EXIT + ENDIF + IF(MOD(ITER,36).EQ.0) THEN + ISTART=ITER+1 + NNADI=NNADI+1 + IF(IMPX.GE.1) WRITE (6,650) NNADI + ENDIF + ENDDO +*---- +* FLUX NORMALIZATION +*---- + FMAX=MAXVAL(EVECT(:LL4F,:NG)) + EVECT(:LL4F,:NG)=EVECT(:LL4F,:NG)/FMAX +*---- +* SOLUTION EDITION +*---- + FKEFF=REAL(1.0D0/EVAL) + IF(IMPX.GE.1) WRITE (6,640) ITER,FKEFF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(IA11Z,IA11Y,IA11X) + DEALLOCATE(WORK,GAF3,GAF2,GAF1,GAR3,GAR2,GAR1,GRAD2,GRAD1,S,S2) + RETURN +*---- +* FORMATS +*---- + 600 FORMAT(1H1/50H NSSEIG: ITERATIVE PROCEDURE BASED ON PRECONDITION, + > 17HED POWER METHOD (,I2,37H ADI ITERATIONS PER OUTER ITERATION)./ + > 9X,16HDIRECT EQUATION.) + 610 FORMAT (10X,3HIN(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2, + > 7H IGDEB=, I13,6H ACCE=,0P,F12.5,12H CONVERGED=,A3) + 620 FORMAT(1X,I3,1P,7E9.1,0P,2F8.3,E14.6,3E10.2,I4/(4X,1P,7E9.1)) + 630 FORMAT(/53H NSSEIG: ***WARNING*** THE MAXIMUM NUMBER OF OUTER IT, + > 20HERATIONS IS REACHED.) + 640 FORMAT(/23H NSSEIG: CONVERGENCE IN,I4,12H ITERATIONS.// + > 42H NSSEIG: EFFECTIVE MULTIPLICATION FACTOR =,1P,E17.10/) + 650 FORMAT(/53H NSSEIG: INCREASING THE NUMBER OF INNER ITERATIONS TO, + 1 I3,36H ADI ITERATIONS PER OUTER ITERATION./) + ! + CONTAINS + SUBROUTINE NSSMPA(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL, + > VOL,MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,SCAT,A11X,A11Y,A11Z, + > EVECT,S2) + ! + ! A*EVECT MULTIPLICATION + ! + INTEGER, INTENT(IN) :: NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG, + > MAT(NEL),IDL(NEL),MUX(LL4F),MUY(LL4F),MUZ(LL4F),IMAX(LL4F), + > IMAY(LL4F),IMAZ(LL4F),IPY(LL4F),IPZ(LL4F) + REAL, INTENT(IN) :: VOL(NEL),SCAT(NMIX,NG,NG) + REAL, INTENT(IN) :: EVECT(LL4F,NG),A11X(NMAX,NG),A11Y(NMAY,NG), + > A11Z(NMAZ,NG) + REAL, INTENT(OUT) :: S2(LL4F,NG) + REAL, ALLOCATABLE, DIMENSION(:) :: GAR1,GAR2 + ! + ALLOCATE(GAR1(LL4F),GAR2(LL4F)) + DO IG=1,NG +* scalar multiplication for a x-oriented matrix. + CALL ALLUM(LL4F,A11X(1,IG),EVECT(1,IG),S2(1,IG),MUX,IMAX,1) + IF(NDIM.GE.2) THEN +* scalar multiplication for a y-oriented matrix. + GAR1(IPY(:LL4F))=EVECT(:LL4F,IG) + GAR2(IPY(:LL4F))=S2(:LL4F,IG) + CALL ALLUM(LL4F,A11Y(1,IG),GAR1(1),GAR2(1),MUY,IMAY,2) + S2(:LL4F,IG)=GAR2(IPY(:LL4F)) + ENDIF + IF(NDIM.EQ.3) THEN +* scalar multiplication for a z-oriented matrix. + GAR1(IPZ(:LL4F))=EVECT(:LL4F,IG) + GAR2(IPZ(:LL4F))=S2(:LL4F,IG) + CALL ALLUM(LL4F,A11Z(1,IG),GAR1(1),GAR2(1),MUZ,IMAZ,2) + S2(:LL4F,IG)=GAR2(IPZ(:LL4F)) + ENDIF + DO JG=1,NG + IF(JG.EQ.IG) CYCLE + DO IEL=1,NEL + IBM=MAT(IEL) + IF(IBM.LE.0) CYCLE + IND=IDL(IEL) + IF(IND.EQ.0) CYCLE + S2(IND,IG)=S2(IND,IG)-VOL(IEL)*SCAT(IBM,IG,JG)* + > EVECT(IND,JG) + ENDDO + ENDDO + ENDDO + DEALLOCATE(GAR2,GAR1) + END SUBROUTINE NSSMPA + ! + SUBROUTINE NSSMPB(LL4F,NEL,NMIX,NG,MAT,IDL,VOL,CHI,SIGF,EVECT, + > S2) + ! + ! B*EVECT MULTIPLICATION + ! + INTEGER, INTENT(IN) :: LL4F,NEL,NMIX,NG,MAT(NEL),IDL(NEL) + REAL, INTENT(IN) :: VOL(NEL),CHI(NMIX,NG),SIGF(NMIX,NG) + REAL, INTENT(IN) :: EVECT(LL4F,NG) + REAL, INTENT(OUT) :: S2(LL4F,NG) + ! + S2(:LL4F,:NG)=0.0D0 + DO IG=1,NG + DO JG=1,NG ! IG <-- JG + DO IEL=1,NEL + IBM=MAT(IEL) + IF(IBM.LE.0) CYCLE + IND=IDL(IEL) + IF(IND.EQ.0) CYCLE + S2(IND,IG)=S2(IND,IG)+VOL(IEL)*CHI(IBM,IG)*SIGF(IBM,JG)* + > EVECT(IND,JG) + ENDDO + ENDDO + ENDDO + END SUBROUTINE NSSMPB + END SUBROUTINE NSSEIG |