PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
test.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module test
5 #include <PB3D_macros.h>
6 #include <wrappers.h>
7  use str_utilities
8  use output_ops
9  use messages
10  use num_vars, only: pi, dp, max_str_ln, iu
11  use input_utilities
12 
13  implicit none
14  private
15 #if ldebug
16  public generic_tests
17 #endif
18 
19 #if ldebug
20 contains
38  integer function generic_tests() result(ierr)
39  use num_vars, only: prog_style
40 
41  character(*), parameter :: rout_name = 'generic_tests'
42 
43  ! initialize ierr
44  ierr = 0
45 
46  call writo('Test smoothed second derivative?')
47  if (get_log(.false.)) then
48  ierr = test_calc_d2_smooth()
49  chckerr('')
50  call pause_prog
51  end if
52 
53  call writo('Test splines?')
54  if (get_log(.false.)) then
55  ierr = test_splines()
56  chckerr('')
57  call pause_prog
58  end if
59 
60  call writo('Test lock system?')
61  if (get_log(.false.)) then
62  ierr = test_lock()
63  chckerr('')
64  call pause_prog
65  end if
66 
67  call writo('Test calculation of RZL?')
68  if (get_log(.false.)) then
69  ierr = test_calc_rzl()
70  chckerr('')
71  call pause_prog
72  end if
73 
74  call writo('Test calculation of toroidal functions?')
75  if (get_log(.false.)) then
76  ierr = test_tor_fun()
77  chckerr('')
78  call pause_prog
79  end if
80 
81  ! select according to program style
82  select case (prog_style)
83  case(1) ! PB3D
84  ! do nothing
85  case(2) ! POST
86  call writo('Test reading of subsets?')
87  if (get_log(.false.)) then
88  ierr = test_read_hdf5_subset()
89  chckerr('')
90  call pause_prog
91  end if
92 
93  call writo('Test calculation of volume integral?')
94  if (get_log(.false.)) then
95  ierr = test_calc_int_vol()
96  chckerr('')
97  call pause_prog
98  end if
99  end select
100  end function generic_tests
101 
105  integer function test_calc_d2_smooth() result(ierr)
106  use num_utilities, only: calc_d2_smooth
107 
108  character(*), parameter :: rout_name = 'test_calc_D2_smooth'
109 
110  ! local variables
111  integer :: kd
112  integer :: nx
113  integer :: fil_n
114  integer :: style
115  integer :: sym_m
116  real(dp), allocatable :: x(:)
117  real(dp), allocatable :: y(:,:)
118  real(dp), allocatable :: d2y(:)
119  logical :: ready
120 
121  ! initialize ierr
122  ierr = 0
123 
124  ! user input
125  call writo('testing Holoborodko''s smoothed second derivatives')
126  call lvl_ud(1)
127 
128  ! set up function
129  call writo('how many points?')
130  nx = get_int(lim_lo=20)
131  allocate(x(nx))
132  allocate(y(nx,0:2))
133  allocate(d2y(nx))
134  do kd = 1,nx
135  x(kd) = (kd-1._dp)/(nx-1._dp) ! 0..1
136  y(kd,0) = sin(2*pi*x(kd))
137  y(kd,1) = 2*pi*cos(2*pi*x(kd))
138  y(kd,2) = -(2*pi)**2*sin(2*pi*x(kd))
139  !!y(kd,0) = x(kd)**3
140  !!y(kd,1) = 3._dp*x(kd)**2
141  !!y(kd,2) = 6._dp*x(kd)
142  end do
143 
144  ready = .false.
145  do while (.not.ready)
146  ! get filter length
147  call writo('filter length?')
148  fil_n = get_int(lim_lo=5,lim_hi=nx)
149  sym_m = (fil_n-1)/2
150 
151  ! get style
152  call writo('style?')
153  call lvl_ud(1)
154  call writo('1: central')
155  call writo('2: left')
156  call lvl_ud(-1)
157  style = get_int(lim_lo=1,lim_hi=2)
158 
159  ! calculate derivative
160  ierr = calc_d2_smooth(fil_n,x,y(:,0),d2y,style)
161  chckerr('')
162 
163  ! plot
164  call print_ex_2d(['y ','D2y'],'',reshape([y(:,2),d2y],[nx,2]),&
165  &x=reshape(x,[nx,1]))
166  select case (style)
167  case (1)
168  call print_ex_2d('y-D2y','',&
169  &y(sym_m+1:nx-sym_m,2)-d2y(sym_m+1:nx-sym_m),&
170  &x=x(sym_m+1:nx-sym_m))
171  case (2)
172  call print_ex_2d('y-D2y','',&
173  &y(2*sym_m+1:nx,2)-d2y(2*sym_m+1:nx),&
174  &x=x(2*sym_m+1:nx))
175  end select
176 
177  call writo('repeat?')
178  ready = .not.get_log(.true.)
179  end do
180 
181  call lvl_ud(-1)
182  end function test_calc_d2_smooth
183 
187  integer function test_splines() result(ierr)
188  use num_utilities, only: spline
189  use num_vars, only: tol_zero
190 
191  character(*), parameter :: rout_name = 'test_splines'
192 
193  ! local variables
194  integer :: kd, id
195  integer :: nx
196  integer :: nx_int
197  integer :: ord
198  integer :: max_deriv
199  integer :: bcs(2) ! boundary conditions
200  real(dp) :: bcs_val(2) ! value of derivative at boundary conditions
201  real(dp) :: extrap_fraction
202  real(dp), allocatable :: x(:)
203  real(dp), allocatable :: x_int(:)
204  real(dp), allocatable :: y(:,:)
205  real(dp), allocatable :: y_int(:)
206  logical :: ready
207  logical :: extrap
208  character(len=5) :: side_str(2)
209 
210  ! initialize ierr
211  ierr = 0
212 
213  ! user input
214  call writo('testing splines')
215  call lvl_ud(1)
216 
217  ! set up function
218  nx = 10
219  allocate(x(nx))
220  allocate(y(nx,0:3))
221  do kd = 1,nx
222  x(kd) = (kd-1._dp)/(nx-1._dp) ! 0..1
223  y(kd,0) = sin(2*pi*x(kd))
224  y(kd,1) = 2*pi*cos(2*pi*x(kd))
225  y(kd,2) = -(2*pi)**2*sin(2*pi*x(kd))
226  y(kd,3) = -(2*pi)**3*cos(2*pi*x(kd))
227  end do
228 
229  ! set up interpolated function
230  nx_int = 400
231  allocate(x_int(nx_int))
232  allocate(y_int(nx_int))
233 
234  ready = .false.
235  do while (.not.ready)
236  ! get order
237  call writo('order?')
238  call lvl_ud(1)
239  call writo('1: linear')
240  call writo('2: akima hermite')
241  call writo('3: cubic')
242  call lvl_ud(-1)
243  ord = get_int(lim_lo=1,lim_hi=3)
244 
245  ! get extrapolation fraction
246  call writo('test extrapolation?')
247  extrap = get_log(.false.)
248  if (extrap) then
249  call writo('Extrapolation fraction?')
250  extrap_fraction = get_real(lim_lo=0._dp,lim_hi=10._dp)
251  else
252  extrap_fraction = -tol_zero
253  end if
254 
255  side_str(1) = 'left'
256  side_str(2) = 'right'
257  if (ord.gt.1) then
258  do kd = 1,2
259  call writo('boundary condition '//trim(side_str(kd))//&
260  &' = ?')
261  call lvl_ud(1)
262  call writo('-1: periodic')
263  call writo(' 0: not-a-knot')
264  call writo(' 1: first derivative prescribed')
265  call writo(' 2: second derivative prescribed')
266  call lvl_ud(-1)
267  bcs(kd) = get_int(lim_lo=-1,lim_hi=2)
268  if (bcs(kd).eq.1 .or. bcs(kd).eq.2) then
269  call writo('value of derivative of order '//&
270  &trim(i2str(bcs(kd))))
271  bcs_val(kd) = get_real()
272  else
273  bcs_val(kd) = 0._dp
274  end if
275  end do
276  else
277  bcs = 0
278  bcs_val = 0._dp
279  end if
280 
281  select case (ord)
282  case (1:2)
283  max_deriv = 1
284  case (3)
285  max_deriv = 3
286  end select
287 
288  ! set up x_int
289  do kd = 1,nx_int
290  x_int(kd) = -extrap_fraction + &
291  &(1+2*extrap_fraction)*(kd-1._dp)/(nx_int-1._dp) ! 0..1 with possible extrapolation fraction
292  end do
293 
294  ! calculate spline and possibly plot
295  do id = 0,max_deriv
296  ierr = spline(x,y(:,0),x_int,y_int,ord=ord,deriv=id,bcs=bcs,&
297  &bcs_val=bcs_val,extrap=extrap) ! use y(0)
298  chckerr('')
299  call plot_y(x,x_int,y(:,id),y_int,id) ! use y(id)
300  end do
301 
302  call writo('repeat?')
303  ready = .not.get_log(.true.)
304  end do
305 
306  call lvl_ud(-1)
307  contains
308  subroutine plot_y(x,x_int,y,y_int,deriv)
309  ! input / output
310  real(dp), intent(in) :: x(:)
311  real(dp), intent(in) :: x_int(:)
312  real(dp), intent(in) :: y(:)
313  real(dp), intent(in) :: y_int(:)
314  integer, intent(in) :: deriv
315 
316  ! local variables
317  integer :: nx, nx_int
318  real(dp) :: minmax_y(2)
319  real(dp), allocatable :: x_plot(:,:)
320  real(dp), allocatable :: y_plot(:,:)
321  character(len=max_str_ln), allocatable :: plot_names(:)
322 
323  ! set up plot variables
324  nx = size(x)
325  nx_int = size(x_int)
326  allocate(x_plot(max(nx,nx_int),2)) ! (x_int,x)
327  allocate(y_plot(max(nx,nx_int),2)) ! (y_int,y)
328 
329  ! x,y
330  x_plot(1:nx_int,1) = x_int
331  x_plot(nx_int+1:max(nx,nx_int),1) = x_int(nx_int)
332  y_plot(1:nx_int,1) = y_int
333  y_plot(nx_int+1:max(nx,nx_int),1) = y_int(nx_int)
334 
335  ! int x,y
336  x_plot(1:nx,2) = x
337  x_plot(nx+1:max(nx,nx_int),2) = x(nx)
338  y_plot(1:nx,2) = y
339  y_plot(nx+1:max(nx,nx_int),2) = y(nx)
340 
341  ! knots
342  minmax_y(1) = minval(y_plot(:,1:2))
343  minmax_y(2) = maxval(y_plot(:,1:2))
344 
345  ! plot
346  allocate(plot_names(2))
347  plot_names(1) = 'D'//trim(i2str(deriv))//'y_int'
348  plot_names(2) = 'D'//trim(i2str(deriv))//'y'
349  call print_ex_2d(plot_names,plot_names(1),y_plot,x=x_plot)
350  end subroutine plot_y
351  end function test_splines
352 
356  integer function test_lock() result(ierr)
357  use mpi_vars, only: lock_type
358  use mpi_utilities, only: lock_req_acc, lock_return_acc, wait_mpi, &
359  &lock_header, lock_wl_change, &
360  &debug_lock
361  use num_vars, only: rank
362 
363  character(*), parameter :: rout_name = 'test_lock'
364 
365  ! local variables
366  type(lock_type) :: lock
367  integer :: wait_time
368  integer(kind=8) :: sleep_time(2)
369  integer :: testcase
370  integer :: n_cycles
371  integer :: id
372  logical :: blocking
373  logical :: red_lat
374  integer, allocatable :: wl(:)
375  character(len=max_str_ln) :: title
376 
377  ! initialize ierr
378  ierr = 0
379 
380  ! debug messages
381  debug_lock = .true.
382 
383  call writo('which test?')
384  call lvl_ud(1)
385  call writo('1. blocking')
386  call writo('2. non-blocking')
387  call writo('3. blocking and non-blocking')
388  testcase = get_int(lim_lo=1,lim_hi=3)
389  call lvl_ud(-1)
390 
391  call writo('How many cycles?')
392  n_cycles = get_int(lim_lo=1)
393 
394  if (testcase.eq.2 .or. testcase.eq.3) then
395  call writo('Do you want to check whether locking and unlocking the &
396  &window multiple times by the master helps reduce latency?')
397  red_lat = get_log(.false.)
398  end if
399 
400  ! set name
401  select case (testcase)
402  case (1)
403  title = 'blocking'
404  case (2)
405  title = 'non-blocking'
406  case (3)
407  title = 'blocking and non-blocking'
408  end select
409 
410  ! create lock
411  ierr = lock%init(100)
412  chckerr('')
413 
414  ! test blocking access
415  call writo('Testing '//trim(title)//' access')
416  call lvl_ud(1)
417  call writo('Every process immediately requests access and then &
418  &waits some seconds between 0 and 3')
419 
420  ! set name
421  select case (testcase)
422  case (1)
423  call writo('Every process is blocking')
424  case (2)
425  call writo('Every process is non-blocking')
426  case (3)
427  call writo('In the first cycle, every process is non-blocking')
428  call writo('Afterwards, each process becomes blocking &
429  &sequentially')
430  end select
431 
432  ! synchronize MPI
433  ierr = wait_mpi()
434  chckerr('')
435 
436  do id = 1,n_cycles
437  select case (testcase)
438  case (1) ! blocking
439  blocking = .true.
440  case (2) ! non-blocking
441  blocking = .false.
442  case (3) ! both
443  if (id.eq.rank+2) then
444  blocking = .true.
445  else
446  blocking = .false.
447  end if
448  end select
449  ierr = lock_req_acc(lock,blocking)
450  chckerr('')
451 
452  wait_time = random_int([0,3])
453  write(*,*) trim(lock_header(lock)),'sleeps ',&
454  &wait_time,' seconds'
455  if (red_lat .and. rank.eq.0 .and. wait_time.gt.0) then
456  call system_clock(sleep_time(1))
457  call system_clock(sleep_time(2))
458  do while (1.e-9_dp*(sleep_time(2)-sleep_time(1)).lt.wait_time)
459  debug_lock = .false.
460  ierr = lock_wl_change(2,lock%blocking,lock,wl)
461  chckerr('')
462  call system_clock(sleep_time(2))
463  end do
464  debug_lock = .true.
465  else
466  call sleep(wait_time)
467  endif
468  write(*,*) trim(lock_header(lock)),'has slept ',&
469  &wait_time,' seconds, returning lock'
470 
471  ierr = lock_return_acc(lock)
472  chckerr('')
473  end do
474 
475  ! clean up
476  ierr = lock%dealloc()
477  chckerr('')
478 
479  ! synchronize MPI
480  ierr = wait_mpi()
481  chckerr('')
482 
483  call lvl_ud(-1)
484  contains
485  ! not very high quality random number between the limits
487  integer function random_int(lim) result(res)
488  ! input / output
489  integer, intent(in) :: lim(2)
490 
491  ! local variables
492  integer :: seed(33)
493  real(dp) :: r_nr
494  integer :: n_steps = 5
495  integer :: id
496 
497  do id = 1,n_steps
498  ! random number for seed
499  call random_number(r_nr)
500 
501  ! set seed
502  seed = int(1000_dp*r_nr*(rank+1))
503  call random_seed(put=seed)
504  end do
505 
506  ! random numbeer
507  call random_number(r_nr)
508  res = lim(1) + floor((lim(2)-lim(1)+1)*r_nr)
509  end function random_int
510  end function test_lock
511 
515  integer function test_calc_rzl() result(ierr)
516  use pb3d_ops, only: reconstruct_pb3d_in
517  use num_vars, only: eq_style
518 
519  character(*), parameter :: rout_name = 'test_calc_RZL'
520 
521  ! initialize ierr
522  ierr = 0
523 
524  ! preliminary things
525  ierr = reconstruct_pb3d_in('in') ! reconstruct miscellaneous PB3D output variables
526  chckerr('')
527 
528  ! select according to equilibrium style
529  select case (eq_style)
530  case(1) ! VMEC
531  ierr = test_calc_rzl_vmec()
532  chckerr('')
533  case(2) ! HELENA
534  ! do nothing
535  end select
536 
537  ! clean up
538  call dealloc_in()
539  contains
541  integer function test_calc_rzl_vmec() result(ierr)
542  use num_vars, only: rank
543  use grid_vars, only: n_r_eq, grid_type
544  use grid_utilities, only: calc_eqd_grid
545  use vmec_utilities, only: calc_trigon_factors, fourier2real
546  use vmec_vars, only: r_v_c, r_v_s, z_v_c, z_v_s, l_v_c, l_v_s, &
547  &is_asym_v
548 
549  character(*), parameter :: rout_name = 'test_calc_RZL_VMEC'
550 
551  ! local variables
552  integer :: id, jd ! counters
553  integer :: dims(3) ! dimensions of grid
554  integer :: norm_id ! normal point for which to do the analysis
555  real(dp) :: grid_lims(2,2) ! limits on grid
556  real(dp), allocatable :: norm(:,:,:) ! copy of normal grid variable
557  real(dp), allocatable :: R(:,:,:), Z(:,:,:), L(:,:,:) ! local R, Z and L
558  type(grid_type) :: grid ! grid
559  integer :: deriv(3) ! local derivative
560  integer :: max_deriv(3) ! maximum derivative
561 
562  ! user output
563  call writo('Going to test the calculation of R, Z and L')
564  call lvl_ud(1)
565 
566  ! initialize ierr
567  ierr = 0
568 
569  if (rank.eq.0) then
570  ! set normal point
571  call writo('For which normal point do you want the analysis?')
572  norm_id = get_int(lim_lo=1,lim_hi=n_r_eq)
573 
574  ! set dimensions [theta,zeta,r]
575  dims = [1,1,1]
576  do id = 2,3
577  call writo('grid size in dimension '//trim(i2str(id))//'?')
578  dims(id-1) = get_int(lim_lo=1)
579 
580  call writo('limits for angle '//trim(i2str(id))//' [pi]?')
581  grid_lims(1,id-1) = get_real()*pi
582  grid_lims(2,id-1) = get_real(lim_lo=grid_lims(1,id-1)/pi)*pi
583  end do
584 
585  ierr = grid%init(dims)
586  chckerr('')
587 
588  ! create grid
589  allocate(norm(dims(1),dims(2),dims(3)))
590  ierr = calc_eqd_grid(grid%theta_E,grid_lims(1,1),&
591  &grid_lims(2,1),1)
592  chckerr('')
593  ierr = calc_eqd_grid(grid%zeta_E,grid_lims(1,2),&
594  &grid_lims(2,2),2)
595  chckerr('')
596  norm = (norm_id-1._dp)/(n_r_eq-1)
597  grid%r_E = norm(1,1,:)
598 
599  ! set maximum derivatives [r,theta,zeta]
600  max_deriv(1) = 0
601  do id = 2,3
602  call writo('maximum derivative in dimension '//&
603  &trim(i2str(id))//'?')
604  max_deriv(id) = get_int(lim_lo=0,lim_hi=3)
605  end do
606 
607  ! calculate R, Z and L
608  allocate(r(dims(1),dims(2),dims(3)))
609  allocate(z(dims(1),dims(2),dims(3)))
610  allocate(l(dims(1),dims(2),dims(3)))
611  ierr = calc_trigon_factors(grid%theta_E,grid%zeta_E,&
612  &grid%trigon_factors)
613  chckerr('')
614  do jd = 0,max_deriv(3)
615  do id = 0,max_deriv(2)
616  deriv = [0,id,jd]
617  call writo('Calculating for deriv = ['//&
618  &trim(i2str(deriv(1)))//','//&
619  &trim(i2str(deriv(2)))//','//&
620  &trim(i2str(deriv(3)))//']')
621  call lvl_ud(1)
622  ierr = fourier2real(r_v_c(:,norm_id:norm_id,deriv(1)),&
623  &r_v_s(:,norm_id:norm_id,deriv(1)),&
624  !&grid%theta_E,grid%zeta_E,&
625  &grid%trigon_factors,&
626  &r,sym=[.true.,is_asym_v],&
627  &deriv=[deriv(2),deriv(3)])
628  chckerr('')
629  ierr = fourier2real(z_v_c(:,norm_id:norm_id,deriv(1)),&
630  &z_v_s(:,norm_id:norm_id,deriv(1)),&
631  !&grid%theta_E,grid%zeta_E,&
632  &grid%trigon_factors,&
633  &z,sym=[is_asym_v,.true.],&
634  &deriv=[deriv(2),deriv(3)])
635  chckerr('')
636  ierr = fourier2real(l_v_c(:,norm_id:norm_id,deriv(1)),&
637  &l_v_s(:,norm_id:norm_id,deriv(1)),&
638  !&grid%theta_E,grid%zeta_E,&
639  &grid%trigon_factors,&
640  &l,sym=[is_asym_v,.true.],&
641  &deriv=[deriv(2),deriv(3)])
642  chckerr('')
643  call plot_hdf5(['R','Z','L'],'TEST_RZL_0_'//&
644  &trim(i2str(id))//'_'//trim(i2str(jd)),&
645  &reshape([r,z,l],[dims,3]),&
646  &x=reshape([grid%theta_E],[dims,1]),&
647  &y=reshape([grid%zeta_E],[dims,1]),&
648  &z=reshape([norm],[dims,1]))
649  call lvl_ud(-1)
650  end do
651  end do
652  end if
653 
654  ! user output
655  call lvl_ud(-1)
656  call writo('Test complete')
657  end function test_calc_rzl_vmec
658  end function test_calc_rzl
659 
663  integer function test_tor_fun() result(ierr)
664  use num_vars, only: rank
665  use grid_utilities, only: calc_eqd_grid
666  use dtorh, only: dtorh1
667 
668  character(*), parameter :: rout_name = 'test_tor_fun'
669 
670  ! local variables
671  integer :: jd, kd ! counters
672  integer :: npt ! number of points
673  real(dp) :: lim(2) ! limits
674  real(dp), allocatable :: x(:) ! abscissa
675  real(dp), allocatable :: f(:,:,:) ! ordinate P and Q
676  integer :: n ! degree
677  integer :: m ! order
678  integer :: max_n ! max. degree that can be calculated
679  integer :: min_n ! max. degree that can be plot
680  logical :: done ! whether done
681  logical :: test_approx ! test approximation
682  logical :: test_rel_diff ! relative difference, not absolute
683  character(len=max_str_ln), allocatable :: plot_name(:) ! names of plot
684  character(len=max_str_ln), allocatable :: fun_names(:,:) ! names of functions that are plot
685 
686  ! initialize ierr
687  ierr = 0
688 
689  done = .false.
690  test_approx = .false.
691 
692  do while (.not.done .and. rank.eq.0)
693  call lvl_ud(1)
694 
695  ! user input
696  call writo('Test approximation at singular point?')
697  test_approx = get_log(.false.)
698  test_rel_diff = .false.
699  if (test_approx) then
700  call writo('Relative difference in stead of absolute?')
701  test_rel_diff = get_log(.false.)
702  end if
703 
704  ! user output
705  call writo('Lower bound to plot?')
706  lim(1) = get_real(lim_lo=1._dp)
707  call writo('Upper bound to plot?')
708  lim(2) = get_real(lim_lo=lim(1))
709  call writo('How many points?')
710  npt = get_int(lim_lo=1)
711  call writo('Max. degree? (subscript n in P_(n-0.5)^m)')
712  n = get_int(lim_lo=0)
713  call writo('Min. degree? (subscript n in P_(n-0.5)^m)')
714  min_n = get_int(lim_lo=0,lim_hi=n)
715  if (.not.test_approx) then
716  call writo('Order? (superscript m in P_(n-0.5)^m)')
717  m = get_int(lim_lo=0)
718  else
719  call writo('Order is 0 when testing approximation')
720  m = 0
721  end if
722 
723  ! set up
724  allocate(x(npt))
725  allocate(f(npt,0:n,2))
726  ierr = calc_eqd_grid(x,lim(1),lim(2))
727  chckerr('')
728  if (test_approx) then
729  allocate(plot_name(min_n:n))
730  allocate(fun_names(min_n:n,2))
731  if (test_rel_diff) then
732  do kd = min_n,n
733  plot_name(kd) = 'tor_Q_rel_diff_'//trim(i2str(kd))
734  fun_names(kd,1) = 'ReldiffQ_{'//trim(i2str(kd))//&
735  &'-0.5}^'//trim(i2str(m))
736  fun_names(kd,2) = 'Q_{'//trim(i2str(kd))//&
737  &'-0.5}^'//trim(i2str(m))
738  end do
739  else
740  do kd = min_n,n
741  plot_name(kd) = 'tor_Q_comp_'//trim(i2str(kd))
742  fun_names(kd,1) = 'ApproxQ_{'//trim(i2str(kd))//&
743  &'-0.5}^'//trim(i2str(m))
744  fun_names(kd,2) = 'Q_{'//trim(i2str(kd))//&
745  &'-0.5}^'//trim(i2str(m))
746  end do
747  end if
748  else
749  allocate(plot_name(2))
750  allocate(fun_names(2,min_n:n))
751  plot_name(1) = 'tor_P'
752  plot_name(2) = 'tor_Q'
753  do kd = min_n,n
754  fun_names(1,kd) = 'P_{'//trim(i2str(kd))//'-0.5}^'//&
755  &trim(i2str(m))
756  fun_names(2,kd) = 'Q_{'//trim(i2str(kd))//'-0.5}^'//&
757  &trim(i2str(m))
758  end do
759  end if
760 
761  ! calculate
762  do kd = 1,npt
763  ierr = dtorh1(x(kd),m,n,f(kd,:,1),f(kd,:,2),max_n)
764  chckerr('')
765  if (max_n.lt.n) call writo('For point x = '//&
766  &trim(r2str(x(kd)))//', the solution did not converge for &
767  &degree > '//trim(i2str(max_n)),warning=.true.)
768  if (test_approx) then
769  f(kd,:,1) = -0.5*log((x(kd)-1._dp)/32)
770  do jd = 1,n
771  f(kd,jd:n,1) = f(kd,jd:n,1) - 2_dp/(2._dp*jd-1._dp)
772  end do
773  if (test_rel_diff) f(kd,:,1) = &
774  &2._dp*(f(kd,:,1)-f(kd,:,2))/(f(kd,:,1)+f(kd,:,2))
775  end if
776  end do
777 
778  ! plot
779  do kd = lbound(plot_name,1),ubound(plot_name,1)
780  if (test_approx) then
781  call print_ex_2d(fun_names(kd,:),trim(plot_name(kd)),&
782  &f(:,kd,:),x=reshape(x,[npt,1]))
783  else
784  call print_ex_2d(fun_names(kd,:),trim(plot_name(kd)),&
785  &f(:,min_n:n,kd),x=reshape(x,[npt,1]))
786  end if
787  call draw_ex(fun_names(kd,:),trim(plot_name(kd)),&
788  &size(fun_names,2),1,.false.)
789  end do
790 
791  ! clean up
792  deallocate(x,f,fun_names,plot_name)
793 
794  ! user input
795  call writo('Test again?')
796  done = .not.get_log(.true.)
797 
798  call lvl_ud(-1)
799  end do
800  end function test_tor_fun
801 
805  integer function test_read_hdf5_subset() result(ierr)
806  use num_vars, only: rank, eq_style, n_procs
807  use mpi_utilities, only: wait_mpi
808  use grid_vars, only: grid_type
809  use grid_utilities, only: copy_grid, calc_xyz_grid, trim_grid
810  use rich_vars, only: rich_lvl
811  use input_utilities, only: dealloc_in
812  use pb3d_ops, only: reconstruct_pb3d_in, reconstruct_pb3d_grid, &
813  &reconstruct_pb3d_eq_1, reconstruct_pb3d_eq_2, &
814  &reconstruct_pb3d_x_1, reconstruct_pb3d_sol
815  use eq_vars, only: eq_1_type, eq_2_type
816  use x_vars, only: x_1_type, modes_type
817  use sol_vars, only: sol_type
818  use x_ops, only: init_modes, setup_modes
819  use vmec_utilities, only: calc_trigon_factors
820 
821  character(*), parameter :: rout_name = 'test_read_HDF5_subset'
822 
823  ! local variables
824  type(modes_type) :: mds ! general modes variables
825  type(grid_type) :: grid ! equilibrium or perturbation grid
826  type(grid_type) :: grid_sub ! grid for subset
827  type(grid_type) :: grid_trim ! trimmed grid for subset
828  type(grid_type) :: grid_eq ! equilibrium grid for XYZ calculation
829  type(eq_1_type) :: eq_1 ! flux equilibrium variables
830  type(eq_2_type) :: eq_2 ! metric equilibrium variables
831  type(x_1_type) :: x_1 ! vectorial perturbation variables
832  type(sol_type) :: sol ! solution variables
833  integer :: rich_lvl_name ! either the Richardson level or zero, to append to names
834  integer :: rich_lvl_name_eq ! either the Richardson level or zero, to append to names
835  integer :: id, jd, kd ! counters
836  integer :: i_sub ! counter
837  integer :: n_tot(3) ! n of equilibrium grid
838  integer :: n_sub(3) ! nr. of subsets
839  integer :: n_sub_norm ! nr. of subsets
840  integer :: norm_ol = 4 ! normal overlap
841  integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
842  integer :: i_lim(2) ! normal limit for this process
843  integer :: plot_dim(3) ! dimensions of plot
844  integer :: plot_offset(3) ! local offset of plot
845  integer :: test_case ! which case to test (1: eq_2, 2: X_1)
846  integer :: i_sec ! which mode for X_1 or sol
847  integer, allocatable :: n_sub_loc(:) ! number of subsets per process
848  integer, allocatable :: n_sub_norm_loc(:) ! number of normal poins per process per subset
849  integer, allocatable :: lims(:,:,:) ! limits (3,2,product(n_sub))
850  integer, allocatable :: lims_loc(:,:,:) ! limits for current process
851  real(dp), allocatable :: xyz(:,:,:,:) ! X, Y and Z on output grid
852  real(dp), allocatable :: xyz_i(:,:,:,:) ! X, Y and Z for slab plot
853  character(len=max_str_ln) :: description ! description for plot
854  character(len=8) :: name_suff(2) ! suffix of name, without and possibly with rank appended
855  character(len=3) :: grid_name ! either eq or X
856  character :: xyz_name(3) ! X, Y and Z
857  logical :: direct_grid_read ! whether grid is directly read with subsets or copied
858  logical :: divide_over_procs ! whether a division over process should be made
859  logical :: pers_out ! whether output is persistent (all procs) or not (only master)
860  logical :: tot_rich ! whether variable is spread out over multiple Richardson levels or not
861 
862  ! initialize ierr
863  ierr = 0
864 
865  ! user output
866  call writo('Checking reading of HDF5 subsets')
867  call lvl_ud(1)
868 
869  ! some preliminary things
870  ierr = reconstruct_pb3d_in('in') ! reconstruct miscellaneous PB3D output variables
871  chckerr('')
872  xyz_name(1) = 'X'
873  xyz_name(2) = 'Y'
874  xyz_name(3) = 'Z'
875 
876  ! get information about which variables to reconstruct from user
877  call writo('Which variable do you want to test?')
878  call lvl_ud(1)
879  call writo('1. eq_2: Jac_FD and kappa_n')
880  call writo('2. X_1: IM(U_0), RE(DU_1)')
881  call writo('3. sol: RE(sol_vec)')
882  call lvl_ud(-1)
883  test_case = get_int(lim_lo=1,lim_hi=3)
884 
885  ! set up whether Richardson level has to be appended to the name and
886  ! whether tot_rich is used or not
887  select case (eq_style)
888  case (1) ! VMEC
889  select case (test_case)
890  case (1) ! eq_2
891  rich_lvl_name = rich_lvl ! append richardson level
892  tot_rich = .true.
893  case (2) ! X_1
894  rich_lvl_name = rich_lvl ! append richardson level
895  tot_rich = .true.
896  case (3) ! sol
897  rich_lvl_name = 0 ! do not append name
898  tot_rich = .false.
899  end select
900  rich_lvl_name_eq = rich_lvl ! append richardson level
901  case (2) ! HELENA
902  rich_lvl_name = 0 ! do not append name
903  rich_lvl_name_eq = 0 ! do not append name
904  tot_rich = .false.
905  end select
906  select case (test_case)
907  case (1) ! eq_2
908  grid_name = 'eq'
909  case (2) ! X_1
910  grid_name = 'X'
911  case (3) ! sol
912  grid_name = 'sol'
913  end select
914 
915  ! reconstruct full grid
916  ierr = reconstruct_pb3d_grid(grid,trim(grid_name),&
917  &rich_lvl=rich_lvl_name,tot_rich=.true.)
918  chckerr('')
919  n_tot = grid%n
920 
921  ! get direct read flag from user
922  call writo('Read the grids directly using subsets?')
923  direct_grid_read = get_log(.true.)
924 
925  ! get flag concerning division over MPI procs
926  if (n_procs.gt.1) then
927  call writo('Divide over MPI processes?')
928  divide_over_procs = get_log(.true.)
929  else
930  divide_over_procs = .false.
931  end if
932 
933  ! get limits from user
934  call writo('Total '//trim(grid_name)//' grid size: ['//&
935  &trim(i2str(n_tot(1)))//','//trim(i2str(n_tot(2)))//','//&
936  &trim(i2str(n_tot(3)))//']')
937  do id = 1,3
938  if (n_tot(id).gt.0) then
939  call writo('how many divisions in dimension '//trim(i2str(id)))
940  n_sub(id) = get_int(lim_lo=0,lim_hi=n_tot(id)-1) + 1
941  else
942  n_sub(id) = 1
943  end if
944  end do
945  allocate(lims(3,2,product(n_sub)))
946 
947  ! also construct equilibrium grid for XYZ calculation, using lims
948  lims(1,:,1) = [0,0]
949  lims(2,:,1) = [0,0]
950  lims(3,:,1) = [-1,-1]
951  ierr = reconstruct_pb3d_grid(grid_eq,'eq',rich_lvl=rich_lvl_name_eq,&
952  &tot_rich=.true.,lim_pos=lims(:,:,1))
953  chckerr('')
954  call print_ex_2d('r_E','r_E_'//trim(i2str(rank)),grid_eq%r_E,&
955  &draw=.false.)
956  call draw_ex(['r_E'],'r_E_'//trim(i2str(rank)),1,1,.false.)
957 
958  ! set limits
959  i_sub = 1
960  do kd = 1,n_sub(3)
961  do jd = 1,n_sub(2)
962  do id = 1,n_sub(1)
963  lims(1,1,i_sub) = (id-1)*(n_tot(1)/n_sub(1)) + &
964  &min(mod(n_tot(1),n_sub(1)),id-1) + 1
965  lims(1,2,i_sub) = id*(n_tot(1)/n_sub(1)) + &
966  &min(mod(n_tot(1),n_sub(1)),id)
967  lims(2,1,i_sub) = (jd-1)*(n_tot(2)/n_sub(2)) + &
968  &min(mod(n_tot(2),n_sub(2)),jd-1) + 1
969  lims(2,2,i_sub) = jd*(n_tot(2)/n_sub(2)) + &
970  &min(mod(n_tot(2),n_sub(2)),jd)
971  lims(3,1,i_sub) = (kd-1)*(n_tot(3)/n_sub(3)) + &
972  &min(mod(n_tot(3),n_sub(3)),kd-1) + 1
973  lims(3,2,i_sub) = kd*(n_tot(3)/n_sub(3)) + &
974  &min(mod(n_tot(3),n_sub(3)),kd)
975  i_sub = i_sub+1
976  end do
977  end do
978  end do
979 
980  ! divide under processes: either just dividing over all processes with
981  ! full normal subranges, or by letting every process work on the same
982  ! subranges, dividing these normally.
983  allocate(n_sub_loc(n_procs))
984  allocate(n_sub_norm_loc(n_procs))
985  if (divide_over_procs) then
986  pers_out = rank.eq.0
987  n_sub_loc = product(n_sub)
988  allocate(lims_loc(3,2,product(n_sub)))
989  lims_loc = lims
990  call writo('All process perform every job together, dividing each &
991  &one normally with overlap '//trim(i2str(norm_ol)))
992  else
993  pers_out = .true.
994  n_sub_loc = product(n_sub)/n_procs
995  n_sub_loc(1:mod(product(n_sub),n_procs)) = &
996  &n_sub_loc(1:mod(product(n_sub),n_procs)) + 1
997  allocate(lims_loc(3,2,n_sub_loc(rank+1)))
998  lims_loc = &
999  &lims(:,:,sum(n_sub_loc(1:rank))+1:sum(n_sub_loc(1:rank+1)))
1000  call writo('Process '//trim(i2str(rank))//' has subsets:',&
1001  &persistent=.true.)
1002  end if
1003  call lvl_ud(1)
1004  do i_sub = 1,size(lims_loc,3)
1005  call writo('subset '//trim(i2str(i_sub)),persistent=pers_out)
1006  call lvl_ud(1)
1007  do id = 1,3
1008  call writo('lim('//trim(i2str(id))//'&
1009  &) = ['//trim(i2str(lims_loc(id,1,i_sub)))//'..'//&
1010  &trim(i2str(lims_loc(id,2,i_sub)))//']',&
1011  &persistent=pers_out)
1012  end do
1013  call lvl_ud(-1)
1014  end do
1015  call lvl_ud(-1)
1016 
1017  ! synchronize MPI
1018  ierr = wait_mpi()
1019  chckerr('')
1020 
1021  ! for X_1 and sol, set up nm and get information about which mode to
1022  ! plot
1023  if (test_case.eq.2 .or. test_case.eq.3) then
1024  ! set up nm in full grids
1025  ierr = reconstruct_pb3d_eq_1(grid_eq,eq_1,'eq_1')
1026  chckerr('')
1027  ierr = init_modes(grid_eq,eq_1)
1028  chckerr('')
1029  ierr = setup_modes(mds,grid_eq,grid)
1030  chckerr('')
1031  call eq_1%dealloc()
1032  ! get user input
1033  call writo('Which perturbation mode do you want to plot?')
1034  i_sec = get_int(lim_lo=1,lim_hi=size(mds%m,2))
1035  end if
1036 
1037  ! synchronize MPI
1038  ierr = wait_mpi()
1039  chckerr('')
1040 
1041  ! loop over all josb
1042  do i_sub = 1,size(lims_loc,3)
1043  ! user output
1044  call writo('Process '//trim(i2str(rank))//', subset '//&
1045  &trim(i2str(i_sub))//':',persistent=.true.)
1046  call lvl_ud(1)
1047 
1048  ! initialize description
1049  description = ''
1050  do id = 1,3
1051  description = trim(description)//' lim('//trim(i2str(id))//&
1052  &') = ['//trim(i2str(lims_loc(id,1,i_sub)))//'..'//&
1053  &trim(i2str(lims_loc(id,2,i_sub)))//']'
1054  if (id.lt.3) description = trim(description)// ','
1055  end do
1056 
1057  ! set up i_lim
1058  n_sub_norm = lims_loc(3,2,i_sub)-lims_loc(3,1,i_sub)+1
1059  if (divide_over_procs) then
1060  n_sub_norm_loc = n_sub_norm/n_procs
1061  n_sub_norm_loc(1:mod(n_sub_norm,n_procs)) = &
1062  &n_sub_norm_loc(1:mod(n_sub_norm,n_procs)) + 1
1063  i_lim = [sum(n_sub_norm_loc(1:rank))+1,&
1064  &sum(n_sub_norm_loc(1:rank+1))]
1065  if (rank+1.lt.n_procs) i_lim(2) = i_lim(2)+norm_ol
1066  description = trim(description)//', normal range = '//&
1067  &trim(i2str(i_lim(1)))//'..'//trim(i2str(i_lim(2)))
1068  else
1069  n_sub_norm_loc = n_sub_norm
1070  i_lim = lims_loc(3,:,i_sub) - lims_loc(3,1,i_sub) + 1
1071  end if
1072 
1073  ! output description
1074  description = adjustl(description)
1075  call writo(trim(description),persistent=.true.)
1076 
1077  ! set up subset grid
1078  if (direct_grid_read) then
1079  ! check whether it coincides with a direct read using subsets
1080  ierr = reconstruct_pb3d_grid(grid_sub,trim(grid_name),&
1081  &rich_lvl=rich_lvl_name,tot_rich=.true.,&
1082  &lim_pos=lims_loc(:,:,i_sub),grid_limits=i_lim)
1083  chckerr('')
1084  else
1085  ierr = copy_grid(grid,grid_sub,&
1086  &lims_b=lims_loc(:,:,i_sub),i_lim=i_lim)
1087  chckerr('')
1088  end if
1089 
1090  ! trim grid
1091  ierr = trim_grid(grid_sub,grid_trim,norm_id)
1092  chckerr('')
1093 
1094  ! set up plot dimensions and local dimensions
1095  plot_dim = grid_trim%n
1096  plot_offset = [0,0,grid_trim%i_min-1]
1097 
1098  ! reconstruct metric equilibrium or vectorial perturbation variables
1099  ! of subset
1100  select case (test_case)
1101  case (1) ! eq_2
1102  ierr = reconstruct_pb3d_eq_2(grid_sub,eq_2,'eq_2',&
1103  &rich_lvl=rich_lvl_name,tot_rich=.true.,&
1104  &lim_pos=lims_loc(:,:,i_sub))
1105  chckerr('')
1106  case (2) ! X_1
1107  ierr = reconstruct_pb3d_x_1(mds,grid_sub,x_1,'X_1',&
1108  &rich_lvl=rich_lvl_name,tot_rich=.true.,&
1109  &lim_pos=lims_loc(:,:,i_sub),&
1110  &lim_sec_x=[i_sec,i_sec])
1111  chckerr('')
1112  case (3) ! sol
1113  ierr = reconstruct_pb3d_sol(mds,grid_sub,sol,'sol',&
1114  &rich_lvl=rich_lvl,&
1115  &lim_pos=lims_loc(3:3,:,i_sub),&
1116  &lim_sec_sol=[i_sec,i_sec])
1117  chckerr('')
1118  end select
1119 
1120  ! set up name suffix
1121  if (divide_over_procs) then
1122  name_suff(1) = '_'//trim(i2str(i_sub))
1123  else
1124  name_suff(1) = '_'//trim(i2str(sum(n_sub_loc(1:rank))+i_sub))
1125  end if
1126  if (divide_over_procs) then
1127  name_suff(2) = trim(name_suff(1))//'_'//trim(i2str(rank))
1128  else
1129  name_suff(2) = trim(name_suff(1))
1130  end if
1131 
1132  ! set up XYZ and XYZ_i
1133  if (test_case.eq.1 .or. test_case.eq.2) then
1134  if (eq_style.eq.1) then ! if VMEC, calculate trigonometric factors of output grid
1135  ierr = calc_trigon_factors(grid_trim%theta_E,&
1136  &grid_trim%zeta_E,grid_trim%trigon_factors)
1137  chckerr('')
1138  end if
1139  allocate(xyz(grid_trim%n(1),grid_trim%n(2),&
1140  &grid_trim%loc_n_r,3))
1141  ierr = calc_xyz_grid(grid_eq,grid_trim,xyz(:,:,:,1),&
1142  &xyz(:,:,:,2),xyz(:,:,:,3))
1143  chckerr('')
1144  allocate(&
1145  &xyz_i(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,3))
1146  do id = 1,grid_trim%n(1)
1147  xyz_i(id,:,:,1) = lims_loc(1,1,i_sub) + id - 2
1148  end do
1149  do jd = 1,grid_trim%n(2)
1150  xyz_i(:,jd,:,2) = lims_loc(2,1,i_sub) + jd - 2
1151  end do
1152  do kd = 1,grid_trim%loc_n_r
1153  xyz_i(:,:,kd,3) = lims_loc(3,1,i_sub) + kd - 2 + &
1154  &grid_trim%i_min - 1
1155  end do
1156  else
1157  allocate(xyz_i(1,grid_trim%loc_n_r,1,1))
1158  do kd = 1,grid_trim%loc_n_r
1159  xyz_i(1,kd,1,1) = lims_loc(3,1,i_sub) + kd - 2 + &
1160  &grid_trim%i_min - 1
1161  end do
1162  end if
1163 
1164  ! plot
1165  select case (test_case)
1166  case (1) ! eq_2
1167  call plot_hdf5('jac_FD',&
1168  &'jac_FD'//trim(name_suff(1)),&
1169  &eq_2%jac_FD(:,:,norm_id(1):norm_id(2),0,0,0),&
1170  &tot_dim=plot_dim,loc_offset=plot_offset,&
1171  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1172  &descr=description)
1173  call plot_hdf5('kappa_n',&
1174  &'kappa_n'//trim(name_suff(1)),&
1175  &eq_2%kappa_n(:,:,norm_id(1):norm_id(2)),&
1176  &tot_dim=plot_dim,loc_offset=plot_offset,&
1177  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1178  &descr=description)
1179  case (2) ! X_1
1180  call plot_hdf5('U_0_IM',&
1181  &'U_0_IM'//trim(name_suff(1)),&
1182  &ip(x_1%U_0(:,:,norm_id(1):norm_id(2),1)),&
1183  &tot_dim=plot_dim,loc_offset=plot_offset,&
1184  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1185  &descr=description)
1186  call plot_hdf5('DU_1_RE',&
1187  &'DU_1_RE'//trim(name_suff(1)),&
1188  &rp(x_1%DU_1(:,:,norm_id(1):norm_id(2),1)),&
1189  &tot_dim=plot_dim,loc_offset=plot_offset,&
1190  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1191  &descr=description)
1192  case (3) ! sol
1193  call print_ex_2d(['sol_RE'],'sol_RE'//trim(name_suff(2)),&
1194  &rp(sol%vec(1,norm_id(1):norm_id(2),:)),&
1195  &x=xyz_i(1,:,:,1),draw=.false.)
1196  call draw_ex(['sol_RE'],'sol_RE'//trim(name_suff(2)),&
1197  &size(sol%vec,3),1,.false.)
1198  end select
1199  if (test_case.eq.1 .or. test_case.eq.2) then
1200  do id = 1,3
1201  call plot_hdf5(xyz_name(id),xyz_name(id)//&
1202  &trim(name_suff(2)),xyz(:,:,:,id),&
1203  &x=xyz_i(:,:,:,1),y=xyz_i(:,:,:,2),z=xyz_i(:,:,:,3))
1204  end do
1205  end if
1206 
1207  ! clean up
1208  call grid_sub%dealloc()
1209  call grid_trim%dealloc()
1210  select case (test_case)
1211  case (1) ! eq_2
1212  call eq_2%dealloc()
1213  case (2) ! X_1
1214  call x_1%dealloc()
1215  case (3) ! sol
1216  call sol%dealloc()
1217  end select
1218  if (test_case.eq.1 .or. test_case.eq.2) deallocate(xyz)
1219  deallocate(xyz_i)
1220  call lvl_ud(-1)
1221  end do
1222 
1223  ! synchronize MPI
1224  ierr = wait_mpi()
1225  chckerr('')
1226 
1227  ! Now for entire region if wanted
1228  call writo('Do you want to plot the entire region at once?')
1229  if (get_log(.false.) .and. rank.eq.0) then
1230  ! reconstruct metric equilibrium or vectorial perturbation variables
1231  select case (test_case)
1232  case (1) ! eq_2
1233  ierr = reconstruct_pb3d_eq_2(grid,eq_2,'eq_2',&
1234  &rich_lvl=rich_lvl_name,tot_rich=.true.)
1235  chckerr('')
1236  case (2) ! X_1
1237  ierr = reconstruct_pb3d_x_1(mds,grid,x_1,'X_1',&
1238  &rich_lvl=rich_lvl_name,tot_rich=.true.,&
1239  &lim_sec_x=[i_sec,i_sec])
1240  chckerr('')
1241  case (3) ! sol
1242  ierr = reconstruct_pb3d_sol(mds,grid,sol,'sol',&
1243  &rich_lvl=rich_lvl,&
1244  &lim_sec_sol=[i_sec,i_sec])
1245  chckerr('')
1246  end select
1247 
1248  ! set up XYZ and XYZ_i
1249  if (test_case.eq.1 .or. test_case.eq.2) then
1250  if (eq_style.eq.1) then ! if VMEC, calculate trigonometric factors of output grid
1251  ierr = calc_trigon_factors(grid%theta_E,&
1252  &grid%zeta_E,grid%trigon_factors)
1253  chckerr('')
1254  end if
1255  allocate(xyz(grid%n(1),grid%n(2),grid%loc_n_r,3))
1256  ierr = calc_xyz_grid(grid_eq,grid,xyz(:,:,:,1),&
1257  &xyz(:,:,:,2),xyz(:,:,:,3))
1258  chckerr('')
1259  allocate(xyz_i(grid%n(1),grid%n(2),grid%loc_n_r,3))
1260  do id = 1,grid%n(1)
1261  xyz_i(id,:,:,1) = id - 1
1262  end do
1263  do jd = 1,grid%n(2)
1264  xyz_i(:,jd,:,2) = jd - 1
1265  end do
1266  do kd = 1,grid%loc_n_r
1267  xyz_i(:,:,kd,3) = kd - 1
1268  end do
1269  else
1270  allocate(xyz_i(1,grid%loc_n_r,1,1))
1271  do kd = 1,grid%loc_n_r
1272  xyz_i(1,kd,1,1) = kd - 1
1273  end do
1274  end if
1275 
1276  ! plot
1277  description = 'total range'
1278  select case (test_case)
1279  case (1) ! eq_2
1280  call plot_hdf5('jac_FD','jac_FD_tot',&
1281  &eq_2%jac_FD(:,:,:,0,0,0),&
1282  &descr=description,&
1283  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3))
1284  call plot_hdf5('kappa_n','kappa_n_tot',eq_2%kappa_n,&
1285  &descr=description,&
1286  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3))
1287  case (2) ! X_1
1288  call plot_hdf5('U_0_IM','U_0_IM_tot',&
1289  &ip(x_1%U_0(:,:,:,1)),&
1290  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1291  &descr=description)
1292  call plot_hdf5('DU_1_RE','DU_1_RE_tot',&
1293  &rp(x_1%DU_1(:,:,:,1)),&
1294  &x=xyz(:,:,:,1),y=xyz(:,:,:,2),z=xyz(:,:,:,3),&
1295  &descr=description)
1296  case (3) ! sol
1297  call print_ex_2d(['sol_RE'],'sol_RE_tot',&
1298  &rp(sol%vec(1,:,:)),x=xyz_i(1,:,:,1),draw=.false.)
1299  call draw_ex(['sol_RE'],'sol_RE_tot',&
1300  &size(sol%vec,3),1,.false.)
1301  end select
1302  if (test_case.eq.1 .or. test_case.eq.2) then
1303  do id = 1,3
1304  call plot_hdf5(xyz_name(id),xyz_name(id)//'_tot',&
1305  &xyz(:,:,:,id),&
1306  &x=xyz_i(:,:,:,1),y=xyz_i(:,:,:,2),z=xyz_i(:,:,:,3))
1307  end do
1308  end if
1309 
1310  ! clean up
1311  call eq_2%dealloc()
1312  if (test_case.eq.1 .or. test_case.eq.2) deallocate(xyz)
1313  deallocate(xyz_i)
1314  call lvl_ud(-1)
1315  end if
1316 
1317  ! clean up
1318  call grid%dealloc()
1319  if (test_case.eq.2 .or. test_case.eq.3) call grid_eq%dealloc()
1320  call dealloc_in()
1321 
1322  ! synchronize MPI
1323  ierr = wait_mpi()
1324  chckerr('')
1325 
1326  ! user output
1327  call lvl_ud(-1)
1328  call writo('Test complete')
1329  end function test_read_hdf5_subset
1330 
1334  integer function test_calc_int_vol() result(ierr)
1335  use num_vars, only: rank
1336  use grid_utilities, only: calc_eqd_grid, calc_int_vol, &
1337  &debug_calc_int_vol
1338 
1339  character(*), parameter :: rout_name = 'test_calc_int_vol'
1340 
1341  ! local variables
1342  real(dp), allocatable :: theta(:,:,:), zeta(:,:,:) ! angles
1343  complex(dp), allocatable :: fun(:,:,:,:) ! function to integrate
1344  real(dp), allocatable :: j(:,:,:) ! Jacobian
1345  real(dp), allocatable :: x(:,:,:), y(:,:,:), z(:,:,:) ! X,Y and Z for plot
1346  real(dp), allocatable :: norm(:) ! normal variable
1347  real(dp), allocatable :: r(:) ! radius
1348  real(dp) :: r_0 ! geometric axis of torus
1349  complex(dp), allocatable :: fun_int(:,:,:) ! integrated function
1350  integer, allocatable :: dims(:,:) ! dimensions of grid
1351  integer :: n_steps ! nr. of steps
1352  integer :: id, jd, kd ! counters
1353  real(dp), allocatable :: eqd_grid_dum(:) ! equidistant grid dummy
1354  character(len=max_str_ln) :: var_name(3), file_name(3) ! name of variable and of file
1355  complex(dp), allocatable :: fun_int_plot(:,:) ! results for plotting
1356 
1357  ! initialize ierr
1358  ierr = 0
1359 
1360  if (rank.eq.0) then
1361  call writo('Checking integral of')
1362  call lvl_ud(1)
1363  call writo('f = 1 - r^2 + i cos(theta)')
1364  call lvl_ud(-1)
1365  call writo('on a torus')
1366 
1367  call writo('Debug mode?')
1368  if (get_log(.false.)) then
1369  debug_calc_int_vol = .true.
1370  else
1371  debug_calc_int_vol = .false.
1372  end if
1373 
1374  call writo('How many different grid sizes?')
1375  n_steps = get_int(lim_lo=1)
1376 
1377  allocate(dims(3,n_steps))
1378  allocate(eqd_grid_dum(n_steps))
1379 
1380  ! get dimensions
1381  if (n_steps.gt.1) then
1382  do kd = 1,size(dims,1)
1383  call writo('Starting dimension '//trim(i2str(kd)))
1384  dims(kd,1) = get_int(lim_lo=4)
1385  call writo('Final dimension '//trim(i2str(kd)))
1386  dims(kd,n_steps) = get_int(lim_lo=dims(kd,1)+n_steps)
1387  ierr = calc_eqd_grid(eqd_grid_dum,1._dp*dims(kd,1),&
1388  &1._dp*dims(kd,n_steps))
1389  chckerr('')
1390  dims(kd,:) = nint(eqd_grid_dum)
1391  end do
1392  call writo('Performing calculations for '//&
1393  &trim(i2str(n_steps))//' grid sizes')
1394  else
1395  do kd = 1,size(dims)
1396  call writo('Dimension '//trim(i2str(kd)))
1397  dims(kd,1) = get_int(lim_lo=1)
1398  if (dims(kd,1).eq.1) then
1399  call lvl_ud(1)
1400  call writo('Make sure the integrand does not depend on &
1401  &dimension '//trim(i2str(kd)))
1402  call writo('A factor such as 2pi can be missing from &
1403  &resulting integral')
1404  call lvl_ud(-1)
1405  end if
1406  end do
1407  end if
1408 
1409  allocate(fun_int(2,1,n_steps))
1410 
1411  do jd = 1,n_steps
1412  ! user output
1413  call writo('grid size = ['//trim(i2str(dims(1,jd)))//' ,'//&
1414  &trim(i2str(dims(2,jd)))//', '//trim(i2str(dims(3,jd)))//&
1415  &']')
1416  call lvl_ud(1)
1417 
1418  ! set up variables
1419  allocate(theta(dims(1,jd),dims(2,jd),dims(3,jd)))
1420  allocate(zeta(dims(1,jd),dims(2,jd),dims(3,jd)))
1421  allocate(j(dims(1,jd),dims(2,jd),dims(3,jd)))
1422  allocate(fun(dims(1,jd),dims(2,jd),dims(3,jd),1))
1423  allocate(x(dims(1,jd),dims(2,jd),dims(3,jd)))
1424  allocate(y(dims(1,jd),dims(2,jd),dims(3,jd)))
1425  allocate(z(dims(1,jd),dims(2,jd),dims(3,jd)))
1426  allocate(norm(dims(3,jd)))
1427  allocate(r(dims(3,jd)))
1428 
1429  r_0 = 2._dp
1430  ierr = calc_eqd_grid(theta,0*pi,2*pi,1)
1431  chckerr('')
1432  ierr = calc_eqd_grid(zeta,0*pi,2*pi,2)
1433  chckerr('')
1434  ierr = calc_eqd_grid(norm,0._dp,1._dp)
1435  chckerr('')
1436  do kd = 1,dims(3,jd)
1437  r(kd) = (kd-1._dp)/(dims(3,jd)-1)
1438  fun(:,:,kd,1) = 1._dp - r(kd)**2 + iu*cos(theta(:,:,kd))
1439  do id = 1,dims(1,jd)
1440  j(id,:,kd) = r(kd)*(r_0+r(kd)*cos(theta(id,:,kd)))
1441  end do
1442  x(:,:,kd) = cos(zeta(:,:,kd))*(r_0+r(kd)*cos(theta(:,:,kd)))
1443  y(:,:,kd) = sin(zeta(:,:,kd))*(r_0+r(kd)*cos(theta(:,:,kd)))
1444  z(:,:,kd) = r(kd)*sin(theta(:,:,kd))
1445  end do
1446 
1447  var_name(1) = 'TEST_J_torus'
1448  var_name(2) = 'TEST_RE_fun_torus'
1449  var_name(3) = 'TEST_IM_fun_torus'
1450  file_name(1) = 'TEST_J_torus'
1451  file_name(2) = 'TEST_RE_fun_torus'
1452  file_name(3) = 'TEST_IM_fun_torus'
1453  if (n_steps.gt.1) then
1454  do kd = 1,3
1455  file_name(kd) = trim(file_name(kd))//'_'//trim(i2str(jd))
1456  end do
1457  end if
1458  call plot_hdf5(var_name(1),file_name(1),j,x=x,y=y,z=z)
1459  call plot_hdf5(var_name(2),file_name(2),rp(fun(:,:,:,1)),&
1460  &x=x,y=y,z=z)
1461  call plot_hdf5(var_name(3),file_name(3),ip(fun(:,:,:,1)),&
1462  &x=x,y=y,z=z)
1463 
1464  ! user output
1465  call writo('Analytically calculating integral')
1466 
1467  ! calculate integral analytically
1468  fun_int(1,1,jd) = r_0*pi**2 + iu*2*pi**2/3
1469 
1470  ! user output
1471  call writo('Numerically calculating integral')
1472 
1473  ! calculate integral numerically
1474  ierr = calc_int_vol(theta,zeta,norm,j,fun,fun_int(2,:,jd))
1475  chckerr('')
1476 
1477  ! output
1478  call writo('Analytical value: '//&
1479  &trim(c2str(fun_int(1,1,jd))))
1480  call writo('Numerical value: '//&
1481  &trim(c2str(fun_int(2,1,jd))))
1482 
1483  ! clean up
1484  deallocate(theta,zeta,j,fun)
1485  deallocate(x,y,z)
1486  deallocate(norm,r)
1487 
1488  call lvl_ud(-1)
1489  end do
1490 
1491  if (n_steps.gt.1) then
1492  allocate(fun_int_plot(n_steps,2))
1493  do jd = 1,n_steps
1494  fun_int_plot(jd,:) = fun_int(:,1,jd)/fun_int(1,1,jd)
1495  end do
1496  call print_ex_2d(['analytical integral, RE',&
1497  &'numerical integral, RE '],'',rp(fun_int_plot))
1498  call print_ex_2d(['analytical integral, IM',&
1499  &'numerical integral, IM '],'',ip(fun_int_plot))
1500  end if
1501  end if
1502  end function test_calc_int_vol
1503 #endif
1504 end module test
pb3d_ops::reconstruct_pb3d_sol
integer function, public reconstruct_pb3d_sol(mds, grid_sol, sol, data_name, rich_lvl, lim_sec_sol, lim_pos)
Reconstructs the solution variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:1359
input_utilities::get_real
real(dp) function, public get_real(lim_lo, lim_hi, ind)
Queries for user input for a real value, where allowable range can be provided as well.
Definition: input_utilities.f90:77
grid_utilities::debug_calc_int_vol
logical, public debug_calc_int_vol
plot debug information for calc_int_vol()
Definition: grid_utilities.f90:25
dtorh
Calculation of toroidal functions and .
Definition: dtorh.f90:8
x_ops
Operations considering perturbation quantities.
Definition: X_ops.f90:4
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition: eq_vars.f90:27
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
input_utilities::get_log
logical function, public get_log(yes, ind)
Queries for a logical value yes or no, where the default answer is also to be provided.
Definition: input_utilities.f90:22
mpi_utilities::lock_header
character(len=max_str_ln) function, public lock_header(lock_loc)
Returns the header for lock debug messages.
Definition: MPI_utilities.f90:1217
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
num_utilities::calc_d2_smooth
integer function, public calc_d2_smooth(fil_N, x, y, D2y, style)
Calculate second derivative with smoothing formula by Holoborodko, .
Definition: num_utilities.f90:2829
rich_vars
Variables concerning Richardson extrapolation.
Definition: rich_vars.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
mpi_utilities::lock_req_acc
integer function, public lock_req_acc(lock, blocking)
Request access to lock of a BL (blocking) or optionally a NB (non-blocking) type.
Definition: MPI_utilities.f90:765
x_ops::init_modes
integer function, public init_modes(grid_eq, eq)
Initializes some variables concerning the mode numbers.
Definition: X_ops.f90:919
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
vmec_utilities::calc_trigon_factors
integer function, public calc_trigon_factors(theta, zeta, trigon_factors)
Calculate the trigonometric cosine and sine factors.
Definition: VMEC_utilities.f90:275
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
grid_vars::n_r_eq
integer, public n_r_eq
nr. of normal points in equilibrium grid
Definition: grid_vars.f90:20
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
grid_utilities::calc_int_vol
integer function, public calc_int_vol(ang_1, ang_2, norm, J, f, f_int)
Calculates volume integral on a 3D grid.
Definition: grid_utilities.f90:1320
num_vars::prog_style
integer, public prog_style
program style (1: PB3D, 2: PB3D_POST)
Definition: num_vars.f90:53
pb3d_ops::reconstruct_pb3d_in
integer function, public reconstruct_pb3d_in(data_name)
Reconstructs the input variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:41
test::test_calc_int_vol
integer function test_calc_int_vol()
Tests calculation of volume integral.
Definition: test.f90:1335
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
test::test_splines
integer function test_splines()
Test spline().
Definition: test.f90:188
mpi_vars
Variables pertaining to MPI.
Definition: MPI_vars.f90:4
sol_vars::sol_type
solution type
Definition: sol_vars.f90:30
mpi_vars::lock_type
lock type
Definition: MPI_vars.f90:63
test
Generic tests.
Definition: test.f90:4
pb3d_ops::reconstruct_pb3d_eq_2
integer function, public reconstruct_pb3d_eq_2(grid_eq, eq, data_name, rich_lvl, tot_rich, lim_pos)
Reconstructs the equilibrium variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:695
plot_y
subroutine plot_y(x, x_int, y, y_int, deriv)
Definition: test.f90:309
test::test_read_hdf5_subset
integer function test_read_hdf5_subset()
Tests reading of HDF5 subset.
Definition: test.f90:806
num_utilities::spline
Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the m...
Definition: num_utilities.f90:276
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
mpi_utilities::debug_lock
logical, public debug_lock
print debug information about lock operations
Definition: MPI_utilities.f90:40
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
vmec_utilities
Numerical utilities related to the output of VMEC.
Definition: VMEC_utilities.f90:4
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition: num_vars.f90:89
x_ops::setup_modes
integer function, public setup_modes(mds, grid_eq, grid, plot_name)
Sets up some variables concerning the mode numbers.
Definition: X_ops.f90:1060
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
mpi_utilities::lock_wl_change
integer function, public lock_wl_change(wl_action, blocking, lock_loc, wl, ranks)
Adds, removes or sets to active a rank from the waiting list for a lock and returns the lock waiting ...
Definition: MPI_utilities.f90:1111
x_vars::modes_type
mode number type
Definition: X_vars.f90:36
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
mpi_utilities::wait_mpi
integer function, public wait_mpi()
Wait for all processes, wrapper to MPI barrier.
Definition: MPI_utilities.f90:744
grid_utilities::copy_grid
integer function, public copy_grid(grid_A, grid_B, lims_B, i_lim)
Copy a grid A to a new grid B, that was not yet initialized.
Definition: grid_utilities.f90:1189
grid_utilities::calc_eqd_grid
Calculate grid of equidistant points, where optionally the last point can be excluded.
Definition: grid_utilities.f90:75
dtorh::dtorh1
integer function, public dtorh1(z, m, nmax, pl, ql, newn, mode, ipre)
Calculates toroidal harmonics of a fixed order m and degrees up to nmax.
Definition: dtorh.f90:75
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
pb3d_ops::reconstruct_pb3d_x_1
integer function, public reconstruct_pb3d_x_1(mds, grid_X, X, data_name, rich_lvl, tot_rich, lim_sec_X, lim_pos)
Reconstructs the vectorial perturbation variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:833
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
vmec_vars::z_v_s
real(dp), dimension(:,:,:), allocatable, public z_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:42
pb3d_ops
Operations on PB3D output.
Definition: PB3D_ops.f90:8
num_vars::pi
real(dp), parameter, public pi
Definition: num_vars.f90:83
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition: grid_utilities.f90:4
pb3d_ops::reconstruct_pb3d_eq_1
integer function, public reconstruct_pb3d_eq_1(grid_eq, eq, data_name, lim_pos)
Reconstructs the equilibrium variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:597
vmec_vars::r_v_c
real(dp), dimension(:,:,:), allocatable, public r_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:39
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
pb3d_ops::reconstruct_pb3d_grid
integer function, public reconstruct_pb3d_grid(grid, data_name, rich_lvl, tot_rich, lim_pos, grid_limits)
Reconstructs grid variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:442
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
input_utilities::get_int
integer function, public get_int(lim_lo, lim_hi, ind)
Queries for user input for an integer value, where allowable range can be provided as well.
Definition: input_utilities.f90:152
input_utilities::pause_prog
subroutine, public pause_prog(ind)
Pauses the running of the program.
Definition: input_utilities.f90:226
vmec_vars
Variables that concern the output of VMEC.
Definition: VMEC_vars.f90:4
vmec_vars::l_v_c
real(dp), dimension(:,:,:), allocatable, public l_v_c
Coeff. of in sine series (HM) and norm. deriv.
Definition: VMEC_vars.f90:43
vmec_vars::r_v_s
real(dp), dimension(:,:,:), allocatable, public r_v_s
Coeff. of in cosine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:40
rich_vars::rich_lvl
integer, public rich_lvl
current level of Richardson extrapolation
Definition: rich_vars.f90:19
x_vars::x_1_type
vectorial perturbation type
Definition: X_vars.f90:51
sol_vars
Variables pertaining to the solution quantities.
Definition: sol_vars.f90:4
test::test_calc_d2_smooth
integer function test_calc_d2_smooth()
Test calc_d2_smooth().
Definition: test.f90:106
mpi_utilities::lock_return_acc
integer function, public lock_return_acc(lock)
Returns access to a lock.
Definition: MPI_utilities.f90:872
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
input_utilities::dealloc_in
subroutine, public dealloc_in()
Cleans up input from equilibrium codes.
Definition: input_utilities.f90:268
str_utilities::c2str
elemental character(len=max_str_ln) function, public c2str(k)
Convert a complex (double) to string.
Definition: str_utilities.f90:66
input_utilities
Numerical utilities related to input.
Definition: input_utilities.f90:4
test::generic_tests
integer function, public generic_tests()
Performs generic tests.
Definition: test.f90:39
vmec_utilities::fourier2real
Inverse Fourier transformation, from VMEC.
Definition: VMEC_utilities.f90:56
vmec_vars::l_v_s
real(dp), dimension(:,:,:), allocatable, public l_v_s
Coeff. of in cosine series (HM) and norm. deriv.
Definition: VMEC_vars.f90:44
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
vmec_vars::z_v_c
real(dp), dimension(:,:,:), allocatable, public z_v_c
Coeff. of in sine series (FM) and norm. deriv.
Definition: VMEC_vars.f90:41
test::test_tor_fun
integer function test_tor_fun()
Test calculation of toroidal functions.
Definition: test.f90:664
grid_utilities::calc_xyz_grid
integer function, public calc_xyz_grid(grid_eq, grid_XYZ, X, Y, Z, L, R)
Calculates , and on a grid grid_XYZ, determined through its E(quilibrium) coordinates.
Definition: grid_utilities.f90:799
grid_utilities::trim_grid
integer function, public trim_grid(grid_in, grid_out, norm_id)
Trim a grid, removing any overlap between the different regions.
Definition: grid_utilities.f90:1636
this documentation was created with Doxygen 1.8.17 by Toon Weyens