Initial commit from Polytechnique Montreal

1 files changed, 300 insertions, 0 deletions

diff --git a/Trivac/src/BIVSPS.f b/Trivac/src/BIVSPS.f
new file mode 100755
index 0000000..d9e1ff4
--- /dev/null
+++ b/Trivac/src/BIVSPS.f
@@ -0,0 +1,300 @@
+*DECK BIVSPS
+      SUBROUTINE BIVSPS(IPTRK,IPMACR,IPSYS,IMPX,NGRP,NEL,NLF,NANI,NBFIS,
+     1 NALBP,LDIFF,MAT,VOL,NBMIX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Recover the cross-section data in LCM object with pointer IPMACR,
+* compute and store the corresponding system matrices for a simplified
+* PN approximation.
+*
+*Copyright:
+* Copyright (C) 2004 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
+* IPTRK   L_TRACK pointer to the bivac tracking information.
+* IPMACR  L_MACROLIB pointer to the cross sections.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IMPX    print parameter (equal to zero for no print).
+* NGRP    number of energy groups.
+* NEL     total number of finite elements.
+* NLF     number of Legendre orders for the flux (even number).
+* NANI    number of Legendre orders for the scattering cross sections.
+* NBFIS   number of fissionable isotopes.
+* NALBP   number of physical albedos per energy group.
+* LDIFF   flag set to .true. to use 1/3D as 'NTOT1' cross sections.
+* MAT     index-number of the mixture type assigned to each volume.
+* VOL     volumes.
+* NBMIX   total number of material mixtures in the macrolib.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPMACR,IPSYS
+      INTEGER IMPX,NGRP,NEL,NLF,NANI,NBFIS,NALBP,MAT(NEL),NBMIX
+      REAL VOL(NEL)
+      LOGICAL LDIFF
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER TEXT12*12,CM*2,HSMG*131
+      LOGICAL LFIS
+      TYPE(C_PTR) JPMACR,KPMACR
+      INTEGER, DIMENSION(:), ALLOCATABLE :: IJJ,NJJ,IPOS,IND
+      REAL, DIMENSION(:), ALLOCATABLE :: WORK
+      REAL, DIMENSION(:,:), ALLOCATABLE :: ALBP,GAMMA,SGD,ZUFIS
+      REAL, DIMENSION(:,:,:), ALLOCATABLE :: CHI,RCAT,RCATI
+*
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IJJ(NBMIX),NJJ(NBMIX),IPOS(NBMIX),IND(NGRP))
+      ALLOCATE(GAMMA(NALBP,NGRP),SGD(NBMIX,2*NLF),WORK(NBMIX*NGRP),
+     1 CHI(NBMIX,NBFIS,NGRP),ZUFIS(NBMIX,NBFIS),RCAT(NGRP,NGRP,NBMIX),
+     2 RCATI(NGRP,NGRP,NBMIX))
+*----
+*  PROCESS PHYSICAL ALBEDO INFORMATION AND CALCULATION OF
+*  MULTIGROUP ALBEDO FUNCTIONS (RAVIART-THOMAS CASE).
+*----
+      IF(NALBP.GT.0) THEN
+         ALLOCATE(ALBP(NALBP,NGRP))
+         CALL LCMGET(IPMACR,'ALBEDO',ALBP)
+         DO IGR=1,NGRP
+            GAMMA(:NALBP,IGR)=0.0
+            DO IALB=1,NALBP
+               IF(ALBP(IALB,IGR).NE.1.0) THEN
+                  GAMMA(IALB,IGR)=1.0/ALB(ALBP(IALB,IGR))
+               ELSE
+                  GAMMA(IALB,IGR)=1.0E20
+               ENDIF
+            ENDDO
+            WRITE(TEXT12,'(9HALBEDO-FU,I3.3)') IGR
+            CALL LCMPUT(IPSYS,TEXT12,NALBP,2,GAMMA(1,IGR))
+         ENDDO
+         DEALLOCATE(ALBP)
+      ENDIF
+*----
+*  PROCESS MACROLIB INFORMATION FOR VARIOUS LEGENDRE ORDERS.
+*----
+      IF(NLF.EQ.0) CALL XABORT('BIVSPS: SPN APPROXIMATION REQUESTED.')
+      JPMACR=LCMGID(IPMACR,'GROUP')
+      DO 112 IL=1,NLF
+      WRITE(CM,'(I2.2)') IL-1
+      RCAT(:NGRP,:NGRP,:NBMIX)=0.0
+      DO 50 IGR=1,NGRP
+*     PROCESS SECONDARY GROUP IGR.
+      KPMACR=LCMGIL(JPMACR,IGR)
+      SGD(:NBMIX,1)=0.0
+      CALL LCMLEN(KPMACR,'SIGW'//CM,LENGT,ITYLCM)
+      IF((LENGT.GT.0).AND.(IL.LE.NANI)) THEN
+         IF(LENGT.GT.NBMIX) CALL XABORT('BIVSPS: INVALID LENGTH FOR'
+     1   //' SIGW'//CM//' CROSS SECTIONS.')
+         CALL LCMGET(KPMACR,'SIGW'//CM,SGD(1,1))
+      ENDIF
+      WRITE(TEXT12,'(4HNTOT,I1)') MIN(IL-1,9)
+      CALL LCMLEN(KPMACR,TEXT12,LENGT,ITYLCM)
+      CALL LCMLEN(KPMACR,'NTOT1',LENGT1,ITYLCM)
+      IF((IL.EQ.1).AND.(LENGT.NE.NBMIX)) CALL XABORT('BIVSPS: NO NTOT0'
+     1 //' CROSS SECTIONS.')
+      IF(MOD(IL-1,2).EQ.0) THEN
+*        macroscopic total cross section in even-parity equations.
+         IF(LENGT.EQ.NBMIX) THEN
+            CALL LCMGET(KPMACR,TEXT12,SGD(1,2))
+         ELSE
+            CALL LCMGET(KPMACR,'NTOT0',SGD(1,2))
+         ENDIF
+      ELSE
+*        macroscopic total cross section in odd-parity equations.
+         IF(LDIFF) THEN
+            CALL LCMLEN(KPMACR,'DIFF',LENGT,ITYLCM)
+            IF(LENGT.EQ.0) CALL XABORT('BIVSPS: DIFFUSION COEFFICIENTS'
+     1      //' EXPECTED IN THE MACROLIB.')
+            IF(LENGT.GT.NBMIX) CALL XABORT('BIVSPS: INVALID LENGTH FOR'
+     1      //' DIFFUSION COEFFICIENTS.')
+            CALL LCMGET(KPMACR,'DIFF',SGD(1,2))
+            DO 5 IBM=1,NBMIX
+            SGD(IBM,2)=1.0/(3.0*SGD(IBM,2))
+    5       CONTINUE
+         ELSE IF(LENGT.EQ.NBMIX) THEN
+            CALL LCMGET(KPMACR,TEXT12,SGD(1,2))
+         ELSE IF(LENGT1.EQ.NBMIX) THEN
+            CALL LCMGET(KPMACR,'NTOT1',SGD(1,2))
+         ELSE
+            CALL LCMGET(KPMACR,'NTOT0',SGD(1,2))
+         ENDIF
+      ENDIF
+      DO 10 IBM=1,NBMIX
+      IF((MOD(IL-1,2).NE.0).AND.LDIFF) THEN
+        RCAT(IGR,IGR,IBM)=SGD(IBM,2)
+      ELSE
+        IF(SGD(IBM,1).GT.SGD(IBM,2)) THEN
+          WRITE(HSMG,'(28HBIVSPS: NEGATIVE XS IN GROUP,I5)') IGR
+          CALL XABORT(HSMG)
+        ENDIF
+        RCAT(IGR,IGR,IBM)=SGD(IBM,2)-SGD(IBM,1)
+      ENDIF
+   10 CONTINUE
+      IF((MOD(IL-1,2).NE.0).AND.LDIFF) GO TO 50
+      CALL LCMLEN(KPMACR,'NJJS'//CM,LENGT,ITYLCM)
+      IF(LENGT.GT.NBMIX) CALL XABORT('BIVSPS: INVALID LENGTH FOR NJJS'
+     1 //CM//' INFORMATION.')
+      IF((LENGT.GT.0).AND.(IL.LE.NANI)) THEN
+         CALL LCMGET(KPMACR,'NJJS'//CM,NJJ)
+         CALL LCMGET(KPMACR,'IJJS'//CM,IJJ)
+         IGMIN=IGR
+         IGMAX=IGR
+         DO 20 IBM=1,NBMIX
+         IGMIN=MIN(IGMIN,IJJ(IBM)-NJJ(IBM)+1)
+         IGMAX=MAX(IGMAX,IJJ(IBM))
+   20    CONTINUE
+         CALL LCMGET(KPMACR,'IPOS'//CM,IPOS)
+         CALL LCMGET(KPMACR,'SCAT'//CM,WORK)
+         DO 40 JGR=IGMAX,IGMIN,-1
+         IF(JGR.EQ.IGR) GO TO 40
+         DO 30 IBM=1,NBMIX
+         IF((JGR.GT.IJJ(IBM)-NJJ(IBM)).AND.(JGR.LE.IJJ(IBM))) THEN
+            RCAT(IGR,JGR,IBM)=-WORK(IPOS(IBM)+IJJ(IBM)-JGR)
+         ENDIF
+   30    CONTINUE
+   40    CONTINUE
+      ENDIF
+   50 CONTINUE
+*----
+*  INVERSION OF THE REMOVAL MATRIX FOR CASES WITH IELEM > 1.
+*----
+      DO 70 IBM=1,NBMIX
+      DO 65 JGR=1,NGRP
+      DO 60 IGR=1,NGRP
+      RCATI(IGR,JGR,IBM)=RCAT(IGR,JGR,IBM)
+   60 CONTINUE
+   65 CONTINUE
+      CALL ALINV(NGRP,RCATI(1,1,IBM),NGRP,IER,IND)
+      IF(IER.NE.0) CALL XABORT('BIVSPS: SINGULAR MATRIX.')
+   70 CONTINUE
+*
+      DO 111 IGR=1,NGRP
+      IGMIN=IGR
+      IGMAX=IGR
+      DO 85 IBM=1,NBMIX
+      DO 80 JGR=1,NGRP
+      IF((RCAT(IGR,JGR,IBM).NE.0.0).OR.(RCATI(IGR,JGR,IBM).NE.0.0)) THEN
+         IGMIN=MIN(IGMIN,JGR)
+         IGMAX=MAX(IGMAX,JGR)
+      ENDIF
+   80 CONTINUE
+   85 CONTINUE
+      DO 110 JGR=IGMIN,IGMAX
+      DO 90 IBM=1,NBMIX
+      WORK(IBM)=RCAT(IGR,JGR,IBM)
+   90 CONTINUE
+      WRITE(TEXT12,'(4HSCAR,A2,2I3.3)') CM,IGR,JGR
+      CALL LCMPUT(IPSYS,TEXT12,NBMIX,2,WORK)
+      DO 100 IBM=1,NBMIX
+      WORK(IBM)=RCATI(IGR,JGR,IBM)
+  100 CONTINUE
+      WRITE(TEXT12,'(4HSCAI,A2,2I3.3)') CM,IGR,JGR
+      CALL LCMPUT(IPSYS,TEXT12,NBMIX,2,WORK)
+  110 CONTINUE
+  111 CONTINUE
+  112 CONTINUE
+*----
+*  COMPUTE AND FACTORIZE THE DIAGONAL SYSTEM MATRICES.
+*----
+      DO 162 IGR=1,NGRP
+      DO 140 IL=1,NLF
+      WRITE(TEXT12,'(4HSCAR,I2.2,2I3.3)') IL-1,IGR,IGR
+      CALL LCMGET(IPSYS,TEXT12,SGD(1,IL))
+      WRITE(TEXT12,'(4HSCAI,I2.2,2I3.3)') IL-1,IGR,IGR
+      CALL LCMGET(IPSYS,TEXT12,SGD(1,NLF+IL))
+  140 CONTINUE
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL BIVASM(TEXT12,0,IPTRK,IPSYS,IMPX,NBMIX,NEL,NLF,2*NLF,NALBP,
+     1 MAT,VOL,GAMMA,SGD)
+*----
+*  PUT A FLAG IN IPSYS TO IDENTIFY NON-ZERO SCATTERING TERMS.
+*----
+      DO 161 IL=1,NLF
+      DO 160 JGR=1,NGRP
+      WRITE(TEXT12,'(4HSCAR,I2.2,2I3.3)') IL-1,IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,LENGT,ITYLCM)
+      IF(LENGT.EQ.NBMIX) THEN
+         WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+         CALL LCMPUT(IPSYS,TEXT12,1,2,0.0)
+      ENDIF
+  160 CONTINUE
+  161 CONTINUE
+  162 CONTINUE
+*----
+*  PROCESS FISSION SPECTRUM TERMS.
+*----
+      KPMACR=LCMGIL(JPMACR,1)
+      CALL LCMLEN(KPMACR,'CHI',LENGT,ITYLCM)
+      IF(LENGT.GT.0) THEN
+         IF(LENGT.NE.NBMIX*NBFIS) CALL XABORT('BIVSPS: INVALID LENGTH '
+     1   //'FOR CHI INFORMATION.')
+         DO 170 IGR=1,NGRP
+         KPMACR=LCMGIL(JPMACR,IGR)
+         CALL LCMGET(KPMACR,'CHI',CHI(1,1,IGR))
+  170    CONTINUE
+      ELSE
+         DO 182 IBM=1,NBMIX
+         DO 181 IFISS=1,NBFIS
+         CHI(IBM,IFISS,1)=1.0
+         DO 180 IGR=2,NGRP
+         CHI(IBM,IFISS,IGR)=0.0
+  180    CONTINUE
+  181    CONTINUE
+  182    CONTINUE
+      ENDIF
+*----
+*  PROCESS FISSION NUSIGF TERMS.
+*----
+      DO 220 IGR=1,NGRP
+*     PROCESS SECONDARY GROUP IGR.
+      LFIS=.FALSE.
+      DO 195 IBM=1,NBMIX
+      DO 190 IFISS=1,NBFIS
+      LFIS=LFIS.OR.(CHI(IBM,IFISS,IGR).NE.0.0)
+  190 CONTINUE
+  195 CONTINUE
+      IF(LFIS) THEN
+         DO 210 JGR=1,NGRP
+         KPMACR=LCMGIL(JPMACR,JGR)
+         CALL LCMLEN(KPMACR,'NUSIGF',LENGT,ITYLCM)
+         IF(LENGT.NE.NBMIX*NBFIS) CALL XABORT('BIVSPS: INVALID LENGTH '
+     1   //'FOR NUSIGF INFORMATION.')
+         IF(LENGT.GT.0) THEN
+            CALL LCMGET(KPMACR,'NUSIGF',ZUFIS)
+            SGD(:NBMIX,:2*NLF)=0.0
+            DO 205 IBM=1,NBMIX
+            DO 200 IFISS=1,NBFIS
+            SGD(IBM,1)=SGD(IBM,1)+CHI(IBM,IFISS,IGR)*ZUFIS(IBM,IFISS)
+  200       CONTINUE
+  205       CONTINUE
+            WRITE(TEXT12,'(4HFISS,2I3.3)') IGR,JGR
+            CALL LCMPUT(IPSYS,TEXT12,NBMIX,2,SGD(1,1))
+            WRITE (TEXT12,'(1HB,2I3.3)') IGR,JGR
+            CALL BIVASM(TEXT12,1,IPTRK,IPSYS,IMPX,NBMIX,NEL,2,4,NALBP,
+     1      MAT,VOL,GAMMA,SGD)
+         ENDIF
+  210    CONTINUE
+      ENDIF
+  220 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IJJ,NJJ,IPOS,IND)
+      DEALLOCATE(GAMMA,SGD,WORK,CHI,ZUFIS,RCAT,RCATI)
+      RETURN
+      END