source: palm/trunk/UTIL/surface_output_processing/surface_output_to_vtk.f90 @ 3494

Last change on this file since 3494 was 3494, checked in by suehring, 3 years ago

Surface output revised and some bugs are fixed + new post-processing tool to convert binary surface output to Paraview readable VTK files

  • Property svn:keywords set to Id
File size: 17.2 KB
Line 
1!> @file surface_output_merge.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018  Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: surface_output_to_vtk.f90 3494 2018-11-06 14:51:27Z suehring $
27! Initial version
28!
29! 3241 2018-09-12 15:02:00Z raasch
30!
31!
32! Authors:
33! --------
34! @author Matthias Suehring and Klaus Ketelsen
35!
36!------------------------------------------------------------------------------!
37! Description:
38! ------------
39!> This routine combines surface output from PALM-subdomains into one file.
40!> Output from every processor element is opened and read and output
41!> from all processor elements are written into one file for each timestep
42!> according to Paraview's VTK standard.
43!> Output is distinguished between instantaneous and time-averaged data.
44!------------------------------------------------------------------------------!
45 PROGRAM surface_output_merge
46
47    IMPLICIT NONE
48   
49    CHARACTER(LEN=4)   ::  char_time            !< string indicating simulated time
50    CHARACTER(LEN=4)   ::  file_suffix = '.bin' !< string which contain the suffix indicating surface data
51   
52    CHARACTER(LEN=10)  ::  char_dum             !< dummy string
53   
54    CHARACTER(LEN=30)  ::  myid_char            !< combined string indicating binary file
55   
56    CHARACTER(LEN=100) ::  path                 !< path to the binary data
57    CHARACTER(LEN=100) ::  run                  !< name of the run
58    CHARACTER(LEN=100) ::  variable_name        !< name of the processed output variable   
59   
60    INTEGER(4)   ::  ftell                      !< intrinsic function, get current position in file
61    INTEGER(4)   ::  fseek                      !< intrinsic function, go to given position in file
62    INTEGER(4)   ::  ndum                       !< return parameter of intrinsic function fseek
63   
64    INTEGER, PARAMETER ::  iwp = 4              !< integer precision
65    INTEGER, PARAMETER ::  wp  = 8              !< float precision
66   
67    INTEGER(iwp) ::  cycle_number               !< cycle number
68    INTEGER(iwp) ::  f                          !< running index over all binary files
69    INTEGER(iwp) ::  file_id_in = 110           !< file unit for input binaray file   
70    INTEGER(iwp) ::  file_id_out = 20           !< file unit for output VTK file       
71    INTEGER(iwp) ::  file_id_out_header = 19    !< file unit for temporary header file
72    INTEGER(iwp) ::  length                     !< length of string on file
73    INTEGER(iwp) ::  n                          !< running index over all points and polygons
74    INTEGER(iwp) ::  npoints_total              !< total number of points
75    INTEGER(iwp) ::  ns_total                   !< total number of polygons
76    INTEGER(iwp) ::  num_pe                     !< number of processors of the run
77   
78    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
79    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  npoints !< number of points/vertices in a binaray file
80    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns      !< number of surface elements in a binaray file
81   
82    LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
83    LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
84   
85    REAL(wp)                              ::  simulated_time !< output time   
86   
87    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< actual surface data
88   
89    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points         !< point / vertex data
90    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons       !< polygon data
91
92!
93!-- Read namelist.
94    CALL surface_output_merge_parin
95!
96!-- Allocate array which contains the file position in each output file,
97!-- in order to skip some reading.
98    ALLOCATE( eof( 0:num_pe-1 )     )
99    ALLOCATE( filepos( 0:num_pe-1 ) )
100    ALLOCATE( npoints( 0:num_pe-1 ) )
101    ALLOCATE( ns( 0:num_pe-1 )      )
102!
103!-- Initialize file position.
104    filepos = 0
105!
106!-- Open a temporary file which contains all necessary header information for the
107!-- VTK format.
108    OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
109           FORM = 'FORMATTED' ) 
110!
111!-- READ grid setup, i.e. the number and position of vertices and surface elements
112!-- and merge all this information into one file. Further, create all the required
113!-- header information of the VTK file.
114!-- Note, PARAVIEW expects one VTK file for each variable and each timestep.
115!-- Hence, header information needs to be duplicated multiple times, which will be
116!-- be done later in a bash script.
117!-- Moreover, Paraview expects consecutive vertex and polygon data, which are
118!-- all distributed over the binaray files. Hence, first read vertex data from
119!-- binary file, write this to the HEADER file, close the binary file and read
120!-- data from the next binary file and so on. This requires several openenings
121!-- and closings of the binary files and temporarily storage of the
122!-- file-positions.
123    DO  f = 0, num_pe - 1
124!
125!--    Create filename of the treated binary file.
126       CALL surface_output_merge_create_file_string
127!
128!--    Open file with surface output for processor f.
129       OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //            &
130              TRIM( myid_char ), FORM='UNFORMATTED' )
131!
132!--    Read number of vertices / points and surface elements
133       READ ( file_id_in + f )  npoints(f)
134       READ ( file_id_in + f )  npoints_total
135       READ ( file_id_in + f )  ns(f)
136       READ ( file_id_in + f )  ns_total
137
138!
139!--    Allocate arrays where all the surfaces and vertices will be stored.
140       ALLOCATE( points(3,1:npoints(f))   )
141       
142!
143!--    Read polygon data and store them in a temporary file.
144       READ ( file_id_in + f )  points
145!
146!--    Obtain current file position. Will be stored for next file opening.
147       filepos(f) = ftell( file_id_in + f )
148!
149!--    Write header information. Only one time required.
150       IF ( f == 0 )  THEN
151          WRITE( file_id_out_header,'(A)' ) "# vtk DataFile Version 3.0"
152          WRITE( file_id_out_header,'(A,F8.2,A)' ) "legacy vtk File generated by PALM,  simulation time = xxx sec"
153          WRITE( file_id_out_header,'(A)') "ASCII"
154
155          WRITE( file_id_out_header,'(A)') "DATASET POLYDATA"
156          WRITE( file_id_out_header,'(A,I5,A)') "POINTS ", npoints_total, " float"
157       ENDIF     
158!
159!--    Write the vertex data into header file.
160       DO  n = 1, npoints(f)
161          WRITE( file_id_out_header, '(8F10.1)' )  points(1:3,n)
162       ENDDO
163!
164!--    Deallocate vertex data and close binary file.
165       DEALLOCATE( points )
166       
167       CLOSE ( file_id_in + f )
168    ENDDO
169!
170!-- Now, treat polygon data.
171    DO  f = 0, num_pe - 1
172!
173!--    Create filename of the treated binary file .   
174       CALL surface_output_merge_create_file_string
175!
176!--    Open file with surface output for processor f.
177       OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //            &
178              TRIM( myid_char ), FORM='UNFORMATTED' )
179!
180!--    Move to last postion.
181       ndum = fseek( file_id_in + f, filepos(f), 0 )
182!
183!--    Allocate array for polygon data
184       ALLOCATE( polygons(5,1:ns(f)) )
185!
186!--    Read polygon data and store them in a temporary file.
187       READ ( file_id_in + f )  polygons
188!
189!--    Obtain current file position after reading the local polygon data.
190!--    Will be used for next file opening.
191       filepos(f) = ftell( file_id_in + f )
192!
193!--    Write further header information. Only one time required.
194       IF ( f == 0 )                                                           &
195          WRITE ( file_id_out_header, '(A,8I10)') "POLYGONS ",                 &
196                                                  ns_total, 5 * ns_total
197!
198!--    Write the polygons into the header file. Note, polygon data is of type
199!--    integer, as it just connects the point-indices which describe the given
200!--    surface element.
201       DO n = 1, ns(f)
202          WRITE ( file_id_out_header, '(8I10)' )  INT (polygons(1:5,n) )
203       ENDDO
204!
205!--    Deallocate array for polygon data and close the file.
206       DEALLOCATE( polygons )
207       
208       CLOSE ( file_id_in + f )
209       
210    ENDDO
211   
212    f = 0
213    CALL surface_output_merge_create_file_string
214!
215!-- Write further header information. Only once required.
216    WRITE ( file_id_out_header, '(A,I10)') "CELL_DATA ", ns_total
217    WRITE ( file_id_out_header, '(A,I10)') "SCALARS cell_scalars float 1 "
218    WRITE ( file_id_out_header, '(A,I10)') "LOOKUP_TABLE default "
219!
220!-- Header creation has now been completed. File can be closed.   
221    CLOSE( file_id_out_header )
222!
223!-- Now, read all the actual data from the surface output. Please note, Paraview
224!-- VTK format expects 1 file per variable and time step.
225!-- The output file is only created once and includes the variable and the
226!-- simulated time.
227!-- In the binaray files, several variables and timesteps are stored, data for a
228!-- given variable, however, is distributed over all binary files. Hence, read
229!-- variable data for a given timestep from the binary file, write this data into
230!-- the target output file, remember file position in binary file and close it,
231!-- open nex binary file and read the variable, and so on, until all variables
232!-- for all timesteps are processed.
233    eof = .FALSE.
234    DO
235       DO  f = 0, num_pe - 1
236!
237!--       Clean up strings-
238          char_time = ''
239          variable_name = ''
240!
241!--       Create filename of the treated binary file.           
242          CALL surface_output_merge_create_file_string
243!
244!--       Open binary file with surface output for processor f.
245          OPEN ( file_id_in + f, FILE = TRIM( path ) // TRIM( run ) //         &
246                 TRIM( myid_char ), FORM='UNFORMATTED' )
247!
248!--       Move to last postion.
249          ndum = fseek( file_id_in + f, filepos(f), 0 )
250!
251!--       Read string length and string indicating the output time step.
252          READ ( file_id_in + f ) length
253          READ ( file_id_in + f ) char_time(1:length)
254!
255!--       If string for the output time indicates the end-of-file, set the eof
256!--       flag and skip the read of the loop.
257          IF ( char_time(1:length) == 'END' )  THEN
258             eof(f) = .TRUE. 
259             CLOSE ( file_id_in + f )
260             CYCLE
261          ENDIF
262!
263!--       Read output time, and variable name.
264          READ ( file_id_in + f ) simulated_time
265          READ ( file_id_in + f ) length
266          READ ( file_id_in + f ) variable_name(1:length)
267!
268!--       For first loop index, open the target output file. First create the
269!--       filename string. Further, copy HEADER file with the given filename
270!--       string. The header information must be given in each VTK file!
271          IF ( f == 0 )  THEN
272             IF ( simulated_time < 10.0_wp )  THEN
273                WRITE( char_dum, '(I1)' )  INT( simulated_time ) 
274             ELSEIF ( simulated_time < 100.0_wp )  THEN
275                WRITE( char_dum, '(I2)' )  INT( simulated_time )
276             ELSEIF ( simulated_time < 1000.0_wp )  THEN
277                WRITE( char_dum, '(I3)' )  INT( simulated_time )
278             ELSEIF ( simulated_time < 10000.0_wp )  THEN
279                WRITE( char_dum, '(I4)' )  INT( simulated_time )   
280             ELSEIF ( simulated_time < 100000.0_wp )  THEN
281                WRITE( char_dum, '(I5)' )  INT( simulated_time )   
282             ELSEIF ( simulated_time < 1000000.0_wp )  THEN
283                WRITE( char_dum, '(I6)' )  INT( simulated_time )
284             ELSEIF ( simulated_time < 10000000.0_wp )  THEN
285                WRITE( char_dum, '(I7)' )  INT( simulated_time )   
286             ELSEIF ( simulated_time < 100000000.0_wp )  THEN
287                WRITE( char_dum, '(I8)' )  INT( simulated_time )   
288             ELSEIF ( simulated_time < 1000000000.0_wp )  THEN
289                WRITE( char_dum, '(I9)' )  INT( simulated_time )
290             ENDIF             
291!
292!--          Copy HEADER file
293             CALL system('cp HEADER ' // TRIM( path ) // TRIM( char_dum ) //   & 
294                                      's_' // TRIM( variable_name ) // '.vtk')
295!--          Open VTK file.
296             OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) //    & 
297                    's_' // TRIM( variable_name ) // '.vtk', FORM='FORMATTED', &
298                    POSITION='APPEND' )           
299          ENDIF
300!
301!--       Allocate and read array for variable data. 
302          ALLOCATE( var(1:ns(f)) )
303       
304          READ( file_id_in + f ) var
305!
306!--       Write variable data into output VTK file.
307          DO n = 1, ns(f)
308             WRITE( file_id_out, * ) var(n)
309          ENDDO
310!
311!--       Remember file position in binary file and close it.
312          filepos(f) = ftell( file_id_in + f )
313       
314          CLOSE ( file_id_in + f )
315!
316!--       Deallocate temporary array for variable data.
317          DEALLOCATE( var )
318       
319       ENDDO
320!
321!--    After data for a variable for one time step is read, close the output
322!--    VTK file and go to next variable or timestep.
323       CLOSE ( file_id_out )
324!
325!--    If all files reached the end-of-file, exit the loop.
326       IF ( ALL( eof ) )  EXIT
327       
328    ENDDO
329!
330!-- Finally, remove HEADER file
331    CALL system( 'rm HEADER' )
332   
333 CONTAINS
334 
335!------------------------------------------------------------------------------!
336! Description:
337! ------------
338!> This subroutine read the namelist file.
339!------------------------------------------------------------------------------!
340    SUBROUTINE surface_output_merge_parin
341       
342       IMPLICIT NONE
343       
344       INTEGER(iwp) ::  file_id_parin = 90
345       
346       NAMELIST /surface_output/  convert_average_data, cycle_number, num_pe,  &
347                                  path, run
348
349!
350!--    Open namelist file.
351       OPEN( file_id_parin, FILE='surface_output_parin',                       &
352             STATUS='OLD', FORM='FORMATTED')
353!
354!--    Read namelist.
355       READ ( file_id_parin, surface_output )
356!
357!--    Close namelist file.
358       CLOSE( file_id_parin )
359       
360    END SUBROUTINE surface_output_merge_parin
361     
362!------------------------------------------------------------------------------!
363! Description:
364! ------------
365!> This subroutine creates the filename string of the treated binary file.
366!------------------------------------------------------------------------------!
367    SUBROUTINE surface_output_merge_create_file_string
368       
369       IMPLICIT NONE
370       
371       CHARACTER(LEN=3) ::  char_av = ''
372       CHARACTER(LEN=4) ::  char_cycle = ''
373       
374!
375!--    Create substring for the cycle number.
376       IF ( cycle_number /= 0 )  THEN
377          IF ( cycle_number < 10 )  THEN
378             WRITE( char_cycle, '(I1)')  cycle_number
379             char_cycle = '.00' // TRIM( char_cycle )
380          ELSEIF ( cycle_number < 100 )  THEN
381             WRITE( char_cycle, '(I2)')  cycle_number
382             char_cycle = '.0' // TRIM( char_cycle )
383          ELSEIF ( cycle_number < 1000 )  THEN
384             WRITE( char_cycle, '(I3)')  cycle_number
385             char_cycle = '.' // TRIM( char_cycle )
386          ENDIF
387       ENDIF
388!
389!--    Create substring for averaged data.
390       IF ( convert_average_data )  char_av = '_av'
391!
392!--    Create substring for the processor id and combine all substrings.
393       IF ( f < 10 )  THEN
394          WRITE( char_dum, '(I1)')  f
395          myid_char = TRIM( char_av ) // '_surf_00000' // TRIM( char_dum ) //  &
396                      TRIM( char_cycle ) // file_suffix
397       ELSEIF ( f < 100     )  THEN
398          WRITE( char_dum, '(I2)')  f
399          myid_char = TRIM( char_av ) // '_surf_0000'  // TRIM( char_dum ) //  &
400                      TRIM( char_cycle ) // file_suffix
401       ELSEIF ( f < 1000    )  THEN
402          WRITE( char_dum, '(I3)')  f
403          myid_char = TRIM( char_av ) // '_surf_000'   // TRIM( char_dum ) //  &
404                      TRIM( char_cycle ) // file_suffix
405       ELSEIF ( f < 10000   )  THEN
406          WRITE( char_dum, '(I4)')  f
407          myid_char = TRIM( char_av ) // '_surf_00'    // TRIM( char_dum ) //  &
408                      TRIM( char_cycle ) // file_suffix
409       ELSEIF ( f < 100000  )  THEN
410          WRITE( char_dum, '(I5)')  f
411          myid_char = TRIM( char_av ) // '_surf_0'     // TRIM( char_dum ) //  &
412                      TRIM( char_cycle ) // file_suffix
413       ELSEIF ( f < 1000000 )  THEN
414          WRITE( char_dum, '(I6)')  f
415          myid_char = TRIM( char_av ) // '_surf_'      // TRIM( char_dum ) //  &
416                      TRIM( char_cycle ) // file_suffix
417       ENDIF
418       
419    END SUBROUTINE surface_output_merge_create_file_string
420   
421 END PROGRAM surface_output_merge
422
423
424
Note: See TracBrowser for help on using the repository browser.