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