diff options
Diffstat (limited to 'Dragon/src/PNFLV.f')
| -rw-r--r-- | Dragon/src/PNFLV.f | 207 |
1 files changed, 207 insertions, 0 deletions
diff --git a/Dragon/src/PNFLV.f b/Dragon/src/PNFLV.f new file mode 100644 index 0000000..d716737 --- /dev/null +++ b/Dragon/src/PNFLV.f @@ -0,0 +1,207 @@ +*DECK PNFLV + SUBROUTINE PNFLV(KPSYS,INCONV,NGIND,IPTRK,IMPX,MAXIT,NGEFF,NREG, + 1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solve N-group transport equation for fluxes using the spherical +* harmonics (PN) method in BIVAC. +* +*Copyright: +* Copyright (C) 2004 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. +* INCONV energy group convergence flag (set to .false. if converged). +* 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. +* 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. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) KPSYS(NGEFF),IPTRK + INTEGER MAXIT,NGEFF,NGIND(NGEFF),IMPX,NREG,NBMIX,NUN, + 1 MAT(NREG),KEYFLX(NREG) + LOGICAL INCONV(NGEFF) + REAL VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (IUNOUT=6,NSTATE=40,EPSINR=1.0E-5,ICL1=3, + 1 ICL2=3) + INTEGER IPAR(NSTATE) + DOUBLE PRECISION F1,F2,R1,R2,DMU + CHARACTER NAMP*12 +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: MU,KN,IPERT + REAL, ALLOCATABLE, DIMENSION(:) :: QFR,SGDI,RR,VV,SYS,OLD1,OLD2, + 1 XX,YY +*---- +* RECOVER PN SPECIFIC PARAMETERS. +*---- + CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR) + IF(NREG.NE.IPAR(1)) CALL XABORT('PNFLV: INVALID VALUE OF NREG.') + IF(NUN.NE.IPAR(2)) CALL XABORT('PNFLV: INVALID VALUE OF NUN.') + ITYPE=IPAR(6) + IELEM=IPAR(8) + ICOL=IPAR(9) + ISPLH=IPAR(10) + L4=IPAR(11) + LX=IPAR(12) + NLF=IPAR(14) + ISPN=IPAR(15) + ISCAT=IPAR(16) + NVD=IPAR(17) +*---- +* RECOVER TRACKING INFORMATION. +*---- + ALLOCATE(MU(L4)) + CALL LCMGET(IPTRK,'MU',MU) + CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) + CALL LCMLEN(IPTRK,'QFR',MAXQF,ITYLCM) + ALLOCATE(KN(MAXKN),QFR(MAXQF)) + CALL LCMGET(IPTRK,'KN',KN) + CALL LCMGET(IPTRK,'QFR',QFR) + IIMAX=MU(L4)*(NLF/2) + ALLOCATE(SYS(IIMAX)) +*---- +* RECOVER THE FINITE ELEMENT UNIT STIFFNESS MATRIX. +*---- + CALL LCMSIX(IPTRK,'BIVCOL',1) + CALL LCMLEN(IPTRK,'T',LC,ITYLCM) + ALLOCATE(RR(LC*LC),VV(LC*(LC-1))) + CALL LCMGET(IPTRK,'R',RR) + CALL LCMGET(IPTRK,'V',VV) + CALL LCMSIX(IPTRK,' ',2) +*---- +* MAIN LOOP OVER ENERGY GROUPS. +*---- + ALLOCATE(OLD1(NUN),OLD2(NUN)) + DO 140 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 140 + IF(IMPX.GT.1) WRITE(IUNOUT,'(/24H PNFLV: PROCESSING GROUP,I5, + 1 6H WITH ,A,1H.)') NGIND(II),'BIVAC/PN' +*---- +* RECOVER CROSS SECTIONS. +*---- + CALL LCMGET(KPSYS(II),'STATE-VECTOR',IPAR) + NAN=IPAR(8) + ALLOCATE(SGDI(NBMIX*NAN)) + DO 10 IL=0,NAN-1 + WRITE(NAMP,'(4HSCAI,I2.2,6H001001)') IL + CALL LCMGET(KPSYS(II),NAMP,SGDI(IL*NBMIX+1)) + 10 CONTINUE +*---- +* INNER ITERATION LOOP FOR ONE-GROUP TRANSPORT EQUATION. +*---- + CALL LCMGET(KPSYS(II),'IA001001',SYS) + OLD2(:NUN)=0.0 + TEST=0.0 + ITER=0 + 30 ITER=ITER+1 + IF(ITER.GT.MAXIT) THEN + WRITE(IUNOUT,'(43H PNFLV: MAXIMUM NUMBER OF ONE-SPEED ITERATI, + 1 11HON REACHED.)') + GO TO 100 + ENDIF + DO 40 I=1,NUN + OLD1(I)=OLD2(I) + OLD2(I)=FUNKNO(I,II) + 40 CONTINUE + IF((ITYPE.EQ.2).OR.((ITYPE.EQ.5).AND.(ISPN.EQ.1))) THEN + ALLOCATE(XX(NREG),YY(NREG)) + CALL LCMGET(IPTRK,'XX',XX) + CALL LCMGET(IPTRK,'YY',YY) + CALL PNFL2E(NREG,IELEM,ICOL,XX,YY,MAT,VOL,NBMIX,NLF,NVD, + 1 NAN,SGDI,L4,KN,QFR,MU,IIMAX,LC,RR,VV,SYS,SUNKNO(1,II), + 2 FUNKNO(1,II),1) + DEALLOCATE(YY,XX) + ELSE IF((ITYPE.EQ.2).OR.((ITYPE.EQ.8).AND.(ISPN.EQ.1))) THEN + NBLOS=LX/3 + ALLOCATE(IPERT(NBLOS)) + CALL LCMGET(IPTRK,'SIDE',SIDE) + CALL LCMGET(IPTRK,'IPERT',IPERT) + CALL PNFH2E(IELEM,ICOL,NBLOS,SIDE,NLF,NVD,L4,IPERT,KN,QFR,MU, + 1 IIMAX,LC,VV,SYS,SUNKNO(1,II),FUNKNO(1,II),1) + DEALLOCATE(IPERT) + ELSE + CALL XABORT('PNFLV: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') + ENDIF + IF(NLF.EQ.2) GO TO 100 +*---- +* VARIATIONAL ACCELERATION. +*---- + DMU=1.0D0 + IF(MOD(ITER-1,ICL1+ICL2).GE.ICL1) THEN + F1=0.0 + F2=0.0 + DO 70 I=1,NUN + R1=OLD2(I)-OLD1(I) + R2=FUNKNO(I,II)-OLD2(I) + F1=F1+R1*(R2-R1) + F2=F2+(R2-R1)*(R2-R1) + 70 CONTINUE + DMU=-F1/F2 + IF(DMU.GT.0.0) THEN + DO 80 I=1,NUN + FUNKNO(I,II)=OLD2(I)+REAL(DMU)*(FUNKNO(I,II)-OLD2(I)) + OLD2(I)=OLD1(I)+REAL(DMU)*(OLD2(I)-OLD1(I)) + 80 CONTINUE + ENDIF + ENDIF +*--- CALCULATE ERROR AND TEST FOR CONVERGENCE. + AAA=0.0 + BBB=0.0 + DO 90 I=1,NREG + IF(KEYFLX(I).EQ.0) GO TO 90 + AAA=MAX(AAA,ABS(FUNKNO(KEYFLX(I),II)-OLD2(KEYFLX(I)))) + BBB=MAX(BBB,ABS(FUNKNO(KEYFLX(I),II))) + 90 CONTINUE + IF(IMPX.GT.2) WRITE(IUNOUT,300) ITER,AAA,BBB,DMU + IF(AAA.LE.EPSINR*BBB) GO TO 100 + IF(ITER.EQ.1) TEST=AAA + IF((ITER.GE.10).AND.(AAA.GT.TEST)) THEN + WRITE(IUNOUT,'(43H PNFLV: UNABLE TO CONVERGE ONE-SPEED ITERAT, + 1 5HIONS.)') + GO TO 100 + ENDIF + GO TO 30 + 100 DEALLOCATE(SGDI) +*---- +* END OF LOOP OVER ENERGY GROUPS. +*---- + 140 CONTINUE + DEALLOCATE(OLD2,OLD1,SYS,VV,RR,QFR,KN,MU) + IF((IMPX.GT.0).AND.(NLF.GT.2)) WRITE(IUNOUT,'(15H PNFLV: NUMBER , + 1 24HOF ONE-SPEED ITERATIONS=,I5,1H.)') ITER + RETURN +* + 300 FORMAT(27H PNFLV: ONE-SPEED ITERATION,I3,8H ERROR=,1P,E11.4, + 1 5H OVER,E11.4,22H ACCELERATION FACTOR=,0P,F7.3) + END |