1 | include(dom.inc)
2 |
3 | #include "seq.h"
4 |
5 | PROGRAM quad_tracing
6 |
7 | ! ==================================================================!
8 | ! !
9 | ! Find the direction of the quadrature !
10 | ! !
11 | ! ==================================================================!
12 |
13 |
14 |
15 | ! use OMP_LIB
16 | MODULE_OMP
17 | IMPLICIT NONE
18 | include 'dom_constants.h'
19 |
20 | DOM_INT, PARAMETER :: FILE_TRACE = 200
21 |
22 | CHARACTER*80 :: path1
23 | CHARACTER*80 :: path2
24 | CHARACTER*80 :: gdatafile
25 | CHARACTER*80 :: c2cfile
26 | CHARACTER*80 :: normfile
27 | CHARACTER*80 :: nodesfile
28 | CHARACTER*80 :: quadrfile
29 | CHARACTER*80 :: c2facfile
30 | CHARACTER*80 :: tracefile
31 | CHARACTER*80 :: cellsfile
32 |
33 | DOM_INT :: numthreads,nthread
34 | DOM_INT :: ntcells, ntnodes, ntfaces,ntfacesmax,ndirs
35 | DOM_INT :: cellnodes
36 | DOM_INT :: icell,i,j,l
37 | DOM_INT :: inear,nodenear
38 | DOM_INT :: idir,unique_cell
39 | DOM_INT :: next_cell,outface
40 | DOM_INT, DIMENSION(10) :: inode
41 | DOM_INT, ALLOCATABLE, DIMENSION(:) :: nfaces
42 | DOM_INT, ALLOCATABLE, DIMENSION(:,:,:) :: neighb
43 |
44 | DOM_REAL :: projection_2">max_projection,projection
45 | DOM_REAL :: xmin,xmax,ymin,ymax,zmin,zmax,xn,yn,zn,xp,yp,zp
46 | DOM_REAL :: r,rmin
47 | DOM_REAL, ALLOCATABLE, DIMENSION(:) :: mu,eta,ksi
48 | DOM_REAL, ALLOCATABLE, DIMENSION(:) :: trace
49 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: trace_priv
50 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: nx,ny,nz,coor
51 |
52 |
53 | LOGICAL :: onecell = .FALSE.
54 | LOGICAL :: first_cell = .TRUE.
55 | ! ------------------!
56 | ! Read choices file !
57 | ! ------------------!
58 |
59 | OPEN (FILE_CHCS , FILE='quad_tracing.choices', FORM='FORMATTED')
60 | READ (FILE_CHCS,*) path1 ! Path for infiles *.in
61 | READ (FILE_CHCS,*) path2 ! Path for infiles *.out
62 | READ (FILE_CHCS,*) unique_cell,xp,yp,zp ! Read the node number : if 0 then, tracking is done for all the domain
63 | if (unique_cell .ne. 0) then
64 | onecell = .TRUE.
65 | print*,'Calculation performed on one point'
66 | else
67 | print*,'Calculation performed on all domains'
68 | endif
69 |
70 | CLOSE(FILE_CHCS)
71 |
72 | ! ---------------!
73 | ! Set file names !
74 | ! ---------------!
75 |
76 | gdatafile = path1(1:len_trim(path1))//'/Global.in'
77 | c2cfile = path1(1:len_trim(path1))//'/Cell2cells.in'
78 | normfile = path1(1:len_trim(path1))//'/Normals.in'
79 | nodesfile = path1(1:len_trim(path1))//'/Nodelist.in'
80 | quadrfile = path1(1:len_trim(path1))//'/Quadrature.in'
81 | c2facfile = path1(1:len_trim(path1))//'/Cell2faces.in'
82 | cellsfile = path1(1:len_trim(path1))//'/Cellnodes.in'
83 | tracefile = path2(1:len_trim(path2))//'/Trace.out'
84 | ! ------------------------------------!
85 | ! Open binary/ascii Global data files !
86 | ! ------------------------------------!
87 |
88 | OPEN(FILE_GDATA, FILE=gdatafile , FORM='UNFORMATTED')
89 |
90 | ! ------------------------!
91 | ! Open binary input files !
92 | ! ------------------------!
93 |
94 | OPEN(FILE_C2C, FILE=c2cfile, FORM='unformatted')
95 | OPEN(FILE_NORM, FILE=normfile, FORM='unformatted')
96 | OPEN(FILE_QUADR, FILE=quadrfile, FORM='unformatted')
97 | OPEN(FILE_NODES, FILE=nodesfile, FORM='unformatted')
98 | OPEN(FILE_CELLS, FILE=cellsfile, FORM='unformatted')
99 | OPEN(FILE_TRACE, FILE=tracefile, FORM='unformatted')
100 |
101 | ! ----------------------!
102 | ! Transform Global data !
103 | ! ----------------------!
104 |
105 | print*,'Read datas'
106 | READ(FILE_GDATA) ntcells, ntnodes, ntfaces,ntfacesmax,
107 | & ndirs
108 | CLOSE(FILE_GDATA)
109 | ! -----------------!
110 | ! Allocate vectors !
111 | ! -----------------!
112 |
113 | ALLOCATE(neighb(ntcells,ntfacesmax,2))
114 | ALLOCATE(nx(ntcells,ntfacesmax))
115 | ALLOCATE(ny(ntcells,ntfacesmax))
116 | ALLOCATE(nz(ntcells,ntfacesmax))
117 | ALLOCATE(nfaces(ntcells))
118 | ALLOCATE(mu(ndirs))
119 | ALLOCATE(eta(ndirs))
120 | ALLOCATE(ksi(ndirs))
121 |
122 | OMP_MAX_THREADS
123 | ! numthreads=omp_get_max_threads()
124 | ! numthreads=1!omp_get_max_threads()
125 | ALLOCATE(trace_priv(ntcells,numthreads))
126 | ALLOCATE(trace(ntcells))
127 | ALLOCATE(coor(3,ntnodes))
128 |
129 | print*,'Vectors sucessfully allocated'
130 | ! ----------------------!
131 | ! Read the datas !
132 | ! ----------------------!
133 | DO icell=1,ntcells
134 | READ(FILE_C2C) i,nfaces(i), &
135 | & (neighb(i,l,1),neighb(i,l,2),l=1,nfaces(i))
136 |
137 | READ(FILE_NORM) i,nfaces(i),(nx(i,l),ny(i,l), &
138 | & nz(i,l),l=1,nfaces(i))
139 |
140 | ENDDO
141 |
142 | DO j=1,ntnodes
143 | READ(FILE_NODES) i,(coor(l,i),l=1,3)
144 | ENDDO
145 |
146 | DO j=1,ndirs
147 | READ(FILE_QUADR) mu(j),eta(j),ksi(j)
148 | ENDDO
149 |
150 | ! ----------------------!
151 | ! Track point if unique !
152 | ! ----------------------!
153 |
154 | if (onecell) then
155 | print*,'Track the nearest node'
156 |
157 | xmin = 10e12
158 | xmax = -10e12
159 | ymin = 10e12
160 | ymax = -10e12
161 | zmin = 10e12
162 | zmax = -10e12
163 | inear = -1
164 | rmin = 10e12
165 |
166 | do i=1,ntnodes
167 |
168 | xn=coor(1,i)
169 | yn=coor(2,i)
170 | zn=coor(3,i)
171 | xmin=min(xmin,xn)
172 | xmax=max(xmax,xn)
173 | ymin=min(ymin,yn)
174 | ymax=max(ymax,yn)
175 | zmin=min(zmin,zn)
176 | zmax=max(zmax,zn)
177 | r=dsqrt((xp-xn)**2+(yp-yn)**2+(zp-zn)**2)
178 | if ( r.lt.rmin ) then
179 | rmin = r
180 | inear = i
181 | endif
182 | enddo
183 | endif
184 | ! Search one cell with the node number 'inear'
185 |
186 | do icell=1,ntcells
187 | READ(FILE_CELLS) i,cellnodes,(inode(l),l=1,cellnodes)
188 | do l=1,cellnodes
189 | if (inode(l).eq.inear) then
190 | nodenear = i
191 | exit
192 | endif
193 | enddo
194 | enddo
195 |
196 | ! --------------------!
197 | ! Loop over the cells !
198 | ! --------------------!
199 | trace_priv(:,:) = 0
200 | print*,'Perform the calculation'
201 |
202 | c$OMP PARALLEL DO
203 | c$OMP& SHARED(ndirs,ntcells,nfaces,nx,ny,nz,mu,eta,ksi,neighb)
204 | c$OMP& SHARED(onecell,nodenear,trace_priv)
205 | c$OMP& PRIVATE(first_cell,icell,idir,next_cell,max_projection,i)
206 | c$OMP& PRIVATE(j,projection,nthread,outface)
207 |
208 | do icell =1,ntcells
209 | OMP_THREAD
210 | ! nthread=omp_get_thread_num()+1
211 | if ((onecell) .and. (icell .ne. nodenear)) cycle
212 | do idir=1,ndirs ! Over all the directions
213 | next_cell = icell
214 | do while(next_cell .ne.0) ! if nex_cell = 0, it's a BC
215 | max_projection = 0
216 | i = next_cell
217 | if(first_cell) then
218 | first_cell = .FALSE.
219 | else
220 | trace_priv(next_cell,nthread) = 1 + &
221 | & trace_priv(next_cell,nthread)
222 | endif
223 | outface = 0
224 |
225 | do j=1,nfaces(i)
226 | projection = nx(i,j)*mu(idir) + ny(i,j)* &
227 | & eta(idir) + nz(i,j)* ksi(idir) !normal(j) * direction(idir)
228 |
229 | if (projection .ge. max_projection) then
230 | outface = j ! Exit face
231 | projection">max_projection = projection
232 | endif
233 |
234 | enddo
235 |
236 | next_cell = neighb(i,outface,1)
237 | enddo
238 | first_cell = .TRUE.
239 | enddo
240 | enddo
241 | !$OMP END PARALLEL DO
242 |
243 | ! If OMP is used we need to make a reduction over all threads
244 | trace(:)=0
245 | do i=1,numthreads
246 | trace(:)= trace_priv(:,i) + trace(:)
247 | enddo
248 |
249 | ! ----------------!
250 | ! Write Trace.out !
251 | ! ----------------!
252 |
253 | do icell=1,ntcells
254 | WRITE(FILE_TRACE) icell,trace(icell)
255 | enddo
256 | print*,'Solution sucessfully written in',TraceFile
257 | ! ----------------!
258 | ! Close all files !
259 | ! ----------------!
260 |
261 | CLOSE(FILE_C2C)
262 | CLOSE(FILE_NORM)
263 | CLOSE(FILE_QUADR)
264 | CLOSE(FILE_NODES)
265 | CLOSE(FILE_CELLS)
266 | CLOSE(FILE_TRACE)
267 |
268 | ! -------------------!
269 | ! Deallocate vectors !
270 | ! -------------------!
271 |
272 | DEALLOCATE(neighb)
273 | DEALLOCATE(nx)
274 | DEALLOCATE(ny)
275 | DEALLOCATE(nz)
276 | DEALLOCATE(nfaces)
277 | DEALLOCATE(mu)
278 | DEALLOCATE(eta)
279 | DEALLOCATE(ksi)
280 | DEALLOCATE(trace)
281 |
282 | END PROGRAM quad_tracing
quad_tracing.F could be called by: