PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
test.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Generic tests.
3!------------------------------------------------------------------------------!
4module 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
20contains
21 !> Performs generic tests.
22 !!
23 !! These tests are run in interactive fashion by running the program with
24 !! the <tt>--test</tt> command-line option.
25 !! \see init_files()
26 !!
27 !! \note
28 !! -# It is probably also recommended to run it with
29 !! <tt>--do_execute_command_line</tt> in order to generate the plots in
30 !! real time.
31 !! -# It has not been tested whether testing routines "pollutes" the
32 !! further simulations. It is therefore best to not use the tests when
33 !! doing a real run.
34 !!
35 !! \ldebug
36 !!
37 !! \return ierr
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
102 !> Test calc_d2_smooth().
103 !!
104 !! \return ierr
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
184 !> Test spline().
185 !!
186 !! \return ierr
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
353 !> Test lock system.
354 !!
355 !! \return ierr
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
486 !> \private
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
512 !> Tests the calculation of R, Z, and L.
513 !!
514 !! \return ierr
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
540 !> \private
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
660 !> Test calculation of toroidal functions.
661 !!
662 !! \return ierr
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
802 !> Tests reading of HDF5 subset.
803 !!
804 !! \return ierr
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
1331 !> Tests calculation of volume integral.
1332 !!
1333 !! \return ierr
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
1504end module test
grid_utilities::calc_eqd_grid
Calculate grid of equidistant points, where optionally the last point can be excluded.
Definition grid_utilities.f90:75
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::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
output_ops::print_ex_2d
Print 2-D output on a file.
Definition output_ops.f90:47
vmec_utilities::fourier2real
Inverse Fourier transformation, from VMEC.
Definition VMEC_utilities.f90:56
dtorh
Calculation of toroidal functions and .
Definition dtorh.f90:8
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
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition eq_vars.f90:27
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition grid_utilities.f90:4
grid_utilities::debug_calc_int_vol
logical, public debug_calc_int_vol
plot debug information for calc_int_vol()
Definition grid_utilities.f90:25
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
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::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_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
grid_vars
Variables pertaining to the different grids used.
Definition grid_vars.f90:4
grid_vars::n_r_eq
integer, public n_r_eq
nr. of normal points in equilibrium grid
Definition grid_vars.f90:20
input_utilities
Numerical utilities related to input.
Definition input_utilities.f90:4
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::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
input_utilities::pause_prog
subroutine, public pause_prog(ind)
Pauses the running of the program.
Definition input_utilities.f90:226
input_utilities::dealloc_in
subroutine, public dealloc_in()
Cleans up input from equilibrium codes.
Definition input_utilities.f90:268
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
messages
Numerical utilities related to giving output.
Definition messages.f90:4
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition messages.f90:254
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition messages.f90:275
mpi_utilities
Numerical utilities related to MPI.
Definition MPI_utilities.f90:20
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
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
mpi_utilities::wait_mpi
integer function, public wait_mpi()
Wait for all processes, wrapper to MPI barrier.
Definition MPI_utilities.f90:744
mpi_utilities::lock_return_acc
integer function, public lock_return_acc(lock)
Returns access to a lock.
Definition MPI_utilities.f90:872
mpi_utilities::debug_lock
logical, public debug_lock
print debug information about lock operations
Definition MPI_utilities.f90:40
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
mpi_vars
Variables pertaining to MPI.
Definition MPI_vars.f90:4
num_utilities
Numerical utilities.
Definition num_utilities.f90:4
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, holoborodko2008diff.
Definition num_utilities.f90:2829
num_vars
Numerical variables used by most other modules.
Definition num_vars.f90:4
num_vars::dp
integer, parameter, public dp
double precision
Definition num_vars.f90:46
num_vars::pi
real(dp), parameter, public pi
Definition num_vars.f90:83
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition num_vars.f90:85
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
num_vars::prog_style
integer, public prog_style
program style (1: PB3D, 2: PB3D_POST)
Definition num_vars.f90:53
num_vars::max_deriv
integer, parameter, public max_deriv
highest derivatives for metric factors in Flux coords.
Definition num_vars.f90:52
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition num_vars.f90:89
num_vars::rank
integer, public rank
MPI rank.
Definition num_vars.f90:68
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition num_vars.f90:133
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
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
pb3d_ops
Operations on PB3D output.
Definition PB3D_ops.f90:8
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
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
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
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
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
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
rich_vars
Variables concerning Richardson extrapolation.
Definition rich_vars.f90:4
rich_vars::rich_lvl
integer, public rich_lvl
current level of Richardson extrapolation
Definition rich_vars.f90:19
sol_vars
Variables pertaining to the solution quantities.
Definition sol_vars.f90:4
str_utilities
Operations on strings.
Definition str_utilities.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
str_utilities::c2str
elemental character(len=max_str_ln) function, public c2str(k)
Convert a complex (double) to string.
Definition str_utilities.f90:66
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition str_utilities.f90:42
test
Generic tests.
Definition test.f90:4
test::generic_tests
integer function, public generic_tests()
Performs generic tests.
Definition test.f90:39
vmec_utilities
Numerical utilities related to the output of VMEC.
Definition VMEC_utilities.f90:4
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
vmec_vars
Variables that concern the output of VMEC.
Definition VMEC_vars.f90:4
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
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
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
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
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
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
x_ops
Operations considering perturbation quantities.
Definition X_ops.f90:4
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_ops::init_modes
integer function, public init_modes(grid_eq, eq)
Initializes some variables concerning the mode numbers.
Definition X_ops.f90:919
x_vars
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
eq_vars::eq_1_type
flux equilibrium type
Definition eq_vars.f90:63
eq_vars::eq_2_type
metric equilibrium type
Definition eq_vars.f90:114
grid_vars::grid_type
Type for grids.
Definition grid_vars.f90:59
mpi_vars::lock_type
lock type
Definition MPI_vars.f90:63
sol_vars::sol_type
solution type
Definition sol_vars.f90:30
x_vars::modes_type
mode number type
Definition X_vars.f90:36
x_vars::x_1_type
vectorial perturbation type
Definition X_vars.f90:51
plot_y
subroutine plot_y(x, x_int, y, y_int, deriv)
Definition test.f90:309
this documentation was created with Doxygen 1.16.1 by Toon Weyens