summaryrefslogtreecommitdiff
path: root/Dragon/src/SYBRX2.f
blob: fa8974c9e5e8839a438b01d1da07c02c373f1fc2 (plain)
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
*DECK SYBRX2
      SUBROUTINE SYBRX2 (IPAS,NPIJ,NPIS,SIGT,SIGW,P,IMPX,NCOUR,IWIGN,
     1 NMCEL,NMERGE,NGEN,IQUAD,XX,YY,NMC,RAYRE,MAIL,RZMAIL,IFR,ALB,
     2 INUM,IGEN)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the global scattering-reduced collision probabilities in a
* 2-D Cartesian or hexagonal assembly using the interface current
* method with Roth approximation.
*
*Copyright:
* Copyright (C) 2002 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
* IPAS    total number of volumes.
* NPIJ    length of cellwise scattering-reduced collision probability
*         matrices.
* NPIS    length of cellwise scattering-reduced collision probability
*         matrices (NPIS=NMC(NGEN+1)).
* SIGT    total macroscopic cross sections.
* SIGW    within group scattering cross sections.
* IMPX    print flag (equal to 0 for no print).
* NCOUR   number of currents surrounding the cells (=4 Cartesian
*         lattice; =6 hexagonal lattice).
* IWIGN   type of cylinderization.
* IQUAD   quadrature parameters.
* NMCEL   total number of cells in the domain.
* IFR     index-number of in-currents.
* ALB     transmission/albedo associated with each in-current.
* NMERGE  total number of merged cells for which specific values
*         of the neutron flux and reactions rates are required.
*         Many cells with different position in the domain can
*         be merged before the neutron flux calculation if they
*         own the same generating cell. This allows some reduction
*         in cpu time and memory (NMERGE.le.NMCEL).
* INUM    index-number of the merged cell associated to each cell.
*         Note: IFR and ALB contains information to rebuild the
*         geometrical 'A' matrix.
* NGEN    total number of generating cells. A generating cell is
*         defined by its material and its position in the domain
*         (NGEN.le.NMERGE).
* XX      X-thickness of the generating cells.
* YY      Y-thickness of the generating cells.
* NMC     offset of the first volume in each generating cell.
* RAYRE   radius of the tubes in each generating cell.
* MAIL    offset of the first tracking information in each
*         generating cell.
* RZMAIL  real tracking information.
* IGEN    index-number of the generating cell associated with each
*         merged cell.
*
*Parameters: output
* P       reduced collision probabilities.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER IPAS,NPIJ,NPIS,IMPX,NCOUR,IWIGN,NMCEL,NMERGE,NGEN,
     1 IQUAD(4),NMC(NGEN+1),MAIL(2,NGEN),IFR(NCOUR*NMCEL),INUM(NMCEL),
     2 IGEN(NMERGE)
      REAL SIGT(IPAS),SIGW(IPAS),P(IPAS,IPAS),XX(NGEN),YY(NGEN),
     1 RAYRE(NPIS),RZMAIL(*),ALB(NCOUR*NMCEL)
*----
*  LOCAL VARIABLES
*----
      REAL PIBB(6)
*----
*  ALLOCATABLE ARRAYS
*----
      REAL, ALLOCATABLE, DIMENSION(:) :: SIGT2,SIGW2,PIJW,PISW,PSJW,PSSW
      REAL, ALLOCATABLE, DIMENSION(:,:) :: PSSB
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(PSSB(NMERGE,2*NMERGE),SIGT2(IPAS),SIGW2(IPAS),PIJW(NPIJ),
     1 PISW(NPIS),PSJW(NPIS),PSSW(NGEN))
