1 | include(dom.inc)
2 |
3 | PROGRAM visualensight
4 |
5 | ! ================================================================!
6 | ! !
7 | ! visualensight.F: Creates EnsightGold files for visualisation !
8 | ! and post processing of the PRISSMA outfiles !
9 | ! !
10 | ! author : J. AMAYA (january 2008) !
11 | ! !
12 | ! ================================================================!
13 |
14 | USE avbp_coor
15 |
16 | USE ensightgold_geometry
17 | USE ensight_variables
18 | USE ensightgold_write
19 |
20 | IMPLICIT NONE
21 |
22 | CHARACTER*80 :: outpath, inpath, geometryfile
23 | CHARACTER*80 :: coorf, connf
24 | CHARACTER*80 :: asciiBoundf, prefix, meshfile
25 | CHARACTER*80 :: GFile, SrFile, HFile, QwFile, QrFile, KpFile
26 | CHARACTER*80 :: TFile, EFile, glfile, ndfile
27 | CHARACTER*80 :: cnfile, cffile, epfile, fcfile
28 | CHARACTER*80 :: filespath, ppfile, wlfile, ccfile
29 |
30 | DOM_INT :: myformat, i, nbnodes, l, j, k
31 | DOM_INT :: i_ncells, i_nnodes, i_nfaces, i_nfacemax
32 | DOM_INT :: iface, icell, nfaces, meshtype, readkp, ios
33 |
34 | DOM_REAL :: xc, yc, zc, xfc, yfc, zfc
35 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: cellnodes
36 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: cellfaces
37 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: nodelist
38 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: facelist
39 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: facenodes
40 | DOM_INT, ALLOCATABLE, DIMENSION(:) :: nodecount
41 | DOM_INT, ALLOCATABLE, DIMENSION(:,:) :: patch
42 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: scalar, vector
43 | DOM_REAL, ALLOCATABLE, DIMENSION(:) :: G, Sr, T, P, Kp
44 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: epsil, Tw
45 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: Qw, H, Qr, X
46 |
47 |
48 | ! ------------------!
49 | ! Read choices file !
50 | ! ------------------!
51 |
52 | OPEN(1, FILE="visualensight.choices", FORM='formatted')
53 |
54 | READ(1,*) meshfile
55 | READ(1,*) meshtype
56 | READ(1,*) inpath
57 | READ(1,*) filespath
58 | READ(1,*) outpath
59 | READ(1,*) myformat
60 |
61 | CLOSE(100)
62 |
63 | ! ----------------!
64 | ! Open data files !
65 | ! ----------------!
66 |
67 | GFile = trim(inpath)//'/G.out'
68 | SrFile = trim(inpath)//'/Sr.out'
69 | HFile = trim(inpath)//'/H.out'
70 | TFile = trim(inpath)//'/T.out'
71 | EFile = trim(inpath)//'/E.out'
72 | QwFile = trim(inpath)//'/Qw.out'
73 | QrFile = trim(inpath)//'/Qr.out'
74 | KpFile = trim(inpath)//'/Kp.out'
75 |
76 | glfile = trim(filespath)//'/Global.in'
77 | ndfile = trim(filespath)//'/Nodelist.in'
78 | cffile = trim(filespath)//'/Cell2faces.in'
79 | cnfile = trim(filespath)//'/Cellnodes.in'
80 | epfile = trim(filespath)//'/Emissivities.in'
81 | fcfile = trim(filespath)//'/Facelist.in'
82 | ppfile = trim(filespath)//'/Properties.in'
83 | wlfile = trim(filespath)//'/CLProperties.in'
84 | ccfile = trim(filespath)//'/Centercells.in'
85 |
86 | OPEN(101 ,FILE=GFile ,FORM='UNFORMATTED')
87 | OPEN(102 ,FILE=SrFile,FORM='UNFORMATTED')
88 | OPEN(103 ,FILE=QrFile,FORM='UNFORMATTED')
89 | OPEN(104 ,FILE=HFile ,FORM='UNFORMATTED')
90 | OPEN(105 ,FILE=QwFile,FORM='UNFORMATTED')
91 | OPEN(108 ,FILE=glfile,FORM='UNFORMATTED')
92 | OPEN(109 ,FILE=ndfile,FORM='UNFORMATTED')
93 | OPEN(110 ,FILE=cnfile,FORM='UNFORMATTED')
94 | OPEN(111 ,FILE=epfile,FORM='UNFORMATTED')
95 | OPEN(112 ,FILE=cffile,FORM='UNFORMATTED')
96 | OPEN(113 ,FILE=fcfile,FORM='UNFORMATTED')
97 | OPEN(114 ,FILE=ppfile,FORM='UNFORMATTED')
98 | OPEN(115 ,FILE=wlfile,FORM='UNFORMATTED')
99 | OPEN(117 ,FILE=ccfile,FORM='UNFORMATTED')
100 |
101 | readkp = 1
102 | OPEN(116,FILE=KpFile,IOSTAT=ios,STATUS='old',FORM='UNFORMATTED')
103 | IF (ios.ne.0) readkp = 0
104 |
105 | ! ---------------------!
106 | ! Get global data info !
107 | ! ---------------------!
108 |
109 | READ(108) i_ncells, i_nnodes, i_nfaces, i_nfacemax
110 | CLOSE(108)
111 |
112 | print*, i_ncells, i_nnodes, i_nfaces, i_nfacemax
113 |
114 | ! ----------------------!
115 | ! Read nodal properties !
116 | ! ----------------------!
117 |
118 | ALLOCATE(nodelist(3,i_nnodes))
119 | print*, " nodelist allocated"
120 |
121 | DO i=1, i_nnodes
122 | READ(109) j, (nodelist(k,j),k=1,3)
123 | ENDDO
124 | CLOSE(109)
125 |
126 | ! ---------------------!
127 | ! Read face properties !
128 | ! ---------------------!
129 |
130 | ALLOCATE(facenodes(6, i_nfaces))
131 | facenodes = 0
132 | print*, " facenodes allocated"
133 |
134 | DO i=1, i_nfaces
135 | READ(113) iface, nbnodes, (facenodes(j,i),j=2,nbnodes+1)
136 | facenodes(1,i) = nbnodes
137 | ENDDO
138 | CLOSE(113)
139 |
140 | ! -------------------------!
141 | ! Read nodal and cell data !
142 | ! -------------------------!
143 | ALLOCATE(nodecount(i_nnodes))
144 | print*, " nodecount allocated"
145 | ALLOCATE(cellnodes(9,i_ncells))
146 | print*, " cellnodes allocated"
147 | ALLOCATE(cellfaces(i_nfacemax+1,i_ncells))
148 | print*, " cellfaces allocated"
149 | ALLOCATE(epsil(i_nfacemax,i_ncells))
150 | print*, " epsil allocated"
151 | ALLOCATE(Tw(i_nfacemax,i_ncells))
152 | print*, " Tw allocated"
153 | ALLOCATE(patch(i_nfacemax,i_ncells))
154 | print*, " Patch allocated"
155 | ALLOCATE(Sr(i_ncells))
156 | print*, " Sr allocated"
157 | ALLOCATE(T(i_ncells))
158 | print*, " T allocated"
159 | ALLOCATE(P(i_ncells))
160 | print*, " P allocated"
161 | ALLOCATE(X(6,i_ncells))
162 | print*, " X allocated"
163 | ALLOCATE(G (i_ncells))
164 | print*, " G allocated"
165 | ALLOCATE(Qr(3,i_ncells))
166 | print*, " Qr allocated"
167 | ALLOCATE(H(i_nfacemax,i_ncells))
168 | print*, " H allocated"
169 | ALLOCATE(Qw(i_nfacemax,i_ncells))
170 | print*, " Qw allocated"
171 | ALLOCATE(Kp(i_ncells))
172 | print*, " Kp allocated"
173 |
174 | Sr = 0.
175 | G = 0.
176 | Qr = 0.
177 | H = 0.
178 | Qw = 0.
179 | Kp = 0.
180 | T = 0.
181 | P = 0.
182 | X = 0.
183 |
184 | cellnodes = 0
185 | epsil = -1.
186 |
187 | DO i=1,i_ncells
188 |
189 | ! -----------------------------------------!
190 | ! Read connections between cells and nodes !
191 | ! -----------------------------------------!
192 | READ(110) icell, nbnodes, (cellnodes(j,icell),j=2,nbnodes+1)
193 | cellnodes(1,icell) = nbnodes
194 |
195 | ! ---------------------!
196 | ! Read this cell faces !
197 | ! ---------------------!
198 | READ(111) icell, nfaces, (epsil(j,icell),j=1,nfaces)
199 | READ(112) icell, nfaces, (cellfaces(j,icell),j=2,nfaces+1)
200 | READ(115) icell, nfaces, (Tw(j,icell),patch(j,icell), &
201 | & j=1,nfaces)
202 | cellfaces(1,icell) = nfaces
203 |
204 | ! print*, " cell: ", icell
205 | ! print*, Tw(:,i)
206 | ! print*, epsil(:,i)
207 | ! print*, patch(:,i)
208 | ! read(*,*)
209 |
210 | ! ---------------------------!
211 | ! Read PRISSMA initial field !
212 | ! ---------------------------!
213 |
214 | READ(114) icell, T(icell), P(icell), (X(j,icell),j=1,6)
215 |
216 | ! --------------------------!
217 | ! Read PRISSMA output files !
218 | ! --------------------------!
219 | READ(117) icell, xc, yc, zc
220 | READ(101) icell, G(icell)
221 | READ(102) icell, Sr(icell)
222 | IF (readkp.eq.1) READ(116) icell, Kp(icell)
223 | READ(103) icell, (Qr(j,icell),j=1,3)
224 |
225 | DO j=1,nfaces
226 | IF (epsil(j,i).ne.-1.) THEN
227 | ! print*, i,"/",i_ncells,":",j,"/",nfaces,":",epsil(j,i)
228 | READ(104) iface, icell, H(iface,icell)
229 | READ(105) iface, icell, Qw(iface,icell)
230 | ENDIF
231 | ENDDO
232 |
233 | ENDDO
234 |
235 | CLOSE(101)
236 | CLOSE(102)
237 | CLOSE(103)
238 | CLOSE(104)
239 | CLOSE(105)
240 | CLOSE(110)
241 | CLOSE(111)
242 | CLOSE(112)
243 | CLOSE(114)
244 | CLOSE(115)
245 | CLOSE(116)
246 | CLOSE(117)
247 |
248 | ! -----------------!
249 | ! Setup parameters !
250 | ! -----------------!
251 |
252 | ensight_overwrite = .TRUE.
253 | geometryfile = trim(outpath)//'/ensight.geo'
254 |
255 | ensight_write_format = myformat
256 |
257 | CALL intialisegeometry
258 | print*, " Geometry initialized"
259 |
260 | IF (meshtype.eq.1) THEN
261 |
262 | coorf = trim(meshfile)//'.coor'
263 | connf = trim(meshfile)//'.conn'
264 | asciiBoundf = trim(meshfile)//'.asciiBound'
265 |
266 | CALL creategeometry(coorf, connf, asciiBoundf, geometryfile)
267 | print*, " Geometry created"
268 |
269 | CALL getcoorinfo
270 |
271 | ensight_nnode = coor_nnode
272 | ensight_vectordims = coor_ndim
273 |
274 | ELSEIF (meshtype.eq.2) THEN
275 |
276 | CALL creategeometry_gambit(meshfile, geometryfile, &
277 | & ensight_nnode,ensight_vectordims)
278 | print*, " Geometry created"
279 |
280 | ENDIF
281 |
282 | ensight_nrec = 1
283 | ensight_nscalars = 15
284 | ensight_nvectors = 1
285 |
286 | ! ------------------------------------!
287 | ! Allocate memory and initialize data !
288 | ! ------------------------------------!
289 |
290 | CALL create_ensight_variables
291 | CALL create_ensight_info
292 | CALL create_ensight_steprec
293 |
294 | print*, " Ensight data initialized"
295 |
296 | ALLOCATE(scalar(ensight_nscalars, ensight_nnode))
297 | ALLOCATE(vector(ensight_nvectors,ensight_nnode* &
298 | & ensight_vectordims))
299 | print*, " Ensight vectors allocated"
300 |
301 | ! ----------------------------!
302 | ! Scatter values to the nodes !
303 | ! ----------------------------!
304 |
305 | scalar = 0.
306 | vector = 0.
307 |
308 | nodecount = 0
309 |
310 | DO j=1, i_ncells
311 | DO i=1, cellnodes(1,j)
312 | k = cellnodes(i+1,j)
313 | nodecount(k) = nodecount(k) + 1
314 | scalar(1,k) = scalar(1,k) + G(j)
315 | scalar(2,k) = scalar(2,k) + Sr(j)
316 | scalar(5,k) = scalar(5,k) + T(j)
317 | scalar(6,k) = scalar(6,k) + X(1,j)
318 | scalar(7,k) = scalar(7,k) + X(2,j)
319 | scalar(8,k) = scalar(8,k) + X(3,j)
320 | scalar(9,k) = scalar(9,k) + X(4,j)
321 | scalar(10,k) = scalar(10,k) + X(5,j)
322 | scalar(11,k) = scalar(11,k) + X(6,j)
323 | scalar(12,k) = scalar(12,k) + Kp(j)
324 |
325 | vector(1,k) = vector(1,k) + Qr(1,j)
326 | vector(1,k+i_nnodes) = vector(1,k+i_nnodes) + Qr(2,j)
327 |
328 | IF (ensight_vectordims.eq.3) THEN
329 | vector(1,k+2*i_nnodes) = vector(1,k+2*i_nnodes) + Qr(3,j)
330 | ENDIF
331 |
332 | ENDDO
333 | ENDDO
334 | print*, " Nodes counted"
335 |
336 | ! --------------------!
337 | ! Divide by nodecount !
338 | ! --------------------!
339 |
340 | scalar(1,:) = scalar(1,:) / real(nodecount(:))
341 | scalar(2,:) = scalar(2,:) / real(nodecount(:))
342 | scalar(5,:) = scalar(5,:) / real(nodecount(:))
343 | scalar(12,:) = scalar(12,:) / real(nodecount(:))
344 |
345 | DO i=6,11
346 | scalar(i,:) = scalar(i,:) / real(nodecount(:))
347 | ENDDO
348 |
349 | print*, " Scalars scattered"
350 |
351 | vector(1,1:i_nnodes) = &
352 | & vector(1,1:i_nnodes) / real(nodecount(:))
353 | vector(1,i_nnodes+1:2*i_nnodes) = &
354 | & vector(1,i_nnodes+1:2*i_nnodes) / real(nodecount(:))
355 | IF (ensight_vectordims.eq.3) THEN
356 | vector(1,2*i_nnodes+1:3*i_nnodes)= &
357 | & vector(1,2*i_nnodes+1:3*i_nnodes) / real(nodecount(:))
358 | ENDIF
359 |
360 | print*, " Vectors scattered"
361 |
362 | ! -------------------------------------!
363 | ! Scatter values to the boundary nodes !
364 | ! -------------------------------------!
365 |
366 | nodecount = 0
367 |
368 | DO j=1, i_ncells
369 | DO i=1, cellfaces(1,j)
370 |
371 | ! print*, " cellfaces:",cellfaces(:,j)
372 | iface = cellfaces(i+1,j)
373 | ! print*, " iface:", iface
374 |
375 | IF (epsil(i,j).ne.-1.) THEN
376 |
377 | ! print*, " epsil=",epsil(i,j)
378 | ! print*, " Tw=",Tw(i,j)
379 | DO k=1, facenodes(1,iface)
380 | l = facenodes(k+1,iface)
381 | nodecount(l) = nodecount(l) + 1
382 | scalar(3,l) = scalar(3,l) + H(i,j)
383 | scalar(4,l) = scalar(4,l) + Qw(i,j)
384 | scalar(13,l) = scalar(13,l) + epsil(i,j)
385 | scalar(14,l) = scalar(14,l) + Tw(i,j)
386 | scalar(15,l) = scalar(15,l) + patch(i,j)
387 | ENDDO
388 |
389 | ENDIF
390 |
391 | ENDDO
392 | ENDDO
393 |
394 | ! ------------------------!
395 | ! Avoid divisions by zero !
396 | ! ------------------------!
397 |
398 | DO i=1,i_nnodes
399 | IF (nodecount(i).eq.0.) nodecount(i) = 1
400 | ENDDO
401 |
402 | ! --------------------!
403 | ! Divide by nodecount !
404 | ! --------------------!
405 |
406 | scalar(3,:) = scalar(3,:) / real(nodecount(:))
407 | scalar(4,:) = scalar(4,:) / real(nodecount(:))
408 | scalar(13,:) = scalar(13,:) / real(nodecount(:))
409 | scalar(14,:) = scalar(14,:) / real(nodecount(:))
410 | scalar(15,:) = scalar(15,:) / real(nodecount(:))
411 |
412 | print*, " Boundary scalars scattered"
413 |
414 | ! -------------------------------!
415 | ! Fill ensight vectors with data !
416 | ! -------------------------------!
417 |
418 | print*, " Filling scalar names"
419 | ensight_scal_name(1) = "Incident_radiation_G"
420 | ensight_scal_name(2) = "Radiation_source_term_Sr"
421 | ensight_scal_name(3) = "Wall_incident_radiation_H"
422 | ensight_scal_name(4) = "Wall_radiative_heatflux_Qw"
423 | ensight_scal_name(5) = "Temperature"
424 | ensight_scal_name(6) = "XH2O"
425 | ensight_scal_name(7) = "XCO2"
426 | ensight_scal_name(8) = "XCO"
427 | ensight_scal_name(9) = "XO2"
428 | ensight_scal_name(10)= "XN2"
429 | ensight_scal_name(11)= "XSOOT"
430 | ensight_scal_name(12)= "Mean_Planck_Absorption"
431 | ensight_scal_name(13)= "Wall_emitance_epsilon"
432 | ensight_scal_name(14)= "Wall_temperature_Tw"
433 | ensight_scal_name(15)= "Patch"
434 |
435 | print*, " Filling vector name"
436 | ensight_vec_name(1) = "Radiated_heat_flux"
437 |
438 | print*, " Filling file names"
439 | ensight_scal_fn(1) = ensight_scal_name(1)
440 | ensight_scal_fn(2) = ensight_scal_name(2)
441 | ensight_scal_fn(3) = ensight_scal_name(3)
442 | ensight_scal_fn(4) = ensight_scal_name(4)
443 | ensight_scal_fn(5) = ensight_scal_name(5)
444 | ensight_scal_fn(6) = ensight_scal_name(6)
445 | ensight_scal_fn(7) = ensight_scal_name(7)
446 | ensight_scal_fn(8) = ensight_scal_name(8)
447 | ensight_scal_fn(9) = ensight_scal_name(9)
448 | ensight_scal_fn(10)= ensight_scal_name(10)
449 | ensight_scal_fn(11)= ensight_scal_name(11)
450 | ensight_scal_fn(12)= ensight_scal_name(12)
451 | ensight_scal_fn(13)= ensight_scal_name(13)
452 | ensight_scal_fn(14)= ensight_scal_name(14)
453 | ensight_scal_fn(15)= ensight_scal_name(15)
454 |
455 | ensight_vec_fn(1) = ensight_vec_name(1)
456 |
457 | ensight_scalars = scalar
458 | ensight_vectors = vector
459 |
460 | print*, " -------------------"
461 | print*, " MAX Wall_Incid_rad = ", MAXVAL(scalar(3,:))
462 | print*, " MIN Wall_Incid_rad = ", MINVAL(scalar(3,:))
463 | print*, " -------------------"
464 | print*, " MAX Heat flux Qw = ", MAXVAL(scalar(4,:))
465 | print*, " MIN Heat flux Qw = ", MINVAL(scalar(4,:))
466 | print*, " -------------------"
467 | print*, " MAX Heat source Sr = ", MAXVAL(scalar(2,:))
468 | print*, " MIN Heat source Sr = ", MINVAL(scalar(2,:))
469 | print*, " -------------------"
470 | print*, " MAX temperature T = ", MAXVAL(scalar(5,:))
471 | print*, " MIN temperature T = ", MINVAL(scalar(5,:))
472 | print*, " -------------------"
473 |
474 | ! -------------------!
475 | ! Write ensight data !
476 | ! -------------------!
477 |
478 | prefix = 'none'
479 | print*, " Calling ensight_write_sol"
480 | CALL ensight_write_sol(outpath," ")
481 | print*, " Calling ensight_write_case"
482 | CALL ensight_write_case(outpath, prefix, 1, 1, 1, 1, .FALSE.)
483 |
484 | END PROGRAM visualensight