PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
SLEPC_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations that use SLEPC (and PETSC) routines.
3!!
4!! \note The routines here require a \b TRIMMED solution grid!
5!------------------------------------------------------------------------------!
6!> \see
7!! References:
8!! \cite Hernandez2005slepc
9!------------------------------------------------------------------------------!
10module 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. !< plot debug information for setup_mats \ldebug
46 logical :: debug_set_bc = .false. !< plot debug information for set_BC \ldebug
47 logical :: test_diff = .false. !< test introduction of numerical diff \ldebug
48 real(dp) :: diff_coeff !< diff coefficient \ldebug
49#endif
50
51contains
52 !> Solve the eigenvalue system using SLEPC.
53 !!
54 !! This subroutine sets up the matrices A ad B of the generalized EV system
55 !! described in \cite Weyens2017PB3D and solves them using the SLEPC suite.
56 !! The solutions closest to a target indicated by \c EV_guess are obtained.
57 !! \see read_input_opts()
58 !!
59 !! \return ierr
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 !< general modes variables
74 type(grid_type), intent(in) :: grid_sol !< solution grid
75 type(x_2_type), intent(in) :: x !< field-averaged perturbation variables (so only first index)
76 type(vac_type), intent(in) :: vac !< vacuum variables
77 type(sol_type), intent(inout) :: sol !< solution variables
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.
225 !> \private
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
328 !> This subroutine starts PETSC and SLEPC with the correct number of
329 !! processors
330 !!
331 !! \return ierr
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 !< nr. of normal points in solution grid
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
398 !> Sets up the matrices \f$\overline{\text{A}}\f$ and
399 !! \f$\overline{\text{B}}\f$ in the EV problem \f$ \overline{\text{A}}
400 !! \vec{X} = \lambda \overline{\text{B}} \vec{X}. \f$
401 !!
402 !! The matrices are set up in the solution grid, so some processes might not
403 !! have any part in the storage thereof (if \c sol_n_procs < \c n_procs),
404 !! but each process sets up the part of the grid that corresponds to their
405 !! own normal range in the solution grid.
406 !!
407 !! \return ierr
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 !< general modes variables
425 type(grid_type), intent(in) :: grid_sol !< solution grid
426 type(x_2_type), intent(in), target :: x !< field-averaged perturbation variables (so only first index)
427 mat, intent(inout) :: a !< matrix A
428 mat, intent(inout) :: b !< matrix 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.
569 !> \private
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
622 !> \private
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.
797 !> \private
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.
834 !> \private
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
899 !> Sets the boundary conditions deep in the plasma (1) and at the edge (2).
900 !!
901 !! Through \c norm_disc_style_sol it can be chosen whether the finite
902 !! differences are:
903 !! 1. central finite differences
904 !! 2. left finite differences
905 !!
906 !! Possibilities for \c BC_style:
907 !! -# Set to zero:
908 !! - An artificial Eigenvalue \c EV_BC is introduced by setting the
909 !! diagonal components of A to EV_BC and of B to 1, and the
910 !! off-diagonal elements to zero.
911 !! -# Minimization of surface energy through asymmetric fin. differences:
912 !! - For symmetric finite differences, the last \c ndps grid points are
913 !! treated asymmetrically in order not go over the edge and the vacuum
914 !! term is added to the edge element.
915 !! - For left differences, this is already standard, so the method
916 !! becomes identical to 2.
917 !! -# Minimization of surface energy through extension of grid:
918 !! - For symmetric finite differences, \c ndps extra grid points are
919 !! introduced after the edge and the vacuum term is added to the edge
920 !! element.
921 !! - For left finite differences, the vacuum term is just added to the
922 !! edge element, so this method becomes idential to 3.
923 !! -# Explicit introduction of the surface energy minimization:
924 !! - The equation due to the minimization of the vacuum term is
925 !! introduced explicitely as an asymmetric finite difference equation
926 !! in the last row.
927 !! -This is done using left finite differences.
928 !!
929 !! Makes use of n_r.
930 !!
931 !! \return ierr
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 !< general modes variables
944 type(grid_type), intent(in) :: grid_sol !< solution grid
945 type(x_2_type), intent(in) :: x !< field-averaged perturbation variables (so only first index)
946 type(vac_type), intent(in) :: vac !< vacuum variables
947 mat, intent(inout) :: a !< Matrix A from A X = lambda B X
948 mat, intent(inout) :: b !< Matrix B from A X = lambda B X
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.
1095 !> \private
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
1164 !> \private
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
1193 !> \private
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.
1225 !> \private
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
1322 !> Sets up EV solver.
1323 !!
1324 !! \return ierr
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 !< field-averaged perturbation variables (so only first index)
1335 mat, intent(in) :: a, b !< matrix A and B
1336 eps, intent(inout) :: solver !< EV 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
1449 !> Sets up guess in solver.
1450 !!
1451 !! \return ierr
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 !< solution variables
1457 mat, intent(in) :: a !< matrix A (or B)
1458 eps, intent(inout) :: solver !< EV 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
1530 !> Get the solution vectors.
1531 !!
1532 !! \return ierr
1533 integer function get_solution(solver) result(ierr)
1534 character(*), parameter :: rout_name = 'get_solution'
1535
1536 ! input / output
1537 eps, intent(inout) :: solver !< EV 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
1560 !> Summarizes solution.
1561 !!
1562 !! \return ierr
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 !< EV solver
1570 petscint, intent(inout) :: max_n_ev !< nr. of EV's saved, up to \c n_conv
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
1629 !> Stores the results.
1630 !!
1631 !! \return ierr
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 !< general modes variables
1647 type(grid_type), intent(in) :: grid_sol !< solution grid
1648 type(sol_type), intent(inout) :: sol !< solution variables
1649 eps, intent(inout) :: solver !< EV solver
1650 petscint, intent(inout) :: max_n_ev !< nr. of EV's saved, up to \c n_conv
1651 mat, intent(inout) :: a !< Matrix A from A X = lambda B X
1652 mat, intent(inout) :: b !< Matrix B from A X = lambda B X
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.
1994 !> \private
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
2050 !> Stop PETSC and SLEPC.
2051 !!
2052 !! \return ierr
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 !< matrices A and B in EV problem
2058 eps, intent(in) :: solver !< EV 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
2079end module
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition MPI_utilities.f90:55
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
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition eq_vars.f90:27
eq_vars::t_0
real(dp), public t_0
derived normalization constant for nondimensionalization
Definition eq_vars.f90:47
files_ops
Operations related to files !
Definition files_ops.f90:4
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
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition grid_utilities.f90:4
grid_utilities::trim_grid
integer function, public trim_grid(grid_in, grid_out, norm_id)
Trim a grid, removing any overlap between the different regions.
Definition grid_utilities.f90:1636
grid_vars
Variables pertaining to the different grids used.
Definition grid_vars.f90:4
grid_vars::n_r_sol
integer, public n_r_sol
nr. of normal points in solution grid
Definition grid_vars.f90:22
input_utilities
Numerical utilities related to input.
Definition input_utilities.f90:4
input_utilities::get_real
real(dp) function, public get_real(lim_lo, lim_hi, ind)
Queries for user input for a real value, where allowable range can be provided as well.
Definition input_utilities.f90:77
input_utilities::get_log
logical function, public get_log(yes, ind)
Queries for a logical value yes or no, where the default answer is also to be provided.
Definition input_utilities.f90:22
messages
Numerical utilities related to giving output.
Definition messages.f90:4
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition messages.f90:254
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition messages.f90:275
mpi_utilities
Numerical utilities related to MPI.
Definition MPI_utilities.f90:20
mpi_utilities::wait_mpi
integer function, public wait_mpi()
Wait for all processes, wrapper to MPI barrier.
Definition MPI_utilities.f90:744
num_utilities
Numerical utilities.
Definition num_utilities.f90:4
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
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
num_utilities::m
integer, dimension(:,:), allocatable, public m
1-D array indices of metric indices
Definition num_utilities.f90:24
num_vars
Numerical variables used by most other modules.
Definition num_vars.f90:4
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::dp
integer, parameter, public dp
double precision
Definition num_vars.f90:46
num_vars::ltest
logical, public ltest
whether or not to call the testing routines
Definition num_vars.f90:112
num_vars::output_name
character(len=3), parameter, public output_name
name of output file
Definition num_vars.f90:55
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::n_procs
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition num_vars.f90:85
num_vars::n_sol_requested
integer, public n_sol_requested
how many solutions requested
Definition num_vars.f90:170
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
num_vars::ev_guess
real(dp), public ev_guess
first guess for eigenvalue
Definition num_vars.f90:117
num_vars::ev_bc
real(dp), public ev_bc
value of artificial Eigenvalue for boundary condition
Definition num_vars.f90:116
num_vars::use_normalization
logical, public use_normalization
whether to use normalization or not
Definition num_vars.f90:115
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
num_vars::prog_name
character(len=4), public prog_name
name of program, used for info
Definition num_vars.f90:54
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition num_vars.f90:89
num_vars::max_it_slepc
integer, public max_it_slepc
maximum nr. of iterations for SLEPC
Definition num_vars.f90:101
num_vars::output_ev_i
integer, parameter, public output_ev_i
file number of output of EV
Definition num_vars.f90:188
num_vars::bc_style
integer, dimension(2), public bc_style
style for BC left and right
Definition num_vars.f90:94
num_vars::rank
integer, public rank
MPI rank.
Definition num_vars.f90:68
num_vars::tol_slepc
real(dp), dimension(:), allocatable, public tol_slepc
tolerance for SLEPC for different Richardson levels
Definition num_vars.f90:118
num_vars::retain_all_sol
logical, public retain_all_sol
retain also faulty solutions
Definition num_vars.f90:152
num_vars::matrix_slepc_style
integer, public matrix_slepc_style
style for matrix storage (1: sparse, 2: shell)
Definition num_vars.f90:96
num_vars::data_dir
character(len=4), public data_dir
directory where to save data for plots
Definition num_vars.f90:155
num_vars::solver_slepc_style
integer, public solver_slepc_style
style for solver (1: Krylov-Schur, 2: GD)
Definition num_vars.f90:97
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
rich_vars
Variables concerning Richardson extrapolation.
Definition rich_vars.f90:4
rich_vars::rich_lvl
integer, public rich_lvl
current level of Richardson extrapolation
Definition rich_vars.f90:19
rich_vars::use_guess
logical, public use_guess
whether a guess is formed from previous level of Richardson
Definition rich_vars.f90:23
slepc_ops
Operations that use SLEPC (and PETSC) routines.
Definition SLEPC_ops.f90:10
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
slepc_ops::setup_guess
integer function, public setup_guess(sol, a, solver)
Sets up guess in solver.
Definition SLEPC_ops.f90:1453
slepc_ops::summarize_solution
integer function, public summarize_solution(solver, max_n_ev)
Summarizes solution.
Definition SLEPC_ops.f90:1564
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
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
slepc_ops::test_diff
logical, public test_diff
test introduction of numerical diff
Definition SLEPC_ops.f90:47
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
slepc_ops::get_solution
integer function, public get_solution(solver)
Get the solution vectors.
Definition SLEPC_ops.f90:1534
slepc_ops::debug_set_bc
logical, public debug_set_bc
plot debug information for set_BC
Definition SLEPC_ops.f90:46
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
slepc_ops::stop_slepc
integer function, public stop_slepc(a, b, solver)
Stop PETSC and SLEPC.
Definition SLEPC_ops.f90:2054
slepc_ops::setup_solver
integer function, public setup_solver(x, a, b, solver)
Sets up EV solver.
Definition SLEPC_ops.f90:1326
slepc_ops::debug_setup_mats
logical, public debug_setup_mats
plot debug information for setup_mats
Definition SLEPC_ops.f90:45
slepc_utilities
Numerical utilities related to SLEPC (and PETSC) operations.
Definition SLEPC_utilities.f90:8
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
sol_vars
Variables pertaining to the solution quantities.
Definition sol_vars.f90:4
str_utilities
Operations on strings.
Definition str_utilities.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition str_utilities.f90:18
str_utilities::c2str
elemental character(len=max_str_ln) function, public c2str(k)
Convert a complex (double) to string.
Definition str_utilities.f90:66
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition str_utilities.f90:42
str_utilities::c2strt
elemental character(len=max_str_ln) function, public c2strt(k)
Convert a complex (double) to string.
Definition str_utilities.f90:88
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition str_utilities.f90:54
vac_vars
Variables pertaining to the vacuum quantities.
Definition vac_vars.f90:4
x_vars
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
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
grid_vars::grid_type
Type for grids.
Definition grid_vars.f90:59
sol_vars::sol_type
solution type
Definition sol_vars.f90:30
vac_vars::vac_type
vacuum type
Definition vac_vars.f90:46
x_vars::modes_type
mode number type
Definition X_vars.f90:36
x_vars::x_1_type
vectorial perturbation type
Definition X_vars.f90:51
x_vars::x_2_type
tensorial perturbation type
Definition X_vars.f90:81
this documentation was created with Doxygen 1.16.1 by Toon Weyens