*
      DO 20 I=1,IPAS
      DO 10 J=1,IPAS
      P(I,J)=0.0
   10 CONTINUE
   20 CONTINUE
      I1=0
      DO 40 IKK=1,NMERGE
      IKG=IGEN(IKK)
      J1=NMC(IKG)
      I2=NMC(IKG+1)-J1
      DO 30 I=1,I2
      SIGT2(J1+I)=SIGT(I1+I)
      SIGW2(J1+I)=SIGW(I1+I)
   30 CONTINUE
      I1=I1+I2
   40 CONTINUE
*
      CALL SYB002 (NGEN,NPIJ,NPIS,SIGT2,SIGW2,IMPX,NCOUR,IWIGN,IQUAD,
     1 XX,YY,NMC,RAYRE,MAIL,RZMAIL,PIJW,PISW,PSJW,PSSW)
*
      IPIJ=0
      DO 80 JKG=1,NGEN
      J2=NMC(JKG+1)-NMC(JKG)
      I1=0
      DO 70 IKK=1,NMERGE
      IKG=IGEN(IKK)
      I2=NMC(IKG+1)-NMC(IKG)
      IF(IKG.EQ.JKG) THEN
         DO 60 J=1,J2
         DO 50 I=1,J2
         P(I1+I,I1+J)=PIJW(IPIJ+(J-1)*J2+I)
   50    CONTINUE
   60    CONTINUE
      ENDIF
      I1=I1+I2
   70 CONTINUE
      IPIJ=IPIJ+J2*J2
   80 CONTINUE
*----
*  COMPUTATION OF PSSB=A*(I-PSS*A)**-1
*----
      DO 100 I=1,NMERGE
      DO 90 J=1,NMERGE
      PSSB(I,J)=0.0
      PSSB(I,NMERGE+J)=0.0
   90 CONTINUE
      PSSB(I,I)=1.0
  100 CONTINUE
      DO 130 ICEL=1,NMCEL
      IKK=INUM(ICEL)
      IKG=IGEN(IKK)
      IS=NCOUR*(ICEL-1)
      IF(NCOUR.EQ.4) THEN
         A=XX(IKG)
         B=YY(IKG)
         DEN1=2.0*(A+B)
         PIBB(1)=B/DEN1
         PIBB(2)=B/DEN1
         PIBB(3)=A/DEN1
         PIBB(4)=A/DEN1
      ELSE
         DO 110 JC=1,NCOUR
         PIBB(JC)=1.0/6.0
  110    CONTINUE
      ENDIF
      ZZZ=PSSW(IKG)
      DO 120 JC=1,NCOUR
      J1=IFR(IS+JC)
      ALBEDO=PIBB(JC)*ALB(IS+JC)
      PSSB(J1,NMERGE+IKK)=PSSB(J1,NMERGE+IKK)+ALBEDO
      PSSB(J1,IKK)=PSSB(J1,IKK)-ZZZ*ALBEDO
  120 CONTINUE
  130 CONTINUE
      CALL ALSB(NMERGE,NMERGE,PSSB,IER,NMERGE)
      IF(IER.NE.0) CALL XABORT('SYBRX2: SINGULAR MATRIX.')
*----
*  COMPUTATION OF PIS*PSSB*PSJ
*----
      I1=0
      DO 170 IKK=1,NMERGE
      IKG=IGEN(IKK)
      I1P=NMC(IKG)
      I2=NMC(IKG+1)-I1P
      DO 160 I=1,I2
      ZZZ=PISW(I1P+I)
      J1=0
      DO 150 JKK=1,NMERGE
      JKG=IGEN(JKK)
      J1P=NMC(JKG)
      J2=NMC(JKG+1)-J1P
      DO 140 J=1,J2
      P(I1+I,J1+J)=P(I1+I,J1+J)+ZZZ*PSSB(JKK,NMERGE+IKK)*PSJW(J1P+J)
  140 CONTINUE
      J1=J1+J2
  150 CONTINUE
  160 CONTINUE
      I1=I1+I2
  170 CONTINUE
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(PSSW,PSJW,PISW,PIJW,SIGW2,SIGT2,PSSB)
      RETURN
      END