PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
SLEPC_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
5 !------------------------------------------------------------------------------!
9 !------------------------------------------------------------------------------!
10 module slepc_ops
11 #include <PB3D_macros.h>
12 #include <wrappers.h>
13 ! for slepc 3.9.0:
14 #include <slepc/finclude/slepceps.h>
15 ! for slepc 3.6.0:
16 !#include <slepc/finclude/slepcepsdef.h>
17 ! for slepc 3.5.3:
18 !#include <finclude/slepcepsdef.h>
19  use str_utilities
20  use output_ops
21  use messages
22  use slepceps
23  use num_vars, only: iu, dp, max_str_ln
24  use grid_vars, only: grid_type
25  use x_vars, only: x_1_type, x_2_type, modes_type
26  use vac_vars, only: vac_type
27  use sol_vars, only: sol_type
28 
29  implicit none
30  private
31  public solve_ev_system_slepc, start_slepc, setup_mats, set_bc, &
32  &setup_solver, setup_guess, get_solution, summarize_solution, &
33  &store_results, stop_slepc
34 #if ldebug
35  public debug_setup_mats, debug_set_bc, test_diff
36 #endif
37 
38  ! global variables
39  petscreal :: step_size ! step size in flux coordinates
40  petscint :: n_r ! n_r of solution
41  petscint :: loc_n_r ! loc_n_r of solution
42  integer, allocatable :: c_tot(:,:,:) ! total c corresponding to symmetric and asymmetric tensorial variables
43  logical :: use_hermitian = .false. ! use hermitian matrices (does not work well currently, see http://lists.mcs.anl.gov/pipermail/petsc-users/2016-October/030781.html)
44 #if ldebug
45  logical :: debug_setup_mats = .false.
46  logical :: debug_set_bc = .false.
47  logical :: test_diff = .false.
48  real(dp) :: diff_coeff
49 #endif
50 
51 contains
60  integer function solve_ev_system_slepc(mds,grid_sol,X,vac,sol) &
61  &result(ierr)
62 
63  use num_vars, only: matrix_slepc_style
64  use rich_vars, only: use_guess
65 #if ldebug
66  use num_vars, only: ltest
67  use input_utilities, only: get_real, get_log
68 #endif
69 
70  character(*), parameter :: rout_name = 'solve_EV_system_SLEPC'
71 
72  ! input / output
73  type(modes_type), intent(in) :: mds
74  type(grid_type), intent(in) :: grid_sol
75  type(x_2_type), intent(in) :: x
76  type(vac_type), intent(in) :: vac
77  type(sol_type), intent(inout) :: sol
78 
79  ! local variables
80  mat, target :: a ! matrix A in EV problem A X = lambda B X
81  mat, target :: b ! matrix B in EV problem A X = lambda B X
82  eps :: solver ! EV solver
83  petscint :: max_n_ev ! how many solutions saved
84  character(len=max_str_ln) :: err_msg ! error message
85 
86  ! initialize ierr
87  ierr = 0
88 
89  ! initialization
90  call writo('Initialize generalized Eigenvalue problem')
91 
92  ! start SLEPC
93  ierr = start_slepc(grid_sol%n(3))
94  chckerr('')
95 
96  ! set up step size
97  step_size = (grid_sol%r_F(grid_sol%n(3))-grid_sol%r_F(1))/&
98  &(grid_sol%n(3)-1) ! constant
99 
100  ! set up the matrix
101  call writo('Set up matrices')
102  call lvl_ud(1)
103  ierr = setup_mats(mds,grid_sol,x,a,b)
104  chckerr('')
105  call lvl_ud(-1)
106 
107 #if ldebug
108  if (debug_setup_mats) then
109  call writo('Testing if A and B are Hermitian by multiplying with &
110  &their Hermitian Transpose and subtracting 1')
111  call lvl_ud(1)
112 
113  ierr = test_mat_hermiticity(a,'A')
114  chckerr('')
115 
116  ierr = test_mat_hermiticity(b,'B')
117  chckerr('')
118 
119  call lvl_ud(-1)
120  end if
121 #endif
122 
123  ! set boundary conditions
124  call writo('Set up boundary conditions')
125  call lvl_ud(1)
126 
127  select case (matrix_slepc_style)
128  case (1) ! sparse
129  ierr = set_bc(mds,grid_sol,x,vac,a,b)
130  chckerr('')
131 #if ldebug
132  if (debug_set_bc) then
133  call writo('Testing if AFTER INSERTING BC''s, A and B &
134  &are Hermitian by multiplying with their &
135  &Hermitian Transpose and subtracting 1')
136  call lvl_ud(1)
137 
138  ierr = test_mat_hermiticity(a,'A_BC')
139  chckerr('')
140 
141  ierr = test_mat_hermiticity(b,'B_BC')
142  chckerr('')
143 
144  call lvl_ud(-1)
145  end if
146 #endif
147  case (2) ! shell
148  ierr = 1
149  err_msg = 'NOT YET IMPLEMENTED FOR SHELL MATRICES'
150  chckerr(err_msg)
151  end select
152 
153  call lvl_ud(-1)
154 
155  ! set up solver
156  call writo('Set up EV solver with defaults')
157  call lvl_ud(1)
158 #if ldebug
159  if (ltest) then
160  call writo('Test spectrum of A or B instead of solving &
161  &generalized Eigenvalue problem?')
162  if (get_log(.false.)) then
163  call writo('Spectrum of A (true) or B (false)?')
164  call lvl_ud(1)
165  if (get_log(.true.)) then ! A
166  call writo('Testing spectrum of A')
167  ierr = setup_solver(x,a,petsc_null_mat,solver)
168  chckerr('')
169  else ! B
170  call writo('Testing spectrum of B')
171  ierr = setup_solver(x,b,petsc_null_mat,solver)
172  chckerr('')
173  end if
174  call lvl_ud(-1)
175  else
176  ierr = setup_solver(x,a,b,solver)
177  chckerr('')
178  end if
179  else
180 #endif
181  ierr = setup_solver(x,a,b,solver)
182  chckerr('')
183 #if ldebug
184  end if
185 #endif
186  call lvl_ud(-1)
187 
188  ! set up guess
189  if (use_guess) then
190  call writo('Set up guess')
191  call lvl_ud(1)
192  ierr = setup_guess(sol,a,solver)
193  chckerr('')
194  call lvl_ud(-1)
195  end if
196 
197  ! get solution
198  call writo('Get solution')
199 
200  ierr = get_solution(solver)
201  chckerr('')
202 
203  ! summarize solution
204  call writo('Summarize solution')
205 
206  ierr = summarize_solution(solver,max_n_ev)
207  chckerr('')
208 
209  ! store results
210  call writo('Store results for '//trim(i2str(max_n_ev))//' least &
211  &stable Eigenvalues')
212 
213  ierr = store_results(mds,grid_sol,sol,solver,max_n_ev,a,b)
214  chckerr('')
215 
216  ! finalize
217  call writo('Finalize SLEPC')
218 
219  ! stop SLEPC
220  ierr = stop_slepc(a,b,solver)
221  chckerr('')
222 #if ldebug
223  contains
224  ! Test Hermiticity explicitely.
226  integer function test_mat_hermiticity(mat,mat_name) result(ierr)
227  use num_vars, only: data_dir
228  !use num_vars, only: rank
229 
230  character(*), parameter :: rout_name = 'test_mat_hermiticity'
231 
232  ! input / output
233  mat, intent(in) :: mat ! matrix to test
234  character(len=*) :: mat_name ! name of matrix
235 
236  ! local variables
237  mat :: mat_t ! Hermitian transpose of mat
238  mat :: mat_loc ! copy of mat, but without block storage (Hermitian transpose does not work for block)
239  petscscalar :: one = 1.0 ! one
240  petscviewer :: file_viewer ! viewer to write matrix to file
241  character(len=max_str_ln) :: file_name ! file name
242 
243  ! initialize ierr
244  ierr = 0
245 
246  ! duplicate mat into mat_loc
247  call matconvert(mat,matmpiaij,mat_initial_matrix,mat_loc,ierr)
248  chckerr('Failed to convert mat into mat_loc')
249 
250  ! test if mat is hermitian
251  call mathermitiantranspose(mat_loc,mat_initial_matrix,mat_t,ierr)
252  chckerr('Hermitian transpose of mat failed')
253  call mataxpy(mat_t,-one,mat_loc,same_nonzero_pattern,ierr)
254  chckerr('mat-mat_t failed')
255 
256  ! visualize the matrices
257  call petscoptionssetvalue(petsc_null_options,'-draw_pause','-1',&
258  &ierr)
259  chckerr('Failed to set option')
260  call writo(trim(mat_name)//' =')
261  call matview(mat,petsc_viewer_stdout_world,ierr)
262  !call MatView(mat,PETSC_VIEWER_DRAW_WORLD,ierr)
263  chckerr('Failed to view')
264  call writo(trim(mat_name)//' - '//trim(mat_name)//'^*T =')
265  call matview(mat_t,petsc_viewer_stdout_world,ierr)
266  chckerr('Failed to view')
267 
268  ! destroy matrices
269  call matdestroy(mat_t,ierr)
270  chckerr('Failed to destroy matrix mat_t')
271  call matdestroy(mat_loc,ierr)
272  chckerr('Failed to destroy matrix mat_loc')
273 
274  !! duplicate mat into mat_loc
275  !call MatDuplicate(mat,MAT_SHARE_NONZERO_PATTERN,mat_loc,ierr)
276  !CHCKERR('failed to duplicate mat into mat_loc')
277 
278  ! write real part to file
279  !file_name = trim(data_dir)//'/'//trim(mat_name)//'_RE'
280  file_name = trim(data_dir)//'/'//trim(mat_name)
281  call petscviewerasciiopen(petsc_comm_world,trim(file_name),&
282  &file_viewer,ierr)
283  chckerr('Unable to open file viewer')
284  !call PetscViewerSetFormat(file_viewer,PETSC_VIEWER_ASCII_DENSE,ierr)
285  call petscviewersetformat(file_viewer,petsc_viewer_ascii_matlab,ierr)
286  chckerr('Unable to set format')
287  !call MatView(mat_loc,file_viewer,ierr)
288  call matview(mat,file_viewer,ierr)
289  chckerr('Unable to write matrix to file')
290  call petscviewerdestroy(file_viewer,ierr)
291  chckerr('Unable to destroy file viewer')
292  !if (rank.eq.0) then
293  !call execute_command_line('tail -n +3 '//trim(file_name)//&
294  !&' > tempfile && mv tempfile '//trim(file_name),&
295  !&EXITSTAT=ierr)
296  !CHCKERR('Failed to execute command')
297  !end if
298 
299  !! write imaginary part to file
300  !file_name = trim(data_dir)//'/'//trim(mat_name)//'_IM'
301  !call MatCopy(mat,mat_loc,MAT_SHARE_NONZERO_PATTERN,ierr) ! mat_loc has same structure as mat
302  !CHCKERR('Failed to copy mat into mat_loc')
303  !call MatImaginaryPart(mat_loc,ierr)
304  !CHCKERR('Failed to get imaginary part')
305  !call PetscViewerASCIIOpen(PETSC_COMM_WORLD,trim(file_name),&
306  !&file_viewer,ierr)
307  !CHCKERR('Unable to open file viewer')
308  !call PetscViewerSetFormat(file_viewer,PETSC_VIEWER_ASCII_DENSE,ierr)
309  !CHCKERR('Unable to set format')
310  !call MatView(mat_loc,file_viewer,ierr)
311  !CHCKERR('Unable to write matrix to file')
312  !call PetscViewerDestroy(file_viewer,ierr)
313  !CHCKERR('Unable to destroy file viewer')
314  !if (rank.eq.0) then
315  !call execute_command_line('tail -n +3 '//trim(file_name)//&
316  !&' > tempfile && mv tempfile '//trim(file_name),&
317  !&EXITSTAT=ierr)
318  !CHCKERR('Failed to execute command')
319  !end if
320 
321  !! destroy matrices
322  !call MatDestroy(mat_loc,ierr)
323  !CHCKERR('Failed to destroy matrix mat_loc')
324  end function test_mat_hermiticity
325 #endif
326  end function solve_ev_system_slepc
327 
332  integer function start_slepc(n_r_sol) result(ierr)
333  use num_vars, only: n_procs, sol_n_procs
334  use files_ops, only: opt_args
335 
336  character(*), parameter :: rout_name = 'start_SLEPC'
337 
338  ! input / output
339  integer, intent(in) :: n_r_sol
340 
341  ! local variables
342  petscbool :: flg ! flag to catch options
343  petscint :: id ! counters
344  character(len=max_str_ln) :: option_name ! options
345  character(len=max_str_ln) :: err_msg ! error message
346 
347  ! initialize ierr
348  ierr = 0
349 
350  ! user message
351  call writo('initialize SLEPC...')
352  call lvl_ud(1)
353 
354  ! use MPI_Comm_world for PETSC_COMM_WORLD
355  petsc_comm_world = mpi_comm_world
356  if (n_procs.gt.n_r_sol) then ! too many processors
357  call writo('using too many processors: '//trim(i2str(n_procs))//&
358  &', while beyond '//trim(i2str(n_r_sol))//&
359  &' does not bring improvement',warning=.true.)
360  end if
361  call slepcinitialize(petsc_null_character,ierr) ! initialize SLEPC
362  chckerr('SLEPC failed to initialize')
363 
364  ! output
365  call writo('slepc started with '//trim(i2str(n_procs))&
366  &//' processes')
367  if (sol_n_procs.lt.n_procs) call writo('(but only '//&
368  &trim(i2str(sol_n_procs))//' processes will be used)')
369 
370  call lvl_ud(-1)
371 
372  ! checking for complex numbers
373  call writo('run tests...')
374 #if defined(PETSC_USE_COMPLEX)
375  ! OK
376 #else
377  err_msg = 'PETSC and SLEPC have to be configured and compiled &
378  &with the option "--with-scalar-type=complex"'
379  call slepcfinalize(ierr)
380  ierr = 1
381  chckerr(err_msg)
382 #endif
383 
384  ! catch options so they don't give a Petsc warning
385  do id = 1,size(opt_args)
386  option_name = opt_args(id)
387  if (trim(option_name).ne.'') then
388  call petscoptionshasname(petsc_null_options,&
389  &petsc_null_character,trim(option_name),flg,ierr) ! Petsc 3.7.6
390  !call PetscOptionsHasName(PETSC_NULL_CHARACTER,&
391  !&trim(option_name),flg,ierr) ! Petsc 3.6.4
392  err_msg = 'Failed to decide whether option has a name'
393  chckerr(err_msg)
394  end if
395  end do
396  end function start_slepc
397 
408  integer function setup_mats(mds,grid_sol,X,A,B) result(ierr)
409  use num_vars, only: ndps => norm_disc_prec_sol, matrix_slepc_style, &
410  &rank, sol_n_procs, bc_style, norm_disc_style_sol
411  use grid_utilities, only: trim_grid
412  use slepc_utilities, only: insert_block_mat
413  use x_vars, only: n_mod_x
414 #if ldebug
415  use num_vars, only: ltest
416  use input_utilities, only: get_real, get_log
417 #endif
418 
419  use num_utilities, only: c
420 
421  character(*), parameter :: rout_name = 'setup_mats'
422 
423  ! input / output
424  type(modes_type), intent(in) :: mds
425  type(grid_type), intent(in) :: grid_sol
426  type(x_2_type), intent(in), target :: x
427  mat, intent(inout) :: a
428  mat, intent(inout) :: b
429 
430  ! local variables
431  type(grid_type) :: grid_sol_trim ! trimmed solution grid
432  integer :: norm_id(2) ! untrimmed normal indices for trimmed solution grid
433  petscint :: kd, id ! counter
434  petscint :: i_lims(2) ! lower and upper limit of grid
435  petscint, allocatable :: d_nz(:) ! nr. of diagonal non-zeros
436  petscint, allocatable :: o_nz(:) ! nr. of off-diagonal non-zeros
437  character(len=max_str_ln) :: err_msg ! error message
438 
439  ! local variables also used in child routines
440  petscint :: bulk_i_lim(2) ! absolute limits of bulk matrix (excluding the BC's)
441 
442  ! initialize ierr
443  ierr = 0
444 
445  ! user output
446  select case (norm_disc_style_sol)
447  case (1) ! central finite differences
448  call writo('Normal discretization with central finite &
449  &differences of order '//trim(i2str(ndps))//&
450  &', stencil width '//trim(i2str(2*ndps+1))) ! ndps left and ndps right
451  case (2) ! left finite differences
452  call writo('Normal discretization with left finite &
453  &differences of order '//trim(i2str(ndps))//&
454  &', stencil width '//trim(i2str(ndps+1))) ! only ndps left
455  case default
456  err_msg = 'No solution normal discretization style associated &
457  &with '//trim(i2str(norm_disc_style_sol))
458  ierr = 1
459  chckerr(err_msg)
460  end select
461 
462 #if ldebug
463  if (ltest) then
464  call writo('Test introduction of numerical friction?')
465  if (get_log(.false.)) then
466  test_diff = .true.
467  call writo('The friction coefficient will be multiplied with &
468  &the central block, and inserted in the off-diagonal parts')
469  call writo('How large should it be?')
470  if (get_log(.false.)) diff_coeff = get_real(lim_lo=0._dp)
471  end if
472  end if
473 #endif
474 
475  ! check nr. of modes
476  if (n_mod_x.ne.x%n_mod(1) .or. n_mod_x.ne.x%n_mod(2)) then
477  ierr = 1
478  err_msg = 'Need square matrix of size [n_mod_X:n_mod_X]'
479  chckerr(err_msg)
480  end if
481 
482  ! setup total c if not yet done
483  if (.not.allocated(c_tot)) then
484  allocate(c_tot(n_mod_x,n_mod_x,2))
485  do kd = 1,n_mod_x
486  do id = 1,n_mod_x
487  c_tot(id,kd,1) = c([id,kd],.true.,n_mod_x)
488  c_tot(id,kd,2) = c([id,kd],.false.,n_mod_x)
489  end do
490  end do
491  end if
492 
493  ! trim grids
494  ierr = trim_grid(grid_sol,grid_sol_trim,norm_id)
495  chckerr('')
496 
497  ! set up loc_n_r and n_r
498  if (rank.lt.sol_n_procs) then
499  loc_n_r = grid_sol_trim%loc_n_r
500  else
501  loc_n_r = 0
502  end if
503  n_r = grid_sol_trim%n(3)
504 
505  ! set up bulk matrix absolute limits
506  ierr = set_bulk_lims(grid_sol_trim,bulk_i_lim)
507  chckerr('')
508  i_lims = [grid_sol_trim%i_min, grid_sol_trim%i_max]
509 
510  ! for BC_style 3 with symmetric finite differences, extend the grid
511  if (norm_disc_style_sol.eq.1 .and. bc_style(2).eq.3) then
512  n_r = n_r + ndps
513  if (rank.eq.sol_n_procs-1) then
514  loc_n_r = loc_n_r + ndps
515  i_lims(2) = i_lims(2) + ndps
516  end if
517  end if
518 
519  ! setup matrix A and B
520  select case (matrix_slepc_style)
521  case (1) ! sparse
522  ! calculate nonzeros
523  call set_nonzeros(i_lims)
524 
525  ! create matrix A
526  call matcreatebaij(petsc_comm_world,n_mod_x,n_mod_x*loc_n_r,&
527  &n_mod_x*loc_n_r,n_mod_x*n_r,n_mod_x*n_r,&
528  &petsc_null_integer,d_nz,petsc_null_integer,o_nz,a,ierr)
529  err_msg = 'MatCreateAIJ failed for matrix A'
530  chckerr(err_msg)
531 
532  ! deallocate tot_nz, d_nz and o_nz
533  deallocate(d_nz,o_nz)
534 
535  ! fill the matrix A
536  ierr = fill_mat(x%PV_0(1,:,norm_id(1):norm_id(2),:),&
537  &x%PV_1(1,:,norm_id(1):norm_id(2),:),&
538  &x%PV_2(1,:,norm_id(1):norm_id(2),:),bulk_i_lim,a)
539  chckerr('')
540 
541  call writo('matrix A set up:')
542  ierr = disp_mat_info(a)
543  chckerr('')
544 
545  ! duplicate A into B
546  call matduplicate(a,mat_share_nonzero_pattern,b,ierr) ! B has same structure as A
547  chckerr('failed to duplicate A into B')
548 
549  ! fill the matrix B
550  ierr = fill_mat(x%KV_0(1,:,norm_id(1):norm_id(2),:),&
551  &x%KV_1(1,:,norm_id(1):norm_id(2),:),&
552  &x%KV_2(1,:,norm_id(1):norm_id(2),:),bulk_i_lim,b)
553  chckerr('')
554 
555  call writo('matrix B set up:')
556  ierr = disp_mat_info(b)
557  chckerr('')
558  case (2) ! shell
559  ierr = 1
560  err_msg = 'NOT YET IMPLEMENTED FOR SHELL MATRICES'
561  chckerr(err_msg)
562  end select
563 
564  ! clean up
565  call grid_sol_trim%dealloc()
566  contains
567  ! Sets the limits of the indices of the bulk matrix, depending on the BC
568  ! style.
570  integer function set_bulk_lims(grid_sol,i_lim) result(ierr)
571  character(*), parameter :: rout_name = 'set_bulk_lims'
572 
573  ! input / output
574  type(grid_type), intent(in) :: grid_sol ! solution grid
575  integer, intent(inout) :: i_lim(2) ! min and max of bulk limits
576 
577  ! initialize ierr
578  ierr = 0
579 
580  ! set up i_min
581  select case (bc_style(1))
582  case (1)
583  i_lim(1) = max(grid_sol%i_min,ndps+1) ! first ndps rows write left BC's
584  case (2:4)
585  err_msg = 'Left BC''s cannot have BC type '//&
586  &trim(i2str(bc_style(1)))
587  ierr = 1
588  chckerr(err_msg)
589  case default
590  err_msg = 'No BC style associated with '//&
591  &trim(i2str(bc_style(1)))
592  ierr = 1
593  chckerr(err_msg)
594  end select
595 
596  ! set up i_max
597  ! Note: for BC_style(2) = 1, one could argue that we
598  ! should use n_r-ndps. This, however, generates an error for
599  ! norm_disc_style_sol = 2 where these elements are not written in
600  ! the main bulk assembly, so that petsc thinks there are no nonzeros
601  ! there. This is avoided by using n_r.
602  select case (bc_style(2))
603  case (1)
604  i_lim(2) = min(grid_sol%i_max,n_r) ! last ndps rows write right BC's, will be overwritten
605  case (2:4)
606  i_lim(2) = min(grid_sol%i_max,n_r) ! right BC is added later
607  case default
608  err_msg = 'No BC style associated with '//&
609  &trim(i2str(bc_style(2)))
610  ierr = 1
611  chckerr(err_msg)
612  end select
613  end function set_bulk_lims
614 
615  ! Fills a matrix according to norm_disc_prec_sol [ADD REF]:
616  ! 1. first order accuracy for all terms
617  ! 2. higher order accuracy for internal first order derivatives
618  ! It is used for both the matrix A and B, corresponding to the plasma
619  ! potential energy and the (perpendicular) kinetic energy.
620  !
621  ! Makes use of mds, n_r and grid_sol_trim
623  integer function fill_mat(V_0,V_1,V_2,bulk_i_lim,mat) result(ierr)
624  use num_utilities, only: con, calc_coeff_fin_diff
625 
626  character(*), parameter :: rout_name = 'fill_mat'
627 
628  ! input / output
629  petscscalar, intent(in) :: v_0(:,:,:) ! either PV_int_0 or KV_int_0 in equilibrium normal grid
630  petscscalar, intent(in) :: v_1(:,:,:) ! either PV_int_1 or KV_int_1 in equilibrium normal grid
631  petscscalar, intent(in) :: v_2(:,:,:) ! either PV_int_2 or KV_int_2 in equilibrium normal grid
632  petscint, intent(in) :: bulk_i_lim(2) ! absolute limits of bulk matrix (excluding the BC's)
633  mat, intent(inout) :: mat ! either A or B
634 
635  ! local variables (not to be used in child routines)
636  character(len=max_str_ln) :: err_msg ! error message
637  petscscalar, allocatable :: loc_block(:,:) ! [n_mod_X:n_mod_X] block matrix for 1 normal point
638  petscreal, allocatable :: ndc(:) ! normal discretization coefficients
639  petscint :: id, jd, kd ! counters
640  petscint :: kd_loc ! kd in local variables
641  petscint :: k, m ! counters
642  petscint :: ndc_ind ! index for finite differences
643  petscint :: st_size ! stencil size
644 #if ldebug
645  petscscalar, allocatable :: loc_block_0_backup(:,:) ! backed up V_0 local block
646 #endif
647 
648  ! for tests
649  petscint :: r_sol_start, r_sol_end ! start block and end block
650 
651  ! initialize ierr
652  ierr = 0
653 
654  ! test whether the matrix range coincides with i_min and i_max
655  call matgetownershiprange(mat,r_sol_start,r_sol_end,ierr) ! starting and ending row r_sol_start and r_sol_end
656  err_msg = 'Couldn''t get ownership range of matrix'
657  chckerr(err_msg)
658  r_sol_start = r_sol_start/n_mod_x ! count per block
659  r_sol_end = r_sol_end/n_mod_x ! count per block
660  if (norm_disc_style_sol.eq.1 .and. bc_style(2).eq.3) then
661  if (rank.eq.sol_n_procs-1) r_sol_end = r_sol_end - ndps
662  end if
663  if (rank.lt.sol_n_procs) then
664  if (grid_sol_trim%i_min.ne.r_sol_start+1) then
665  ierr = 1
666  err_msg = 'start of matrix in this process does not &
667  &coincide with tabulated value'
668  chckerr(err_msg)
669  end if
670  if (grid_sol_trim%i_max.ne.r_sol_end) then
671  ierr = 1
672  err_msg = 'end of matrix in this process does not &
673  &coincide with tabulated value'
674  chckerr(err_msg)
675  end if
676  end if
677 
678  ! allocate local block
679  allocate(loc_block(n_mod_x,n_mod_x))
680 
681  ! iterate over all rows of this rank
682  do kd = bulk_i_lim(1)-1,bulk_i_lim(2)-1 ! (indices start with 0 here)
683  ! set up local kd
684  kd_loc = kd+2-grid_sol_trim%i_min
685 
686  ! set up ndc
687  select case (norm_disc_style_sol)
688  case (1) ! central differences
689  if ((bc_style(2).eq.2 .or. bc_style(2).eq.4) &
690  &.and. kd.gt.n_r-1-ndps) then ! asymmetric finite differences: shifted indices
691  ndc_ind = 1+2*ndps - (n_r-1-kd)
692  else
693  ndc_ind = 1+ndps
694  end if
695  st_size = 2*ndps
696  case (2) ! left differences
697  ndc_ind = 1+ndps
698  st_size = ndps
699  end select
700  ierr = calc_coeff_fin_diff(1,st_size+1,ndc_ind,ndc)
701  chckerr('')
702  ndc = ndc/step_size ! scale by step size
703 
704  ! fill the blocks
705 
706  ! -------------!
707  ! BLOCKS ~ V_0 !
708  ! -------------!
709  ! fill local block
710  do m = 1,n_mod_x
711  do k = 1,n_mod_x
712  loc_block(k,m) = &
713  &con(sum(v_0(:,kd_loc,c_tot(k,m,1))),&
714  &[k,m],.true.) ! symmetric matrices need con()
715  end do
716  end do
717  loc_block = loc_block*step_size ! include step size of normal integral
718 
719 #if ldebug
720  if (test_diff) then
721  ! back up the V_0 block
722  allocate(loc_block_0_backup(n_mod_x,n_mod_x))
723  loc_block_0_backup = loc_block
724  end if
725 #endif
726 
727  ! add block to kd + (0,0)
728  ierr = insert_block_mat(mds,loc_block,mat,kd,[0,0],n_r)
729  chckerr('')
730 
731  ! -------------!
732  ! BLOCKS ~ V_1 !
733  ! -------------!
734  ! fill local block
735  do m = 1,n_mod_x
736  do k = 1,n_mod_x
737  loc_block(k,m) = sum(v_1(:,kd_loc,c_tot(k,m,2))) ! asymetric matrices don't need con()
738  end do
739  end do
740  loc_block = loc_block*step_size ! include step size of normal integral
741 
742  ! add block to kd +
743  ! (0,1-ndc_ind:1-st_size+jd) + Hermitian conjugate
744  do jd = 1,1+st_size
745  ierr = insert_block_mat(mds,loc_block*ndc(jd),mat,&
746  &kd,[0,jd-ndc_ind],n_r,transp=.true.) ! so that jd = ndc_ind falls on [0,0]
747  chckerr('')
748  end do
749 
750  ! -------------!
751  ! BLOCKS ~ V_2 !
752  ! -------------!
753  ! fill local block
754  do m = 1,n_mod_x
755  do k = 1,n_mod_x
756  loc_block(k,m) = &
757  &con(sum(v_2(:,kd_loc,c_tot(k,m,1))),[k,m],.true.) ! symmetric matrices need con()
758  end do
759  end do
760  loc_block = loc_block*step_size ! include step size of normal integral
761 
762 #if ldebug
763  if (test_diff) then
764  loc_block = loc_block + diff_coeff*loc_block_0_backup
765  deallocate(loc_block_0_backup)
766  end if
767 #endif
768 
769  ! add block to kd +
770  ! (1-ndc_ind:1-st_size+id,1-ndc_ind:1-st_size+jd)
771  do jd = 1,1+st_size
772  do id = 1,1+st_size
773  ierr = insert_block_mat(mds,loc_block*ndc(id)*ndc(jd),&
774  &mat,kd,[id-ndc_ind,jd-ndc_ind],n_r) ! so that id,jd = ndc_ind,ndc_ind falls on [0,0]
775  chckerr('')
776  end do
777  end do
778  end do
779 
780  ! assemble the matrix
781  call matassemblybegin(mat,mat_final_assembly,ierr)
782  chckerr('Coulnd''t begin assembly of matrix')
783  call matassemblyend(mat,mat_final_assembly,ierr)
784  chckerr('Coulnd''t end assembly of matrix')
785 
786  ! Hermitian matrix
787  if (use_hermitian) then
788  call matsetoption(mat,mat_hermitian,petsc_true,ierr)
789  chckerr('Coulnd''t set option Hermitian')
790  end if
791 
792  ! clean up
793  deallocate(loc_block)
794  end function fill_mat
795 
796  ! Display information about matrix.
798  integer function disp_mat_info(mat) result(ierr)
799  character(*), parameter :: rout_name = 'disp_mat_info'
800 
801  ! input / output
802  mat, intent(inout) :: mat ! either A or B
803 
804  ! local variables
805  real(dp) :: mat_info(mat_info_size) ! information about matrix
806 
807  ! initialize ierr
808  ierr = 0
809 
810  call lvl_ud(1)
811 
812  select case (matrix_slepc_style)
813  case (1) ! sparse
814  call matgetinfo(mat,mat_global_sum,mat_info,ierr)
815  chckerr('')
816  call writo('memory usage: '//&
817  &trim(r2strt(mat_info(mat_info_memory)*1.e-6_dp))//&
818  &' MB')
819  call writo('nonzero''s allocated: '//&
820  &trim(r2strt(mat_info(mat_info_nz_allocated))))
821  if (mat_info(mat_info_nz_unneeded).gt.0._dp) then
822  call writo('of which unused: '//&
823  &trim(r2strt(mat_info(mat_info_nz_unneeded))))
824  end if
825  case (2) ! shell
826  call writo('shell matrix')
827  end select
828 
829  call lvl_ud(-1)
830  end function disp_mat_info
831 
832  ! Sets nonzero elements d_nz and o_nz.
833  ! Makes use of ndps.
835  subroutine set_nonzeros(i_lims)
836  ! input / output
837  petscint, intent(in) :: i_lims(2) ! lower and upper limit of grid
838 
839  ! local variables
840  petscint, allocatable :: col_lims(:,:) ! column limits of nonzeros
841  petscint, allocatable :: tot_nz(:) ! nr. of total non-zeros
842  petscint :: st_size(2) ! half stencil size left and right
843 
844  ! initialize the numbers of non-zeros in diagonal and off-diagonal
845  allocate(tot_nz(loc_n_r)); tot_nz = 0
846  allocate(d_nz(loc_n_r)); d_nz = 0
847  allocate(o_nz(loc_n_r)); o_nz = 0
848 
849  ! initialize helper variable
850  allocate(col_lims(loc_n_r,2)); col_lims = 0
851 
852  if (rank.lt.sol_n_procs) then
853  ! set (half) stencil size
854  select case (norm_disc_style_sol)
855  case (1) ! central finite differences
856  st_size = [ndps,ndps]
857  case (2) ! left finite differences
858  st_size = [ndps,0]
859  end select
860 
861  ! Set column limits of nonzero columns per row
862  ! At row i the largest nonzero column is the one due to the
863  ! row at i+st_size(1), as this one will work back at most
864  ! st_size(1) places. The last column that this row will
865  ! contribute to is st_size(2) further, which makes for i +
866  ! sum(st_size). For the lower bound a similar argument can be
867  ! constructed.
868  do kd = 1,loc_n_r
869  col_lims(kd,1) = max(kd+i_lims(1)-1-sum(st_size),1)
870  col_lims(kd,2) = min(kd+i_lims(1)-1+sum(st_size),n_r)
871  end do
872 
873  ! calculate number of total nonzero entries:
874  ! Number of nonzeros columns.
875  tot_nz = col_lims(:,2) - col_lims(:,1) + 1
876 
877  ! calculate number of nonzero off-diagonal entries
878  ! For row i, this is equal to the number of elements with
879  ! column index greater than the maximum row index (i_lims(2))
880  ! plus the number of elements with column index lower than the
881  ! minimum row index (i_lims(1)).
882  do kd = 1,loc_n_r
883  o_nz(kd) = max(0,i_lims(1)-col_lims(kd,1)) + &
884  &max(0,col_lims(kd,2)-i_lims(2))
885  end do
886 
887  ! calculate number of nonzero diagonal entries
888  d_nz = tot_nz-o_nz
889  else
890  d_nz = 0
891  o_nz = 0
892  end if
893 
894  ! clean up
895  deallocate(tot_nz)
896  end subroutine set_nonzeros
897  end function setup_mats
898 
932  integer function set_bc(mds,grid_sol,X,vac,A,B) result(ierr)
933  use num_vars, only: ndps => norm_disc_prec_sol, bc_style, ev_bc, &
934  &norm_disc_style_sol
935  use x_vars, only: n_mod_x
936  use mpi_utilities, only: get_ser_var, wait_mpi
937  use num_utilities, only: con, c
938  use slepc_utilities, only: insert_block_mat
939 
940  character(*), parameter :: rout_name = 'set_BC'
941 
942  ! input / output
943  type(modes_type), intent(in) :: mds
944  type(grid_type), intent(in) :: grid_sol
945  type(x_2_type), intent(in) :: x
946  type(vac_type), intent(in) :: vac
947  mat, intent(inout) :: a
948  mat, intent(inout) :: b
949 
950  ! local variables
951  petscint :: kd ! counter
952  petscint :: n_min, n_max ! absolute limits excluding the BC's
953  character(len=max_str_ln) :: err_msg ! error message
954 
955  ! initialize ierr
956  ierr = 0
957 
958  call writo('Preparing variables')
959  call lvl_ud(1)
960 
961  ! check nr. of modes
962  if (n_mod_x.ne.x%n_mod(1) .or. n_mod_x.ne.x%n_mod(2)) then
963  ierr = 1
964  err_msg = 'Need square matrix of size [n_mod_X:n_mod_X]'
965  chckerr(err_msg)
966  end if
967 
968  call lvl_ud(-1)
969 
970  ! set up n_min, depending on BC style
971  select case (bc_style(1))
972  case (1)
973  n_min = ndps ! dirichlet BC requires stencil
974  case (2:4)
975  err_msg = 'Left BC''s cannot have BC type '//&
976  &trim(i2str(bc_style(1)))
977  ierr = 1
978  chckerr(err_msg)
979  case default
980  err_msg = 'No BC style associated with '//&
981  &trim(i2str(bc_style(1)))
982  ierr = 1
983  chckerr(err_msg)
984  end select
985 
986  ! set up n_max, depending on BC style
987  select case (bc_style(2))
988  case (1)
989  select case (norm_disc_style_sol)
990  case (1) ! central differences
991  n_max = ndps ! dirichlet BC requires stencil
992  case (2) ! left differences
993  n_max = 1 ! dirichlet BC only on last point
994  end select
995  case (2)
996  n_max = 1 ! only 1 element carries vacuum contribution
997  case (3)
998  n_max = 1 ! only 1 element carries vacuum contribution
999  case (4)
1000  n_max = 1 ! only last element carries BC
1001  case default
1002  err_msg = 'No BC style associated with '//&
1003  &trim(i2str(bc_style(1)))
1004  ierr = 1
1005  chckerr(err_msg)
1006  end select
1007 
1008  call writo('Setting BC deep within plasma')
1009  call lvl_ud(1)
1010  call writo('Using artificial eigenvalue EV_BC = '//trim(r2strt(ev_bc)))
1011 
1012  ! iterate over all positions where to set left BC
1013  do kd = 1,n_min
1014  if (grid_sol%i_min.le.kd .and. grid_sol%i_max.ge.kd) then ! decide which process does the setting
1015  select case (bc_style(1))
1016  case (1)
1017  ierr = set_bc_1(kd-1,a,b,.false.) ! indices start at 0
1018  chckerr('')
1019  case (2:4)
1020  err_msg = 'Left BC''s cannot have BC type '//&
1021  &trim(i2str(bc_style(1)))
1022  ierr = 1
1023  chckerr(err_msg)
1024  case default
1025  err_msg = 'No BC style associated with '//&
1026  &trim(i2str(bc_style(1)))
1027  ierr = 1
1028  chckerr(err_msg)
1029  end select
1030  end if
1031  end do
1032 
1033  ! assemble the matrices (cannot mix overwrite and not)
1034  call matassemblybegin(a,mat_final_assembly,ierr)
1035  chckerr('Coulnd''t begin assembly of matrix')
1036  call matassemblybegin(b,mat_final_assembly,ierr)
1037  chckerr('Coulnd''t begin assembly of matrix')
1038  call matassemblyend(a,mat_final_assembly,ierr)
1039  chckerr('Coulnd''t end assembly of matrix')
1040  call matassemblyend(b,mat_final_assembly,ierr)
1041  chckerr('Coulnd''t end assembly of matrix')
1042 
1043  call lvl_ud(-1)
1044 
1045  ! wait to get messages consistent
1046  ierr = wait_mpi()
1047  chckerr('')
1048 
1049  call writo('Setting BC at plasma edge')
1050  call lvl_ud(1)
1051 
1052  ! iterate over all positions where to set right BC
1053  do kd = grid_sol%n(3),grid_sol%n(3)-n_max+1,-1
1054  if (grid_sol%i_min.le.kd .and. grid_sol%i_max.ge.kd) then ! decide which process does the setting
1055  select case (bc_style(2))
1056  case (1)
1057  ierr = set_bc_1(kd-1,a,b,.true.) ! indices start at 0
1058  chckerr('')
1059  case (2)
1060  ierr = set_bc_2(kd-1,a) ! indices start at 0
1061  chckerr('')
1062  case (3)
1063  ierr = set_bc_3(kd-1,a) ! indices start at 0
1064  chckerr('')
1065  case (4)
1066  ierr = set_bc_4(kd-1,kd-grid_sol%i_min+1,x,a,b,&
1067  &grid_sol%n(3)) ! indices start at 0
1068  chckerr('')
1069  case default
1070  err_msg = 'No BC style associated with '//&
1071  &trim(i2str(bc_style(2)))
1072  ierr = 1
1073  chckerr(err_msg)
1074  end select
1075  end if
1076  end do
1077 
1078  ! assemble the matrices
1079  call matassemblybegin(a,mat_final_assembly,ierr)
1080  chckerr('Coulnd''t begin assembly of matrix')
1081  call matassemblybegin(b,mat_final_assembly,ierr)
1082  chckerr('Coulnd''t begin assembly of matrix')
1083  call matassemblyend(a,mat_final_assembly,ierr)
1084  chckerr('Coulnd''t end assembly of matrix')
1085  call matassemblyend(b,mat_final_assembly,ierr)
1086  chckerr('Coulnd''t end assembly of matrix')
1087 
1088  call lvl_ud(-1)
1089  contains
1090  ! set BC style 1:
1091  ! Set to zero
1092  ! A(ind,ind) = EV_BC, B(ind,ind) = 1,
1093  ! A(ind+1..ind+p,ind+1..ind+p) = 0, B(ind+1..ind+p,ind+1..ind+p) = 0,
1094  ! where ind indicates the row where the BC is centered.
1096  integer function set_bc_1(r_id,A,B,BC_right) result(ierr)
1097  character(*), parameter :: rout_name = 'set_BC_1'
1098 
1099  ! input / output
1100  integer, intent(in) :: r_id ! position at which to set BC
1101  mat, intent(inout) :: a, b ! Matrices A and B from A X = lambda B X
1102  logical :: bc_right ! if BC is at right (so there are vacuum terms)
1103 
1104  ! local variables
1105  integer :: kd ! counter
1106  petscscalar, allocatable :: loc_block(:,:) ! [n_mod_X:n_mod_X] block matrix for 1 normal point
1107  petscint :: pmone ! plus or minus one
1108  petscint :: st_size ! stencil size
1109 
1110  ! initialize ierr
1111  ierr = 0
1112 
1113  ! user output
1114  call writo('Boundary style at row '//trim(i2str(r_id+1))//&
1115  &': Eigenvector set to zero',persistent=.true.)
1116 
1117  ! initialize local blocks
1118  allocate(loc_block(n_mod_x,n_mod_x))
1119  loc_block = 0.0_dp
1120  do kd = 1,n_mod_x
1121  loc_block(kd,kd) = 1.0_dp
1122  end do
1123 
1124  ! set up plus minus one
1125  if (bc_right) then
1126  pmone = 1
1127  else
1128  pmone = -1
1129  end if
1130 
1131  ! set up stencil size
1132  select case (norm_disc_style_sol)
1133  case (1) ! central differences
1134  st_size = 2*ndps
1135  case (2) ! left differences
1136  st_size = ndps
1137  end select
1138 
1139  ! set block r_id + (0,0)
1140  ierr = insert_block_mat(mds,ev_bc*loc_block,a,r_id,[0,0],n_r,&
1141  &overwrite=.true.,ind_insert=.true.)
1142  chckerr('')
1143  ierr = insert_block_mat(mds,loc_block,b,r_id,[0,0],n_r,&
1144  &overwrite=.true.,ind_insert=.true.)
1145  chckerr('')
1146 
1147  ! iterate over range 2
1148  do kd = 1,st_size
1149  ! set block r_id +/- (0,kd) and Hermitian conjugate
1150  ierr = insert_block_mat(mds,0*loc_block,a,r_id,[0,-pmone*kd],&
1151  &n_r,overwrite=.true.,transp=.true.,ind_insert=.true.)
1152  chckerr('')
1153  ierr = insert_block_mat(mds,0*loc_block,b,r_id,[0,-pmone*kd],&
1154  &n_r,overwrite=.true.,transp=.true.,ind_insert=.true.)
1155  chckerr('')
1156  end do
1157 
1158  ! clean up
1159  deallocate(loc_block)
1160  end function set_bc_1
1161 
1162  ! set BC style 2:
1163  ! Minimization of surface energy through asymmetric fin. differences
1165  integer function set_bc_2(r_id,A) result(ierr)
1166  character(*), parameter :: rout_name = 'set_BC_2'
1167 
1168  ! input / output
1169  integer, intent(in) :: r_id ! position at which to set BC
1170  mat, intent(inout) :: a ! Matrices A from A X = lambda B X
1171 
1172  ! initialize ierr
1173  ierr = 0
1174 
1175  ! user output
1176  call writo('Boundary style at row '//trim(i2str(r_id+1))//&
1177  &': Minimization of surface energy through asymmetric finite &
1178  & differences',persistent=.true.)
1179 
1180  ! -------------!
1181  ! BLOCKS ~ vac !
1182  ! -------------!
1183  ! add block to r_id + (0,0)
1184  ierr = 2
1185  chckerr('Vacuum has not been implemented yet!')
1186  ierr = insert_block_mat(mds,vac%res,a,r_id,[0,0],n_r,&
1187  &ind_insert=.true.)
1188  chckerr('')
1189  end function set_bc_2
1190 
1191  ! set BC style 3:
1192  ! Minimization of surface energy through extension of grid
1194  integer function set_bc_3(r_id,A) result(ierr)
1195  character(*), parameter :: rout_name = 'set_BC_3'
1196 
1197  ! input / output
1198  integer, intent(in) :: r_id ! position at which to set BC
1199  mat, intent(inout) :: a ! Matrices A from A X = lambda B X
1200 
1201  ! initialize ierr
1202  ierr = 0
1203 
1204  ! user output
1205  call writo('Boundary style at row '//trim(i2str(r_id+1))//&
1206  &': Minimization of surface energy through extension of grid',&
1207  &persistent=.true.)
1208 
1209  ! -------------!
1210  ! BLOCKS ~ vac !
1211  ! -------------!
1212  ! add block to r_id + (0,0)
1213  ierr = 2
1214  chckerr('Vacuum has not been implemented yet!')
1215  ierr = insert_block_mat(mds,vac%res,a,r_id,[0,0],n_r,&
1216  &ind_insert=.true.)
1217  chckerr('')
1218  end function set_bc_3
1219 
1220  ! set BC style 4:
1221  ! Explicit introduction of the surface energy minimization
1222  ! V1^T X + V2 X' + delta_vac X = 0 at surface
1223  !
1224  ! Makes use of mds.
1226  integer function set_bc_4(r_id,r_id_loc,X,A,B,n_r) result(ierr)
1227  use num_vars, only: norm_disc_style_sol
1228  use num_utilities, only: calc_coeff_fin_diff
1229 
1230  character(*), parameter :: rout_name = 'set_BC_4'
1231 
1232  ! input / output
1233  integer, intent(in) :: r_id ! global position at which to set BC (starting at 0)
1234  integer, intent(in) :: r_id_loc ! local index in perturbation tables
1235  type(x_2_type), intent(in) :: x ! field-averaged perturbation variables (so only first index)
1236  mat, intent(inout) :: a, b ! Matrices A and B from A X = lambda B X
1237  integer, intent(in) :: n_r ! number of grid points of solution grid
1238 
1239  ! local variables
1240  petscint :: jd ! counter
1241  petscint :: k, m ! counters
1242  petscint :: st_size ! stencil size
1243  petscreal, allocatable :: ndc(:) ! normal discretization coefficients
1244  petscscalar, allocatable :: loc_block(:,:,:) ! [n_mod_X:n_mod_X] block matrix for 1 normal point
1245 
1246  ! initialize ierr
1247  ierr = 0
1248 
1249  ! user output
1250  call writo('Boundary style at row '//trim(i2str(r_id+1))//&
1251  &': Explicit introduction of surface energy minimization',&
1252  &persistent=.true.)
1253 
1254  ! initialize local blocks
1255  allocate(loc_block(n_mod_x,n_mod_x,2))
1256 
1257  ! set up ndc
1258  select case (norm_disc_style_sol)
1259  case (1) ! central differences
1260  st_size = 2*ndps
1261  case (2) ! left differences
1262  st_size = ndps
1263  end select
1264  ierr = calc_coeff_fin_diff(1,st_size+1,st_size+1,ndc)
1265  chckerr('')
1266  ndc = ndc/step_size ! scale by step size
1267 
1268  ! ----------------------!
1269  ! BLOCKS ~ V1^T and vac !
1270  ! ----------------------!
1271  ! fill local blocks for A and B
1272  do m = 1,n_mod_x
1273  do k = 1,n_mod_x
1274  loc_block(k,m,1) = conjg(sum(&
1275  &x%PV_1(1,:,r_id_loc,c([m,k],.false.,n_mod_x)))) ! asymetric matrices don't need con()
1276  loc_block(k,m,2) = conjg(sum(&
1277  &x%KV_1(1,:,r_id_loc,c([m,k],.false.,n_mod_x)))) ! asymetric matrices don't need con()
1278  end do
1279  end do
1280  loc_block(:,:,1) = loc_block(:,:,1) + vac%res
1281 
1282  ! set block r_id + (0,0)
1283  ierr = insert_block_mat(mds,loc_block(:,:,1)*ndc(1+st_size),a,r_id,&
1284  &[0,0],n_r,overwrite=.true.,ind_insert=.true.)
1285  chckerr('')
1286  ! set block r_id + (0,0)
1287  ierr = insert_block_mat(mds,loc_block(:,:,2)*ndc(1+st_size),b,r_id,&
1288  &[0,0],n_r,overwrite=.true.,ind_insert=.true.)
1289  chckerr('')
1290 
1291  ! -------------!
1292  ! BLOCKS ~ V_2 !
1293  ! -------------!
1294  ! fill local blocks for A and B
1295  do m = 1,n_mod_x
1296  do k = 1,n_mod_x
1297  loc_block(k,m,1) = &
1298  &con(sum(x%PV_2(1,:,r_id_loc,c([k,m],.true.,n_mod_x))),&
1299  &[k,m],.true.) ! symmetric matrices need con()
1300  loc_block(k,m,2) = &
1301  &con(sum(x%KV_2(1,:,r_id_loc,c([k,m],.true.,n_mod_x))),&
1302  &[k,m],.true.) ! symmetric matrices need con()
1303  end do
1304  end do
1305 
1306  ! set block r_id + (0,-st_size:0)
1307  do jd = -st_size,0
1308  ierr = insert_block_mat(mds,loc_block(:,:,1)*&
1309  &ndc(jd+st_size+1)*ndc(1+st_size),a,&
1310  &r_id,[0,jd],n_r,overwrite=.true.,&
1311  &ind_insert=.true.)
1312  chckerr('')
1313  ierr = insert_block_mat(mds,loc_block(:,:,2)*&
1314  &ndc(jd+st_size+1)*ndc(1+st_size),b,&
1315  &r_id,[0,jd],n_r,overwrite=.true.,&
1316  &ind_insert=.true.)
1317  chckerr('')
1318  end do
1319  end function set_bc_4
1320  end function set_bc
1321 
1325  integer function setup_solver(X,A,B,solver) result(ierr)
1326  use num_vars, only: n_sol_requested, tol_slepc, max_it_slepc, &
1327  &ev_guess, solver_slepc_style
1328  use rich_vars, only: rich_lvl
1329  use x_vars, only: n_mod_x
1330 
1331  character(*), parameter :: rout_name = 'setup_solver'
1332 
1333  ! input / output
1334  type(x_2_type), intent(in) :: x
1335  mat, intent(in) :: a, b
1336  eps, intent(inout) :: solver
1337 
1338  ! local variables
1339  petscscalar :: one = 1.0 ! one
1340  st :: solver_st ! spectral transformation associated to the solver
1341  petscint :: n_sol ! how many solutions can be requested (normally n_sol_requested)
1342  character(len=max_str_ln) :: err_msg ! error message
1343 
1344  ! initialize ierr
1345  ierr = 0
1346 
1347  ! check nr. of modes
1348  if (n_mod_x.ne.x%n_mod(1) .or. n_mod_x.ne.x%n_mod(2)) then
1349  ierr = 1
1350  err_msg = 'Need square matrix of size [n_mod_X:n_mod_X]'
1351  chckerr(err_msg)
1352  end if
1353 
1354  !call PetscOptionsSetValue(PETSC_NULL_OPTION,'-eps_view','-1',ierr)
1355  call epscreate(petsc_comm_world,solver,ierr)
1356  chckerr('EPSCreate failed')
1357 
1358  call epssetoperators(solver,a,b,ierr) ! generalized EV problem A X = lambda B X
1359  chckerr('EPSetOperators failed')
1360 
1361  if (use_hermitian) then
1362  call writo('Set problem type to Generalized Hermitian')
1363  call epssetproblemtype(solver,eps_ghep,ierr)
1364  chckerr('Set problem type failed')
1365  end if
1366 
1367  ! set some options depending on which default style
1368  select case (solver_slepc_style)
1369  case (1) ! Krylov-Schur with shift-invert
1370  ! Krylov-Schur is the most robust
1371  call writo('set algorithm to Krylov-Schur')
1372  call epssettype(solver,epskrylovschur,ierr)
1373  err_msg = 'Failed to set type to Krylov-Schur'
1374  chckerr(err_msg)
1375 
1376  ! get spectral transformation
1377  call epsgetst(solver,solver_st,ierr)
1378  err_msg = 'Failed to get ST from solver'
1379  chckerr(err_msg)
1380 
1381  ! shift-invert is the best method for generalized EV problems,
1382  ! because a matrix needs to be inverted anyway
1383  call writo('use shift-invert')
1384  call stsettype(solver_st,stsinvert,ierr)
1385  err_msg = 'Failed to set type to STSINVERT'
1386  chckerr(err_msg)
1387  case (2) ! generalized Davidson with preconditioner
1388  ! GD might be faster than Krylov-Schur
1389  call writo('set algorithm to Generalized Davidson')
1390  call epssettype(solver,epsgd,ierr)
1391  err_msg = 'Failed to set type to Generalized Davidson'
1392  chckerr(err_msg)
1393 
1394  ! warning
1395  call lvl_ud(1)
1396  call writo('This is often less robust than the Krylov-Schur &
1397  &method.',alert=.true.)
1398  call lvl_ud(1)
1399  call writo('If it converges to a strange result, try using &
1400  &that')
1401  call lvl_ud(-1)
1402  call lvl_ud(-1)
1403 
1404  ! get spectral transformation
1405  call epsgetst(solver,solver_st,ierr)
1406  err_msg = 'Failed to get ST from solver'
1407  chckerr(err_msg)
1408 
1409  ! GD works only with a preconditioner
1410  call writo('Use preconditioner')
1411  call stsettype(solver_st,stprecond,ierr)
1412  err_msg = 'Failed to set type to STPRECOND'
1413  chckerr(err_msg)
1414  case default
1415  call writo('SLEPC solver style '//&
1416  &trim(i2str(solver_slepc_style))//' not recognized',&
1417  &alert=.true.)
1418  end select
1419 
1420  ! set the guess as target
1421  call writo('set the target eigenvalue to '//trim(r2str(ev_guess)))
1422  call epssettarget(solver,ev_guess*one,ierr)
1423  err_msg = 'Failed to set the target of solver'
1424  chckerr(err_msg)
1425 
1426  ! request n_sol_requested Eigenpairs
1427  if (n_sol_requested.gt.n_r*n_mod_x) then
1428  call writo('max. nr. of solutions requested capped to problem &
1429  &dimension ('//trim(i2str(n_r*n_mod_x))//')',&
1430  &warning=.true.)
1431  call writo('Increase either n_r or number of pol. modes or &
1432  &decrease n_sol_requested')
1433  n_sol = n_r*n_mod_x
1434  else
1435  n_sol = n_sol_requested
1436  end if
1437  call epssetdimensions(solver,n_sol,petsc_decide,&
1438  &petsc_decide,ierr)
1439  chckerr('EPSSetDimensions failed')
1440  call epssettolerances(solver,tol_slepc(rich_lvl),max_it_slepc,ierr)
1441  chckerr('EPSSetTolerances failed')
1442 
1443  ! user output
1444  call writo('set tolerance to '//trim(r2str(tol_slepc(rich_lvl))))
1445  call writo('set maximum nr. of iterations to '//&
1446  &trim(i2str(max_it_slepc)))
1447  end function setup_solver
1448 
1452  integer function setup_guess(sol,A,solver) result(ierr)
1453  character(*), parameter :: rout_name = 'setup_guess'
1454 
1455  ! input / output
1456  type(sol_type), intent(in) :: sol
1457  mat, intent(in) :: a
1458  eps, intent(inout) :: solver
1459 
1460  ! local variables
1461  vec, allocatable :: guess_vec(:) ! guess to solution EV parallel vector
1462  petscint :: kd ! counter
1463  petscint :: n_ev_prev ! nr. of previous solutions
1464  petscscalar, pointer :: guess_vec_ptr(:) ! pointer to local guess_vec
1465  character(len=max_str_ln) :: err_msg ! error message
1466 
1467  ! initialize ierr
1468  ierr = 0
1469 
1470  call lvl_ud(1)
1471 
1472  ! set guess for EV if sol vec is allocated
1473  if (allocated(sol%val)) then
1474  ! set n_EV_prev
1475  n_ev_prev = size(sol%val)
1476 
1477  if (n_ev_prev.gt.0) then
1478  ! user output
1479  call writo('Set '//trim(i2str(n_ev_prev))//&
1480  &' vector(s) as guess')
1481 
1482  ! allocate guess vectors
1483  allocate(guess_vec(n_ev_prev))
1484 
1485  ! create vecctor guess_vec and set values
1486  do kd = 1,n_ev_prev
1487  ! create the vectors
1488  call matcreatevecs(a,guess_vec(kd),petsc_null_vec,ierr) ! get compatible parallel vectors to matrix A
1489  chckerr('Failed to create vector')
1490 
1491  ! get pointer
1492  call vecgetarrayf90(guess_vec(kd),guess_vec_ptr,ierr)
1493  chckerr('Failed to get pointer')
1494 
1495  ! copy the values
1496  guess_vec_ptr(1:size(sol%vec(:,:,kd))) = &
1497  &reshape([sol%vec(:,:,kd)],&
1498  &[size(sol%vec(:,:,kd))])
1499  guess_vec_ptr(size(sol%vec(:,:,kd))+1:size(guess_vec_ptr)) &
1500  &= guess_vec_ptr(size(sol%vec(:,:,kd))) ! some BC's have a grid extension
1501 
1502  ! return pointer
1503  call vecrestorearrayf90(guess_vec(kd),guess_vec_ptr,ierr)
1504  chckerr('Failed to restore pointer')
1505 
1506  !! visualize guess
1507  !call VecView(guess_vec(kd),PETSC_VIEWER_STDOUT_WORLD,ierr)
1508  !CHCKERR('Cannot view vector')
1509  end do
1510 
1511  ! set guess
1512  call epssetinitialspace(solver,n_ev_prev,guess_vec,ierr)
1513  chckerr('Failed to set guess')
1514 
1515  ! destroy guess vector
1516  do kd = 1,size(guess_vec)
1517  call vecdestroy(guess_vec(kd),ierr) ! destroy guess vector
1518  err_msg = 'Failed to destroy guess vector '//trim(i2str(kd))
1519  chckerr(err_msg)
1520  end do
1521  deallocate(guess_vec)
1522  else
1523  call writo('No vectors saved to use as guess')
1524  end if
1525  end if
1526 
1527  call lvl_ud(-1)
1528  end function setup_guess
1529 
1533  integer function get_solution(solver) result(ierr)
1534  character(*), parameter :: rout_name = 'get_solution'
1535 
1536  ! input / output
1537  eps, intent(inout) :: solver
1538 
1539  ! local variables
1540  character(len=max_str_ln) :: err_msg ! error message
1541 
1542  ! initialize ierr
1543  ierr = 0
1544 
1545  call lvl_ud(1)
1546 
1547  ! set run-time options
1548  call epssetfromoptions(solver,ierr)
1549  chckerr('EPSetFromOptions failed')
1550 
1551  ! solve EV problem
1552  call epssolve(solver,ierr)
1553  err_msg = 'EPS couldn''t find a solution. Maybe you should increase &
1554  &the number of parallel points.'
1555  chckerr(err_msg)
1556 
1557  call lvl_ud(-1)
1558  end function get_solution
1559 
1563  integer function summarize_solution(solver,max_n_EV) result (ierr)
1564  use num_vars, only: n_sol_requested
1565 
1566  character(*), parameter :: rout_name = 'summarize_solution'
1567 
1568  ! input / output
1569  eps, intent(inout) :: solver
1570  petscint, intent(inout) :: max_n_ev
1571 
1572  ! local variables
1573  petscint :: n_it ! nr. of iterations
1574  petscint :: n_conv ! nr. of converged solutions
1575  petscint :: n_ev, ncv, mpd ! nr. of requested EV, max. dim of subspace and for projected problem
1576  epstype :: eps_type ! string with name of type of EPS solver
1577  petscreal :: tol ! tolerance of EPS solver
1578  petscint :: max_it ! maximum number of iterations
1579 
1580  ! initialize ierr
1581  ierr = 0
1582 
1583  call lvl_ud(1)
1584 
1585  ! user output
1586  call epsgettype(solver,eps_type,ierr) ! EPS solver type
1587  chckerr('EPSGetType failed')
1588  call epsgettolerances(solver,tol,max_it,ierr) ! tolerance and max. nr. of iterations
1589  chckerr('EPSGetTolerances failed')
1590  call writo(trim(eps_type)//' solver with tolerance '//&
1591  &trim(r2strt(tol))//' and maximum '//trim(i2str(max_it))//&
1592  &' iterations')
1593  call epsgetconverged(solver,n_conv,ierr) ! nr. of converged solutions
1594  chckerr('EPSGetConverged failed')
1595  call epsgetiterationnumber(solver,n_it,ierr) ! nr. of iterations
1596  chckerr('EPSGetIterationNumber failed')
1597  call epsgetdimensions(solver,n_ev,ncv,mpd,ierr) ! nr. of requested EV, ...
1598  chckerr('EPSGetDimensions failed')
1599  call writo('number of iterations: '//trim(i2str(n_it)))
1600  call writo('number of converged solutions: '//trim(i2str(n_conv)))
1601  call writo('number of requested solutions: '//trim(i2str(n_ev)))
1602  call writo('maximum dimension of the subspace to be used by solver: '//&
1603  &trim(i2str(ncv)))
1604  call writo('maximum dimension allowed for projected problem : '//&
1605  &trim(i2str(mpd)))
1606 
1607  ! set maximum nr of solutions to be saved
1608  if (n_sol_requested.gt.n_conv) then
1609  if (n_conv.eq.0) then
1610  call writo('no solutions found',warning=.true.)
1611  else
1612  call writo('max. nr. of solutions found only '//&
1613  &trim(i2str(n_conv)),warning=.true.)
1614  end if
1615  max_n_ev = n_conv
1616  else
1617  max_n_ev = n_sol_requested
1618  end if
1619 
1620  ! check whether solutions found
1621  if (max_n_ev.le.0) then
1622  ierr = 1
1623  chckerr('No solutions found')
1624  end if
1625 
1626  call lvl_ud(-1)
1627  end function summarize_solution
1628 
1632  integer function store_results(mds,grid_sol,sol,solver,max_n_EV,A,B) &
1633  &result(ierr)
1634 
1635  use eq_vars, only: t_0
1636  use x_vars, only: n_mod_x
1637  use num_vars, only: use_normalization, ev_bc, prog_name, output_name, &
1638  &n_procs, rank, tol_slepc, eq_style, output_ev_i, retain_all_sol
1639  use mpi_utilities, only: get_ser_var
1640  use rich_vars, only: rich_lvl
1641  use grid_utilities, only: trim_grid
1642 
1643  character(*), parameter :: rout_name = 'store_results'
1644 
1645  ! input / output
1646  type(modes_type), intent(in) :: mds
1647  type(grid_type), intent(in) :: grid_sol
1648  type(sol_type), intent(inout) :: sol
1649  eps, intent(inout) :: solver
1650  petscint, intent(inout) :: max_n_ev
1651  mat, intent(inout) :: a
1652  mat, intent(inout) :: b
1653 
1654  ! local variables
1655  type(grid_type) :: grid_sol_trim ! trimmed solution grid
1656  mat :: err_mat ! A - omega^2 B
1657  vec :: sol_vec ! solution EV parallel vector
1658  vec :: err_vec ! AX - omega^2 BX
1659  vec :: e_vec(2) ! AX and BX
1660  petscint :: id, jd ! counters
1661  petscint :: id_tot ! counters
1662  petscint :: one = 1 ! one
1663  petscint, allocatable :: sol_vec_max_loc(:) ! location of sol_vec_max of this process
1664  petscreal :: error ! error of EPS solver
1665  petscreal :: error_est ! error estimate of EPS solver
1666  petscreal :: err_norm ! norm of error
1667  petscreal :: tol_complex ! tolerance for an EV to be considered complex
1668  petscreal :: tol_ev_bc = 1.e-3_dp ! tolerance for an EV to be considered due to the BC
1669  petscreal, parameter :: infinity = 1.e40_dp ! beyond this value, modes are not saved
1670  petscscalar :: sol_val_loc ! local X val
1671  petscscalar :: err_val ! X*(A-omega^2B)X
1672  petscscalar :: e_val(2) ! X*AX and X*BX
1673  petscscalar, allocatable :: sol_vec_loc(:,:) ! local solution vector, possibly in extended grid
1674  petscscalar, allocatable :: sol_vec_max(:) ! max of sol_vec of all processes, including phase at max
1675  character(len=max_str_ln) :: full_output_name ! full name
1676  character(len=max_str_ln) :: format_val ! format
1677  character(len=max_str_ln) :: format_head ! header
1678  character(len=max_str_ln) :: ev_err_str ! String with information about EV error
1679  character(len=max_str_ln) :: err_msg ! error message
1680  character(len=2*max_str_ln) :: temp_output_str ! temporary output string
1681  integer :: n_digits ! nr. of digits for the integer number
1682  integer :: n_err(3) ! how many errors there were
1683 
1684  ! initialize ierr
1685  ierr = 0
1686 
1687  call lvl_ud(1)
1688 
1689  ! trim grid
1690  ierr = trim_grid(grid_sol,grid_sol_trim)
1691  chckerr('')
1692 
1693  ! create solution variables
1694  if (allocated(sol%val)) call sol%dealloc()
1695  call sol%init(mds,grid_sol_trim,max_n_ev)
1696 
1697  ! create solution vector
1698  call veccreatempiwitharray(petsc_comm_world,one,loc_n_r*n_mod_x,&
1699  &n_r*n_mod_x,petsc_null_scalar,sol_vec,ierr)
1700  chckerr('Failed to create MPI vector with arrays')
1701 
1702  ! set up EV error string and format string:
1703  ! 1: index of EV
1704  ! 2: Real part
1705  ! 3: Imaginary part
1706  ! 4: relative error ||Ax - EV Bx||_2/||EV x||_2 [1]
1707  n_err = 0
1708  ev_err_str = ''
1709  if (rank.eq.0) then
1710  n_digits = min(ceiling(log10(1._dp*max_n_ev)),16)
1711  format_val = '(A2,I'//trim(i2str(n_digits))//'," ",ES23.16," ",&
1712  &ES23.16," ",ES23.16)'
1713  format_head = '(A'//trim(i2str(n_digits+3))//'," ",A23," ",A23," ",&
1714  &A23)'
1715  end if
1716 
1717  ! open output file for the log
1718  if (rank.eq.0) then
1719  full_output_name = prog_name//'_'//trim(output_name)//'_EV_R_'//&
1720  &trim(i2str(rich_lvl))//'.txt'
1721  open(output_ev_i,file=full_output_name,status='replace',&
1722  &iostat=ierr)
1723  chckerr('Cannot open EV output file')
1724  end if
1725 
1726  ! output to file
1727  if (rank.eq.0) then
1728  if (use_normalization) then
1729  write(unit=output_ev_i,fmt='(A)',iostat=ierr) &
1730  &'# Eigenvalues normalized to the squared Alfven &
1731  &frequency omega_A^2 = '
1732  chckerr('Failed to write')
1733  select case (eq_style)
1734  case (1) ! VMEC
1735  write(unit=output_ev_i,fmt='(A)',iostat=ierr) &
1736  &'# ('//trim(r2str(1._dp/t_0))//' Hz)^2 = '//&
1737  &trim(r2str(1._dp/t_0**2))//' Hz^2'
1738  chckerr('Failed to write')
1739  case (2) ! HELENA
1740  write(unit=output_ev_i,fmt='(A)',iostat=ierr) &
1741  &'# (MISHKA normalization)'
1742  chckerr('Failed to write')
1743  end select
1744  else
1745  write(unit=output_ev_i,fmt='(A)',iostat=ierr) '# Eigenvalues'
1746  chckerr('Failed to write')
1747  end if
1748  write(temp_output_str,format_head) &
1749  &'# I ', &
1750  &'real part ', 'imaginary part ', &
1751  &'relative precision '
1752  write(unit=output_ev_i,fmt='(A)',iostat=ierr) trim(temp_output_str)
1753  chckerr('Failed to write')
1754  end if
1755 
1756  ! initialize helper variables
1757  allocate(sol_vec_max(n_procs))
1758  allocate(sol_vec_max_loc(2))
1759  allocate(sol_vec_loc(n_mod_x,loc_n_r))
1760 
1761  ! start id
1762  id = 1
1763  id_tot = 1
1764 
1765  ! store them
1766  do while (id.le.max_n_ev)
1767  call vecplacearray(sol_vec,sol_vec_loc,ierr) ! place array local sol_vec in solution vec
1768  chckerr('Failed to place array')
1769  call epsgeteigenpair(solver,id_tot-1,sol%val(id),petsc_null_scalar,&
1770  &sol_vec,petsc_null_vec,ierr) ! get solution EV vector and value (starts at index 0)
1771  chckerr('EPSGetEigenpair failed')
1772  sol%vec(:,:,id) = sol_vec_loc(:,1:size(sol%vec,2))
1773  call epscomputeerror(solver,id_tot-1,eps_error_relative,error,ierr) ! get error (starts at index 0) (petsc 3.6.1)
1774  !call EPSComputeRelativeError(solver,id_tot-1,error,ierr) ! get error (starts at index 0) (petsc 3.5.3)
1775  chckerr('EPSComputeError failed')
1776 
1777  ! set up local solution val
1778  sol_val_loc = sol%val(id)
1779 
1780  ! tests
1781  tol_complex = 10*tol_slepc(rich_lvl)
1782  if (abs(ip(sol%val(id))/rp(sol%val(id))).gt.tol_complex) then ! test for unphysical complex solution
1783  ev_err_str = '# WARNING: Unphysical complex Eigenvalue!'
1784  n_err(1) = n_err(1)+1
1785  call remove_ev(id,max_n_ev)
1786  else if (abs(rp(sol%val(id)-ev_bc)/ev_bc).lt.tol_ev_bc) then ! test for artificial EV due to BC
1787  ev_err_str = '# WARNING: Eigenvalue probably due to BC''s &
1788  &(artifically set to '//trim(r2strt(ev_bc))//')'
1789  n_err(2) = n_err(2)+1
1790  call remove_ev(id,max_n_ev)
1791  else if (abs(sol%val(id)).gt.infinity) then ! test for infinity
1792  ev_err_str = '# WARNING: The next Eigenvalues are larger &
1793  &than '//trim(r2strt(infinity))//' and are discarded'
1794  n_err(3) = 1
1795  call remove_ev(id,max_n_ev,remove_next=.true.)
1796  exit
1797  end if
1798 
1799  if (rank.eq.0) then
1800  ! output EV to file
1801  if (ev_err_str.ne.'') then
1802  write(temp_output_str,format_val) '#', id,rp(sol_val_loc),&
1803  &ip(sol_val_loc),error ! if error, comment the next line
1804  else
1805  write(temp_output_str,format_val) '', id,rp(sol_val_loc),&
1806  &ip(sol_val_loc),error
1807  end if
1808  write(unit=output_ev_i,fmt='(A)',iostat=ierr) &
1809  &trim(temp_output_str)
1810  chckerr('Failed to write')
1811 
1812  ! if error, print explanation
1813  if (ev_err_str.ne.'') then
1814  write(unit=output_ev_i,fmt='(A)',iostat=ierr) &
1815  &trim(ev_err_str)
1816  chckerr('Failed to write')
1817  end if
1818  end if
1819 
1820  ! normalize Eigenvectors to make output more easily comparable
1821  if (ev_err_str.eq.'') then
1822  ! find local maximum
1823  sol_vec_max_loc = maxloc(abs(sol%vec(:,:,id)))
1824  ierr = get_ser_var(&
1825  &[sol%vec(sol_vec_max_loc(1),sol_vec_max_loc(2),id)],&
1826  &sol_vec_max,scatter=.true.)
1827  chckerr('')
1828  ! find global maximum
1829  sol_vec_max_loc(1) = maxloc(abs(sol_vec_max),1)
1830  sol%vec(:,:,id) = sol%vec(:,:,id) / &
1831  &sol_vec_max(sol_vec_max_loc(1))
1832 
1833  call epsgeterrorestimate(solver,id_tot-1,error_est,ierr) ! get error estimate
1834  chckerr('EPSGetErrorEstimate failed')
1835  ! user message
1836  call writo('Checking whether A x - omega^2 B x = 0 for EV '//&
1837  &trim(i2str(id))//': '//trim(c2strt(sol%val(id))))
1838  call lvl_ud(1)
1839 
1840  ! set up error matrix A - omega^2 B
1841  call matduplicate(a,mat_share_nonzero_pattern,err_mat,ierr)
1842  err_msg = 'failed to duplicate mat into err_mat'
1843  chckerr(err_msg)
1844  call matcopy(a,err_mat,mat_share_nonzero_pattern,ierr) ! err_mat has same structure as A
1845  chckerr('Failed to copy mat into mat_loc')
1846  call mataxpy(err_mat,-sol%val(id),b,different_nonzero_pattern,&
1847  &ierr) ! for some reason, SAME_NONZERO_PATTERN does not work
1848  chckerr('Failed to perform AXPY')
1849 
1850  ! set up error vector
1851  call vecduplicate(sol_vec,err_vec,ierr)
1852  chckerr('Failed to duplicate vector')
1853 
1854  ! calculate Ax-lambdaBx
1855  call matmult(err_mat,sol_vec,err_vec,ierr)
1856  chckerr('Failed to multiply')
1857 
1858  ! multiply with X* and give output
1859  call vecdot(err_vec,sol_vec,err_val,ierr)
1860  chckerr('Failed to do dot product')
1861  call writo('X*(A-omega^2B)X = '//trim(c2strt(err_val)))
1862 
1863  ! calculate absolute norm
1864  call vecnorm(err_vec,norm_2,err_norm,ierr)
1865  chckerr('Failed to calculate norm')
1866 
1867  ! get relative norm
1868  err_norm = err_norm/abs(sol%val(id))
1869 
1870  ! visualize solution and error
1871  call writo('error: '//trim(r2str(err_norm))//', given: '//&
1872  &trim(r2str(error))//', estimate: '//trim(r2str(error_est)))
1873  !!call VecView(err_vec,PETSC_VIEWER_STDOUT_WORLD,ierr)
1874  !!CHCKERR('Cannot view vector')
1875 
1876  ! set up Energy vector
1877  call vecduplicate(sol_vec,e_vec(1),ierr)
1878  chckerr('Failed to duplicate vector')
1879  call vecduplicate(sol_vec,e_vec(2),ierr)
1880  chckerr('Failed to duplicate vector')
1881 
1882  ! calculate Ax and Bx
1883  call matmult(a,sol_vec,e_vec(1),ierr)
1884  chckerr('Failed to multiply')
1885  call matmult(b,sol_vec,e_vec(2),ierr)
1886  chckerr('Failed to multiply')
1887 
1888  ! multiply with X* and give output
1889  call vecdot(e_vec(1),sol_vec,e_val(1),ierr)
1890  chckerr('Failed to do dot product')
1891  call writo('step_size = '//trim(r2str(step_size)))
1892  call writo('E_pot = X*AX*step_size = '//&
1893  &trim(c2str(e_val(1)*step_size)))
1894  call vecdot(e_vec(2),sol_vec,e_val(2),ierr)
1895  chckerr('Failed to do dot product')
1896  call writo('E_kin = X*BX*step_size = '//&
1897  &trim(c2str(e_val(2)*step_size)))
1898  call writo('X*AX/X*BX = '//trim(c2str(e_val(1)/e_val(2))))
1899  call writo('omega^2 = '//trim(c2str(sol%val(id))))
1900 
1901  ! if stable EV, give warning
1902  if (rp(sol%val(id)).gt.0._dp) then
1903  call writo('The eigenvalue is positive, so the mode could &
1904  &be stable.',warning=.true.)
1905  call lvl_ud(1)
1906  call writo('But this could also be due to value of "ncv" &
1907  &that is too low.')
1908  call writo('Have a look at the POST output and/or try &
1909  &increasing "ncv".')
1910  call lvl_ud(-1)
1911  end if
1912 
1913  ! increment counter
1914  id = id + 1
1915 
1916  call lvl_ud(-1)
1917 
1918  ! clean up
1919  call vecdestroy(err_vec,ierr) ! destroy error vector
1920  chckerr('Failed to destroy err_vec')
1921  call matdestroy(err_mat,ierr) ! destroy error matrix
1922  chckerr('Failed to destroy err_mat')
1923  do jd = 1,size(e_vec)
1924  call vecdestroy(e_vec(jd),ierr) ! destroy energy vector
1925  err_msg = 'Failed to destroy E_vec '//trim(i2str(jd))
1926  chckerr(err_msg)
1927  end do
1928  else
1929  ! reinitialize error string if error
1930  ev_err_str = ''
1931  end if
1932 
1933  ! reset vector
1934  call vecresetarray(sol_vec,ierr)
1935  chckerr('Failed to reset array')
1936 
1937  ! increment total id
1938  id_tot = id_tot+1
1939  end do
1940 
1941  ! deallocate helper variables
1942  deallocate(sol_vec_max_loc)
1943  deallocate(sol_vec_max)
1944 
1945  ! close output file if master
1946  if (rank.eq.0) close(output_ev_i)
1947 
1948  ! user output
1949  if (rank.eq.0) call writo(trim(i2str(max_n_ev))//&
1950  &' Eigenvalues were written in the file "'//&
1951  &trim(full_output_name)//'"')
1952  if (sum(n_err).ne.0) then
1953  call writo('there were some Eigenvalues that were removed:')
1954  call lvl_ud(1)
1955 
1956  if (n_err(1).gt.0) call writo(trim(i2str(n_err(1)))//&
1957  &' unphysical complex Eigenvalues')
1958  if (n_err(2).gt.0) call writo(trim(i2str(n_err(2)))//&
1959  &' Eigenvalues due to the boundary conditions')
1960  if (n_err(3).gt.0) call writo('The final Eigenvalues became &
1961  &infinite')
1962 
1963  call lvl_ud(-1)
1964  if (retain_all_sol) then
1965  call writo('(But this was overriden using "retain_all_sol")')
1966  else
1967  call writo('(Override this behavior using "retain_all_sol")')
1968  end if
1969  end if
1970 
1971  if (max_n_ev.gt.0) then
1972  call writo('basic statistics:')
1973  call lvl_ud(1)
1974 
1975  call writo('min: '//trim(c2strt(sol%val(1))))
1976  call writo('max: '//trim(c2strt(sol%val(max_n_ev))))
1977 
1978  call lvl_ud(-1)
1979  end if
1980 
1981  !!call EPSErrorView(solver,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr)
1982  !!CHCKERR('Failed to print solution')
1983 
1984  call vecdestroy(sol_vec,ierr) ! destroy solution vector
1985  chckerr('Failed to destroy sol_vec')
1986 
1987  call lvl_ud(-1)
1988 
1989  ! clean up
1990  call grid_sol_trim%dealloc()
1991  contains
1992  ! Removes id'th Eigenvalue and -vector. Optionally, all the following
1993  ! can be removed as well.
1995  subroutine remove_ev(id,max_id,remove_next)
1996  ! input / output
1997  petscint, intent(inout) :: id ! id of faulty values
1998  petscint, intent(inout) :: max_id ! maximum id
1999  petscbool, intent(in), optional :: remove_next ! whether all next values have to be removed as well
2000 
2001  ! local variables
2002  petscscalar, allocatable :: sol_val_loc(:) ! local copy of sol_val
2003  petscscalar, allocatable :: sol_vec_loc(:,:,:) ! local copy of sol_vec
2004  petscbool :: remove_next_loc = .false. ! local copy of remove_next
2005 
2006  ! only remove if faulty solutions are not optionally retained
2007  if (retain_all_sol) then
2008  id = id+1 ! increment the solution
2009  else
2010  ! set up local remove_next
2011  if (present(remove_next)) remove_next_loc = remove_next
2012 
2013  ! save old arrays
2014  allocate(sol_val_loc(max_id))
2015  allocate(sol_vec_loc(n_mod_x,grid_sol_trim%loc_n_r,max_id))
2016  sol_val_loc = sol%val
2017  sol_vec_loc = sol%vec
2018 
2019  ! reallocate arrays
2020  deallocate(sol%val,sol%vec)
2021  if (remove_next_loc) then
2022  allocate(sol%val(id-1))
2023  allocate(sol%vec(n_mod_x,grid_sol_trim%loc_n_r,id-1))
2024  else
2025  allocate(sol%val(max_id-1))
2026  allocate(sol%vec(n_mod_x,grid_sol_trim%loc_n_r,max_id-1))
2027  end if
2028 
2029  ! copy other values
2030  sol%val(1:id-1) = sol_val_loc(1:id-1)
2031  sol%vec(:,:,1:id-1) = sol_vec_loc(:,:,1:id-1)
2032  if (.not.remove_next_loc) then
2033  sol%val(id:max_id-1) = sol_val_loc(id+1:max_id)
2034  sol%vec(:,:,id:max_id-1) = sol_vec_loc(:,:,id+1:max_id)
2035  end if
2036 
2037  ! adapt max_id
2038  if (remove_next_loc) then
2039  max_id = id-1
2040  else
2041  max_id = max_id-1
2042  end if
2043 
2044  ! clean up
2045  deallocate(sol_val_loc,sol_vec_loc)
2046  end if
2047  end subroutine remove_ev
2048  end function store_results
2049 
2053  integer function stop_slepc(A,B,solver) result(ierr)
2054  character(*), parameter :: rout_name = 'stop_SLEPC'
2055 
2056  ! input / output
2057  mat, intent(in) :: a, b
2058  eps, intent(in) :: solver
2059 
2060  ! initialize ierr
2061  ierr = 0
2062 
2063  call lvl_ud(1)
2064 
2065  ! destroy objects
2066  call epsdestroy(solver,ierr)
2067  chckerr('Failed to destroy EPS solver')
2068  call matdestroy(a,ierr)
2069  chckerr('Failed to destroy matrix A')
2070  call matdestroy(b,ierr)
2071  chckerr('Failed to destroy matrix B')
2072 
2073  ! stop SLEPC
2074  call slepcfinalize(ierr)
2075  chckerr('Failed to Finalize SLEPC')
2076 
2077  call lvl_ud(-1)
2078  end function stop_slepc
2079 end module
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
input_utilities::get_real
real(dp) function, public get_real(lim_lo, lim_hi, ind)
Queries for user input for a real value, where allowable range can be provided as well.
Definition: input_utilities.f90:77
num_vars::output_ev_i
integer, parameter, public output_ev_i
file number of output of EV
Definition: num_vars.f90:188
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
num_utilities::con
Either takes the complex conjugate of a square matrix element A, defined on a 3-D grid,...
Definition: num_utilities.f90:221
num_vars::use_normalization
logical, public use_normalization
whether to use normalization or not
Definition: num_vars.f90:115
num_vars::bc_style
integer, dimension(2), public bc_style
style for BC left and right
Definition: num_vars.f90:94
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
input_utilities::get_log
logical function, public get_log(yes, ind)
Queries for a logical value yes or no, where the default answer is also to be provided.
Definition: input_utilities.f90:22
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
files_ops::opt_args
character(len=max_str_ln), dimension(:), allocatable, public opt_args
optional arguments that can be passed using –[name]
Definition: files_ops.f90:19
rich_vars
Variables concerning Richardson extrapolation.
Definition: rich_vars.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
num_vars::retain_all_sol
logical, public retain_all_sol
retain also faulty solutions
Definition: num_vars.f90:152
num_vars::n_sol_requested
integer, public n_sol_requested
how many solutions requested
Definition: num_vars.f90:170
slepc_utilities
Numerical utilities related to SLEPC (and PETSC) operations.
Definition: SLEPC_utilities.f90:8
eq_vars::t_0
real(dp), public t_0
derived normalization constant for nondimensionalization
Definition: eq_vars.f90:47
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
slepc_ops::debug_setup_mats
logical, public debug_setup_mats
plot debug information for setup_mats
Definition: SLEPC_ops.f90:45
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
num_vars::norm_disc_prec_sol
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition: num_vars.f90:122
num_vars::ev_guess
real(dp), public ev_guess
first guess for eigenvalue
Definition: num_vars.f90:117
slepc_ops::start_slepc
integer function, public start_slepc(n_r_sol)
This subroutine starts PETSC and SLEPC with the correct number of processors.
Definition: SLEPC_ops.f90:333
num_vars::data_dir
character(len=4), public data_dir
directory where to save data for plots
Definition: num_vars.f90:155
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition: str_utilities.f90:54
num_vars::prog_name
character(len=4), public prog_name
name of program, used for info
Definition: num_vars.f90:54
num_vars::ltest
logical, public ltest
whether or not to call the testing routines
Definition: num_vars.f90:112
slepc_ops::store_results
integer function, public store_results(mds, grid_sol, sol, solver, max_n_EV, A, B)
Stores the results.
Definition: SLEPC_ops.f90:1634
vac_vars
Variables pertaining to the vacuum quantities.
Definition: vac_vars.f90:4
sol_vars::sol_type
solution type
Definition: sol_vars.f90:30
num_vars::max_it_slepc
integer, public max_it_slepc
maximum nr. of iterations for SLEPC
Definition: num_vars.f90:101
num_vars::norm_disc_style_sol
integer, public norm_disc_style_sol
style for normal discretization for solution (1: central fin. diff., 2: left fin. diff....
Definition: num_vars.f90:123
slepc_ops::setup_mats
integer function, public setup_mats(mds, grid_sol, X, A, B)
Sets up the matrices and in the EV problem .
Definition: SLEPC_ops.f90:409
num_vars::ev_bc
real(dp), public ev_bc
value of artificial Eigenvalue for boundary condition
Definition: num_vars.f90:116
num_utilities::calc_coeff_fin_diff
integer function, public calc_coeff_fin_diff(deriv, nr, ind, coeff)
Calculate the coefficients for finite differences.
Definition: num_utilities.f90:2449
x_vars::x_2_type
tensorial perturbation type
Definition: X_vars.f90:81
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition: num_vars.f90:89
num_vars::tol_slepc
real(dp), dimension(:), allocatable, public tol_slepc
tolerance for SLEPC for different Richardson levels
Definition: num_vars.f90:118
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
rich_vars::use_guess
logical, public use_guess
whether a guess is formed from previous level of Richardson
Definition: rich_vars.f90:23
mpi_utilities::wait_mpi
integer function, public wait_mpi()
Wait for all processes, wrapper to MPI barrier.
Definition: MPI_utilities.f90:744
slepc_ops
Operations that use SLEPC (and PETSC) routines.
Definition: SLEPC_ops.f90:10
slepc_ops::test_diff
logical, public test_diff
test introduction of numerical diff
Definition: SLEPC_ops.f90:47
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
slepc_ops::get_solution
integer function, public get_solution(solver)
Get the solution vectors.
Definition: SLEPC_ops.f90:1534
slepc_ops::solve_ev_system_slepc
integer function, public solve_ev_system_slepc(mds, grid_sol, X, vac, sol)
Solve the eigenvalue system using SLEPC.
Definition: SLEPC_ops.f90:62
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition: grid_utilities.f90:4
str_utilities::c2strt
elemental character(len=max_str_ln) function, public c2strt(k)
Convert a complex (double) to string.
Definition: str_utilities.f90:88
files_ops
Operations related to files !
Definition: files_ops.f90:4
vac_vars::vac_type
vacuum type
Definition: vac_vars.f90:46
slepc_utilities::insert_block_mat
integer function, public insert_block_mat(mds, block, mat, r_id, ind, n_r, transp, overwrite, ind_insert)
Insert a block pertaining to 1 normal point into a matrix.
Definition: SLEPC_utilities.f90:83
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
slepc_ops::setup_solver
integer function, public setup_solver(X, A, B, solver)
Sets up EV solver.
Definition: SLEPC_ops.f90:1326
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
slepc_ops::debug_set_bc
logical, public debug_set_bc
plot debug information for set_BC
Definition: SLEPC_ops.f90:46
slepc_ops::setup_guess
integer function, public setup_guess(sol, A, solver)
Sets up guess in solver.
Definition: SLEPC_ops.f90:1453
num_vars::output_name
character(len=3), parameter, public output_name
name of output file
Definition: num_vars.f90:55
rich_vars::rich_lvl
integer, public rich_lvl
current level of Richardson extrapolation
Definition: rich_vars.f90:19
x_vars::x_1_type
vectorial perturbation type
Definition: X_vars.f90:51
sol_vars
Variables pertaining to the solution quantities.
Definition: sol_vars.f90:4
slepc_ops::summarize_solution
integer function, public summarize_solution(solver, max_n_EV)
Summarizes solution.
Definition: SLEPC_ops.f90:1564
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
num_vars::solver_slepc_style
integer, public solver_slepc_style
style for solver (1: Krylov-Schur, 2: GD)
Definition: num_vars.f90:97
num_vars::sol_n_procs
integer, public sol_n_procs
nr. of MPI processes for solution with SLEPC
Definition: num_vars.f90:70
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
str_utilities::c2str
elemental character(len=max_str_ln) function, public c2str(k)
Convert a complex (double) to string.
Definition: str_utilities.f90:66
input_utilities
Numerical utilities related to input.
Definition: input_utilities.f90:4
num_vars::matrix_slepc_style
integer, public matrix_slepc_style
style for matrix storage (1: sparse, 2: shell)
Definition: num_vars.f90:96
slepc_ops::stop_slepc
integer function, public stop_slepc(A, B, solver)
Stop PETSC and SLEPC.
Definition: SLEPC_ops.f90:2054
slepc_ops::set_bc
integer function, public set_bc(mds, grid_sol, X, vac, A, B)
Sets the boundary conditions deep in the plasma (1) and at the edge (2).
Definition: SLEPC_ops.f90:933
grid_utilities::trim_grid
integer function, public trim_grid(grid_in, grid_out, norm_id)
Trim a grid, removing any overlap between the different regions.
Definition: grid_utilities.f90:1636
this documentation was created with Doxygen 1.8.17 by Toon Weyens