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