Initial commit from Polytechnique Montreal

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