diff options
Diffstat (limited to 'Dragon/src/PIJRHL.f')
| -rw-r--r-- | Dragon/src/PIJRHL.f | 194 |
1 files changed, 194 insertions, 0 deletions
diff --git a/Dragon/src/PIJRHL.f b/Dragon/src/PIJRHL.f new file mode 100644 index 0000000..9e2040b --- /dev/null +++ b/Dragon/src/PIJRHL.f @@ -0,0 +1,194 @@ +*DECK PIJRHL + SUBROUTINE PIJRHL(IPRT,NREG,NSOUT,SIGTAL,PROB) +* +*----------------------------------------------------------------------- +* +*Purpose: +* HELIOS type normalization of collision probs (CP). +* +*Copyright: +* Copyright (C) 1994 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, E. Varin +* +*Parameters: input +* IPRT print level. +* NREG number of zones for geometry. +* NSOUT number of surfaces for geometry. +* SIGTAL albedo-sigt vector. +* +*Parameters: input/output +* PROB CP matrix for all types. +* +*References: +* R. Roy and G. Marleau, +* Normalization Techniques for CP Matrices, +* CONF/PHYSOR-90, Marseille/France, V 2, P IX-40 (1990). +* \\\\ +* E.A. Villarino, R.J.J. Stamm'ler, A.A. Ferri and J.J. Casal +* HELIOS: Angularly Dependent Collision Probabilities. +* Nucl.Sci.Eng. 112,16-31, 1992. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IPRT,NREG,NSOUT + REAL SIGTAL(-NSOUT:NREG) + DOUBLE PRECISION PROB(*) +*---- +* LOCAL VARIABLES +*---- + INTEGER IUNOUT,NITMAX,NIT,IPRINT,IR,JR,IP,IPRB,IND,I,J,CPTLB, + > CPTAC,CTOT,NSURC,NSURM,NVOLC,NVOLM + LOGICAL NOTCON + DOUBLE PRECISION NOM,DENOM,DMU,WFSPAD,WFSP,EPSCON,R1,R2,TOTCON, + > TMPCON + CHARACTER HSMG*131 + PARAMETER (IUNOUT=6, IPRINT=10, EPSCON=1.0E-6, NITMAX=20) +*---- +* ALLOCATABLE ARRAYS +*---- + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: CHI + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: WEIG +* +*----- INTRINSIC FUNCTION FOR POSITION IN CONDENSE PIJ MATRIX +* + IND(I,J)=(MAX(I+NSOUT+1,J+NSOUT+1)* + > (MAX(I+NSOUT+1,J+NSOUT+1)-1))/2 + > +MIN(I+NSOUT+1,J+NSOUT+1) +*---- +* SCRATCH STORAGE ALLOCATION +* WEIG : ADDITIVE WEIGHT +*---- + ALLOCATE(WEIG(-NSOUT:NREG,3),CHI(-NSOUT:NREG)) +* + NOTCON= .FALSE. + CPTLB = 3 + CPTAC = 3 + CTOT = CPTAC+CPTLB +* +* INITIALISATION OF WEIGHTS + DO 60 IR=-NSOUT, NREG + WEIG(IR,1)=0.0D0 + WEIG(IR,2)=0.5D0 + WEIG(IR,3)=0.5D0 + 60 CONTINUE + DO 50 IR=-NSOUT, NREG + CHI(IR)= 1.0D0 + IF( IR.GE.0.AND.SIGTAL(IR).EQ.0.0D0 )THEN + CHI(IR)= 0.0D0 + ENDIF + 50 CONTINUE +* +* MAIN ITERATION LOOP + IF(IPRT.GT.2) WRITE(IUNOUT,'(A24)') + > 'ITER. MU ERROR ' + DO 110 NIT=1,NITMAX +* + DO 220 IR= -NSOUT, NREG + WFSPAD = PROB(IND(IR,0)) + > + CHI(IR)*PROB(IND(IR,IR))*WEIG(IR,3) + WFSP = CHI(IR)*PROB(IND(IR,IR)) + DO 200 JR=-NSOUT, NREG + WFSPAD = WFSPAD - CHI(JR)*WEIG(JR,3)*PROB(IND(IR,JR)) + WFSP = WFSP + CHI(JR)*PROB(IND(IR,JR)) + 200 CONTINUE + IF(WFSP.NE.0.0) WEIG(IR,3) = WFSPAD / WFSP + 220 CONTINUE +* +* ACCELERATION TECHNIQUE + IF( MOD(NIT-1,CTOT).GE.CPTAC )THEN + NOM = 0.0D0 + DENOM = 0.0D0 + DO 10 IR=-NSOUT, NREG + R1= WEIG(IR,2) - WEIG(IR,1) + R2= WEIG(IR,3) - WEIG(IR,2) + NOM = NOM + R1*(R2-R1) + DENOM = DENOM + (R2-R1)*(R2-R1) + 10 CONTINUE + IF(DENOM.EQ.0.0D0) THEN + DMU = 1.0D0 + ELSE + DMU = - NOM / DENOM + ENDIF + IF( DMU.GT.50.0D0 .OR. DMU.LT.0.0D0 ) THEN + WRITE(HSMG,'(37HPIJRHL: PROBLEM OF ACCELERATION (DMU=,1P, + > E11.4,2H).)') DMU + CALL XABORT(HSMG) + ENDIF + DO 20 IR=-NSOUT, NREG + WEIG(IR,3) = WEIG(IR,2) + DMU * + > (WEIG(IR,3) - WEIG(IR,2)) + WEIG(IR,2) = WEIG(IR,1) + DMU * + > (WEIG(IR,2) - WEIG(IR,1)) + 20 CONTINUE + ELSE + DMU = 1.0D0 + ENDIF +* +* CALCULATIONS OF SQUARE DISTANCE BETWEEN 2 ITERATIONS +* AND UPDATING THE SOLUTION + TOTCON = 0.0D0 + DO 100 IR=-NSOUT, NREG + TMPCON=ABS(WEIG(IR,3)-WEIG(IR,2))/WEIG(IR,3) + TOTCON=MAX(TMPCON,TOTCON) + WEIG(IR,1)= WEIG(IR,2) + WEIG(IR,2)= WEIG(IR,3) + 100 CONTINUE + IF( IPRT.GT.2 ) WRITE(IUNOUT,'(I3,F9.5,E15.7)') NIT,DMU,TOTCON +* +* CONVERGENCE TEST + IF( TOTCON.LT.EPSCON )GO TO 120 +* + 110 CONTINUE + NOTCON=.TRUE. + WRITE(IUNOUT,'(35H PIJRHL: WEIGHTS NOT CONVERGED )') + 120 CONTINUE +* +* RENORMALIZE "PIJ" SYMMETRIC MATRIX + IPRB = 0 + DO 240 IR = -NSOUT, NREG + DO 230 JR= -NSOUT, IR + IPRB= IPRB+1 + IF( IR.NE.0.AND.JR.NE.0 )THEN + PROB(IPRB)=PROB(IPRB)*(WEIG(IR,1)+WEIG(JR,1)) + ENDIF + 230 CONTINUE + 240 CONTINUE +* +* PRINT WEIGHT FACTORS IF THERE IS A PROBLEM... + IF( NOTCON .OR. IPRT.GE.IPRINT )THEN + WRITE(IUNOUT,'(30H0 SURFACE WEIGHTS FACTORS /)') + NSURC = -1 + DO 300 IP = 1, (9 +NSOUT) / 10 + NSURM= MAX( -NSOUT, NSURC-9 ) + WRITE(IUNOUT,'(10X,10( A5, I6)/)') + > (' SUR ',-IR,IR= NSURC, NSURM, -1) + WRITE(IUNOUT,'(10H WEIGHT ,10F11.5)') + > (WEIG(IR,1),IR=NSURC,NSURM,-1) + NSURC = NSURC - 10 + 300 CONTINUE + WRITE(IUNOUT,'(30H0 VOLUME WEIGHTS FACTORS /)') + NVOLC = 1 + DO 310 IP = 1, (9 + NREG) / 10 + NVOLM= MIN( NREG, NVOLC+9 ) + WRITE(IUNOUT,'(10X,10( A5 , I6)/)') + > (' VOL ',IR,IR=NVOLC,NVOLM, 1) + WRITE(IUNOUT,'(10H WEIGHT ,10F11.5)') + > (WEIG(IR,1),IR=NVOLC,NVOLM, 1) + NVOLC = NVOLC + 10 + 310 CONTINUE + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(CHI,WEIG) + RETURN + END |