summaryrefslogtreecommitdiff
path: root/Dragon/src/BIVSOU.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/BIVSOU.f')
-rw-r--r--Dragon/src/BIVSOU.f277
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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:22:52 +0000