diff options
Diffstat (limited to 'Trivac/src/GPTFLU.f')
| -rwxr-xr-x | Trivac/src/GPTFLU.f | 398 |
1 files changed, 398 insertions, 0 deletions
diff --git a/Trivac/src/GPTFLU.f b/Trivac/src/GPTFLU.f new file mode 100755 index 0000000..03a8ce0 --- /dev/null +++ b/Trivac/src/GPTFLU.f @@ -0,0 +1,398 @@ +*DECK GPTFLU + SUBROUTINE GPTFLU(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Computes generalized adjoints. +* GPTFLU = Generalized Perturbation Theory FLUx calculation +* +*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): A. Hebert, E. Varin and R. Chambon +* +*Parameters: input/ouput +* NENTRY number of linked lists or files used by the module +* HENTRY character*12 name of each linked list or file +* HENTRY(1): create or modification type(L_FLUX) (GPT solution) +* HENTRY(2): read-only type(L_SOURCE) => GPT fixed source +* HENTRY(3): read-only type(L_FLUX) => unperturbed solution +* HENTRY(4): read-only type(L_SYSTEM) => reference matrices +* HENTRY(5): read-only type(L_TRACK) => TRIVAC tracking. +* IENTRY =1 linked list; =2 xsm file; =3 sequential binary file; +* =4 sequential ascii file +* JENTRY =0 the linked list or file is created +* =1 the linked list or file is open for modifications; +* =2 the linked list or file is open in read-only mode +* KENTRY =file unit number; =linked list address otherwise. +* +*Comments: +* The GPTFLU: calling specifications are: +* FLUX\_GPT := GPTFLU: [ FLUX\_GPT ] GPT FLUX0 SYST TRACK :: (gptflu\_data) ; +* where +* FLUX\_GPT : name of the \emph{lcm} object (type L\_FLUX) containing the GPT +* solution. If FLUX\_GPT} appears on the RHS, the solution previously stored +* in FLUX\_GPT} is used to initialize the new iterative process; otherwise, +* a uniform unknown vector is used. +* GPT : name of the \emph{lcm} object (type L\_GPT) containing the +* fixed sources. +* FLUX0 : name of the \emph{lcm} object (type L\_FLUX) containing the +* unperturbed flux used to decontaminate the GPT solution. +* SYST : name of the \emph{lcm} object (type L\_SYSTEM) containing the +* unperturbed system matrices. +* TRACK : name of the \emph{lcm} object (type L\_TRACK) containing the +* \emph{tracking}. +* gptflu\_data}] : structure containing the data to module GPTFLU:} +* +*----------------------------------------------------------------------- +* + USE GANLIB + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY) + CHARACTER HENTRY(NENTRY)*12 + TYPE(C_PTR) KENTRY(NENTRY) +*---- +* LOCAL VARIABLES +*---- + INTEGER IPRINT,NITMA,I,J,IGR,ITYP,LENGT,SIGNA(3),ITY,NSTART + DOUBLE PRECISION DFLOTT + REAL FLOTT + CHARACTER TEXT12*12,CMODUL*12 + LOGICAL LFLU + INTEGER NEL,NUN,NGRP,LL4,ISRC + CHARACTER TITLE*72 +*---- +* STATE-VECTOR VARIABLES +*---- + INTEGER NSTATE + PARAMETER (NSTATE=40) + INTEGER FLUPRM(NSTATE),SYSPRM(NSTATE),TRKPRM(NSTATE), + 1 GPTPRM(NSTATE) +*---- +* Generalized Adjoint calculation +*---- + INTEGER SRCFRM,SRCTO,MAXOUT,ICL1,ICL2,NADI,IMPH,NLF,MAXINR + REAL FKEFF,FKEFF2,EPSOUT,EPSINR,EPSCON(5) + LOGICAL ADJ,REC,RECP + INTEGER, DIMENSION(:), ALLOCATABLE :: MAT,IDL + REAL, DIMENSION(:), ALLOCATABLE :: VOL + REAL, DIMENSION(:,:), ALLOCATABLE :: EVECT,ADECT,EASS,SOUR + TYPE(C_PTR) IPFLUP,IPFLU,IPGPT,IPTRK,IPSYS,JPFLU1,JPFLU2,JPGPT, + 1 KPGPT,JPFLUP,KPFLUP +*---- +* VALIDITY OF OBJECTS +*---- + IF(NENTRY.LT.4) CALL XABORT('GPTFLU: 5 OBJECTS EXPECTED.') + IPFLUP=C_NULL_PTR + IPFLU =C_NULL_PTR + IPGPT =C_NULL_PTR + IPTRK =C_NULL_PTR + IPSYS =C_NULL_PTR + RECP=(JENTRY(1).EQ.1) + LFLU=.FALSE. + DO 2 I=1,NENTRY + TEXT12=HENTRY(I) + IF((IENTRY(I).NE.1).AND.(IENTRY(I).NE.2))CALL XABORT('GPTFLU:' + 1 //' LINKED LIST OR XSM FILE EXPECTED ('//TEXT12//')') + IF((JENTRY(I).EQ.0).AND.(I.EQ.1)) THEN + TEXT12='L_FLUX' + READ(TEXT12,'(3A4)') (SIGNA(J),J=1,3) + CALL LCMPUT(KENTRY(I),'SIGNATURE',3,3,SIGNA) + ELSE + CALL LCMGET(KENTRY(I),'SIGNATURE',SIGNA) + WRITE(TEXT12,'(3A4)') (SIGNA(J),J=1,3) + IF(JENTRY(I).EQ.0) CALL XABORT('GPTFLU:'//TEXT12//' IS ' + 1 //'NOT ON RHS') + ENDIF + IF(TEXT12.EQ.'L_FLUX') THEN + IF(LFLU) THEN + IPFLU=KENTRY(I) + ELSE + IPFLUP=KENTRY(I) + LFLU=.TRUE. + ENDIF + ELSEIF(TEXT12.EQ.'L_SOURCE') THEN + IPGPT=KENTRY(I) + ELSEIF(TEXT12.EQ.'L_TRACK') THEN + IPTRK=KENTRY(I) + ELSEIF(TEXT12.EQ.'L_SYSTEM') THEN + IPSYS=KENTRY(I) + ELSE + CALL XABORT('GPTFLU: NOT GOOD TYPE OF OBJECT') + ENDIF + 2 CONTINUE + IF(.NOT.C_ASSOCIATED(IPGPT)) + 1 CALL XABORT('GPTFLU: MISSING GPT SOURCE OBJECT.') + IF(.NOT.C_ASSOCIATED(IPFLU)) + 1 CALL XABORT('GPTFLU: MISSING FLUX OBJECT.') + IF(.NOT.C_ASSOCIATED(IPSYS)) + 1 CALL XABORT('GPTFLU: MISSING SYSTEM OBJECT.') + IF(.NOT.C_ASSOCIATED(IPTRK)) + 1 CALL XABORT('GPTFLU: MISSING TRACK OBJECT.') +*---- +* VARIABLE INITIALISATION +*---- + CALL LCMGET(IPFLU,'STATE-VECTOR',FLUPRM) + NGRP = FLUPRM(1) + NUN = FLUPRM(2) + MAXINR = FLUPRM(11) + MAXOUT = FLUPRM(12) + CALL LCMGET(IPFLU,'EPS-CONVERGE',EPSCON) + EPSINR=EPSCON(1) + EPSOUT=EPSCON(2) + CALL LCMGET(IPTRK,'STATE-VECTOR',TRKPRM) + NEL = TRKPRM(1) + IF(NUN.NE.TRKPRM(2)) CALL XABORT('GPTFLU: TRACKING AND UNPERTURB' + + //'ED FLUX HAVE DIFFERENT NUMBER OF UNKNOWS') + CALL LCMLEN(IPTRK,'TITLE',LENGT,ITYP) + IF(LENGT.GT.0) THEN + CALL LCMGTC(IPTRK,'TITLE',72,TITLE) + ELSE + TITLE='*** NO TITLE PROVIDED ***' + ENDIF + CALL LCMGTC(IPTRK,'TRACK-TYPE',12,CMODUL) + IF(CMODUL.NE.'TRIVAC') CALL XABORT('GPTFLU: TRIVAC TRACKING EXPE' + + //'CTED.') + LL4 = TRKPRM(11) + NLF = TRKPRM(30) + CALL LCMGET(IPSYS,'STATE-VECTOR',SYSPRM) + IF( SYSPRM(1).NE.NGRP )CALL XABORT('GPTFLU: L_SYSTEM AND L_FLUX' + + //'FOR UNPERTURBED STATE HAVE DIFFERENT NUMBER OF GROUPS') + IF( SYSPRM(2).NE.LL4 )CALL XABORT('GPTFLU: UNPERTURBED SYSTEM A' + + //'ND TRACKING OBJECTS HAVE DIFFERENT NUMBER OF LINEAR ORDER') + ITY = SYSPRM(4) + IF(ITY.EQ.13) LL4=LL4*NLF/2 + CALL LCMGET(IPGPT,'STATE-VECTOR',GPTPRM) + IF( GPTPRM(1).NE.NGRP )CALL XABORT('GPTFLU: L_SOURCE AND L_FLUX ' + + //'HAVE DIFFERENT NUMBER OF GROUPS') + IF( GPTPRM(2).NE.NUN )CALL XABORT('GPTFLU: L_SOURCE AND L_FLUX H' + + //'AVE DIFFERENT NUMBER OF UNKNOWS') +*---- +* READ USER INPUT: +*---- + IPRINT=0 + IMPH=0 + IF(RECP) THEN +* RECOVER EXISTING OPTIONS. + CALL LCMGET(IPFLU,'STATE-VECTOR',FLUPRM) + ICL1=FLUPRM(8) + ICL2=FLUPRM(9) + MAXINR=FLUPRM(11) + MAXOUT=FLUPRM(12) + NADI=FLUPRM(13) + NSTART=FLUPRM(16) + CALL LCMGET(IPFLU,'EPS-CONVERGE',EPSCON) + EPSINR=EPSCON(1) + EPSOUT=EPSCON(2) + ELSE +* DEFAULT OPTIONS. + ICL1=3 + ICL2=3 + NADI=TRKPRM(33) + MAXINR=0 + MAXOUT=200 + NSTART=0 + EPSINR=1.0E-2 + EPSOUT=1.0E-4 + ENDIF + IF(GPTPRM(3).LE.GPTPRM(4)) THEN + ADJ=.FALSE. + ELSE + ADJ=.TRUE. + ENDIF + REC=.FALSE. + 505 CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + 506 IF(ITYP.NE.3) CALL XABORT('GPTFLU: CHARACTER DATA EXPECTED.') + IF(TEXT12.EQ.'EDIT') THEN + CALL REDGET(ITYP,IPRINT,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: *IPRINT* MUST BE INTEGER') + GO TO 505 + ELSEIF((TEXT12.EQ.'VAR1').OR.(TEXT12.EQ.'ACCE')) THEN + CALL REDGET(ITYP,ICL1,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + CALL REDGET(ITYP,ICL2,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + GO TO 505 + ELSEIF(TEXT12.EQ.'GMRES') THEN + CALL REDGET(ITYP,NSTART,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + IF(NSTART.LT.0) CALL XABORT('GPTFLU: POSITIVE VALUE EXPECTED.') + GO TO 505 + ELSEIF(TEXT12.EQ.'IMPLICIT') THEN + ADJ=.TRUE. + GO TO 505 + ELSEIF(TEXT12.EQ.'EXPLICIT') THEN + ADJ=.FALSE. + GO TO 505 + ELSEIF(TEXT12.EQ.'RCVR-LAST') THEN + REC=.TRUE. + GO TO 505 + ELSEIF(TEXT12.EQ.'EXTE') THEN + 507 CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.EQ.1) THEN + MAXOUT=NITMA + ELSE IF(ITYP.EQ.2) THEN + EPSOUT=FLOTT + ELSE + GO TO 506 + ENDIF + GO TO 507 + ELSEIF(TEXT12.EQ.'ADI') THEN + CALL REDGET(ITYP,NADI,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + GO TO 505 + ELSEIF(TEXT12.EQ.'THER') THEN + MAXINR = NGRP*2 + 508 CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.EQ.1) THEN + MAXINR=NITMA + ELSEIF(ITYP.EQ.2) THEN + EPSINR=FLOTT + ELSE + GO TO 506 + ENDIF + GO TO 508 + ELSEIF(TEXT12.EQ.'HIST') THEN + CALL REDGET(ITYP,IMPH,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + GO TO 505 + ENDIF + IF(TEXT12.EQ.'FROM-TO') THEN + CALL REDGET(ITYP,SRCFRM,FLOTT,TEXT12,DFLOTT) + IF((ITYP.EQ.3).AND.(TEXT12.EQ.'ALL')) THEN + SRCFRM=1 + IF(ADJ) THEN + SRCTO=GPTPRM(4) + ELSE + SRCTO=GPTPRM(3) + ENDIF + ELSE + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + CALL REDGET(ITYP,SRCTO,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('GPTFLU: INTEGER DATA EXPECTED.') + IF(ADJ) THEN + IF(SRCTO.GT.GPTPRM(4)) WRITE(6,*) 'THE NUMBER OF THE '// + 1 'SOURCE ',SRCTO,' IS GREATER THAN THE NUMBER OF CONST'// + 2 'RAINTS +1',GPTPRM(4) + ELSE + IF(SRCTO.GT.GPTPRM(3)) WRITE(6,*) 'THE NUMBER OF THE '// + 1 'SOURCE ',SRCTO,' IS GREATER THAN THE NUMBER OF VARIA'// + 2 'BLES',GPTPRM(3) + ENDIF + IF(SRCFRM.GT.SRCTO) CALL XABORT('GPTFLU:SCRFRM .GT. SCRTO') + ENDIF + ELSEIF(TEXT12.EQ.';') THEN + GO TO 1000 + ELSE + WRITE(6,*) 'Your keyword is : ',TEXT12 + CALL XABORT('GPTFLU:"FROM-TO" or ";" EXPECTED') + ENDIF +*---- +* RECOVER TRACKING INFORMATION +*---- + ALLOCATE(MAT(NEL),VOL(NEL),IDL(NEL)) + CALL LCMGET(IPTRK,'MATCOD',MAT) + CALL LCMGET(IPTRK,'VOLUME',VOL) + CALL LCMGET(IPTRK,'KEYFLX',IDL) +*---- +* RECOVER UNPERTURBED K-EFFECTIVE AND FLUXES. +*---- + ALLOCATE(EVECT(NUN,NGRP),ADECT(NUN,NGRP),EASS(NUN,NGRP)) + CALL LCMGET(IPFLU,'K-EFFECTIVE',FKEFF) + JPFLU1=LCMGID(IPFLU,'FLUX') + JPFLU2=LCMGID(IPFLU,'AFLUX') + DO 510 IGR=1,NGRP + CALL LCMGDL(JPFLU1,IGR,EVECT(1,IGR)) + CALL LCMGDL(JPFLU2,IGR,ADECT(1,IGR)) + 510 CONTINUE + ALLOCATE(SOUR(NUN,NGRP)) +*---- +* RECOVER FIXED SOURCE AND SET INITIAL VALUE OF GPFLUX +*---- + IF(ADJ) THEN + JPFLUP=LCMLID(IPFLUP,'ADFLUX',SRCTO) + ELSE + JPFLUP=LCMLID(IPFLUP,'DFLUX',SRCTO) + ENDIF + DO 590 ISRC=SRCFRM,SRCTO + IF(ADJ) THEN + JPGPT=LCMGID(IPGPT,'ASOUR') + ELSE + JPGPT=LCMGID(IPGPT,'DSOUR') + ENDIF + CALL LCMLEL(JPGPT,ISRC,LENGT,ITYP) + IF(LENGT.EQ.0) GO TO 590 + KPGPT=LCMGIL(JPGPT,ISRC) + DO 520 IGR=1,NGRP + CALL LCMGDL(KPGPT,IGR,SOUR(1,IGR)) + 520 CONTINUE + IF(REC.AND.(IMPH.EQ.0)) THEN + CALL LCMLEL(JPFLUP,ISRC,LENGT,ITYP) + IF(LENGT.EQ.0) THEN + WRITE(TEXT12,'(I4,3H-TH)') ISRC + CALL XABORT('GPTFLU: '//TEXT12//' GENERALIZED ADJOINT CANN' + 1 //'OT BE RECOVERED.') + ENDIF + KPFLUP=LCMGIL(JPFLUP,ISRC) + DO 530 IGR=1,NGRP + CALL LCMGDL(KPFLUP,IGR,EASS(1,IGR)) + 530 CONTINUE + ELSE + EASS(:NUN,:NGRP)=1.0 + ENDIF +*---- +* ADJOINT NEUTRON FLUX CALCULATION +*---- + IF(IPRINT.GE.1) WRITE(6,*) 'GPTFLU: ISRC=',ISRC + IF(ADJ) THEN + IF(IPRINT.GE.2) WRITE(6,*) 'implicit' + CALL GPTAFL(IPTRK,IPSYS,IPFLUP,LL4,ITY,NUN,NGRP,ICL1,ICL2, + 1 NSTART,IPRINT,IMPH,TITLE,EPSOUT,MAXINR,EPSINR,NADI,MAXOUT, + 2 FKEFF,EVECT,ADECT,FKEFF2,EASS,SOUR) + ELSE + IF(IPRINT.GE.2) WRITE(6,*) 'explicit' + CALL GPTDFL(IPTRK,IPSYS,IPFLUP,LL4,ITY,NUN,NGRP,ICL1,ICL2, + 1 NSTART,IPRINT,IMPH,TITLE,EPSOUT,MAXINR,EPSINR,NADI,MAXOUT, + 2 FKEFF,EVECT,ADECT,FKEFF2,EASS,SOUR) + ENDIF + CALL LCMPUT(IPFLUP,'K-EFFECTIVE',1,2,FKEFF2) + KPFLUP=LCMLIL(JPFLUP,ISRC,NGRP) + DO 550 IGR=1,NGRP + CALL FLDTRI(IPTRK,NEL,NUN,EASS(1,IGR),MAT,VOL,IDL) + CALL LCMPDL(KPFLUP,IGR,NUN,2,EASS(1,IGR)) + 550 CONTINUE + 590 CONTINUE + DEALLOCATE(EASS,ADECT,EVECT,SOUR,IDL,VOL,MAT) + GO TO 505 +*---- +* END +*---- + 1000 CALL LCMPTC(IPFLUP,'TRACK-TYPE',12,CMODUL) + IF(ADJ) THEN + FLUPRM(3)=1000 + ELSE + FLUPRM(3)=100 + ENDIF + FLUPRM(5)=SRCTO + FLUPRM(6)=1 + FLUPRM(16)=NSTART + IF(IPRINT.GT.0) WRITE(6,2020) (FLUPRM(I),I=1,5),FLUPRM(16) + CALL LCMPUT(IPFLUP,'STATE-VECTOR',NSTATE,1,FLUPRM) + RETURN +* + 2020 FORMAT(/8H OPTIONS/8H -------/ + 1 7H NGRP ,I8,28H (NUMBER OF ENERGY GROUPS)/ + 2 7H NUN ,I8,40H (NUMBER OF UNKNOWNS PER ENERGY GROUP)/ + 3 7H IADJ ,I8,30H (=100/1000: DIRECT/ADJOINT)/ + 4 7H NMOD ,I8,13H (NOT USED)/ + 5 7H SRCTO ,I8,48H (NUMBER OF FIXED-SOURCE EIGENVALUE EQUATIONS)/ + 6 7H NSTART,I8,46H (NUMBER OF GMRES ITERATIONS BEFORE RESTART)) + END |