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
|
*DECK MCGFFAL
SUBROUTINE MCGFFAL(SUBSCH,K,KPN,M,N,H,NOM,NZON,WEIGHT,XST,SP,SM,
1 DSP,DSM,NREG,NMU,NANI,NFUNLX,TRHAR,KEYCUR,IMU,B,F,
2 PHIV,DPHIV)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solution of transport equation on a finite track.
* Linear-discontinuous-characteristics approximation.
* Ray-tracing (anisotropic scattering case,'source term isolation' off).
*
*Copyright:
* Copyright (C) 2015 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
* SUBSCH track coefficients calculation subroutine.
* K total number of volumes for which specific values
* of the neutron flux and reactions rates are required.
* KPN total number of unknowns in vectors F.
* M number of material mixtures.
* N number of elements in the current track.
* H vector containing the lenght of the different segments of this
* track.
* NOM vector containing the region number of the different segments
* of this track.
* NZON index-number of the mixture type assigned to each volume.
* WEIGHT track weight.
* XST macroscopic total cross section.
* SP total source vector for + direction.
* SM total source vector for - direction.
* DSP linear component of the total source vector for + direction.
* DSM linear component of the total source vector for - direction.
* NREG number of volumes.
* NMU order of the polar quadrature set.
* NANI scattering anisotropy (=1 for isotropic scattering).
* NFUNLX number of moments of the spherical harmonics.
* NMOD third dimension of TRHAR.
* TRHAR spherical harmonics components for this angle in the plane.
* KEYCUR position of current elements in PHI vector.
* IMU azimuthal angle corresponding to this track.
*
*Parameters: input/output
* F vector containing the zonal scalar flux (surface components).
* PHIV vector containing the zonal scalar flux (component 1).
* DPHIV vector containing the zonal scalar flux (components 2 and 3).
*
*Parameters: scratch
* B undefined.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER K,KPN,M,N,NOM(N),NZON(K),NMU,NFUNLX,NREG,
1 KEYCUR(K-NREG),IMU,NANI
REAL XST(0:M),TRHAR(NMU,NFUNLX,2)
DOUBLE PRECISION WEIGHT,H(N),SP(N),SM(N),DSP(N),DSM(N),B(0:5,N),
1 F(KPN),PHIV(NFUNLX,NREG),DPHIV(2*NFUNLX,NREG)
EXTERNAL SUBSCH
*----
* LOCAL VARIABLES
*---
REAL ETA,XI
DOUBLE PRECISION F0,DF0,SI,SJ,DSI,DSJ,RM,RP,DSIG,DH
INTEGER I,NOMI,IND1,INDN,JF,J,NOMJ
*----
* Calculation of the coefficients for this track.
*----
* MCGSCAL: Linear discontinuous-Characteristics Scheme with
* Tabulated Exponentials
* MCGDDFL: Diamond-Differencing DD1 Scheme
* MCGSCEL: Linear discontinuous-Characteristics Scheme with
* Exact Exponentials
IF(NANI.LE.0) CALL XABORT('MCGFFAL: INVALID VALUE OF NANI.')
CALL SUBSCH(N,K,M,NOM,NZON,H,XST,B)
*----
* Summation along the track in both directions
*----
* incoming flux in + direction
IND1=KEYCUR(-NOM(1))
RP=SP(1)
* incoming flux in - direction
INDN=KEYCUR(-NOM(N))
RM=SM(N)
* track angles in 3D
ETA=TRHAR(IMU,3,1)
XI=TRHAR(IMU,2,1)
DO I=2,N-1
* + direction
NOMI=NOM(I)
DSIG=XST(NZON(NOMI))
DH=H(I)
SI=SP(I)
DSI=DSP(I)
F0=B(1,I)*RP+B(2,I)*SI+B(3,I)*DSI
DF0=B(4,I)*RP-B(3,I)*SI/(DH*DH)+B(5,I)*DSIG*DSI
RP=B(0,I)*RP+B(1,I)*SI-B(3,I)*DSIG*DSI
DO JF=1,NFUNLX
PHIV(JF,NOMI)=PHIV(JF,NOMI)+WEIGHT*F0*TRHAR(IMU,JF,1)
DPHIV(JF,NOMI)=DPHIV(JF,NOMI)+WEIGHT*DF0*ETA*
> TRHAR(IMU,JF,1)
DPHIV(NFUNLX+JF,NOMI)=DPHIV(NFUNLX+JF,NOMI)+WEIGHT*DF0*XI*
> TRHAR(IMU,JF,1)
ENDDO
* - direction
J=N+1-I
NOMJ=NOM(J)
DSIG=XST(NZON(NOMJ))
DH=H(J)
SJ=SM(J)
DSJ=DSM(J)
F0=B(1,J)*RM+B(2,J)*SJ+B(3,J)*DSJ
DF0=B(4,J)*RM-B(3,J)*SJ/(DH*DH)+B(5,J)*DSIG*DSJ
RM=B(0,J)*RM+B(1,J)*SJ-B(3,J)*DSIG*DSJ
DO JF=1,NFUNLX
PHIV(JF,NOMJ)=PHIV(JF,NOMJ)+WEIGHT*F0*TRHAR(IMU,JF,2)
DPHIV(JF,NOMJ)=DPHIV(JF,NOMJ)-WEIGHT*DF0*ETA*
> TRHAR(IMU,JF,2)
DPHIV(NFUNLX+JF,NOMJ)=DPHIV(NFUNLX+JF,NOMJ)-WEIGHT*DF0*XI*
> TRHAR(IMU,JF,2)
ENDDO
ENDDO
* outgoing flux in + direction
F(INDN)=F(INDN)+WEIGHT*RP
* outgoing flux in - direction
F(IND1)=F(IND1)+WEIGHT*RM
*
RETURN
END
|
