Home

Context Navigation

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

Last change on this file since 3602 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 sward $
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |