PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
vac_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
22 !------------------------------------------------------------------------------!
23 module vac_ops
24 #include <PB3D_macros.h>
25 #include <wrappers.h>
26  use strumpackdensepackage
27  use str_utilities
28  use messages
29  use output_ops
30  use num_vars, only: dp, pi, max_str_ln, iu
31  use grid_vars, only: grid_type
32  use eq_vars, only: eq_1_type, eq_2_type
33  use x_vars, only: x_2_type, modes_type
34  use vac_vars, only: copy_vac, &
35  &vac_type
36 
37  implicit none
38  private
39  public store_vac, calc_vac_res, print_output_vac, vac_pot_plot
40 #if ldebug
41  public debug_calc_gh, debug_calc_vac_res, debug_vac_pot_plot
42 #endif
43 
44  ! global variables
45 #if ldebug
46  logical :: debug_calc_gh = .false.
47  logical :: debug_calc_vac_res = .false.
48  logical :: debug_vac_pot_plot = .false.
49 #endif
50 
51 contains
112  integer function store_vac(grid,eq_1,eq_2,vac) result(ierr)
113  use pb3d_ops, only: reconstruct_pb3d_vac
114  use num_vars, only: rich_restart_lvl, eq_job_nr, eq_jobs_lims, &
115  &eq_style, rank, n_procs, use_pol_flux_f
116  use rich_vars, only: rich_lvl
117  use mpi_utilities, only: broadcast_var
118  use num_utilities, only: c
119  use grid_vars, only: n_r_eq
120  use x_vars, only: prim_x
121 
122  character(*), parameter :: rout_name = 'store_vac'
123 
124  ! input / output
125  type(grid_type), intent(in) :: grid
126  type(eq_1_type), intent(in) :: eq_1
127  type(eq_2_type), intent(in) :: eq_2
128  type(vac_type), intent(inout) :: vac
129 
130  ! local variables
131  integer :: id, jd ! counter
132  integer :: r_id ! normal index
133  real(dp) :: jq ! iota or q
134  character(len=max_str_ln) :: err_msg ! error message
135 
136  ! initialize ierr
137  ierr = 0
138 
139  ! set jq
140  if (use_pol_flux_f) then
141  jq = eq_1%q_saf_FD(grid%loc_n_r,0)
142  else
143  jq = eq_1%rot_t_FD(grid%loc_n_r,0)
144  ierr = 1
145  err_msg = 'TOR. FLUX HAS NOT BEEN TESTED AND IS MOST PROBABLY NOT &
146  &CORRECT YET'
147  chckerr(err_msg)
148  end if
149  ierr = broadcast_var(jq,n_procs-1)
150  chckerr('')
151 
152  ! call the approriate procedure
153  select case (eq_style)
154  case (1) ! VMEC
155  ierr = store_vac_vmec()
156  chckerr('')
157  case (2) ! HELENA
158  ierr = store_vac_hel()
159  chckerr('')
160  end select
161  contains
162  ! VMEC version
164  integer function store_vac_vmec() result(ierr)
165  use num_vars, only: jump_to_sol
166  use grid_vars, only: n_alpha
167  use rich_vars, only: n_par_x
168  use eq_utilities, only: calc_inv_met
169  use vac_utilities, only: interlaced_vac_copy
170 
171  character(*), parameter :: rout_name = 'store_vac_VMEC'
172 
173  ! local variables
174  integer :: par_id(3) ! total parallel index
175  real(dp), allocatable :: norm_com_c(:,:,:) ! Cylindrical components of norm
176  real(dp), allocatable :: t_cv_loc(:,:,:,:,:,:,:) ! local T_CV
177  type(vac_type) :: vac_old ! old vacuum variables
178  logical :: interlaced_vac_copy_needed ! interlaced copy vacuum needed
179 
180  ! initialize ierr
181  ierr = 0
182 
183  ! user output
184  call writo('Start storing vacuum quantities')
185 
186  call lvl_ud(1)
187 
188  ! if start of Richardson level that is not the first, copy previous
189  ! results
190  if (eq_job_nr.eq.1) then
191  if (rich_lvl.eq.rich_restart_lvl) then ! restarting
192  if (jump_to_sol) then ! jumping straight to solution
193  ! nothing needs to be done here
194  return
195  else ! not jumping to solution
196  if (rich_restart_lvl.eq.1) then ! not actually restarting, just starting
197  ! nothing extra needed
198  interlaced_vac_copy_needed = .false.
199  else ! restarting
200  ! reconstruct old vacuum
201  ierr = reconstruct_pb3d_vac(vac_old,'vac',&
202  &rich_lvl=rich_lvl-1)
203  chckerr('')
204 
205  ! interlaced copy necessary
206  interlaced_vac_copy_needed = .true.
207  end if
208  end if
209  else ! not restarting
210  ! copy old vacuum temporarily
211  ierr = copy_vac(vac,vac_old)
212  chckerr('')
213 
214  ! deallocate
215  call vac%dealloc()
216 
217  ! interlaced copy necessary
218  interlaced_vac_copy_needed = .true.
219  end if
220 
221  ! allocate
222  ierr = vac%init(1,n_par_x*n_alpha,prim_x,[n_par_x,n_alpha],jq)
223  chckerr('')
224 
225  ! interlaced copy if needed
226  if (interlaced_vac_copy_needed) then
227  ierr = interlaced_vac_copy(vac_old,vac)
228  chckerr('')
229  end if
230  end if
231 
232  ! add results from current equilibrium job to the variables
233  if (rank.eq.n_procs-1) then
234  ! grid point to save is last local
235  r_id = grid%loc_n_r
236 
237  ! calculate Cartesian components of J nabla psi
238  ! = J_F (Dpsi_pol/Dr_V) / 2pi (T_C^V)^T (e^C)
239  ! = J_F T_F^V(1,2) (T_C^V)^T (e^C)
240  allocate(norm_com_c(grid%n(1),grid%n(2),3))
241  norm_com_c = 0._dp
242  allocate(t_cv_loc(size(eq_2%T_VC,1),size(eq_2%T_VC,2),1:1,&
243  &size(eq_2%T_VC,4),0:0,0:0,0:0))
244  ierr = calc_inv_met(t_cv_loc,eq_2%T_VC(:,:,r_id:r_id,:,:,:,:),&
245  &[0,0,0])
246  chckerr('')
247  do id = 1,3
248  norm_com_c(:,:,id) = eq_2%jac_FD(:,:,r_id,0,0,0) * &
249  &eq_2%T_EF(:,:,r_id,c([1,2],.false.),0,0,0)*&
250  &t_cv_loc(:,:,1,c([id,1],.false.),0,0,0) ! Cylindrical contravariant components
251  end do
252  deallocate(t_cv_loc)
253 
254  ! iterate over all field lines
255  do jd = 1,n_alpha
256  ! set total parallel id
257  if (rich_lvl.eq.1) then
258  par_id = [eq_jobs_lims(:,eq_job_nr),1]
259  else
260  par_id = [2*eq_jobs_lims(:,eq_job_nr),2]
261  end if
262  par_id(1:2) = par_id(1:2) + (jd-1)*n_par_x
263 
264  ! save X, Y and Z
265  vac%x_vec(par_id(1):par_id(2):par_id(3),1) = &
266  &eq_2%R_E(:,jd,r_id,0,0,0) * cos(grid%zeta_E(:,jd,r_id))
267  vac%x_vec(par_id(1):par_id(2):par_id(3),2) = &
268  &eq_2%R_E(:,jd,r_id,0,0,0) * sin(grid%zeta_E(:,jd,r_id))
269  vac%x_vec(par_id(1):par_id(2):par_id(3),3) = &
270  &eq_2%Z_E(:,jd,r_id,0,0,0)
271 
272  ! save norm
273  vac%norm(par_id(1):par_id(2):par_id(3),1) = &
274  &norm_com_c(:,jd,1)*cos(grid%zeta_E(:,jd,r_id)) - &
275  &norm_com_c(:,jd,2)*sin(grid%zeta_E(:,jd,r_id)) / &
276  &eq_2%R_E(:,jd,r_id,0,0,0) ! ~ e^X
277  vac%norm(par_id(1):par_id(2):par_id(3),2) = &
278  &norm_com_c(:,jd,1)*sin(grid%zeta_E(:,jd,r_id)) + &
279  &norm_com_c(:,jd,2)*cos(grid%zeta_E(:,jd,r_id)) / &
280  &eq_2%R_E(:,jd,r_id,0,0,0) ! ~ e^X
281  vac%norm(par_id(1):par_id(2):par_id(3),3) = &
282  &norm_com_c(:,jd,3) ! ~ e^Z
283 
284  ! save metric factors and H factor
285  vac%h_fac(par_id(1):par_id(2):par_id(3),1) = &
286  &eq_2%g_FD(:,jd,r_id,c([1,1],.true.),0,0,0) ! g_11
287  vac%h_fac(par_id(1):par_id(2):par_id(3),2) = &
288  &eq_2%g_FD(:,jd,r_id,c([1,3],.true.),0,0,0) ! g_13
289  vac%h_fac(par_id(1):par_id(2):par_id(3),3) = &
290  &eq_2%g_FD(:,jd,r_id,c([3,3],.true.),0,0,0) ! g_33
291  vac%h_fac(par_id(1):par_id(2):par_id(3),4) = 0.5_dp* &
292  &eq_2%jac_FD(:,jd,r_id,0,0,0)**2 * &
293  &(eq_2%h_FD(:,jd,r_id,c([2,2],.true.),0,0,0) / &
294  &eq_2%g_FD(:,jd,r_id,c([1,1],.true.),0,0,0))**1.5 * ( &
295  &eq_2%jac_FD(:,jd,r_id,0,1,0)/&
296  &eq_2%jac_FD(:,jd,r_id,0,0,0) + &
297  &0.5_dp*eq_2%h_FD(:,jd,r_id,c([2,2],.true.),0,1,0)/&
298  &eq_2%h_FD(:,jd,r_id,c([2,2],.true.),0,0,0) - &
299  &0.5_dp*eq_2%g_FD(:,jd,r_id,c([1,1],.true.),0,1,0)/&
300  &eq_2%g_FD(:,jd,r_id,c([1,1],.true.),0,0,0) &
301  &) ! J/2 h^22/g_11 d/d2 (J sqrt(h^22/g_11)))
302  end do
303  !call print_ex_2D(['norm'],'norm_'//trim(i2str(rich_lvl)),&
304  !&vac%norm,persistent=.true.)
305  !call print_ex_2D(['x_vec'],'x_vec'//trim(i2str(rich_lvl)),&
306  !&vac%x_vec,persistent=.true.)
307  end if
308 
309  ! broadcast to other processes
310  do id = 1,size(vac%x_vec,2)
311  ierr = broadcast_var(vac%x_vec(:,id),n_procs-1)
312  chckerr('')
313  ierr = broadcast_var(vac%norm(:,id),n_procs-1)
314  chckerr('')
315  end do
316  do id = 1,size(vac%h_fac,2)
317  ierr = broadcast_var(vac%h_fac(:,id),n_procs-1)
318  chckerr('')
319  end do
320 
321  call lvl_ud(-1)
322 
323  call writo('Done storing vacuum quantities')
324  end function store_vac_vmec
325  ! HELENA version
327  integer function store_vac_hel() result(ierr)
328  use helena_vars, only: r_h, z_h, nchi, chi_h, ias
329  use grid_utilities, only: nufft
330 
331  character(*), parameter :: rout_name = 'store_vac_HEL'
332 
333  ! local variables
334  integer :: n_theta ! total number of points, including overlap
335  real(dp), allocatable :: x_vec_f(:,:) ! Fourier components of x_vec
336  character(len=max_str_ln) :: err_msg ! error message
337 
338  ! initialize ierr
339  ierr = 0
340 
341  ierr = 2
342  chckerr('Vacuum has not been implemented yet!')
343 
344  ! user output
345  call writo('Start storing vacuum quantities')
346 
347  call lvl_ud(1)
348 
349  ! if restarting a Richardson level, reconstruct. If not, only
350  ! calculate for the first level
351  if (rich_lvl.eq.rich_restart_lvl) then
352  ! check for start from zero versus restart
353  if (rich_restart_lvl.eq.1) then ! starting from zero
354  ! test
355  if (eq_job_nr.ne.1) then
356  ierr = 1
357  err_msg = 'HELENA should only have 1 equilibrium job &
358  &for Richardson level 1'
359  chckerr(err_msg)
360  end if
361 
362  if (ias.eq.0) then ! top-bottom symmmetric
363  n_theta = 2*(nchi-1)+1
364  else
365  n_theta = nchi
366  end if
367 
368  ! allocate
369  ierr = vac%init(eq_style,n_theta,prim_x,[n_theta,1],jq)
370  chckerr('')
371  ! save points
372  if (rank.eq.n_procs-1) then
373  ! grid point to save is last total
374  r_id = n_r_eq
375 
376  ! save only R and Z
377  if (ias.eq.0) then
378  vac%x_vec(1:nchi,1) = r_h(1:nchi,r_id)
379  vac%x_vec(1:nchi,2) = z_h(1:nchi,r_id)
380  vac%ang(1:nchi,1) = chi_h(1:nchi)
381  vac%x_vec(nchi+1:n_theta,1) = r_h(nchi-1:1:-1,r_id)
382  vac%x_vec(nchi+1:n_theta,2) = -z_h(nchi-1:1:-1,r_id)
383  vac%ang(nchi+1:n_theta,1) = chi_h(2:nchi) + pi
384  else
385  vac%x_vec(:,1) = r_h(:,r_id)
386  vac%x_vec(:,2) = z_h(:,r_id)
387  vac%ang(:,1) = chi_h
388  end if
389  !call print_ex_2D(['R','Z'],'',vac%x_vec,&
390  !&x=vac%ang(:,1:1),persistent=.true.)
391  !call print_ex_2D('cs','',vac%x_vec(:,2),&
392  !&x=vac%x_vec(:,1),persistent=.true.)
393 
394  ! calculate Cylindrical components of - J nabla psi
395  ! = - R (Z_theta nabla R - R_theta nabla Z)
396  ! and its poloidal derivative
397  vac%norm = 0._dp
398  vac%dnorm = 0._dp
399  ! nufft of Z
400  ierr = nufft(vac%ang(1:n_theta-1,1),&
401  &vac%x_vec(1:n_theta-1,2),x_vec_f)
402  chckerr('')
403  ! dZ/dtheta
404  do id = 1,size(x_vec_f,1)
405  vac%norm(1:n_theta-1,1) = &
406  &vac%norm(1:n_theta-1,1) + (id-1) * &
407  &(-x_vec_f(id,1)*sin((id-1)*&
408  &vac%ang(1:n_theta-1,1))&
409  &+x_vec_f(id,2)*cos((id-1)*&
410  &vac%ang(1:n_theta-1,1)))
411  end do
412  ! nufft of -R
413  ierr = nufft(vac%ang(1:n_theta-1,1),&
414  &-vac%x_vec(1:n_theta-1,1),x_vec_f)
415  chckerr('')
416  ! d(-R)/dtheta
417  do id = 1,size(x_vec_f,1)
418  vac%norm(1:n_theta-1,2) = &
419  &vac%norm(1:n_theta-1,2) + (id-1) * &
420  &(-x_vec_f(id,1)*sin((id-1)*&
421  &vac%ang(1:n_theta-1,1))&
422  &+x_vec_f(id,2)*cos((id-1)*&
423  &vac%ang(1:n_theta-1,1)))
424  end do
425  vac%norm(n_theta,:) = vac%norm(1,:)
426  ! multiply by -R
427  do id = 1,2
428  vac%norm(:,id) = - vac%x_vec(:,1)*vac%norm(:,id)
429  end do
430  ! nufft of norm
431  do jd = 1,2
432  ierr = nufft(vac%ang(1:n_theta-1,1),&
433  &vac%norm(1:n_theta-1,jd),x_vec_f)
434  chckerr('')
435  do id = 1,size(x_vec_f,1)
436  vac%dnorm(1:n_theta-1,jd) = &
437  &vac%dnorm(1:n_theta-1,jd) + (id-1) * &
438  &(-x_vec_f(id,1)*sin((id-1)*&
439  &vac%ang(1:n_theta-1,1))&
440  &+x_vec_f(id,2)*cos((id-1)*&
441  &vac%ang(1:n_theta-1,1)))
442  end do
443  end do
444  vac%dnorm(n_theta,:) = vac%dnorm(1,:)
445  end if
446 
447  ! broadcast to other processes
448  do id = 1,2
449  ierr = broadcast_var(vac%x_vec(:,id),n_procs-1)
450  chckerr('')
451  ierr = broadcast_var(vac%norm(:,id),n_procs-1)
452  chckerr('')
453  ierr = broadcast_var(vac%dnorm(:,id),n_procs-1)
454  chckerr('')
455  end do
456  ierr = broadcast_var(vac%ang(:,1),n_procs-1)
457  chckerr('')
458  else ! restarting
459  ! reconstruct old vacuum
460  ierr = reconstruct_pb3d_vac(vac,'vac')
461  chckerr('')
462  end if
463  end if
464 
465  call lvl_ud(-1)
466 
467  call writo('Done storing vacuum quantities')
468  end function store_vac_hel
469  end function store_vac
470 
499  integer function calc_gh(vac,n_r_in,lims_r_in,x_vec_in,G_in,H_in) &
500  &result(ierr)
501 
502  use num_vars, only: rank
503  use vac_vars, only: in_context
504  use vac_utilities, only: vec_dis2loc
505 #if ldebug
506  use num_vars, only: n_procs
507  use vac_utilities, only: mat_dis2loc
508 #endif
509 
510  character(*), parameter :: rout_name = 'calc_GH'
511 
512  ! input / output
513  type(vac_type), intent(inout), target :: vac
514  integer, intent(in), optional :: n_r_in
515  integer, intent(in), optional, target :: lims_r_in(:,:)
516  real(dp), intent(in), optional, target :: x_vec_in(:,:)
517  real(dp), intent(in), optional, target :: g_in(:,:)
518  real(dp), intent(in), optional, target :: h_in(:,:)
519 
520  ! local variables
521  integer :: n_r ! number of rows
522  integer :: id, jd ! counters
523  integer, pointer :: lims_r(:,:) ! row limits
524  real(dp), pointer :: g(:,:) ! G
525  real(dp), pointer :: h(:,:) ! H
526  real(dp), pointer :: x_vec(:,:) ! x_vec used in row
527  character(len=max_str_ln) :: err_msg ! error message
528  logical :: ext_in ! external influence point
529 #if ldebug
530  integer :: rd ! global counter for row
531  integer :: rdl ! local counter for row
532  integer :: i_rd ! index of subrow
533  integer :: desc_res(blacsctxtsize) ! descriptor for the result
534  real(dp), allocatable :: loc_ser(:,:) ! dummy serial variable
535  real(dp), allocatable :: r_sph(:) ! spherical r with r^2 = R^2+Z^2
536  real(dp), allocatable :: phi(:,:) ! test phi
537  real(dp), allocatable :: dphi(:,:) ! test dphi/dn
538  real(dp), allocatable :: res(:,:,:) ! H phi, G dphi
539  character(len=max_str_ln) :: plot_title(3) ! plot title
540  character(len=max_str_ln) :: plot_name ! plot name
541 #endif
542 
543  ! initialize ierr
544  ierr = 0
545 
546  ! user output
547  call writo('Calculate the G and H matrices')
548  call lvl_ud(1)
549 
550  ! test
551  if (present(n_r_in) .and. .not.present(x_vec_in) .or. &
552  &present(n_r_in) .and. .not.present(lims_r_in) .or. &
553  &present(n_r_in) .and. .not.present(g_in) .or. &
554  &present(n_r_in) .and. .not.present(h_in)) then
555  ierr = 1
556  chckerr('Need all optional variables')
557  end if
558 
559  ! set limits and variables
560  if (present(lims_r_in)) then
561  n_r = n_r_in
562  lims_r => lims_r_in
563  x_vec => x_vec_in
564  g => g_in
565  h => h_in
566  ext_in = .true.
567  else
568  n_r = vac%n_bnd
569  allocate(lims_r(2,size(vac%lims_r,2)))
570  allocate(x_vec(vac%n_bnd,size(vac%x_vec,2)))
571  lims_r = vac%lims_r
572  x_vec = vac%x_vec
573  g => vac%G
574  h => vac%H
575  ext_in = .false.
576  end if
577 
578  ! call specific procedure for vacuum style
579  select case (vac%style)
580  case (1) ! field-line 3-D
581  ierr = calc_gh_1(vac,n_r,lims_r,x_vec,g,h,ext_in)
582  chckerr('')
583  case (2) ! axisymmetric
584  ierr = calc_gh_2(vac,n_r,lims_r,x_vec,g,h,ext_in)
585  chckerr('')
586  case default
587  ierr = 1
588  err_msg = 'No vacuum style '//trim(i2str(vac%style))//&
589  &' possible'
590  chckerr(err_msg)
591  end select
592 
593  ! clean up
594  if (.not.present(lims_r_in)) then
595  deallocate(lims_r)
596  deallocate(x_vec)
597  end if
598  nullify(lims_r)
599  nullify(x_vec)
600  nullify(g)
601  nullify(h)
602 
603 #if ldebug
604  if (debug_calc_gh .and. in_context(vac%ctxt_HG)) then
605  ! initialize
606  call descinit(desc_res,vac%n_bnd,1,vac%bs,vac%bs,0,0,&
607  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
608  chckerr('descinit failed for res')
609 
610  ! plot H and G in HDF5
611  allocate(loc_ser(vac%n_bnd,vac%n_bnd))
612  ierr = mat_dis2loc(vac%ctxt_HG,vac%H,vac%lims_r,vac%lims_c,loc_ser,&
613  &proc=n_procs-1)
614  chckerr('')
615  if (rank.eq.n_procs-1) then
616  call plot_hdf5('H','H',reshape(loc_ser,[vac%n_bnd,vac%n_bnd,1]))
617  call plot_diff_hdf5(reshape(loc_ser,[vac%n_bnd,vac%n_bnd,1]),&
618  &reshape(transpose(loc_ser),[vac%n_bnd,vac%n_bnd,1]),&
619  &'H_transp')
620  end if
621  ierr = mat_dis2loc(vac%ctxt_HG,vac%G,vac%lims_r,vac%lims_c,loc_ser,&
622  &proc=n_procs-1)
623  chckerr('')
624  if (rank.eq.n_procs-1) then
625  call plot_hdf5('G','G',reshape(loc_ser,[vac%n_bnd,vac%n_bnd,1]))
626  call plot_diff_hdf5(reshape(loc_ser,[vac%n_bnd,vac%n_bnd,1]),&
627  &reshape(transpose(loc_ser),[vac%n_bnd,vac%n_bnd,1]),&
628  &'G_transp')
629  end if
630  deallocate(loc_ser)
631 
632  ! the rest is only for axisymmetric vacua.
633  if (vac%style.eq.2) then
634  ! test whether G and H hold for test potentials phi
635  if (vac%lims_c(1,1).eq.1) then ! this process owns (part of) first column
636  ! set up variables
637  allocate(res(vac%n_loc(1),3,2))
638  allocate(r_sph(vac%n_loc(1)))
639  allocate(phi(vac%n_loc(1),2))
640  allocate(dphi(vac%n_loc(1),2))
641 
642  ! set up rhs for H (phi) and G (dphi)
643  subrows: do i_rd = 1,size(vac%lims_r,2)
644  row: do rd = vac%lims_r(1,i_rd),vac%lims_r(2,i_rd)
645  ! set local row index
646  rdl = sum(vac%lims_r(2,1:i_rd-1)-&
647  &vac%lims_r(1,1:i_rd-1)+1) + &
648  &rd-vac%lims_r(1,i_rd)+1
649 
650  ! spherical radius
651  r_sph(rdl) = sqrt(sum(vac%x_vec(rd,:)**2))
652 
653  ! test 1
654  phi(rdl,1) = vac%x_vec(rd,1)**vac%prim_X
655  dphi(rdl,1) = - vac%prim_X*vac%norm(rd,1)/&
656  &vac%x_vec(rd,1) * phi(rdl,1)
657 
658  ! test 2
659  phi(rdl,2) = (vac%x_vec(rd,1)/&
660  &(r_sph(rdl)+abs(vac%x_vec(rd,2))))**vac%prim_X
661  dphi(rdl,2) = vac%prim_X*(vac%norm(rd,2)-&
662  &vac%norm(rd,1)*vac%x_vec(rd,2)/&
663  &vac%x_vec(rd,1)) / &
664  &r_sph(rdl)*sign(1._dp,vac%x_vec(rd,2)) * &
665  &phi(rdl,2)
666  end do row
667  end do subrows
668  else
669  allocate(res(0,3,2))
670  allocate(r_sph(0))
671  allocate(phi(0,2))
672  allocate(dphi(0,2))
673  end if
674 
675  ! output
676  allocate(loc_ser(vac%n_bnd,3))
677  ierr = vec_dis2loc(vac%ctxt_HG,r_sph,vac%lims_r,loc_ser(:,1),&
678  &proc=n_procs-1)
679  chckerr('')
680  if (rank.eq.n_procs-1) then
681  plot_title(1) = 'spherical r'
682  plot_name = 'r_sph'
683  call print_ex_2d(plot_title(1),plot_name,loc_ser(:,1),&
684  &x=vac%ang(:,1),draw=.false.)
685  call draw_ex([plot_title],plot_name,1,1,.false.)
686  end if
687  do jd = 1,size(phi,2)
688  ierr = vec_dis2loc(vac%ctxt_HG,phi(:,jd),vac%lims_r,&
689  &loc_ser(:,jd),proc=n_procs-1)
690  chckerr('')
691  end do
692  if (rank.eq.n_procs-1) then
693  plot_title(1) = 'test potential Ï•'
694  plot_name = 'phi'
695  call print_ex_2d([plot_title(1)],plot_name,&
696  &loc_ser(:,1:size(phi,2)),x=vac%ang(:,1:1),draw=.false.)
697  call draw_ex([plot_title(1)],plot_name,size(phi,2),1,&
698  &.false.)
699  end if
700  do jd = 1,size(dphi,2)
701  ierr = vec_dis2loc(vac%ctxt_HG,dphi(:,jd),vac%lims_r,&
702  &loc_ser(:,jd),proc=n_procs-1)
703  chckerr('')
704  end do
705  if (rank.eq.n_procs-1) then
706  plot_title(1) = 'normal derivative of test potential dÏ•/dn'
707  plot_name = 'dphi'
708  call print_ex_2d([plot_title(1)],plot_name,&
709  &loc_ser(:,1:size(dphi,2)),x=vac%ang(:,1:1),&
710  &draw=.false.)
711  call draw_ex([plot_title(1)],plot_name,size(dphi,2),1,&
712  &.false.)
713  end if
714  deallocate(loc_ser)
715 
716  ! loop over test potentials
717  do jd = 1,size(phi,2)
718  ! calculate H phi
719  call pdgemv('N',vac%n_bnd,vac%n_bnd,1._dp,vac%H,1,1,&
720  &vac%desc_H,phi(:,jd),1,1,desc_res,1,0._dp,res(:,1,jd),&
721  &1,1,desc_res,1)
722 
723  ! calculate - G dphi
724  call pdgemv('N',vac%n_bnd,vac%n_bnd,-1._dp,vac%G,1,1,&
725  &vac%desc_G,dphi(:,jd),1,1,desc_res,1,0._dp,&
726  &res(:,2,jd),1,1,desc_res,1)
727 
728  ! add them together
729  res(:,3,jd) = sum(res(:,1:2,jd),2)
730  end do
731 
732  ! output of tests
733  allocate(loc_ser(vac%n_bnd,3))
734  do id = 1,size(phi,2)
735  do jd = 1,3
736  ierr = vec_dis2loc(vac%ctxt_HG,res(:,jd,id),vac%lims_r,&
737  &loc_ser(:,jd),proc=n_procs-1)
738  chckerr('')
739  end do
740  if (rank.eq.n_procs-1) then
741  select case (id)
742  case (1)
743  plot_title(1) = 'H R^n'
744  plot_title(2) = '-G n Z_θ R^n'
745  plot_title(3) = '(H - G n Z_θ) R^n'
746  plot_name = 'test_vac_1'
747  case (2)
748  plot_title(1) = 'H (R/(r+|Z|))^n'
749  plot_title(2) = &
750  &'- G |Z|/Z n dr/dθ (R/(r+|Z|))^n'
751  plot_title(3) = '(H - G |Z|/Z n dr/dθ) &
752  &(R/(r+|Z|))^n'
753  plot_name = 'test_vac_2'
754  end select
755  call print_ex_2d(plot_title,trim(plot_name)//'_all',&
756  &loc_ser(:,:),x=vac%ang(:,1:1),&
757  &draw=.false.)
758  call draw_ex(plot_title,trim(plot_name)//'_all',3,1,&
759  &.false.)
760  call print_ex_2d(plot_title(3),plot_name,&
761  &loc_ser(:,3),x=vac%ang(:,1),&
762  &draw=.false.)
763  call draw_ex(plot_title(3:3),plot_name,1,1,.false.)
764  end if
765  end do
766  deallocate(loc_ser)
767  end if
768  end if
769 #endif
770 
771  call lvl_ud(-1)
772  end function calc_gh
773 
777  integer function calc_gh_1(vac,n_r_in,lims_r_in,x_vec_in,G_in,H_in,ext_in) &
778  &result(ierr)
779 
780  use vac_utilities, only: calc_gh_int_1
781  use grid_vars, only: min_par_x, max_par_x, min_alpha, max_alpha, n_alpha
782  use rich_vars, only: n_par_x
783  use num_vars, only: tol_zero
784  use vac_vars, only: in_context
785  use vac_utilities, only: vec_dis2loc
786 
787  character(*), parameter :: rout_name = 'calc_GH_1'
788 
789  ! input / output
790  type(vac_type), intent(inout) :: vac
791  integer, intent(in) :: n_r_in
792  integer, intent(in) :: lims_r_in(:,:)
793  real(dp), intent(in) :: x_vec_in(:,:)
794  real(dp), intent(in), pointer :: g_in(:,:)
795  real(dp), intent(in), pointer :: h_in(:,:)
796  logical, intent(in) :: ext_in
797 
798  ! local variables
799  integer :: kd ! counter
800  integer :: i_rd, i_cd ! index of subrow and column
801  integer :: rd, cd ! global counters for row and column
802  integer :: rdl, cdl ! local counters for row and column, used for G and H
803  integer :: lims_c_loc(2) ! local column limits
804  integer :: desc_res(blacsctxtsize) ! descriptor for the result
805  real(dp) :: dpar_x ! step size in par_X
806  real(dp) :: dalpha ! step size in alpha
807  real(dp) :: tol_loc ! local tolerance on rho for singular points
808  real(dp) :: g_loc(2) ! local G
809  real(dp) :: h_loc(2) ! local H
810  real(dp), allocatable :: res(:,:) ! sums of rows of H and vector of ones
811  real(dp), allocatable :: loc_res(:) ! local res
812 
813  ! initialize ierr
814  ierr = 0
815 
816  call writo('Using field-line 3-D Boundary Element Method')
817 
818  ! initialize variables
819  dpar_x = (max_par_x-min_par_x)/(n_par_x-1)*pi
820  if (n_alpha.gt.1) then
821  dalpha = (max_alpha-min_alpha)/(n_alpha-1)*pi
822  else
823  dalpha = 0._dp
824  end if
825 
826  ! set up local tolerance
827  tol_loc = tol_zero*10 ! tolerance for singular interval
828 
829  ! iterate over all subrows: position of singular point
830  subrows: do i_rd = 1,size(lims_r_in,2)
831  ! iterate over all rows of this subrow
832  row: do rd = lims_r_in(1,i_rd),lims_r_in(2,i_rd)
833  ! set local row index
834  rdl = sum(lims_r_in(2,1:i_rd-1)-lims_r_in(1,1:i_rd-1)+1) + &
835  &rd-lims_r_in(1,i_rd)+1
836 
837  ! iterate over all subcolumns: subintegrals
838  subcols: do i_cd = 1,size(vac%lims_c,2)
839  ! set up local limits, including possible overlap
840  lims_c_loc(1) = max(1,vac%lims_c(1,i_cd)-1)
841  lims_c_loc(2) = min(vac%n_bnd,vac%lims_c(2,i_cd)+1)
842 
843  ! iterate over all pairs of columns of this subcolumn
844  ! (upper limit is one lower as lims_c_loc
845  col: do cd = lims_c_loc(1),lims_c_loc(2)-1
846  ! set local column index
847  cdl = sum(vac%lims_c(2,1:i_cd-1)-&
848  &vac%lims_c(1,1:i_cd-1)+1) + &
849  &cd-vac%lims_c(1,i_cd)+1
850 
851  ! calculate contributions to H and G
852  call calc_gh_int_1(g_loc,h_loc,&
853  &vac%x_vec(cd:cd+1,:),x_vec_in(rd,:),&
854  &vac%norm(cd:cd+1,:),vac%h_fac(rd,:),&
855  &[dpar_x,dalpha],tol_loc**2)
856 
857  ! loop over left and right side of interval
858  do kd = 0,1
859  if (cd+kd.ge.vac%lims_c(1,i_cd) .and. & ! lower local bound
860  &cd+kd.le.vac%lims_c(2,i_cd)) then ! upper local bound
861  ! cycle if first or last global point
862  if (.not.on_field_line(cd,kd,vac%n_ang(1))) &
863  &cycle
864  g_in(rdl,cdl+kd) = g_in(rdl,cdl+kd) + &
865  &g_loc(1+kd)
866  h_in(rdl,cdl+kd) = h_in(rdl,cdl+kd) + &
867  &h_loc(1+kd)
868  end if
869  end do
870  end do col
871  end do subcols
872  end do row
873  end do subrows
874 
875  ! if in context and not external, indirctly calculate contribution due
876  ! to fundament solution through constant potential
877  if (in_context(vac%ctxt_HG) .and. .not.ext_in) then
878  if (vac%lims_c(1,1).eq.1) then ! this process owns (part of) first column
879  allocate(res(vac%n_loc(1),2))
880  else
881  allocate(res(0,2))
882  end if
883  call descinit(desc_res,vac%n_bnd,1,vac%bs,vac%bs,0,0,&
884  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
885  chckerr('descinit failed for res')
886 
887  ! set diagonal of H to zero
888  ! iterate over all subrows: position of singular point
889  subrows2: do i_rd = 1,size(lims_r_in,2)
890  ! iterate over all rows of this subrow
891  row2: do rd = lims_r_in(1,i_rd),lims_r_in(2,i_rd)
892  ! set local row index
893  rdl = sum(lims_r_in(2,1:i_rd-1)-&
894  &lims_r_in(1,1:i_rd-1)+1) + &
895  &rd-lims_r_in(1,i_rd)+1
896 
897  ! iterate over all subcolumns: subintegrals
898  subcols2: do i_cd = 1,size(vac%lims_c,2)
899  cd = rd
900  if (vac%lims_c(1,i_cd).le.cd .and. &
901  &cd.le.vac%lims_c(2,i_cd)) then ! this subcolumn includes the diagonal
902  ! set local column index
903  cdl = sum(vac%lims_c(2,1:i_cd-1)-&
904  &vac%lims_c(1,1:i_cd-1)+1) + &
905  &cd-vac%lims_c(1,i_cd)+1
906 
907  ! set diagonal of H to zero
908  h_in(rdl,cdl) = 0._dp
909  end if
910  end do subcols2
911  end do row2
912  end do subrows2
913 
914  ! set res(1) to ones
915  res(:,1) = 1._dp
916 
917  ! calculate sum of rows of H
918  call pdgemv('N',vac%n_bnd,vac%n_bnd,1._dp,vac%H,1,1,&
919  &vac%desc_H,res(:,1),1,1,desc_res,1,0._dp,res(:,2),&
920  &1,1,desc_res,1)
921 
922  ! add 4pi
923  res(:,2) = res(:,2) + 4*pi
924 
925  ! distribute to all processes
926  allocate(loc_res(vac%n_bnd))
927  ierr = vec_dis2loc(vac%ctxt_HG,res(:,2),vac%lims_r,loc_res,&
928  &proc=-1)
929  chckerr('')
930 
931  ! Put in diagonal of H
932  ! iterate over all subrows: position of singular point
933  subrows3: do i_rd = 1,size(lims_r_in,2)
934  ! iterate over all rows of this subrow
935  row3: do rd = lims_r_in(1,i_rd),lims_r_in(2,i_rd)
936  ! set local row index
937  rdl = sum(lims_r_in(2,1:i_rd-1)-&
938  &lims_r_in(1,1:i_rd-1)+1) + &
939  &rd-lims_r_in(1,i_rd)+1
940 
941  ! iterate over all subcolumns: subintegrals
942  subcols3: do i_cd = 1,size(vac%lims_c,2)
943  cd = rd
944  if (vac%lims_c(1,i_cd).le.cd .and. &
945  &cd.le.vac%lims_c(2,i_cd)) then ! this subcolumn includes the diagonal
946  ! set local column index
947  cdl = sum(vac%lims_c(2,1:i_cd-1)-&
948  &vac%lims_c(1,1:i_cd-1)+1) + &
949  &cd-vac%lims_c(1,i_cd)+1
950 
951  ! set diagonal of H to inverse of sum
952  h_in(rdl,cdl) = - loc_res(rd)
953  end if
954  end do subcols3
955  end do row3
956  end do subrows3
957 
958  ! clean up
959  deallocate(res)
960  deallocate(loc_res)
961  end if
962  contains
963  ! checks whether an interval is on the field line
964  logical function on_field_line(cd,kd,n_ang) result(res)
965  ! input / output
966  integer, intent(in) :: cd ! position of interval on field line
967  integer, intent(in) :: kd ! left or right point of interval
968  integer, intent(in) :: n_ang ! number of points on field line
969 
970  res = .true.
971  if (mod(cd,n_ang).eq.1 .and. kd.eq.0) res = .false. ! start of field line
972  if (mod(cd,n_ang).eq.0 .and. kd.eq.1) res = .false. ! end of field line
973  end function on_field_line
974  end function calc_gh_1
975 
981  integer function calc_gh_2(vac,n_r_in,lims_r_in,x_vec_in,G_in,H_in,ext_in) &
982  &result(ierr)
983 
984  use dtorh, only: dtorh1
985  use vac_utilities, only: calc_gh_int_2
986  use mpi_utilities, only: get_ser_var
987  use num_ops, only: calc_zero_hh, calc_zero_zhang
988  use num_vars, only: rank
989 
990  character(*), parameter :: rout_name = 'calc_GH_2'
991 
992  ! input / output
993  type(vac_type), intent(inout) :: vac
994  integer, intent(in) :: n_r_in
995  integer, intent(in) :: lims_r_in(:,:)
996  real(dp), intent(in) :: x_vec_in(:,:)
997  real(dp), intent(in), pointer :: g_in(:,:)
998  real(dp), intent(in), pointer :: h_in(:,:)
999  logical, intent(in) :: ext_in
1000 
1001  ! local variables
1002  integer :: kd ! counter
1003  integer :: i_rd, i_cd ! index of subrow and column
1004  integer :: rd, cd ! global counters for row and column
1005  integer :: rdl, cdl ! local counters for row and column, used for G and H
1006  integer :: max_n ! maximum degree
1007  integer :: lims_c_loc(2) ! local column limits
1008  real(dp) :: b_coeff(2) ! sum_k=1^n 2/2k-1 for n and n-1
1009  real(dp) :: tol_loc ! local tolerance on rho for singular points
1010  real(dp) :: del_r(2,2) ! unit vectors of next and previous subinterval
1011  real(dp) :: g_loc(2) ! local G
1012  real(dp) :: h_loc(2) ! local H
1013  real(dp) :: beta_comp ! complementary beta
1014  real(dp) :: ang ! local angle
1015  real(dp), parameter :: tol_q = 1.e-6_dp ! tolerance on Q approximation
1016  real(dp), allocatable :: beta(:) ! beta
1017  real(dp), allocatable :: loc_ser(:) ! dummy serial variable
1018  real(dp), allocatable :: rho2(:) ! square of distance in projected poloidal plane for one column
1019  real(dp), allocatable :: arg(:) ! argument 1 + rho^2/(RiRj)(
1020  real(dp), allocatable :: aij(:) ! Aij helper variable for H
1021  real(dp), allocatable :: ql(:,:,:) ! ql for all rho2's
1022  real(dp), allocatable :: pl_loc(:), ql_loc(:) ! local toroidal functions
1023  logical, allocatable :: was_calc(:,:) ! was already calculated
1024  logical, allocatable :: sing_col(:) ! singular column, used only with external influence
1025  logical :: was_calc_sym ! symmetrical was_calc
1026  character(len=max_str_ln) :: err_msg ! error message
1027 #if ldebug
1028  real(dp) :: rho2_lims(2) ! limits on rho2
1029  real(dp) :: gamma_lims(2) ! limits on gamma
1030  integer :: beta_lims(2) ! limits to plot beta
1031  !integer :: lims_cl(2) ! vac%lims_c in local coordinates
1032  !character(len=max_str_ln) :: plot_title ! plot title
1033  !character(len=max_str_ln) :: plot_name ! plot name
1034 #endif
1035 
1036  ! initialize ierr
1037  ierr = 0
1038 
1039  call writo('Using axisymmetric Boundary Element Method')
1040 
1041  ! set up b_coeff and local tolerance
1042  b_coeff = 0._dp
1043  do kd = 1,vac%prim_X
1044  b_coeff(2) = b_coeff(2) + 2._dp/(2*kd-1)
1045  end do
1046  b_coeff(1) = b_coeff(2) - 2._dp/(2*vac%prim_X-1)
1047  err_msg = calc_zero_zhang(tol_loc,fun_tol,&
1048  &[1.e-10_dp,& ! very low lower limit
1049  &exp(-2*b_coeff(2))]) ! upper limit for -1/2ln(x) = b
1050  if (trim(err_msg).ne.'') then
1051  ierr = 1
1052  chckerr(trim(err_msg))
1053  end if
1054  tol_loc = 8*minval(vac%x_vec(:,1))*tol_loc ! tolerance on rho
1055  ierr = get_ser_var([tol_loc],loc_ser,scatter=.true.)
1056  chckerr('')
1057  tol_loc = maxval(loc_ser)
1058  deallocate(loc_ser)
1059  call writo('Tolerance on rho for analytical integral: '//&
1060  &trim(r2str(tol_loc)))
1061 
1062  ! allocate helper variables
1063  allocate(was_calc(1:n_r_in,1:vac%n_bnd))
1064  was_calc = .false. ! will remain false if external influence
1065  allocate(ql(1:n_r_in,1:vac%n_bnd,2))
1066  allocate(pl_loc(0:max(vac%prim_X,1)))
1067  allocate(ql_loc(0:max(vac%prim_X,1)))
1068 #if ldebug
1069  rho2_lims = [huge(1._dp),-huge(1._dp)]
1070  gamma_lims = [huge(1._dp),-huge(1._dp)]
1071 #endif
1072 
1073  ! iterate over all subrows: position of singular point
1074  subrows: do i_rd = 1,size(lims_r_in,2)
1075  ! initialize helper variables
1076  allocate(beta(lims_r_in(1,i_rd):lims_r_in(2,i_rd))) ! global index
1077 #if ldebug
1078  beta_lims = [lims_r_in(2,i_rd)+1,lims_r_in(1,i_rd)-1]
1079 #endif
1080 
1081  ! iterate over all rows of this subrow
1082  row: do rd = lims_r_in(1,i_rd),lims_r_in(2,i_rd)
1083  ! set local row index
1084  rdl = sum(lims_r_in(2,1:i_rd-1)-lims_r_in(1,1:i_rd-1)+1) + &
1085  &rd-lims_r_in(1,i_rd)+1
1086 
1087  ! iterate over all subcolumns: subintegrals
1088  subcols: do i_cd = 1,size(vac%lims_c,2)
1089  ! set up local limits, including possible overlap
1090  lims_c_loc(1) = max(1,vac%lims_c(1,i_cd)-1)
1091  lims_c_loc(2) = min(vac%n_bnd,vac%lims_c(2,i_cd)+1)
1092 
1093  ! initialize helper variables
1094  allocate(rho2(lims_c_loc(1):lims_c_loc(2))) ! global index
1095  allocate(arg(lims_c_loc(1):lims_c_loc(2))) ! global index
1096  allocate(aij(lims_c_loc(1):lims_c_loc(2))) ! global index
1097 
1098  ! initialize singular column flags
1099  if (ext_in) then
1100  if (allocated(sing_col)) deallocate(sing_col)
1101  allocate(sing_col(&
1102  &vac%lims_c(1,i_cd):vac%lims_c(2,i_cd)))
1103  sing_col = .false.
1104  end if
1105 
1106  ! iterate over all columns of this subcolumn, possibly with
1107  ! an overlap left and/or right to set up the Aij
1108  col: do cd = lims_c_loc(1),lims_c_loc(2)
1109  ! check rho2 and set arg to calculate tor. functions
1110  rho2(cd) = sum((vac%x_vec(cd,:)-x_vec_in(rd,:))**2)
1111  arg(cd) = 1._dp + &
1112  &rho2(cd)/(2*vac%x_vec(cd,1)*x_vec_in(rd,1))
1113  was_calc_sym = .false.
1114  if (.not.ext_in) was_calc_sym = was_calc(cd,rd)
1115  if (rho2(cd).le.tol_loc**2) then
1116  ! don't set it
1117  else if (was_calc_sym) then ! check symmetric element
1118  ql(rd,cd,:) = ql(cd,rd,:) ! set it
1119  else ! not yet calculated by this process
1120 #if ldebug
1121  ! update extrema
1122  if (rho2(cd).lt.rho2_lims(1)) &
1123  &rho2_lims(1) = rho2(cd)
1124  if (rho2(cd).gt.rho2_lims(2)) &
1125  &rho2_lims(2) = rho2(cd)
1126  if (arg(cd).lt.gamma_lims(1)) &
1127  &gamma_lims(1) = arg(cd)
1128  if (arg(cd).gt.gamma_lims(2)) &
1129  &gamma_lims(2) = arg(cd)
1130 #endif
1131  err_msg = 'The solutions did not converge for &
1132  &rho^2 = '//trim(r2str(rho2(cd)))//&
1133  &' and n = '//trim(i2str(vac%prim_X))
1134  if (vac%prim_X.gt.0) then ! positive n
1135  ierr = dtorh1(arg(cd),0,vac%prim_X,pl_loc,&
1136  &ql_loc,max_n)
1137  chckerr('')
1138  if (max_n.lt.vac%prim_X) then
1139  ierr = 1
1140  chckerr(err_msg)
1141  end if
1142  ql(rd,cd,:) = ql_loc(vac%prim_X-1:vac%prim_X)
1143  else if (vac%prim_X.eq.0) then
1144  ierr = dtorh1(arg(cd),0,1,pl_loc,ql_loc,max_n)
1145  chckerr('')
1146  if (max_n.lt.vac%prim_X) then
1147  ierr = 1
1148  chckerr(err_msg)
1149  end if
1150  ql(rd,cd,1) = ql_loc(1) ! Q_-3/2 = Q_1/2
1151  ql(rd,cd,2) = ql_loc(0) ! Q_-1/2
1152  else ! negative n
1153  ierr = 1
1154  err_msg = 'negative n not yet implemented' ! SHOULD JUST BE Q_{n-1/2} = Q_{-n-1/2}
1155  chckerr('')
1156  end if
1157  if (.not.ext_in) was_calc(rd,cd) = .true.
1158  end if
1159 
1160  ! set up Aij
1161  if (rho2(cd).le.tol_loc**2) then
1162  if (ext_in) then
1163  if (cd.ge.vac%lims_c(1,i_cd) .and. &
1164  &cd.le.vac%lims_c(2,i_cd)) &
1165  &sing_col(cd) = .true.
1166  !!Aij(cd) = 0._dp
1167  else
1168  aij(cd) = 0.5_dp*(vac%prim_X - 0.5_dp) * &
1169  &(x_vec_in(rd,1) * &
1170  &(vac%norm(rd,1)*vac%dnorm(rd,2)-&
1171  &vac%norm(rd,2)*vac%dnorm(rd,1))/&
1172  &sum(vac%norm(rd,:)**2) &
1173  &+vac%norm(rd,1)/x_vec_in(rd,1))
1174  end if
1175  else
1176  aij(cd) = 2._dp*x_vec_in(rd,1)*vac%x_vec(cd,1) / &
1177  &(4._dp*x_vec_in(rd,1)*vac%x_vec(cd,1)+&
1178  &rho2(cd)) *(vac%prim_X - 0.5_dp) * &
1179  &(- 2._dp/rho2(cd) * sum(vac%norm(cd,:)*&
1180  &(vac%x_vec(cd,:)-x_vec_in(rd,:))) + &
1181  &vac%norm(cd,1)/vac%x_vec(cd,1))
1182  end if
1183  end do col
1184 
1185 #if ldebug
1186  if (debug_calc_gh) then
1187  !! output
1188  !plot_title = 'for row '//trim(i2str(rd))//&
1189  !&' and starting column '//trim(i2str(lims_c_loc(1)))
1190  !plot_name = '_'//trim(i2str(rd))//'_'//&
1191  !&trim(i2str(lims_c_loc(1)))
1192  !call print_ex_2D('rho^2 '//trim(plot_title),&
1193  !&'rho2'//trim(plot_name),rho2,&
1194  !&x=vac%ang(lims_c_loc(1):lims_c_loc(2),1),&
1195  !&draw=.false.)
1196  !call draw_ex(['rho^2 '//trim(plot_title)],&
1197  !&'rho2'//trim(plot_name),1,1,.false.)
1198  !call print_ex_2D('γ '//trim(plot_title),&
1199  !&'gamma'//trim(plot_name),arg,&
1200  !&x=vac%ang(lims_c_loc(1):lims_c_loc(2),1),&
1201  !&draw=.false.)
1202  !call draw_ex(['γ '//trim(plot_title)],&
1203  !&'gamma'//trim(plot_name),1,1,.false.)
1204  !call print_ex_2D('A_{ij} '//trim(plot_title),&
1205  !&'Aij'//trim(plot_name),arg,&
1206  !&x=vac%ang(lims_c_loc(1):lims_c_loc(2),1),&
1207  !&draw=.false.)
1208  !call draw_ex(['A_{ij} '//trim(plot_title)],&
1209  !&'Aij'//trim(plot_name),1,1,.false.)
1210  end if
1211 #endif
1212 
1213  ! iterate over all pairs of columns of this subcolumn
1214  ! (upper limit is one lower as previous)
1215  col2: do cd = lims_c_loc(1),lims_c_loc(2)-1
1216  ! set local column index
1217  cdl = sum(vac%lims_c(2,1:i_cd-1)-&
1218  &vac%lims_c(1,1:i_cd-1)+1) + &
1219  &cd-vac%lims_c(1,i_cd)+1
1220 
1221  ! check for singular columns if external influence
1222  if (ext_in) then
1223  if (any(sing_col)) then
1224  do kd = 0,1
1225  if (cd+kd.ge.vac%lims_c(1,i_cd) .and. &
1226  &cd+kd.le.vac%lims_c(2,i_cd)) then
1227  g_in(rdl,cdl+kd) = 0._dp
1228  h_in(rdl,cdl+kd) = 4*pi
1229  end if
1230  end do
1231  cycle col2 ! pass to next column
1232  end if
1233  end if
1234 
1235  ! set ang
1236  if (ext_in) then
1237  ang = 0._dp
1238  else
1239  ang = vac%ang(rd,1)
1240  end if
1241 
1242  ! calculate contributions to H and G
1243  call calc_gh_int_2(g_loc,h_loc,&
1244  &vac%ang(cd:cd+1,1),ang,&
1245  &vac%x_vec(cd:cd+1,:),x_vec_in(rd,:),&
1246  &vac%norm(cd:cd+1,:),vac%norm(rd,:),&
1247  &aij(cd:cd+1),ql(rd,cd:cd+1,:),&
1248  &tol_loc**2,b_coeff,vac%prim_X)
1249 
1250  do kd = 0,1
1251  if (cd+kd.ge.vac%lims_c(1,i_cd) .and. &
1252  &cd+kd.le.vac%lims_c(2,i_cd)) then
1253  g_in(rdl,cdl+kd) = g_in(rdl,cdl+kd) + &
1254  &g_loc(1+kd)
1255  h_in(rdl,cdl+kd) = h_in(rdl,cdl+kd) + &
1256  &h_loc(1+kd)
1257  end if
1258  end do
1259  end do col2
1260 
1261 #if ldebug
1262  if (debug_calc_gh) then
1263  !! local limits
1264  !lims_cl = sum(vac%lims_c(2,1:i_cd-1)-&
1265  !&vac%lims_c(1,1:i_cd-1)+1) + &
1266  !&vac%lims_c(:,i_cd)-vac%lims_c(1,i_cd)+1
1267  !! output
1268  !call print_ex_2D('G '//trim(plot_title),&
1269  !&'G'//trim(plot_name),&
1270  !&G_in(rdl,lims_cl(1):lims_cl(2)),x&
1271  !&=vac%ang(vac%lims_c(1,i_cd):vac%lims_c(2,i_cd),1),&
1272  !&draw=.false.)
1273  !call draw_ex(['G '//trim(plot_title)],&
1274  !&'G'//trim(plot_name),1,1,.false.)
1275  !call print_ex_2D('H '//trim(plot_title),&
1276  !&'H'//trim(plot_name),&
1277  !&H_in(rdl,lims_cl(1):lims_cl(2)),x&
1278  !&=vac%ang(vac%lims_c(1,i_cd):vac%lims_c(2,i_cd),1),&
1279  !&draw=.false.)
1280  !call draw_ex(['H '//trim(plot_title)],&
1281  !&'H'//trim(plot_name),1,1,.false.)
1282  end if
1283 #endif
1284 
1285  ! calculate contribution due to fundamental solution
1286  ! (only if not external)
1287  if (.not.ext_in .and. vac%lims_c(1,i_cd).le.rd .and. &
1288  &rd.le.vac%lims_c(2,i_cd)) then ! this subcolumn includes the diagonal
1289  ! set local column index
1290  cdl = sum(vac%lims_c(2,1:i_cd-1)-&
1291  &vac%lims_c(1,1:i_cd-1)+1) + &
1292  &rd-vac%lims_c(1,i_cd)+1
1293 
1294  if (rd.eq.1 .or. rd.eq.vac%n_bnd) then ! first or last global point
1295  ! loop around to get previous direction vector
1296  del_r(1,:) = vac%x_vec(vac%n_bnd,:)-&
1297  &vac%x_vec(vac%n_bnd-1,:)
1298  del_r(2,:) = vac%x_vec(2,:)-vac%x_vec(1,:)
1299  else ! all points in between
1300  del_r(1,:) = vac%x_vec(rd,:)-vac%x_vec(rd-1,:)
1301  del_r(2,:) = vac%x_vec(rd+1,:)-vac%x_vec(rd,:)
1302  end if
1303  do kd = 1,2
1304  del_r(kd,:) = del_r(kd,:)/sqrt(sum(del_r(kd,:)**2))
1305  end do
1306  beta_comp = acos(sum(product(del_r,1)))
1307  if (del_r(1,1)*del_r(2,2) .lt. &
1308  &del_r(1,2)*del_r(2,1)) beta_comp = - beta_comp ! concave
1309  beta(rd) = beta_comp + pi
1310 #if ldebug
1311  if (debug_calc_gh) then
1312  if (rd.lt.beta_lims(1)) beta_lims(1) = rd
1313  if (rd.gt.beta_lims(2)) beta_lims(2) = rd
1314  end if
1315 #endif
1316 
1317  h_in(rdl,cdl) = h_in(rdl,cdl) - 2*beta(rd)
1318  end if
1319 
1320  ! clean up
1321  deallocate(arg,rho2,aij)
1322  end do subcols
1323  end do row
1324 
1325 #if ldebug
1326  if (debug_calc_gh) then
1327  !! output
1328  !if (beta_lims(2).ge.beta_lims(1)) then
1329  !plot_title = 'beta/Ï€ for starting row '//&
1330  !&trim(i2str(beta_lims(1)))
1331  !plot_name = 'beta_'//trim(i2str(beta_lims(1)))
1332  !call print_ex_2D(plot_title,plot_name,&
1333  !&beta(beta_lims(1):beta_lims(2))/pi,&
1334  !&x=vac%ang(beta_lims(1):beta_lims(2),1),draw=.false.)
1335  !call draw_ex([plot_title],plot_name,1,1,.false.)
1336  !end if
1337  end if
1338 #endif
1339 
1340  ! clean up
1341  deallocate(beta)
1342  end do subrows
1343 
1344 #if ldebug
1345  if (debug_calc_gh) then
1346  ! output
1347  ierr = get_ser_var([rho2_lims(1)],loc_ser)
1348  chckerr('')
1349  if (rank.eq.0) rho2_lims(1) = minval(loc_ser)
1350  deallocate(loc_ser)
1351  ierr = get_ser_var([rho2_lims(2)],loc_ser)
1352  chckerr('')
1353  if (rank.eq.0) rho2_lims(2) = maxval(loc_ser)
1354  deallocate(loc_ser)
1355  ierr = get_ser_var([gamma_lims(1)],loc_ser)
1356  chckerr('')
1357  if (rank.eq.0) gamma_lims(1) = minval(loc_ser)
1358  deallocate(loc_ser)
1359  ierr = get_ser_var([gamma_lims(2)],loc_ser)
1360  chckerr('')
1361  if (rank.eq.0) gamma_lims(2) = maxval(loc_ser)
1362  deallocate(loc_ser)
1363  if (rank.eq.0) then
1364  call writo('limits in rho: '//trim(r2str(sqrt(rho2_lims(1))))//&
1365  &'..'//trim(r2str(sqrt(rho2_lims(2)))))
1366  call writo('limits in gamma: '//trim(r2str(gamma_lims(1)))//&
1367  &'..'//trim(r2str(gamma_lims(2))))
1368  end if
1369  end if
1370 #endif
1371  contains
1372  ! function that returns f = -1/2 ln(x)-b - x/T.
1373  !
1374  ! It uses b_coeff (b) and tol_Q (T) from the parent function where T is
1375  ! a hard-coded positive value << 1.
1376  !
1378  function fun_tol(val)
1379  ! input / output
1380  real(dp), intent(in) :: val
1381  real(dp) :: fun_tol
1382 
1383  fun_tol = -0.5_dp*log(val)-b_coeff(2)-val/tol_q
1384  end function fun_tol
1385  end function calc_gh_2
1386 
1416  integer function calc_vac_res(mds,vac) result(ierr)
1417  use grid_vars, only: min_par_x, max_par_x, min_alpha, max_alpha, n_alpha
1418  use mpi
1419  use rich_vars, only: rich_lvl, n_par_x
1420  use num_vars, only: jump_to_sol, rich_restart_lvl, eq_style, &
1421  &use_pol_flux_f, rank, n_procs
1422  use pb3d_ops, only: reconstruct_pb3d_vac
1423  use x_vars, only: n_mod_x
1424  use vac_vars, only: set_loc_lims, in_context
1425  use vac_utilities, only: mat_dis2loc
1426  use eq_vars, only: vac_perm
1427  use num_utilities, only: lcm
1428  use mpi_utilities, only: broadcast_var
1429 
1430  character(*), parameter :: rout_name = 'calc_vac_res'
1431 
1432  ! input / output
1433  type(modes_type), intent(in) :: mds
1434  type(vac_type), intent(inout) :: vac
1435 
1436  ! local variables
1437  integer :: jd ! counter
1438  integer :: rd, cd ! global counter for row and col
1439  integer :: rdl, cdl ! local counter for row and col
1440  integer :: i_rd, i_cd ! index of subrow and subcol
1441  integer :: rich_lvl_name ! either the Richardson level or zero, to append to names
1442  integer :: sec_x_loc ! local sec_X
1443  integer :: n_loc(2) ! n_loc for n_mod_X instead of n_bnd
1444  integer :: lims_cl(2) ! lims_c_loc in local coordinates
1445  integer :: par_id ! parallel index
1446  integer :: alpha_id ! field line label
1447  integer, allocatable :: lims_r_loc(:,:) ! local row limits for n_loc
1448  integer, allocatable :: lims_c_loc(:,:) ! local column limits for n_loc
1449  integer, target :: desc_phiep(blacsctxtsize) ! descriptor for Phi and EP
1450  integer, target :: desc_gep(blacsctxtsize) ! descriptor for GEP
1451  integer, target :: desc_res2(blacsctxtsize) ! descriptor for res2
1452  real(dp) :: arg_loc ! local argument
1453  real(dp) :: dpar_x ! step size in par_X
1454  real(dp) :: dalpha ! step size in alpha
1455  real(dp) :: i_loc ! local integration rule
1456  real(dp) :: par_loc ! local parallel angle
1457  real(dp) :: alpha_loc ! local field line label
1458  real(dp), allocatable, target :: ep(:,:) ! EP
1459  real(dp), allocatable, target :: gep(:,:) ! GEP
1460  real(dp), allocatable, target :: phi(:,:) ! intermediate Phi
1461  real(dp), allocatable, target :: res2(:,:) ! vacuum response, real version
1462  real(dp), allocatable, target :: res2_loc(:,:) ! local version of res2
1463  logical :: rank_o ! output rank
1464 #if ldebug
1465  integer :: lwork ! size of work array
1466  real(dp), allocatable :: res2_ev(:,:) ! eigenvalues of res2
1467  real(dp), allocatable :: tau(:) ! array tau in upper Hessenberg form
1468  real(dp), allocatable :: work(:) ! work array
1469  real(dp), allocatable :: x_plot(:,:) ! X of plot
1470  real(dp), allocatable :: y_plot(:,:) ! Y of plot
1471  complex(dp), allocatable :: res_loc(:,:) ! local vac%res
1472  complex(dp), allocatable :: res_ev(:) ! eigenvalues of res2
1473  complex(dp), allocatable :: work_c(:) ! work array
1474  complex(dp), allocatable :: tau_c(:) ! array tau in upper Hessenberg form
1475 #endif
1476 
1477  ! initialize ierr
1478  ierr = 0
1479 
1480  ! decide whether to append Richardson level to name
1481  select case (eq_style)
1482  case (1) ! VMEC
1483  rich_lvl_name = rich_lvl ! append richardson level
1484  case (2) ! HELENA
1485  rich_lvl_name = 0 ! do not append name
1486  end select
1487 
1488  if (rich_lvl.eq.rich_restart_lvl .and. jump_to_sol) then
1489  ! reconstruct old vacuum
1490  ierr = reconstruct_pb3d_vac(vac,'vac',rich_lvl=rich_lvl_name)
1491  chckerr('')
1492  return
1493  end if
1494 
1495  call writo('Start calculation of vacuum quantities')
1496 
1497  call lvl_ud(1)
1498 
1499  ! set ouptut rank persistency
1500  if (rank.eq.n_procs-1) then
1501  rank_o = .true.
1502  else
1503  rank_o = .false.
1504  end if
1505 
1506  ! calculate matrices G and H
1507  ierr = calc_gh(vac)
1508  chckerr('')
1509 
1510  ! Only for processes that are in the blacs context
1511  if (in_context(vac%ctxt_HG)) then
1512  ! user output
1513  call writo('Use STRUMPack to solve system',persistent=rank_o)
1514  call lvl_ud(1)
1515 
1516  ! set step sizes if 3-D vacuum
1517  if (vac%style.eq.1) then
1518  dpar_x = (max_par_x-min_par_x)/(n_par_x-1)*pi
1519  if (n_alpha.gt.1) then
1520  dalpha = (max_alpha-min_alpha)/(n_alpha-1)*pi
1521  else
1522  dalpha = 0._dp
1523  end if
1524  end if
1525 
1526  ! set secondary mode numbers
1527  allocate(vac%sec_X(n_mod_x))
1528  if (use_pol_flux_f) then
1529  vac%sec_X = mds%m(size(mds%m,1),:)
1530  else
1531  vac%sec_X = mds%n(size(mds%n,1),:)
1532  end if
1533 
1534  ! set sizes
1535  do jd = 1,2
1536  n_loc(jd) = numroc(2*n_mod_x,vac%bs,vac%ind_p(jd),0,vac%n_p(jd))
1537  end do
1538  call set_loc_lims(n_loc(1),vac%bs,vac%ind_p(1),vac%n_p(1),&
1539  &lims_r_loc)
1540  call set_loc_lims(n_loc(2),vac%bs,vac%ind_p(2),vac%n_p(2),&
1541  &lims_c_loc)
1542 
1543  ! initialize descriptors
1544  call descinit(desc_phiep,vac%n_bnd,2*n_mod_x,vac%bs,vac%bs,0,0,&
1545  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
1546  chckerr('descinit failed for EP')
1547  call descinit(desc_gep,vac%n_bnd,2*n_mod_x,vac%bs,vac%bs,0,0,&
1548  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
1549  chckerr('descinit failed for GEP')
1550  call descinit(desc_res2,2*n_mod_x,2*n_mod_x,vac%bs,vac%bs,0,0,&
1551  &vac%ctxt_HG,max(1,n_loc(1)),ierr)
1552  chckerr('descinit failed for res2')
1553 
1554  ! allocate auxililary matrices
1555  allocate(ep(vac%n_loc(1),n_loc(2)))
1556  allocate(gep(vac%n_loc(1),n_loc(2)))
1557  allocate(phi(vac%n_loc(1),n_loc(2)))
1558  allocate(res2(n_loc(1),n_loc(2)))
1559  allocate(res2_loc(2*n_mod_x,2*n_mod_x))
1560 
1561  ! set EP
1562  ! iterate over all subcolumns of matrix with 2n_mod_X columns
1563  subcols: do i_cd = 1,size(lims_c_loc,2)
1564  col: do cd = lims_c_loc(1,i_cd),lims_c_loc(2,i_cd)
1565  ! set local column index
1566  cdl = sum(lims_c_loc(2,1:i_cd-1)-&
1567  &lims_c_loc(1,1:i_cd-1)+1) + &
1568  &cd-lims_c_loc(1,i_cd)+1
1569 
1570  ! set local secondary mode number
1571  sec_x_loc = vac%sec_X(mod(cd-1,n_mod_x)+1)
1572 
1573  ! iterate over all subrows of matrix with n_bnd rows
1574  subrows: do i_rd = 1,size(vac%lims_r,2)
1575  row: do rd = vac%lims_r(1,i_rd),vac%lims_r(2,i_rd)
1576  ! set local row index
1577  rdl = sum(vac%lims_r(2,1:i_rd-1)-&
1578  &vac%lims_r(1,1:i_rd-1)+1) + &
1579  &rd-vac%lims_r(1,i_rd)+1
1580 
1581  ! set parallel index and field line label
1582  par_id = mod(rd-1,vac%n_ang(1))+1
1583  alpha_id = (rd-1)/vac%n_ang(1) + 1
1584 
1585  ! set up local argument
1586  select case (vac%style)
1587  case (1) ! field-line 3-D
1588  par_loc = pi*min_par_x + &
1589  &(par_id-1)*dpar_x
1590  alpha_loc = pi*min_alpha + &
1591  &(alpha_id-1)*dalpha
1592  if (use_pol_flux_f) then
1593  arg_loc = vac%prim_X*alpha_loc + &
1594  &(vac%prim_X*vac%jq-sec_x_loc)*&
1595  &par_loc
1596  else
1597  arg_loc = vac%prim_X*alpha_loc + &
1598  &(sec_x_loc-vac%prim_X*vac%jq)*&
1599  &par_loc
1600  end if
1601  case (2) ! axisymmetric
1602  arg_loc = -sec_x_loc*&
1603  &vac%ang(par_id,alpha_id)
1604  end select
1605 
1606  ! set up argument
1607 
1608  ! save EP
1609  if (cd.le.n_mod_x) then ! real part of rhs
1610  ep(rdl,cdl) = cos(arg_loc) ! rp(i exp)
1611  else
1612  ep(rdl,cdl) = sin(arg_loc) ! ip(i exp)
1613  end if
1614  if (use_pol_flux_f) then
1615  ep(rdl,cdl) = (vac%prim_X*vac%jq-sec_x_loc)*&
1616  &ep(rdl,cdl)
1617  else
1618  ep(rdl,cdl) = (sec_x_loc-vac%prim_X*vac%jq)*&
1619  &ep(rdl,cdl)
1620  end if
1621  end do row
1622  end do subrows
1623  end do col
1624  end do subcols
1625  call plot_hdf5('EP','EP',reshape(ep,[vac%n_bnd,n_loc(2),1]))
1626 
1627  ! solve for Phi
1628  ierr = solve_phi_bem(vac,ep,phi,[vac%n_bnd,2*n_mod_x],&
1629  &[vac%n_loc(1),n_loc(2)],lims_c_loc,desc_phiep)
1630  chckerr('')
1631 
1632  ! calculate GEP
1633  call pdgemm('N','N',vac%n_bnd,2*n_mod_x,vac%n_bnd,1._dp,vac%G,1,1,&
1634  &vac%desc_G,ep,1,1,desc_phiep,0._dp,gep,1,1,desc_gep)
1635 
1636  call lvl_ud(-1)
1637 
1638  ! user output
1639  call writo('Combine results into vacuum response',persistent=rank_o)
1640  call lvl_ud(1)
1641 
1642  ! Convert EP into IEP where I contains the integration rule:
1643  ! (θ_{i+1}-θ_{i-1})/2 for i = 2..n-1
1644  ! I_i = (θ_2-θ_1}/2 for i = 1
1645  ! (θ_{n}-θ_{n-1}}/2 for i = n
1646  ! with n = vac%n_bnd.
1647  ! For vacua of type 1, the step size is constant. Also, there is an
1648  ! additional step size in the alpha direction.
1649  subrows3: do i_rd = 1,size(vac%lims_r,2)
1650  row3: do rd = vac%lims_r(1,i_rd),vac%lims_r(2,i_rd)
1651  ! set local row index
1652  rdl = sum(vac%lims_r(2,1:i_rd-1)-&
1653  &vac%lims_r(1,1:i_rd-1)+1) + &
1654  &rd-vac%lims_r(1,i_rd)+1
1655 
1656  subcols3: do i_cd = 1,size(lims_c_loc,2)
1657  ! set local column limits
1658  lims_cl = sum(lims_c_loc(2,1:i_cd-1)-&
1659  &lims_c_loc(1,1:i_cd-1)+1) + &
1660  &lims_c_loc(:,i_cd)-lims_c_loc(1,i_cd)+1
1661 
1662  select case (vac%style)
1663  case (1) ! field-line 3-D
1664  i_loc = dpar_x*dalpha
1665  if ((rd-1)/vac%n_ang(1).eq.0 .or. & ! first point on field line
1666  &(rd-1)/vac%n_ang(1).eq.vac%n_ang(1)-1) & ! last point on field line
1667  &i_loc = i_loc*0.5_dp
1668  case (2) ! axisymmetric
1669  i_loc = 0.5_dp*&
1670  &(vac%ang(min(rd+1,vac%n_bnd),1)-&
1671  &vac%ang(max(rd-1,1),1))
1672  end select
1673 
1674  ep(rdl,lims_cl(1):lims_cl(2)) = &
1675  &ep(rdl,lims_cl(1):lims_cl(2)) * i_loc
1676  end do subcols3
1677  end do row3
1678  end do subrows3
1679 
1680  ! calculate (EP)^T Phi
1681  ! Note that this means that the second half of the rows in EP^T are
1682  ! still missing a minus sign if they are to represent EP^†. This
1683  ! will be taken into account when res2_loc is converted to vac%res.
1684  call pdgemm('T','N',2*n_mod_x,2*n_mod_x,vac%n_bnd,1._dp,ep,1,1,&
1685  &desc_phiep,phi,1,1,desc_phiep,0._dp,res2,1,1,desc_res2)
1686 
1687  ! get serial version on last process
1688  ierr = mat_dis2loc(vac%ctxt_HG,res2,lims_r_loc,lims_c_loc,res2_loc,&
1689  &proc=n_procs-1)
1690  chckerr('')
1691 
1692 #if ldebug
1693  ! user output
1694  if (debug_calc_vac_res) then
1695  call writo('Calculate the eigenvalues of res_loc to &
1696  &investigate whether it is positive definite',&
1697  &persistent=rank_o)
1698  call lvl_ud(1)
1699  end if
1700 #endif
1701 
1702  ! convert to complex numbers and save in vac%res
1703  if (rank.eq.n_procs-1) then
1704  vac%res = res2_loc(1:n_mod_x,1:n_mod_x)-&
1705  &(-res2_loc(n_mod_x+1:2*n_mod_x,n_mod_x+1:2*n_mod_x)) ! minus sign for lower rows
1706  vac%res = vac%res + iu * &
1707  &(res2_loc(1:n_mod_x,n_mod_x+1:2*n_mod_x)-&
1708  &res2_loc(n_mod_x+1:2*n_mod_x,1:n_mod_x)) ! minus sign for lower rows
1709  vac%res = vac%res/vac_perm
1710 #if ldebug
1711  allocate(x_plot(n_mod_x,n_mod_x))
1712  allocate(y_plot(n_mod_x,n_mod_x))
1713  do jd = 1,n_mod_x
1714  x_plot(jd,:) = vac%sec_X(jd)
1715  y_plot(:, jd) = vac%sec_X(jd)
1716  end do
1717  call print_ex_2d(['real vacuum response'],'vac_res_Re',&
1718  &rp(vac%res),x=x_plot,draw=.false.)
1719  call draw_ex(['real vacuum response'],'vac_res_Re',n_mod_x,1,&
1720  &.false.)
1721  call print_ex_2d(['imag vacuum response'],'vac_res_Im',&
1722  &ip(vac%res),x=x_plot,draw=.false.)
1723  call draw_ex(['imag vacuum response'],'vac_res_Im',n_mod_x,1,&
1724  &.false.)
1725  call plot_hdf5(['real part','imag part','imag diff'],'vac_res',&
1726  &reshape([rp(vac%res),ip(vac%res),&
1727  &ip(vac%res+transpose(vac%res))],[n_mod_x,n_mod_x,1,3]),&
1728  &x=reshape(x_plot,[n_mod_x,n_mod_x,1,1]),&
1729  &y=reshape(y_plot,[n_mod_x,n_mod_x,1,1]))
1730  deallocate(x_plot,y_plot)
1731 
1732  if (debug_calc_vac_res) then
1733  ! copy vacuum response
1734  allocate(res_loc(n_mod_x,n_mod_x))
1735  res_loc = vac%res
1736 
1737  ! get size of work matrix required for zgehrd
1738  allocate(work_c(1))
1739  allocate(tau_c(n_mod_x-1))
1740  call zgehrd(n_mod_x,1,n_mod_x,res_loc,n_mod_x,tau_c,work_c,&
1741  &-1,ierr)
1742  chckerr('zgehrd failed')
1743  lwork = nint(rp(work_c(1)))
1744  deallocate(work_c)
1745  allocate(work_c(lwork))
1746 
1747  ! convert to upper Heisenberg form
1748  call zgehrd(n_mod_x,1,n_mod_x,res_loc,n_mod_x,tau_c,work_c,&
1749  &lwork,ierr)
1750  chckerr('zgehrd failed')
1751  deallocate(work_c)
1752 
1753  ! get size of work matrix required for zhseqr
1754  allocate(work_c(1))
1755  allocate(res_ev(n_mod_x))
1756  call zhseqr('E','N',n_mod_x,1,n_mod_x,res_loc,n_mod_x,&
1757  &res_ev,0._dp,n_mod_x,work_c,-1,ierr)
1758  chckerr('zhseqr failed')
1759  lwork = nint(rp(work_c(1)))
1760  deallocate(work_c)
1761  allocate(work_c(lwork))
1762 
1763  ! calculate eigenvalues
1764  call zhseqr('E','N',n_mod_x,1,n_mod_x,res_loc,n_mod_x,&
1765  &res_ev,0._dp,n_mod_x,work_c,lwork,ierr)
1766  chckerr('zhseqr failed')
1767  call print_ex_2d(['real part','imag part'],'res_ev',&
1768  &reshape([rp(res_ev),ip(res_ev)],[n_mod_x,2]),&
1769  &draw=.false.)
1770  call draw_ex(['real part','imag part'],'res_ev',2,1,&
1771  &.false.)
1772  deallocate(work_c)
1773  deallocate(res_ev)
1774  end if
1775 #endif
1776  end if
1777 
1778 #if ldebug
1779  ! calculate the Eigenvalues of the original real res2
1780  if (debug_calc_vac_res) then
1781  ! get size of work matrix required for pdgehrd
1782  allocate(work(1))
1783  allocate(tau(2*n_mod_x-1))
1784  call pdgehrd(2*n_mod_x,1,2*n_mod_x,res2,1,1,desc_res2,tau,work,&
1785  &-1,ierr)
1786  chckerr('pdgehrd failed')
1787  lwork = nint(work(1))
1788  deallocate(work)
1789  allocate(work(lwork))
1790 
1791  ! convert to upper Heisenberg form
1792  call pdgehrd(2*n_mod_x,1,2*n_mod_x,res2,1,1,desc_res2,tau,&
1793  &work,lwork,ierr)
1794  chckerr('pdgehrd failed')
1795  deallocate(work)
1796  deallocate(tau)
1797 
1798  ! calculate eigenvalues
1799  lwork = 2*n_mod_x / vac%bs
1800  if (lwork*vac%bs.lt.2*n_mod_x) lwork = lwork + 1
1801  lwork = 7*lwork / lcm(vac%n_p(1),vac%n_p(2))
1802  lwork = 3*2*n_mod_x + max(2*2*n_mod_x+2*n_loc(2) , lwork)
1803  allocate(work(lwork))
1804  allocate(res2_ev(2*n_mod_x,2))
1805  call pdlahqr(.false.,.false.,2*n_mod_x,1,2*n_mod_x,res2,&
1806  &desc_res2,res2_ev(:,1),res2_ev(:,2),1,2*n_mod_x,[0._dp],&
1807  &desc_res2,work,lwork,[0],1,ierr)
1808  chckerr('pdlahqr failed')
1809  call print_ex_2d(['real part','imag part'],'res2_ev',&
1810  &res2_ev,draw=.false.)
1811  call draw_ex(['real part','imag part'],'res2_ev',2,1,.false.)
1812  deallocate(work)
1813  deallocate(res2_ev)
1814 
1815  ! user output
1816  call lvl_ud(-1)
1817  call writo('They should all be positive',persistent=rank_o)
1818  end if
1819 #endif
1820 
1821  call lvl_ud(-1)
1822  end if
1823 
1824  call lvl_ud(-1)
1825 
1826  call writo('Done calculating vacuum quantities')
1827  end function calc_vac_res
1828 
1856  integer function vac_pot_plot(grid_sol,vac,sol,X_id) result(ierr)
1857  use num_vars, only: min_zeta_plot, max_zeta_plot, min_rvac_plot, &
1858  &max_rvac_plot, min_zvac_plot, max_zvac_plot, n_vac_plot, &
1859  &n_zeta_plot, n_procs, rank, use_pol_flux_f
1860  use sol_vars, only: sol_type
1861  use mpi_utilities, only: broadcast_var
1862  use iso_c_binding
1863  use vac_vars, only: in_context
1864  use vac_utilities, only: mat_dis2loc
1865  use x_vars, only: prim_x, n_mod_x
1866  use vac_vars, only: set_loc_lims
1867 
1868  ! input / output
1869  type(grid_type), intent(in) :: grid_sol
1870  type(vac_type), intent(inout) :: vac
1871  type(sol_type), intent(in) :: sol
1872  integer, intent(in) :: x_id
1873 
1874  ! local variables
1875  integer :: id, jd, kd ! counters
1876  integer :: rd, cd ! global counter for row and col
1877  integer :: i_rd, i_cd ! index of subrow and subcol
1878  integer :: rdl, cdl ! local counter for row and col
1879  integer :: desc_rhs(blacsctxtsize) ! descriptor for RHS
1880  integer :: desc_gh_loc(blacsctxtsize) ! descriptor for local G and H
1881  integer :: desc_phi(blacsctxtsize) ! descriptor for Phi
1882  integer :: n_row_plot ! n_loc(1) for n_vac_plot rows
1883  integer :: n_col_2 ! n_loc(2) for 2 columns (real and imaginary part)
1884  integer :: lims_cl(2) ! lims_c_loc in local coordinates
1885  integer, allocatable :: lims_r_loc(:,:) ! local column limits for n_row_plot
1886  integer, allocatable :: lims_c_loc(:,:) ! local column limits for n_col_2
1887  integer :: sec_x_loc ! local sec_X
1888  integer :: fac_x(2) ! (n,-m) for poloidal or (-m,n) for toroidal flux
1889  real(dp) :: zeta_loc ! local zeta
1890  real(dp), allocatable :: x_vec(:,:) ! x_vec of influence points
1891  real(dp), allocatable :: g_loc(:,:) ! local G
1892  real(dp), allocatable :: h_loc(:,:) ! local H
1893  real(dp), allocatable :: rhs(:,:,:) ! 2 right-hand sides, real and complex
1894  real(dp), allocatable :: phi(:,:) ! Phi
1895  real(dp), allocatable :: phi_ser(:,:) ! serial version of Phi
1896  real(dp), allocatable :: phi_2d(:,:) ! 2-D serial version of Phi
1897  complex(dp), allocatable :: sol_vec(:) ! solution vector at edge
1898  complex(dp) :: rhs_complex ! complex local RHS
1899  character(len=max_str_ln) :: err_msg ! error message
1900 #if ldebug
1901  real(dp), allocatable :: rhs_loc(:,:) ! local version of RHS
1902  !character(len=max_str_ln) :: plot_title ! plot title
1903  !character(len=max_str_ln) :: plot_name ! plot name
1904 #endif
1905 
1906  character(*), parameter :: rout_name = 'vac_pot_plot'
1907 
1908  ! initialize ierr
1909  ierr = 0
1910 
1911  ! distribute solution vector
1912  allocate(sol_vec(n_mod_x))
1913  if (rank.eq.n_procs-1) sol_vec = sol%vec(:,grid_sol%loc_n_r,x_id)
1914  ierr = broadcast_var(sol_vec,n_procs-1)
1915  chckerr('')
1916 
1917  ! calculate matrices G and H
1918  ierr = calc_gh(vac)
1919  chckerr('')
1920 
1921  ! Only for processes that are in the blacs context
1922  if (in_context(vac%ctxt_HG)) then
1923  ! set sizes
1924  n_row_plot = numroc(product(n_vac_plot),vac%bs,vac%ind_p(1),0,&
1925  &vac%n_p(1))
1926  n_col_2 = numroc(2,vac%bs,vac%ind_p(2),0,vac%n_p(2))
1927  call set_loc_lims(n_row_plot,vac%bs,vac%ind_p(1),vac%n_p(1),&
1928  &lims_r_loc)
1929  call set_loc_lims(n_col_2,vac%bs,vac%ind_p(2),vac%n_p(2),&
1930  &lims_c_loc)
1931 
1932  ! allocate RHS and set descriptor
1933  allocate(rhs(vac%n_loc(1),n_col_2,2))
1934  rhs = 0._dp
1935  call descinit(desc_rhs,vac%n_bnd,2,vac%bs,vac%bs,0,0,&
1936  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
1937  chckerr('descinit failed for RHS')
1938 
1939  ! initialize the right hand sides as the normal derivative of Phi at
1940  ! the plasma edge
1941  subcols: do i_cd = 1,size(lims_c_loc,2)
1942  col: do cd = lims_c_loc(1,i_cd),lims_c_loc(2,i_cd)
1943  ! set local column index
1944  cdl = sum(lims_c_loc(2,1:i_cd-1)-&
1945  &lims_c_loc(1,1:i_cd-1)+1) + &
1946  &cd-lims_c_loc(1,i_cd)+1
1947 
1948  ! iterate over all subrows of matrix with n_bnd rows
1949  subrows: do i_rd = 1,size(vac%lims_r,2)
1950  row: do rd = vac%lims_r(1,i_rd),vac%lims_r(2,i_rd)
1951  ! set local row index
1952  rdl = sum(vac%lims_r(2,1:i_rd-1)-&
1953  &vac%lims_r(1,1:i_rd-1)+1) + &
1954  &rd-vac%lims_r(1,i_rd)+1
1955 
1956  ! sum over all modes to get parallel derivative of X
1957  do jd = 1,n_mod_x
1958  sec_x_loc = vac%sec_X(jd)
1959  if (use_pol_flux_f) then
1960  fac_x = [prim_x,-sec_x_loc]
1961  else
1962  ierr = 1
1963  err_msg = 'Currently must use poloidal flux'
1964  chckerr(err_msg)
1965  !fac_X = [-prim_X,sec_X_loc]
1966  end if
1967 
1968  rhs_complex = iu*(fac_x(1)*vac%jq+fac_x(2))*&
1969  &exp(iu*fac_x(2)*vac%ang(rd,1))*sol_vec(jd)
1970  if (cd.eq.1) rhs(rdl,cdl,:) = &
1971  &rhs(rdl,cdl,:) + rp(rhs_complex)
1972  if (cd.eq.2) rhs(rdl,cdl,:) = &
1973  &rhs(rdl,cdl,:) + ip(rhs_complex)
1974  end do
1975  end do row
1976  end do subrows
1977  end do col
1978  end do subcols
1979 
1980  ! solve for Phi and save in second index of RHS
1981  ierr = solve_phi_bem(vac,rhs(:,:,1),rhs(:,:,2),[vac%n_bnd,2],&
1982  &[vac%n_loc(1),n_col_2],lims_c_loc,desc_rhs)
1983  chckerr('')
1984 
1985 #if ldebug
1986  ! get serial version on last process
1987  allocate(rhs_loc(vac%n_bnd,2))
1988  do jd = 1,2
1989  ierr = mat_dis2loc(vac%ctxt_HG,rhs(:,:,jd),vac%lims_r,&
1990  &lims_c_loc,rhs_loc,proc=n_procs-1)
1991  chckerr('')
1992  if (rank.eq.n_procs-1) then
1993  call print_ex_2d(['dphi','phi '],'RHS_'//trim(i2str(jd)),&
1994  &rhs_loc,x=vac%ang(:,1:1),draw=.false.)
1995  call draw_ex(['dphi','phi '],'RHS_'//trim(i2str(jd)),2,1,&
1996  &.false.)
1997  end if
1998  end do
1999  deallocate(rhs_loc)
2000 #endif
2001 
2002  ! convert RHS into IRHS where I contains the integration rule:
2003  ! (θ_{i+1}-θ_{i-1})/2 for i = 2..n-1
2004  ! I_i = (θ_2-θ_1}/2 for i = 1
2005  ! (θ_{n}-θ_{n-1}}/2 for i = n
2006  ! with n = vac%n_bnd.
2007  subrows2: do i_rd = 1,size(vac%lims_r,2)
2008  row2: do rd = vac%lims_r(1,i_rd),vac%lims_r(2,i_rd)
2009  ! set local row index
2010  rdl = sum(vac%lims_r(2,1:i_rd-1)-&
2011  &vac%lims_r(1,1:i_rd-1)+1) + &
2012  &rd-vac%lims_r(1,i_rd)+1
2013 
2014  subcols2: do i_cd = 1,size(lims_c_loc,2)
2015  ! local limits
2016  lims_cl = sum(lims_c_loc(2,1:i_cd-1)-&
2017  &lims_c_loc(1,1:i_cd-1)+1) + &
2018  &lims_c_loc(:,i_cd)-lims_c_loc(1,i_cd)+1
2019 
2020  rhs(rdl,lims_cl(1):lims_cl(2),:) = &
2021  &rhs(rdl,lims_cl(1):lims_cl(2),:) * 0.5_dp*&
2022  &(vac%ang(min(rd+1,vac%n_bnd),1)-&
2023  &vac%ang(max(rd-1,1),1))
2024  end do subcols2
2025  end do row2
2026  end do subrows2
2027 
2028  ! allocate G, H and x_vec and set descriptor
2029  allocate(g_loc(n_row_plot,vac%n_loc(2)))
2030  allocate(h_loc(n_row_plot,vac%n_loc(2)))
2031  allocate(x_vec(product(n_vac_plot),2))
2032  g_loc = 0._dp
2033  h_loc = 0._dp
2034  call descinit(desc_gh_loc,product(n_vac_plot),vac%n_bnd,vac%bs,&
2035  &vac%bs,0,0,vac%ctxt_HG,max(1,n_row_plot),ierr)
2036  chckerr('descinit failed for GH_loc')
2037 
2038  ! loop over all toroidal points
2039  do kd = 1,n_zeta_plot
2040  ! set local zeta
2041  zeta_loc = min_zeta_plot + (max_zeta_plot-min_zeta_plot) * &
2042  &(kd-1._dp)/max(1._dp,n_zeta_plot-1._dp)
2043 
2044  ! set up variables to calculate local G and H
2045  row3: do rd = 1,product(n_vac_plot)
2046  ! set up 2-D index for x_vec
2047  id = mod(rd-1,n_vac_plot(1))+1
2048  jd = (rd-1)/n_vac_plot(1)+1
2049 
2050  ! set up x_vec
2051  x_vec(rd,1) = min_rvac_plot + &
2052  &(max_rvac_plot-min_rvac_plot)*&
2053  &(id-1._dp)/(n_vac_plot(1)-1._dp)
2054  x_vec(rd,2) = min_zvac_plot + &
2055  &(max_zvac_plot-min_zvac_plot)*&
2056  &(jd-1._dp)/(n_vac_plot(2)-1._dp)
2057  end do row3
2058  !call print_ex_2D(['R','Z'],'x_vec_'//trim(i2str(rank)),x_vec,draw=.false.)
2059  !call draw_ex(['R','Z'],'x_vec_'//trim(i2str(rank)),2,1,.false.)
2060 
2061  ! calc local G and H
2062  ierr = calc_gh(vac,n_r_in=product(n_vac_plot),&
2063  &lims_r_in=lims_r_loc,x_vec_in=x_vec,g_in=g_loc,h_in=h_loc)
2064  chckerr('')
2065 
2066 #if ldebug
2067  !if (debug_vac_pot_plot) then
2068  !subrows4: do i_rd = 1,size(lims_r_loc,2)
2069  !row4: do rd = lims_r_loc(1,i_rd),lims_r_loc(2,i_rd)
2070  !subcols4: do i_cd = 1,size(vac%lims_c,2)
2071  !! output
2072  !plot_title = 'for row '//trim(i2str(rd))//&
2073  !&' and starting column '//&
2074  !&trim(i2str(vac%lims_c(1,i_cd)))
2075  !plot_name = '_'//trim(i2str(kd))//'_'//&
2076  !&trim(i2str(rd))//'_'//&
2077  !&trim(i2str(vac%lims_c(1,i_cd)))
2078 
2079  !! local limits
2080  !lims_cl = sum(vac%lims_c(2,1:i_cd-1)-&
2081  !&vac%lims_c(1,1:i_cd-1)+1) + &
2082  !&vac%lims_c(:,i_cd)-vac%lims_c(1,i_cd)+1
2083 
2084  !! output
2085  !call print_ex_2D('G '//trim(plot_title),&
2086  !&'G_loc'//trim(plot_name),&
2087  !&G_loc(rdl,lims_cl(1):lims_cl(2)),x=vac%ang&
2088  !&(vac%lims_c(1,i_cd):vac%lims_c(2,i_cd),1),&
2089  !&draw=.false.)
2090  !call draw_ex(['G '//trim(plot_title)],&
2091  !&'G_loc'//trim(plot_name),1,1,.false.)
2092  !call print_ex_2D('H '//trim(plot_title),&
2093  !&'H_loc'//trim(plot_name),&
2094  !&H_loc(rdl,lims_cl(1):lims_cl(2)),x=vac%ang&
2095  !&(vac%lims_c(1,i_cd):vac%lims_c(2,i_cd),1),&
2096  !&draw=.false.)
2097  !call draw_ex(['H '//trim(plot_title)],&
2098  !&'H_loc'//trim(plot_name),1,1,.false.)
2099  !end do subcols4
2100  !end do row4
2101  !end do subrows4
2102  !end if
2103 #endif
2104 
2105  ! multiply G with RHS(1) and save in Phi
2106  allocate(phi(n_row_plot,n_col_2))
2107  call descinit(desc_phi,product(n_vac_plot),2,vac%bs,&
2108  &vac%bs,0,0,vac%ctxt_HG,max(1,n_row_plot),ierr)
2109  chckerr('descinit failed for Phi')
2110  call pdgemm('N','N',product(n_vac_plot),2,vac%n_bnd,1._dp,&
2111  &g_loc,1,1,desc_gh_loc,rhs(:,:,1),1,1,desc_rhs,0._dp,&
2112  &phi,1,1,desc_phi)
2113 
2114  ! add -H multiplied by RHS(2)
2115  call pdgemm('N','N',product(n_vac_plot),2,vac%n_bnd,-1._dp,&
2116  &h_loc,1,1,desc_gh_loc,rhs(:,:,2),1,1,desc_rhs,1._dp,&
2117  &phi,1,1,desc_phi)
2118 
2119  ! get serial version on last process
2120  allocate(phi_ser(product(n_vac_plot),2))
2121  ierr = mat_dis2loc(vac%ctxt_HG,phi,lims_r_loc,&
2122  &lims_c_loc,phi_ser,proc=n_procs-1)
2123  chckerr('')
2124 #if ldebug
2125  if (rank.eq.n_procs-1) then
2126  call print_ex_2d(['dphi','phi '],'Phi'//trim(i2str(kd)),&
2127  &phi_ser,draw=.false.)
2128  call draw_ex(['dphi','phi '],'Phi'//trim(i2str(kd)),2,1,&
2129  &.false.)
2130  end if
2131 #endif
2132 
2133  ! multiply with exponential of toroidal angle and take real part
2134  ! (phi_C = - zeta_F)
2135  if (rank.eq.n_procs-1) then
2136  allocate(phi_2d(n_vac_plot(1),n_vac_plot(2)))
2137  phi_2d = -1./(4*pi)*reshape(&
2138  &phi_ser(:,1)*cos(vac%prim_X*zeta_loc) + &
2139  &phi_ser(:,2)*sin(vac%prim_X*zeta_loc),n_vac_plot)
2140  call plot_hdf5('Phi','Phi',reshape(phi_2d,[n_vac_plot,1]),&
2141  &x=cos(vac%prim_X*zeta_loc)*&
2142  &reshape(x_vec(:,1),[n_vac_plot,1]),&
2143  &y=sin(vac%prim_X*zeta_loc)*&
2144  &reshape(x_vec(:,1),[n_vac_plot,1]),&
2145  &z=reshape(x_vec(:,2),[n_vac_plot,1]))
2146  end if
2147  deallocate(x_vec)
2148  end do
2149  end if
2150  end function vac_pot_plot
2151 
2161  integer function solve_phi_bem(vac,R,Phi,n_RPhi,n_loc_RPhi,lims_c_RPhi,&
2162  &desc_RPhi) result(ierr)
2163 
2164  use vac_vars, only: in_context
2165 #if ldebug
2166  use grid_vars, only: min_par_x, max_par_x
2167  use rich_vars, only: n_par_x
2168 #endif
2169 
2170  character(*), parameter :: rout_name = 'solve_Phi_BEM'
2171 
2172  ! input / output
2173  type(vac_type), intent(in), target :: vac
2174  real(dp), intent(inout), target :: phi(:,:)
2175  real(dp), intent(in) :: r(:,:)
2176  integer, intent(in) :: n_rphi(2)
2177  integer, intent(in) :: n_loc_rphi(2)
2178  integer, intent(in) :: lims_c_rphi(:,:)
2179  integer, intent(in), target :: desc_rphi(blacsctxtsize)
2180 
2181  ! local variables
2182  integer, target :: desc_gr(blacsctxtsize) ! descriptor for GR
2183  real(dp) :: sdp_err ! error
2184  real(dp), allocatable, target :: gr(:,:) ! GR
2185  type(c_ptr) :: cdesch, ch ! C pointers to descriptor of H and H itself
2186  type(c_ptr) :: cdescgr, cgr ! C pointers to descriptor of GR and GR itself
2187  type(c_ptr) :: cdescphi, cphi ! C pointers to descriptor of Phi and Phi itself
2188  type(strumpackdensepackage_f90_double) :: sdp_loc ! Strumpack object for local solve H Phi = GR
2189 #if ldebug
2190  integer :: lims_rl(2) ! vac%lims_r in local coordinates
2191  integer :: id ! counter
2192  integer :: cd ! global counter for row and col
2193  integer :: cdl ! local counter for col
2194  integer :: i_rd, i_cd ! index of subrow and subcol
2195  character(len=max_str_ln) :: plot_title ! plot title
2196  character(len=max_str_ln) :: plot_name ! plot name
2197  real(dp), allocatable :: ang(:) ! angle
2198 #endif
2199 
2200  ! initialize ierr
2201  ierr = 0
2202 
2203  ! Only for processes that are in the blacs context
2204  if (in_context(vac%ctxt_HG)) then
2205  ! allocate auxililary matrices
2206  allocate(gr(vac%n_loc(1),n_loc_rphi(2)))
2207 
2208  ! initialize descriptors
2209  call descinit(desc_gr,vac%n_bnd,n_rphi(2),vac%bs,vac%bs,0,0,&
2210  &vac%ctxt_HG,max(1,vac%n_loc(1)),ierr)
2211  chckerr('descinit failed for GR')
2212 
2213  ! calculate GR
2214  call pdgemm('N','N',vac%n_bnd,n_rphi(2),vac%n_bnd,1._dp,vac%G,1,1,&
2215  &vac%desc_G,r,1,1,desc_rphi,0._dp,gr,1,1,desc_gr)
2216 
2217  ! set C pointers
2218  ch = c_loc(vac%H)
2219  cdesch = c_loc(vac%desc_H)
2220  cgr = c_loc(gr)
2221  cdescgr = c_loc(desc_gr)
2222  cphi = c_loc(phi)
2223  cdescphi = c_loc(desc_rphi)
2224 
2225  ! initialize the solver and set parameters
2226  call sdp_f90_double_init(sdp_loc,vac%MPI_Comm)
2227  sdp_loc%use_HSS=1
2228  sdp_loc%levels_HSS=4
2229  sdp_loc%min_rand_HSS=10
2230  sdp_loc%lim_rand_HSS=5
2231  sdp_loc%inc_rand_HSS=10
2232  sdp_loc%max_rand_HSS=100
2233  sdp_loc%tol_HSS=1d-12
2234  sdp_loc%steps_IR=10
2235  sdp_loc%tol_IR=1d-10
2236 
2237  ! compression
2238  call sdp_f90_double_compress_a(sdp_loc,ch,cdesch)
2239 
2240 #if ldebug
2241  ! accuracy checking
2242  sdp_err = sdp_f90_double_check_compression(sdp_loc,ch,cdesch)
2243 #endif
2244 
2245  ! factorization
2246  call sdp_f90_double_factor(sdp_loc,ch,cdesch)
2247 
2248  ! solve H Phi = G R
2249  call sdp_f90_double_solve(sdp_loc,cphi,cdescphi,cgr,cdescgr)
2250 
2251 #if ldebug
2252  ! accuracy checking
2253  sdp_err = sdp_f90_double_check_solution(sdp_loc,ch,cdesch,cphi,&
2254  &cdescphi,cgr,cdescgr)
2255 #endif
2256 
2257  ! iterative refinement
2258  sdp_err = sdp_f90_double_refine(sdp_loc,ch,cdesch,cphi,cdescphi,&
2259  &cgr,cdescgr)
2260 
2261 #if ldebug
2262  ! statistics
2263  call sdp_f90_double_print_statistics(sdp_loc)
2264 
2265  ! user output
2266  subcols: do i_cd = 1,size(lims_c_rphi,2)
2267  col: do cd = lims_c_rphi(1,i_cd),lims_c_rphi(2,i_cd)
2268  ! set local column index
2269  cdl = sum(lims_c_rphi(2,1:i_cd-1)-&
2270  &lims_c_rphi(1,1:i_cd-1)+1) + &
2271  &cd-lims_c_rphi(1,i_cd)+1
2272 
2273  subrows2: do i_rd = 1,size(vac%lims_r,2)
2274  ! set local row limits
2275  lims_rl = sum(vac%lims_r(2,1:i_rd-1)-&
2276  &vac%lims_r(1,1:i_rd-1)+1) + &
2277  &vac%lims_r(:,i_rd)-vac%lims_r(1,i_rd)+1
2278 
2279  ! set angle
2280  allocate(ang(vac%lims_r(2,i_rd)-vac%lims_r(1,i_rd)+1))
2281  select case (vac%style)
2282  case (1) ! field-line 3-D
2283  do id = 1,vac%n_ang(1)
2284  ang(id:&
2285  &id+vac%n_ang(1)*(vac%n_ang(2)-1):&
2286  &vac%n_ang(1)) = &
2287  &pi * (min_par_x + (id-1)*&
2288  &(max_par_x-min_par_x)/(n_par_x-1))
2289  end do
2290  case (2) ! axisymmetric
2291  ang = reshape(vac%ang,[vac%n_bnd])
2292  end select
2293 
2294  ! output
2295  plot_title = 'for column '//trim(i2str(cd))//&
2296  &' and starting row '//&
2297  &trim(i2str(vac%lims_r(1,i_rd)))
2298  plot_name = '_'//trim(i2str(cd))//'_'//&
2299  &trim(i2str(vac%lims_r(1,i_rd)))
2300  call print_ex_2d('Phi '//trim(plot_title),&
2301  &'Phi'//trim(plot_name),&
2302  &phi(lims_rl(1):lims_rl(2),cdl),&
2303  &x=ang(vac%lims_r(1,i_rd):vac%lims_r(2,i_rd)),&
2304  &draw=.false.)
2305  call draw_ex(['Phi '//trim(plot_title)],&
2306  &'Phi'//trim(plot_name),1,1,.false.)
2307 
2308  ! clean up
2309  deallocate(ang)
2310  end do subrows2
2311  end do col
2312  end do subcols
2313 #endif
2314 
2315  ! destroy Strumpack object
2316  call sdp_f90_double_destroy(sdp_loc)
2317  end if
2318  end function solve_phi_bem
2319 
2346  integer function print_output_vac(vac,data_name,rich_lvl) result(ierr)
2347  use x_vars, only: n_mod_x
2348  use num_vars, only: pb3d_name
2349  use hdf5_ops, only: print_hdf5_arrs
2350  use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
2351  &max_dim_var_1d
2352 
2353  character(*), parameter :: rout_name = 'print_output_vac'
2354 
2355  ! input / output
2356  type(vac_type), intent(in) :: vac
2357  character(len=*), intent(in) :: data_name
2358  integer, intent(in), optional :: rich_lvl
2359 
2360  ! local variables
2361  type(var_1d_type), allocatable, target :: vac_1d(:) ! 1D equivalent of vac variables
2362  type(var_1d_type), pointer :: vac_1d_loc => null() ! local element in vac_1D
2363  integer :: id ! counter
2364 
2365  ! initialize ierr
2366  ierr = 0
2367 
2368  ! user output
2369  call writo('Write vacuum variables to output file')
2370  call lvl_ud(1)
2371 
2372  ! Set up the 1D equivalents of the perturbation variables
2373  allocate(vac_1d(max_dim_var_1d))
2374 
2375  ! set up variables vac_1D
2376  id = 1
2377 
2378  ! misc_vac
2379  vac_1d_loc => vac_1d(id); id = id+1
2380  vac_1d_loc%var_name = 'misc_vac'
2381  allocate(vac_1d_loc%tot_i_min(1),vac_1d_loc%tot_i_max(1))
2382  allocate(vac_1d_loc%loc_i_min(1),vac_1d_loc%loc_i_max(1))
2383  vac_1d_loc%loc_i_min = [1]
2384  vac_1d_loc%loc_i_max = [6]
2385  vac_1d_loc%tot_i_min = vac_1d_loc%loc_i_min
2386  vac_1d_loc%tot_i_max = vac_1d_loc%loc_i_max
2387  allocate(vac_1d_loc%p(6))
2388  vac_1d_loc%p = [vac%style*1._dp,vac%n_bnd*1._dp,vac%prim_X*1._dp,&
2389  &vac%n_ang*1._dp,vac%jq]
2390 
2391  ! sec_X
2392  vac_1d_loc => vac_1d(id); id = id+1
2393  vac_1d_loc%var_name = 'sec_X'
2394  allocate(vac_1d_loc%tot_i_min(1),vac_1d_loc%tot_i_max(1))
2395  allocate(vac_1d_loc%loc_i_min(1),vac_1d_loc%loc_i_max(1))
2396  vac_1d_loc%tot_i_min = [1]
2397  vac_1d_loc%tot_i_max = n_mod_x
2398  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2399  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2400  allocate(vac_1d_loc%p(n_mod_x))
2401  vac_1d_loc%p = vac%sec_X*1._dp
2402 
2403  ! norm
2404  vac_1d_loc => vac_1d(id); id = id+1
2405  vac_1d_loc%var_name = 'norm'
2406  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2407  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2408  vac_1d_loc%tot_i_min = [1,1]
2409  vac_1d_loc%tot_i_max = shape(vac%norm)
2410  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2411  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2412  allocate(vac_1d_loc%p(size(vac%norm)))
2413  vac_1d_loc%p = reshape(vac%norm,[size(vac%norm)])
2414 
2415  ! x_vec
2416  vac_1d_loc => vac_1d(id); id = id+1
2417  vac_1d_loc%var_name = 'x_vec'
2418  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2419  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2420  vac_1d_loc%tot_i_min = [1,1]
2421  vac_1d_loc%tot_i_max = shape(vac%x_vec)
2422  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2423  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2424  allocate(vac_1d_loc%p(size(vac%x_vec)))
2425  vac_1d_loc%p = reshape(vac%x_vec,[size(vac%x_vec)])
2426 
2427  ! RE_res
2428  vac_1d_loc => vac_1d(id); id = id+1
2429  vac_1d_loc%var_name = 'RE_res'
2430  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2431  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2432  vac_1d_loc%tot_i_min = [1,1]
2433  vac_1d_loc%tot_i_max = [n_mod_x,n_mod_x]
2434  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2435  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2436  allocate(vac_1d_loc%p(size(vac%res)))
2437  vac_1d_loc%p = reshape(rp(vac%res),[size(vac%res)])
2438 
2439  ! IM_res
2440  vac_1d_loc => vac_1d(id); id = id+1
2441  vac_1d_loc%var_name = 'IM_res'
2442  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2443  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2444  vac_1d_loc%tot_i_min = [1,1]
2445  vac_1d_loc%tot_i_max = [n_mod_x,n_mod_x]
2446  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2447  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2448  allocate(vac_1d_loc%p(size(vac%res)))
2449  vac_1d_loc%p = reshape(ip(vac%res),[size(vac%res)])
2450 
2451  ! copy variables specific to style
2452  select case (vac%style)
2453  case (1) ! field-line 3-D
2454  ! h_fac
2455  vac_1d_loc => vac_1d(id); id = id+1
2456  vac_1d_loc%var_name = 'h_fac'
2457  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2458  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2459  vac_1d_loc%tot_i_min = [1,1]
2460  vac_1d_loc%tot_i_max = shape(vac%h_fac)
2461  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2462  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2463  allocate(vac_1d_loc%p(size(vac%h_fac)))
2464  vac_1d_loc%p = reshape(vac%h_fac,[size(vac%h_fac)])
2465  case (2) ! axisymmetric
2466  ! ang
2467  vac_1d_loc => vac_1d(id); id = id+1
2468  vac_1d_loc%var_name = 'ang'
2469  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2470  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2471  vac_1d_loc%tot_i_min = [1,1]
2472  vac_1d_loc%tot_i_max = shape(vac%ang)
2473  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2474  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2475  allocate(vac_1d_loc%p(size(vac%ang)))
2476  vac_1d_loc%p = reshape(vac%ang,[size(vac%ang)])
2477 
2478  ! dnorm
2479  vac_1d_loc => vac_1d(id); id = id+1
2480  vac_1d_loc%var_name = 'dnorm'
2481  allocate(vac_1d_loc%tot_i_min(2),vac_1d_loc%tot_i_max(2))
2482  allocate(vac_1d_loc%loc_i_min(2),vac_1d_loc%loc_i_max(2))
2483  vac_1d_loc%tot_i_min = [1,1]
2484  vac_1d_loc%tot_i_max = shape(vac%dnorm)
2485  vac_1d_loc%loc_i_min = vac_1d_loc%tot_i_min
2486  vac_1d_loc%loc_i_max = vac_1d_loc%tot_i_max
2487  allocate(vac_1d_loc%p(size(vac%dnorm)))
2488  vac_1d_loc%p = reshape(vac%dnorm,[size(vac%dnorm)])
2489  end select
2490 
2491  ! write
2492  ierr = print_hdf5_arrs(vac_1d(1:id-1),pb3d_name,trim(data_name),&
2493  &rich_lvl=rich_lvl,ind_print=.true.)
2494  chckerr('')
2495 
2496  ! clean up
2497  call dealloc_var_1d(vac_1d)
2498  nullify(vac_1d_loc)
2499 
2500  ! user output
2501  call lvl_ud(-1)
2502  end function print_output_vac
2503 end module vac_ops
num_utilities::c
integer function, public c(ij, sym, n, lim_n)
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
Definition: num_utilities.f90:2556
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition: MPI_utilities.f90:55
on_field_line
logical function on_field_line(cd, kd, n_ang)
Definition: vac_ops.f90:965
vac_utilities
Numerical utilities related to the vacuum response.
Definition: vac_utilities.f90:6
vac_ops::calc_gh_1
integer function calc_gh_1(vac, n_r_in, lims_r_in, x_vec_in, G_in, H_in, ext_in)
Calculates matrices G and H in 3-D configuration.
Definition: vac_ops.f90:779
dtorh
Calculation of toroidal functions and .
Definition: dtorh.f90:8
output_ops::plot_diff_hdf5
subroutine, public plot_diff_hdf5(A, B, file_name, tot_dim, loc_offset, descr, output_message)
Takes two input vectors and plots these as well as the relative and absolute difference in a HDF5 fil...
Definition: output_ops.f90:1765
hdf5_vars::dealloc_var_1d
Deallocates 1D variable.
Definition: HDF5_vars.f90:68
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
eq_vars::vac_perm
real(dp), public vac_perm
either usual mu_0 (default) or normalized
Definition: eq_vars.f90:48
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
vac_vars::copy_vac
integer function, public copy_vac(vac, vac_copy)
Copy a vacuum type.
Definition: vac_vars.f90:339
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
vac_ops::debug_vac_pot_plot
logical, public debug_vac_pot_plot
plot debug information for vac_pot_plot()
Definition: vac_ops.f90:48
vac_ops::debug_calc_gh
logical, public debug_calc_gh
plot debug information for calc_GH()
Definition: vac_ops.f90:46
vac_utilities::interlaced_vac_copy
integer function, public interlaced_vac_copy(vac_old, vac)
Copies vacuum variables of a previous level into the appropriate, interlaced place of a next level.
Definition: vac_utilities.f90:355
rich_vars
Variables concerning Richardson extrapolation.
Definition: rich_vars.f90:4
num_ops
Numerical operations.
Definition: num_ops.f90:4
vac_ops::vac_pot_plot
integer function, public vac_pot_plot(grid_sol, vac, sol, X_id)
Calculate vacuum potential at some positions that are not resonant.
Definition: vac_ops.f90:1857
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition: HDF5_ops.f90:27
grid_vars::max_par_x
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
Definition: grid_vars.f90:25
num_vars::n_zeta_plot
integer, public n_zeta_plot
nr. of toroidal points in plot
Definition: num_vars.f90:157
eq_utilities::calc_inv_met
Calculate from and where and , according to .
Definition: eq_utilities.f90:61
vac_ops::calc_gh
integer function calc_gh(vac, n_r_in, lims_r_in, x_vec_in, G_in, H_in)
Calculate the matrices G and H.
Definition: vac_ops.f90:501
num_vars::eq_job_nr
integer, public eq_job_nr
nr. of eq job
Definition: num_vars.f90:79
vac_ops::debug_calc_vac_res
logical, public debug_calc_vac_res
plot debug information for calc_vac_res()
Definition: vac_ops.f90:47
num_vars::min_zvac_plot
real(dp), public min_zvac_plot
min. of Z in which to plot vacuum
Definition: num_vars.f90:167
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
hdf5_vars
Variables pertaining to HDF5 and XDMF.
Definition: HDF5_vars.f90:4
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
vac_ops::store_vac
integer function, public store_vac(grid, eq_1, eq_2, vac)
Stores calculation of the vacuum response by storing the necessary variables.
Definition: vac_ops.f90:113
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
pb3d_ops::reconstruct_pb3d_vac
integer function, public reconstruct_pb3d_vac(vac, data_name, rich_lvl)
Reconstructs the vacuum variables from PB3D HDF5 output.
Definition: PB3D_ops.f90:1228
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
eq_utilities
Numerical utilities related to equilibrium variables.
Definition: eq_utilities.f90:4
num_vars::min_zeta_plot
real(dp), public min_zeta_plot
min. of zeta_plot
Definition: num_vars.f90:161
vac_vars
Variables pertaining to the vacuum quantities.
Definition: vac_vars.f90:4
sol_vars::sol_type
solution type
Definition: sol_vars.f90:30
helena_vars::r_h
real(dp), dimension(:,:), allocatable, public r_h
major radius (xout)
Definition: HELENA_vars.f90:33
vac_ops
Operations and variables pertaining to the vacuum response.
Definition: vac_ops.f90:23
num_vars::max_zvac_plot
real(dp), public max_zvac_plot
max. of Z in which to plot vacuum
Definition: num_vars.f90:168
hdf5_vars::max_dim_var_1d
integer, parameter, public max_dim_var_1d
maximum dimension of var_1D
Definition: HDF5_vars.f90:21
vac_vars::set_loc_lims
subroutine, public set_loc_lims(n_loc, bs, ind_p, n_p, lims)
Calculates the limits in local index.
Definition: vac_vars.f90:261
grid_vars::min_par_x
real(dp), public min_par_x
min. of parallel coordinate [ ] in field-aligned grid
Definition: grid_vars.f90:24
x_vars::prim_x
integer, public prim_x
n_X (pol. flux) or m_X (tor. flux)
Definition: X_vars.f90:126
vac_utilities::vec_dis2loc
integer function, public vec_dis2loc(ctxt, vec_dis, lims, vec_loc, proc)
Gathers a distributed vector on a single process.
Definition: vac_utilities.f90:223
num_vars::rich_restart_lvl
integer, public rich_restart_lvl
starting Richardson level (0: none [default])
Definition: num_vars.f90:173
helena_vars::z_h
real(dp), dimension(:,:), allocatable, public z_h
height (yout)
Definition: HELENA_vars.f90:34
hdf5_vars::var_1d_type
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition: HDF5_vars.f90:48
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
x_vars::x_2_type
tensorial perturbation type
Definition: X_vars.f90:81
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
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition: num_vars.f90:89
hdf5_ops::print_hdf5_arrs
integer function, public print_hdf5_arrs(vars, PB3D_name, head_name, rich_lvl, disp_info, ind_print, remove_previous_arrs)
Prints a series of arrays, in the form of an array of pointers, to an HDF5 file.
Definition: HDF5_ops.f90:1132
num_utilities::lcm
recursive integer function, public lcm(u, v)
Returns common multiple using the Euclid's algorithm.
Definition: num_utilities.f90:2657
rich_vars::n_par_x
integer, public n_par_x
nr. of parallel points in field-aligned grid
Definition: rich_vars.f90:20
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
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
num_ops::calc_zero_hh
Finds the zero of a function using Householder iteration.
Definition: num_ops.f90:35
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
vac_vars::in_context
logical function, public in_context(ctxt)
Indicates whether current process is in the context.
Definition: vac_vars.f90:378
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
vac_ops::calc_gh_2
integer function calc_gh_2(vac, n_r_in, lims_r_in, x_vec_in, G_in, H_in, ext_in)
Calculates matrices G and H in axisymmetric configurations.
Definition: vac_ops.f90:983
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
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
helena_vars
Variables that have to do with HELENA quantities.
Definition: HELENA_vars.f90:4
num_vars::max_rvac_plot
real(dp), public max_rvac_plot
max. of R in which to plot vacuum
Definition: num_vars.f90:166
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
vac_utilities::calc_gh_int_1
subroutine, public calc_gh_int_1(G, H, x_s, x_in, norm_s, h_fac_in, step_size, tol)
Calculate G_ij and H_ij on an interval for field line aligned configurations.
Definition: vac_utilities.f90:33
grid_vars::max_alpha
real(dp), public max_alpha
max. of field-line label [ ] in field-aligned grid
Definition: grid_vars.f90:27
vac_utilities::mat_dis2loc
integer function, public mat_dis2loc(ctxt, mat_dis, lims_r, lims_c, mat_loc, proc)
Gathers a distributed vector on a single process.
Definition: vac_utilities.f90:277
vac_utilities::calc_gh_int_2
subroutine, public calc_gh_int_2(G, H, t_s, t_in, x_s, x_in, norm_s, norm_in, Aij, ql, tol, b_coeff, n_tor)
Calculate G_ij and H_ij on an interval for axisymmetric configurations.
Definition: vac_utilities.f90:99
num_vars::pb3d_name
character(len=max_str_ln), public pb3d_name
name of PB3D output file
Definition: num_vars.f90:139
vac_vars::vac_type
vacuum type
Definition: vac_vars.f90:46
num_vars::jump_to_sol
logical, public jump_to_sol
jump to solution
Definition: num_vars.f90:141
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
helena_vars::ias
integer, public ias
0 if top-bottom symmetric, 1 if not
Definition: HELENA_vars.f90:36
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
x_vars::n_mod_x
integer, public n_mod_x
size of m_X (pol. flux) or n_X (tor. flux)
Definition: X_vars.f90:129
num_vars::min_rvac_plot
real(dp), public min_rvac_plot
min. of R in which to plot vacuum
Definition: num_vars.f90:165
vac_ops::solve_phi_bem
integer function solve_phi_bem(vac, R, Phi, n_RPhi, n_loc_RPhi, lims_c_RPhi, desc_RPhi)
Calculates potential Phi on boundary of BEM in terms of some right-hand side .
Definition: vac_ops.f90:2163
num_ops::calc_zero_zhang
character(len=max_str_ln) function, public calc_zero_zhang(zero, fun, x_int_in)
Finds the zero of a function using Zhang's method, which is simpler than Brent's method.
Definition: num_ops.f90:462
num_vars::use_pol_flux_f
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
Definition: num_vars.f90:114
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
num_vars::eq_jobs_lims
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
Definition: num_vars.f90:77
grid_vars::min_alpha
real(dp), public min_alpha
min. of field-line label [ ] in field-aligned grid
Definition: grid_vars.f90:26
helena_vars::chi_h
real(dp), dimension(:), allocatable, public chi_h
poloidal angle
Definition: HELENA_vars.f90:23
grid_vars::n_alpha
integer, public n_alpha
nr. of field-lines
Definition: grid_vars.f90:23
vac_ops::calc_vac_res
integer function, public calc_vac_res(mds, vac)
Calculates the vacuum response.
Definition: vac_ops.f90:1417
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
num_vars::max_zeta_plot
real(dp), public max_zeta_plot
max. of zeta_plot
Definition: num_vars.f90:162
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
grid_utilities::nufft
integer function, public nufft(x, f, f_F, plot_name)
calculates the cosine and sine mode numbers of a function defined on a non-regular grid.
Definition: grid_utilities.f90:2901
num_vars::n_vac_plot
integer, dimension(2), public n_vac_plot
nr. of points in R and Z in vacuum
Definition: num_vars.f90:158
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
helena_vars::nchi
integer, public nchi
nr. of poloidal points
Definition: HELENA_vars.f90:35
vac_ops::print_output_vac
integer function, public print_output_vac(vac, data_name, rich_lvl)
Print vacuum quantities of a certain order to an output file.
Definition: vac_ops.f90:2347
mpi_utilities::broadcast_var
Wrapper function to broadcast a single variable using MPI.
Definition: MPI_utilities.f90:87
this documentation was created with Doxygen 1.8.17 by Toon Weyens