Home

Context Navigation

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

Last change on this file since 3775 was 3755, checked in by suehring, 6 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 gronemeier $
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |