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

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

Change format descriptor in VTK postprocessor for surface output; switch back to non-standard Fortran functions since standard-conform C-functions lead to unexpected segmentation faults for some unknown reason

  • Property svn:keywords set to Id
File size: 20.7 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 3755 2019-02-19 17:03:19Z eckhard $
27! - Change format description to avoid *** for larger domains
28! - switch back from to non-standard Fortran ftell and fseek since posix
29!   interface do not work any more for some unknown reason
30!
31! 3523 2018-11-13 16:09:31Z suehring
32! Implement interface for posix conform C-systemcalls in order to replace
33! non-standard FORTRAN functions ftell and fseek.
34!
35! 3496 2018-11-06 15:59:50Z suehring
36! Use subroutine call for fseek instead of function call. gfortran has some
37! problems with this.
38!
39! 3494 2018-11-06 14:51:27Z suehring
40! Initial version
41!
42! Authors:
43! --------
44! @author Matthias Suehring and Klaus Ketelsen
45!
46!------------------------------------------------------------------------------!
47! Description:
48! ------------
49!> Interface to posix systemcalls. The module allows to call the following
50!> C system calls from FORTRAN: ftell, lseek.
51!------------------------------------------------------------------------------!
52 MODULE posix_interface
53 
54    INTERFACE
55       INTEGER (C_SIZE_T) FUNCTION C_LSEEK (fd, offset, whence)                &
56                                   BIND(C, NAME='lseek')
57         
58          USE ISO_C_BINDING
59         
60          IMPLICIT NONE
61         
62          INTEGER(KIND=C_INT),    VALUE ::  fd
63          INTEGER(KIND=C_SIZE_T), VALUE ::  offset
64          INTEGER(KIND=C_INT),    VALUE ::  whence
65         
66       END FUNCTION C_LSEEK
67       
68    END INTERFACE
69   
70    INTERFACE
71       INTEGER (C_SIZE_T) FUNCTION C_FTELL ( fd )                              &
72                                   BIND(C, NAME='ftell')
73         
74          USE ISO_C_BINDING
75         
76          IMPLICIT NONE
77         
78          INTEGER(KIND=C_INT),    VALUE ::  fd
79
80       END FUNCTION C_FTELL
81       
82    END INTERFACE
83   
84    PUBLIC posix_lseek, posix_ftell
85   
86    CONTAINS
87!
88!------------------------------------------------------------------------------!
89! Description:
90! ------------
91!> Interface for the C-routine lseek.
92!------------------------------------------------------------------------------!
93    SUBROUTINE posix_lseek( fid, offset )
94   
95       USE ISO_C_BINDING
96   
97       IMPLICIT NONE
98
99       INTEGER,INTENT(IN)                     ::  fid    !< file unit
100       INTEGER(KIND=C_SIZE_T),INTENT(IN)      ::  offset !< file offset from the beginning
101
102       INTEGER(KIND=C_INT)                    ::  my_fid !< file unit
103       INTEGER(KIND=C_SIZE_T)                 ::  retval !< return value
104       INTEGER(KIND=C_INT)                    ::  whence !< mode, here start from the beginning
105
106       my_fid = fid
107       whence = 0
108
109       retval = C_LSEEK( fid, offset, whence )
110       
111    END SUBROUTINE posix_lseek
112!
113!------------------------------------------------------------------------------!
114! Description:
115! ------------
116!> Interface for the C-routine ftell.
117!------------------------------------------------------------------------------!
118    SUBROUTINE posix_ftell( fid, filepos )
119   
120       USE ISO_C_BINDING
121   
122       IMPLICIT NONE
123
124       INTEGER,INTENT(IN)      ::  fid                     !< file unit
125       INTEGER(KIND=C_SIZE_T), INTENT(INOUT)   ::  filepos !< file position
126
127       filepos = C_FTELL( fid )
128       
129    END SUBROUTINE posix_ftell
130
131 END MODULE posix_interface
132!
133!------------------------------------------------------------------------------!
134! Description:
135! ------------
136!> This routine combines surface output from PALM-subdomains into one file.
137!> Output from every processor element is opened and read and output
138!> from all processor elements are written into one file for each timestep
139!> according to Paraview's VTK standard.
140!> Output is distinguished between instantaneous and time-averaged data.
141!------------------------------------------------------------------------------!
142 PROGRAM surface_output_to_vtk
143 
144    USE, INTRINSIC ::  ISO_C_BINDING
145 
146    USE posix_interface,                                                       &
147        ONLY:  posix_ftell, posix_lseek
148
149    IMPLICIT NONE
150   
151    CHARACTER(LEN=4)   ::  char_time              !< string indicating simulated time
152    CHARACTER(LEN=4)   ::  file_suffix = '.bin'   !< string which contain the suffix indicating surface data
153   
154    CHARACTER(LEN=10)  ::  char_dum               !< dummy string
155   
156    CHARACTER(LEN=30)  ::  myid_char              !< combined string indicating binary file
157   
158
159    CHARACTER(LEN=100) ::  path                   !< path to the binary data
160    CHARACTER(LEN=100) ::  run                    !< name of the run
161    CHARACTER(LEN=100) ::  variable_name          !< name of the processed output variable   
162
163    INTEGER(4)   ::  ftell                      !< intrinsic function, get current position in file
164    INTEGER(4)   ::  ndum                       !< return parameter of intrinsic function fseek   
165       
166    INTEGER, PARAMETER ::  iwp = 4                !< integer precision
167    INTEGER, PARAMETER ::  wp  = 8                !< float precision
168    INTEGER, PARAMETER ::  OFFSET_KIND = C_SIZE_T !< unsigned integer for the C-interface
169
170    INTEGER(iwp) ::  cycle_number                 !< cycle number
171    INTEGER(iwp) ::  f                            !< running index over all binary files
172    INTEGER(iwp) ::  file_id_in = 18              !< file unit for input binaray file   
173    INTEGER(iwp) ::  file_id_out = 20             !< file unit for output VTK file       
174    INTEGER(iwp) ::  file_id_out_header = 19      !< file unit for temporary header file
175    INTEGER(iwp) ::  length                       !< length of string on file
176    INTEGER(iwp) ::  n                            !< running index over all points and polygons
177    INTEGER(iwp) ::  npoints_total                !< total number of points
178    INTEGER(iwp) ::  ns_total                     !< total number of polygons
179    INTEGER(iwp) ::  num_pe                       !< number of processors of the run
180
181!     INTEGER(OFFSET_KIND),DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
182    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  filepos !< current fileposition in binary file
183   
184    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  npoints !< number of points/vertices in a binaray file
185    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns      !< number of surface elements in a binaray file
186   
187    LOGICAL ::  convert_average_data = .FALSE.          !< namelist parameter to decide whether average or instantaneous data should be converted
188    LOGICAL, DIMENSION(:), ALLOCATABLE      ::  eof     !< flag indicating that end of binary file is reached
189   
190    REAL(wp)                              ::  simulated_time !< output time   
191   
192    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< actual surface data
193   
194    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points         !< point / vertex data
195    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons       !< polygon data
196   
197    logical :: flag
198
199!
200!-- Read namelist.
201    CALL surface_output_parin
202!
203!-- Allocate array which contains the file position in each output file,
204!-- in order to skip some reading.
205    ALLOCATE( eof(0:num_pe-1)     )
206    ALLOCATE( filepos(0:num_pe-1) )
207    ALLOCATE( npoints(0:num_pe-1) )
208    ALLOCATE( ns(0:num_pe-1)      )
209!
210!-- Initialize file position.
211    filepos = 0
212!
213!-- Open a temporary file which contains all necessary header information for the
214!-- VTK format.
215    OPEN ( file_id_out_header, FILE = 'HEADER', STATUS = 'REPLACE',            &
216           FORM = 'FORMATTED' ) 
217!
218!-- READ grid setup, i.e. the number and position of vertices and surface elements
219!-- and merge all this information into one file. Further, create all the required
220!-- header information of the VTK file.
221!-- Note, PARAVIEW expects one VTK file for each variable and each timestep.
222!-- Hence, header information needs to be duplicated multiple times, which will be
223!-- be done later in a bash script.
224!-- Moreover, Paraview expects consecutive vertex and polygon data, which are
225!-- all distributed over the binaray files. Hence, first read vertex data from
226!-- binary file, write this to the HEADER file, close the binary file and read
227!-- data from the next binary file and so on. This requires several openenings
228!-- and closings of the binary files and temporarily storage of the
229!-- file-positions.
230    DO  f = 0, num_pe - 1
231!
232!--    Create filename of the treated binary file.
233       CALL surface_output_create_file_string
234!
235!--    Open file with surface output for processor f.
236       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
237              TRIM( myid_char ), FORM = 'UNFORMATTED' )
238!
239!--    Read number of vertices / points and surface elements
240       READ ( file_id_in )  npoints(f)
241       READ ( file_id_in )  npoints_total
242       READ ( file_id_in )  ns(f)
243       READ ( file_id_in )  ns_total
244
245!
246!--    Allocate arrays where all the surfaces and vertices will be stored.
247       ALLOCATE( points(3,1:npoints(f))   )
248       
249!
250!--    Read polygon data and store them in a temporary file.
251       READ ( file_id_in )  points
252!
253!--    Obtain current file position. Will be stored for next file opening.
254!        CALL posix_ftell( file_id_in, filepos(f) )
255       filepos(f) = ftell( file_id_in )
256!
257!--    Write header information. Only one time required.
258       IF ( f == 0 )  THEN
259          WRITE( file_id_out_header,'(A)' ) "# vtk DataFile Version 3.0"
260          WRITE( file_id_out_header,'(A,F8.2,A)' ) "legacy vtk File generated by PALM,  simulation time = xxx sec"
261          WRITE( file_id_out_header,'(A)') "ASCII"
262
263          WRITE( file_id_out_header,'(A)') "DATASET POLYDATA"
264          WRITE( file_id_out_header,'(A,I5,A)') "POINTS ", npoints_total, " float"
265       ENDIF     
266!
267!--    Write the vertex data into header file.
268       DO  n = 1, npoints(f)
269          WRITE( file_id_out_header, '(8F15.1)' )  points(1:3,n)
270       ENDDO
271!
272!--    Deallocate vertex data and close binary file.
273       DEALLOCATE( points )
274       
275       CLOSE ( file_id_in )
276    ENDDO
277!
278!-- Now, treat polygon data.
279    DO  f = 0, num_pe - 1
280!
281!--    Create filename of the treated binary file .   
282       CALL surface_output_create_file_string
283!
284!--    Open file with surface output for processor f.
285       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
286              TRIM( myid_char ), FORM = 'UNFORMATTED' )
287!
288!--    Move to last postion.
289!        CALL posix_lseek( file_id_in, filepos(f) )
290       CALL FSEEK( file_id_in, filepos(f), 0, ndum )
291!
292!--    Allocate array for polygon data
293       ALLOCATE( polygons(5,1:ns(f)) )
294!
295!--    Read polygon data and store them in a temporary file.
296       READ ( file_id_in )  polygons
297!
298!--    Obtain current file position after reading the local polygon data.
299!--    Will be used for next file opening.
300!        CALL posix_ftell( file_id_in, filepos(f) )
301       filepos(f) = ftell( file_id_in )
302!
303!--    Write further header information. Only one time required.
304       IF ( f == 0 )                                                           &
305          WRITE ( file_id_out_header, '(A,2I10)') "POLYGONS ",                 &
306                                                  ns_total, 5 * ns_total
307!
308!--    Write the polygons into the header file. Note, polygon data is of type
309!--    integer, as it just connects the point-indices which describe the given
310!--    surface element.
311       DO n = 1, ns(f)
312          WRITE ( file_id_out_header, '(5I18)' )  INT( polygons(1:5,n) )
313       ENDDO
314!
315!--    Deallocate array for polygon data and close the file.
316       DEALLOCATE( polygons )
317       
318       CLOSE ( file_id_in )
319       
320    ENDDO
321   
322    f = 0
323    CALL surface_output_create_file_string
324!
325!-- Write further header information. Only once required.
326    WRITE ( file_id_out_header, '(A,I10)') "CELL_DATA ", ns_total
327    WRITE ( file_id_out_header, '(A,I10)') "SCALARS cell_scalars float 1 "
328    WRITE ( file_id_out_header, '(A,I10)') "LOOKUP_TABLE default "
329!
330!-- Header creation has now been completed. File can be closed.   
331    CLOSE( file_id_out_header )
332!
333!-- Now, read all the actual data from the surface output. Please note, Paraview
334!-- VTK format expects 1 file per variable and time step.
335!-- The output file is only created once and includes the variable and the
336!-- simulated time.
337!-- In the binaray files, several variables and timesteps are stored, data for a
338!-- given variable, however, is distributed over all binary files. Hence, read
339!-- variable data for a given timestep from the binary file, write this data into
340!-- the target output file, remember file position in binary file and close it,
341!-- open nex binary file and read the variable, and so on, until all variables
342!-- for all timesteps are processed.
343    eof = .FALSE.
344    DO
345       DO  f = 0, num_pe - 1
346!
347!--       Clean up strings-
348          char_time = ''
349          variable_name = ''
350!
351!--       Create filename of the treated binary file.           
352          CALL surface_output_create_file_string
353!
354!--       Open binary file with surface output for processor f.
355          OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //         &
356                 TRIM( myid_char ), FORM = 'UNFORMATTED' )
357!
358!--       Move to last postion.
359!           CALL posix_lseek( file_id_in, filepos(f) )
360          CALL FSEEK( file_id_in, filepos(f), 0, ndum )
361!
362!--       Read string length and string indicating the output time step.
363          READ ( file_id_in ) length
364          READ ( file_id_in ) char_time(1:length)
365!
366!--       If string for the output time indicates the end-of-file, set the eof
367!--       flag and skip the read of the loop.
368          IF ( char_time(1:length) == 'END' )  THEN
369             eof(f) = .TRUE. 
370             CLOSE ( file_id_in )
371             CYCLE
372          ENDIF
373!
374!--       Read output time, and variable name.
375          READ ( file_id_in ) simulated_time
376          READ ( file_id_in ) length
377          READ ( file_id_in ) variable_name(1:length)
378!
379!--       For first loop index, open the target output file. First create the
380!--       filename string. Further, copy HEADER file with the given filename
381!--       string. The header information must be given in each VTK file!
382          IF ( f == 0 )  THEN
383             IF ( simulated_time < 10.0_wp )  THEN
384                WRITE( char_dum, '(I1)' )  INT( simulated_time ) 
385             ELSEIF ( simulated_time < 100.0_wp )  THEN
386                WRITE( char_dum, '(I2)' )  INT( simulated_time )
387             ELSEIF ( simulated_time < 1000.0_wp )  THEN
388                WRITE( char_dum, '(I3)' )  INT( simulated_time )
389             ELSEIF ( simulated_time < 10000.0_wp )  THEN
390                WRITE( char_dum, '(I4)' )  INT( simulated_time )   
391             ELSEIF ( simulated_time < 100000.0_wp )  THEN
392                WRITE( char_dum, '(I5)' )  INT( simulated_time )   
393             ELSEIF ( simulated_time < 1000000.0_wp )  THEN
394                WRITE( char_dum, '(I6)' )  INT( simulated_time )
395             ELSEIF ( simulated_time < 10000000.0_wp )  THEN
396                WRITE( char_dum, '(I7)' )  INT( simulated_time )   
397             ELSEIF ( simulated_time < 100000000.0_wp )  THEN
398                WRITE( char_dum, '(I8)' )  INT( simulated_time )   
399             ELSEIF ( simulated_time < 1000000000.0_wp )  THEN
400                WRITE( char_dum, '(I9)' )  INT( simulated_time )
401             ENDIF             
402!
403!--          Copy HEADER file
404             CALL system('cp HEADER ' // TRIM( path ) // TRIM( char_dum ) //   & 
405                                      's_' // TRIM( variable_name ) // '.vtk')
406!--          Open VTK file.
407             OPEN ( file_id_out, FILE = TRIM( path ) // TRIM( char_dum ) //    & 
408                    's_' // TRIM( variable_name ) // '.vtk', FORM='FORMATTED', &
409                    POSITION = 'APPEND' )           
410          ENDIF
411!
412!--       Allocate and read array for variable data. 
413          ALLOCATE( var(1:ns(f)) )
414       
415          READ( file_id_in ) var
416!
417!--       Write variable data into output VTK file.
418          DO n = 1, ns(f)
419             WRITE( file_id_out, * ) var(n)
420          ENDDO
421!
422!--       Remember file position in binary file and close it.
423!           CALL posix_ftell( file_id_in, filepos(f) )
424          filepos(f) = ftell( file_id_in )
425       
426          CLOSE ( file_id_in )
427!
428!--       Deallocate temporary array for variable data.
429          DEALLOCATE( var )
430       
431       ENDDO
432!
433!--    After data for a variable for one time step is read, close the output
434!--    VTK file and go to next variable or timestep.
435       CLOSE ( file_id_out )
436!
437!--    If all files reached the end-of-file, exit the loop.
438       IF ( ALL( eof ) )  EXIT
439       
440    ENDDO
441!
442!-- Finally, remove HEADER file
443    CALL system( 'rm HEADER' )
444   
445 CONTAINS
446 
447!------------------------------------------------------------------------------!
448! Description:
449! ------------
450!> This subroutine read the namelist file.
451!------------------------------------------------------------------------------!
452    SUBROUTINE surface_output_parin
453       
454       IMPLICIT NONE
455       
456       INTEGER(iwp) ::  file_id_parin = 90
457       
458       NAMELIST /surface_output/  convert_average_data, cycle_number, num_pe,  &
459                                  path, run
460
461!
462!--    Open namelist file.
463       OPEN( file_id_parin, FILE='surface_output_parin',                       &
464             STATUS='OLD', FORM='FORMATTED')
465!
466!--    Read namelist.
467       READ ( file_id_parin, surface_output )
468!
469!--    Close namelist file.
470       CLOSE( file_id_parin )
471       
472    END SUBROUTINE surface_output_parin
473     
474!------------------------------------------------------------------------------!
475! Description:
476! ------------
477!> This subroutine creates the filename string of the treated binary file.
478!------------------------------------------------------------------------------!
479    SUBROUTINE surface_output_create_file_string
480       
481       IMPLICIT NONE
482       
483       CHARACTER(LEN=3) ::  char_av = ''
484       CHARACTER(LEN=4) ::  char_cycle = ''
485       
486!
487!--    Create substring for the cycle number.
488       IF ( cycle_number /= 0 )  THEN
489          IF ( cycle_number < 10 )  THEN
490             WRITE( char_cycle, '(I1)')  cycle_number
491             char_cycle = '.00' // TRIM( char_cycle )
492          ELSEIF ( cycle_number < 100 )  THEN
493             WRITE( char_cycle, '(I2)')  cycle_number
494             char_cycle = '.0' // TRIM( char_cycle )
495          ELSEIF ( cycle_number < 1000 )  THEN
496             WRITE( char_cycle, '(I3)')  cycle_number
497             char_cycle = '.' // TRIM( char_cycle )
498          ENDIF
499       ENDIF
500!
501!--    Create substring for averaged data.
502       IF ( convert_average_data )  char_av = '_av'
503!
504!--    Create substring for the processor id and combine all substrings.
505       IF ( f < 10 )  THEN
506          WRITE( char_dum, '(I1)')  f
507          myid_char = TRIM( char_av ) // '_surf_00000' // TRIM( char_dum ) //  &
508                      TRIM( char_cycle ) // file_suffix
509       ELSEIF ( f < 100     )  THEN
510          WRITE( char_dum, '(I2)')  f
511          myid_char = TRIM( char_av ) // '_surf_0000'  // TRIM( char_dum ) //  &
512                      TRIM( char_cycle ) // file_suffix
513       ELSEIF ( f < 1000    )  THEN
514          WRITE( char_dum, '(I3)')  f
515          myid_char = TRIM( char_av ) // '_surf_000'   // TRIM( char_dum ) //  &
516                      TRIM( char_cycle ) // file_suffix
517       ELSEIF ( f < 10000   )  THEN
518          WRITE( char_dum, '(I4)')  f
519          myid_char = TRIM( char_av ) // '_surf_00'    // TRIM( char_dum ) //  &
520                      TRIM( char_cycle ) // file_suffix
521       ELSEIF ( f < 100000  )  THEN
522          WRITE( char_dum, '(I5)')  f
523          myid_char = TRIM( char_av ) // '_surf_0'     // TRIM( char_dum ) //  &
524                      TRIM( char_cycle ) // file_suffix
525       ELSEIF ( f < 1000000 )  THEN
526          WRITE( char_dum, '(I6)')  f
527          myid_char = TRIM( char_av ) // '_surf_'      // TRIM( char_dum ) //  &
528                      TRIM( char_cycle ) // file_suffix
529       ENDIF
530       
531    END SUBROUTINE surface_output_create_file_string
532   
533 END PROGRAM surface_output_to_vtk
534
535
536
Note: See TracBrowser for help on using the repository browser.