diff options
Diffstat (limited to 'Dragon/src/EDIWP1.f')
| -rw-r--r-- | Dragon/src/EDIWP1.f | 233 |
1 files changed, 233 insertions, 0 deletions
diff --git a/Dragon/src/EDIWP1.f b/Dragon/src/EDIWP1.f new file mode 100644 index 0000000..488dc75 --- /dev/null +++ b/Dragon/src/EDIWP1.f @@ -0,0 +1,233 @@ +*DECK EDIWP1 + SUBROUTINE EDIWP1(IPFLUX,NW,NGROUP,NUN,NREGIO,NDIM,IADJ,NLIN, + > NFUNL,NGCOND,NMERGE,KEYANI,VOLUME,IGCOND,IMERGE,FLUXES,AFLUXE) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Evaluate the PN weighting spectra for an homogenization based on +* spherical harmonic moments of the flux. +* +*Copyright: +* Copyright (C) 2019 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 +* IPFLUX pointer to the flux LCM object. +* NW order of the spherical harmonic expansion for the flux. +* NGROUP number of energy groups. +* NUN number of unknowns in flux array. +* NREGIO number of regions. +* NDIM number of dimensions. +* IADJ type of flux weighting: +* =0: direct flux weighting; +* =1: direct-adjoint flux weighting. +* NLIN number of polynomial components in flux. +* NFUNL number of spherical harmonic components in flux. +* NGCOND number of merged energy groups. +* NMERGE number of merged regions. +* KEYANI position of spherical harmonic components in unknown vector. +* VOLUME volumes. +* IGCOND limit condensed groups. +* IMERGE region merging matrix. +* +*Parameters: input/output +* FLUXES weighting function for PN fluxes. +* AFLUXE weighting function for PN adjoint fluxes. +* +*Reference: +* Jean-Francois Vidal et al., APOLLO3 homogenization techniques for +* transport core calculations - application to the ASTRID CFV core, +* Nuclear Engineering and Technology 49 (2017) 1379 - 1387. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPFLUX + INTEGER NW,NREGIO,NGROUP,NUN,NDIM,IADJ,NLIN,NFUNL,NGCOND, + > NMERGE,KEYANI(NREGIO,NLIN,NFUNL),IGCOND(NGCOND), + > IMERGE(NREGIO) + REAL VOLUME(NREGIO),FLUXES(NREGIO,NGROUP,NW), + > AFLUXE(NREGIO,NGROUP,NW) +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) JPFLUX,JPFLUA + DOUBLE PRECISION DVOL +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, ALLOCATABLE, DIMENSION(:) :: WORKF,WORKA + REAL, ALLOCATABLE, DIMENSION(:,:) :: FDEN,ADEN + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: FLUANI,AFLANI + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: SVOL + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: DFLUA,DAFLA +*---- +* INITIALIZATION +*---- + IF(NFUNL.EQ.1) CALL XABORT('EDIWP1: ANIS.GE.2 EXPECTED IN TRACKI' + > //'NG.') + JPFLUX=LCMGID(IPFLUX,'FLUX') + ALLOCATE(WORKF(NUN)) + IF(IADJ.EQ.1) THEN + JPFLUA=LCMGID(IPFLUX,'AFLUX') + ALLOCATE(WORKA(NUN)) + ENDIF +*---- +* PROCESS TRIVIAL 1D CASE +*---- + IF(NDIM.EQ.1) THEN + DO IL=1,NW + DO IGR=1,NGROUP + IF(IADJ.EQ.0) THEN + CALL LCMGDL(JPFLUX,IGR,WORKF) + DO IREG=1,NREGIO + FLUXES(IREG,IGR,IL)=WORKF(KEYANI(IREG,1,IL+1)) + AFLUXE(IREG,IGR,IL)=1.0 + ENDDO + ELSE IF(IADJ.EQ.1) THEN + CALL LCMGDL(JPFLUX,IGR,WORKF) + CALL LCMGDL(JPFLUA,IGR,WORKA) + DO IREG=1,NREGIO + FLUXES(IREG,IGR,IL)=WORKF(KEYANI(IREG,1,IL+1)) + AFLUXE(IREG,IGR,IL)=WORKA(KEYANI(IREG,1,IL+1)) + ENDDO + ENDIF + DO IREG=1,NREGIO + FLUXES(IREG,IGR,IL)=MAX(ABS(FLUXES(IREG,IGR,IL)),1.0E-10) + AFLUXE(IREG,IGR,IL)=MAX(ABS(AFLUXE(IREG,IGR,IL)),1.0E-10) + ENDDO + ENDDO + ENDDO + DEALLOCATE(WORKF) + IF(IADJ.EQ.1) DEALLOCATE(WORKA) + RETURN + ENDIF +*---- +* RECOVER PN MOMENTS OF THE FLUX +*---- + IOF=1 + DO IL=1,NW + IOF0=IOF + NTRM=1 + IF(NDIM.EQ.2) THEN + NTRM=IL+1 + ELSE IF(NDIM.EQ.3) THEN + NTRM=2*IL+1 + ENDIF + ALLOCATE(FLUANI(NREGIO,NGROUP,NTRM),AFLANI(NREGIO,NGROUP,NTRM)) + DO IGR=1,NGROUP + IF(IADJ.EQ.0) THEN + CALL LCMGDL(JPFLUX,IGR,WORKF) + DO IREG=1,NREGIO + IOF=IOF0 + ID=0 + DO IM=-IL,IL + IF((NDIM.EQ.2).AND.(MOD(IL+IM,2).EQ.1)) CYCLE + IOF=IOF+1 + IF(IOF.GT.NFUNL) CALL XABORT('EDIWP1: KEYANI OVERFLOW.') + ID=ID+1 + FLUANI(IREG,IGR,ID)=WORKF(KEYANI(IREG,1,IOF)) + AFLANI(IREG,IGR,ID)=1.0 + ENDDO + ENDDO + ELSE IF(IADJ.EQ.1) THEN + CALL LCMGDL(JPFLUX,IGR,WORKF) + CALL LCMGDL(JPFLUA,IGR,WORKA) + DO IREG=1,NREGIO + IOF=IOF0 + ID=0 + DO IM=-IL,IL + IF((NDIM.EQ.2).AND.(MOD(IL+IM,2).EQ.1)) CYCLE + IOF=IOF+1 + IF(IOF.GT.NFUNL) CALL XABORT('EDIWP1: KEYANI OVERFLOW.') + ID=ID+1 + FLUANI(IREG,IGR,ID)=WORKF(KEYANI(IREG,1,IOF)) + AFLANI(IREG,IGR,ID)=WORKA(KEYANI(IREG,1,IOF)) + ENDDO + ENDDO + ENDIF + ENDDO +*---- +* CONDENSATION AND HOMOGENIZATION OF SPHERICAL HARMONIC MOMENTS +*---- + ALLOCATE(DFLUA(NMERGE,NGCOND,NTRM),DAFLA(NMERGE,NGCOND,NTRM), + > SVOL(NMERGE)) + DFLUA(:NMERGE,:NGCOND,:NTRM)=0.0D0 + DAFLA(:NMERGE,:NGCOND,:NTRM)=0.0D0 + IGRFIN=0 + DO IGRC=1,NGCOND + IGRDEB=IGRFIN+1 + IGRFIN=IGCOND(IGRC) + DO IGR=IGRDEB,IGRFIN + SVOL(:NMERGE)=0.0D0 + DO IREG=1,NREGIO + IRA=IMERGE(IREG) + IF(IRA.EQ.0) CYCLE + DVOL=VOLUME(IREG) + SVOL(IRA)=SVOL(IRA)+DVOL + DO ID=1,NTRM + DFLUA(IRA,IGRC,ID)=DFLUA(IRA,IGRC,ID)+ + > FLUANI(IREG,IGR,ID)*DVOL + DAFLA(IRA,IGRC,ID)=DAFLA(IRA,IGRC,ID)+ + > FLUANI(IREG,IGR,ID)*AFLANI(IREG,IGR,ID)*DVOL + ENDDO + ENDDO + DO IRA=1,NMERGE + DO ID=1,NTRM + DFLUA(IRA,IGRC,ID)=DFLUA(IRA,IGRC,ID)/SVOL(IRA) + DAFLA(IRA,IGRC,ID)=DAFLA(IRA,IGRC,ID)/ + > (DFLUA(IRA,IGRC,ID)*SVOL(IRA)) + ENDDO + ENDDO + ENDDO + ENDDO +*---- +* USE APOLLO3 FORMULA +*---- + ALLOCATE(FDEN(NREGIO,NGROUP),ADEN(NREGIO,NGROUP)) + FLUXES(:NREGIO,:NGROUP,IL)=0.0D0 + AFLUXE(:NREGIO,:NGROUP,IL)=0.0D0 + FDEN(:NREGIO,:NGROUP)=0.0D0 + ADEN(:NREGIO,:NGROUP)=0.0D0 + IGRFIN=0 + DO IGRC=1,NGCOND + IGRDEB=IGRFIN+1 + IGRFIN=IGCOND(IGRC) + DO IGR=IGRDEB,IGRFIN + DO IREG=1,NREGIO + IRA=IMERGE(IREG) + IF(IRA.EQ.0) CYCLE + DO ID=1,NTRM + FLUXES(IREG,IGR,IL)=FLUXES(IREG,IGR,IL)+ + > REAL(FLUANI(IREG,IGR,ID)*DFLUA(IRA,IGRC,ID)) + AFLUXE(IREG,IGR,IL)=AFLUXE(IREG,IGR,IL)+ + > REAL(AFLANI(IREG,IGR,ID)*DAFLA(IRA,IGRC,ID)) + FDEN(IREG,IGR)=FDEN(IREG,IGR)+REAL(DFLUA(IRA,IGRC,ID)) + ADEN(IREG,IGR)=ADEN(IREG,IGR)+REAL(DAFLA(IRA,IGRC,ID)) + ENDDO + ENDDO + ENDDO + ENDDO + DO IGR=1,NGROUP + DO IREG=1,NREGIO + FLUXES(IREG,IGR,IL)=FLUXES(IREG,IGR,IL)/FDEN(IREG,IGR) + AFLUXE(IREG,IGR,IL)=AFLUXE(IREG,IGR,IL)/ADEN(IREG,IGR) + FLUXES(IREG,IGR,IL)=MAX(ABS(FLUXES(IREG,IGR,IL)),1.0E-10) + AFLUXE(IREG,IGR,IL)=MAX(ABS(AFLUXE(IREG,IGR,IL)),1.0E-10) + ENDDO + ENDDO + DEALLOCATE(ADEN,FDEN,SVOL,DAFLA,DFLUA,AFLANI,FLUANI) + ENDDO + DEALLOCATE(WORKF) + IF(IADJ.EQ.1) DEALLOCATE(WORKA) + RETURN + END |