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
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
|
!
!-----------------------------------------------------------------------
!
!Purpose:
! To generate the standard tracking lines (isotropic tracking) for a
! geometry using the SALT algorithm.
!
!Copyright:
! Copyright (C) 2014 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
! IFTEMP pointer to a temporary tracking data structure in creation
! mode.
! IPRINT print level.
! IGTRK flag to generate the tracking file. In the case where IGTRK=1,
! the tracking is performed and used to evaluate the track
! normalisation factor and the tracking file is generated.
! When IGTRK=0, the tracking is still performed and used to
! evaluate the track normalisation factor but the tracking file
! is not generated.
! NFREG number of regions.
! NBANGL number of angles.
! NQUAD number of quarter (in 2-D).
! RENO track normalisation option. A value RENO=-1 implies
! a direction dependent normalization of the tracks for
! the volume. A value reno=0, implies a global normalisation.
! NBDR number of directions for track normalization (no normalization
! when reno=1).
! IFMT tracking file format:
! IFMT=0 for short file;
! IFMT=1 long file required by TLM: module.
! DENUSR user defined track density.
! DANGLT angle cosines.
! DDENWT angular density for each angle.
! GG geometry basic information.
!
!Parameters: output
! NBTDIR number of tracks directions considered.
! MAXSGL maximum number of segments in a line.
! NTLINE total number of lines generated.
! DVNOR ratio of analytic to tracked volume.
!
!-----------------------------------------------------------------------
!
SUBROUTINE SALTLS(IFTEMP,IPRINT,IGTRK,NFREG,NBANGL,NQUAD,RENO,NBDR,IFMT,DENUSR, &
DANGLT,DDENWT,GG,NBTDIR,MAXSGL,NTLINE,DVNOR)
USE PRECISION_AND_KINDS, ONLY : PDB,SMALL,PI,TWOPI,HALFPI,INFINITY
USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC
USE SAL_TRACKING_TYPES, ONLY : NMAX2,MINLEN,NNN,ITRAC2,RTRAC2,DPIECE,CNT,CNT0,NBTRAC,LGMORE, &
LGMORE,IERR,DD0,EX0,EY0,DELX,AX,AY,DINIT,ANGTAB,ELMTAB,NB_MAX,IMPX
USE SAL_AUX_MOD, ONLY : SAL231,SAL232,SAL235
USE SAL_TRAJECTORY_MOD, ONLY : SALTRA
IMPLICIT NONE
!----
! subroutine arguments
!----
INTEGER :: IFTEMP,IPRINT,IGTRK,NFREG,NBANGL,NQUAD,RENO,NBDR,IFMT,MAXSGL,NBTDIR,NTLINE,OK
REAL(PDB) :: DENUSR,DANGLT(2,NQUAD*NBANGL),DDENWT(NQUAD,NBANGL),DVNOR(NFREG,NBDR)
TYPE(T_G_BASIC) :: GG
!----
! local parameters
!----
INTEGER :: II0,IPHI,MQ,MQUAD,NPIECE,INODE,IANGL,II,IQUAD,FIRST,ICURR, &
JCURR,LASTI,NTSEG,IAVERR
LOGICAL :: LGON
REAL :: EPS0,X
REAL(PDB) :: WT,WR,DEL0,DEL2,DNEW,DOLD,ANGLE,DELM,KEEP_DELX, &
MOV_DELX,SUMM,NORM,XFACT,EPSILON_PDB,DCERR,DAVERR, &
DMVERR,DSVERR,TORIG(2)
REAL(PDB), DIMENSION(2) :: THETA0
REAL(PDB), ALLOCATABLE, DIMENSION(:) :: VOLN,SURFN,CURRN
REAL(PDB), ALLOCATABLE, DIMENSION(:,:) :: FACNRM ! aux for normalisation
REAL, PARAMETER :: EPS3 = 1E-3
INTEGER, PARAMETER :: FOUT =6
!
IMPX=IPRINT
!----
! recompute weights
!----
NBTDIR=0
SUMM=0._PDB
DO IANGL=1,NBANGL
DO IQUAD=1,NQUAD
IF(DDENWT(IQUAD,IANGL).GT.0._PDB) NBTDIR=NBTDIR+1
SUMM=SUMM+DDENWT(IQUAD,IANGL)
ENDDO
ENDDO
NORM=0.5_PDB/SUMM
NB_MAX=1
ALLOCATE(ANGTAB(2),ELMTAB(2))
!----
! isotropic tracking loop
!----
! define quadrature
! get limit intervals for radial quadrature
! minimum radial interval to contain one trajectory
EPSILON_PDB=SQRT(EPSILON(X))
EPS0=REAL(10.*EPSILON_PDB)
DEL0=REAL(EPS3,PDB)/DENUSR
DEL2=-1._PDB
! start a set of tracks
! CNT0 = address of beginning of trajectory - 1
! CNT = address for trajectory data - 1
! NBTRAC = number of trajectories in a record
NBTRAC=0
MAXSGL=0
IPHI=0
ALLOCATE(VOLN(GG%NB_NODE),SURFN(GG%NB_SURF2), &
CURRN(GG%NB_SURF2),FACNRM(GG%NB_NODE,NBANGL*NQUAD),STAT=OK)
IF(OK/=0) CALL XABORT('SALTLS: NOT ENOUGH MEMORY R')
FACNRM(:GG%NB_NODE,:NBANGL*NQUAD)=0._PDB
DO IANGL=1,NBANGL
DO IQUAD=1,NQUAD
IPHI=IPHI+1
! KEEP IPHI
EX0=DANGLT(1,(IANGL-1)*NQUAD+IQUAD)
EY0=DANGLT(2,(IANGL-1)*NQUAD+IQUAD)
ANGLE=DACOS(EX0)
! get theta- and theta+ from angle (to be used in SAL235 to
! decide whether to include projections of tangents to arcs):
THETA0(1)=ANGLE+HALFPI
THETA0(2)=THETA0(1)+PI
IF(THETA0(2)>TWOPI)THETA0(2)=THETA0(2)-TWOPI
! get projection of points onto axis orthogonal to tracking:
! only first and last points from macro perimeter
CALL SAL235(NPIECE,THETA0,EX0,EY0,GG%IPAR,GG%RPAR,GG%PERIM_MAC2,GG%NPERIM_MAC2)
! integrate on each piece
DOLD=DPIECE(1)
DNEW=DPIECE(2)
DELM=DNEW-DOLD
IF(DELM>DEL0)THEN
! compute nber of intervals for step =< 1/denusr
MQUAD=1+INT(DELM*DENUSR)
DELM=DELM/MQUAD
DO MQ=1,MQUAD
DELX=DOLD+0.5_PDB*DELM
KEEP_DELX=DELX
WR=REAL(DELM)
! initialize entering distance
DD0=-INFINITY
LGON=.TRUE.
DO WHILE (LGON)
CNT0=0
CNT=CNT0+NNN
IERR=0
MOV_DELX=0.
DO WHILE (IERR==0)
! compute one trajectory:
AX=DELX*EY0
AY=-DELX*EX0
CALL SALTRA(DANGLT,GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%ISURF2_ELEM,GG%IPAR, &
GG%RPAR,GG%PPERIM_NODE,GG%PERIM_NODE)
IF(IERR==0)THEN
! if IERR=0,trajectory has entered into the element joint, moves
! DELX -> DELX+EPSILON_PDB*N
MOV_DELX=MOV_DELX+EPS0
DELX=KEEP_DELX+MOV_DELX
CNT=CNT0+NNN
ENDIF
ENDDO
! a trajectory has been completed: store angle order nber, total weight and wr
!!! corrected by A. Hebert in ev3874
!!! ITRAC2(CNT0+4)=IPHI
ITRAC2(CNT0+7)=IPHI
RTRAC2(CNT0+7)=DDENWT(IQUAD,IANGL)*WR*NORM
RTRAC2(CNT0+8)=WR
II0=CNT0+1
IF(IPRINT > 3) CALL SAL231(RTRAC2(II0:),ITRAC2(II0:),DELX,EX0,EY0,ANGLE)
DO II=1,ITRAC2(II0)
IF(RTRAC2(II0+II+NNN-1) <= 0.0) CALL XABORT('SALTLS: INVALID SEGMENT LENGTH')
ENDDO
! compute volumes
CALL SAL232(ITRAC2(II0:),RTRAC2(II0:),FACNRM,GG,SURFN,CURRN)
! next line
IF(CNT+NNN+MINLEN>=NMAX2) CALL XABORT('SALTLS: BUFFER OVERFLOW')
NBTRAC=NBTRAC+1
WT=RTRAC2(II0+7-1)
WR=RTRAC2(II0+8-1)
FIRST=1
LASTI=ITRAC2(II0)
IF(GG%NB_SURF2/=0)THEN
ICURR=ITRAC2(II0+5-1)
JCURR=ITRAC2(II0+6-1)
ELSE
ICURR=0
JCURR=0
ENDIF
NTSEG=LASTI-FIRST+3
MAXSGL=MAX(MAXSGL,NTSEG)
IF(IGTRK == 1) THEN
IF(IFMT == 1) THEN
TORIG(1)=AX+DANGLT(1,IPHI)*DINIT ; TORIG(2)=AY+DANGLT(2,IPHI)*DINIT ;
WRITE(IFTEMP) 1,NTSEG,WT,IPHI, &
-ICURR,(ITRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),-JCURR, &
0.5D0,(RTRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),0.5D0, &
NBTRAC,1,MQ,1,TORIG(1),TORIG(2)
ELSE
WRITE(IFTEMP) 1,NTSEG,WT,IPHI, &
-ICURR,(ITRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),-JCURR, &
0.5D0,(RTRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),0.5D0
ENDIF
IF(IPRINT>5) WRITE(FOUT,'(2X,''TRAJ# DELX '',''IERR = '',I6,3X,1P,E12.4,I5)') &
NBTRAC,DELX,IERR
ENDIF
LGON=LGMORE
ENDDO
! end of one interval: move into beginning of next
DOLD=DOLD+DELM
ENDDO
! end of one piece
ENDIF
! end of trajectories for this angle
IF(IPRINT > 5) THEN
WRITE(FOUT,'('' ANGLE EX EY NTRA = '',1P,3E12.4,I8,/)') ANGLE,EX0,EY0,NBTRAC
ENDIF
ENDDO
ENDDO
NTLINE=NBTRAC
!----
! Compute merged normalization factors
!----
IF(RENO/=1) THEN
DO INODE=1,GG%NB_NODE
VOLN(INODE)=0._PDB
DO IANGL=1,NBANGL
DO IQUAD=1,NQUAD
VOLN(INODE)=VOLN(INODE)+2._PDB*FACNRM(INODE,IANGL+(IQUAD-1)*NBANGL)* &
DDENWT(IQUAD,IANGL)*NORM
ENDDO
ENDDO
ENDDO
DMVERR=0.0D0
DSVERR=0.0D0
DAVERR=0.0D0
IAVERR=0
DO INODE=1,GG%NB_NODE
DO IPHI=1,NBANGL*NQUAD
IF(RENO==0) THEN
IF(ABS(VOLN(INODE))>SMALL) THEN
FACNRM(INODE,IPHI)=GG%VOL_NODE(INODE)/VOLN(INODE)
ENDIF
ELSE IF(RENO==-1) THEN
IF(ABS(FACNRM(INODE,IPHI))>SMALL) THEN
FACNRM(INODE,IPHI)=GG%VOL_NODE(INODE)/FACNRM(INODE,IPHI)
ENDIF
ENDIF
ENDDO
IF(ABS(VOLN(INODE))>SMALL) THEN
IAVERR=IAVERR+1
DCERR=100.0D0*(1.0D0-GG%VOL_NODE(INODE)/VOLN(INODE))
DMVERR=MAX(DMVERR,ABS(DCERR))
DSVERR=DSVERR+DCERR*DCERR
DAVERR=DAVERR+DCERR
ENDIF
ENDDO
DSVERR=SQRT(DSVERR/DBLE(IAVERR))
DAVERR=DAVERR/DBLE(IAVERR)
IF(IPRINT > 0) WRITE(FOUT,6005) DSVERR,DMVERR,DAVERR
IF(IPRINT > 5) THEN
DO IPHI=1,NBANGL*NQUAD
WRITE(*,*) 'iphi=',IPHI
WRITE(*,*) 'facnrm : ',(FACNRM(INODE,IPHI),INODE=1,MIN(10,GG%NB_NODE))
ENDDO
ENDIF
IF(GG%NB_NODE > NFREG) CALL XABORT('SALTLS: bug.')
DVNOR(:,:)=0._PDB
DO INODE=1,GG%NB_NODE
IF(RENO==0) THEN
DVNOR(INODE,1)=DVNOR(INODE,1)+FACNRM(INODE,1)*GG%VOL_NODE(INODE)
ELSE IF(RENO==-1) THEN
DO IPHI=1,NBANGL*NQUAD
XFACT=FACNRM(INODE,IPHI)
DVNOR(INODE,IPHI+1)=DVNOR(INODE,IPHI+1)+XFACT*GG%VOL_NODE(INODE)
ENDDO
ENDIF
ENDDO
DO IPHI=1,NBDR
DVNOR(:,IPHI)=DVNOR(:,IPHI)/GG%VOL_NODE(:)
ENDDO
IF(IPRINT > 4) THEN
DO IPHI=1,NBDR
WRITE(*,*) 'iphi=',IPHI
WRITE(*,*) 'dvnor : ',(DVNOR(INODE,IPHI),INODE=1,MIN(10,NFREG))
ENDDO
ENDIF
ENDIF
DEALLOCATE(FACNRM,CURRN,SURFN,VOLN)
DEALLOCATE(ELMTAB,ANGTAB)
RETURN
6005 FORMAT(' SALTLS: Global RMS, maximum and average errors (%) ', &
'on region volumes :',3(2X,F10.5))
END SUBROUTINE SALTLS
|
