Initial commit from Polytechnique Montreal

1 files changed, 277 insertions, 0 deletions

diff --git a/Dragon/src/BIVSOU.f b/Dragon/src/BIVSOU.f
new file mode 100644
index 0000000..32c9f2e
--- /dev/null
+++ b/Dragon/src/BIVSOU.f
@@ -0,0 +1,277 @@
+*DECK BIVSOU
+      SUBROUTINE BIVSOU(MAX1,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNKNO,
+     > NGRP,MATCOD,VOL,FUNKNO,SUNKNO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the source for the solution of diffusion or PN equations.
+* BIVAC-specific version.
+*
+*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
+* MAX1    first dimension of FUNKNO and SUNKNO arrays.
+* IG      secondary group.
+* IPTRK   pointer to the tracking LCM object.
+* KPMACR  pointer to the secondary-group related macrolib information.
+* NANIS   maximum cross section Legendre order.
+* NREG    number of regions.
+* NMAT    number of mixtures.
+* NUNKNO  number of unknowns per energy group including spherical
+*         harmonic terms, interface currents and fundamental
+*         currents.
+* NGRP    number of energy groups.
+* MATCOD  mixture indices.
+* VOL     volumes.
+* FUNKNO  fluxes.
+*
+*Parameters: output
+* SUNKNO  sources.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,KPMACR
+      INTEGER MAX1,IG,NANIS,NREG,NMAT,NUNKNO,NGRP,MATCOD(NREG)
+      REAL VOL(NREG),FUNKNO(MAX1,NGRP),SUNKNO(MAX1,NGRP)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40)
+      INTEGER JPAR(NSTATE),IJ1(25),IJ2(25)
+      CHARACTER CAN(0:9)*2
+      LOGICAL CYLIND
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,KN,IDL
+      REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT,XX,DD
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: RR,RS
+*----
+*  DATA STATEMENTS
+*----
+      DATA CAN /'00','01','02','03','04','05','06','07','08','09'/
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IJJ(0:NMAT),NJJ(0:NMAT),IPOS(0:NMAT))
+      ALLOCATE(XSCAT(0:NMAT*NGRP))
+*----
+*  RECOVER BIVAC SPECIFIC PARAMETERS.
+*----
+      CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR)
+      IF(JPAR(1).NE.NREG) CALL XABORT('BIVSOU: INCONSISTENT NREG.')
+      IF(JPAR(2).NE.NUNKNO) CALL XABORT('BIVSOU: INCONSISTENT NUNKNO.')
+      ITYPE=JPAR(6)
+      IELEM=JPAR(8)
+      ICOL=JPAR(9)
+      L4=JPAR(11)
+      LX=JPAR(12)
+      NLF=JPAR(14)
+      ISCAT=JPAR(16)
+      CYLIND=(ITYPE.EQ.3).OR.(ITYPE.EQ.6)
+      IF(ICOL.EQ.4) THEN
+        CALL XABORT('BIVSOU: COLLOCATION NODAL NOT IMPLEMENTED.')
+      ELSE IF((ITYPE.NE.2).AND.(ITYPE.NE.5)) THEN
+        CALL XABORT('BIVSOU: CARTESIAN 1D OR 2D GEOMETRY EXPECTED.')
+      ENDIF
+      CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+      ALLOCATE(XX(NREG),DD(NREG),KN(MAXKN),IDL(NREG))
+      CALL LCMGET(IPTRK,'XX',XX)
+      CALL LCMGET(IPTRK,'DD',DD)
+      CALL LCMGET(IPTRK,'KN',KN)
+      CALL LCMGET(IPTRK,'KEYFLX',IDL)
+*----
+*  RECOVER THE FINITE ELEMENT UNIT STIFFNESS MATRIX.
+*----
+      LL=0
+      IF((NLF.GT.0).OR.(IELEM.LT.0)) THEN
+         CALL LCMSIX(IPTRK,'BIVCOL',1)
+         CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+         ALLOCATE(RR(LC,LC),RS(LC,LC))
+         CALL LCMGET(IPTRK,'R',RR)
+         CALL LCMGET(IPTRK,'RS',RS)
+         CALL LCMSIX(IPTRK,' ',2)
+*----
+*  COMPUTE VECTORS IJ1 AND IJ2
+*----
+         LL=LC*LC
+         DO 10 I=1,LL
+         IJ1(I)=1+MOD(I-1,LC)
+         IJ2(I)=1+(I-IJ1(I))/LC
+   10    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE SOURCE
+*----
+      IF(NLF.EQ.0) THEN
+*----
+*  ++++ DIFFUSION THEORY ++++
+*----
+        CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
+        CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+        CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
+        CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
+        IF(IELEM.GT.0) THEN
+*----
+*  CARTESIAN 2D DUAL (RAVIART-THOMAS) CASE.
+*----
+          NUM1=0
+          DO 30 IR=1,NREG
+          IBM=MATCOD(IR)
+          IF(IBM.LE.0) GO TO 30
+          IF(VOL(IR).EQ.0.0) GO TO 26
+          DO 25 I0=1,IELEM*IELEM
+          IND=KN(NUM1+1)+I0-1
+          JG=IJJ(IBM)
+          DO 20 JND=1,NJJ(IBM)
+          IF(JG.NE.IG) THEN
+            SUNKNO(IND,IG)=SUNKNO(IND,IG)+FUNKNO(IND,JG)*VOL(IR)*
+     >      XSCAT(IPOS(IBM)+JND-1)
+          ENDIF
+          JG=JG-1
+   20     CONTINUE
+   25     CONTINUE
+   26     NUM1=NUM1+5
+   30     CONTINUE
+        ELSE IF(IELEM.LT.0) THEN
+*----
+*  CARTESIAN 2D PRIM (LAGRANGE) CASE.
+*----
+          NUM1=0
+          DO 170 IR=1,NREG
+          IBM=MATCOD(IR)
+          IF(IBM.LE.0) GO TO 170
+          IF(VOL(IR).EQ.0.0) GO TO 160
+          DO 140 I=1,LL
+          IND1=KN(NUM1+I)
+          IF(IND1.EQ.0) GO TO 140
+          I1=IJ1(I)
+          I2=IJ2(I)
+          DO 130 J=1,LL
+          IND2=KN(NUM1+J)
+          IF(IND2.EQ.0) GO TO 130
+          IF(CYLIND) THEN
+            AUXX=(RR(I1,IJ1(J))+RS(I1,IJ1(J))*XX(IR)/DD(IR))*
+     >      RR(I2,IJ2(J))*VOL(IR)
+          ELSE
+            AUXX=RR(I1,IJ1(J))*RR(I2,IJ2(J))*VOL(IR)
+          ENDIF
+          JG=IJJ(IBM)
+          DO 120 JND=1,NJJ(IBM)
+          IF(JG.NE.IG) THEN
+            SUNKNO(IND1,IG)=SUNKNO(IND1,IG)+AUXX*FUNKNO(IND2,JG)*
+     >      XSCAT(IPOS(IBM)+JND-1)
+          ENDIF
+          JG=JG-1
+  120     CONTINUE
+  130     CONTINUE
+  140     CONTINUE
+          ! append the integrated volumic sources
+          JG=IJJ(IBM)
+          DO 150 JND=1,NJJ(IBM)
+          SUNKNO(IDL(IR),IG)=SUNKNO(IDL(IR),IG)+FUNKNO(IDL(IR),JG)*
+     >    VOL(IR)*XSCAT(IPOS(IBM)+JND-1)
+          JG=JG-1
+  150     CONTINUE
+          !
+  160     NUM1=NUM1+LL
+  170     CONTINUE
+        ENDIF
+      ELSE
+*----
+*  ++++ SPN THEORY ++++
+*----
+        DO 330 IL=0,MIN(ABS(ISCAT)-1,NANIS)
+        FACT=REAL(2*IL+1)
+        CALL LCMGET(KPMACR,'NJJS'//CAN(IL),NJJ(1))
+        CALL LCMGET(KPMACR,'IJJS'//CAN(IL),IJJ(1))
+        CALL LCMGET(KPMACR,'IPOS'//CAN(IL),IPOS(1))
+        CALL LCMGET(KPMACR,'SCAT'//CAN(IL),XSCAT(1))
+        NUM1=0
+        DO 320 IR=1,NREG
+        IBM=MATCOD(IR)
+        IF(IBM.LE.0) GO TO 320
+        IF(VOL(IR).EQ.0.0) GO TO 310
+        IF(MOD(IL,2).EQ.0) THEN
+          DO 255 I0=1,IELEM*IELEM
+          IND=(IL/2)*L4+KN(NUM1+1)+I0-1
+          JG=IJJ(IBM)
+          DO 250 JND=1,NJJ(IBM)
+          IF(JG.NE.IG) THEN
+            SUNKNO(IND,IG)=SUNKNO(IND,IG)+FACT*FUNKNO(IND,JG)*
+     >      VOL(IR)*XSCAT(IPOS(IBM)+JND-1)
+          ENDIF
+          JG=JG-1
+  250     CONTINUE
+  255     CONTINUE
+        ELSE
+          DO 305 I0=1,IELEM
+          DO 275 IC=1,2
+          IIC=1+(IC-1)*IELEM
+          IND1=(IL/2)*L4+ABS(KN(NUM1+1+IC))+I0-1
+          S1=REAL(SIGN(1,KN(NUM1+1+IC)))
+          DO 270 JC=1,2
+          JJC=1+(JC-1)*IELEM
+          IND2=(IL/2)*L4+ABS(KN(NUM1+1+JC))+I0-1
+          IF((KN(NUM1+1+IC).NE.0).AND.(KN(NUM1+1+JC).NE.0)) THEN
+            S2=REAL(SIGN(1,KN(NUM1+1+JC)))
+            AUXX=S1*S2*FACT*RR(IIC,JJC)*VOL(IR)
+            JG=IJJ(IBM)
+            DO 260 JND=1,NJJ(IBM)
+            IF(JG.NE.IG) THEN
+              SUNKNO(IND1,IG)=SUNKNO(IND1,IG)-AUXX*FUNKNO(IND2,JG)*
+     1        XSCAT(IPOS(IBM)+JND-1)
+            ENDIF
+            JG=JG-1
+  260       CONTINUE
+          ENDIF
+  270     CONTINUE
+  275     CONTINUE
+          DO 300 IC=3,4
+          IIC=1+(IC-3)*IELEM
+          IND1=(IL/2)*L4+ABS(KN(NUM1+1+IC))+I0-1
+          S1=REAL(SIGN(1,KN(NUM1+1+IC)))
+          DO 290 JC=3,4
+          JJC=1+(JC-3)*IELEM
+          IND2=(IL/2)*L4+ABS(KN(NUM1+1+JC))+I0-1
+          IF((KN(NUM1+1+IC).NE.0).AND.(KN(NUM1+1+JC).NE.0)) THEN
+            S2=REAL(SIGN(1,KN(NUM1+1+JC)))
+            AUXX=S1*S2*FACT*RR(IIC,JJC)*VOL(IR)
+            JG=IJJ(IBM)
+            DO 280 JND=1,NJJ(IBM)
+            IF(JG.NE.IG) THEN
+              SUNKNO(IND1,IG)=SUNKNO(IND1,IG)-AUXX*FUNKNO(IND2,JG)*
+     1        XSCAT(IPOS(IBM)+JND-1)
+            ENDIF
+            JG=JG-1
+  280       CONTINUE
+          ENDIF
+  290     CONTINUE
+  300     CONTINUE
+  305     CONTINUE
+        ENDIF
+  310   NUM1=NUM1+5
+  320   CONTINUE
+  330   CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      IF((NLF.GT.0).OR.(IELEM.LT.0)) DEALLOCATE(RS,RR)
+      DEALLOCATE(IDL,KN,DD,XX)
+      DEALLOCATE(XSCAT)
+      DEALLOCATE(IPOS,NJJ,IJJ)
+      RETURN
+      END