1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
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
|
