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

Last change on this file since 3698 was 3523, checked in by suehring, 6 years ago

Implement interface for posix conform C-systemcalls in order to replace non-standard FORTRAN functions ftell and fseek

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