diff options
| author | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
|---|---|---|
| committer | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
| commit | 7dfcc480ba1e19bd3232349fc733caef94034292 (patch) | |
| tree | 03ee104eb8846d5cc1a981d267687a729185d3f3 /Trivac/src/BIVASM.f | |
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Trivac/src/BIVASM.f')
| -rwxr-xr-x | Trivac/src/BIVASM.f | 229 |
1 files changed, 229 insertions, 0 deletions
diff --git a/Trivac/src/BIVASM.f b/Trivac/src/BIVASM.f new file mode 100755 index 0000000..6d86410 --- /dev/null +++ b/Trivac/src/BIVASM.f @@ -0,0 +1,229 @@ +*DECK BIVASM
+ SUBROUTINE BIVASM(HNAMT,ITY,IPTRK,IPSYS,IMPX,NBMIX,NEL,NLF,NDIM,
+ 1 NALBP,MAT,VOL,GAMMA,SGD)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Assembling of a single-group system matrix for BIVAC.
+*
+*Copyright:
+* Copyright (C) 2008 Ecole Polytechnique de Montreal.
+*
+*Author(s): A. Hebert
+*
+*Parameters: input/output
+* HNAMT name of the matrix.
+* ITY type of assembly: =0: leakage-removal matrix assembly;
+* =1: cross section matrix assembly.
+* IPTRK L_TRACK pointer to the BIVAC tracking information.
+* IPSYS L_SYSTEM pointer to system matrices.
+* IMPX print parameter (equal to zero for no print).
+* NBMIX total number of material mixtures.
+* NEL total number of finite elements.
+* NLF number of Legendre orders for the flux (even number). Equal
+* to zero for diffusion theory.
+* NDIM second dimension of matrix SGD.
+* NALBP number of physical albedos.
+* MAT mixture index assigned to each volume.
+* VOL volume of each element.
+* GAMMA physical albedo functions.
+* SGD nuclear properties per material mixture.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ CHARACTER HNAMT*(*)
+ TYPE(C_PTR) IPTRK,IPSYS
+ INTEGER ITY,IMPX,NBMIX,NEL,NLF,NDIM,NALBP,MAT(NEL)
+ REAL VOL(NEL),GAMMA(NALBP),SGD(NBMIX,NDIM)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (NSTATE=40)
+ LOGICAL CYLIND
+ CHARACTER TEXT11*11
+ INTEGER ITP(NSTATE)
+ INTEGER, DIMENSION(:), ALLOCATABLE :: KN,IQFR,MU,IPERT
+ REAL, DIMENSION(:), ALLOCATABLE :: XX,YY,DD,QFR
+ REAL, DIMENSION(:,:), ALLOCATABLE :: R,RS,Q,QS,V,H,RH,QH,RT,QT
+ REAL, DIMENSION(:), POINTER :: SYS,ASS
+ TYPE(C_PTR) SYS_PTR,ASS_PTR
+*----
+* RECOVER BIVAC SPECIFIC TRACKING INFORMATION
+*----
+ CALL LCMGET(IPTRK,'STATE-VECTOR',ITP)
+ ITYPE=ITP(6)
+ CYLIND=(ITYPE.EQ.3).OR.(ITYPE.EQ.6)
+ IELEM=ITP(8)
+ ICOL=ITP(9)
+ ISPLH=ITP(10)
+ LL4=ITP(11)
+ LX=ITP(12)
+ LY=ITP(13)
+ NVD=ITP(17)
+ CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
+ CALL LCMLEN(IPTRK,'QFR',MAXQF,ITYLCM)
+ ALLOCATE(XX(LX*LY),YY(LX*LY),DD(LX*LY),KN(MAXKN),QFR(MAXQF),
+ 1 IQFR(MAXQF),MU(LL4))
+ IF(ITYPE.EQ.8) THEN
+ CALL LCMGET(IPTRK,'SIDE',SIDE)
+ ELSE
+ CALL LCMGET(IPTRK,'XX',XX)
+ CALL LCMGET(IPTRK,'YY',YY)
+ CALL LCMGET(IPTRK,'DD',DD)
+ ENDIF
+ CALL LCMGET(IPTRK,'KN',KN)
+ CALL LCMGET(IPTRK,'QFR',QFR)
+ CALL LCMGET(IPTRK,'IQFR',IQFR)
+ CALL LCMGET(IPTRK,'MU',MU)
+*----
+* APPLY PHYSICAL ALBEDO FUNCTIONS
+*----
+ IF(NALBP.GT.0) THEN
+ DO IQW=1,MAXQF
+ IALB=IQFR(IQW)
+ IF(IALB.NE.0) QFR(IQW)=QFR(IQW)*GAMMA(IALB)
+ ENDDO
+ ENDIF
+*
+ TEXT11=HNAMT
+ IF(IMPX.GT.0) WRITE(6,'(/36H BIVASM: ASSEMBLY OF SYMMETRIC MATRI,
+ 1 3HX '',A11,38H'' IN COMPRESSED DIAGONAL STORAGE MODE.)') TEXT11
+*----
+* ASSEMBLY OF THE SYSTEM MATRICES
+*----
+ CALL KDRCPU(TK1)
+ IIMAX=MU(LL4)
+ IF(NLF.NE.0) IIMAX=IIMAX*NLF/2
+ SYS_PTR=LCMARA(IIMAX)
+ CALL C_F_POINTER(SYS_PTR,SYS,(/ IIMAX /))
+ SYS(:IIMAX)=0.0
+*
+ IF((IELEM.LT.0).AND.(ITYPE.NE.8)) THEN
+* MESH CORNER FINITE DIFFERENCES OR LAGRANGIAN FINITE ELEMENTS
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+ ALLOCATE(R(LC,LC),RS(LC,LC),Q(LC,LC),QS(LC,LC))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'RS',RS)
+ CALL LCMGET(IPTRK,'Q',Q)
+ CALL LCMGET(IPTRK,'QS',QS)
+ CALL LCMSIX(IPTRK,' ',2)
+ CALL BIVA01(ITY,MAXKN,SGD,CYLIND,NEL,LL4,NBMIX,IIMAX,
+ 1 XX,YY,DD,MAT,KN,QFR,VOL,MU,LC,R,RS,Q,QS,SYS)
+ DEALLOCATE(R,RS,Q,QS)
+ ELSE IF((IELEM.GT.0).AND.(ITYPE.NE.8).AND.(NLF.GT.0)) THEN
+* MIXED-DUAL FINITE ELEMENTS (SIMPLIFIED PN THEORY)
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+ ALLOCATE(R(LC,LC),V(LC,LC-1))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'V',V)
+ CALL LCMSIX(IPTRK,' ',2)
+ CALL PNDM2E(ITY,NEL,LL4,IELEM,ICOL,MAT,VOL,NBMIX,NLF,NVD,
+ 1 NDIM/2,SGD(1,1),SGD(1,1+NDIM/2),XX,YY,KN,QFR,MU,IIMAX,LC,
+ 2 R,V,SYS)
+ DEALLOCATE(R,V)
+ ELSE IF((IELEM.GT.0).AND.(ITYPE.NE.8)) THEN
+* MIXED-DUAL FINITE ELEMENTS (DIFFUSION THEORY).
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+ ALLOCATE(R(LC,LC),V(LC,LC-1))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'V',V)
+ CALL LCMSIX(IPTRK,' ',2)
+ CALL BIVA02(ITY,SGD,CYLIND,IELEM,ICOL,NEL,LL4,NBMIX,IIMAX,XX,
+ 1 YY,DD,MAT,KN,QFR,VOL,MU,LC,R,V,SYS)
+ DEALLOCATE(R,V)
+ ELSE IF((IELEM.LT.0).AND.(ITYPE.EQ.8)) THEN
+* MESH CORNER FINITE DIFFERENCES FOR HEXAGONS
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ ALLOCATE(R(2,2),RH(6,6),QH(6,6),RT(3,3),QT(3,3))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'RH',RH)
+ CALL LCMGET(IPTRK,'QH',QH)
+ CALL LCMGET(IPTRK,'RT',RT)
+ CALL LCMGET(IPTRK,'QT',QT)
+ CALL LCMSIX(IPTRK,' ',2)
+ IF(ISPLH.EQ.1) THEN
+ NELEM=MAXKN/7
+ ELSE
+ NELEM=MAXKN/4
+ ENDIF
+ CALL BIVA03(ITY,MAXKN,MAXQF,SGD,NEL,LL4,ISPLH,NELEM,NBMIX,
+ 1 IIMAX,SIDE,MAT,KN,QFR,VOL,MU,R,RH,QH,RT,QT,SYS)
+ DEALLOCATE(R,RH,QH,RT,QT)
+ ELSE IF((IELEM.GT.0).AND.(ITYPE.EQ.8).AND.(ICOL.EQ.4)) THEN
+* MESH CENTERED FINITE DIFFERENCES FOR HEXAGONS
+ CALL BIVA04(ITY,MAXKN,MAXQF,SGD,NEL,LL4,ISPLH,NBMIX,IIMAX,
+ 1 SIDE,MAT,KN,QFR,VOL,MU,SYS)
+ ELSE IF((IELEM.GT.0).AND.(ITYPE.EQ.8).AND.(NLF.GT.0)) THEN
+* THOMAS-RAVIART-SCHNEIDER METHOD FOR HEXAGONS (SIMPLIFIED PN
+* THEORY)
+ LXH=LX/(3*ISPLH**2)
+ NBLOS=LXH*ISPLH**2
+ ALLOCATE(IPERT(NBLOS))
+ CALL LCMGET(IPTRK,'IPERT',IPERT)
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+ ALLOCATE(R(LC,LC),V(LC,LC-1),H(LC,LC-1))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'V',V)
+ CALL LCMGET(IPTRK,'H',H)
+ CALL LCMSIX(IPTRK,' ',2)
+ CALL PNDH2E(ITY,IELEM,ICOL,NBLOS,LL4,NBMIX,IIMAX,SIDE,MAT,
+ 1 IPERT,SGD(1,1),KN,QFR,NLF,NVD,NDIM/2,MU,LC,R,V,H,SYS)
+ DEALLOCATE(R,V,H,IPERT)
+ ELSE IF((IELEM.GT.0).AND.(ITYPE.EQ.8)) THEN
+* THOMAS-RAVIART-SCHNEIDER METHOD FOR HEXAGONS (DIFFUSION THEORY)
+ LXH=LX/(3*ISPLH**2)
+ NBLOS=LXH*ISPLH**2
+ ALLOCATE(IPERT(NBLOS))
+ CALL LCMGET(IPTRK,'IPERT',IPERT)
+ CALL LCMSIX(IPTRK,'BIVCOL',1)
+ CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
+ ALLOCATE(R(LC,LC),V(LC,LC-1),H(LC,LC-1))
+ CALL LCMGET(IPTRK,'R',R)
+ CALL LCMGET(IPTRK,'V',V)
+ CALL LCMGET(IPTRK,'H',H)
+ CALL LCMSIX(IPTRK,' ',2)
+ CALL BIVA05(ITY,SGD,IELEM,NBLOS,LL4,NBMIX,IIMAX,SIDE,MAT,IPERT,
+ 1 KN,QFR,MU,LC,R,V,H,SYS)
+ DEALLOCATE(R,V,H,IPERT)
+ ENDIF
+ CALL LCMPPD(IPSYS,TEXT11,IIMAX,2,SYS_PTR)
+ CALL KDRCPU(TK2)
+ IF(IMPX.GT.0) WRITE(6,'(/35H BIVASM: CPU TIME FOR SYSTEM MATRIX,
+ 1 11H ASSEMBLY =,F9.2,3H S.)') TK2-TK1
+*----
+* MATRIX FACTORIZATIONS
+*----
+ IF((ITY.EQ.0).OR.(TEXT11.EQ.'RM')) THEN
+ CALL KDRCPU(TK1)
+ ASS_PTR=LCMARA(IIMAX)
+ CALL C_F_POINTER(ASS_PTR,ASS,(/ IIMAX /))
+ CALL LCMGET(IPSYS,TEXT11,ASS)
+ IF(NLF.EQ.0) THEN
+ CALL ALLDLF(LL4,ASS(1),MU)
+ ELSE
+ IOF=1
+ DO 50 IL=0,NLF-2,2
+ CALL ALLDLF(LL4,ASS(IOF),MU)
+ IOF=IOF+MU(LL4)
+ 50 CONTINUE
+ ENDIF
+ CALL LCMPPD(IPSYS,'I'//TEXT11,IIMAX,2,ASS_PTR)
+ CALL KDRCPU(TK2)
+ IF(IMPX.GT.1) WRITE(6,'(/34H BIVASM: CPU TIME FOR LDLT FACTORI,
+ 1 18HZATION OF MATRIX '',A11,2H''=,F9.2,3H S.)') TEXT11,TK2-TK1
+ ENDIF
+*----
+* RELEASE BIVAC SPECIFIC TRACKING INFORMATION
+*----
+ DEALLOCATE(MU,IQFR,QFR,KN,DD,XX,YY)
+ RETURN
+ END
|