diff options
Diffstat (limited to 'Donjon/src/DRESOU.f')
| -rw-r--r-- | Donjon/src/DRESOU.f | 167 |
1 files changed, 167 insertions, 0 deletions
diff --git a/Donjon/src/DRESOU.f b/Donjon/src/DRESOU.f new file mode 100644 index 0000000..e82d4e3 --- /dev/null +++ b/Donjon/src/DRESOU.f @@ -0,0 +1,167 @@ +*DECK DRESOU + SUBROUTINE DRESOU(IPRINT,IPGPT,IPMAC1,IPMAC2,IPFLX,IPGRAD,NG,NREG, + 1 NMIL,NUN,MATCOD,KEYFLX,VOL,LNO,RMSD) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the GPT sources corresponding to the gradient of the RMS power +* distribution. Case with no direct effect. +* +*Copyright: +* Copyright (C) 2012 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 +* IPRINT print parameter +* IPGPT pointer to the L_GPT data structure. +* IPMAC1 pointer to the actual macrolib structure. +* IPMAC2 pointer to the reference macrolib structure. +* IPFLX pointer to the multigroup flux. +* IPGRAD pointer to the L_OPTIMIZE object. +* NG number of energy groups. +* NREG number of regions. +* NMIL number of material mixtures. +* NUN number of unknowns per energy group. +* MATCOD material mixture indices per region. +* KEYFLX position of averaged fluxes in unknown vector. +* VOL volumes. +* LNO flag set to .true. to exit after calculation of RMS. +* +*Parameters: output +* RMSD RMS error on power distribution. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPGPT,IPMAC1,IPMAC2,IPFLX,IPGRAD + INTEGER IPRINT,NG,NREG,NMIL,NUN,MATCOD(NREG),KEYFLX(NREG) + REAL VOL(NREG) + DOUBLE PRECISION RMSD + LOGICAL LNO +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) JPMAC1,JPMAC2,KPMAC1,KPMAC2,JPFLX,JPGPT,KPGPT + DOUBLE PRECISION SOUT2,SOUTOT,PW1TOT,PW2TOT,DSUM,AIL,BIL + REAL, ALLOCATABLE, DIMENSION(:) :: POW1,H1,F1,POW2,H2,F2,SUNK, + 1 FLUX + CHARACTER HSMG*131 +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(POW1(NMIL),H1(NMIL),F1(NMIL),POW2(NMIL),H2(NMIL), + 1 F2(NMIL)) +*---- +* COMPUTE THE ACTUAL AND REFERENCE POWER DISTRIBUTION +*---- + POW1(:NMIL)=0.0 + POW2(:NMIL)=0.0 + JPMAC1=LCMGID(IPMAC1,'GROUP') + JPMAC2=LCMGID(IPMAC2,'GROUP') + DO IG=1,NG + KPMAC1=LCMGIL(JPMAC1,IG) + KPMAC2=LCMGIL(JPMAC2,IG) + CALL LCMLEN(KPMAC1,'FLUX-INTG',ILG,ITYLCM) + IF(ILG.EQ.0) CALL XABORT('DRESOU: MISSING ACTUAL FLUX.') + CALL LCMLEN(KPMAC2,'FLUX-INTG',ILG,ITYLCM) + IF(ILG.EQ.0) CALL XABORT('DRESOU: MISSING REFERENCE FLUX.') + CALL LCMLEN(KPMAC1,'H-FACTOR',ILG,ITYLCM) + IF(ILG.EQ.0) CALL XABORT('DRESOU: MISSING ACTUAL H-FACTOR.') + CALL LCMLEN(KPMAC2,'H-FACTOR',ILG,ITYLCM) + IF(ILG.EQ.0) CALL XABORT('DRESOU: MISSING REFERENCE H-FACTOR.') + CALL LCMGET(KPMAC1,'FLUX-INTG',F1) + CALL LCMGET(KPMAC2,'FLUX-INTG',F2) + CALL LCMGET(KPMAC1,'H-FACTOR',H1) + CALL LCMGET(KPMAC2,'H-FACTOR',H2) + DO IBM=1,NMIL + POW1(IBM)=POW1(IBM)+F1(IBM)*H1(IBM) + POW2(IBM)=POW2(IBM)+F2(IBM)*H2(IBM) + ENDDO + ENDDO +*---- +* COMPUTE THE RMS FUNCTIONAL +*---- + PW1TOT=0.0D0 + PW2TOT=0.0D0 + DO IBM=1,NMIL + PW1TOT=PW1TOT+POW1(IBM) + PW2TOT=PW2TOT+POW2(IBM) + ENDDO + RMSD=0.0D0 + DO IBM=1,NMIL + RMSD=RMSD+(POW1(IBM)/PW1TOT-POW2(IBM)/PW2TOT)**2 + ENDDO + CALL LCMPUT(IPGRAD,'FOBJ-CST-VAL',1,4,RMSD) + IF(IPRINT.GT.0) WRITE(6,100) RMS + IF(LNO) GO TO 10 +*---- +* COMPUTE THE GRADIENT OF THE RMS FUNCTIONAL +*---- + ALLOCATE(SUNK(NUN)) + JPFLX=LCMGID(IPFLX,'FLUX') + JPGPT=LCMLID(IPGPT,'ASOUR',1) + KPGPT=LCMLIL(JPGPT,1,NG) + DO IG=1,NG + SUNK(:NUN)=0.0 + KPMAC1=LCMGIL(JPMAC1,IG) + CALL LCMGET(KPMAC1,'H-FACTOR',H1) + DO IR=1,NREG + IUNK=KEYFLX(IR) + IF(IUNK.EQ.0) CYCLE + IBM=MATCOD(IR) + IF(IBM.EQ.0) CYCLE + SOUT2=0.0D0 + DO JBM=1,NMIL + SOUTOT=0.0D0 + IF(IBM.EQ.JBM) SOUTOT=1.0D0 + SOUTOT=SOUTOT-POW1(JBM)/PW1TOT + SOUT2=SOUT2+SOUTOT*(POW1(JBM)/PW1TOT-POW2(JBM)/PW2TOT) + ENDDO + SUNK(IUNK)=2.0*VOL(IR)*H1(IBM)*REAL(SOUT2/PW1TOT) + ENDDO + CALL LCMPDL(KPGPT,IG,NUN,2,SUNK) + ENDDO +*---- +* CHECK SOURCE ORTHOGONALITY +*---- + ALLOCATE(FLUX(NUN)) + AIL=0.0D0 + BIL=0.0D0 + DO IG=1,NG + CALL LCMGDL(KPGPT,IG,SUNK) + CALL LCMGDL(JPFLX,IG,FLUX) + DO IUNK=1,NUN + GAZ=FLUX(IUNK)*SUNK(IUNK) + DAZ=FLUX(IUNK)**2 + AIL=AIL+GAZ + BIL=BIL+DAZ + ENDDO + ENDDO + DSUM=ABS(AIL)/ABS(BIL)/REAL(NUN) + IF(IPRINT.GT.0) THEN + WRITE(6,'(/21H DRESOU: DOT PRODUCT=,1P,E11.4)') DSUM + ENDIF + IF(ABS(DSUM).GT.1.0E-5) THEN + WRITE(HSMG,'(36HDRESOU: NON ORTHOGONAL SOURCE (DSUM=,1P,E11.3, + 1 2H).)') DSUM + CALL XABORT(HSMG) + ENDIF + DEALLOCATE(FLUX,SUNK) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + 10 DEALLOCATE(F2,H2,POW2,F1,H1,POW1) + RETURN +* + 100 FORMAT(/41H DRESOU: RMS ERROR ON POWER DISTRIBUTION=,1P,E11.4) + END |