PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
output_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
4 !------------------------------------------------------------------------------!
5 module output_ops
6 #include <PB3D_macros.h>
7  use str_utilities
8  use messages
9  use num_vars, only: dp, max_str_ln, no_plots, iu, plot_dir, data_dir, &
10  &script_dir, plot_size
11  use files_utilities, only: nextunit
12  use grid_vars, only: grid_type
13 
14 
15  implicit none
16  private
17  public print_ex_2d, print_ex_3d, draw_ex, plot_hdf5, plot_diff_hdf5
18 
19  ! global variables
20  integer :: temp_id(2) = 1
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;'
27 #if ldebug
28  character(len=0) :: err_output_str = ''
29 #else
30  character(len=14) :: err_output_str = ' 2> /dev/null'
31 #endif
32 
33  ! interfaces
34 
47  interface print_ex_2d
49  module procedure print_ex_2d_ind
51  module procedure print_ex_2d_arr
52  end interface
53 
65  interface print_ex_3d
67  module procedure print_ex_3d_ind
69  module procedure print_ex_3d_arr
70  end interface
71 
72 
137  interface plot_hdf5
139  module procedure plot_hdf5_ind
141  module procedure plot_hdf5_arr
142  end interface
143 
144 contains
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
149  character(len=*), intent(in) :: file_name_i
150  real(dp), intent(in) :: y(1:)
151  real(dp), intent(in), optional :: x(1:)
152  logical, intent(in), optional :: draw
153  logical, intent(in), optional :: persistent
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
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(:)
176  character(len=*), intent(in) :: file_name_i
177  real(dp), intent(in) :: y(1:,1:)
178  real(dp), intent(in), optional :: x(1:,1:)
179  logical, intent(in), optional :: draw
180  logical, intent(in), optional :: persistent
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 
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
271  character(len=*), intent(in) :: file_name_i
272  real(dp), intent(in) :: z(1:,1:)
273  real(dp), intent(in), optional :: x(1:,1:)
274  real(dp), intent(in), optional :: y(1:,1:)
275  logical, intent(in), optional :: draw
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
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(:)
314  character(len=*), intent(in) :: file_name_i
315  real(dp), intent(in) :: z(1:,1:,1:)
316  real(dp), intent(in), optional :: x(1:,1:,1:)
317  real(dp), intent(in), optional :: y(1:,1:,1:)
318  logical, intent(in), optional :: draw
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 
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)
435  use hdf5_ops, only: open_hdf5_file, add_hdf5_item, print_hdf5_top, &
436  &print_hdf5_geom, print_hdf5_3d_data_item, print_hdf5_att, &
437  &print_hdf5_grid, close_hdf5_file, merge_hdf5_3d_data_items
438  use hdf5_vars, only: dealloc_xml_str, &
439  &xml_str_type, hdf5_file_type
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(:)
445  character(len=*), intent(in) :: file_name
446  real(dp), intent(in), target :: vars(:,:,:,:)
447  integer, intent(in), optional :: tot_dim(4)
448  integer, intent(in), optional :: loc_offset(4)
449  real(dp), intent(in), target, optional :: X(:,:,:,:)
450  real(dp), intent(in), target, optional :: Y(:,:,:,:)
451  real(dp), intent(in), target, optional :: Z(:,:,:,:)
452  integer, intent(in), optional :: col_id
453  integer, intent(in), optional :: col
454  integer, intent(in), optional :: sym_type
455  logical, intent(in), optional :: cont_plot
456  character(len=*), intent(in), optional :: descr
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
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
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
952  character(len=*), intent(in) :: file_name
953  real(dp), intent(in) :: var(:,:,:)
954  integer, intent(in), optional :: tot_dim(3)
955  integer, intent(in), optional :: loc_offset(3)
956  real(dp), intent(in), target, optional :: X(:,:,:)
957  real(dp), intent(in), target, optional :: Y(:,:,:)
958  real(dp), intent(in), target, optional :: Z(:,:,:)
959  logical, intent(in), optional :: cont_plot
960  character(len=*), intent(in), optional :: descr
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 
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(:)
1084  character(len=*), intent(in) :: draw_name
1085  integer, intent(in) :: nplt
1086  integer, intent(in) :: draw_dim
1087  logical, intent(in) :: plot_on_screen
1088  integer, intent(in), optional :: ex_plot_style
1089  character(len=*), intent(in), optional :: data_name
1090  character(len=*), intent(in), optional :: draw_ops(:)
1091  character(len=*), intent(in), optional :: extra_ops
1092  logical, intent(in), optional :: is_animated
1093  real(dp), intent(in), optional :: ranges(:,:)
1094  integer, intent(in), optional :: delay
1095  logical, intent(in), optional :: persistent
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
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 
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) 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 .eq..and..eq. if (size(ranges,1)2 size(ranges,2)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 .eq. if (n_draw_ops0) 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 .eq. if (istat0) then
1685  ! set up loc_data
1686  allocate(loc_data(2*nplt))
1687 
1688  ! read the data file
1689  istat = 0
1690 .eq. do while (istat0)
1691  read(data_i,*,IOSTAT=istat) loc_data_char ! read first character of data
1692 .eq. if (istat0) then ! read succesful
1693 .ne. if (loc_data_char'#') 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 .gt. if (n_draw_ops0) 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 .gt. if (n_draw_ops0) 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 .ne..or..ne..or. if (size(A,1)size(B,1) size(A,2)size(B,2) &
1797 .ne. &size(A,3)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 .eq..and..eq. if (tot_dim_loc(1)size(A,1) tot_dim_loc(2)size(A,2) &
1809 .and..eq. & tot_dim_loc(3)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 .not. if (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)
1914  use num_vars, only: do_execute_command_line, prog_name, &
1915  &shell_commands_name, shell_commands_i
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 .eq. if (istat0) 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
1966 end module output_ops
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition: MPI_utilities.f90:55
num_vars::ex_plot_style
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
Definition: num_vars.f90:175
output_ops::plot_diff_hdf5
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...
Definition: output_ops.f90:1765
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
num_vars::script_dir
character(len=7), public script_dir
directory where to save scripts for plots
Definition: num_vars.f90:154
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
hdf5_ops::print_hdf5_geom
subroutine, public print_hdf5_geom(geom_id, geom_type, geom_dataitems, reset, ind_plot)
Prints an HDF5 geometry.
Definition: HDF5_ops.f90:805
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
hdf5_ops::print_hdf5_3d_data_item
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
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
hdf5_ops::print_hdf5_top
subroutine, public print_hdf5_top(top_id, top_type, top_n_elem, ind_plot)
Prints an HDF5 topology.
Definition: HDF5_ops.f90:755
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition: HDF5_ops.f90:27
num_vars::plot_size
integer, dimension(2), public plot_size
size of plot in inches
Definition: num_vars.f90:172
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
hdf5_vars
Variables pertaining to HDF5 and XDMF.
Definition: HDF5_vars.f90:4
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
hdf5_ops::add_hdf5_item
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
num_vars::data_dir
character(len=4), public data_dir
directory where to save data for plots
Definition: num_vars.f90:155
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
hdf5_vars::dealloc_xml_str
Deallocates XML_str_type.
Definition: HDF5_vars.f90:60
hdf5_ops::open_hdf5_file
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
output_ops::draw_ex
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.
Definition: output_ops.f90:1079
hdf5_ops::print_hdf5_att
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
output_ops::plot_hdf5
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Definition: output_ops.f90:137
files_utilities::nextunit
integer function, public nextunit(unit)
Search for available new unit.
Definition: files_utilities.f90:27
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
hdf5_ops::merge_hdf5_3d_data_items
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
output_ops::use_execute_command_line
subroutine use_execute_command_line(command, exitstat, cmdstat, cmdmsg)
Executes command line, or displays a message if disabled.
Definition: output_ops.f90:1914
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
hdf5_ops::close_hdf5_file
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
str_utilities::merge_strings
character((len(input_strings)+2) *size(input_strings)) function, public merge_strings(input_strings)
Merge array of strings.
Definition: str_utilities.f90:152
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
output_ops::print_ex_3d
Print 3-D output on a file.
Definition: output_ops.f90:65
num_vars::no_plots
logical, public no_plots
no plots made
Definition: num_vars.f90:140
hdf5_vars::hdf5_file_type
HDF5 data type, containing the information about HDF5 files.
Definition: HDF5_vars.f90:40
num_vars::plot_dir
character(len=5), public plot_dir
directory where to save plots
Definition: num_vars.f90:153
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
hdf5_ops::print_hdf5_grid
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
files_utilities
Numerical utilities related to files.
Definition: files_utilities.f90:4
hdf5_vars::xml_str_type
XML strings used in XDMF.
Definition: HDF5_vars.f90:33
this documentation was created with Doxygen 1.8.17 by Toon Weyens