diff options
Diffstat (limited to 'Dragon/src/SNFLUX.f')
| -rw-r--r-- | Dragon/src/SNFLUX.f | 1129 |
1 files changed, 1129 insertions, 0 deletions
diff --git a/Dragon/src/SNFLUX.f b/Dragon/src/SNFLUX.f new file mode 100644 index 0000000..09a832d --- /dev/null +++ b/Dragon/src/SNFLUX.f @@ -0,0 +1,1129 @@ +*DECK SNFLUX + SUBROUTINE SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG, + 1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,ITER,NBS,KPSOU1,KPSOU2, + 2 FLUXC) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solve a single non-accelerated scattering iteration of the N-group +* transport equation for fluxes using the discrete ordinates (SN) +* method. +* +*Copyright: +* Copyright (C) 2007 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. +* ITER number of previous calls to SNFLUX. +* NBS +* KPSOU1 +* KPSOU2 +* +*Parameters: input/output +* FUNKNO unknown vector. +* FLUXC flux at the cutoff energy. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) KPSYS(NGEFF),IPTRK,KPSOU1(NGEFF),KPSOU2(NGEFF) + INTEGER NGEFF,NGIND(NGEFF),IMPX,NREG,NBMIX,NUN,MAT(NREG), + 1 KEYFLX(NREG),ITER,NBS(NGEFF) + LOGICAL INCONV(NGEFF) + REAL VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF) + REAL,OPTIONAL :: FLUXC(NREG) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (IUNOUT=6,NSTATE=40) + INTEGER IPAR(NSTATE),NCODE(6),BSINFO(2),EELEM,ESCHM,P,Q, + > LOZSWP(3,6) + REAL ZCODE(6),SIDE,LAMBDA0 + LOGICAL LFIXUP,LDSA,LSHOOT,LADPT(4) +*---- +* ALLOCATABLE ARRAYS +*--- + INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYSPN + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: COORDMAP + INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: ISBSM + REAL, ALLOCATABLE, DIMENSION(:,:) :: QEXT,OLD,SGAR,BS,EMOMTR, + 1 SIGMATR,MTR + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: SGAS,ESTOPW +* + TYPE(C_PTR) U_PTR,W_PTR,PL_PTR,JOP_PTR,UPQ_PTR,WPQ_PTR,ALPHA_PTR, + 1 PLZ_PTR,SURF_PTR,XXX_PTR,DU_PTR,DE_PTR,MRM_PTR,MRMX_PTR,MRMY_PTR, + 2 MRMZ_PTR,DB_PTR,DA_PTR,DAL_PTR,DZ_PTR,DC_PTR,WX_PTR,WE_PTR, + 3 CST_PTR,MN_PTR,DN_PTR,IL_PTR + INTEGER, POINTER, DIMENSION(:) :: JOP,MRM,MRMX,MRMY,MRMZ,ISLG,IL + REAL, POINTER, DIMENSION(:) :: U,W,PL,UPQ,WPQ,ALPHA,PLZ,SURF,DU, + 1 DE,XXX,DB,DA,DAL,DZ,DC,DELTAE,WX,WE,CST,MN,DN +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(QEXT(NUN,NGEFF),KEYSPN(NREG),DELTAE(NGEFF),ISLG(NGEFF)) +*---- +* RECOVER SN SPECIFIC PARAMETERS. +*---- + CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR) + ITYPE=IPAR(6) + NSCT=IPAR(7) + IELEM=IPAR(8) + ISCHM=IPAR(10) + L4=IPAR(11) + LX=IPAR(12) + LY=IPAR(13) + LZ=IPAR(14) + NLF=IPAR(15) + ISCAT=IPAR(16) + LFIXUP=IPAR(18).EQ.1 + LDSA=IPAR(19).EQ.1 + NSDSA=IPAR(21) + ISPLH=IPAR(26) + NKBA=IPAR(28) + IGAV=IPAR(29) + LSHOOT=.TRUE. + IF(IPAR(30).EQ.0) LSHOOT=.FALSE. + IBFP=IPAR(31) + EELEM=IPAR(35) + ESCHM=IPAR(36) + LADPT(1:4)=.FALSE. + IF(ESCHM.EQ.3) LADPT(1)=.TRUE. + IF(ISCHM.EQ.3) LADPT(2:4)=.TRUE. +*---- +* TEST FOR DSA and SAVE OLD FLUX. +*---- + LDSA=.FALSE. + IF(MOD(ITER,(NSDSA+1))==0) LDSA=.TRUE. + + IF(LDSA)THEN + ALLOCATE(OLD(NUN,NGEFF)) + OLD(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF) + ENDIF +*---- +* RECOVER TOTAL AND WITHIN-GROUP SCATTERING MULTIGROUP CROSS SECTIONS. +*---- + ALLOCATE(SGAR(0:NBMIX,NGEFF),SGAS(0:NBMIX,ISCAT,NGEFF)) + ALLOCATE(ESTOPW(0:NBMIX,2,NGEFF),EMOMTR(0:NBMIX,NGEFF)) + NANI=1 + DO 10 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 10 + CALL LCMLEN(KPSYS(II),'DRAGON-TXSC',ILONG,ITYLCM) + IF(ILONG.NE.NBMIX+1) CALL XABORT('SNFLUX: INVALID TXSC LENGTH.') + CALL LCMLEN(KPSYS(II),'DRAGON-S0XSC',ILONG,ITYLCM) + NANI=MAX(NANI,ILONG/(NBMIX+1)) + IF(NANI.GT.ISCAT) CALL XABORT('SNFLUX: INVALID S0XSC LENGTH.') + CALL LCMGET(KPSYS(II),'DRAGON-TXSC',SGAR(0,II)) + CALL LCMGET(KPSYS(II),'DRAGON-S0XSC',SGAS(0,1,II)) +*---- +* TEST FOR FOKKER-PLANCK TREATMENT. +*---- + IF(IBFP.GT.0) THEN + CALL LCMLEN(KPSYS(II),'DRAGON-ESTOP',IFP,ITYLCM) + IF(IFP.NE.(NBMIX+1)*2) CALL XABORT('SNFLUX: INVALID ESTOPW LEN' + 1 //'GTH.') + CALL LCMGET(KPSYS(II),'DRAGON-ESTOP',ESTOPW(0,1,II)) + CALL LCMGET(KPSYS(II),'DRAGON-DELTE',DELTAE(II)) + CALL LCMGET(KPSYS(II),'DRAGON-ISLG',ISLG(II)) + IF(ISLG(II).EQ.1) FLUXC(:)=0.0 + CALL LCMLEN(KPSYS(II),'DRAGON-EMOMT',IFP,ITYLCM) + IF(IFP.NE.NBMIX+1) CALL XABORT('SNFLUX: INVALID EMOMTR LENGTH.') + CALL LCMGET(KPSYS(II),'DRAGON-EMOMT',EMOMTR(0,II)) + ENDIF +*---- +* PRINT ZEROTH MOMENT OF SOURCES. +*---- + IF(IMPX.GT.3) THEN + WRITE(IUNOUT,500) NGIND(II) + WRITE(IUNOUT,'(1P,6(5X,E15.7))') (SUNKNO(KEYFLX(I),II),I=1,NREG) + ENDIF + 10 CONTINUE +*---- +* RECOVER INFORMATION ABOUT BOUNDARY SOURCES +*---- + ISBS=0 + MAXL=0 + IF(SUM(NBS).NE.0) THEN + ISBS=1 + IF(ITYPE.EQ.2) THEN + MAXL=1 + ALLOCATE(ISBSM(2,NLF,NGEFF)) + ELSE IF(ITYPE.EQ.5) THEN + MAXL=MAX(LX,LY) + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + ALLOCATE(ISBSM(4,NPQ,NGEFF)) + ELSE IF(ITYPE.EQ.7) THEN + MAXL=MAX(LX*LY,LX*LZ,LY*LZ) + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + ALLOCATE(ISBSM(6,NPQ,NGEFF)) + ELSE + CALL XABORT('SNFLUX: BOUNDARY SOURCE NOT IMPLEMENTED FOR THAT' + 1 //' GEOMETRY') + ALLOCATE(ISBSM(0,0,0),BS(0,0)) + ENDIF + ALLOCATE(BS(MAXL,SUM(NBS))) + BS=0.0 + ISBSM=0 + JJ=1 + DO 11 II=1,NGEFF + IF(NBS(II).NE.0) THEN + DO N=1,NBS(II) + CALL LCMGDL(KPSOU1(II),N,BS(1,JJ)) + CALL LCMGDL(KPSOU2(II),N,BSINFO) + ISBSM(BSINFO(1),BSINFO(2),II)=JJ + JJ=JJ+1 + ENDDO + ENDIF + 11 CONTINUE + ELSE + ALLOCATE(ISBSM(0,0,0),BS(0,0)) + ENDIF +*---- +* COMPUTE THE FLUX. +*---- + IF((ITYPE.EQ.2).AND.(IBFP.EQ.0)) THEN +*------------ +* 1D SLAB +*------------ + ! EXTRACTING PARAMETERS + IF(IELEM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW' + 1 //'(1)') + CALL LCMGPD(IPTRK,'U',U_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL C_F_POINTER(U_PTR,U,(/ NLF /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NLF*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NLF*NSCT /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! LOOP FOR GROUPS + DO 40 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 40 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-slab' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 30 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 30 + DO 20 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 20 + DO 15 IEL=1,IELEM + IND=(IR-1)*NSCT*IELEM+IELEM*(P-1)+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 15 CONTINUE + 20 CONTINUE + 30 CONTINUE + + ! ONE-SPEED FLUX CALCULATION + CALL SNFBC1(NREG,NBMIX,IELEM,NLF,NSCT,U,MAT,VOL,SGAR(0,II), + 1 NCODE,ZCODE,QEXT(1,II),LFIXUP,LSHOOT,ISBS,SUM(NBS), + 2 ISBSM(1,1,II),BS,WX,CST,LADPT(2),NUN,FUNKNO(1,II),MN,DN) + + 40 CONTINUE ! END OF ENERGY LOOP + + ELSE IF(ITYPE.EQ.2) THEN +*------------ +* 1D SLAB BOLTZMANN-FOKKER-PLANCK +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM*EELEM + IF((IELEM*NLF+NM*NSCT)*NREG.GT.NUN) THEN + CALL XABORT('SNFLUX: QEXT OVERFLOW(1a)') + ENDIF + IOF=NM*NSCT*NREG+1 + CALL LCMGPD(IPTRK,'U',U_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'WE',WE_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL C_F_POINTER(U_PTR,U,(/ NLF /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NLF*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NLF*NSCT /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! ANGULAR FOKKER-PLANCK OPERATOR + ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NLF,NLF)) + SIGMATR=0.0 + DO P=1,NSCT + SIGMATR(P,P)=-IL(P)*(IL(P)+1) + ENDDO + MTR=MATMUL(MATMUL(RESHAPE(MN,[NLF,NSCT]),SIGMATR), + 1 RESHAPE(DN,[NSCT,NLF])) + LAMBDA0=0.0 + DO M=1,NLF + IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) + ENDDO + DO M=1,NLF + MTR(M,M)=MTR(M,M)+LAMBDA0 + ENDDO + SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT,NLF]),MTR), + 1 RESHAPE(MN,[NLF,NSCT])) + DEALLOCATE(MTR) + DO II=1,NGEFF + DO IBM=1,NBMIX + SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) + ENDDO + ENDDO + + ! LOOP FOR GROUPS + DO 80 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 80 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/1D-slab' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 70 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 70 + DO 60 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 60 + DO 50 IEL=1,NM + IND=(IR-1)*NSCT*NM+NM*(P-1)+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + DO Q=1,NSCT + JND=(IR-1)*NSCT*NM+NM*(Q-1)+IEL + QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) + 1 *FUNKNO(JND,II) + ENDDO + 50 CONTINUE + 60 CONTINUE + 70 CONTINUE + + ! ONE-SPEED FLUX CALCULATION + CALL SNFE1D(NREG,NBMIX,IELEM,EELEM,NM,NLF,NSCT,U,MAT,VOL, + 1 SGAR(0,II),ESTOPW(0,1,II),NCODE,ZCODE,DELTAE(II),QEXT(1,II), + 2 LFIXUP,LSHOOT,FUNKNO(1,II),ISBS,SUM(NBS),ISBSM(1,1,II),BS,WX,WE, + 3 CST,LADPT(1:2),IBFP,NUN,MN,DN) + + ! SOLUTION FLUX AT THE CUTOFF ENERGY + IF(ISLG(II).EQ.1) THEN + LXNI=EELEM*NLF + IF(LSHOOT) LXNI=0 + DO 71 IR=1,NREG + IBM=MAT(IR) + FLUXC(IR)=0.0 + DO 72 IP=1,NLF + IND=NM*NSCT*NREG+LXNI+(IR-1)*NLF*IELEM+(IP-1)*IELEM+1 + JND=(IP-1)*NSCT+1 + FLUXC(IR)=FLUXC(IR)+FUNKNO(IND,II)*DN(JND) + 72 CONTINUE + 71 CONTINUE + ENDIF + + 80 CONTINUE ! END OF ENERGY LOOP + DEALLOCATE(SIGMATR) + + ELSE IF(ITYPE.EQ.3) THEN +*------------ +* 1D TUBE/CYLINDRICAL +*------------ + ! EXTRACTING PARAMETERS + IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(1).') + IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(2)') + M2=NLF/2 + CALL LCMLEN(IPTRK,'UPQ',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'JOP',JOP_PTR) + CALL LCMGPD(IPTRK,'U',U_PTR) + CALL LCMGPD(IPTRK,'UPQ',UPQ_PTR) + CALL LCMGPD(IPTRK,'WPQ',WPQ_PTR) + CALL LCMGPD(IPTRK,'ALPHA',ALPHA_PTR) + CALL LCMGPD(IPTRK,'PLZ',PLZ_PTR) + CALL LCMGPD(IPTRK,'PL',PL_PTR) + CALL LCMGPD(IPTRK,'SURF',SURF_PTR) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + CALL C_F_POINTER(JOP_PTR,JOP,(/ M2 /)) + CALL C_F_POINTER(U_PTR,U,(/ NPQ /)) + CALL C_F_POINTER(UPQ_PTR,UPQ,(/ NPQ /)) + CALL C_F_POINTER(WPQ_PTR,WPQ,(/ NPQ /)) + CALL C_F_POINTER(ALPHA_PTR,ALPHA,(/ NPQ /)) + CALL C_F_POINTER(PLZ_PTR,PLZ,(/ NSCT*M2 /)) + CALL C_F_POINTER(PL_PTR,PL,(/ NSCT*NPQ /)) + CALL C_F_POINTER(SURF_PTR,SURF,(/ NREG+1 /)) + + ! LOOP FOR GROUPS + DO 120 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 120 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-cyl' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 110 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 110 + IOF=0 + DO 100 P=0,NANI-1 + DO 90 IM=0,P + IF(MOD(P+IM,2).EQ.1) GO TO 90 + IOF=IOF+1 + IND=(IR-1)*NSCT+IOF + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,P+1,II)*FUNKNO(IND,II) + 90 CONTINUE + 100 CONTINUE + 110 CONTINUE + + ! ONE-SPEED FLUX CALCULATION + CURR=0.0 + CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, + 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR, + 2 FUNKNO(1,II)) + IF(ZCODE(2).NE.0.0) THEN + CA=CURR + CB=1.0 + CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, + 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CB, + 2 FUNKNO(1,II)) + CURR=ZCODE(2)*CA/(1.0+ZCODE(2)*(CA-CB)) + CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, + 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR, + 2 FUNKNO(1,II)) + ENDIF + + 120 CONTINUE ! END OF ENERGY LOOP + + ELSE IF(ITYPE.EQ.4) THEN +*------------ +* 1D SPHERE +*------------ + ! EXTRACTING PARAMETERS + IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(2).') + NSCT=ISCAT + IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(3)') + CALL LCMGPD(IPTRK,'U',U_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'ALPHA',ALPHA_PTR) + CALL LCMGPD(IPTRK,'PLZ',PLZ_PTR) + CALL LCMGPD(IPTRK,'PL',PL_PTR) + CALL LCMGPD(IPTRK,'SURF',SURF_PTR) + CALL LCMGPD(IPTRK,'XXX',XXX_PTR) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + CALL C_F_POINTER(U_PTR,U,(/ NLF /)) + CALL C_F_POINTER(W_PTR,W,(/ NLF /)) + CALL C_F_POINTER(ALPHA_PTR,ALPHA,(/ NLF /)) + CALL C_F_POINTER(PLZ_PTR,PLZ,(/ NSCT /)) + CALL C_F_POINTER(PL_PTR,PL,(/ NSCT*NLF /)) + CALL C_F_POINTER(SURF_PTR,SURF,(/ NREG+1 /)) + CALL C_F_POINTER(XXX_PTR,XXX,(/ NREG+1 /)) + + ! LOOP FOR GROUPS + DO 150 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 150 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-sph' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 140 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 140 + DO 130 P=0,NANI-1 + IND=(IR-1)*NSCT+P+1 + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,P+1,II)*FUNKNO(IND,II) + 130 CONTINUE + 140 CONTINUE + + ! ONE-SPEED FLUX CALCULATION + CURR=0.0 + CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL,SURF, + 1 XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR,FUNKNO(1,II)) + IF(ZCODE(2).NE.0.0) THEN + CA=CURR + CB=1.0 + CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL, + 1 SURF,XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CB,FUNKNO(1,II)) + CURR=ZCODE(2)*CA/(1.0+ZCODE(2)*(CA-CB)) + CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL, + 1 SURF,XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR,FUNKNO(1,II)) + ENDIF + + 150 CONTINUE ! END OF ENERGY LOOP + + ELSE IF((ITYPE.EQ.5).AND.(IBFP.EQ.0)) THEN +*------------ +* 2D CARTESIAN +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM**2 + NMX=IELEM + IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVE' + 1 //'RFLOW(4a)') + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRM',MRM_PTR) + CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) + CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! LOOP FOR GROUPS + DO 200 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 200 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-car' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 190 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 190 + IOF=0 + DO 180 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 180 + DO 160 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 160 CONTINUE + 180 CONTINUE + 190 CONTINUE + 200 CONTINUE ! END OF ENERGY LOOP + IF(NKBA.EQ.0)THEN + CALL SNFBC2(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM,NM,NMX, + 1 NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, + 2 DU,DE,W,MRM,MRMY,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3), + 3 ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) + ELSE IF(NKBA.GT.0)THEN + IF(ISBS.EQ.1) CALL XABORT('SNFLUX: BOUNDARY SOURCES NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + CALL SNFKC2(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM,NM, + 1 NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, + 2 DU,DE,W,MRM,MRMY,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3), + 3 ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) + ENDIF + ELSE IF(ITYPE.EQ.5) THEN +*------------ +* 2D CARTESIAN BOLTZMANN-FOKKER-PLANCK +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM**2*EELEM + NME=IELEM**2 + NMX=IELEM*EELEM + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + IF((NSCT*NM+NPQ*NME)*NREG.GT.NUN) THEN + CALL XABORT('SNFLUX: QEXT OVERFLOW(4b)') + ENDIF + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRM',MRM_PTR) + CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'WE',WE_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) + CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! ANGULAR FOKKER-PLANCK OPERATOR + ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NPQ,NPQ)) + SIGMATR=0.0 + DO P=1,NSCT + SIGMATR(P,P)=-IL(P)*(IL(P)+1) + ENDDO + MTR=MATMUL(MATMUL(RESHAPE(MN,[NPQ,NSCT]),SIGMATR), + 1 RESHAPE(DN,[NSCT,NPQ])) + LAMBDA0=0.0 + DO M=1,NPQ + IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) + ENDDO + DO M=1,NPQ + MTR(M,M)=MTR(M,M)+LAMBDA0 + ENDDO + SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT, NPQ]),MTR), + 1 RESHAPE(MN,[NPQ,NSCT])) + DEALLOCATE(MTR) + DO II=1,NGEFF + DO IBM=1,NBMIX + SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) + ENDDO + ENDDO + + ! LOOP FOR GROUPS + DO 250 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 250 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/2D-car' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 240 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 240 + IOF=0 + DO 230 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 230 + DO 210 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + DO Q=1,NSCT + JND=(IR-1)*NSCT*NM+(Q-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) + 1 *FUNKNO(JND,II) + ENDDO + 210 CONTINUE + 230 CONTINUE + 240 CONTINUE + 250 CONTINUE ! END OF ENERGY LOOP + DEALLOCATE(SIGMATR) + + ! FLUX CALCULATION + CALL SNFE2D(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM, + 1 EELEM,NM,NME,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,ESTOPW, + 2 NCODE,ZCODE,DELTAE,QEXT,LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,FUNKNO, + 3 ISLG,FLUXC,ISBS,SUM(NBS),ISBSM,BS,MAXL,WX,WX,WE,CST,LADPT(1:3), + 4 IBFP,MN,DN) + + ELSE IF(ITYPE.EQ.6) THEN +*------------ +* TUBE 2D (R-Z) +*------------ + ! EXTRACTING PARAMETERS + IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(3).') + IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(5)') + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRM',MRM_PTR) + CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'DAL',DAL_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) + CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DAL_PTR,DAL,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! LOOP FOR GROUPS + DO 290 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 290 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-rz' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 280 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 280 + IOF=0 + DO 270 P=1,NSCT + IF(IL(P).GT.MIN(ISCAT,NANI)-1) GO TO 270 + IND=(IR-1)*NSCT+(P-1)*IELEM*IELEM+1 + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 270 CONTINUE + 280 CONTINUE + + ! ONE-SPEED FLUX CALCULATION + CALL SNFC12(LX,LY,NBMIX,NPQ,NSCT,MAT,VOL,SGAR(0,II),NCODE,ZCODE, + 1 QEXT(1,II),LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,DAL,FUNKNO(1,II), + 2 FUNKNO(L4+1,II),FUNKNO(L4+IELEM*LY*NPQ+1,II),MN,DN) + + 290 CONTINUE ! END OF ENERGY GROUP + + ELSE IF((ITYPE.EQ.7).AND.(IBFP.EQ.0)) THEN +*---- +* 3D CARTESIAN CASE +*---- + ! EXTRACTING PARAMETERS + NM=IELEM**3 + NMX=IELEM**2 + IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QE' + 1 //'XT OVERFLOW(6a)') + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'DZ',DZ_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRMX',MRMX_PTR) + CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) + CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) + CALL LCMGPD(IPTRK,'DC',DC_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRMX_PTR,MRMX,(/ NPQ /)) + CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) + CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) + CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! LOOP FOR GROUPS + DO 340 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 340 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/3D-car' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 330 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 330 + IOF=0 + DO 320 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 320 + DO 300 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 300 CONTINUE + 320 CONTINUE + 330 CONTINUE + 340 CONTINUE ! END OF ENERGY LOOP + + ! FLUX CALCULATION + IF(NKBA.EQ.0)THEN + CALL SNFBC3(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ,IELEM,NM,NMX, + 1 NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, + 2 DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, + 3 LADPT(2:4),ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) + ELSE IF(NKBA.GT.0)THEN + IF(ISBS.EQ.1) CALL XABORT('SNFLUX: BOUNDARY SOURCES NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + CALL SNFKC3(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ,IELEM, + 1 NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT, + 2 LFIXUP,DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, + 3 LADPT(2:4),ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) + ENDIF + ELSE IF(ITYPE.EQ.7) THEN +*------------ +* 3D CARTESIAN BOLTZMANN-FOKKER-PLANCK +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM**3*EELEM + NME=IELEM**3 + NMX=IELEM**2*EELEM + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + IF((NSCT*NM+NPQ*NME)*NREG.GT.NUN) THEN + CALL XABORT('SNFLUX: QEXT OVERFLOW(6b)') + ENDIF + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'DZ',DZ_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRMX',MRMX_PTR) + CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) + CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) + CALL LCMGPD(IPTRK,'DC',DC_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'WE',WE_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRMX_PTR,MRMX,(/ NPQ /)) + CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) + CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) + CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + + ! ANGULAR FOKKER-PLANCK OPERATOR + ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NPQ,NPQ)) + SIGMATR=0.0 + DO P=1,NSCT + SIGMATR(P,P)=-IL(P)*(IL(P)+1) + ENDDO + MTR=MATMUL(MATMUL(RESHAPE(MN,[NPQ, NSCT]),SIGMATR), + 1 RESHAPE(DN,[NSCT,NPQ])) + LAMBDA0=0.0 + DO M=1,NPQ + IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) + ENDDO + DO M=1,NPQ + MTR(M,M)=MTR(M,M)+LAMBDA0 + ENDDO + SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT,NPQ]),MTR), + 1 RESHAPE(MN,[NPQ,NSCT])) + DEALLOCATE(MTR) + DO II=1,NGEFF + DO IBM=1,NBMIX + SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) + ENDDO + ENDDO + + ! LOOP FOR GROUPS + DO 390 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 390 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/3D-car' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 380 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 380 + IOF=0 + DO 370 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 370 + DO 350 IEL=1,NM + IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + DO Q=1,NSCT + JND=(IR-1)*NSCT*NM+(Q-1)*NM+IEL + QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) + 1 *FUNKNO(JND,II) + ENDDO + 350 CONTINUE + 370 CONTINUE + 380 CONTINUE + 390 CONTINUE ! END OF ENERGY LOOP + DEALLOCATE(SIGMATR) + + ! FLUX CALCULATION + CALL SNFE3D(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ, + 1 IELEM,EELEM,NM,NME,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR, + 2 ESTOPW,NCODE,ZCODE,DELTAE,QEXT,LFIXUP,DU,DE,DZ,W,MRMX,MRMY,MRMZ, + 3 DC,DB,DA,FUNKNO,ISLG,FLUXC,ISBS,SUM(NBS),ISBSM,BS,MAXL, + 4 WX,WX,WX,WE,CST,LADPT,IBFP,MN,DN) + + ELSE IF(ITYPE.EQ.8) THEN +*------------ +* 2D HEXAGONAL +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM**2 + NMX=IELEM + IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVE' + 1 //'RFLOW(7)') + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'SIDE',SIDE) + CALL LCMGET(IPTRK,'LOZSWP',LOZSWP) + NHEX=LX/(3*ISPLH**2) + ALLOCATE(COORDMAP(3,NHEX)) + CALL LCMGET(IPTRK,'COORDMAP',COORDMAP) + + ! LOOP FOR GROUPS + DO 440 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 440 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-hex' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 430 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 430 + IOF=0 + DO 420 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 420 + DO 400 IEL=1,NM + IND=(((IR-1)*NSCT+(P-1))*NM)+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 400 CONTINUE + 420 CONTINUE + 430 CONTINUE + 440 CONTINUE ! END OF ENERGY LOOP + + ! FLUX CALCULATION + IF(NKBA.EQ.0)THEN + CALL SNFBH2(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,ISPLH,SIDE, + 1 IELEM,NM,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,QEXT,LFIXUP,DU, + 2 DE,W,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3),LOZSWP,COORDMAP, + 3 FUNKNO) + ELSE IF(NKBA.GE.1)THEN + IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + CALL SNFKH2(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,ISPLH,SIDE, + 1 IELEM,NM,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,QEXT,LFIXUP,DU, + 2 DE,W,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3),LOZSWP,COORDMAP, + 3 FUNKNO) + + ENDIF + DEALLOCATE(COORDMAP) + + ELSE IF(ITYPE.EQ.9) THEN +*------------ +* 3D HEXAGONAL +*------------ + ! EXTRACTING PARAMETERS + NM=IELEM**3 + NMX=IELEM**2 + IF(IELEM*IELEM*IELEM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QE' + 1 //'XT OVERFLOW(8)') + CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) + CALL LCMGPD(IPTRK,'DU',DU_PTR) + CALL LCMGPD(IPTRK,'DE',DE_PTR) + CALL LCMGPD(IPTRK,'DZ',DZ_PTR) + CALL LCMGPD(IPTRK,'W',W_PTR) + CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) + CALL LCMGPD(IPTRK,'DC',DC_PTR) + CALL LCMGPD(IPTRK,'DB',DB_PTR) + CALL LCMGPD(IPTRK,'DA',DA_PTR) + CALL LCMGPD(IPTRK,'WX',WX_PTR) + CALL LCMGPD(IPTRK,'CST',CST_PTR) + CALL LCMGPD(IPTRK,'IL',IL_PTR) + CALL LCMGPD(IPTRK,'MN',MN_PTR) + CALL LCMGPD(IPTRK,'DN',DN_PTR) + CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) + CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) + CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) + CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) + CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) + CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) + CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) + CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) + CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) + CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) + CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) + CALL LCMGET(IPTRK,'NCODE',NCODE) + CALL LCMGET(IPTRK,'ZCODE',ZCODE) + CALL LCMGET(IPTRK,'SIDE',SIDE) + CALL LCMGET(IPTRK,'LOZSWP',LOZSWP) + ! Number of hexagons in one plane only + NHEX=LX/(3*ISPLH**2) + ALLOCATE(COORDMAP(3,NHEX)) + CALL LCMGET(IPTRK,'COORDMAP',COORDMAP) + + ! LOOP FOR GROUPS + DO 490 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 490 + IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/3D-hex' + + ! FIXED VOLUMIC SOURCES + QEXT(:NUN,II)=SUNKNO(:NUN,II) + + ! IN-GROUP SCATTERING + DO 480 IR=1,NREG + IBM=MAT(IR) + IF(IBM.EQ.0) GO TO 480 + IOF=0 + DO 470 P=1,NSCT + IF(IL(P).GT.NANI-1) GO TO 470 + DO 450 IEL=1,IELEM*IELEM*IELEM + IND=(((IR-1)*NSCT+(P-1))*IELEM**3)+IEL + QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) + 450 CONTINUE + 470 CONTINUE + 480 CONTINUE + 490 CONTINUE ! END OF ENERGY LOOP + + ! FLUX CALCULATION + IF(NKBA.EQ.0)THEN + CALL SNFBH3(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH,SIDE, + 1 IELEM,NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE, + 2 QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, + 3 LADPT(2:4),LOZSWP,COORDMAP,FUNKNO) + ELSE IF(NKBA.GE.1)THEN + IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' + 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' + 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') + CALL SNFKH3(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH, + 1 SIDE,IELEM,NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE, + 2 ZCODE,QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX, + 3 CST,LADPT(2:4),LOZSWP,COORDMAP,FUNKNO) + ENDIF + DEALLOCATE(COORDMAP) + ELSE + CALL XABORT('SNFLUX: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') + ENDIF + DEALLOCATE(ESTOPW,SGAS,SGAR) +*---- +* PRINT COMPLETE UNKNOWN VECTOR. +*---- + DO 495 II=1,NGEFF + IF(.NOT.INCONV(II)) GO TO 495 + IF(IMPX.GT.5) THEN + WRITE(IUNOUT,520) NGIND(II) + WRITE(IUNOUT,'(1P,4(5X,E15.7))') (FUNKNO(:,II)) + ENDIF + 495 CONTINUE +*---- +* DIFFUSION SYNTHETIC ACCELERATION. +*---- + IF(LDSA) THEN + ISOLVSA=IPAR(33) + + CALL LCMSIX(IPTRK,'DSA',1) + CALL LCMGET(IPTRK,'KEYFLX',KEYSPN) + CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR) + IF(NREG.NE.IPAR(1)) CALL XABORT('SNFLUX: INVALID NREG (DSA).') + NUNSA=IPAR(2) + ITYPE=IPAR(6) + + IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN + IELEMSA=IPAR(9) + ELSE + IF(ISOLVSA.EQ.1)THEN + IELEMSA=IPAR(8) + ELSEIF(ISOLVSA.EQ.2)THEN + IELEMSA=IPAR(9) + ENDIF + ENDIF + + IMPY=MAX(0,IMPX-1) + CALL LCMSIX(IPTRK,' ',2) + + CALL SNDSA(KPSYS,INCONV,NGIND,IPTRK,IMPY,NGEFF,NREG,NBMIX, + 1 NUN,MAT,VOL,KEYFLX,KEYSPN,NUNSA,IELEMSA,ZCODE,OLD,FUNKNO, + 2 NHEX) + + DEALLOCATE(OLD) + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(DELTAE,KEYSPN,QEXT,BS,ISBSM) + RETURN +* + 500 FORMAT(//41H SNFLUX: N E U T R O N S O U R C E S (,I5,3H ):) + 510 FORMAT(/25H SNFLUX: PROCESSING GROUP,I5,6H WITH ,A,1H.) + 520 FORMAT(//41H SNFLUX: A F T E R T R A N S P O R T (,I5,3H ):) + END |