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
|
*DECK BIVFIS
SUBROUTINE BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL,
> XSCHI,XSNUF,FUNKNO,SUNKNO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the fission source for a BIVAC tracking.
*
*Copyright:
* Copyright (C) 2025 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
* IPTRK pointer to the tracking LCM object.
* NANIS maximum cross section Legendre order.
* NREG number of regions.
* NMAT number of mixtures.
* NIFIS number of fissile isotopes.
* NUNKNO number of unknowns per energy group.
* NGRP number of energy groups.
* MATCOD mixture indices.
* VOL volumes.
* XSCHI fission spectra.
* XSNUP nu times the fission cross sections.
* FUNKNO fluxes.
*
*Parameters: output
* SUNKNO sources.
*
*-----------------------------------------------------------------------
*
USE GANLIB
USE DOORS_MOD
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK
INTEGER NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD(NREG)
REAL VOL(NREG),XSCHI(0:NMAT,NIFIS,NGRP),XSNUF(0:NMAT,NIFIS,NGRP),
1 FUNKNO(NUNKNO,NGRP),SUNKNO(NUNKNO,NGRP)
*----
* LOCAL VARIABLES
*----
PARAMETER(NSTATE=40)
INTEGER JPAR(NSTATE),IJ1(25),IJ2(25)
LOGICAL CYLIND
CHARACTER*12 CXDOOR
*----
* ALLOCATABLE ARRAYS
*----
INTEGER, ALLOCATABLE, DIMENSION(:) :: KN,IDL
REAL, ALLOCATABLE, DIMENSION(:) :: XX,DD,FXSOR
REAL, ALLOCATABLE, DIMENSION(:,:) :: RR,RS
*----
* RECOVER BIVAC SPECIFIC PARAMETERS.
*----
CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR)
IF(JPAR(1).NE.NREG) CALL XABORT('BIVFIS: INCONSISTENT NREG.')
IF(JPAR(2).NE.NUNKNO) CALL XABORT('BIVFIS: INCONSISTENT NUNKNO.')
ITYPE=JPAR(6)
IELEM=JPAR(8)
ICOL=JPAR(9)
ISPLH=JPAR(10)
L4=JPAR(11)
LX=JPAR(12)
NLF=JPAR(14)
ISCAT=JPAR(16)
CYLIND=(ITYPE.EQ.3).OR.(ITYPE.EQ.6)
IF(ICOL.EQ.4) THEN
CALL XABORT('BIVFIS: COLLOCATION NODAL NOT IMPLEMENTED.')
ELSE IF((ITYPE.NE.2).AND.(ITYPE.NE.5)) THEN
CALL XABORT('BIVFIS: CARTESIAN 1D OR 2D GEOMETRY EXPECTED.')
ENDIF
CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
ALLOCATE(XX(NREG),DD(NREG),KN(MAXKN),IDL(NREG))
CALL LCMGET(IPTRK,'XX',XX)
CALL LCMGET(IPTRK,'DD',DD)
CALL LCMGET(IPTRK,'KN',KN)
CALL LCMGET(IPTRK,'KEYFLX',IDL)
IF(IELEM.LT.0) THEN
! Lagrangian finite element method
*----
* RECOVER THE FINITE ELEMENT UNIT STIFFNESS MATRIX.
*----
CALL LCMSIX(IPTRK,'BIVCOL',1)
CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
ALLOCATE(RR(LC,LC),RS(LC,LC))
CALL LCMGET(IPTRK,'R',RR)
CALL LCMGET(IPTRK,'RS',RS)
CALL LCMSIX(IPTRK,' ',2)
*----
* COMPUTE VECTORS IJ1 AND IJ2
*----
LL=LC*LC
DO I=1,LL
IJ1(I)=1+MOD(I-1,LC)
IJ2(I)=1+(I-IJ1(I))/LC
ENDDO
*----
* COMPUTE THE SOURCE
*----
NUM1=0
DO IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) CYCLE
IF(VOL(IR).EQ.0.0) GO TO 10
DO I=1,LL
IND1=KN(NUM1+I)
IF(IND1.EQ.0) CYCLE
I1=IJ1(I)
I2=IJ2(I)
DO J=1,LL
IND2=KN(NUM1+J)
IF(IND2.EQ.0) CYCLE
IF(CYLIND) THEN
AUXX=(RR(I1,IJ1(J))+RS(I1,IJ1(J))*XX(IR)/DD(IR))*
> RR(I2,IJ2(J))*VOL(IR)
ELSE
AUXX=RR(I1,IJ1(J))*RR(I2,IJ2(J))*VOL(IR)
ENDIF
DO IG=1,NGRP
DO JG=1,NGRP
DO IS=1,NIFIS
SIGG=XSCHI(IBM,IS,IG)*XSNUF(IBM,IS,JG)
SUNKNO(IND1,IG)=SUNKNO(IND1,IG)+AUXX*SIGG*
> FUNKNO(IND2,JG)
ENDDO ! IS
ENDDO ! JG
ENDDO ! IG
ENDDO ! J
ENDDO ! I
10 NUM1=NUM1+LL
ENDDO ! IR
DEALLOCATE(RS,RR)
! append the integrated volumic sources
ALLOCATE(FXSOR(NUNKNO))
DO IS=1,NIFIS
FXSOR(:NUNKNO)=0.0
DO IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.LE.0) CYCLE
DO IG=1,NGRP
FXSOR(IDL(IR))=FXSOR(IDL(IR))+VOL(IR)*XSNUF(IBM,IS,IG)*
> FUNKNO(IDL(IR),IG)
ENDDO ! IG
DO IG=1,NGRP
SUNKNO(IDL(IR),IG)=SUNKNO(IDL(IR),IG)+XSCHI(IBM,IS,IG)*
> FXSOR(IDL(IR))
ENDDO ! IG
ENDDO ! IR
ENDDO ! IS
DEALLOCATE(FXSOR)
ELSE
! Raviart-Thomas finite element method
CXDOOR='BIVAC'
ALLOCATE(FXSOR(NUNKNO))
DO IS=1,NIFIS
FXSOR(:NUNKNO)=0.0
DO IG=1,NGRP
CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG),
> FXSOR,FUNKNO(1,IG))
ENDDO
DO IR=1,NREG
IBM=MATCOD(IR)
IF(IBM.EQ.0) CYCLE
DO IE=1,IELEM**2
IND=IDL(IR)+IE-1
IF(IND.EQ.0) CYCLE
DO IG=1,NGRP
SUNKNO(IND,IG)=SUNKNO(IND,IG)+XSCHI(IBM,IS,IG)*
> FXSOR(IND)
ENDDO
ENDDO
ENDDO
ENDDO ! IS
DEALLOCATE(FXSOR)
ENDIF
DEALLOCATE(IDL,KN,DD,XX)
RETURN
END
|
