diff options
Diffstat (limited to 'Dragon/src/MCGDS2.f')
| -rw-r--r-- | Dragon/src/MCGDS2.f | 232 |
1 files changed, 232 insertions, 0 deletions
diff --git a/Dragon/src/MCGDS2.f b/Dragon/src/MCGDS2.f new file mode 100644 index 0000000..16ac09f --- /dev/null +++ b/Dragon/src/MCGDS2.f @@ -0,0 +1,232 @@ +*DECK MCGDS2 + SUBROUTINE MCGDS2(SUBDSC,LC,M,N,H,NOM,NZON,TR,SC,W,NFI,DIAGF, + 1 DIAGQ,CA,CQ,PREV,NEXT,DINV2,A2,B2) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Calculation of contribution in second-order ACA coefficients on one +* track. +* +*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): I. Suslov and R. Le Tellier +* +*Parameters: input +* SUBDSC ACA coefficients calculation subroutine. +* LC dimension of vector MCU. +* M number of material mixtures. +* N number of elements for this track. +* H tracking widths. +* NOM integer tracking elements. +* NZON index-number of the mixture type assigned to each volume. +* TR macroscopic total cross section. +* SC macroscopic P0 scattering cross section. +* W weight associated with this track. +* NFI total number of volumes and surfaces for which specific values +* of the neutron flux and reactions rates are required. +*Parameters: input/output +* CA undefined. +* CQ undefined. +* DIAGQ undefined. +* DIAGF undefined. +* +*Parameters: scratch +* PREV undefined. +* NEXT undefined. +* DINV2 undefined. +* A2 undefined. +* B2 undefined. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +* +*--- +* SUBROUTINE ARGUMENTS +*--- + INTEGER LC,M,N,NFI,NZON(NFI),NOM(N),PREV(N),NEXT(N) + DOUBLE PRECISION W,H(N),CA(LC),DIAGF(NFI),DINV2(N),A2(N),B2(N) + REAL TR(0:M),SC(0:M),DIAGQ(NFI),CQ(LC) + EXTERNAL SUBDSC +*--- +* LOCAL VARIABLES +*--- + INTEGER IBCV + DOUBLE PRECISION DMINV,DMAX + PARAMETER(DMINV=2.D-2,IBCV=-7) + DOUBLE PRECISION WW,AAW,CAW,AQW,CQW,CAWP,CQWP, + 1 B,A,DINV,DH,DHP,BN,AN,DINVN,BP,AP,DINVP + INTEGER I,I1,NOMI,NZI,NOMIN,NZIN,NOMIP,NZIP,ICN,ICP +* + DMAX=1.D0/DMINV + WW=W +*--- +* CALCULATE COEFFICIENTS OF THIS TRACK +*--- +* MCGDS2A: Tabulated Exponentials +* MCGDS2E: Exact Exponentials + CALL SUBDSC(N,M,NFI,NOM,NZON,H,TR,SC,DINV2,B2,A2) +*---- +* CONSTRUCTION OF ACA MATRICES +*--- +* ------------------- +* Left outer boundary +* ------------------- + I1=2 + ICN=NEXT(1) + ICP=PREV(1) + NZIP=IBCV + NOMI=NOM(1) + NZI=NZON(NOMI) + NOMIN=NOM(I1) + NZIN=NZON(NOMIN) + IF (NZI.NE.IBCV) THEN +* Other than void boundary condition (treated as white reflection) + DIAGF(NOMI)=DIAGF(NOMI)+WW + DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*(DINV2(I1)-1.D0)) + IF(ICN.GT.0) THEN + CA(ICN)=CA(ICN)-WW*A2(I1) + CQ(ICN)=CQ(ICN)+REAL(W*B2(I1)) + ENDIF + ENDIF +* ------------ +* Volume Cells +* ------------ + DHP=0.0 + AP=0.0 + BP=0.0 + AN=0.0 + BN=0.0 + DO I=2,N-1 + ICN=NEXT(I) + ICP=PREV(I) + NOMIP=NOMI + NZIP=NZI + NOMI=NOMIN + NZI=NZIN + NOMIN=NOM(I+1) + NZIN=NZON(NOMIN) + DINV=DINV2(I) + B=B2(I) + A=A2(I) + DH=0.0 + IF (NZIN.EQ.IBCV) THEN +* next cell is a fixed boundary condition + DH=1.D0/(1.D0+DINV) + ELSEIF (NZIN.GE.0) THEN +* next cell is a volume + I1=I+1 + DINVN=DINV2(I1) + BN=B2(I1) + AN=A2(I1) + IF (ABS(DINV+DINVN).LT.DMINV) THEN + DH=DMAX + ELSE + DH=1.D0/(DINV+DINVN) + ENDIF + ENDIF + IF (NZIP.EQ.IBCV) THEN +* previous cell is a fixed boundary condition + DHP=1.D0/(1.D0+DINV) + ELSEIF (NZIP.GE.0) THEN +* previous cell is a volume + I1=I-1 + DINVP=DINV2(I1) + BP=B2(I1) + AP=A2(I1) + IF (ABS(DINV+DINVP).LT.DMINV) THEN + DHP=DMAX + ELSE + DHP=1.D0/(DINV+DINVP) + ENDIF + ENDIF +* assembling coefficients + IF ((NZIN.LT.0).AND.(NZIN.NE.IBCV)) THEN +* next cell is a surface with reflective boundary condition + AAW=0.D0 + AQW=0.D0 + CAW=0.D0 + CQW=1.D0 + ELSE +* next cell is a volume or a fixed boundary condition + AAW=DH*A + AQW=DH*B + IF (NZIN.GE.0) THEN +* next cell is a volume + CAW=DH*AN + CQW=DH*BN + ELSE +* next cell is a fixed boundary condition + CAW=0.D0 + CQW=0.D0 + ENDIF + ENDIF + IF ((NZIP.LT.0).AND.(NZIP.NE.IBCV)) THEN +* previous cell is a surface with reflective boundary condition + CAWP=0.D0 + CQWP=1.D0 + ELSE +* previous cell is a volume or a fixed boundary condition + AAW=AAW+DHP*A + AQW=AQW+DHP*B + IF (NZIP.GE.0) THEN +* previous cell is a volume + CAWP=DHP*AP + CQWP=DHP*BP + ELSE +* previous cell is a fixed boundary condition + CAWP=0.D0 + CQWP=0.D0 + ENDIF + ENDIF +* assembling matrices + DIAGF(NOMI)=DIAGF(NOMI)+AAW*WW + DIAGQ(NOMI)=DIAGQ(NOMI)-REAL(W*AQW) + IF(ICN.GT.0) THEN +* next cell is a volume different from this one + CA(ICN)=CA(ICN)-CAW*WW + CQ(ICN)=CQ(ICN)+REAL(W*CQW) + ELSE +* next cell is a voided boundary or a volume identical to this one + DIAGF(NOMI)=DIAGF(NOMI)-CAW*WW + DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*CQW) + ENDIF + IF(ICP.GT.0) THEN +* previous cell is a volume different from this one + CA(ICP)=CA(ICP)-CAWP*WW + CQ(ICP)=CQ(ICP)+REAL(W*CQWP) + ELSE +* previous cell is a voided boundary or a volume identical to this one + DIAGF(NOMI)=DIAGF(NOMI)-CAWP*WW + DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*CQWP) + ENDIF + ENDDO +* -------------------- +* Right outer boundary +* -------------------- + ICN=NEXT(N) + ICP=PREV(N) + NOMIP=NOMI + NZIP=NZI + NOMI=NOMIN + NZI=NZIN + NZIN=IBCV + IF (NZI.NE.IBCV) THEN +* Other than void boundary condition (treated as white reflection) + I1=N-1 + DIAGF(NOMI)=DIAGF(NOMI)+WW + DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*(DINV2(I1)-1.D0)) + IF(ICP.GT.0) THEN + CA(ICP)=CA(ICP)-WW*A2(I1) + CQ(ICP)=CQ(ICP)+REAL(W*B2(I1)) + ENDIF + ENDIF +* + RETURN + END |