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
230
231
232
233
234
235
236
237
238
239
240
241
|
*DECK PSOUSN
SUBROUTINE PSOUSN(NUNF,NUNS,IG,IPTRK,JPTRK,KPMACR,NANIS,NREG,NMAT,
> NGRP1,NGRP2,MATCOD,FLUX,SOURCE)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the source from companion particle for the solution of SN
* equations.
*
*Copyright:
* Copyright (C) 2007 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
* NUNF first dimension of FLUX arrays.
* NUNS first dimension of SOURCE arrays.
* IG secondary group.
* IPTRK pointer to the tracking LCM object (from main particle).
* JPTRK pointer to the tracking LCM object (from companion particle).
* KPMACR pointer to the secondary-group related macrolib information.
* NANIS maximum cross section Legendre order.
* NREG number of regions.
* NMAT number of mixtures.
* NGRP1 number of primary energy groups.
* NGRP2 number of secondary energy groups.
* MATCOD mixture indices.
* FLUX fluxes.
*
*Parameters: output
* SOURCE sources.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK,JPTRK,KPMACR
INTEGER NUNF,NUNS,IG,NANIS,NREG,NMAT,NGRP1,NGRP2,MATCOD(NREG)
REAL FLUX(NUNF,NGRP1),SOURCE(NUNS,NGRP2)
*----
* LOCAL VARIABLES
*----
PARAMETER(NSTATE=40,PI4=12.5663706144)
INTEGER IPAR(NSTATE),JPAR(NSTATE),P,P2,IELEM,EELEM,IELEM2,EELEM2,
1 NM,NM2,EEL,IEL,IND,IND2,EL,EL2
CHARACTER CANIL*2
*----
* ALLOCATABLE ARRAYS
*----
INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,MAP
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: M_INDEXES
REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT
TYPE(C_PTR) IL_PTR,IM_PTR,IL2_PTR,IM2_PTR
INTEGER, POINTER, DIMENSION(:) :: IL,IM,IL2,IM2
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(IJJ(0:NMAT),NJJ(0:NMAT),IPOS(0:NMAT))
ALLOCATE(XSCAT(0:NMAT*NGRP1))
*----
* RECOVER SNT SPECIFIC PARAMETERS.
*----
CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR)
IF(IPAR(1).NE.NREG) CALL XABORT('PSOUSN: INCONSISTENT NREG.')
ITYPE=IPAR(6)
NSCT=IPAR(7)
IELEM=IPAR(8)
ISCAT=IPAR(16)
EELEM=IPAR(35)
CALL LCMGPD(IPTRK,'IL',IL_PTR)
CALL LCMGPD(IPTRK,'IM',IM_PTR)
CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /))
CALL C_F_POINTER(IM_PTR,IM,(/ NSCT /))
CALL LCMGET(JPTRK,'STATE-VECTOR',JPAR)
IF(JPAR(1).NE.NREG) CALL XABORT('PSOUSN: INCONSISTENT NREG.')
ITYPE2=IPAR(6)
NSCT2=JPAR(7)
IELEM2=JPAR(8)
ISCAT2=JPAR(16)
EELEM2=JPAR(35)
CALL LCMGPD(JPTRK,'IL',IL2_PTR)
CALL LCMGPD(JPTRK,'IM',IM2_PTR)
CALL C_F_POINTER(IL2_PTR,IL2,(/ NSCT2 /))
CALL C_F_POINTER(IM2_PTR,IM2,(/ NSCT2 /))
IF(ITYPE.NE.ITYPE2.OR.NSCT.NE.NSCT2.OR.ISCAT.NE.ISCAT2)
1 CALL XABORT('PSOUSN: INCONSISTENCE OF ANGULAR DISCRETISATION'
2 //'BETWEEN THE PARTICLE AND ITS COMPANION PARTICLE.')
*----
* MAPPING BETWEEN SPACE-ENERGY MOMENTS
*----
NMX=IELEM
NMX2=IELEM2
IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) THEN
NMX=IELEM**2
NMX2=IELEM2**2
ELSEIF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN
NMX=IELEM**3
NMX2=IELEM2**3
ENDIF
ALLOCATE(M_INDEXES(MAX(NMX,NMX2),MAX(EELEM,EELEM2)))
ALLOCATE(MAP(NMX*EELEM))
M_INDEXES=0
MAP=0
DO IEL=1,NMX2
DO EEL=1,EELEM2
M_INDEXES(IEL,EEL)=EELEM2*(IEL-1)+EEL
ENDDO
ENDDO
DO IEL=1,NMX
DO EEL=1,EELEM
IF(IEL.LE.NMX2.AND.EEL.LE.EELEM2) THEN
IND=EELEM*(IEL-1)+EEL
MAP(IND)=M_INDEXES(IEL,EEL)
ENDIF
ENDDO
ENDDO
DEALLOCATE(M_INDEXES)
*----
* CONSTRUCT THE SOURCE.
*----
IJJ(0)=0
NJJ(0)=0
IPOS(0)=0
XSCAT(0)=0.0
IOF0=0
DO 100 P=1,NSCT
ILP = IL(P)
IF(ILP.GT.NANIS-1) GO TO 100
WRITE(CANIL,'(I2.2)') ILP
CALL LCMGET(KPMACR,'NJJS'//CANIL,NJJ(1))
CALL LCMGET(KPMACR,'IJJS'//CANIL,IJJ(1))
CALL LCMGET(KPMACR,'IPOS'//CANIL,IPOS(1))
CALL LCMGET(KPMACR,'SCAT'//CANIL,XSCAT(1))
IF((ITYPE.EQ.2).OR.(ITYPE.EQ.4)) THEN
*----
* SLAB OR SPHERICAL 1D CASE.
*----
NM=IELEM*EELEM
NM2=IELEM2*EELEM2
DO 20 IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) GO TO 20
DO 15 IEL=1,NM
IF(MAP(IEL).EQ.0) CONTINUE
IND=(IR-1)*NSCT*NM+NM*(P-1)+IEL
IND2=(IR-1)*NSCT*NM2+NM2*(P-1)+MAP(IEL)
JG=IJJ(IBM)
DO 10 JND=1,NJJ(IBM)
SOURCE(IND,IG)=SOURCE(IND,IG)+FLUX(IND2,JG)*
> XSCAT(IPOS(IBM)+JND-1)
JG=JG-1
10 CONTINUE
15 CONTINUE
20 CONTINUE
ELSE IF(ITYPE.EQ.3) THEN
*----
* CYLINDRICAL 1D CASE.
*----
DO 50 P2=0,(P-1)
IF(MOD(P-1+P2,2).EQ.1) GO TO 50
IOF0=IOF0+1
DO 40 IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) GO TO 40
IND=(IR-1)*NSCT+IOF0
JG=IJJ(IBM)
DO 30 JND=1,NJJ(IBM)
SOURCE(IND,IG)=SOURCE(IND,IG)+FLUX(IND,JG)*
> XSCAT(IPOS(IBM)+JND-1)
JG=JG-1
30 CONTINUE
40 CONTINUE
50 CONTINUE
ELSE IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) THEN
*----
* 2D CASES (CARTESIAN OR R-Z).
*----
NM=IELEM*IELEM*EELEM
NM2=IELEM2*IELEM2*EELEM2
DO 70 IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) GO TO 70
DO 65 IEL=1,IELEM**2
DO 64 EEL=1,EELEM
IF(IEL.GT.IELEM2**2.OR.EEL.GT.EELEM2) CONTINUE
EL=EELEM*(IEL-1)+EEL
EL2=EELEM2*(IEL-1)+EEL
IND=(IR-1)*NSCT*NM+NM*(P-1)+EL
IND2=(IR-1)*NSCT*NM2+NM2*(P-1)+EL2
JG=IJJ(IBM)
DO 60 JND=1,NJJ(IBM)
SOURCE(IND,IG)=SOURCE(IND,IG)+FLUX(IND2,JG)*
> XSCAT(IPOS(IBM)+JND-1)
JG=JG-1
60 CONTINUE
64 CONTINUE
65 CONTINUE
70 CONTINUE
ELSE IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN
*----
* 3D CARTESIAN CASE
*----
NM=IELEM*IELEM*IELEM*EELEM
NM2=IELEM2*IELEM2*IELEM2*EELEM2
DO 90 IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) GO TO 90
DO 85 IEL=1,IELEM**3
DO 84 EEL=1,EELEM
IF(IEL.GT.IELEM2**3.OR.EEL.GT.EELEM2) CONTINUE
EL=EELEM*(IEL-1)+EEL
EL2=EELEM2*(IEL-1)+EEL
IND=(IR-1)*NSCT*NM+NM*(P-1)+EL
IND2=(IR-1)*NSCT*NM2+NM2*(P-1)+EL2
JG=IJJ(IBM)
DO 80 JND=1,NJJ(IBM)
SOURCE(IND,IG)=SOURCE(IND,IG)+FLUX(IND2,JG)*
> XSCAT(IPOS(IBM)+JND-1)
JG=JG-1
80 CONTINUE
84 CONTINUE
85 CONTINUE
90 CONTINUE
ELSE
CALL XABORT('PSOUSN: TYPE OF DISCRETIZATION NOT IMPLEMENTED.')
ENDIF
100 CONTINUE
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(XSCAT,MAP)
DEALLOCATE(IPOS,NJJ,IJJ)
RETURN
END
|
