Initial commit from Polytechnique Montreal

1 files changed, 303 insertions, 0 deletions

diff --git a/Trivac/src/TRIVAA.f b/Trivac/src/TRIVAA.f
new file mode 100755
index 0000000..acee907
--- /dev/null
+++ b/Trivac/src/TRIVAA.f
@@ -0,0 +1,303 @@
+*DECK TRIVAA
+      SUBROUTINE TRIVAA(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* TRIVAC type (3-D and ADI) system matrix assembly operator.
+*
+*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
+*
+*Parameters: input/output
+* NENTRY  number of LCM objects or files used by the operator.
+* HENTRY  name of each LCM object or file:
+*         HENTRY(1): create or modification type(L_SYSTEM);
+*         HENTRY(2): read-only type(L_MACROLIB) (unperturbed);
+*         HENTRY(3): read-only type(L_TRACK);
+*         HENTRY(4): optional read-only type(L_MACROLIB) (perturbed).
+* IENTRY  type of each LCM object or file:
+*         =1 LCM memory object; =2 XSM file; =3 sequential binary file;
+*         =4 sequential ascii file.
+* JENTRY  access of each LCM object or file:
+*         =0 the LCM object or file is created;
+*         =1 the LCM object or file is open for modifications;
+*         =2 the LCM object or file is open in read-only mode.
+* KENTRY  LCM object address or file unit number.
+*
+*Comments:
+* The TRIVAA: calling specifications are:
+* SYST := TRIVAA: [ SYST ] MACRO  TRACK [ DMACRO ] :: (trivaa\_data) ;
+* where
+*   SYST   : name of the \emph{lcm} object (type L\_SYSTEM) containing the 
+*     system matrices. If SYST appears on the RHS, the system matrices 
+*     previously stored in SYST are kept.
+*   MACRO  : name of the \emph{lcm} object (type L\_MACROLIB) containing the 
+*     macroscopic cross sections and diffusion coefficients.
+*   TRACK  : name of the \emph{lcm} object (type L\_TRIVAC) containing the 
+*     TRIVAC \emph{tracking}.
+*   DMACRO : name of the \emph{lcm} object  (type L\_MACROLIB) containing 
+*     derivatives or perturbations of the macroscopic cross sections and 
+*     diffusion coefficients. If DMACRO is given, only the derivatives or 
+*     perturbations of the system matrices are computed.
+*   trivaa\_data : structure containing the data to module TRIVAA:
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY)
+      CHARACTER HENTRY(NENTRY)*12
+      TYPE(C_PTR) KENTRY(NENTRY)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40)
+      CHARACTER TEXT4*4,TEXT11*12,TEXT12*12,HSMG*131,TITLE*72,CNAM*12
+      DOUBLE PRECISION DFLOTT
+      INTEGER IGP(NSTATE),IPAR(NSTATE),ITR(NSTATE)
+      LOGICAL LDIFF
+      TYPE(C_PTR) IPSYS,JPSYS,KPSYS,IPMACR,JPMACR,KPMACR,IPTRK,IPMACP
+      INTEGER, DIMENSION(:), ALLOCATABLE :: MAT
+      REAL, DIMENSION(:), ALLOCATABLE :: VOL,UN,VII
+*----
+*  PARAMETER VALIDATION.
+*----
+      IF(NENTRY.LE.2) CALL XABORT('TRIVAA: THREE PARAMETERS EXPECTED.')
+      IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2)) CALL XABORT('TRIVAA: L'
+     1 //'CM OBJECT EXPECTED AT LHS.')
+      IF((JENTRY(1).NE.0).AND.(JENTRY(1).NE.1)) CALL XABORT('TRIVAA: E'
+     1 //'NTRY IN CREATE OR MODIFICATION MODE EXPECTED.')
+      IF((JENTRY(3).NE.2).OR.((IENTRY(3).NE.1).AND.(IENTRY(3).NE.2)))
+     1 CALL XABORT('TRIVAA: LCM OBJECT IN READ-ONLY MODE EXPECTED AT F'
+     2 //'IRST RHS.')
+      IF((JENTRY(2).NE.2).OR.((IENTRY(2).NE.1).AND.(IENTRY(2).NE.2)))
+     1 CALL XABORT('TRIVAA: LCM OBJECT IN READ-ONLY MODE EXPECTED AT S'
+     2 //'ECOND RHS.')
+      CALL LCMGTC(KENTRY(3),'SIGNATURE',12,TEXT11)
+      IF(TEXT11.NE.'L_TRACK') THEN
+         TEXT12=HENTRY(3)
+         CALL XABORT('TRIVAA: SIGNATURE OF '//TEXT12//' IS '//TEXT11//
+     1   '. L_TRACK EXPECTED.')
+      ENDIF
+      CALL LCMGTC(KENTRY(3),'TRACK-TYPE',12,TEXT11)
+      IF(TEXT11.NE.'TRIVAC') THEN
+         TEXT12=HENTRY(3)
+         CALL XABORT('TRIVAA: TRACK-TYPE OF '//TEXT12//' IS '//TEXT11
+     1   //'. TRIVAC EXPECTED.')
+      ENDIF
+      TEXT11='L_SYSTEM'
+      IPSYS=KENTRY(1)
+      CALL LCMPTC(IPSYS,'SIGNATURE',12,TEXT11)
+      IPMACR=KENTRY(2)
+      IPTRK=KENTRY(3)
+      TEXT12=HENTRY(2)
+      CALL LCMPTC(IPSYS,'LINK.MACRO',12,TEXT12)
+      TEXT12=HENTRY(3)
+      CALL LCMPTC(IPSYS,'LINK.TRACK',12,TEXT12)
+*----
+*  RECOVER GENERAL TRACKING INFORMATION.
+*----
+      CALL LCMGET(IPTRK,'STATE-VECTOR',IGP)
+      NEL=IGP(1)
+      NLF=IGP(30)
+      ISCAT=IGP(32)
+      LDIFF=(ISCAT.LT.0)
+      ISCAT=ABS(ISCAT)
+      IF((NLF.NE.0).AND.(IGP(31).NE.1)) CALL XABORT('TRIVAA: ONLY SPN '
+     1 //'DISCRETIZATIONS ARE ALLOWED.')
+      ITY=2
+      IF(IGP(12).EQ.2) ITY=3
+      ALLOCATE(MAT(NEL),VOL(NEL))
+      CALL LCMGET(IPTRK,'MATCOD',MAT)
+      CALL LCMGET(IPTRK,'VOLUME',VOL)
+      CALL LCMLEN(IPTRK,'TITLE',LENGT,ITYLCM)
+      IF(LENGT.GT.0) THEN
+         CALL LCMGTC(IPTRK,'TITLE',72,TITLE)
+      ELSE
+         TITLE='*** NO TITLE PROVIDED ***'
+      ENDIF
+*----
+*  RECOVER MACROLIB PARAMETERS.
+*----
+      CALL LCMGTC(IPMACR,'SIGNATURE',12,TEXT11)
+      IF(TEXT11.NE.'L_MACROLIB') THEN
+         TEXT12=HENTRY(2)
+         CALL XABORT('TRIVAA: SIGNATURE OF '//TEXT12//' IS '//TEXT11//
+     1   '. L_MACROLIB EXPECTED.')
+      ENDIF
+      CALL LCMGET(IPMACR,'STATE-VECTOR',IPAR)
+      NGRP=IPAR(1)
+      NBMIX=IPAR(2)
+      NANI=IPAR(3)
+      NBFIS=IPAR(4)
+      NALBP=IPAR(8)
+      IF(IGP(4).GT.NBMIX) THEN
+         WRITE(HSMG,'(46HTRIVAA: THE NUMBER OF MIXTURES IN THE TRACKING,
+     1   2H (,I5,51H) IS GREATER THAN THE NUMBER OF MIXTURES IN THE MAC,
+     2   7HROLIB (,I5,2H).)') IGP(4),NBMIX
+         CALL XABORT(HSMG)
+      ENDIF
+*
+      IMPX=1
+      IASM=0
+      IPR=0
+      IUNIT=0
+      IOVEL=0
+      NSTEP=0
+   10 CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT)
+      IF(INDIC.EQ.10) GO TO 30
+      IF(INDIC.NE.3) CALL XABORT('TRIVAA: CHARACTER DATA EXPECTED.')
+      IF(TEXT4.EQ.'EDIT') THEN
+         CALL REDGET(INDIC,IMPX,FLOTT,TEXT4,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('TRIVAA: INTEGER DATA EXPECTED(1).')
+      ELSE IF(TEXT4.EQ.'SKIP') THEN
+*        OPTION TO SKIP THE SYSTEM MATRIX ASSEMBLY (DO NOT SKIP THE
+*        LDLT FACTORIZATION).
+         IASM=1
+      ELSE IF(TEXT4.EQ.'DERI') THEN
+         IPR=1
+         WRITE(6,'(/43H TRIVAA: USE DERIVATIVE OF SYSTEM MATRICES.)')
+      ELSE IF(TEXT4.EQ.'PERT') THEN
+         IPR=2
+         WRITE(6,'(/41H TRIVAA: PERTURBATION OF SYSTEM MATRICES.)')
+      ELSE IF(TEXT4.EQ.'UNIT') THEN
+*       COMPUTE THE UNITARY WEIGHTING MATRIX.
+        IUNIT=1
+        ALLOCATE(UN(NBMIX))
+        UN(:NBMIX)=1.0
+        CALL TRIDIG('RM',IPTRK,IPSYS,IMPX,NBMIX,NEL,0,MAT,VOL,UN)
+        DEALLOCATE(UN)
+      ELSE IF(TEXT4.EQ.'OVEL') THEN
+*        COMPUTE THE RECIPROCAL NEUTRON VELOCITIES MATRIX.
+         IOVEL=1
+         JPMACR=LCMGID(IPMACR,'GROUP')
+         ALLOCATE(VII(NBMIX))
+         DO 25 IGR=1,NGRP
+         KPMACR=LCMGIL(JPMACR,IGR)
+         CALL LCMLEN(KPMACR,'OVERV',LENGT,ITYLCM)
+         IF(LENGT.EQ.0) THEN
+            CALL XABORT('TRIVAA: NO ''VELOCITY'' INFORMATION.')
+         ELSE IF(LENGT.GT.NBMIX) THEN
+            CALL XABORT('TRIVAA: INVALID LENGTH FOR ''VELOCITY'' IN'
+     1      //'FORMATION.')
+         ENDIF
+         CALL LCMGET(KPMACR,'OVERV',VII)
+         WRITE (CNAM,'(1HV,2I3.3)') IGR,IGR
+         CALL TRIDIG(CNAM,IPTRK,IPSYS,IMPX,NBMIX,NEL,0,MAT,VOL,VII)
+   25    CONTINUE
+         DEALLOCATE(VII)
+      ELSE IF(TEXT4.EQ.';') THEN
+         GO TO 30
+      ELSE
+         CALL XABORT('TRIVAA: '//TEXT4//' IS AN INVALID KEY WORD.')
+      ENDIF
+      GO TO 10
+*----
+*  2-MACROLIBS PERTURBATION CALCULATION.
+*----
+   30 IF(IPR.GT.0) THEN
+         IF(NENTRY.LE.3) CALL XABORT('TRIVAA: 4 PARAMETERS EXPECTED WIT'
+     1   //'H DERI OR PERT OPTIONS.')
+         IF((JENTRY(4).NE.2).OR.((IENTRY(4).NE.1).AND.(IENTRY(4).NE.2)))
+     1   CALL XABORT('TRIVAA: LINKED LIST OR XSM FILE IN READ-ONLY MODE'
+     2   //' EXPECTED AT THIRD RHS.')
+         IPMACP=KENTRY(4)
+         CALL LCMGTC(IPMACP,'SIGNATURE',12,TEXT11)
+         IF(TEXT11.NE.'L_MACROLIB') THEN
+            TEXT12=HENTRY(4)
+            CALL XABORT('TRIVAA: SIGNATURE OF '//TEXT12//' IS '
+     1      //TEXT11//'. L_MACROLIB EXPECTED.')
+         ENDIF
+         CALL LCMGET(IPMACP,'STATE-VECTOR',IPAR)
+         NSTEP=IPAR(11)
+         IF((IPAR(1).NE.NGRP).OR.(IPAR(2).GT.NBMIX)) THEN
+            WRITE(HSMG,'(43HTRIVAA: INCONSISTENT PERTURBATION MACROLIB ,
+     1      1H'',A12,8H''. NGRP=,2I5,7H NBMIX=,2I9)') HENTRY(4),IPAR(1),
+     2      NGRP,IPAR(2),NBMIX
+            CALL XABORT(HSMG)
+         ENDIF
+      ENDIF
+*----
+*  SET THE STATE VECTOR FOR THE L_SYSTEM OBJECT
+*----
+      ITR(:NSTATE)=0
+      ITR(1)=NGRP
+      ITR(2)=IGP(11)
+      ITR(3)=0
+      ITR(4)=ITY
+      IF((NLF.GT.0).AND.(ITY.GE.3)) ITR(4)=10+ITR(4)
+      IF(IUNIT.EQ.1) ITR(5)=1
+      ITR(6)=NSTEP
+      ITR(7)=NBMIX
+      NAN=MIN(ISCAT,NANI)
+      ITR(8)=NLF
+      ITR(9)=IPR
+      CALL LCMPUT(IPSYS,'STATE-VECTOR',NSTATE,1,ITR)
+*----
+*  SYSTEM MATRIX ASSEMBLY.
+*----
+      IF((IASM.EQ.0).AND.(IPR.EQ.0)) THEN
+         IF(NLF.EQ.0) THEN
+*           DIFFUSION THEORY.
+            CALL TRISYS(IPTRK,IPMACR,IPMACR,IPSYS,IMPX,NGRP,NEL,NBFIS,
+     1      NALBP,IPR,MAT,VOL,NBMIX)
+         ELSE
+*           SIMPLIFIED PN THEORY.
+            CALL TRISPS(IPTRK,IPMACR,IPMACR,IPSYS,IMPX,NGRP,NEL,NLF,
+     1      NAN,NBFIS,NALBP,LDIFF,IPR,MAT,VOL,NBMIX)
+         ENDIF
+      ELSE IF((IASM.EQ.1).AND.(IPR.EQ.0)) THEN
+*        PERFORM FACTORIZATION WITHOUT ASSEMBLY.
+         DO 40 I=1,NGRP
+         WRITE(TEXT11,'(1HA,2I3.3)') I,I
+         CALL MTLDLF(TEXT11,IPTRK,IPSYS,ITY,IMPX)
+   40    CONTINUE
+      ELSE IF((IPR.GT.0).AND.(NSTEP.EQ.0)) THEN
+*        ASSEMBLY OF PERTURBED SYSTEM MATRICES (NO STEP DIRECTORIES).
+         IF(NLF.EQ.0) THEN
+*           DIFFUSION THEORY.
+            CALL TRISYS(IPTRK,IPMACR,IPMACP,IPSYS,IMPX,NGRP,NEL,NBFIS,
+     1      NALBP,IPR,MAT,VOL,NBMIX)
+         ELSE
+*           SIMPLIFIED PN THEORY.
+            CALL TRISPS(IPTRK,IPMACR,IPMACP,IPSYS,IMPX,NGRP,NEL,NLF,
+     1      NAN,NBFIS,NALBP,LDIFF,IPR,MAT,VOL,NBMIX)
+         ENDIF
+      ELSE IF(NSTEP.GT.0) THEN
+*        ASSEMBLY OF PERTURBED SYSTEM MATRICES (WITH STEP DIRECTORIES).
+         JPMACR=LCMGID(IPMACP,'STEP')
+         JPSYS=LCMLID(IPSYS,'STEP',NSTEP)
+         DO 50 ISTEP=1,NSTEP
+         KPMACR=LCMGIL(JPMACR,ISTEP)
+         KPSYS=LCMDIL(JPSYS,ISTEP)
+         CALL LCMPUT(KPSYS,'STATE-VECTOR',NSTATE,1,ITR)
+         IF(NLF.EQ.0) THEN
+*           DIFFUSION THEORY.
+            CALL TRISYS(IPTRK,IPMACR,KPMACR,KPSYS,IMPX,NGRP,NEL,NBFIS,
+     1      NALBP,IPR,MAT,VOL,NBMIX)
+         ELSE
+*           SIMPLIFIED PN THEORY.
+            CALL TRISPS(IPTRK,IPMACR,KPMACR,KPSYS,IMPX,NGRP,NEL,NLF,
+     1      NAN,NBFIS,NALBP,LDIFF,IPR,MAT,VOL,NBMIX)
+         ENDIF
+   50    CONTINUE
+      ELSE
+         CALL XABORT('TRIVAA: INVALID REQUEST.')
+      ENDIF
+*
+      IF(IMPX.GE.3) CALL LCMLIB(IPSYS)
+*----
+*  RELEASE GENERAL TRACKING INFORMATION.
+*----
+      DEALLOCATE(VOL,MAT)
+      RETURN
+      END