PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
output_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations concerning giving output, on the screen as well as in output
3!! files.
4!------------------------------------------------------------------------------!
6#include <PB3D_macros.h>
8 use messages
11 use files_utilities, only: nextunit
12 use grid_vars, only: grid_type
13
14
15 implicit none
16 private
18
19 ! global variables
20 integer :: temp_id(2) = 1 !< \private will be appended to temporary data and script files
21 character(len=9) :: line_clrs(20) = [&
22 &'"#7297E6"','"#67EB84"','"#F97A6D"','"#F9C96D"','"#1D4599"',&
23 &'"#11AD34"','"#E62B17"','"#E69F17"','"#2F3F60"','"#2F6C3D"',&
24 &'"#8F463F"','"#8F743F"','"#031A49"','"#025214"','"#6D0D03"',&
25 &'"#6D4903"','"#A9BDE6"','"#A6EBB5"','"#F9B7B0"','"#F9E0B0"'] !< \private line colors, from <https://github.com/Gnuplotting/gnuplot-palettes>
26 character(len=max_str_ln) :: line_style = 'lt 1 lw 1 pt 7 ps 0.5;' !< \private line style without little space in line connecting points
27#if ldebug
28 character(len=0) :: err_output_str = '' !< \private string with error output
29#else
30 character(len=14) :: err_output_str = ' 2> /dev/null' !< \private string with error output (/dev/null)
31#endif
32
33 ! interfaces
34
35 !> \public Print 2-D output on a file.
36 !!
37 !! The variables \c var_name and \c file_name hold the name of the plot and
38 !! of the file in which the plot data is to be saved, respectively. \c y is
39 !! the an array containing the function which is stored and \c x is an
40 !! optional vector with the x-values. The logical \c draw can optionally
41 !! disable calling the external drawing procedure for output on screen
42 !! [default], without modifying the plot file.
43 !!
44 !! The first index of \c y (and \c x) contains the points of a current plot.
45 !!
46 !! The second index indicates various plots (one or more)
47 interface print_ex_2d
48 !> \public
49 module procedure print_ex_2d_ind
50 !> \public
51 module procedure print_ex_2d_arr
52 end interface
53
54 !> \public Print 3-D output on a file.
55 !!
56 !! The variables \c var_name and \c file_name hold the name of the plot and
57 !! of the file in which the plot data is to be saved, respectively. \c z is
58 !! the an array containing the function which is stored and \c x and\c y are
59 !! optional vectors with the x and y-values. The logical \c draw can
60 !! optionally disable calling the external drawing procedure for output on
61 !! screen [default], without modifying the plot file.
62 !!
63 !! The first index of \c z (and \c x, \c y) contains the points of a current
64 !! The plot second index indicates various plots (one or more)
65 interface print_ex_3d
66 !> \public
67 module procedure print_ex_3d_ind
68 !> \public
69 module procedure print_ex_3d_arr
70 end interface
71
72
73 !> \public Prints variables \c vars with names \c var_names in an HDF5 file
74 !! with name c file_name and accompanying XDMF file.
75 !!
76 !! For XDMF collections (\cite xdmf), only the first value for \c var_names
77 !! is used, so it should have a size of one.
78 !!
79 !! The plot is generally 3-D, but if one of the dimensions provided is equal
80 !! to 1, it is checked whether there is poloidal or toroidal axisymmetry and
81 !! if so, the plot becomes 2-D. This can be forced using the optional input
82 !! argument \c sym_type.
83 !!
84 !! Optionally, the (curvilinear) grid can be provided through \c X, \c Y and
85 !! \c Z. If not, the grid is assumed to be Cartesian with discrete indices
86 !! where \c X corresponds to the first dimensions, \c Y to the second and \c
87 !! Z to the third.
88 !!
89 !! Additionally, the total grid size and local offset can be provided in \c
90 !! tot_dim and \c loc_offset to run this routine in parallel automatically.
91 !! Optionally, one of the dimensions (\c col_id, default 4) can be
92 !! associated to a collection dimension using \c col different from 1:
93 !! - \c col = 1: no collection, just plots of different variables
94 !! - \c col = 2: time collection
95 !! - \c col = 3: spatial collection
96 !! - \c col = 4: vector field
97 !!
98 !! Furthermore, using the variable \c cont_plot, a plot can be
99 !! (over-)written in multiple writes. By this is meant that there should be
100 !! an initial plot, with collection type 1, 2, 3 or 4, which can then be
101 !! followed by an arbitrary number of additional writes. As these additional
102 !! writes currently cannot modify the plot structure, nor the XDMF
103 !! variables, their collection dimension should be complete from the start.
104 !!
105 !! This has no implications for single plots but means that for collection
106 !! types 2 to 4 all the elements in the collection have to be present,
107 !! though they do not necessary need to have been completely written in the
108 !! other dimensions.
109 !!
110 !! Subsequent writes with \c cont_plot can then, for instance, write parts
111 !! of the data that had not yet been written, or overwrite ones that had.
112 !! This can be useful for post-processing where the memory requirements are
113 !! large so that the work has to be split.
114 !!
115 !! \note
116 !! -# For a vector field, the number of variables has to be 2 for 2-D plots
117 !! or 3 for 3-D plots. This should be rewritten in the future, so that
118 !! collections can be used for vectors as well.
119 !! -# In order to merge collections in their collection dimension, the XDMF
120 !! files can always easily be joined.
121 !! -# If necessary, a lock system should be used when multiple processes
122 !! are writing the same file, including continued writes.
123 !! -# To plot this with VisIt, use:
124 !! - for temporal collections: pseudocolor using the variable name
125 !! (other names are ignored).
126 !! - for spatial collections: subset of blocks or pseudocolor using the
127 !! variable name (other names are ignored).
128 !! - for vector plot: Vector plot.
129 !! - for without collections: pseudocolor using the variable names.
130 !! -# Currently all of possibly multiple processes that write data have to
131 !! cover the entire range of the variables in the dimension indicated by \c
132 !! col_id. (This could be implemented by changing how n_plot is defined and
133 !! selectively letting each processer write in the main loop at its
134 !! corresponding indices.)
135 !! -# To project the data to 2-D in VisIt, use the projection tool under
136 !! Operators > Transform
137 interface plot_hdf5
138 !> \public
139 module procedure plot_hdf5_ind
140 !> \public
141 module procedure plot_hdf5_arr
142 end interface
143
144contains
145 !> \private individual version
146 subroutine print_ex_2d_ind(var_name,file_name_i,y,x,draw,persistent)
147 ! input / output
148 character(len=*), intent(in) :: var_name !< name of variable in legend
149 character(len=*), intent(in) :: file_name_i !< name of input file
150 real(dp), intent(in) :: y(1:) !< ordinate
151 real(dp), intent(in), optional :: x(1:) !< absicca
152 logical, intent(in), optional :: draw !< whether to draw the plot as well
153 logical, intent(in), optional :: persistent !< keep on-screen plot open
154
155 ! local variables
156 integer :: npoints
157
158 ! set npoints, ny
159 npoints = size(y)
160
161 ! call multiplot version
162 if (present(x)) then
163 call print_ex_2d_arr([var_name],file_name_i,reshape(y,[npoints,1]),&
164 &x=reshape(x,[npoints,1]),draw=draw,persistent=persistent)
165 else
166 call print_ex_2d_arr([var_name],file_name_i,reshape(y,[npoints,1]),&
167 &draw=draw,persistent=persistent)
168 end if
169 end subroutine print_ex_2d_ind
170 !> \private array version
171 subroutine print_ex_2d_arr(var_names,file_name_i,y,x,draw,persistent)
172 use num_vars, only: rank
173
174 ! input / output
175 character(len=*), intent(in) :: var_names(:) !< names of variables in legend
176 character(len=*), intent(in) :: file_name_i !< name of input file
177 real(dp), intent(in) :: y(1:,1:) !< ordinate
178 real(dp), intent(in), optional :: x(1:,1:) !< absicca
179 logical, intent(in), optional :: draw !< whether to draw the plot as well
180 logical, intent(in), optional :: persistent !< keep on-screen plot open
181
182 ! local variables
183 integer :: file_i
184 integer :: iplt, ipnt, nplt, npnt
185 real(dp), allocatable :: x_fin(:,:)
186 integer :: istat
187 character(len=max_str_ln) :: file_name
188 character(len=max_str_ln) :: write_format
189 logical :: persistent_loc ! local persistent
190
191 ! set nplt, npnt
192 npnt = size(y,1)
193 nplt = size(y,2)
194 if (present(x)) then
195 if (size(x,1).ne.size(y,1) .or. &
196 &(size(x,2).ne.1 .and. size(x,2).ne.size(y,2))) then
197 write(*,*) 'shape(x)', shape(x)
198 write(*,*) 'shape(y)', shape(y)
199 call writo('In print_ex_2D, the size of x and y has to be the &
200 &same... Skipping plot',persistent=.true.,warning=.true.)
201 end if
202 end if
203
204 ! set other variables
205 persistent_loc = .false.
206 if (present(persistent)) persistent_loc = persistent
207
208 ! set x
209 allocate(x_fin(npnt,nplt))
210 if (present(x)) then
211 if (size(x,2).eq.1) then ! assume copy
212 do iplt = 1,nplt
213 x_fin(:,iplt) = x(:,1)
214 end do
215 else
216 x_fin = x
217 end if
218 else
219 do ipnt = 1,npnt
220 x_fin(ipnt,:) = ipnt
221 end do
222 end if
223
224 ! set default file name if empty
225 if (trim(file_name_i).eq.'') then
226 file_name = 'temp_data_print_ex_2D_'//trim(i2str(rank))//'_'//&
227 trim(i2str(temp_id(1)))
228 temp_id(1) = temp_id(1)+1
229 else
230 file_name = trim(file_name_i)
231 end if
232
233 ! open output file
234 open(unit=nextunit(file_i),file=data_dir//'/'//trim(file_name)//'.dat',&
235 &iostat=istat)
236 chckstt
237
238 ! write to output file
239 write(file_i,'(1X,A)',iostat=istat) &
240 &'# '//trim(merge_strings(var_names))//':'
241 write_format="("//trim(i2str(nplt*2))//"(ES23.16,' '))"
242 do ipnt = 1,npnt
243 write(file_i,fmt=trim(write_format),iostat=istat) &
244 &(x_fin(ipnt,iplt), iplt = 1,nplt), &
245 &(y(ipnt,iplt), iplt = 1,nplt)
246 enddo
247 write(file_i,'(1X,A)',iostat=istat) ''
248
249 ! close output file
250 close(file_i)
251
252 ! bypass plots if no_plots
253 if (no_plots) return
254
255 ! if draw is present and equal to .false., cancel calling draw_ex
256 if (present(draw)) then
257 if (.not.draw) return
258 end if
259 call draw_ex(var_names,file_name,nplt,1,.true.,persistent=persistent)
260
261 if (trim(file_name_i).eq.'' .and. .not.persistent_loc) then
262 call use_execute_command_line('rm '//data_dir//'/'//&
263 &trim(file_name)//'.dat')
264 end if
265 end subroutine print_ex_2d_arr
266
267 !> \private individual version
268 subroutine print_ex_3d_ind(var_name,file_name_i,z,y,x,draw)
269 ! input / output
270 character(len=*), intent(in) :: var_name !< name of variable in legend
271 character(len=*), intent(in) :: file_name_i !< name of input file
272 real(dp), intent(in) :: z(1:,1:) !< ordinate
273 real(dp), intent(in), optional :: x(1:,1:) !< absicca
274 real(dp), intent(in), optional :: y(1:,1:) !< absicca
275 logical, intent(in), optional :: draw !< whether to draw the plot as well
276
277 ! local variables
278 integer :: npoints(2)
279
280 ! set npoints, ny
281 npoints = [size(z,1),size(z,2)]
282
283 ! call multiplot version
284 if (present(y)) then
285 if (present(x)) then
286 call print_ex_3d_arr([var_name],file_name_i,&
287 &reshape(z,[size(z,1),size(z,2),1]),&
288 &y=reshape(y,[size(z,1),size(z,2),1]),&
289 &x=reshape(x,[size(z,1),size(z,2),1]),draw=draw)
290 else
291 call print_ex_3d_arr([var_name],file_name_i,&
292 &reshape(z,[size(z,1),size(z,2),1]),&
293 &y=reshape(y,[size(z,1),size(z,2),1]),&
294 &draw=draw)
295 end if
296 else
297 if (present(x)) then
298 call print_ex_3d_arr([var_name],file_name_i,&
299 &reshape(z,[size(z,1),size(z,2),1]),&
300 &x=reshape(x,[size(z,1),size(z,2),1]),&
301 &draw=draw)
302 else
303 call print_ex_3d_arr([var_name],file_name_i,&
304 &reshape(z,[size(z,1),size(z,2),1]),draw=draw)
305 end if
306 end if
307 end subroutine print_ex_3d_ind
308 !> \private array version
309 subroutine print_ex_3d_arr(var_names,file_name_i,z,x,y,draw)
310 use num_vars, only: rank
311
312 ! input / output
313 character(len=*), intent(in) :: var_names(:) !< names of variables in legend
314 character(len=*), intent(in) :: file_name_i !< name of input file
315 real(dp), intent(in) :: z(1:,1:,1:) !< ordinate
316 real(dp), intent(in), optional :: x(1:,1:,1:) !< absicca
317 real(dp), intent(in), optional :: y(1:,1:,1:) !< absicca
318 logical, intent(in), optional :: draw !< whether to draw the plot as well
319
320 ! local variables
321 integer :: file_i
322 integer :: iplt, ipntx, ipnty, nplt, npntx, npnty
323 real(dp), allocatable :: x_fin(:,:,:)
324 real(dp), allocatable :: y_fin(:,:,:)
325 integer :: istat
326 character(len=max_str_ln) :: file_name
327 character(len=max_str_ln) :: write_format
328
329 ! set nplt, npnt
330 npntx = size(z,1)
331 npnty = size(z,2)
332 nplt = size(z,3)
333 if (present(x)) then
334 if (size(x,1).ne.size(z,1) .or. size(x,2).ne.size(z,2) .or. &
335 &(size(x,3).ne.1 .and. size(x,3).ne.size(z,3))) then
336 write(*,*) 'shape(x)', shape(x)
337 write(*,*) 'shape(z)', shape(z)
338 call writo('In print_ex, the size of x and z has to be the &
339 &same... Skipping plot',persistent=.true.,warning=.true.)
340 return
341 end if
342 end if
343 if (present(y)) then
344 if (size(y,1).ne.size(z,1) .or. size(y,2).ne.size(z,2) .or. &
345 &(size(y,3).ne.1 .and. size(y,3).ne.size(z,3))) then
346 write(*,*) 'shape(y)', shape(y)
347 write(*,*) 'shape(z)', shape(z)
348 call writo('In print_ex, the size of y and z has to be the &
349 &same... Skipping plot',persistent=.true.,warning=.true.)
350 return
351 end if
352 end if
353
354 ! set x
355 allocate(x_fin(npntx,npnty,nplt))
356 if (present (x)) then
357 if (size(x,3).eq.1) then ! assume copy
358 do iplt = 1,nplt
359 x_fin(:,:,iplt) = x(:,:,1)
360 end do
361 else
362 x_fin = x
363 end if
364 else
365 do ipntx = 1,npntx
366 x_fin(ipntx,:,:) = ipntx
367 end do
368 end if
369 ! set y
370 allocate(y_fin(npntx,npnty,nplt))
371 if (present (y)) then
372 if (size(y,3).eq.1) then ! assume copy
373 do iplt = 1,nplt
374 y_fin(:,:,iplt) = y(:,:,1)
375 end do
376 else
377 y_fin = y
378 end if
379 else
380 do ipnty = 1,npnty
381 y_fin(:,ipnty,:) = ipnty
382 end do
383 end if
384
385 ! set default file name if empty
386 if (trim(file_name_i).eq.'') then
387 file_name = 'temp_data_print_ex_3D_'//trim(i2str(rank))//'_'//&
388 trim(i2str(temp_id(1)))
389 temp_id(1) = temp_id(1)+1
390 else
391 file_name = trim(file_name_i)
392 end if
393
394 ! open output file
395 open(nextunit(file_i),file=data_dir//'/'//trim(file_name)//'.dat',&
396 &iostat=istat)
397 chckstt
398
399 ! write to output file
400 write(file_i,'(A)',iostat=istat) &
401 &'# '//trim(merge_strings(var_names))//':'
402 write_format="("//trim(i2str(nplt*3))//"(ES23.16,' '))"
403 do ipntx = 1,npntx
404 do ipnty = 1,npnty
405 write(file_i,fmt=trim(write_format),iostat=istat) &
406 &(x_fin(ipntx,ipnty,iplt), iplt = 1,nplt), &
407 &(y_fin(ipntx,ipnty,iplt), iplt = 1,nplt), &
408 &(z(ipntx,ipnty,iplt), iplt = 1,nplt)
409 end do
410 write(file_i,'(A)',iostat=istat) ''
411 enddo
412 write(file_i,'(A)',iostat=istat) ''
413
414 ! close output file
415 close(file_i)
416
417 ! bypass plots if no_plots
418 if (no_plots) return
419
420 ! if draw is present and equal to .false., cancel calling draw_ex
421 if (present(draw)) then
422 if (.not.draw) return
423 end if
424 call draw_ex(var_names,file_name,nplt,2,.true.)
425
426 if (trim(file_name_i).eq.'') then
427 call use_execute_command_line('rm '//data_dir//'/'//&
428 &trim(file_name)//'.dat')
429 end if
430 end subroutine print_ex_3d_arr
431
432 !> \private array version
433 subroutine plot_hdf5_arr(var_names,file_name,vars,tot_dim,loc_offset,&
434 &X,Y,Z,col_id,col,sym_type,cont_plot,descr)
438 use hdf5_vars, only: dealloc_xml_str, &
440 use num_vars, only: n_procs
441 use mpi_utilities, only: get_ser_var
442
443 ! input / output
444 character(len=*), intent(in) :: var_names(:) !< names of variable to be plot
445 character(len=*), intent(in) :: file_name !< file name
446 real(dp), intent(in), target :: vars(:,:,:,:) !< variables to plot
447 integer, intent(in), optional :: tot_dim(4) !< total dimensions of the arrays
448 integer, intent(in), optional :: loc_offset(4) !< offset of local dimensions
449 real(dp), intent(in), target, optional :: X(:,:,:,:) !< curvilinear grid X points
450 real(dp), intent(in), target, optional :: Y(:,:,:,:) !< curvilinear grid Y points
451 real(dp), intent(in), target, optional :: Z(:,:,:,:) !< curvilinear grid Z points
452 integer, intent(in), optional :: col_id !< index of time dimension
453 integer, intent(in), optional :: col !< whether a collection is made
454 integer, intent(in), optional :: sym_type !< type of symmetry (1: no symmetry, 2: toroidal, 3: poloidal)
455 logical, intent(in), optional :: cont_plot !< continued plot
456 character(len=*), intent(in), optional :: descr !< description
457
458 ! local variables
459 type(hdf5_file_type) :: file_info ! file info
460 integer :: istat ! status
461 integer :: col_id_loc ! local copy of col_id
462 integer :: col_loc ! local copy of col
463 integer :: n_plot ! nr. of plots
464 integer :: id, jd ! counter
465 integer :: sym_type_loc ! local sym_type
466 integer :: sym_pol, sym_tor ! used to determine symmetry
467 integer :: att_id ! index in att
468 integer, allocatable :: tot_sym_pol(:), tot_sym_tor(:) ! sym_pol and sym_tor for all processes
469 integer :: tot_dim_loc(4) ! local copy of tot_dim
470 integer :: loc_offset_loc(4) ! local copy of loc_offset
471 integer :: tot_dim_3D(3) ! tot_dim except collection
472 integer :: loc_dim_3D(3) ! loc_dim except collection
473 integer :: loc_offset_3D(3) ! loc_offset except collection
474 type(xml_str_type) :: col_grid ! grid with collection
475 type(xml_str_type), allocatable :: grids(:) ! the grids in the time collection
476 type(xml_str_type) :: top ! topology
477 type(xml_str_type), allocatable :: dat(:) ! data items for geometry and attribute
478 type(xml_str_type) :: dat_vec ! data items for merged attribute
479 type(xml_str_type) :: geom ! geometry
480 type(xml_str_type) :: att(1) ! attribute
481 logical :: col_mask(4) ! to select out the collection dimension
482 logical :: ind_plot ! individual plot
483 logical :: cont_plot_loc ! local cont_plot
484 logical :: first_vc, last_vc ! first or last vector components (for the others, some XDMF operations can be skipped)
485 real(dp), allocatable :: sym_ang(:,:,:) ! angle to be checked for symmetry
486 real(dp) :: tol_sym = 1.e-8_dp ! tolerance for symmetry determination
487 real(dp), pointer :: var_3D(:,:,:) => null() ! pointer to vars
488 real(dp), pointer :: X_3D(:,:,:) => null() ! pointer to X
489 real(dp), pointer :: Y_3D(:,:,:) => null() ! pointer to Y
490 real(dp), pointer :: Z_3D(:,:,:) => null() ! pointer to Z
491 character(len=max_str_ln), allocatable :: grd_names(:) ! grid names
492 character(len=max_str_ln), allocatable :: att_names(:) ! attribute names
493
494 ! set up local col_id and col
495 col_id_loc = 4 ! default collection dimension: last index
496 if (present(col_id)) col_id_loc = col_id
497 col_loc = 1 ! default no spatial collection
498 if (present(col)) col_loc = col
499
500 ! set up nr. of plots
501 n_plot = size(vars,col_id_loc)
502
503 ! set up local tot_dim and loc_offset
504 tot_dim_loc = shape(vars)
505 if (present(tot_dim)) tot_dim_loc = tot_dim
506 loc_offset_loc = 0
507 if (present(loc_offset)) loc_offset_loc = loc_offset
508
509 ! tests
510 if (tot_dim_loc(col_id_loc).ne.n_plot) then
511 istat = 1
512 call writo('In plot_HDF5, all the processes need to have the full &
513 &range in the dimension given by col_id',persistent=.true.,&
514 &warning=.true.)
515 chckstt
516 end if
517 if (n_plot.eq.1 .and. col_loc.ne.1) then
518 istat = 1
519 call writo('In plot_HDF5, if single plot, the collection type &
520 &needs to be one',persistent=.true.,warning=.true.)
521 chckstt
522 end if
523
524 ! set 3D dimensions
525 col_mask = .false.
526 col_mask(col_id_loc) = .true.
527 tot_dim_3d = pack(tot_dim_loc,.not.col_mask)
528 loc_dim_3d = pack(shape(vars),.not.col_mask)
529 loc_offset_3d = pack(loc_offset_loc,.not.col_mask)
530
531 ! set up individual plot
532 if (tot_dim_3d(1).eq.loc_dim_3d(1) .and. &
533 &tot_dim_3d(2).eq.loc_dim_3d(2) .and. &
534 &tot_dim_3d(3).eq.loc_dim_3d(3)) then
535 ind_plot = .true.
536 else
537 ind_plot = .false.
538 end if
539
540 ! set up continued plot
541 cont_plot_loc = .false.
542 if (present(cont_plot)) cont_plot_loc = cont_plot
543
544 ! default symmetry type
545 sym_type_loc = 1
546
547 ! Find symmetry type by checking whether Y/X is constant (toroidal
548 ! symmetry) or Z^2/(X^2+Y^2) is constant (poloidal symmetry), for all
549 ! plots.
550 if (present(sym_type)) then
551 sym_type_loc = sym_type
552 else if (minval(tot_dim_3d).eq.1) then ! possibly symmetry
553 ! allocate helper variable
554 allocate(sym_ang(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
555 ! initialize sym_pol and sym_tor
556 sym_pol = 0
557 sym_tor = 0
558 ! loop over all plots
559 do id = 1,n_plot
560 ! assign pointers
561 call assign_pointers(id)
562 ! check poloidal angle
563 sym_ang = atan2(y_3d,x_3d)
564 if (maxval(sym_ang)-minval(sym_ang).lt.tol_sym .and. &
565 &(maxval(x_3d).ge.0._dp .neqv. &
566 &minval(x_3d).lt.0._dp).and.& ! X has to be either positive or negative
567 &(maxval(y_3d).ge.0._dp .neqv. &
568 &minval(y_3d).lt.0._dp)) & ! Y has to be either positive or negative
569 &sym_pol = sym_pol+1 ! poloidal symmetry for this plot
570 ! check toroidal angle
571 sym_ang = atan2(sqrt(z_3d**2),sqrt(x_3d**2+y_3d**2))
572 if (maxval(sym_ang)-minval(sym_ang).lt.tol_sym) &
573 &sym_tor = sym_tor+1 ! toroidal symmetry for this plot
574 end do
575 ! get total sym_pol and sym_tor
576 if (ind_plot) then ! so that below test succeeds
577 allocate(tot_sym_pol(n_procs),tot_sym_tor(n_procs))
578 tot_sym_pol = sym_pol
579 tot_sym_tor = sym_tor
580 else ! get from all the processes
581 istat = get_ser_var([sym_pol],tot_sym_pol,scatter=.true.)
582 chckstt
583 istat = get_ser_var([sym_tor],tot_sym_tor,scatter=.true.)
584 chckstt
585 end if
586
587 ! check total results
588 if (sum(tot_sym_pol).eq.n_procs*n_plot) then ! poloidal symmetry for all plots
589 sym_type_loc = 2
590 else if (sum(tot_sym_tor).eq.n_procs*n_plot) then ! toroidal symmetry for all plots
591 sym_type_loc = 3
592 end if
593 ! deallocate helper variables
594 deallocate(sym_ang)
595 end if
596
597 ! tests
598 if (col_loc.eq.4) then
599 if (sym_type_loc.eq.1) then
600 if (n_plot.ne.3) then
601 istat = 1
602 call writo('For vector field plots, need 3 dimensions',&
603 &warning=.true.)
604 chckstt
605 end if
606 else
607 if (n_plot.ne.2) then
608 istat = 1
609 call writo('For symmetric vector field plots, need 2 &
610 &dimensions',warning=.true.)
611 chckstt
612 end if
613 end if
614 end if
615
616 ! set up and local var names if not a continued plot
617 if (.not.cont_plot_loc) then
618 ! set up local var_names
619 allocate(grd_names(n_plot))
620 allocate(att_names(n_plot))
621 if (col_loc.eq.1) then ! without collection
622 if (n_plot.eq.1) then ! just one plot: attribute name is important
623 if (size(var_names).eq.n_plot) then ! the right number of variable names provided
624 att_names = var_names
625 else if (size(var_names).gt.n_plot) then ! too many variable names provided
626 att_names = var_names(1:n_plot)
627 call writo('Too many variable names provided',&
628 &persistent=.true.,warning=.true.)
629 else ! not enough variable names provided
630 att_names(1:size(var_names)) = var_names
631 do id = size(var_names)+1,n_plot
632 att_names(id) = 'unnamed variable '//trim(i2str(id))
633 end do
634 call writo('Not enough variable names provided',&
635 &persistent=.true.,warning=.true.)
636 end if
637 grd_names = 'default_grid_name'
638 else ! multiple plots: grid name is important
639 if (size(var_names).eq.n_plot) then ! the right number of variable names provided
640 grd_names = var_names
641 else if (size(var_names).gt.n_plot) then ! too many variable names provided
642 grd_names = var_names(1:n_plot)
643 call writo('Too many variable names provided',&
644 &persistent=.true.,warning=.true.)
645 else ! not enough variable names provided
646 grd_names(1:size(var_names)) = var_names
647 do id = size(var_names)+1,n_plot
648 grd_names(id) = 'unnamed variable '//trim(i2str(id))
649 end do
650 call writo('Not enough variable names provided',&
651 &persistent=.true.,warning=.true.)
652 end if
653 att_names = 'default_att_name'
654 end if
655 else ! collections: attribute name is important
656 att_names = var_names(1)
657 if (size(var_names).gt.1) call writo('For collections, only &
658 &the first variable name is used',persistent=.true.,&
659 &warning=.true.)
660 grd_names = 'default_grid_name'
661 end if
662 end if
663
664 ! open HDF5 file
665 istat = open_hdf5_file(file_info,file_name,sym_type=sym_type_loc,&
666 &descr=descr,ind_plot=ind_plot,cont_plot=cont_plot_loc)
667 chckstt
668
669 ! create grid for collection if not continued plot
670 if (.not.cont_plot_loc) allocate(grids(n_plot))
671
672 ! allocate data item arrays
673 if (sym_type_loc.eq.1) then ! 3D grid
674 allocate(dat(3))
675 else ! 2D grid
676 allocate(dat(2))
677 end if
678
679 ! loop over all plots
680 do id = 1,n_plot
681 ! set up whether last vector component
682 ! (for scalar plots, this is always true)
683 first_vc = (id.eq.1 .or. col_loc.ne.4)
684 last_vc = (id.eq.n_plot .or. col_loc.ne.4)
685
686 ! print topology if not continued plot and last vector components
687 if (.not.cont_plot_loc .and. last_vc) then
688 if (sym_type_loc.eq.1) then ! 3D grid
689 call print_hdf5_top(top,2,tot_dim_3d,ind_plot=ind_plot)
690 else ! 2D grid
691 call print_hdf5_top(top,1,tot_dim_3d,ind_plot=ind_plot)
692 end if
693 end if
694
695 ! assign pointers
696 call assign_pointers(id)
697
698 ! print data item for X, Y and Z (no symmetry), R = sqrt(X^2+Y^2)
699 ! and Z (poloidal symmetry) or X and Y (toroidal symmetry) if first
700 ! vector component (afterwards, dat is overwritten)
701 if (first_vc) then
702 select case (sym_type_loc)
703 case (1) ! no symmetry
704 istat = print_hdf5_3d_data_item(dat(1),file_info,&
705 &'X_'//trim(i2str(id)),x_3d,tot_dim_3d,loc_dim_3d,&
706 &loc_offset_3d,ind_plot=ind_plot,&
707 &cont_plot=cont_plot_loc)
708 chckstt
709 istat = print_hdf5_3d_data_item(dat(2),file_info,&
710 &'Y_'//trim(i2str(id)),y_3d,tot_dim_3d,loc_dim_3d,&
711 &loc_offset_3d,ind_plot=ind_plot,&
712 &cont_plot=cont_plot_loc)
713 chckstt
714 istat = print_hdf5_3d_data_item(dat(3),file_info,&
715 &'Z_'//trim(i2str(id)),z_3d,tot_dim_3d,loc_dim_3d,&
716 &loc_offset_3d,ind_plot=ind_plot,&
717 &cont_plot=cont_plot_loc)
718 chckstt
719 case (2) ! poloidal symmetry
720 istat = print_hdf5_3d_data_item(dat(1),file_info,&
721 &'R_'//trim(i2str(id)),sqrt(x_3d**2+y_3d**2),&
722 &tot_dim_3d,loc_dim_3d,loc_offset_3d,&
723 &ind_plot=ind_plot,cont_plot=cont_plot_loc)
724 chckstt
725 istat = print_hdf5_3d_data_item(dat(2),file_info,&
726 &'Z_'//trim(i2str(id)),z_3d,tot_dim_3d,loc_dim_3d,&
727 &loc_offset_3d,ind_plot=ind_plot,&
728 &cont_plot=cont_plot_loc)
729 chckstt
730 case (3) ! toroidal symmetry
731 istat = print_hdf5_3d_data_item(dat(1),file_info,&
732 &'X_'//trim(i2str(id)),x_3d,tot_dim_3d,loc_dim_3d,&
733 &loc_offset_3d,ind_plot=ind_plot,&
734 &cont_plot=cont_plot_loc)
735 chckstt
736 istat = print_hdf5_3d_data_item(dat(2),file_info,&
737 &'Y_'//trim(i2str(id)),y_3d,tot_dim_3d,loc_dim_3d,&
738 &loc_offset_3d,ind_plot=ind_plot,&
739 &cont_plot=cont_plot_loc)
740 chckstt
741 case default ! no symmetry
742 istat = 1
743 call writo('symmetry type '//&
744 &trim(i2str(sym_type_loc))//' not recognized',&
745 &persistent=.true.,warning=.true.)
746 chckstt
747 end select
748 end if
749
750 ! print geometry with X, Y and Z data item if not continued plot and
751 ! first vector component
752 if (.not.cont_plot_loc .and. first_vc) then
753 if (sym_type_loc.eq.1) then ! no symmetry so 3D geometry
754 call print_hdf5_geom(geom,2,dat,reset=.true.,&
755 &ind_plot=ind_plot)
756 else ! symmetry so 2D geometry
757 call print_hdf5_geom(geom,1,dat,reset=.true.,&
758 &ind_plot=ind_plot)
759 end if
760 end if
761
762 ! print data item for plot variable
763 if (col_loc.eq.4) then
764 att_id = id
765 else
766 att_id = 1
767 end if
768 istat = print_hdf5_3d_data_item(dat(att_id),file_info,'var_'//&
769 &trim(i2str(id)),var_3d,tot_dim_3d,loc_dim_3d,&
770 &loc_offset_3d,ind_plot=ind_plot,cont_plot=cont_plot_loc)
771 chckstt
772
773 ! print attribute with this data item if not continued plot
774 if (.not.cont_plot_loc) then
775 if (col_loc.eq.4) then
776 if (last_vc) then
777 call merge_hdf5_3d_data_items(dat_vec,dat,'var_vec',&
778 &tot_dim_3d,reset=.true.,ind_plot=ind_plot,&
779 &cont_plot=cont_plot_loc)
780 call print_hdf5_att(att(1),dat_vec,&
781 &att_names(1),1,2,reset=.true.,ind_plot=ind_plot)
782 end if
783 else
784 call print_hdf5_att(att(1),dat(1),&
785 &att_names(id),1,1,reset=.true.,ind_plot=ind_plot)
786 end if
787 end if
788
789 ! create a grid with the topology, the geometry, the attribute and
790 ! time if time collection, if not continued plot and last vector
791 ! component
792 if (.not.cont_plot_loc .and. last_vc) then
793 if (col_loc.eq.2) then ! time collection
794 istat = print_hdf5_grid(grids(id),grd_names(id),1,&
795 &grid_time=id*1._dp,grid_top=top,grid_geom=geom,&
796 &grid_atts=att,reset=.true.,ind_plot=ind_plot)
797 chckstt
798 else ! no time collection
799 istat = print_hdf5_grid(grids(id),grd_names(id),1,&
800 &grid_top=top,grid_geom=geom,grid_atts=att,&
801 &reset=.true.,ind_plot=ind_plot)
802 chckstt
803 end if
804 end if
805 end do
806
807 ! either create collection or just use individual grids if not continued
808 ! plot
809 if (.not.cont_plot_loc) then
810 if (col_loc.eq.1 .or. col_loc.eq.4) then
811 ! add individual grids to HDF5 file and reset them
812 do id = 1,n_plot
813 if (col_loc.eq.4 .and. id.lt.n_plot) cycle ! only plot last for vector
814 istat = add_hdf5_item(file_info,grids(id),reset=.true.,&
815 &ind_plot=ind_plot)
816 chckstt
817 end do
818 else
819 ! create grid collection from individual grids and reset them
820 istat = print_hdf5_grid(col_grid,'domain of mesh',col_loc,&
821 &grid_grids=grids,reset=.true.,ind_plot=ind_plot)
822 chckstt
823
824 ! add collection grid to HDF5 file and reset it
825 istat = add_hdf5_item(file_info,col_grid,reset=.true.,&
826 &ind_plot=ind_plot)
827 chckstt
828 end if
829 end if
830
831 ! close HDF5 file
832 istat = close_hdf5_file(file_info,ind_plot=ind_plot,&
833 &cont_plot=cont_plot_loc)
834 chckstt
835
836 ! clean up
837 nullify(var_3d)
838 nullify(x_3d,y_3d,z_3d)
839 call dealloc_xml_str(dat)
840 if (.not.cont_plot_loc) then
841 call dealloc_xml_str(grids)
842 call dealloc_xml_str(top)
843 call dealloc_xml_str(geom)
844 call dealloc_xml_str(att(1))
845 call dealloc_xml_str(col_grid)
846 end if
847 contains
848 !> \private assigns the 3D subarray pointer variables
849 subroutine assign_pointers(id)
850 ! input / output
851 integer :: id ! index at which to assing pointer
852
853 ! local variables
854 integer :: id_loc ! local id
855
856 ! X
857 if (present(x)) then
858 id_loc = id
859 if (col_id_loc.eq.1) then
860 if (size(x,1).eq.1) id_loc = 1
861 x_3d => x(id_loc,:,:,:)
862 else if (col_id_loc.eq.2) then
863 if (size(x,2).eq.1) id_loc = 1
864 x_3d => x(:,id_loc,:,:)
865 else if (col_id_loc.eq.3) then
866 if (size(x,3).eq.1) id_loc = 1
867 x_3d => x(:,:,id_loc,:)
868 else if (col_id_loc.eq.4) then
869 if (size(x,4).eq.1) id_loc = 1
870 x_3d => x(:,:,:,id_loc)
871 else
872 istat = 1
873 chckstt
874 end if
875 else
876 allocate(x_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
877 do jd = 1,loc_dim_3d(1)
878 x_3d(jd,:,:) = loc_offset_3d(1) + jd - 1
879 end do
880 end if
881 ! Y
882 if (present(y)) then
883 id_loc = id
884 if (col_id_loc.eq.1) then
885 if (size(y,1).eq.1) id_loc = 1
886 y_3d => y(id_loc,:,:,:)
887 else if (col_id_loc.eq.2) then
888 if (size(y,2).eq.1) id_loc = 1
889 y_3d => y(:,id_loc,:,:)
890 else if (col_id_loc.eq.3) then
891 if (size(y,3).eq.1) id_loc = 1
892 y_3d => y(:,:,id_loc,:)
893 else if (col_id_loc.eq.4) then
894 if (size(y,4).eq.1) id_loc = 1
895 y_3d => y(:,:,:,id_loc)
896 else
897 istat = 1
898 chckstt
899 end if
900 else
901 allocate(y_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
902 do jd = 1,loc_dim_3d(2)
903 y_3d(:,jd,:) = loc_offset_3d(2) + jd - 1
904 end do
905 end if
906 ! Z
907 if (present(z)) then
908 id_loc = id
909 if (col_id_loc.eq.1) then
910 if (size(z,1).eq.1) id_loc = 1
911 z_3d => z(id_loc,:,:,:)
912 else if (col_id_loc.eq.2) then
913 if (size(z,2).eq.1) id_loc = 1
914 z_3d => z(:,id_loc,:,:)
915 else if (col_id_loc.eq.3) then
916 if (size(z,3).eq.1) id_loc = 1
917 z_3d => z(:,:,id_loc,:)
918 else if (col_id_loc.eq.4) then
919 if (size(z,4).eq.1) id_loc = 1
920 z_3d => z(:,:,:,id_loc)
921 else
922 istat = 1
923 chckstt
924 end if
925 else
926 allocate(z_3d(loc_dim_3d(1),loc_dim_3d(2),loc_dim_3d(3)))
927 do jd = 1,loc_dim_3d(3)
928 z_3d(:,:,jd) = loc_offset_3d(3) + jd - 1
929 end do
930 end if
931 ! variable
932 if (col_id_loc.eq.1) then
933 var_3d => vars(id,:,:,:)
934 else if (col_id_loc.eq.2) then
935 var_3d => vars(:,id,:,:)
936 else if (col_id_loc.eq.3) then
937 var_3d => vars(:,:,id,:)
938 else if (col_id_loc.eq.4) then
939 var_3d => vars(:,:,:,id)
940 else
941 istat = 1
942 chckstt
943 end if
944 end subroutine assign_pointers
945 end subroutine plot_hdf5_arr
946 !> \private individual version
947 subroutine plot_hdf5_ind(var_name,file_name,var,tot_dim,loc_offset,&
948 &X,Y,Z,cont_plot,descr)
949
950 ! input / output
951 character(len=*), intent(in) :: var_name !< name of variable to be plot
952 character(len=*), intent(in) :: file_name !< file name
953 real(dp), intent(in) :: var(:,:,:) !< variable to plot
954 integer, intent(in), optional :: tot_dim(3) !< total dimensions of the arrays
955 integer, intent(in), optional :: loc_offset(3) !< offset of local dimensions
956 real(dp), intent(in), target, optional :: X(:,:,:) !< curvilinear grid X points
957 real(dp), intent(in), target, optional :: Y(:,:,:) !< curvilinear grid Y points
958 real(dp), intent(in), target, optional :: Z(:,:,:) !< curvilinear grid Z points
959 logical, intent(in), optional :: cont_plot !< continued plot
960 character(len=*), intent(in), optional :: descr !< description
961
962 ! local variables
963 integer :: tot_dim_loc(4), loc_offset_loc(4) ! local versions of total dimensions and local offset
964
965 ! set up local tot_dim and loc_offset
966 tot_dim_loc = [shape(var),1]
967 if (present(tot_dim)) tot_dim_loc = [tot_dim,1]
968 loc_offset_loc = 0
969 if (present(loc_offset)) loc_offset_loc = [loc_offset,0]
970
971 if (present(x)) then
972 if (present(y)) then
973 if (present(z)) then
974 call plot_hdf5_arr([var_name],file_name,&
975 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
976 &tot_dim_loc,loc_offset_loc,&
977 &x=reshape(x,[size(x,1),size(x,2),size(x,3),1]),&
978 &y=reshape(y,[size(y,1),size(y,2),size(y,3),1]),&
979 &z=reshape(z,[size(z,1),size(z,2),size(z,3),1]),&
980 &col=1,cont_plot=cont_plot,descr=descr)
981 else
982 call plot_hdf5_arr([var_name],file_name,&
983 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
984 &tot_dim_loc,loc_offset_loc,&
985 &x=reshape(x,[size(x,1),size(x,2),size(x,3),1]),&
986 &y=reshape(y,[size(y,1),size(y,2),size(y,3),1]),&
987 &col=1,cont_plot=cont_plot,descr=descr)
988 end if
989 else
990 if (present(z)) then
991 call plot_hdf5_arr([var_name],file_name,&
992 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
993 &tot_dim_loc,loc_offset_loc,&
994 &x=reshape(x,[size(x,1),size(x,2),size(x,3),1]),&
995 &z=reshape(z,[size(z,1),size(z,2),size(z,3),1]),&
996 &col=1,cont_plot=cont_plot,descr=descr)
997 else
998 call plot_hdf5_arr([var_name],file_name,&
999 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
1000 &tot_dim_loc,loc_offset_loc,&
1001 &x=reshape(x,[size(x,1),size(x,2),size(x,3),1]),&
1002 &col=1,cont_plot=cont_plot,descr=descr)
1003 end if
1004 end if
1005 else
1006 if (present(y)) then
1007 if (present(z)) then
1008 call plot_hdf5_arr([var_name],file_name,&
1009 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
1010 &tot_dim_loc,loc_offset_loc,&
1011 &y=reshape(y,[size(y,1),size(y,2),size(y,3),1]),&
1012 &z=reshape(z,[size(z,1),size(z,2),size(z,3),1]),&
1013 &col=1,cont_plot=cont_plot,descr=descr)
1014 else
1015 call plot_hdf5_arr([var_name],file_name,&
1016 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
1017 &tot_dim_loc,loc_offset_loc,&
1018 &y=reshape(y,[size(y,1),size(y,2),size(y,3),1]),&
1019 &col=1,cont_plot=cont_plot,descr=descr)
1020 end if
1021 else
1022 if (present(z)) then
1023 call plot_hdf5_arr([var_name],file_name,&
1024 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
1025 &tot_dim_loc,loc_offset_loc,&
1026 &col=1,z=reshape(z,[size(z,1),size(z,2),size(z,3),1]),&
1027 &cont_plot=cont_plot,descr=descr)
1028 else
1029 call plot_hdf5_arr([var_name],file_name,&
1030 &reshape(var,[size(var,1),size(var,2),size(var,3),1]),&
1031 &tot_dim_loc,loc_offset_loc,&
1032 &col=1,cont_plot=cont_plot,descr=descr)
1033 end if
1034 end if
1035 end if
1036 end subroutine plot_hdf5_ind
1037
1038 !> Use external program to draw a plot.
1039 !!
1040 !!
1041 !! The external programs used are determined by \c ex_plot_style:
1042 !! - \c ex_plot_style = 1: Use GnuPlot for 2-D and 3-D
1043 !! - \c ex_plot_style = 1: Use Bokeh in 2-D and Mayavi in 3-D
1044 !!
1045 !! The output is saved in a file <tt>[draw_name].[ext]</tt> where the
1046 !! extension depends on the plotting style and the dimension, or on the
1047 !! screen, depending on \c plot_on_screen.
1048 !!
1049 !! If not specified otherwise through data_name, it is assumed that the
1050 !! datafile to be plot is situated in <tt>[draw_name].dat</tt>.
1051 !!
1052 !! The title(s) of the plot(s) are provided through \c var_names and the
1053 !! number of plots through \c nplt. If there are less names than the number
1054 !! of plots, the first one is taken and appended by the plot number.
1055 !!
1056 !! Furthermore, \c draw_dim determines whether the plot is to be 2-D, 3-D or
1057 !! 2-D slices in 3-D.
1058 !!
1059 !! Also, an optional command \c extra_ops can be provided, to provide overal
1060 !! options, as well as draw_ops that specifies the line style for the plots
1061 !! from the file. If \c less draw_ops are provided than plots, they will be
1062 !! cycled.
1063 !!
1064 !! Finally, there is an option to provide animated plots, but is only
1065 !! available in 2-D. Also, for external plot style 1, no on-screen view is
1066 !! possible. Animations can take optional ranges arguments as well as delay.
1067 !! \note \c About draw_dim:
1068 !! - \c draw_dim = 1: 2-D plot; should be called with the output of
1069 !! output_ops.print_ex_2d(), \c nplt should be correctly set.
1070 !! - \c draw_dim = 2: 3-D plot; should be called with the output of
1071 !! print_ex_3D. For GNUPlot, \c nplt is ignored, but not for Mayavi, as it
1072 !! needs this information to be able to reconstruct 2D arrays for X, Y and
1073 !! Z.
1074 !! - \c draw_dim = 3: 2-D plot in 3-D slices; should be called with the
1075 !! output of output_ops.print_ex_2d(), \c nplt should be correctly set.
1076 subroutine draw_ex(var_names,draw_name,nplt,draw_dim,plot_on_screen,&
1077 &ex_plot_style,data_name,draw_ops,extra_ops,is_animated,ranges,delay,&
1078 &persistent)
1079
1080 use num_vars, only: ex_plot_style_global => ex_plot_style, rank
1081
1082 ! input / output
1083 character(len=*), intent(in) :: var_names(:) !< name of variables
1084 character(len=*), intent(in) :: draw_name !< name of drawing
1085 integer, intent(in) :: nplt !< number of plots
1086 integer, intent(in) :: draw_dim !< 1: 2-D, 2: 3-D, 3: decoupled 3-D
1087 logical, intent(in) :: plot_on_screen !< True if on screen, false if in file
1088 integer, intent(in), optional :: ex_plot_style !< alternative external plot style
1089 character(len=*), intent(in), optional :: data_name !< name of data file
1090 character(len=*), intent(in), optional :: draw_ops(:) !< drawing options
1091 character(len=*), intent(in), optional :: extra_ops !< extra options
1092 logical, intent(in), optional :: is_animated !< plot is animated
1093 real(dp), intent(in), optional :: ranges(:,:) !< x and y range of animated plot
1094 integer, intent(in), optional :: delay !< time delay between animated plot frames
1095 logical, intent(in), optional :: persistent !< keep on-screen plot open
1096
1097 ! local variables
1098 character(len=4) :: ex_ext ! output extension of external program
1099 character(len=5*max_str_ln) :: cmdmsg ! error message of C system command
1100 character(len=max_str_ln) :: script_name ! name of script, including path
1101 character(len=max_str_ln) :: ex_prog_name ! name of external program
1102 character(len=max_str_ln) :: data_name_loc ! local data name
1103 character(len=max_str_ln), allocatable :: var_names_loc(:) ! local variables names
1104 logical :: is_animated_loc ! local is_animated
1105 logical :: run_shell_necessary ! whether shell needs to be run
1106 logical :: persistent_loc ! local persistent
1107 real(dp), allocatable :: ranges_loc(:,:) ! local copy of ranges
1108 integer :: delay_loc ! local copy of delay
1109 integer :: n_draw_ops ! number of draw options provided
1110 integer :: istat ! status of opening a file
1111 integer :: iplt ! counter
1112 integer :: cmd_i ! file number for script file
1113 integer :: cmdstat ! status of C system command
1114 integer :: ex_plot_style_loc ! local ex_plot_style
1115
1116 ! skip if no plots
1117 if (no_plots) return
1118
1119 ! set up local data_name, variable names, is_animated and ex_plot_style
1120 data_name_loc = draw_name
1121 if (present(data_name)) data_name_loc = data_name
1122 allocate(var_names_loc(nplt))
1123 if (size(var_names).eq.nplt) then
1124 var_names_loc = var_names
1125 else
1126 do iplt = 1,nplt
1127 var_names_loc(iplt) = trim(var_names(1))//' ('//&
1128 &trim(i2str(iplt))//' / '//trim(i2str(nplt))//')'
1129 end do
1130 end if
1131 is_animated_loc = .false.
1132 if (present(is_animated)) is_animated_loc = is_animated
1133 ex_plot_style_loc = ex_plot_style_global
1134 if (present(ex_plot_style)) ex_plot_style_loc = ex_plot_style
1135 persistent_loc = .false.
1136 if (present(persistent)) persistent_loc = persistent
1137
1138 ! run shell to produce the plot by default
1139 run_shell_necessary = .true.
1140
1141 ! test whether animation is in 2D and set some variables
1142 if (is_animated_loc) then
1143 if (draw_dim.ne.1) then
1144 call writo('Animations are only possible in 2D',warning=.true.,&
1145 &persistent=.true.)
1146 return
1147 else
1148 ! calculate delay [1/100 s]
1149 if (present(delay)) then
1150 delay_loc = delay
1151 else
1152 delay_loc = max(250/nplt,10)
1153 end if
1154 end if
1155 end if
1156
1157 ! create the external program command
1158 if (plot_on_screen .and. &
1159 &(ex_plot_style_loc.ne.1 .and. is_animated_loc)) then ! not for GNUPlot animation
1160 script_name = trim(script_dir)//'/'//'temp_script_draw_ex_'//&
1161 &trim(i2str(rank))//'_'//trim(i2str(temp_id(2)))
1162 temp_id(2) = temp_id(2)+1
1163 else
1164 script_name = trim(script_dir)//'/'//trim(draw_name)
1165 end if
1166 select case (ex_plot_style_loc)
1167 case (1) ! GNUPlot
1168 script_name = trim(script_name)//'.gnu'
1169 ex_ext = 'pdf'
1170 ex_prog_name = 'gnuplot'
1171 case (2) ! Bokeh or Mayavi
1172 script_name = trim(script_name)//'.py'
1173 if (draw_dim.eq.1) then ! 2D: Bokeh
1174 ex_ext = 'html'
1175 else ! 3D and 2D slices in 3D: Mayavi
1176 ex_ext = 'png'
1177 end if
1178 ex_prog_name = 'python'
1179 end select
1180
1181 ! open script file
1182 open(nextunit(cmd_i),file=trim(script_name),iostat=istat)
1183 chckstt
1184
1185 ! set number of draw ops
1186 if (present(draw_ops)) then
1187 n_draw_ops = size(draw_ops)
1188 else
1189 n_draw_ops = 0
1190 end if
1191
1192 select case (ex_plot_style_loc)
1193 case (1) ! GNUPlot
1194 if (is_animated_loc) then
1195 call draw_ex_animated_gnuplot()
1196 else
1197 call draw_ex_gnuplot()
1198 end if
1199 case (2) ! Bokeh
1200 select case (draw_dim)
1201 case (1) ! 2D
1202 !!run_shell_necessary = plot_on_screen ! only need to run if plot on screen (?)
1203 if (is_animated_loc) then
1204 call draw_ex_animated_bokeh()
1205 else
1206 call draw_ex_bokeh()
1207 end if
1208 case (2:3) ! 3D, 2D slices in 3D
1209 call draw_ex_mayavi()
1210 end select
1211 end select
1212
1213 ! closing the command
1214 close(cmd_i)
1215
1216 ! run shell if necessary
1217 istat = 0
1218 cmdstat = 0
1219 if (run_shell_necessary) &
1220 &call use_execute_command_line(trim(ex_prog_name)//' "'//&
1221 &trim(script_name)//'"'//err_output_str,exitstat=istat,&
1222 &cmdstat=cmdstat,cmdmsg=cmdmsg)
1223
1224 if (istat.ne.0) then
1225 call writo('Failed to plot '//trim(draw_name)//'.'//&
1226 &trim(ex_ext),persistent=.true.,warning=.true.)
1227 else
1228 if (cmdstat.ne.0) then
1229 call writo('Failed to plot '//trim(draw_name)//'.'//&
1230 &trim(ex_ext),persistent=.true.,warning=.true.)
1231 call lvl_ud(1)
1232 call writo('System message: "'//trim(cmdmsg)//'"',&
1233 &persistent=.true.)
1234 if (.not.plot_on_screen) call writo(&
1235 &'Try running "'//trim(ex_prog_name)//' "'//&
1236 &trim(script_name)//'"'//'" manually',persistent=.true.)
1237 call lvl_ud(-1)
1238 else
1239 if (plot_on_screen .and. run_shell_necessary .and. &
1240 &.not.persistent_loc) then
1241 call use_execute_command_line('rm '//&
1242 &trim(script_name),exitstat=istat,&
1243 &cmdstat=cmdstat,cmdmsg=cmdmsg)
1244 ! ignore errors
1245 else
1246 call writo('Plot in output file "'//&
1247 &trim(plot_dir)//'/'//trim(draw_name)//&
1248 &'.'//trim(ex_ext)//'"',persistent=.true.)
1249 end if
1250 end if
1251 end if
1252 contains
1253 !> \private GNUPlot version: wxt terminal or pdf output
1254 subroutine draw_ex_gnuplot
1255 ! initialize the script
1256 write(cmd_i,"(A)",iostat=istat) 'set grid'
1257 write(cmd_i,"(A)",iostat=istat) 'set border 4095 front linetype -1 &
1258 &linewidth 1.0'
1259 if (plot_on_screen) then
1260 if (persistent_loc) then
1261 write(cmd_i,"(A)",iostat=istat) 'set terminal wxt persist'
1262 else
1263 write(cmd_i,"(A)",iostat=istat) 'set terminal wxt'
1264 end if
1265 else
1266 write(cmd_i,"(A)",iostat=istat) 'set terminal pdf size '//&
1267 &trim(i2str(plot_size(1)))//','//trim(i2str(plot_size(2)))
1268 write(cmd_i,"(A)",iostat=istat) 'set output "'//&
1269 &trim(plot_dir)//'/'//trim(draw_name)//'.pdf"'
1270 end if
1271
1272 ! no legend if too many plots
1273 if (nplt.gt.10) write(cmd_i,"(A)",iostat=istat) 'set nokey;'
1274
1275 ! write extra options
1276 if (present(extra_ops)) write(cmd_i,"(A)",iostat=istat) &
1277 &trim(extra_ops)
1278
1279 ! set up line styles
1280 if (n_draw_ops.eq.0) then
1281 do iplt = 1,min(nplt,size(line_clrs))
1282 write(cmd_i,"(A)",iostat=istat) 'set style line '//&
1283 &trim(i2str(iplt))//' lc rgb '//&
1284 &line_clrs(iplt)//' '//trim(line_style)
1285 end do
1286 end if
1287
1288 ! individual plots
1289 select case (draw_dim)
1290 case (1) ! 2D
1291 write(cmd_i,"(A)",iostat=istat) 'plot \'
1292 do iplt = 1,nplt
1293 write(cmd_i,"(A)",iostat=istat,advance="no") &
1294 &' "'//trim(data_dir)//'/'//trim(data_name_loc)//&
1295 &'.dat" using '//trim(i2str(iplt))//':'//&
1296 &trim(i2str(nplt+iplt))//' title "'//&
1297 &trim(var_names_loc(iplt))//'" '//&
1298 &trim(loc_draw_op())
1299 if (iplt.eq.nplt) then
1300 write(cmd_i,"(A)",iostat=istat) ''
1301 else
1302 write(cmd_i,"(A)",iostat=istat) ', \'
1303 end if
1304 end do
1305 case (2) ! 3D
1306 write(cmd_i,"(A)",iostat=istat) 'splot \'
1307 do iplt = 1,nplt
1308 write(cmd_i,"(A)",iostat=istat,advance="no") &
1309 &' "'//trim(data_dir)//'/'//trim(data_name_loc)//&
1310 &'.dat" using '//trim(i2str(iplt))//':'//&
1311 &trim(i2str(nplt+iplt))//':'//&
1312 &trim(i2str(2*nplt+iplt))//' title "'//&
1313 &trim(var_names_loc(iplt))//&
1314 &'" '//trim(loc_draw_op())
1315 if (iplt.eq.nplt) then
1316 write(cmd_i,"(A)",iostat=istat) ''
1317 else
1318 write(cmd_i,"(A)",iostat=istat) ', \'
1319 end if
1320 end do
1321 case (3) ! 2D slices in 3D
1322 write(cmd_i,"(A)",iostat=istat) 'splot \'
1323 do iplt = 1,nplt
1324 write(cmd_i,"(A)",iostat=istat,advance="no") &
1325 &' "'//trim(data_dir)//'/'//trim(data_name_loc)//&
1326 &'.dat" using ('//trim(i2str(iplt))//'):'//&
1327 &trim(i2str(iplt))//':'//trim(i2str(nplt+iplt))//&
1328 &' title "'//trim(var_names_loc(iplt))//'" '//&
1329 &trim(loc_draw_op())
1330 if (iplt.eq.nplt) then
1331 write(cmd_i,"(A)",iostat=istat) ''
1332 else
1333 write(cmd_i,"(A)",iostat=istat) ', \'
1334 end if
1335 end do
1336 case default
1337 call writo('No draw_dim associated with '//&
1338 &trim(i2str(draw_dim)),persistent=.true.)
1339 istat = 1
1340 chckstt
1341 end select
1342
1343 ! finishing the command
1344 write(cmd_i,"(A)",iostat=istat) ''
1345 if (plot_on_screen) write(cmd_i,"(A)",iostat=istat) 'pause -1'
1346 end subroutine draw_ex_gnuplot
1347
1348 !> \private Bokeh version: 2D html output
1349 ! Interactive checkbox from https://github.com/bokeh/bokeh/issues/3715.
1350 ! (Could be replaced by http://bokeh.pydata.org/en/latest/docs/
1351 ! user_guide/interaction/legends.html#userguide-interaction-legends)
1352 subroutine draw_ex_bokeh
1353 ! initialize the script
1354 write(cmd_i,"(A)",iostat=istat) 'from numpy import genfromtxt'
1355 write(cmd_i,"(A)",iostat=istat) 'from bokeh.plotting import &
1356 &figure, output_file, save, show'
1357 write(cmd_i,"(A)",iostat=istat) 'from bokeh.layouts import &
1358 &widgetbox, row'
1359 write(cmd_i,"(A)",iostat=istat) 'from bokeh.models.widgets import &
1360 &CheckboxGroup'
1361 write(cmd_i,"(A)",iostat=istat) 'from bokeh.models import CustomJS'
1362 write(cmd_i,"(A)",iostat=istat) ''
1363 write(cmd_i,"(A)",iostat=istat) 'TOOLS="hover,pan,wheel_zoom,&
1364 &box_zoom,reset,save"'
1365 write(cmd_i,"(A)",iostat=istat) ''
1366 write(cmd_i,"(A)",iostat=istat) 'output_file("'//trim(plot_dir)//&
1367 &'/'//trim(draw_name)//'.html", title="'//trim(draw_name)//&
1368 &'", mode="cdn") '
1369 write(cmd_i,"(A)",iostat=istat) ''
1370 write(cmd_i,"(A)",iostat=istat) 'data = genfromtxt("'//&
1371 &trim(data_dir)//'/'//trim(data_name_loc)//'.dat")'
1372 write(cmd_i,"(A)",iostat=istat) ''
1373 write(cmd_i,"(A)",iostat=istat) 'p = figure(toolbar_location=&
1374 &"above",tools=TOOLS,active_drag="box_zoom",active_scroll=&
1375 &"wheel_zoom",plot_width=600, plot_height=600)'
1376 write(cmd_i,"(A)",iostat=istat) ''
1377
1378 ! write extra options
1379 if (present(extra_ops)) write(cmd_i,"(A)",iostat=istat) &
1380 &trim(extra_ops)
1381
1382 ! individual plots
1383 do iplt = 1,nplt
1384 ! plot the lines
1385 write(cmd_i,"(A)",iostat=istat) 'l'//trim(i2str(iplt))//&
1386 &' = p.line(data[:,'//trim(i2str(iplt-1))//'],data[:,'//&
1387 &trim(i2str(nplt+iplt-1))//'],'//trim(loc_draw_op(1))//')'
1388 ! plot the circles
1389 write(cmd_i,"(A)",iostat=istat) 'c'//trim(i2str(iplt))//&
1390 &' = p.circle(data[:,'//&
1391 &trim(i2str(iplt-1))//'],data[:,'//&
1392 &trim(i2str(nplt+iplt-1))//'],'//trim(loc_draw_op(2))//')'
1393 end do
1394 write(cmd_i,"(A)",iostat=istat) ''
1395
1396 ! create checkbox
1397 write(cmd_i,"(A)",iostat=istat) 'labels=list()'
1398 do iplt = 1,nplt
1399 write(cmd_i,"(A)",iostat=istat) 'labels.append("'//&
1400 &trim(var_names_loc(iplt))//'")'
1401 end do
1402 write(cmd_i,"(A)",iostat=istat) 'active=list()'
1403 do iplt = 1,nplt
1404 write(cmd_i,"(A)",iostat=istat) 'active.append('//&
1405 &trim(i2str(iplt-1))//')'
1406 end do
1407 write(cmd_i,"(A)",iostat=istat) 'checkbox = CheckboxGroup(&
1408 &labels=labels, active=active)'
1409 write(cmd_i,"(A)",iostat=istat) ''
1410
1411 ! create callback
1412 write(cmd_i,"(A)",iostat=istat) 'args=dict()'
1413 do iplt = 1,nplt
1414 write(cmd_i,"(A)",iostat=istat) 'args["l'//&
1415 &trim(i2str(iplt))//'"] = l'//trim(i2str(iplt))
1416 write(cmd_i,"(A)",iostat=istat) 'args["c'//&
1417 &trim(i2str(iplt))//'"] = c'//trim(i2str(iplt))
1418 end do
1419
1420 write(cmd_i,"(A)",iostat=istat) 'code="""'
1421 write(cmd_i,"(A)",iostat=istat) '//console.log(cb_obj.active);'
1422 do iplt = 1,nplt
1423 write(cmd_i,"(A)",iostat=istat) 'l'//trim(i2str(iplt))//&
1424 &'.visible = false;'
1425 write(cmd_i,"(A)",iostat=istat) 'c'//trim(i2str(iplt))//&
1426 &'.visible = false;'
1427 end do
1428 write(cmd_i,"(A)",iostat=istat) 'for (i in cb_obj.active) {'
1429 write(cmd_i,"(A)",iostat=istat) &
1430 &' //console.log(cb_obj.active[i]);'
1431 do iplt = 1,nplt
1432 write(cmd_i,"(A)",iostat=istat) ' if (cb_obj.active[i] == '&
1433 &//trim(i2str(iplt-1))//') {'
1434 write(cmd_i,"(A)",iostat=istat) ' l'//trim(i2str(iplt))//&
1435 &'.visible = true;'
1436 write(cmd_i,"(A)",iostat=istat) ' c'//trim(i2str(iplt))//&
1437 &'.visible = true;'
1438 write(cmd_i,"(A)",iostat=istat) ' }'
1439 end do
1440 write(cmd_i,"(A)",iostat=istat) '}'
1441 write(cmd_i,"(A)",iostat=istat) '"""'
1442
1443 write(cmd_i,"(A)",iostat=istat) 'checkbox.callback = CustomJS(&
1444 &args=args, code=code)'
1445
1446 ! create layout
1447 write(cmd_i,"(A)",iostat=istat) &
1448 &'layout = row(p, widgetbox(checkbox))'
1449 write(cmd_i,"(A)",iostat=istat) ''
1450
1451 ! finishing the command
1452 write(cmd_i,"(A)",iostat=istat) ''
1453 if (plot_on_screen) then
1454 write(cmd_i,"(A)",iostat=istat) 'show(layout)'
1455 else
1456 write(cmd_i,"(A)",iostat=istat) 'save(layout)'
1457 end if
1458 end subroutine draw_ex_bokeh
1459
1460 !> \private Mayavi version: 3D png output
1461 subroutine draw_ex_mayavi
1462 ! initialize the script
1463 write(cmd_i,"(A)",iostat=istat) 'from numpy import genfromtxt, &
1464 &array, size, zeros, amin, amax'
1465 write(cmd_i,"(A)",iostat=istat) 'from mayavi.mlab import outline, &
1466 &mesh, points3d, savefig, colorbar, axes, savefig, show'
1467 write(cmd_i,"(A)",iostat=istat) ''
1468 write(cmd_i,"(A)",iostat=istat) 'nplt = '//trim(i2str(nplt))
1469
1470 ! calculate the number of y points
1471 write(cmd_i,"(A)",iostat=istat) 'npnty = -1'
1472 write(cmd_i,"(A)",iostat=istat) 'with open("'//trim(data_dir)//&
1473 &'/'//trim(data_name_loc)//'.dat") as f:'
1474 write(cmd_i,"(A)",iostat=istat) ' for line in f:'
1475 write(cmd_i,"(A)",iostat=istat) ' if (npnty < 0):'
1476 write(cmd_i,"(A)",iostat=istat) ' if &
1477 &(line.strip()[0] == "#"):'
1478 write(cmd_i,"(A)",iostat=istat) ' npnty = 0'
1479 write(cmd_i,"(A)",iostat=istat) ' else:'
1480 write(cmd_i,"(A)",iostat=istat) ' if (not line.strip()):'
1481 write(cmd_i,"(A)",iostat=istat) ' break'
1482 write(cmd_i,"(A)",iostat=istat) ' else:'
1483 write(cmd_i,"(A)",iostat=istat) ' npnty += 1'
1484
1485 ! get data
1486 write(cmd_i,"(A)",iostat=istat) 'data = genfromtxt("'//&
1487 &trim(data_dir)//'/'//trim(data_name_loc)//'.dat")'
1488 write(cmd_i,"(A)",iostat=istat) 'dims = &
1489 &array([size(data,0)/npnty,npnty,nplt])'
1490 write(cmd_i,"(A)",iostat=istat) ''
1491 write(cmd_i,"(A)",iostat=istat) 'X = zeros(dims)'
1492 write(cmd_i,"(A)",iostat=istat) 'Y = zeros(dims)'
1493 write(cmd_i,"(A)",iostat=istat) 'Z = zeros(dims)'
1494 write(cmd_i,"(A)",iostat=istat) 'for i in range(0,dims[0]):'
1495 write(cmd_i,"(A)",iostat=istat) ' for j in range(0,dims[2]):'
1496 write(cmd_i,"(A)",iostat=istat) ' X[i,:,j] = &
1497 &data[npnty*i:npnty*(i+1),j]'
1498 write(cmd_i,"(A)",iostat=istat) ' Y[i,:,j] = &
1499 &data[npnty*i:npnty*(i+1),nplt+j]'
1500 write(cmd_i,"(A)",iostat=istat) ' Z[i,:,j] = &
1501 &data[npnty*i:npnty*(i+1),nplt*2+j]'
1502 write(cmd_i,"(A)",iostat=istat) ''
1503 write(cmd_i,"(A)",iostat=istat) 'minZ = amin(Z)'
1504 write(cmd_i,"(A)",iostat=istat) 'maxZ = amax(Z)'
1505 write(cmd_i,"(A)",iostat=istat) ''
1506
1507 ! write extra options
1508 if (present(extra_ops)) write(cmd_i,"(A)",iostat=istat) &
1509 &trim(extra_ops)
1510
1511 ! plot
1512 select case (draw_dim)
1513 case (2) ! 3D
1514 write(cmd_i,"(A)",iostat=istat) 'for j in range(0,dims[2]):'
1515 write(cmd_i,"(A)",iostat=istat) ' &
1516 &mesh(X[:,:,j],Y[:,:,j],Z[:,:,j],vmin=minZ,vmax=maxZ)'
1517 case (3) ! 2D slices in 3D
1518 write(cmd_i,"(A)",iostat=istat) 'for j in range(0,dims[2]):'
1519 write(cmd_i,"(A)",iostat=istat) ' &
1520 &points3d(X[:,:,j],Y[:,:,j],Z[:,:,j],mode="sphere",&
1521 &scale_factor=0.05,vmin=minZ,vmax=maxZ)'
1522 end select
1523 write(cmd_i,"(A)",iostat=istat) 'outline()'
1524 write(cmd_i,"(A)",iostat=istat) 'colorbar()'
1525 write(cmd_i,"(A)",iostat=istat) 'axes()'
1526
1527 ! pause on screen or output to file
1528 if (plot_on_screen) then
1529 write(cmd_i,"(A)",iostat=istat) ''
1530 write(cmd_i,"(A)",iostat=istat) 'show()'
1531 else
1532 write(cmd_i,"(A)",iostat=istat) ''
1533 write(cmd_i,"(A)",iostat=istat) '#show()'
1534 write(cmd_i,"(A)",iostat=istat) 'savefig("'//trim(plot_dir)//&
1535 &'/'//trim(draw_name)//'.png"'//',magnification=4)'
1536 end if
1537 end subroutine draw_ex_mayavi
1538
1539 !> \private GNUPlot animated version: gif output
1540 subroutine draw_ex_animated_gnuplot
1541 ! initialize the script
1542 write(cmd_i,"(A)",iostat=istat) 'set grid'
1543 write(cmd_i,"(A)",iostat=istat) 'set border 4095 front linetype -1 &
1544 &linewidth 1.0'
1545 write(cmd_i,"(A)",iostat=istat) 'set terminal gif animate &
1546 &delay '//trim(i2str(delay_loc))//' size '//&
1547 &trim(i2str(128*plot_size(1)))//','//&
1548 &trim(i2str(128*plot_size(2)))
1549 write(cmd_i,"(A)",iostat=istat) 'set output "'//trim(plot_dir)//&
1550 &'/'//trim(draw_name)//'.gif"'
1551
1552 ! find and set ranges
1553 allocate(ranges_loc(2,2))
1554 if (present(ranges)) then
1555 if (size(ranges,1).eq.2 .and. size(ranges,2).eq.2) then
1556 ranges_loc = ranges
1557 else
1558 call writo('invalid ranges given &
1559 &to draw_ex_animated',persistent=.true.,&
1560 &warning=.true.)
1561 end if
1562 else
1563 ! initialize ranges
1564 ranges_loc(:,1) = huge(1._dp) ! minimum value
1565 ranges_loc(:,2) = -huge(1._dp) ! maximum value
1566
1567 call get_ranges(ranges_loc)
1568 end if
1569
1570 ! set ranges
1571 write(cmd_i,"(A)",iostat=istat) 'set xrange ['//&
1572 &trim(r2str(ranges_loc(1,1)))//':'//&
1573 &trim(r2str(ranges_loc(1,2)))//'];'
1574 write(cmd_i,"(A)",iostat=istat) 'set yrange ['//&
1575 &trim(r2str(ranges_loc(2,1)))//':'//&
1576 &trim(r2str(ranges_loc(2,2)))//'];'
1577
1578 ! write extra options
1579 if (present(extra_ops)) write(cmd_i,"(A)",iostat=istat) &
1580 &trim(extra_ops)
1581
1582 ! set up line styles
1583 if (n_draw_ops.eq.0) then
1584 do iplt = 1,min(nplt,size(line_clrs))
1585 write(cmd_i,"(A)",iostat=istat) 'set style line '//&
1586 &trim(i2str(iplt))//' lc rgb '//&
1587 &line_clrs(iplt)//' '//trim(line_style)
1588 end do
1589 end if
1590
1591 ! individual plots
1592 do iplt = 1,nplt
1593 write(cmd_i,"(A)",iostat=istat) 'plot "'//trim(data_dir)//'/'//&
1594 &trim(data_name_loc)//'.dat" using '//trim(i2str(iplt))//&
1595 &':'//trim(i2str(nplt+iplt))//' title "'//&
1596 &trim(var_names_loc(iplt))//'" '//&
1597 &trim(loc_draw_op())
1598 end do
1599
1600 ! finishing the command
1601 write(cmd_i,"(A)",iostat=istat) ''
1602 end subroutine draw_ex_animated_gnuplot
1603
1604 !> \private Bokeh animated version: html output
1605 subroutine draw_ex_animated_bokeh
1606 ! initialize the script
1607 write(cmd_i,"(A)",iostat=istat) '# Note: it is necessary to first &
1608 &run "bokeh serve" before calling python'
1609 write(cmd_i,"(A)",iostat=istat) ''
1610 write(cmd_i,"(A)",iostat=istat) 'from numpy import genfromtxt'
1611 write(cmd_i,"(A)",iostat=istat) 'from bokeh.plotting import &
1612 &figure, output_file, save, show, curdoc'
1613 write(cmd_i,"(A)",iostat=istat) 'from bokeh.client import &
1614 &push_session'
1615 write(cmd_i,"(A)",iostat=istat) 'from bokeh.driving import repeat'
1616 write(cmd_i,"(A)",iostat=istat) ''
1617 write(cmd_i,"(A)",iostat=istat) 'TOOLS="crosshair,pan,wheel_zoom,&
1618 &box_zoom,reset,save"'
1619 write(cmd_i,"(A)",iostat=istat) ''
1620 write(cmd_i,"(A)",iostat=istat) 'output_file("'//trim(plot_dir)//&
1621 &'/'//trim(draw_name)//'.html", title="'//trim(draw_name)//&
1622 &'", mode="cdn") '
1623 write(cmd_i,"(A)",iostat=istat) ''
1624 write(cmd_i,"(A)",iostat=istat) 'data = genfromtxt("'//&
1625 &trim(data_dir)//'/'//trim(data_name_loc)//'.dat")'
1626 write(cmd_i,"(A)",iostat=istat) ''
1627 write(cmd_i,"(A)",iostat=istat) 'p = figure(toolbar_location=&
1628 &"above",'//'tools=TOOLS,active_drag="pan",active_scroll=&
1629 &"wheel_zoom",'//'plot_width=600, plot_height=600)'
1630 write(cmd_i,"(A)",iostat=istat) ''
1631
1632 ! write extra options
1633 if (present(extra_ops)) write(cmd_i,"(A)",iostat=istat) &
1634 &trim(extra_ops)
1635
1636 ! individual plot for first plot
1637 write(cmd_i,"(A)",iostat=istat) 'line = p.line(data[:,0],data[:,'//&
1638 &trim(i2str(nplt))//'],'//trim(loc_draw_op(1))//',legend="'//&
1639 &trim(var_names_loc(1))//'")'
1640 ! plot the circles
1641 write(cmd_i,"(A)",iostat=istat) 'circle = p.circle(data[:,'//&
1642 &trim(i2str(0))//'],data[:,'//trim(i2str(nplt))//'],'//&
1643 &trim(loc_draw_op(2))//')'
1644
1645 ! finishing the command (needs to use show)
1646 write(cmd_i,"(A)",iostat=istat) ''
1647 write(cmd_i,"(A)",iostat=istat) 'session = push_session(curdoc())'
1648 write(cmd_i,"(A)",iostat=istat) 'session.show(p)'
1649 write(cmd_i,"(A)",iostat=istat) ''
1650 write(cmd_i,"(A)",iostat=istat) 'seq = range(0,'//&
1651 &trim(i2str(nplt))//')'
1652 write(cmd_i,"(A)",iostat=istat) '@repeat(seq)'
1653 write(cmd_i,"(A)",iostat=istat) 'def update(i):'
1654 write(cmd_i,"(A)",iostat=istat) ' line.data_source.data["x"] = &
1655 &data[:,i]'
1656 write(cmd_i,"(A)",iostat=istat) ' line.data_source.data["y"] = &
1657 &data[:,'//trim(i2str(nplt))//'+i]'
1658 write(cmd_i,"(A)",iostat=istat) ' circle.data_source.data["x"] = &
1659 &data[:,i]'
1660 write(cmd_i,"(A)",iostat=istat) ' circle.data_source.data["y"] = &
1661 &data[:,'//&
1662 &trim(i2str(nplt))//'+i]'
1663 write(cmd_i,"(A)",iostat=istat) ''
1664 write(cmd_i,"(A)",iostat=istat) 'curdoc().add_periodic_callback(&
1665 &update, '//trim(i2str(delay_loc*10))//')'
1666 write(cmd_i,"(A)",iostat=istat) 'session.loop_until_closed()'
1667 end subroutine draw_ex_animated_bokeh
1668
1669 !> \private gets ranges for animated plot
1670 subroutine get_ranges(ranges)
1671 ! input / output
1672 real(dp), intent(inout) :: ranges(:,:) ! x and y range, and z range (if 3D) of plot
1673
1674 ! local variables
1675 integer :: data_i ! file number of data file
1676 real(dp), allocatable :: loc_data(:) ! one line of data
1677 character(len=1) :: loc_data_char ! local first char
1678
1679
1680 ! open data file
1681 open(nextunit(data_i),file=data_dir//'/'//trim(data_name_loc)//&
1682 &'.dat',status='old',iostat=istat)
1683
1684 if (istat.eq.0) then
1685 ! set up loc_data
1686 allocate(loc_data(2*nplt))
1687
1688 ! read the data file
1689 istat = 0
1690 do while (istat.eq.0)
1691 read(data_i,*,iostat=istat) loc_data_char ! read first character of data
1692 if (istat.eq.0) then ! read succesful
1693 if (loc_data_char.ne.'#') then ! exclude comment lines
1694 backspace(unit=data_i) ! go back one line
1695 read(data_i,*,iostat=istat) loc_data ! read data again, but now ful
1696 ranges(1,1) = &
1697 &min(ranges(1,1),minval(loc_data(1:nplt)))
1698 ranges(1,2) = &
1699 &max(ranges(1,2),maxval(loc_data(1:nplt)))
1700 ranges(2,1) = &
1701 &min(ranges(2,1),&
1702 &minval(loc_data(nplt+1:2*nplt)))
1703 ranges(2,2) = &
1704 &max(ranges(2,2),&
1705 &maxval(loc_data(nplt+1:2*nplt)))
1706 end if
1707 end if
1708 end do
1709 !write(*,*) 'ranges(1,:) = ', ranges(1,:)
1710 !write(*,*) 'ranges(2,:) = ', ranges(2,:)
1711 !write(*,*) 'ranges(3,:) = ', ranges(3,:)
1712
1713 ! close data file
1714 close(data_i,iostat=istat)
1715 end if
1716 end subroutine get_ranges
1717
1718 !> \private sets local draw options either from pre-defined line styles or through
1719 ! user specified option, for GNUPlot or Bokeh
1720 ! As for Bokeh, there are two default draw options possible (for the
1721 ! lines and for the points), there is an option to select between them.
1722 function loc_draw_op(Bokeh_style)
1723 ! input / output
1724 character(len=max_str_ln) :: loc_draw_op ! local drawing option
1725 integer, intent(in), optional :: bokeh_style ! lines(1) or points (2)
1726
1727 select case (ex_plot_style_loc)
1728 case (1) ! GNUPlot
1729 if (n_draw_ops.gt.0) then
1730 loc_draw_op = trim(draw_ops(mod(iplt-1,n_draw_ops)+1))
1731 else
1732 loc_draw_op = 'with linespoints linestyle '//&
1733 &trim(i2str(mod(iplt-1,size(line_clrs))+1))
1734 end if
1735 case (2) ! Bokeh / Mayavi
1736 if (n_draw_ops.gt.0) then
1737 loc_draw_op = trim(draw_ops(mod(iplt-1,n_draw_ops)+1))
1738 else
1739 select case (bokeh_style)
1740 case (1) ! lines
1741 loc_draw_op = 'line_color='//&
1742 &line_clrs(mod(iplt-1,size(line_clrs))+1)//&
1743 &',line_width=2'
1744 case (2) ! points
1745 loc_draw_op = 'color='//&
1746 &line_clrs(mod(iplt-1,size(line_clrs))+1)//&
1747 &',size=5'
1748 case default
1749 loc_draw_op = ''
1750 end select
1751 end if
1752 end select
1753 end function loc_draw_op
1754 end subroutine draw_ex
1755
1756 !> Takes two input vectors and plots these as well as the relative and
1757 !! absolute difference in a HDF5 file.
1758 !!
1759 !! This is similar to a basic version of output_ops.plot_hdf5().
1760 !!
1761 !! Optionally, an output message can be displayed on screen with the maximum
1762 !! relative and absolute error.
1763 subroutine plot_diff_hdf5(A,B,file_name,tot_dim,loc_offset,descr,&
1764 &output_message)
1765
1766 ! input / output
1767 real(dp), intent(in) :: a(:,:,:) !< vector A
1768 real(dp), intent(in) :: b(:,:,:) !< vector B
1769 character(len=*), intent(in) :: file_name !< name of plot
1770 integer, intent(in), optional :: tot_dim(3) !< total dimensions of the arrays
1771 integer, intent(in), optional :: loc_offset(3) !< offset of local dimensions
1772 character(len=*), intent(in), optional :: descr !< description
1773 logical, intent(in), optional :: output_message !< whether to display a message or not
1774
1775 ! local variables
1776 real(dp), allocatable :: plot_var(:,:,:,:) ! variable containing plot
1777 integer :: tot_dim_loc(4) ! local version of tot_dim
1778 integer :: loc_offset_loc(4) ! local version of loc_offset
1779 character(len=max_str_ln) :: var_names(5) ! names of variables in plot
1780 logical :: output_message_loc ! local version of output_message
1781 real(dp) :: lim_lo ! lower limit of errors
1782 real(dp) :: lim_hi ! upper limit of errors
1783 real(dp) :: err_av ! average error
1784 real(dp), allocatable :: tot_lim(:) ! total lower or upper limit
1785 real(dp), allocatable :: tot_err(:) ! total sum of errors
1786 integer :: istat ! status
1787 logical :: ind_plot ! individual plot or not
1788
1789 ! set up local tot_dim and loc_offset
1790 tot_dim_loc = [shape(a),5]
1791 if (present(tot_dim)) tot_dim_loc = [tot_dim,5]
1792 loc_offset_loc = 0
1793 if (present(loc_offset)) loc_offset_loc = [loc_offset,5]
1794
1795 ! tests
1796 if (size(a,1).ne.size(b,1) .or. size(a,2).ne.size(b,2) .or. &
1797 &size(a,3).ne.size(b,3)) then
1798 call writo('in plot_diff_HDF5, A and B need to have the correct &
1799 &size',persistent=.true.,warning=.true.)
1800 return
1801 end if
1802
1803 ! set up local ouput message
1804 output_message_loc = .false.
1805 if (present(output_message)) output_message_loc = output_message
1806
1807 ind_plot = .false.
1808 if (tot_dim_loc(1).eq.size(a,1) .and. tot_dim_loc(2).eq.size(a,2) &
1809 &.and. tot_dim_loc(3).eq.size(a,3)) ind_plot = .true.
1810
1811 ! set up plot_var
1812 allocate(plot_var(size(a,1),size(a,2),size(a,3),5))
1813 plot_var(:,:,:,1) = a
1814 plot_var(:,:,:,2) = b
1815 plot_var(:,:,:,3) = diff(a,b,shape(a),rel=.true.) ! rel. diff.
1816 plot_var(:,:,:,4) = log10(abs(diff(a,b,shape(a),rel=.true.))) ! log of abs. value of rel. diff.
1817 plot_var(:,:,:,5) = diff(a,b,shape(a),rel=.false.) ! abs. diff.
1818
1819 ! set up var_names
1820 var_names(1) = 'v1'
1821 var_names(2) = 'v2'
1822 var_names(3) = 'rel v1 - v2'
1823 var_names(4) = 'log abs rel v1 - v2'
1824 var_names(5) = 'abs v1 - v2'
1825
1826 ! plot
1827 call plot_hdf5_arr(var_names,file_name,plot_var,tot_dim=tot_dim_loc,&
1828 &loc_offset=loc_offset_loc,col=1,descr=descr)
1829
1830 ! output message if requested
1831 if (output_message_loc) then
1832 call writo('Information about errors:',persistent=.true.)
1833 call lvl_ud(1)
1834 ! relative error
1835 call stats(tot_dim_loc(1:3),plot_var(:,:,:,3),lim_lo,lim_hi,err_av)
1836 call writo(trim(r2strt(lim_lo))//' < rel. err. < '//&
1837 &trim(r2strt(lim_hi))//', average value: '//&
1838 &trim(r2strt(err_av)),persistent=.true.)
1839 ! log of absolute relative error
1840 call stats(tot_dim_loc(1:3),plot_var(:,:,:,4),lim_lo,lim_hi,err_av)
1841 call writo(trim(r2strt(lim_lo))//' < log(abs(rel. err.)) < '//&
1842 &trim(r2strt(lim_hi))//', average value: '//&
1843 &trim(r2strt(err_av)),persistent=.true.)
1844 ! absolute error
1845 call stats(tot_dim_loc(1:3),plot_var(:,:,:,5),lim_lo,lim_hi,err_av)
1846 call writo(trim(r2strt(lim_lo))//' < abs. err. < '//&
1847 &trim(r2strt(lim_hi))//', average value: '//&
1848 &trim(r2strt(err_av)),persistent=.true.)
1849 call lvl_ud(-1)
1850 end if
1851 contains
1852 ! returns relative or absolute difference between inputs A and B
1853 !> \private
1854 function diff(A,B,dims,rel) result(C)
1855 use num_vars, only: tol_zero
1856
1857 ! input / output
1858 real(dp), intent(in) :: a(:,:,:) ! input A
1859 real(dp), intent(in) :: b(:,:,:) ! input B
1860 integer, intent(in) :: dims(3) ! dimensions of A and B
1861 logical, intent(in) :: rel ! .true. if relative and .false. if absolute error
1862 real(dp) :: c(dims(1),dims(2),dims(3)) ! output C
1863
1864 ! return output
1865 if (rel) then
1866 c = 2*(a-b)/max(tol_zero,(abs(a)+abs(b)))
1867 else
1868 c = abs(a-b)
1869 end if
1870 end function diff
1871
1872 ! returns limits and average value
1873 !> \private
1874 subroutine stats(tot_dims,var,lim_lo,lim_hi,err_av)
1875 use mpi_utilities, only: get_ser_var
1876
1877 ! input / output
1878 integer, intent(in) :: tot_dims(3) ! total dimensions of grid
1879 real(dp), intent(in) :: var(:,:,:) ! input variable for which to calculate statistics
1880 real(dp), intent(inout) :: lim_lo ! lower limit of errors
1881 real(dp), intent(inout) :: lim_hi ! upper limit of errors
1882 real(dp), intent(inout) :: err_av ! average error
1883
1884 ! local variables
1885 real(dp) :: sum_err ! sum of errors
1886
1887 ! calculate limits on absolute error
1888 lim_lo = minval(var)
1889 lim_hi = maxval(var)
1890 sum_err = sum(var)
1891
1892 ! get most stringent limits from all processes
1893 if (.not.ind_plot) then
1894 istat = get_ser_var([lim_lo],tot_lim)
1895 chckstt
1896 lim_lo = minval(tot_lim)
1897 istat = get_ser_var([lim_hi],tot_lim)
1898 chckstt
1899 lim_hi = maxval(tot_lim)
1900 istat = get_ser_var([sum_err],tot_err)
1901 chckstt
1902 sum_err = sum(tot_err)
1903 end if
1904
1905 ! calculate average error
1906 err_av = sum_err/product(tot_dims)
1907 end subroutine
1908 end subroutine plot_diff_hdf5
1909
1910 !> Executes command line, or displays a message if disabled.
1911 !!
1912 !! It also keeps a log of all shell commands executed.
1913 subroutine use_execute_command_line(command,exitstat,cmdstat,cmdmsg)
1916#if ( lwith_intel && !lwith_gnu)
1917 use ifport
1918#endif
1919
1920 ! input / output
1921 character(len=*), intent(in) :: command !< command to execute
1922 integer, intent(inout), optional :: exitstat !< exit status
1923 integer, intent(inout), optional :: cmdstat !< command status
1924 character(len=*), intent(inout), optional :: cmdmsg !< command message
1925
1926 ! local variables
1927 integer :: istat ! status
1928 character(len=max_str_ln) :: full_name ! full name
1929
1930 ! set up the full shell commands file name
1931 full_name = prog_name//'_'//trim(shell_commands_name)//'.sh'
1932
1933 ! write the command to the file
1934 open(shell_commands_i,file=trim(full_name),status='old',&
1935 &position='append',iostat=istat)
1936 if (istat.eq.0) then
1937 write(shell_commands_i,'(A)',iostat=istat) trim(command)
1938 close(shell_commands_i,iostat=istat)
1939 else
1940 call writo('Failed to write to shell commands log file "'&
1941 &//trim(full_name)//'"',warning=.true.)
1942 end if
1943
1944 ! initialize stati
1945 if (present(exitstat)) exitstat = 0
1946 if (present(cmdstat)) cmdstat = 0
1947 if (present(cmdmsg)) cmdmsg = ''
1948
1949 ! execute command line
1950 if (do_execute_command_line) then
1951#if ( lwith_intel && !lwith_gnu)
1952 exitstat = system(command)
1953#else
1954 call execute_command_line(command,exitstat=exitstat,&
1955 &cmdstat=cmdstat,cmdmsg=cmdmsg)
1956#endif
1957 else
1958 call writo('Command added to log file "'//trim(full_name)//'":')
1959 call lvl_ud(1)
1960 call writo(command)
1961 call lvl_ud(-1)
1962 call writo('Run with "--do_execute_command_line" to execute &
1963 &command line instead')
1964 end if
1965 end subroutine
1966end module output_ops
Deallocates XML_str_type.
Definition HDF5_vars.f90:60
Gather parallel variable in serial version on group master.
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Print 2-D output on a file.
Print 3-D output on a file.
Numerical utilities related to files.
integer function, public nextunit(unit)
Search for available new unit.
Variables pertaining to the different grids used.
Definition grid_vars.f90:4
Operations on HDF5 and XDMF variables.
Definition HDF5_ops.f90:27
integer function, public print_hdf5_grid(grid_id, grid_name, grid_type, grid_time, grid_top, grid_geom, grid_atts, grid_grids, reset, ind_plot)
Prints an HDF5 grid.
Definition HDF5_ops.f90:886
subroutine, public print_hdf5_att(att_id, att_dataitem, att_name, att_center, att_type, reset, ind_plot)
Prints an HDF5 attribute.
Definition HDF5_ops.f90:696
integer function, public add_hdf5_item(file_info, xdmf_item, reset, ind_plot)
Add an XDMF item to a HDF5 file.
Definition HDF5_ops.f90:340
integer function, public open_hdf5_file(file_info, file_name, sym_type, descr, ind_plot, cont_plot)
Opens an HDF5 file and accompanying xmf file and returns the handles.
Definition HDF5_ops.f90:78
subroutine, public merge_hdf5_3d_data_items(merged_id, dataitem_ids, var_name, dim_tot, reset, ind_plot, cont_plot)
Joins dataitems in a vector.
Definition HDF5_ops.f90:603
subroutine, public print_hdf5_geom(geom_id, geom_type, geom_dataitems, reset, ind_plot)
Prints an HDF5 geometry.
Definition HDF5_ops.f90:805
integer function, public close_hdf5_file(file_info, ind_plot, cont_plot)
Closes an HDF5 file and writes the accompanying xmf file.
Definition HDF5_ops.f90:264
integer function, public print_hdf5_3d_data_item(dataitem_id, file_info, var_name, var, dim_tot, loc_dim, loc_offset, init_val, ind_plot, cont_plot)
Prints an HDF5 data item.
Definition HDF5_ops.f90:396
subroutine, public print_hdf5_top(top_id, top_type, top_n_elem, ind_plot)
Prints an HDF5 topology.
Definition HDF5_ops.f90:755
Variables pertaining to HDF5 and XDMF.
Definition HDF5_vars.f90:4
Numerical utilities related to giving output.
Definition messages.f90:4
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition messages.f90:254
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition messages.f90:275
Numerical utilities related to MPI.
Numerical variables used by most other modules.
Definition num_vars.f90:4
integer, parameter, public dp
double precision
Definition num_vars.f90:46
character(len=7), public script_dir
directory where to save scripts for plots
Definition num_vars.f90:154
integer, parameter, public shell_commands_i
file number of shell commands file
Definition num_vars.f90:186
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
complex(dp), parameter, public iu
complex unit
Definition num_vars.f90:85
integer, dimension(2), public plot_size
size of plot in inches
Definition num_vars.f90:172
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
character(len=5), public plot_dir
directory where to save plots
Definition num_vars.f90:153
character(len=4), public prog_name
name of program, used for info
Definition num_vars.f90:54
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
Definition num_vars.f90:175
integer, public rank
MPI rank.
Definition num_vars.f90:68
logical, public do_execute_command_line
call "execute_command_line" inside program
Definition num_vars.f90:148
real(dp), public tol_zero
tolerance for zeros
Definition num_vars.f90:133
character(len=14), parameter, public shell_commands_name
name of shell commands file
Definition num_vars.f90:56
character(len=4), public data_dir
directory where to save data for plots
Definition num_vars.f90:155
logical, public no_plots
no plots made
Definition num_vars.f90:140
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
subroutine, public plot_diff_hdf5(a, b, file_name, tot_dim, loc_offset, descr, output_message)
Takes two input vectors and plots these as well as the relative and absolute difference in a HDF5 fil...
subroutine, public draw_ex(var_names, draw_name, nplt, draw_dim, plot_on_screen, ex_plot_style, data_name, draw_ops, extra_ops, is_animated, ranges, delay, persistent)
Use external program to draw a plot.
Operations on strings.
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
character((len(input_strings)+2) *size(input_strings)) function, public merge_strings(input_strings)
Merge array of strings.
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Type for grids.
Definition grid_vars.f90:59
HDF5 data type, containing the information about HDF5 files.
Definition HDF5_vars.f90:40
XML strings used in XDMF.
Definition HDF5_vars.f90:33