PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
X_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations considering perturbation quantities.
3!------------------------------------------------------------------------------!
4module x_ops
5#include <PB3D_macros.h>
6#include <wrappers.h>
7 use str_utilities
8 use output_ops
9 use messages
10 use num_vars, only: dp, iu, max_str_ln, max_name_ln, pi
11 use grid_vars, onlY: grid_type
12 use eq_vars, only: eq_1_type, eq_2_type
13 use x_vars, only: x_1_type, x_2_type, modes_type
14
15 implicit none
16 private
17 public calc_x, calc_magn_ints, print_output_x, resonance_plot, &
18 &calc_res_surf, check_x_modes, init_modes, setup_modes, divide_x_jobs, &
19 &redistribute_output_x
20#if ldebug
21 public debug_check_x_modes_2, print_debug_x_1, print_debug_x_2
22#endif
23
24 ! global variables
25#if ldebug
26 logical :: debug_check_x_modes_2 = .false. !< plot debug information for check_x_modes_2() \ldebug
27 logical :: debug_setup_modes = .false. !< plot debug information for setup_modes() \ldebug
28#endif
29
30 ! interfaces
31
32 !> \public Calculates either vectorial or tensorial perturbation variables.
33 !!
34 !! Optionally, the secondary mode number can be specified (\c m if poloidal
35 !! flux is used as normal coordinate and c n if toroidal flux). By default,
36 !! they are taken from the global \c X_vars variables.
37 !!
38 !! \return ierr
39 interface calc_x
40 !> \public
41 module procedure calc_x_1
42 !> \public
43 module procedure calc_x_2
44 end interface
45
46 !> \public Redistribute the perturbation variables.
47 !!
48 !! \see redistribute_output_grid()
49 !!
50 !! \return ierr
51 interface redistribute_output_x
52 !> \public
53 module procedure redistribute_output_x_1
54 !> \public
55 module procedure redistribute_output_x_2
56 end interface
57
58 !> \public Print either vectorial or tensorial perturbation quantities of a
59 !! certain order to an output file.
60 !!
61 !! - vectorial:
62 !! - U
63 !! - DU
64 !! - tensorial:
65 !! - PV_int
66 !! - KV_int
67 !!
68 !! (The non-integrated variables are not saved because they are heavy and
69 !! not requested.)
70 !!
71 !! If \c rich_lvl is provided, \c "_R_[rich_lvl]" is appended to the data
72 !! name if it is \c >0.
73 !!
74 !! Optionally, it can be specified that this is a divided
75 !! parallel grid, corresponding to the variable \c eq_jobs_lims with index
76 !! \c eq_job_nr. In this case, the total grid size is adjusted to the one
77 !! specified by \c eq_jobs_lims and the grid is written as a subset.
78 !! \note
79 !! -# The equilibrium quantities are outputted in Flux coordinates.
80 !! -# The tensorial perturbation type can also be used for field-aligned
81 !! variables, in which case the first index is assumed to have dimension 1
82 !! only. This can be triggered using \c is_field_averaged.
83 !!
84 !! \return ierr
85 interface print_output_x
86 !> \public
87 module procedure print_output_x_1
88 !> \public
89 module procedure print_output_x_2
90 end interface
91
92contains
93 !> \private vectorial version
94 integer function calc_x_1(mds,grid_eq,grid_X,eq_1,eq_2,X,lim_sec_X) &
95 &result(ierr)
96
97 character(*), parameter :: rout_name = 'calc_X_1'
98
99 ! input / output
100 type(modes_type), intent(in) :: mds !< general modes variables
101 type(grid_type), intent(in) :: grid_eq !< equilibrium grid variables
102 type(grid_type), intent(in) :: grid_x !< perturbation grid variables
103 type(eq_1_type), intent(in) :: eq_1 !< flux equilibrium
104 type(eq_2_type), intent(in) :: eq_2 !< metric equilibrium
105 type(x_1_type), intent(inout) :: x !< vectorial perturbation variables
106 integer, intent(in), optional :: lim_sec_x(2) !< limits of \c m_X (pol. flux) or \ n_X (tor. flux)
107
108 ! initialize ierr
109 ierr = 0
110
111 ! user output
112 call writo('Calculate vectorial perturbation variables')
113 call lvl_ud(1)
114
115 ! create perturbation with modes of current X job
116 call x%init(mds,grid_x,lim_sec_x)
117
118 ! calculate U and DU
119 call writo('Calculating U and DU...')
120 call lvl_ud(1)
121 ierr = calc_u(grid_eq,grid_x,eq_1,eq_2,x)
122 chckerr('')
123 call lvl_ud(-1)
124
125 ! user output
126 call lvl_ud(-1)
127 end function calc_x_1
128 !> \private tensorial version
129 integer function calc_x_2(mds,grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,&
130 &lim_sec_X) result(ierr)
131
132 character(*), parameter :: rout_name = 'calc_X_2'
133
134 ! input / output
135 type(modes_type), intent(in) :: mds !< general modes variables
136 type(grid_type), intent(in) :: grid_eq !< equilibrium grid variables
137 type(grid_type), intent(in) :: grid_x !< perturbation grid variables
138 type(eq_1_type), intent(in) :: eq_1 !< flux equilibrium
139 type(eq_2_type), intent(in) :: eq_2 !< metric equilibrium
140 type(x_1_type), intent(inout) :: x_a !< vectorial perturbation variables (dimension 1)
141 type(x_1_type), intent(inout) :: x_b !< vectorial perturbation variables (dimension 2)
142 type(x_2_type), intent(inout) :: x !< tensorial perturbation variables
143 integer, intent(in), optional :: lim_sec_x(2,2) !< limits of \c m_X (pol flux) or \c n_X (tor flux) for both dimensions
144
145 ! initialize ierr
146 ierr = 0
147
148 ! user output
149 call writo('Calculate tensorial perturbation variables')
150 call lvl_ud(1)
151
152 ! create perturbation with modes of current X job
153 call x%init(mds,grid_x,lim_sec_x)
154
155 ! Calculate PV_i for all (k,m) pairs and n_r (equilibrium)
156 ! values of the normal coordinate
157 call writo('Calculating PV...')
158 call lvl_ud(1)
159 ierr = calc_pv(grid_eq,grid_x,eq_1,eq_2,x_a,x_b,x,lim_sec_x)
160 chckerr('')
161 call lvl_ud(-1)
162
163 ! Calculate KV_i for all (k,m) pairs and n_r (equilibrium)
164 ! values of the normal coordinate
165 call writo('Calculating KV...')
166 call lvl_ud(1)
167 ierr = calc_kv(grid_eq,grid_x,eq_1,eq_2,x_a,x_b,x,lim_sec_x)
168 chckerr('')
169 call lvl_ud(-1)
170
171 ! user output
172 call lvl_ud(-1)
173 end function calc_x_2
174
175 !> \private flux version
176 integer function redistribute_output_x_1(mds,grid,grid_out,X,X_out) &
177 &result(ierr)
178 use mpi_utilities, only: redistribute_var
179
180 character(*), parameter :: rout_name = 'redistribute_output_X_1'
181
182 ! input / output
183 type(modes_type), intent(in) :: mds !< general modes variables
184 type(grid_type), intent(in) :: grid !< perturbation grid variables
185 type(grid_type), intent(in) :: grid_out !< redistributed perturbation grid variables
186 type(x_1_type), intent(in) :: x !< vectorial perturbation variables
187 type(x_1_type), intent(inout) :: x_out !< vectorial perturbation variables in redistributed grid
188
189 ! local variables
190 integer :: ld ! counter
191 integer :: lims(2), lims_dis(2) ! limits and distributed limits, taking into account the angular extent
192 integer :: siz(3), siz_dis(3) ! size for geometric part of variable
193 real(dp), allocatable :: temp_var(:) ! temporary variable
194
195 ! initialize ierr
196 ierr = 0
197
198 ! user output
199 call writo('Redistribute vectorial perturbation variables')
200 call lvl_ud(1)
201
202 ! test
203 if (grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2)) then
204 ierr = 1
205 chckerr('grid and grid_out are not compatible')
206 end if
207
208 ! create redistributed vectorial perturbation variables
209 call x_out%init(mds,grid_out,lim_sec_x=x%lim_sec_X)
210
211 ! set up limits taking into account angular extent and temporary var
212 lims(1) = product(grid%n(1:2))*(grid%i_min-1)+1
213 lims(2) = product(grid%n(1:2))*grid%i_max
214 lims_dis(1) = product(grid%n(1:2))*(grid_out%i_min-1)+1
215 lims_dis(2) = product(grid%n(1:2))*grid_out%i_max
216 siz = [grid%n(1:2),grid%loc_n_r]
217 siz_dis = [grid%n(1:2),grid_out%loc_n_r]
218 allocate(temp_var(product(siz_dis)))
219
220 ! for all mode numbers
221 do ld = 1,size(x%U_0,4)
222 ! U_0
223 ierr = redistribute_var(reshape(rp(x%U_0(:,:,:,ld)),&
224 &[product(siz)]),temp_var,lims,lims_dis)
225 chckerr('')
226 x_out%U_0(:,:,:,ld) = &
227 &reshape(temp_var(1:product(siz_dis)),siz_dis)
228 ierr = redistribute_var(reshape(ip(x%U_0(:,:,:,ld)),&
229 &[product(siz)]),temp_var,lims,lims_dis)
230 chckerr('')
231 x_out%U_0(:,:,:,ld) = x_out%U_0(:,:,:,ld) + iu*&
232 &reshape(temp_var(1:product(siz_dis)),siz_dis)
233
234 ! U_1
235 ierr = redistribute_var(reshape(rp(x%U_1(:,:,:,ld)),&
236 &[product(siz)]),temp_var,lims,lims_dis)
237 chckerr('')
238 x_out%U_1(:,:,:,ld) = &
239 &reshape(temp_var(1:product(siz_dis)),siz_dis)
240 ierr = redistribute_var(reshape(ip(x%U_1(:,:,:,ld)),&
241 &[product(siz)]),temp_var,lims,lims_dis)
242 chckerr('')
243 x_out%U_1(:,:,:,ld) = x_out%U_1(:,:,:,ld) + iu*&
244 &reshape(temp_var(1:product(siz_dis)),siz_dis)
245
246 ! DU_0
247 ierr = redistribute_var(reshape(rp(x%DU_0(:,:,:,ld)),&
248 &[product(siz)]),temp_var,lims,lims_dis)
249 chckerr('')
250 x_out%DU_0(:,:,:,ld) = &
251 &reshape(temp_var(1:product(siz_dis)),siz_dis)
252 ierr = redistribute_var(reshape(ip(x%DU_0(:,:,:,ld)),&
253 &[product(siz)]),temp_var,lims,lims_dis)
254 chckerr('')
255 x_out%DU_0(:,:,:,ld) = x_out%DU_0(:,:,:,ld) + iu*&
256 &reshape(temp_var(1:product(siz_dis)),siz_dis)
257
258 ! DU_1
259 ierr = redistribute_var(reshape(rp(x%DU_1(:,:,:,ld)),&
260 &[product(siz)]),temp_var,lims,lims_dis)
261 chckerr('')
262 x_out%DU_1(:,:,:,ld) = &
263 &reshape(temp_var(1:product(siz_dis)),siz_dis)
264 ierr = redistribute_var(reshape(ip(x%DU_1(:,:,:,ld)),&
265 &[product(siz)]),temp_var,lims,lims_dis)
266 chckerr('')
267 x_out%DU_1(:,:,:,ld) = x_out%DU_1(:,:,:,ld) + iu*&
268 &reshape(temp_var(1:product(siz_dis)),siz_dis)
269 end do
270
271 ! user output
272 call lvl_ud(-1)
273 end function redistribute_output_x_1
274 !> \private metric version
275 integer function redistribute_output_x_2(mds,grid,grid_out,X,X_out) &
276 &result(ierr)
277 use mpi_utilities, only: redistribute_var
278
279 character(*), parameter :: rout_name = 'redistribute_output_X_2'
280
281 ! input / output
282 type(modes_type), intent(in) :: mds !< general modes variables
283 type(grid_type), intent(in) :: grid !< perturbation grid variables
284 type(grid_type), intent(in) :: grid_out !< redistributed perturbation grid variables
285 type(x_2_type), intent(in) :: x !< tensorial perturbation variables
286 type(x_2_type), intent(inout) :: x_out !< tensorial perturbation variables in redistributed grid
287
288 ! local variables
289 integer :: ld ! counters
290 integer :: n_loc(2) ! local n
291 integer :: lims(2), lims_dis(2) ! limits and distributed limits, taking into account the angular extent
292 integer :: siz(3), siz_dis(3) ! size for geometric part of variable
293 logical :: is_field_averaged ! whether field-averaged, so only one dimension for first index
294 real(dp), allocatable :: temp_var(:) ! temporary variable
295
296 ! initialize ierr
297 ierr = 0
298
299 ! user output
300 call writo('Redistribute tensorial perturbation variables')
301 call lvl_ud(1)
302
303 ! test
304 is_field_averaged = size(x%PV_0,1).eq.1
305 if ((grid%n(1).ne.grid_out%n(1) .or. grid%n(2).ne.grid_out%n(2)) .and. &
306 &.not.is_field_averaged) then ! for field-averaged grids, don't do tests
307 ierr = 1
308 chckerr('grid and grid_out are not compatible')
309 end if
310
311 ! create redistributed tensorial perturbation variables
312 call x_out%init(mds,grid_out,lim_sec_x=x%lim_sec_X,&
313 &is_field_averaged=is_field_averaged)
314
315 ! set up limits taking into account angular extent and temporary var
316 n_loc = grid%n(1:2)
317 if (is_field_averaged) n_loc(1) = 1
318 lims(1) = product(n_loc)*(grid%i_min-1)+1
319 lims(2) = product(n_loc)*grid%i_max
320 lims_dis(1) = product(n_loc)*(grid_out%i_min-1)+1
321 lims_dis(2) = product(n_loc)*grid_out%i_max
322 siz = [n_loc,grid%loc_n_r]
323 siz_dis = [n_loc,grid_out%loc_n_r]
324 allocate(temp_var(product(siz_dis)))
325
326 ! for all mode number combinations, symmetric
327 do ld = 1,size(x%PV_0,4)
328 ! PV_0
329 ierr = redistribute_var(reshape(rp(x%PV_0(1:n_loc(1),:,:,ld)),&
330 &[product(siz)]),temp_var,lims,lims_dis)
331 chckerr('')
332 x_out%PV_0(1:n_loc(1),:,:,ld) = &
333 &reshape(temp_var(1:product(siz_dis)),siz_dis)
334 ierr = redistribute_var(reshape(ip(x%PV_0(1:n_loc(1),:,:,ld)),&
335 &[product(siz)]),temp_var,lims,lims_dis)
336 chckerr('')
337 x_out%PV_0(1:n_loc(1),:,:,ld) = x_out%PV_0(1:n_loc(1),:,:,ld) + iu*&
338 &reshape(temp_var(1:product(siz_dis)),siz_dis)
339
340 ! PV_2
341 ierr = redistribute_var(reshape(rp(x%PV_2(1:n_loc(1),:,:,ld)),&
342 &[product(siz)]),temp_var,lims,lims_dis)
343 chckerr('')
344 x_out%PV_2(1:n_loc(1),:,:,ld) = &
345 &reshape(temp_var(1:product(siz_dis)),siz_dis)
346 ierr = redistribute_var(reshape(ip(x%PV_2(1:n_loc(1),:,:,ld)),&
347 &[product(siz)]),temp_var,lims,lims_dis)
348 chckerr('')
349 x_out%PV_2(1:n_loc(1),:,:,ld) = x_out%PV_2(1:n_loc(1),:,:,ld) + iu*&
350 &reshape(temp_var(1:product(siz_dis)),siz_dis)
351
352 ! KV_0
353 ierr = redistribute_var(reshape(rp(x%KV_0(1:n_loc(1),:,:,ld)),&
354 &[product(siz)]),temp_var,lims,lims_dis)
355 chckerr('')
356 x_out%KV_0(1:n_loc(1),:,:,ld) = &
357 &reshape(temp_var(1:product(siz_dis)),siz_dis)
358 ierr = redistribute_var(reshape(ip(x%KV_0(1:n_loc(1),:,:,ld)),&
359 &[product(siz)]),temp_var,lims,lims_dis)
360 chckerr('')
361 x_out%KV_0(1:n_loc(1),:,:,ld) = x_out%KV_0(1:n_loc(1),:,:,ld) + iu*&
362 &reshape(temp_var(1:product(siz_dis)),siz_dis)
363
364 ! KV_2
365 ierr = redistribute_var(reshape(rp(x%KV_2(1:n_loc(1),:,:,ld)),&
366 &[product(siz)]),temp_var,lims,lims_dis)
367 chckerr('')
368 x_out%KV_2(1:n_loc(1),:,:,ld) = &
369 &reshape(temp_var(1:product(siz_dis)),siz_dis)
370 ierr = redistribute_var(reshape(ip(x%KV_2(1:n_loc(1),:,:,ld)),&
371 &[product(siz)]),temp_var,lims,lims_dis)
372 chckerr('')
373 x_out%KV_2(1:n_loc(1),:,:,ld) = x_out%KV_2(1:n_loc(1),:,:,ld) + iu*&
374 &reshape(temp_var(1:product(siz_dis)),siz_dis)
375 end do
376
377 ! for all mode number combinations, asymmetric
378 do ld = 1,size(x%PV_1,4)
379 ! PV_1
380 ierr = redistribute_var(reshape(rp(x%PV_1(1:n_loc(1),:,:,ld)),&
381 &[product(siz)]),temp_var,lims,lims_dis)
382 chckerr('')
383 x_out%PV_1(1:n_loc(1),:,:,ld) = &
384 &reshape(temp_var(1:product(siz_dis)),siz_dis)
385 ierr = redistribute_var(reshape(ip(x%PV_1(1:n_loc(1),:,:,ld)),&
386 &[product(siz)]),temp_var,lims,lims_dis)
387 chckerr('')
388 x_out%PV_1(1:n_loc(1),:,:,ld) = x_out%PV_1(1:n_loc(1),:,:,ld) + iu*&
389 &reshape(temp_var(1:product(siz_dis)),siz_dis)
390
391 ! KV_1
392 ierr = redistribute_var(reshape(rp(x%KV_1(1:n_loc(1),:,:,ld)),&
393 &[product(siz)]),temp_var,lims,lims_dis)
394 chckerr('')
395 x_out%KV_1(1:n_loc(1),:,:,ld) = &
396 &reshape(temp_var(1:product(siz_dis)),siz_dis)
397 ierr = redistribute_var(reshape(ip(x%KV_1(1:n_loc(1),:,:,ld)),&
398 &[product(siz)]),temp_var,lims,lims_dis)
399 chckerr('')
400 x_out%KV_1(1:n_loc(1),:,:,ld) = x_out%KV_1(1:n_loc(1),:,:,ld) + iu*&
401 &reshape(temp_var(1:product(siz_dis)),siz_dis)
402 end do
403
404 ! user output
405 call lvl_ud(-1)
406 end function redistribute_output_x_2
407
408 !> \private vectorial version.
409 integer function print_output_x_1(grid_X,X,data_name,rich_lvl,par_div,&
410 &lim_sec_X) result(ierr)
411
412 use num_vars, only: pb3d_name, eq_jobs_lims, eq_job_nr
413 use grid_utilities, only: trim_grid
414 use hdf5_ops, only: print_hdf5_arrs
415 use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
416 &max_dim_var_1d
417 use x_vars, only: n_mod_x
418
419 character(*), parameter :: rout_name = 'print_output_X_1'
420
421 ! input / output
422 type(grid_type), intent(in) :: grid_x !< perturbation grid variables
423 type(x_1_type), intent(in) :: x !< vectorial perturbation variables
424 character(len=*), intent(in) :: data_name !< name under which to store
425 integer, intent(in), optional :: rich_lvl !< Richardson level to print
426 logical, intent(in), optional :: par_div !< is a parallely divided grid
427 integer, intent(in), optional :: lim_sec_x(2) !< limits of \c m_X (pol. flux) or \c n_X (tor. flux)
428
429 ! local variables
430 type(grid_type) :: grid_trim ! trimmed grid
431 integer :: n_tot(3) ! total n
432 integer :: par_id(2) ! parallel interval
433 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
434 type(var_1d_type), allocatable, target :: x_1d(:) ! 1D equivalent of X variables
435 type(var_1d_type), pointer :: x_1d_loc => null() ! local element in X_1D
436 integer :: id ! counters
437 logical :: par_div_loc ! local par_div
438 integer :: lim_sec_x_loc(2) ! local lim_sec_X
439 integer :: loc_size ! local size
440
441 ! initialize ierr
442 ierr = 0
443
444 ! user output
445 call writo('Write vectorial perturbation variables to output file')
446 call lvl_ud(1)
447
448 ! trim grid
449 ierr = trim_grid(grid_x,grid_trim,norm_id)
450 chckerr('')
451
452 ! set local par_div
453 par_div_loc = .false.
454 if (present(par_div)) par_div_loc = par_div
455
456 ! set total n and parallel interval
457 n_tot = grid_trim%n
458 par_id = [1,n_tot(1)]
459 if (grid_trim%n(1).gt.0 .and. par_div_loc) then ! total grid includes all equilibrium jobs
460 n_tot(1) = maxval(eq_jobs_lims)-minval(eq_jobs_lims)+1
461 par_id = eq_jobs_lims(:,eq_job_nr)
462 end if
463
464 ! set local size and lim_sec_X
465 loc_size = size(x%U_0(:,:,norm_id(1):norm_id(2),:))
466 lim_sec_x_loc = [1,n_mod_x]
467 if (present(lim_sec_x)) lim_sec_x_loc = lim_sec_x
468
469 ! Set up the 1D equivalents of the perturbation variables
470 allocate(x_1d(max_dim_var_1d))
471
472 ! set up variables X_1D
473 id = 1
474
475 ! RE_U_0
476 x_1d_loc => x_1d(id); id = id+1
477 x_1d_loc%var_name = 'RE_U_0'
478 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
479 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
480 x_1d_loc%tot_i_min = [1,1,1,1]
481 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
482 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
483 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
484 &lim_sec_x_loc(2)]
485 allocate(x_1d_loc%p(loc_size))
486 x_1d_loc%p = reshape(rp(x%U_0(:,:,norm_id(1):norm_id(2),:)),&
487 &[loc_size])
488
489 ! IM_U_0
490 x_1d_loc => x_1d(id); id = id+1
491 x_1d_loc%var_name = 'IM_U_0'
492 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
493 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
494 x_1d_loc%tot_i_min = [1,1,1,1]
495 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
496 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
497 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
498 &lim_sec_x_loc(2)]
499 allocate(x_1d_loc%p(loc_size))
500 x_1d_loc%p = reshape(ip(x%U_0(:,:,norm_id(1):norm_id(2),:)),&
501 &[loc_size])
502
503 ! RE_U_1
504 x_1d_loc => x_1d(id); id = id+1
505 x_1d_loc%var_name = 'RE_U_1'
506 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
507 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
508 x_1d_loc%tot_i_min = [1,1,1,1]
509 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
510 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
511 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
512 &lim_sec_x_loc(2)]
513 allocate(x_1d_loc%p(loc_size))
514 x_1d_loc%p = reshape(rp(x%U_1(:,:,norm_id(1):norm_id(2),:)),&
515 &[loc_size])
516
517 ! IM_U_1
518 x_1d_loc => x_1d(id); id = id+1
519 x_1d_loc%var_name = 'IM_U_1'
520 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
521 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
522 x_1d_loc%tot_i_min = [1,1,1,1]
523 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
524 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
525 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
526 &lim_sec_x_loc(2)]
527 allocate(x_1d_loc%p(loc_size))
528 x_1d_loc%p = reshape(ip(x%U_1(:,:,norm_id(1):norm_id(2),:)),&
529 &[loc_size])
530
531 ! RE_DU_0
532 x_1d_loc => x_1d(id); id = id+1
533 x_1d_loc%var_name = 'RE_DU_0'
534 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
535 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
536 x_1d_loc%tot_i_min = [1,1,1,1]
537 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
538 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
539 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
540 &lim_sec_x_loc(2)]
541 allocate(x_1d_loc%p(loc_size))
542 x_1d_loc%p = reshape(rp(x%DU_0(:,:,norm_id(1):norm_id(2),:)),&
543 &[loc_size])
544
545 ! IM_DU_0
546 x_1d_loc => x_1d(id); id = id+1
547 x_1d_loc%var_name = 'IM_DU_0'
548 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
549 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
550 x_1d_loc%tot_i_min = [1,1,1,1]
551 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
552 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
553 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
554 &lim_sec_x_loc(2)]
555 allocate(x_1d_loc%p(loc_size))
556 x_1d_loc%p = reshape(ip(x%DU_0(:,:,norm_id(1):norm_id(2),:)),&
557 &[loc_size])
558
559 ! RE_DU_1
560 x_1d_loc => x_1d(id); id = id+1
561 x_1d_loc%var_name = 'RE_DU_1'
562 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
563 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
564 x_1d_loc%tot_i_min = [1,1,1,1]
565 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
566 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
567 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
568 &lim_sec_x_loc(2)]
569 allocate(x_1d_loc%p(loc_size))
570 x_1d_loc%p = reshape(rp(x%DU_1(:,:,norm_id(1):norm_id(2),:)),&
571 &[loc_size])
572
573 ! IM_DU_1
574 x_1d_loc => x_1d(id); id = id+1
575 x_1d_loc%var_name = 'IM_DU_1'
576 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
577 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
578 x_1d_loc%tot_i_min = [1,1,1,1]
579 x_1d_loc%tot_i_max = [n_tot,n_mod_x]
580 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,lim_sec_x_loc(1)]
581 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
582 &lim_sec_x_loc(2)]
583 allocate(x_1d_loc%p(loc_size))
584 x_1d_loc%p = reshape(ip(x%DU_1(:,:,norm_id(1):norm_id(2),:)),&
585 &[loc_size])
586
587 ! write
588 ierr = print_hdf5_arrs(x_1d(1:id-1),pb3d_name,trim(data_name),&
589 &rich_lvl=rich_lvl,ind_print=.not.grid_trim%divided)
590 chckerr('')
591
592 ! clean up
593 call grid_trim%dealloc()
594 call dealloc_var_1d(x_1d)
595 nullify(x_1d_loc)
596
597 ! user output
598 call lvl_ud(-1)
599 end function print_output_x_1
600 !> \private tensorial version
601 integer function print_output_x_2(grid_X,X,data_name,rich_lvl,par_div,&
602 &lim_sec_X,is_field_averaged) result(ierr)
603 use num_vars, only: pb3d_name, eq_jobs_lims, eq_job_nr
604 use grid_utilities, only: trim_grid
605 use hdf5_ops, only: print_hdf5_arrs
606 use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
607 &max_dim_var_1d
608 use x_utilities, only: is_necessary_x, get_sec_x_range
609 use x_vars, only: set_nn_mod
610 use num_utilities, only: c
611
612 character(*), parameter :: rout_name = 'print_output_X_2'
613
614 ! input / output
615 type(grid_type), intent(in) :: grid_x !< perturbation grid variables
616 type(x_2_type), intent(in) :: x !< tensorial perturbation variables
617 character(len=*), intent(in) :: data_name !< name under which to store
618 integer, intent(in), optional :: rich_lvl !< Richardson level to print
619 logical, intent(in), optional :: par_div !< is a parallely divided grid
620 integer, intent(in), optional :: lim_sec_x(2,2) !< limits of \c m_X (pol. flux) or \c n_X (tor. flux)
621 logical, intent(in), optional :: is_field_averaged !< whether field-averaged, so only one dimension for first index
622
623 ! local variables
624 type(grid_type) :: grid_trim ! trimmed grid
625 integer :: n_tot(3) ! total n
626 integer :: par_id(2) ! parallel interval
627 integer :: par_id_loc(2) ! local parallel interval
628 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
629 type(var_1d_type), allocatable, target :: x_1d(:) ! 1D equivalent of X variables
630 type(var_1d_type), pointer :: x_1d_loc => null() ! local element in X_1D
631 logical :: print_this(2) ! whether symmetric and asymmetric variables need to be printed
632 integer :: nn_mod_tot(2) ! total nr. of modes for symmetric and asymmetric variables
633 integer :: nn_mod_loc(2) ! local nr. of modes for symmetric and asymmetric variables
634 integer :: id ! counter
635 logical :: par_div_loc ! local par_div
636 integer :: loc_size(2) ! local size for symmetric and asymmetric variables
637 integer :: m, k ! counters
638 integer :: sxr_loc(2,2) ! local secondary X limits for symmetric and asymmetric variables
639 integer :: sxr_tot(2,2) ! total secondary X limits for symmetric and asymmetric variables
640
641 ! initialize ierr
642 ierr = 0
643
644 ! user output
645 call writo('Write tensorial perturbation variables to output file')
646 call lvl_ud(1)
647
648 ! trim grid
649 ierr = trim_grid(grid_x,grid_trim,norm_id)
650 chckerr('')
651
652 ! set local par_div
653 par_div_loc = .false.
654 if (present(par_div)) par_div_loc = par_div
655
656 ! set total n, parallel interval and total n_mod_X
657 n_tot = grid_trim%n
658 par_id = [1,n_tot(1)]
659 par_id_loc = [1,n_tot(1)]
660 if (grid_trim%n(1).gt.0 .and. par_div_loc) then ! total grid includes all equilibrium jobs
661 n_tot(1) = maxval(eq_jobs_lims)-minval(eq_jobs_lims)+1
662 par_id = eq_jobs_lims(:,eq_job_nr)
663 par_id_loc = par_id - par_id(1)+1
664 else if (present(is_field_averaged)) then
665 if (is_field_averaged) then ! only first point
666 par_id = [1,1]
667 par_id_loc = [1,1]
668 n_tot(1) = 1
669 end if
670 end if
671 nn_mod_tot(1) = set_nn_mod(.true.)
672 nn_mod_tot(2) = set_nn_mod(.false.)
673
674 ! set dimensions, parallel limits and local and total nn_mod_X
675
676 ! Set up the 1D equivalents of the perturbation variables
677 allocate(x_1d(max_dim_var_1d))
678
679 ! set up variables X_1D
680 id = 1
681
682 ! loop over modes of dimension 2
683 do m = 1,x%n_mod(2)
684 ! get contiguous range of modes of this m
685 call get_sec_x_range(sxr_loc(:,1),sxr_tot(:,1),m,.true.,lim_sec_x)
686 call get_sec_x_range(sxr_loc(:,2),sxr_tot(:,2),m,.false.,lim_sec_x)
687 nn_mod_loc = sxr_loc(2,:)-sxr_loc(1,:)+1
688 print_this = .false.
689 do k = 1,2
690 if (sxr_loc(1,k).le.sxr_loc(2,k)) print_this(k) = .true. ! a bound is found
691 end do
692
693 ! set local size
694 loc_size(1) = (par_id(2)-par_id(1)+1)*n_tot(2)*&
695 &(norm_id(2)-norm_id(1)+1)*nn_mod_loc(1)
696 loc_size(2) = (par_id(2)-par_id(1)+1)*n_tot(2)*&
697 &(norm_id(2)-norm_id(1)+1)*nn_mod_loc(2)
698
699 if (print_this(1)) then
700 ! RE_PV_0
701 x_1d_loc => x_1d(id); id = id+1
702 x_1d_loc%var_name = 'RE_PV_0'
703 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
704 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
705 x_1d_loc%tot_i_min = [1,1,1,1]
706 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
707 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
708 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
709 &sxr_tot(2,1)]
710 allocate(x_1d_loc%p(loc_size(1)))
711 x_1d_loc%p = reshape(rp(x%PV_0(par_id_loc(1):par_id_loc(2),:,&
712 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
713 &[loc_size(1)])
714
715 ! IM_PV_0
716 x_1d_loc => x_1d(id); id = id+1
717 x_1d_loc%var_name = 'IM_PV_0'
718 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
719 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
720 x_1d_loc%tot_i_min = [1,1,1,1]
721 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
722 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
723 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
724 &sxr_tot(2,1)]
725 allocate(x_1d_loc%p(loc_size(1)))
726 x_1d_loc%p = reshape(ip(x%PV_0(par_id_loc(1):par_id_loc(2),:,&
727 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
728 &[loc_size(1)])
729
730 ! RE_PV_2
731 x_1d_loc => x_1d(id); id = id+1
732 x_1d_loc%var_name = 'RE_PV_2'
733 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
734 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
735 x_1d_loc%tot_i_min = [1,1,1,1]
736 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
737 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
738 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
739 &sxr_tot(2,1)]
740 allocate(x_1d_loc%p(loc_size(1)))
741 x_1d_loc%p = reshape(rp(x%PV_2(par_id_loc(1):par_id_loc(2),:,&
742 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
743 &[loc_size(1)])
744
745 ! IM_PV_2
746 x_1d_loc => x_1d(id); id = id+1
747 x_1d_loc%var_name = 'IM_PV_2'
748 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
749 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
750 x_1d_loc%tot_i_min = [1,1,1,1]
751 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
752 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
753 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
754 &sxr_tot(2,1)]
755 allocate(x_1d_loc%p(loc_size(1)))
756 x_1d_loc%p = reshape(ip(x%PV_2(par_id_loc(1):par_id_loc(2),:,&
757 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
758 &[loc_size(1)])
759
760 ! RE_KV_0
761 x_1d_loc => x_1d(id); id = id+1
762 x_1d_loc%var_name = 'RE_KV_0'
763 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
764 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
765 x_1d_loc%tot_i_min = [1,1,1,1]
766 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
767 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
768 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
769 &sxr_tot(2,1)]
770 allocate(x_1d_loc%p(loc_size(1)))
771 x_1d_loc%p = reshape(rp(x%KV_0(par_id_loc(1):par_id_loc(2),:,&
772 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
773 &[loc_size(1)])
774
775 ! IM_KV_0
776 x_1d_loc => x_1d(id); id = id+1
777 x_1d_loc%var_name = 'IM_KV_0'
778 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
779 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
780 x_1d_loc%tot_i_min = [1,1,1,1]
781 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
782 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
783 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
784 &sxr_tot(2,1)]
785 allocate(x_1d_loc%p(loc_size(1)))
786 x_1d_loc%p = reshape(ip(x%KV_0(par_id_loc(1):par_id_loc(2),:,&
787 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
788 &[loc_size(1)])
789
790 ! RE_KV_2
791 x_1d_loc => x_1d(id); id = id+1
792 x_1d_loc%var_name = 'RE_KV_2'
793 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
794 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
795 x_1d_loc%tot_i_min = [1,1,1,1]
796 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
797 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
798 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
799 &sxr_tot(2,1)]
800 allocate(x_1d_loc%p(loc_size(1)))
801 x_1d_loc%p = reshape(rp(x%KV_2(par_id_loc(1):par_id_loc(2),:,&
802 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
803 &[loc_size(1)])
804
805 ! IM_KV_2
806 x_1d_loc => x_1d(id); id = id+1
807 x_1d_loc%var_name = 'IM_KV_2'
808 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
809 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
810 x_1d_loc%tot_i_min = [1,1,1,1]
811 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(1)]
812 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,1)]
813 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
814 &sxr_tot(2,1)]
815 allocate(x_1d_loc%p(loc_size(1)))
816 x_1d_loc%p = reshape(ip(x%KV_2(par_id_loc(1):par_id_loc(2),:,&
817 &norm_id(1):norm_id(2),sxr_loc(1,1):sxr_loc(2,1))),&
818 &[loc_size(1)])
819 end if
820
821 if (print_this(2)) then
822 ! RE_PV_1
823 x_1d_loc => x_1d(id); id = id+1
824 x_1d_loc%var_name = 'RE_PV_1'
825 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
826 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
827 x_1d_loc%tot_i_min = [1,1,1,1]
828 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
829 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
830 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
831 &sxr_tot(2,2)]
832 allocate(x_1d_loc%p(loc_size(2)))
833 x_1d_loc%p = reshape(rp(x%PV_1(par_id_loc(1):par_id_loc(2),:,&
834 &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
835 &[loc_size(2)])
836
837 ! IM_PV_1
838 x_1d_loc => x_1d(id); id = id+1
839 x_1d_loc%var_name = 'IM_PV_1'
840 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
841 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
842 x_1d_loc%tot_i_min = [1,1,1,1]
843 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
844 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
845 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
846 &sxr_tot(2,2)]
847 allocate(x_1d_loc%p(loc_size(2)))
848 x_1d_loc%p = reshape(ip(x%PV_1(par_id_loc(1):par_id_loc(2),:,&
849 &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
850 &[loc_size(2)])
851
852 ! RE_KV_1
853 x_1d_loc => x_1d(id); id = id+1
854 x_1d_loc%var_name = 'RE_KV_1'
855 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
856 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
857 x_1d_loc%tot_i_min = [1,1,1,1]
858 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
859 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
860 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
861 &sxr_tot(2,2)]
862 allocate(x_1d_loc%p(loc_size(2)))
863 x_1d_loc%p = reshape(rp(x%KV_1(par_id_loc(1):par_id_loc(2),:,&
864 &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
865 &[loc_size(2)])
866
867 ! IM_KV_1
868 x_1d_loc => x_1d(id); id = id+1
869 x_1d_loc%var_name = 'IM_KV_1'
870 allocate(x_1d_loc%tot_i_min(4),x_1d_loc%tot_i_max(4))
871 allocate(x_1d_loc%loc_i_min(4),x_1d_loc%loc_i_max(4))
872 x_1d_loc%tot_i_min = [1,1,1,1]
873 x_1d_loc%tot_i_max = [n_tot,nn_mod_tot(2)]
874 x_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min,sxr_tot(1,2)]
875 x_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max,&
876 &sxr_tot(2,2)]
877 allocate(x_1d_loc%p(loc_size(2)))
878 x_1d_loc%p = reshape(ip(x%KV_1(par_id_loc(1):par_id_loc(2),:,&
879 &norm_id(1):norm_id(2),sxr_loc(1,2):sxr_loc(2,2))),&
880 &[loc_size(2)])
881 end if
882 end do
883
884 ! write
885 ierr = print_hdf5_arrs(x_1d(1:id-1),pb3d_name,trim(data_name),&
886 &rich_lvl=rich_lvl,ind_print=.not.grid_trim%divided)
887 chckerr('')
888
889 ! clean up
890 call grid_trim%dealloc()
891 call dealloc_var_1d(x_1d)
892 nullify(x_1d_loc)
893
894 ! user output
895 call lvl_ud(-1)
896 end function print_output_x_2
897
898 !> Initializes some variables concerning the mode numbers.
899 !!
900 !! Setup minimum and maximum of mode numbers at every flux surface in
901 !! equilibrium coordinates
902 !! - \c min_n_X, \c max_n_X
903 !! - \c min_m_X, \c max_m_X,
904 !!
905 !! It functions depending on the \c X_style used: 1 (prescribed) or 2
906 !! (fast). For the fast style, at every flux surface the range of modes is
907 !! sought that resonates most:
908 !! - \f$m = nq \pm \frac{N}{2}\f$ for poloidal flux
909 !! - \f$n = \iota m \pm \frac{N}{2}\f$ for toroidal flux
910 !!
911 !! where \f$N\f$ = \c n_mod_X.
912 !!
913 !! However, at the same time, both \c m and \c n have to be larger, in
914 !! absolute value, than \c min_sec_X. Therefore, the range of width \c
915 !! n_mod_X can be shifted upwards.
916 !!
917 !! \return ierr
918 integer function init_modes(grid_eq,eq) result(ierr)
919 use num_vars, only: use_pol_flux_f, x_style
920 use x_vars, only: prim_x, min_sec_x, max_sec_x, n_mod_x, min_n_x, &
921 &max_n_x, min_m_x, max_m_x, min_nm_x
922 use mpi_utilities, only: get_ser_var
923 use grid_utilities, only: trim_grid
924
925 character(*), parameter :: rout_name = 'init_modes'
926
927 ! input / output
928 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
929 type(eq_1_type), intent(in) :: eq !< flux equilibrium
930
931 ! local variables
932 type(grid_type) :: grid_eq_trim ! trimmed version of equilibrium
933 real(dp), allocatable :: jq_tot(:) ! saf. fac. or rot. transf. and derivs. in Flux coords.
934 integer :: norm_eq_id(2) ! untrimmed normal indices for trimmed grids
935
936 ! initialize ierr
937 ierr = 0
938
939 ! user output
940 call writo('Initializing perturbation modes variables')
941 call lvl_ud(1)
942
943 ! get trimmed grids
944 ierr = trim_grid(grid_eq,grid_eq_trim,norm_eq_id)
945 chckerr('')
946
947 ! (re)allocate variables
948 if (allocated(min_n_x)) deallocate(min_n_x)
949 if (allocated(max_n_x)) deallocate(max_n_x)
950 if (allocated(min_m_x)) deallocate(min_m_x)
951 if (allocated(max_m_x)) deallocate(max_m_x)
952 allocate(min_n_x(grid_eq_trim%n(3)),max_n_x(grid_eq_trim%n(3)))
953 allocate(min_m_x(grid_eq_trim%n(3)),max_m_x(grid_eq_trim%n(3)))
954
955 ! get serial version of safety factor or rot. transform
956 allocate(jq_tot(grid_eq_trim%n(3)))
957 if (use_pol_flux_f) then
958 if (grid_eq%divided) then
959 ierr = get_ser_var(eq%q_saf_FD(norm_eq_id(1):norm_eq_id(2),0),&
960 &jq_tot,scatter=.true.)
961 chckerr('')
962 else
963 jq_tot = eq%q_saf_FD(:,0)
964 end if
965 else
966 if (grid_eq%divided) then
967 ierr = get_ser_var(eq%rot_t_FD(norm_eq_id(1):norm_eq_id(2),0),&
968 &jq_tot,scatter=.true.)
969 chckerr('')
970 else
971 jq_tot = eq%rot_t_FD(:,0)
972 end if
973 end if
974
975 ! set the limits depending on the X style
976 select case (x_style)
977 case (1) ! prescribed
978 if (use_pol_flux_f) then
979 min_n_x = prim_x
980 max_n_x = prim_x
981 min_m_x = min_sec_x
982 max_m_x = max_sec_x
983 else
984 min_n_x = min_sec_x
985 max_n_x = max_sec_x
986 min_m_x = prim_x
987 max_m_x = prim_x
988 end if
989 case (2) ! fast
990 if (use_pol_flux_f) then
991 min_n_x = prim_x
992 max_n_x = prim_x
993 if (prim_x*jq_tot(1).gt.0) then ! nq > 0: m > 0
994 min_m_x = nint(prim_x*jq_tot-(n_mod_x-1)*0.5)
995 min_m_x = max(min_nm_x,min_m_x)
996 max_m_x = min_m_x + n_mod_x - 1
997 else ! nq < 0: m < 0
998 max_m_x = nint(prim_x*jq_tot+(n_mod_x-1)*0.5)
999 max_m_x = min(-min_nm_x,max_m_x)
1000 min_m_x = max_m_x - n_mod_x + 1
1001 end if
1002 else
1003 if (prim_x*jq_tot(1).gt.0) then ! m iota > 0: n > 0
1004 min_n_x = nint(prim_x*jq_tot-(n_mod_x-1)*0.5)
1005 min_n_x = max(min_nm_x,min_n_x)
1006 max_n_x = min_n_x + n_mod_x - 1
1007 else ! m iota < 0: n < 0
1008 max_n_x = nint(prim_x*jq_tot+(n_mod_x-1)*0.5)
1009 max_n_x = min(-min_nm_x,max_n_x)
1010 min_n_x = max_n_x - n_mod_x + 1
1011 end if
1012 min_m_x = prim_x
1013 max_m_x = prim_x
1014 end if
1015 end select
1016
1017 ! clean up
1018 call grid_eq_trim%dealloc()
1019
1020 call lvl_ud(-1)
1021 call writo('Perturbation modes variables initialized')
1022 end function init_modes
1023
1024 !> Sets up some variables concerning the mode numbers.
1025 !!
1026 !! Apart from the mode numbers
1027 !! - \c n, \c m,
1028 !!
1029 !! also the variables of the secondary modes
1030 !! - sec,
1031 !!
1032 !! are set up in the coordinates of the grid passed. The normal values of
1033 !! this grid are saved as well, in
1034 !! - r_F.
1035 !!
1036 !! All variables are tabulated in the full normal grid. The mode indices,
1037 !! however, are chosen so that there is maximum overlap between different
1038 !! normal positions. This is trivial for \c X_style 1 (prescribed), but for
1039 !! \c X_style 2 (fast), this is best explained by an example:
1040 !!
1041 !! kd m1 m2 m3
1042 !! ------------
1043 !! 1 10 11 12 <- start
1044 !! 2 10 11 12 <- no change in limits
1045 !! 2 13 11 12 <- limits shift up by one
1046 !! 3 13 11 12 <- no change in limits
1047 !! 4 13 14 12 <- limits shift up by one
1048 !! 5 16 14 15 <- limits shift up by two
1049 !! 6 13 14 15 <- limits shift down by one
1050 !!
1051 !! This procedure makes use of the global variables
1052 !! - \c min_m_X, max_m_X, min_n_X, max_m_X
1053 !!
1054 !! that have to be set up using init_nm_x().
1055 !!
1056 !! Optionally, \c n and \c m can be plot.
1057 !!
1058 !! \return ierr
1059 integer function setup_modes(mds,grid_eq,grid,plot_name) result(ierr)
1060 use num_vars, only: use_pol_flux_f, rank
1061 use x_vars, only: n_mod_x, min_n_x, min_m_x, max_n_x, max_m_x
1062 use grid_utilities, only: trim_grid
1063 use eq_vars, only: max_flux_f
1064 use num_utilities, only: spline
1065
1066 character(*), parameter :: rout_name = 'setup_modes'
1067
1068 ! input / output
1069 type(modes_type), intent(inout), target :: mds !< modes variables
1070 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
1071 type(grid_type), intent(in) :: grid !< grid at which to calculate modes
1072 character(len=*), intent(in), optional :: plot_name !< name to be used when plotting \c n and \c m
1073
1074 ! local variables
1075 real(dp), allocatable :: min_nm_x(:) ! interpolated minimum mode number
1076 real(dp), allocatable :: x_plot(:,:) ! abscissa of plot
1077 integer, pointer :: nm_x(:,:) ! either n or m
1078 integer :: id, ld, kd ! counters
1079 integer :: ld_loc ! shfited ld
1080 integer :: n_mod_tot ! total number of modes
1081 integer :: ind_id ! current size of ind_tot
1082 integer :: ld_shift ! shift in table indices
1083 integer :: mod_x_range ! total mode range
1084 integer :: delta_ld ! change in total mode numbers
1085 integer, allocatable :: ind_cur(:) ! current indices
1086 integer, allocatable :: ind_tot(:,:) ! total index information, will be later cut to sec
1087 character(len=max_str_ln), allocatable :: plot_titles(:) ! title for plots
1088 character(len=max_str_ln) :: plot_name_tot ! file name for plots
1089 character(len=max_str_ln) :: err_msg ! error message
1090#if ldebug
1091 real(dp), allocatable :: y_plot(:,:) ! ordinate of plot
1092#endif
1093
1094 ! initialize ierr
1095 ierr = 0
1096
1097 ! user output
1098 call writo('Setting up general modes variables')
1099 call lvl_ud(1)
1100
1101 ! set up normal tabulation values
1102 allocate(ind_tot(-grid%n(3)*n_mod_x:grid%n(3)*n_mod_x,4)) ! not quite absolute maximum but should never be reached
1103 allocate(ind_cur(n_mod_x)) ! indices of total modes currently being treated
1104
1105 ! calculate n and m
1106 allocate(mds%n(grid%n(3),n_mod_x))
1107 allocate(mds%m(grid%n(3),n_mod_x))
1108
1109 ! auxiliary variable
1110 allocate(min_nm_x(grid%n(3)))
1111
1112 ! iterate over n (id=1) and m (id=2)
1113 do id = 1,2
1114 ! point nm_X, and interpolate mode minimum (linear, to preserve
1115 ! form)
1116 if (id.eq.1) then
1117 nm_x => mds%n
1118 mod_x_range = max_n_x(1)-min_n_x(1) ! assumes that this is true for all normal points
1119 ierr = spline(grid_eq%r_F,min_n_x*1._dp,grid%r_F,min_nm_x,ord=1)
1120 chckerr('')
1121 else
1122 nm_x => mds%m
1123 mod_x_range = max_m_x(1)-min_m_x(1) ! assumes that this is true for all normal points
1124 ierr = spline(grid_eq%r_F,min_m_x*1._dp,grid%r_F,min_nm_x,ord=1)
1125 chckerr('')
1126 end if
1127
1128 !!! For non-monotonous tests:
1129 !!!if (id .eq. 2) then
1130 !!!min_nm_X(nint(grid%n(3)*0.7):grid%n(3)) = &
1131 !!!&2*min_nm_X(nint(grid%n(3)*0.7)) - &
1132 !!!&min_nm_X(nint(grid%n(3)*0.7):grid%n(3))
1133 !!!end if
1134
1135 ! initialize variables for first normal grid point
1136 ld_shift = 0
1137 ind_id = 0
1138 do ld = 1,n_mod_x
1139 nm_x(1,ld) = nint(min_nm_x(1)) + mod_x_range*(ld-1)/(n_mod_x-1)
1140 ind_tot(ld,:) = [nm_x(1,ld),1,0,ld]
1141 ind_id = ind_id + 1
1142#if ldebug
1143 if (debug_setup_modes) write(*,*) 'starting with', ld, ':', &
1144 &ind_tot(ld,:)
1145#endif
1146 end do
1147 ind_cur = [(ld, ld=1,n_mod_x)]
1148
1149 ! iterate over all next normal grid points
1150 do kd = 2,grid%n(3)
1151 ! update ld delta and shift
1152 delta_ld = nint(min_nm_x(kd)) - nint(min_nm_x(kd-1))
1153 ld_shift = ld_shift + delta_ld
1154
1155 ! set n or m
1156 do ld = 1,n_mod_x
1157 ! shift ld and wrap back to [1..n_mod_X]
1158 ld_loc = modulo(ld+ld_shift-1,n_mod_x)+1
1159 nm_x(kd,ld_loc) = nint(min_nm_x(kd)) + &
1160 &mod_x_range*(ld-1)/(n_mod_x-1)
1161 end do
1162
1163 ! if there was a shift, set new total index information
1164 if (abs(delta_ld).gt.0) then
1165#if ldebug
1166 if (debug_setup_modes) write(*,*) 'kd', kd, 'delta', &
1167 &delta_ld
1168#endif
1169 do ld = 1,abs(delta_ld)
1170 if (delta_ld.gt.0) then ! mode number has increased
1171 ind_tot(ind_cur(1),3) = kd - 1 ! upper limit in normal range
1172 ind_tot(ind_id+ld,:) = [ &
1173 &ind_tot(ind_cur(n_mod_x),1) + 1, & ! total mode number
1174 &kd, & ! lower limit in normal range
1175 &0, & ! initalize upper limit normal range
1176 &modulo(ind_tot(ind_id,4)+ld-1,n_mod_x) + 1] ! index in tables, shifted and wrapped around
1177#if ldebug
1178 if (debug_setup_modes) then
1179 write(*,*) ' FINISHES', ind_cur(1), &
1180 &':', ind_tot(ind_cur(1),:)
1181 write(*,*) ' STARTS', ind_id+ld, ':', &
1182 &ind_tot(ind_id+ld,:)
1183 end if
1184#endif
1185 ind_cur = [ind_cur(2:n_mod_x),ind_id+ld] ! add ind_id+ld at top
1186 else ! mode number has decreased
1187 ind_tot(ind_cur(n_mod_x),3) = kd - 1 ! upper limit normal range
1188 ind_tot(ind_id+ld,:) = [ &
1189 &ind_tot(ind_cur(1),1) - 1, & ! total mode number
1190 &kd, & ! lower limit in normal range
1191 &0, & ! initalize upper limit normal range
1192 &modulo(ind_tot(ind_id,4)-ld,n_mod_x) + 1] ! index in tables, shifted and wrapped around
1193#if ldebug
1194 if (debug_setup_modes) then
1195 write(*,*) ' FINISHES', ind_cur(n_mod_x), &
1196 &':', ind_tot(ind_cur(n_mod_x),:)
1197 write(*,*) ' STARTS', ind_id+ld, ':', &
1198 &ind_tot(ind_id+ld,:)
1199 end if
1200#endif
1201 ind_cur = [ind_id+ld,ind_cur(1:n_mod_x-1)] ! add ind_id+ld at bottom
1202 end if
1203 end do
1204 ind_id = ind_id+abs(delta_ld)
1205#if ldebug
1206 if (debug_setup_modes) write(*,*) ' current indices:', &
1207 &ind_cur
1208#endif
1209 end if
1210 end do
1211
1212 ! close last total modes
1213 do ld = 1,size(ind_cur)
1214 ind_tot(ind_cur(ld),3) = grid%n(3)
1215 end do
1216
1217 ! save in index information
1218 if ((use_pol_flux_f .and. id.eq.2) .or. &
1219 &(.not.use_pol_flux_f .and. id.eq.1)) then
1220 allocate(mds%sec(ind_id,4))
1221 mds%sec = ind_tot(1:ind_id,:)
1222 end if
1223 end do
1224
1225 ! master plots output if requested
1226 if (rank.eq.0 .and. present(plot_name)) then
1227 call writo('Plotting mode numbers')
1228 call lvl_ud(1)
1229
1230 ! total number of modes
1231 n_mod_tot = size(mds%sec,1)
1232
1233 ! set up x values
1234 allocate(x_plot(grid%n(3),n_mod_x))
1235 do ld = 1,n_mod_x
1236 x_plot(:,ld) = grid%r_F
1237 end do
1238 x_plot = x_plot*2*pi/max_flux_f
1239 allocate(plot_titles(1))
1240
1241 ! plot poloidal modes
1242 plot_titles(1) = 'poloidal mode numbers'
1243 plot_name_tot = 'modes_m_'//trim(plot_name)
1244 call print_ex_2d(plot_titles(1:1),plot_name_tot,mds%m*1._dp,&
1245 &x=x_plot,draw=.false.)
1246 call draw_ex(plot_titles(1:1),plot_name_tot,n_mod_x,1,.false.)
1247
1248 ! plot toroidal modes
1249 plot_titles(1) = 'toroidal mode numbers'
1250 plot_name_tot = 'modes_n_'//trim(plot_name)
1251 call print_ex_2d(plot_titles(1:1),plot_name_tot,mds%n*1._dp,&
1252 &x=x_plot,draw=.false.)
1253 call draw_ex(plot_titles(1:1),plot_name_tot,n_mod_x,1,.false.)
1254
1255 ! clean up
1256 deallocate(x_plot)
1257 deallocate(plot_titles)
1258
1259#if ldebug
1260 ! plot secondary limits
1261 allocate(x_plot(n_mod_tot,n_mod_tot+2))
1262 allocate(y_plot(n_mod_tot,n_mod_tot+2))
1263 allocate(plot_titles(n_mod_tot+2))
1264 do ld = 1,n_mod_tot
1265 x_plot(:,ld) = ld
1266 y_plot(:,ld) = mds%sec(ld,2) ! lower boundary
1267 y_plot(n_mod_tot,ld) = mds%sec(ld,3) ! upper boundary
1268 plot_titles(ld) = 'mode '//trim(i2str(mds%sec(ld,1)))
1269 end do
1270 x_plot(:,n_mod_tot+1) = x_plot(1,1:n_mod_tot)
1271 y_plot(:,n_mod_tot+1) = mds%sec(:,1) ! mode number
1272 plot_titles(n_mod_tot+1) = 'mode number'
1273 x_plot(:,n_mod_tot+2) = x_plot(1,1:n_mod_tot)
1274 y_plot(:,n_mod_tot+2) = mds%sec(:,4) ! index
1275 plot_titles(n_mod_tot+2) = 'index'
1276 plot_name_tot = 'modes_sec_'//trim(plot_name)
1277 call print_ex_2d(plot_titles,plot_name_tot,y_plot,x=x_plot,&
1278 &draw=.false.)
1279 call draw_ex(plot_titles,plot_name_tot,n_mod_tot+2,1,.false.)
1280 plot_name_tot = 'modes_sec_flux_'//trim(plot_name)
1281 do ld = 1,n_mod_tot
1282 y_plot(:,ld) = grid%r_F(mds%sec(ld,2))*2*pi/max_flux_f ! lower boundary
1283 y_plot(n_mod_tot,ld) = grid%r_F(mds%sec(ld,3))*2*pi/max_flux_f ! upper boundary
1284 end do
1285 call print_ex_2d(plot_titles(1:n_mod_tot),plot_name_tot,&
1286 &y_plot(:,1:n_mod_tot),x=x_plot(:,1:n_mod_tot),draw=.false.)
1287 call draw_ex(plot_titles(1:n_mod_tot),plot_name_tot,n_mod_tot,1,&
1288 &.false.)
1289 deallocate(x_plot)
1290 deallocate(y_plot)
1291 deallocate(plot_titles)
1292#endif
1293
1294 call lvl_ud(-1)
1295 end if
1296
1297 ! test whether all modes have at least one normal position
1298 if (minval(mds%sec(:,3)-mds%sec(:,2)+1).lt.1) then
1299 ierr = 1
1300 call writo('Mode '//&
1301 &trim(i2str(minloc(mds%sec(:,3)-mds%sec(:,2),1)))//&
1302 &' is not present in any flux surface')
1303 err_msg = 'Aument number of modes or choose finer equilibrium'
1304 chckerr(err_msg)
1305 end if
1306
1307 ! clean up
1308 nullify(nm_x)
1309
1310 call lvl_ud(-1)
1311 call writo('General mode variables set up')
1312 end function setup_modes
1313
1314 !> Checks whether the high-n approximation is valid:
1315 !!
1316 !! This depends on the \c X_style:
1317 !!
1318 !! For X_style 1 (prescribed): every mode should resonate at least somewhere
1319 !! in the whole normal range:
1320 !! - \f$\frac{\left|nq-m\right|}{\left|n\right|} < T\f$ and
1321 !! \f$\frac{\left|nq-m\right|}{\left|m\right|} < T\f$
1322 !! for poloidal flux
1323 !! - \f$\frac{\left|q-\iota m\right|}{\left|m\right|} < T\f$ and
1324 !! \f$\frac{\left|n-\iota m\right|}{\left|n\right|} < T\f$
1325 !! for toroidal flux
1326 !!
1327 !! where \f$T\f$ = \c tol << 1.
1328 !!
1329 !! This condition is determined by the sign of \f$q\f$ (or \f$\iota\f$) and
1330 !! given by:
1331 !! - \f$\max\left(q_\text{min}-T,\frac{q_\text{min}}{1+T}\right) <
1332 !! \frac{m}{n} <
1333 !! \min\left(q_\text{max}+T,\frac{q_\text{max}}{1-T}\right)\f$,
1334 !! \f$q>0\f$
1335 !! - \f$\max\left(q_\text{min}-T,\frac{q_\text{min}}{1-T}\right) <
1336 !! \frac{m}{n} <
1337 !! \min\left(q_\text{max}+T,\frac{q_\text{max}}{1+T}\right)\f$,
1338 !! \f$q<0\f$
1339 !!
1340 !! for poloidal flux.
1341 !!
1342 !! For toroidal flux, \f$q\f$ should be replaced by \f$\iota\f$ and
1343 !! \f$\frac{m}{n}\f$ by \f$\frac{n}{m}\f$.
1344 !!
1345 !! For \c X_style 2 (fast): the resonance has been taken care of, but it
1346 !! remains to be checked whether the number of modes is efficient.
1347 !!
1348 !! \return ierr
1349 integer function check_x_modes(grid_eq,eq) result(ierr)
1350 use num_vars, only: x_style, rank, use_pol_flux_f, n_procs
1351 use mpi_utilities, only: get_ser_var
1352
1353 character(*), parameter :: rout_name = 'check_X_modes'
1354
1355 ! input / output
1356 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
1357 type(eq_1_type), intent(in) :: eq !< flux equilibrium
1358
1359 ! local variables
1360 character(len=max_str_ln) :: err_msg ! error message
1361
1362 ! initialize ierr
1363 ierr = 0
1364
1365 call writo('Checking mode numbers')
1366 call lvl_ud(1)
1367
1368 ! check the modes depending on the X style
1369 select case (x_style)
1370 case (1) ! prescribed
1371 ierr = check_x_modes_1(eq)
1372 chckerr('')
1373 case (2) ! fast
1374 ierr = check_x_modes_2(grid_eq,eq)
1375 chckerr('')
1376 end select
1377
1378 call lvl_ud(-1)
1379 call writo('Mode numbers checked')
1380 contains
1381 ! Version for X style 1: Checks whether there exists a range in which
1382 ! each of the modes resonates, with a certain tolerance tol_norm.
1383 !> \private
1384 integer function check_x_modes_1(eq) result(ierr)
1385 use x_vars, only: prim_x, min_sec_x, n_mod_x
1386 use num_vars, only: tol_norm
1387
1388 character(*), parameter :: rout_name = 'check_X_modes_1'
1389
1390 ! input / output
1391 type(eq_1_type), intent(in) :: eq ! flux equilibrium
1392
1393 ! local variables
1394 integer :: id ! counter
1395 integer :: pmone ! plus or minus one
1396 real(dp) :: min_jq, max_jq ! min. and max. values of q or iota
1397 real(dp) :: lim_lo, lim_hi ! lower and upper limit on n/m (or m/n)
1398 real(dp), allocatable :: temp_var(:) ! temporary variable
1399 character(len=max_str_ln) :: jq_name ! either safety factor or rotational transform
1400 character(len=1) :: mode_name ! either n or m
1401 logical :: all_modes_fine ! all modes are fine
1402
1403 ! initialize ierr
1404 ierr = 0
1405
1406 ! user output
1407 call writo('The tolerance used is '//trim(r2strt(tol_norm))//'...')
1408
1409 ! set min_jq and max_jq in Flux coordinate system
1410 if (use_pol_flux_f) then
1411 min_jq = minval(eq%q_saf_FD(:,0))
1412 max_jq = maxval(eq%q_saf_FD(:,0))
1413 else
1414 min_jq = minval(eq%rot_t_FD(:,0))
1415 max_jq = maxval(eq%rot_t_FD(:,0))
1416 end if
1417
1418 ! combine all processes
1419 ierr = get_ser_var([min_jq],temp_var)
1420 chckerr('')
1421 if (rank.eq.0) min_jq = minval(temp_var)
1422 ierr = get_ser_var([max_jq],temp_var)
1423 chckerr('')
1424 if (rank.eq.0) max_jq = maxval(temp_var)
1425
1426 ! output if master
1427 if (rank.eq.0) then
1428 ! set up jq name and all_modes_fine
1429 if (use_pol_flux_f) then
1430 jq_name = 'safety factor'
1431 mode_name = 'm'
1432 else
1433 jq_name = 'rotational transform'
1434 mode_name = 'n'
1435 end if
1436 all_modes_fine = .true.
1437
1438 ! set up plus minus one, according to the sign of jq
1439 if (min_jq.lt.0 .and. max_jq.lt.0) then
1440 pmone = -1
1441 else if (min_jq.gt.0 .and. max_jq.gt.0) then
1442 pmone = 1
1443 else
1444 err_msg = trim(jq_name)//' cannot change sign'
1445 ierr = 1
1446 chckerr(err_msg)
1447 end if
1448
1449 ! calculate upper and lower limits
1450 lim_lo = max(min_jq-tol_norm,min_jq/(1+pmone*tol_norm))
1451 lim_hi = min(max_jq+tol_norm,max_jq/(1-pmone*tol_norm))
1452
1453 ! for every mode (n,m) check whether m/n is inside the range of
1454 ! q values or n/m inside the range of iota values
1455 do id = 1,n_mod_x
1456 ! check if limits are met
1457 if (use_pol_flux_f) then
1458 if ((min_sec_x+id-1._dp)/prim_x.lt.lim_lo .or. &
1459 &(min_sec_x+id-1._dp)/prim_x.gt.lim_hi) then
1460 call writo('for (n,m) = ('//trim(i2str(prim_x))//&
1461 &','//trim(i2str(min_sec_x+id-1))//&
1462 &'), there is no range in the plasma where the &
1463 &ratio |n q - m| << |n|,|m| is met',alert=.true.)
1464 all_modes_fine = .false.
1465 end if
1466 else
1467 if ((min_sec_x+id-1._dp)/prim_x.lt.lim_lo .or. &
1468 &(min_sec_x+id-1._dp)/prim_x.gt.lim_hi) then
1469 call writo('for (n,m) = ('//&
1470 &trim(i2str(min_sec_x+id-1))//','//&
1471 &trim(i2str(prim_x))//'), there is no range in &
1472 &the plasma where the ratio |n - iota m| << &
1473 &|m|,|n| is met',alert=.true.)
1474 all_modes_fine = .false.
1475 end if
1476 end if
1477 end do
1478
1479 ! output message
1480 if (all_modes_fine) then
1481 call writo('The modes are all within the allowed range &
1482 &of '//trim(i2str(ceiling(prim_x*lim_lo)))//' < '//&
1483 &mode_name//' < '//trim(i2str(floor(prim_x*lim_hi)))//&
1484 &'...')
1485 else
1486 call writo('Not all modes are within the allowed range &
1487 &of '//trim(i2str(ceiling(prim_x*lim_lo)))//' < '//&
1488 &mode_name//' < '//trim(i2str(floor(prim_x*lim_hi)))//&
1489 &'...')
1490 end if
1491 end if
1492 end function check_x_modes_1
1493
1494 ! version for X style 2: Check how efficient the chosen number of modes
1495 ! is.
1496 !> \private
1497 integer function check_x_modes_2(grid_eq,eq) result(ierr)
1498 use x_vars, only: min_n_x, min_m_x, n_mod_x
1499 use grid_utilities, only: trim_grid
1500
1501 character(*), parameter :: rout_name = 'check_X_modes_2'
1502
1503 ! input / output
1504 type(grid_type), intent(in) :: grid_eq ! equilibrium grid
1505 type(eq_1_type), intent(in) :: eq ! flux equilibrium
1506
1507 ! local variables
1508 type(grid_type) :: grid_eq_trim ! trimmed equilibrium grid
1509 integer :: norm_id(2) ! untrimmed normal indices for trimmed grid
1510 integer :: jd, kd, ld ! counters
1511 real(dp), allocatable :: max_frac(:,:) ! maximum fraction
1512 real(dp), allocatable :: max_frac_sum(:) ! sum of maximum fraction for all processes
1513 real(dp), allocatable :: max_frac_max(:) ! maximum maximum fraction for all processes
1514 real(dp), allocatable :: max_frac_temp(:) ! auxilliary variable
1515 real(dp), allocatable :: n(:,:), m(:,:) ! n and m
1516 real(dp), allocatable :: fac_n(:), fac_m(:) ! factors to be multiplied with n m m
1517 character(len=max_str_ln) :: frac_name ! name of fraction
1518#if ldebug
1519 real(dp), allocatable :: x_plot(:,:) ! for plotting
1520 character(len=max_str_ln) :: plot_title ! title for plots
1521 character(len=max_str_ln) :: plot_name ! file name for plots
1522#endif
1523
1524 ! initialize ierr
1525 ierr = 0
1526
1527 ! set up helper variables
1528 allocate(max_frac(grid_eq%loc_n_r,n_mod_x))
1529 allocate(fac_n(grid_eq%loc_n_r))
1530 allocate(fac_m(grid_eq%loc_n_r))
1531 allocate(n(grid_eq%loc_n_r,n_mod_x))
1532 allocate(m(grid_eq%loc_n_r,n_mod_x))
1533 if (use_pol_flux_f) then
1534 do jd = 1,n_mod_x
1535 n(:,jd) = min_n_x(grid_eq%i_min:grid_eq%i_max)
1536 m(:,jd) = min_m_x(grid_eq%i_min:grid_eq%i_max)+jd-1
1537 end do
1538 fac_n = eq%q_saf_FD(:,0)
1539 fac_m = 1._dp
1540 else
1541 do jd = 1,n_mod_x
1542 n(:,jd) = min_n_x(grid_eq%i_min:grid_eq%i_max)+jd-1
1543 m(:,jd) = min_m_x(grid_eq%i_min:grid_eq%i_max)
1544 end do
1545 fac_n = 1._dp
1546 fac_m = eq%rot_t_FD(:,0)
1547 end if
1548
1549 ! loop over all modes
1550 do ld = 1,n_mod_x
1551 do kd = 1,grid_eq%loc_n_r
1552 max_frac(kd,ld) = &
1553 &abs(n(kd,ld)*fac_n(kd)-m(kd,ld)*fac_m(kd))/&
1554 &min(abs(n(kd,ld)),abs(m(kd,ld)))
1555 end do
1556 end do
1557
1558 ! trim grid
1559 ierr = trim_grid(grid_eq,grid_eq_trim,norm_id=norm_id)
1560 chckerr('')
1561
1562 ! gather all fractions
1563 if (rank.eq.0) allocate(max_frac_temp(n_procs))
1564 if (rank.eq.0) allocate(max_frac_sum(n_mod_x))
1565 if (rank.eq.0) allocate(max_frac_max(n_mod_x))
1566 do ld = 1,n_mod_x
1567 ierr = get_ser_var([sum(max_frac(norm_id(1):norm_id(2),ld))],&
1568 &max_frac_temp)
1569 chckerr('')
1570 if (rank.eq.0) max_frac_sum(ld) = sum(max_frac_temp)
1571 ierr = get_ser_var(&
1572 &[maxval(max_frac(norm_id(1):norm_id(2),ld))],max_frac_temp)
1573 chckerr('')
1574 if (rank.eq.0) max_frac_max(ld) = maxval(max_frac_temp)
1575 end do
1576
1577 ! output if master
1578 if (rank.eq.0) then
1579 if (use_pol_flux_f) then
1580 frac_name = '|nq-m|/|n| or |nq-m|/|m|'
1581 else
1582 frac_name = '|n-iota m|/|n| or |n-iota m|/|m|'
1583 end if
1584 call writo('The fraction '//trim(frac_name)//' is')
1585 call lvl_ud(1)
1586 do ld = 1,n_mod_x
1587 call writo('for mode '//trim(i2str(ld))//' maximally '//&
1588 &trim(r2strt(max_frac_max(ld)))//&
1589 &', average '//&
1590 &trim(r2strt(max_frac_sum(ld)/grid_eq%n(3))))
1591 end do
1592 if (n_mod_x.gt.1) call writo('so for all modes maximally '//&
1593 &trim(r2strt(maxval(max_frac_max)))//', average '//&
1594 &trim(r2strt(sum(max_frac_sum)/(grid_eq%n(3)*n_mod_x))))
1595 call lvl_ud(-1)
1596 end if
1597#if ldebug
1598 if (debug_check_x_modes_2) then
1599 call writo('Plotting the fraction for all modes')
1600 allocate(x_plot(grid_eq%n(3),n_mod_x))
1601 do ld = 1,n_mod_x
1602 x_plot(:,ld) = grid_eq%r_F
1603 end do
1604 plot_name = 'TEST_max_frac'
1605 plot_title = 'maximum fraction'
1606 call print_ex_2d([plot_title],plot_name,max_frac,x=x_plot,&
1607 &draw=.false.)
1608 call draw_ex([plot_title],plot_name,n_mod_x,1,.false.)
1609 end if
1610#endif
1611
1612 ! clean up
1613 call grid_eq_trim%dealloc()
1614 end function check_x_modes_2
1615 end function check_x_modes
1616
1617 !> Calculates resonating flux surfaces for the perturbation modes.
1618 !!
1619 !! The output consists of mode number, resonating normal position and the
1620 !! fraction \f$\frac{m}{n}\f$ or \f$\frac{n}{m}\f$. for those modes for
1621 !! which a solution is found that is within the plasma range.
1622 !!
1623 !! It contains three pieces of information:
1624 !! - <tt>(:,1)</tt>: the mode index
1625 !! - <tt>(:,2)</tt>: the radial position in Flux coordinates
1626 !! - <tt>(:,3)</tt>: the fraction \f$\frac{m}{n}\f$ or \f$\frac{n}{m}\f$
1627 !!
1628 !! for every single mode in \c sec of \c mds, which can be tabulated in an
1629 !! arbitrary grid, not necessarily the equilibrium one.
1630 !!
1631 !! Optionally, the total safety factor or rotational transform can be
1632 !! returned to the master.
1633 !!
1634 !! Also, information can be displayed with info.
1635 !!
1636 !! \return ierr
1637 integer function calc_res_surf(mds,grid_eq,eq,res_surf,info,jq) result(ierr)
1638 use x_vars, only: prim_x
1639 use num_vars, only: use_pol_flux_f, norm_disc_prec_eq
1640 use eq_vars, only: max_flux_f, max_flux_f
1641 use num_ops, only: calc_zero_zhang
1642 use grid_utilities, only: trim_grid
1643 use mpi_utilities, only: get_ser_var
1644
1645 character(*), parameter :: rout_name = 'calc_res_surf'
1646
1647 ! input / output
1648 type(modes_type), intent(in) :: mds !< general modes variables
1649 type(grid_type), intent(in) :: grid_eq !< equilibrium grid_eq
1650 type(eq_1_type), intent(in) :: eq !< flux equilibrium
1651 real(dp), intent(inout), allocatable :: res_surf(:,:) !< resonant surface
1652 logical, intent(in), optional :: info !< if info is displayed
1653 real(dp), intent(inout), optional, allocatable :: jq(:) !< either safety factor or rotational transform
1654
1655 ! local variables (also used in child functions)
1656 real(dp) :: nmfrac_fun ! fraction m/n or n/m to determine resonant flux surface
1657
1658 ! local variables (not to be used in child functions)
1659 integer :: norm_id(2) ! untrimmed normal indices for trimmed grid
1660 integer :: ld ! counter
1661 integer :: ld_loc ! local ld
1662 logical :: info_loc ! local version of info
1663 real(dp), allocatable :: res_surf_loc(:,:) ! local copy of res_surf
1664 real(dp), allocatable :: jq_tot(:) ! saf. fac. or rot. transf. and derivs. in Flux coords.
1665 real(dp), allocatable :: jq_loc(:) ! local version of jq
1666 real(dp) :: norm_factor ! normalization factor to make normal coordinate 0..1
1667 type(grid_type) :: grid_eq_trim ! trimmed version of equilibrium grid
1668 integer :: n_loc, m_loc ! local n and m
1669 character(len=max_str_ln) :: err_msg ! error message
1670
1671 ! initialize ierr
1672 ierr = 0
1673
1674 ! set local info
1675 info_loc = .false.
1676 if (present(info)) info_loc = info
1677
1678 ! get trimmed grid
1679 ierr = trim_grid(grid_eq,grid_eq_trim,norm_id)
1680 chckerr('')
1681
1682 ! get serial version of safety factor or rot. transform
1683 allocate(jq_tot(grid_eq_trim%n(3)))
1684 if (use_pol_flux_f) then
1685 if (grid_eq%divided) then
1686 ierr = get_ser_var(eq%q_saf_FD(norm_id(1):norm_id(2),0),&
1687 &jq_loc,scatter=.true.)
1688 chckerr('')
1689 jq_tot = jq_loc
1690 else
1691 jq_tot = eq%q_saf_FD(:,0)
1692 end if
1693 else
1694 if (grid_eq%divided) then
1695 ierr = get_ser_var(eq%rot_t_FD(norm_id(1):norm_id(2),0),&
1696 &jq_loc,scatter=.true.)
1697 chckerr('')
1698 jq_tot = jq_loc
1699 else
1700 jq_tot = eq%rot_t_FD(:,0)
1701 end if
1702 end if
1703 if (present(jq)) then
1704 allocate(jq(grid_eq_trim%n(3)))
1705 jq = jq_tot
1706 end if
1707
1708 ! initialize local res_surf
1709 allocate(res_surf_loc(1:size(mds%sec,1),3))
1710
1711 ! calculate normalization factor max_flux / 2pi
1712 norm_factor = max_flux_f/(2*pi)
1713
1714 ! loop over all modes (and shift the index in m_loc and n_loc by 1)
1715 ld_loc = 1
1716 do ld = 1,size(mds%sec,1)
1717 ! Find place where q = m/n or iota = n/m in Flux coordinates by
1718 ! solving q-m/n = 0 or iota-n/m=0, using the functin jq_fun.
1719
1720 ! set up local m, n and nmfrac for function
1721 if (use_pol_flux_f) then
1722 n_loc = prim_x
1723 m_loc = mds%sec(ld,1)
1724 nmfrac_fun = 1.0_dp*m_loc/n_loc
1725 else
1726 n_loc = mds%sec(ld,1)
1727 m_loc = prim_x
1728 nmfrac_fun = 1.0_dp*n_loc/m_loc
1729 end if
1730
1731 ! calculate zero using Zhang
1732 if (use_pol_flux_f) then
1733 res_surf_loc(ld_loc,1) = m_loc
1734 else
1735 res_surf_loc(ld_loc,1) = n_loc
1736 end if
1737 res_surf_loc(ld_loc,3) = nmfrac_fun
1738 err_msg = calc_zero_zhang(res_surf_loc(ld_loc,2),jq_fun,&
1739 &[minval(grid_eq_trim%r_F),maxval(grid_eq_trim%r_F)])
1740
1741 ! intercept error
1742 if (err_msg.ne.'') then
1743 if (info_loc) then
1744 call writo('Mode (n,m) = ('//trim(i2str(n_loc))//','//&
1745 &trim(i2str(m_loc))//') is not found')
1746 call lvl_ud(1)
1747 call writo(trim(err_msg))
1748 call lvl_ud(-1)
1749 end if
1750 else if (res_surf_loc(ld_loc,2).lt.minval(grid_eq_trim%r_F) .or. &
1751 &res_surf_loc(ld_loc,2).gt.maxval(grid_eq_trim%r_F)) then
1752 if (info_loc) call writo('Mode (n,m) = ('//&
1753 &trim(i2str(n_loc))//','//trim(i2str(m_loc))//&
1754 &') does not resonate in plasma')
1755 else
1756 if (info_loc) call writo('Mode (n,m) = ('//&
1757 &trim(i2str(n_loc))//','//trim(i2str(m_loc))//&
1758 &') resonates in plasma at normalized flux surface '//&
1759 &trim(r2str(res_surf_loc(ld_loc,2)/norm_factor)))
1760 ld_loc = ld_loc + 1 ! advance ld_loc
1761 end if
1762 end do
1763
1764 ! set res_surf from local copy
1765 allocate(res_surf(ld_loc-1,3))
1766 res_surf = res_surf_loc(1:ld_loc-1,:)
1767
1768 ! clean up
1769 call grid_eq_trim%dealloc()
1770 contains
1771 ! Returns q-m/n or iota-n/m in Flux coordinates, used to solve for q =
1772 ! m/n or iota = n/m.
1773 !> \private
1774 real(dp) function jq_fun(pt) result(res)
1775 use num_utilities, only: spline
1776
1777 ! input / output
1778 real(dp), intent(in) :: pt ! normal position at which to evaluate
1779
1780 ! local variables
1781 integer :: i_min, i_max ! index of minimum and maximum value of x
1782 real(dp) :: x_min, x_max ! minimum and maximum value of x
1783 real(dp) :: res_loc(1) ! local copy of res
1784
1785 ! initialize res
1786 res = 0
1787
1788 ! set up min. and max index and value
1789 x_min = minval(grid_eq_trim%r_F)
1790 x_max = maxval(grid_eq_trim%r_F)
1791 i_min = minloc(grid_eq_trim%r_F,1)
1792 i_max = maxloc(grid_eq_trim%r_F,1)
1793
1794 ! check whether to interpolate or extrapolate
1795 if (pt.ge.x_min .and. pt.le.x_max) then ! point requested within range
1796 ! interpolate
1797 ierr = spline(grid_eq_trim%r_F,jq_tot-nmfrac_fun,[pt],res_loc,&
1798 &ord=norm_disc_prec_eq)
1799 chckerr('')
1800 ! copy
1801 res = res_loc(1)
1802 end if
1803 end function jq_fun
1804 end function calc_res_surf
1805
1806 !> plot \f$q\f$-profile or \f$\iota\f$-profile in flux coordinates with
1807 !! \f$nq-m = 0\f$ or \f$n-\iota m = 0\f$ indicated if requested.
1808 !!
1809 !! The plot will be done in the grid in which \c mds is tabulated, which is
1810 !! not necessarily the equilibrium one.
1811 !!
1812 !! \return ierr
1813 integer function resonance_plot(mds,grid_eq,eq) result(ierr)
1814 use num_vars, only: use_pol_flux_f, no_plots, n_theta_plot, &
1815 &n_zeta_plot, rank, eq_style, min_theta_plot, max_theta_plot, &
1816 &min_zeta_plot, max_zeta_plot
1817 use eq_vars, only: max_flux_f
1818 use grid_utilities, only: calc_xyz_grid, calc_eqd_grid, coord_f2e, &
1819 &trim_grid
1820 use x_vars, only: prim_x
1821 use vmec_utilities, only: calc_trigon_factors
1822
1823 character(*), parameter :: rout_name = 'resonance_plot'
1824
1825 ! input / output
1826 type(modes_type), intent(in) :: mds !< general modes variables
1827 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
1828 type(eq_1_type), intent(in) :: eq !< flux equilibrium
1829
1830 ! local variables (not to be used in child functions)
1831 integer :: ld ! counter
1832 integer :: n_mod_loc ! local n_mod
1833 real(dp), allocatable :: res_surf(:,:) ! resonant surfaces
1834 real(dp), allocatable :: x_plot_loc(:,:) ! for plotting
1835 real(dp), allocatable :: y_plot_loc(:,:) ! for plotting
1836 character(len=max_str_ln) :: plot_title, file_name ! name of plot, of file
1837 real(dp), allocatable :: jq(:) ! saf. fac. or rot. transf. in Flux coords.
1838 integer :: n_r ! total number of normal points
1839 integer :: plot_dim(4) ! plot dimensions (total = local because only masters)
1840 type(grid_type) :: grid_trim ! trimmed version of grid
1841 type(grid_type) :: grid_plot ! grid for plotting
1842 real(dp), allocatable :: r_plot_e(:) ! normal E coordinates of plot (needed to calculate X, Y and Z)
1843 real(dp), allocatable :: theta_plot(:,:,:), zeta_plot(:,:,:) ! pol. and tor. angle of plot
1844 real(dp), allocatable :: x_plot(:,:,:,:), y_plot(:,:,:,:), &
1845 &Z_plot(:,:,:,:) ! X, Y and Z of plot of all surfaces
1846 real(dp), allocatable :: vars(:,:,:,:) ! variable to plot
1847 character(len=max_str_ln), allocatable :: plot_titles(:) ! name of plots
1848
1849 ! initialize ierr
1850 ierr = 0
1851
1852 ! bypass plots if no_plots
1853 if (no_plots) return
1854
1855 ! get trimmed grid
1856 ierr = trim_grid(grid_eq,grid_trim)
1857 chckerr('')
1858
1859 ! initialize variables
1860 n_r = grid_trim%n(3)
1861
1862 ! print user output
1863 if (use_pol_flux_f) then
1864 call writo('Plotting safety factor q and resonant surfaces &
1865 &q = m/n...')
1866 plot_title = 'safety factor q'
1867 file_name = 'q_saf'
1868 else
1869 call writo('Plotting rotational transform iota and resonant &
1870 &surfaces iota = n/m...')
1871 plot_title = 'rotational transform iota'
1872 file_name = 'rot_t'
1873 end if
1874
1875 call lvl_ud(1)
1876
1877 call writo('Calculate resonant surfaces')
1878 call lvl_ud(1)
1879
1880 ! find resonating surfaces
1881 ierr = calc_res_surf(mds,grid_eq,eq,res_surf,info=.true.,jq=jq)
1882 chckerr('')
1883
1884 call lvl_ud(-1)
1885
1886 ! only master
1887 if (rank.eq.0) then
1888 ! set local n_mod
1889 n_mod_loc = size(res_surf,1)
1890
1891 ! set up plot titles
1892 allocate(plot_titles(n_mod_loc))
1893 if (use_pol_flux_f) then ! n is fixed and m = m/n n
1894 do ld = 1,n_mod_loc
1895 plot_titles(ld) = 'resonance for m,n = '//&
1896 &trim(i2str(nint(res_surf(ld,3)*prim_x)))//','//&
1897 &trim(i2str(prim_x))
1898 end do
1899 else ! m is fixed and n = n/m m
1900 do ld = 1,n_mod_loc
1901 plot_titles(ld) = 'resonance for m,n = '//&
1902 &trim(i2str(prim_x))//','//&
1903 &trim(i2str(nint(res_surf(ld,3)*prim_x)))
1904 end do
1905 end if
1906
1907 ! initialize x_plot_loc and y_plot_loc
1908 allocate(x_plot_loc(n_r,n_mod_loc+1)); x_plot_loc = 0
1909 allocate(y_plot_loc(n_r,n_mod_loc+1)); y_plot_loc = 0
1910
1911 ! set x_plot_loc and y_plot_loc for first column
1912 x_plot_loc(:,1) = grid_trim%r_F
1913 y_plot_loc(:,1) = jq(:)
1914
1915 ! set x_plot_loc and y_plot_loc for other columns
1916 do ld = 1,n_mod_loc
1917 x_plot_loc(:,ld+1) = res_surf(ld,2)
1918 y_plot_loc(n_r,ld+1) = res_surf(ld,3)
1919 end do
1920
1921 ! only if resonant surfaces
1922 if (size(res_surf,1).gt.0) then
1923 ! user message
1924 call writo('Plot resonant flux surfaces using HDF5')
1925
1926 call lvl_ud(1)
1927
1928 ! set up pol. and tor. angle for plot
1929 allocate(theta_plot(n_theta_plot,n_zeta_plot,1))
1930 allocate(zeta_plot(n_theta_plot,n_zeta_plot,1))
1931 ierr = calc_eqd_grid(theta_plot,min_theta_plot*pi,&
1932 &max_theta_plot*pi,1)
1933 chckerr('')
1934 ierr = calc_eqd_grid(zeta_plot,min_zeta_plot*pi,&
1935 &max_zeta_plot*pi,2)
1936 chckerr('')
1937
1938 ! set up vars
1939 allocate(vars(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1940 do ld = 1,n_mod_loc
1941 vars(:,:,:,ld) = y_plot_loc(n_r,ld+1)
1942 end do
1943
1944 ! set dimensions for single flux surface
1945 plot_dim = [n_theta_plot,n_zeta_plot,1,n_mod_loc]
1946
1947 ! calculate normal vars in Equilibrium coords.
1948 allocate(r_plot_e(n_mod_loc))
1949 ierr = coord_f2e(grid_eq,x_plot_loc(n_r,2:n_mod_loc+1),&
1950 &r_plot_e,r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
1951 chckerr('')
1952
1953 ! create plot grid for single flux surface
1954 ierr = grid_plot%init(plot_dim(1:3))
1955 chckerr('')
1956 grid_plot%theta_E = theta_plot
1957 grid_plot%zeta_E = zeta_plot
1958
1959 ! if VMEC, calculate trigonometric factors of plot grid
1960 if (eq_style.eq.1) then
1961 ierr = calc_trigon_factors(grid_plot%theta_E,&
1962 &grid_plot%zeta_E,grid_plot%trigon_factors)
1963 chckerr('')
1964 end if
1965
1966 ! set up plot X, Y and Z
1967 allocate(x_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1968 allocate(y_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1969 allocate(z_plot(n_theta_plot,n_zeta_plot,1,n_mod_loc))
1970
1971 ! loop over all resonant surfaces to calculate X, Y and Z values
1972 do ld = 1,n_mod_loc
1973 ! set loc_r_E of plot grid
1974 grid_plot%loc_r_E = r_plot_e(ld)
1975
1976 ! calculate X, Y and Z
1977 ierr = calc_xyz_grid(grid_eq,grid_plot,x_plot(:,:,:,ld),&
1978 &y_plot(:,:,:,ld),z_plot(:,:,:,ld))
1979 chckerr('')
1980 end do
1981
1982 ! plot using HDF5
1983 call plot_hdf5(plot_titles,file_name,vars,x=x_plot,y=y_plot,&
1984 &z=z_plot,col=1,descr='resonant surfaces')
1985
1986 ! deallocate local variables
1987 deallocate(vars)
1988 deallocate(theta_plot,zeta_plot,r_plot_e)
1989 deallocate(x_plot,y_plot,z_plot)
1990 call grid_plot%dealloc()
1991
1992 call lvl_ud(-1)
1993 end if
1994
1995 call writo('Plot safety factor with resonant flux surfaces using &
1996 &external program')
1997
1998 call lvl_ud(1)
1999
2000 ! rescale x_plot_loc by max_flux_F/2*pi
2001 x_plot_loc = x_plot_loc*2*pi/max_flux_f
2002
2003 ! plot using external program
2004 call print_ex_2d([plot_title,plot_titles],file_name,y_plot_loc,&
2005 &x=x_plot_loc,draw=.false.)
2006 call draw_ex([plot_title,plot_titles],file_name,n_mod_loc+1,1,&
2007 &.false.)
2008
2009 call lvl_ud(-1)
2010
2011 ! only if resonant surfaces
2012 if (size(res_surf,1).gt.0) then
2013 call writo('Plot 2D map of resonances using HDF5')
2014
2015 call lvl_ud(1)
2016
2017 ! initialize plot dimensions
2018 plot_dim = [1,n_mod_loc,1,0]
2019
2020 ! set up X, Y and Z
2021 allocate(x_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2022 allocate(y_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2023 allocate(z_plot(plot_dim(1),plot_dim(2),plot_dim(3),1))
2024
2025 do ld = 1,n_mod_loc
2026 x_plot(1,ld,1,1) = x_plot_loc(1,ld+1)
2027 y_plot(1,ld,1,1) = res_surf(ld,1)
2028 z_plot(1,ld,1,1) = 1._dp
2029 end do
2030
2031 ! plot 2D (to be used with plots of harmonics in sol_ops)
2032 call plot_hdf5('resonating surfaces','res_surf',&
2033 &z_plot(:,:,:,1),x=x_plot(:,:,:,1),y=y_plot(:,:,:,1))
2034
2035 ! deallocate local variables
2036 deallocate(x_plot,y_plot,z_plot)
2037
2038 call lvl_ud(-1)
2039 end if
2040 end if
2041
2042 ! clean up
2043 call grid_trim%dealloc()
2044
2045 call lvl_ud(-1)
2046 end function resonance_plot
2047
2048 !> Calculate \f$U_m^0\f$, \f$U_m^1\f$ or \f$U_n^0\f$, \f$U_n^1\f$.
2049 !!
2050 !! This is done at eq \c loc_n_r values of the normal coordinate, \c n_par
2051 !! values of the parallel coordinate and \c size_X values of the poloidal
2052 !! mode number, or of the toroidal mode number, as well as the poloidal
2053 !! derivatives
2054 !!
2055 !! Use is made of variables \c Ti that are set up in the equilibrium grid
2056 !! and are afterwards converted to \c Ti_X in the perturbation grid, which
2057 !! needs to have the same angular coordinates as the equilibrium grid:
2058 !! \f[\begin{aligned}
2059 !! T_1 &= \frac{B_\alpha}{B_\theta} \\
2060 !! T_2 &= \Theta^\alpha + q' \theta \\
2061 !! T_3 &= \frac{B_\alpha}{B_\theta} q' +
2062 !! \frac{\mathcal{J}^2}{B_\theta} \mu_0 p' \\
2063 !! T_4 &= \frac{B_\alpha}{B_\theta} q' \theta -
2064 !! \frac{B_\psi}{B_\theta} \\
2065 !! T_5 &= \frac{B_\alpha}{B_\theta}
2066 !! \frac{\partial \Theta^\theta}{\partial \theta} -
2067 !! \frac{\partial\Theta^\theta}{\partial \alpha} \\
2068 !! T_6 &= \frac{B_\alpha}{B_\theta} \Theta^\theta \ ,
2069 !! \end{aligned}\f]
2070 !! with \f$T_i\f$ = \c Ti.
2071 !!
2072 !! which is valid for poloidal Flux coordinates.
2073 !!
2074 !! For toroidal Flux coordinates, \f$q'\f$ has to be replaced by
2075 !! \f$-\iota'\f$ and \f$\theta\f$ by \f$\zeta\f$.
2076 !!
2077 !! The interpolated \c Ti_X are then used to calculate \f$U\f$:
2078 !! \f[\begin{aligned}
2079 !! U_0 =& -T_2 &\text{(order 1)}\\
2080 !! &+ \frac{i}{n}\left(T_3 + i(nq-m) T_4\right) &\text{(order 2)}\\
2081 !! &+ \left(\frac{i}{n}\right)^2 i(nq-m)
2082 !! \left(-T_5 - (nq-m) T_6\right) &\text{(order 3)}\\
2083 !! U_1 =& \frac{i}{n} &\text{(order 1)}\\
2084 !! &+ \left(\frac{i}{n}\right)^2 i(nq-m)(-T_1) &\text{(order 2)} .
2085 !! \end{aligned}\f]
2086 !!
2087 !! This is valid for poloidal Flux coordinates and where \f$n\f$ is to be
2088 !! replaced by \f$m\f$ and \f$(nq-m)\f$ by \f$(n-\iota m)\f$ for toroidal
2089 !! Flux coordinates.
2090 !!
2091 !! For VMEC, these factors are also derived in the parallel coordinate.
2092 !!
2093 !! \see See \cite weyens2014theory.
2094 !!
2095 !! \return ierr
2096 integer function calc_u(grid_eq,grid_X,eq_1,eq_2,X) result(ierr)
2097 use num_vars, only: use_pol_flux_f, eq_style, u_style, &
2098 &norm_disc_prec_x, x_grid_style
2099 use num_utilities, only: c, spline
2100 use input_utilities, only: get_log, pause_prog
2101 use eq_vars, only: vac_perm
2102#if ldebug
2103 use num_vars, only: ltest
2104#endif
2105
2106 character(*), parameter :: rout_name = 'calc_U'
2107
2108 ! input / output
2109 type(grid_type), intent(in) :: grid_eq !< equilibrium grid variables
2110 type(grid_type), intent(in) :: grid_x !< perturbation grid variables
2111 type(eq_1_type), intent(in), target :: eq_1 !< flux equilibrium
2112 type(eq_2_type), intent(in), target :: eq_2 !< metric equilibrium
2113 type(x_1_type), intent(inout) :: x !< vectorial perturbation variables
2114
2115 ! local variables
2116 integer :: id, jd, kd, ld ! counters
2117 integer :: t_size ! 2 for VMEC and 1 for HELENA
2118
2119 ! Jacobian
2120 real(dp), pointer :: j(:,:,:) ! Jacobian
2121 real(dp), pointer :: d3j(:,:,:) ! D_theta Jacobian
2122 ! lower metric factors
2123 real(dp), pointer :: g13(:,:,:) ! g_alpha,theta
2124 real(dp), pointer :: d3g13(:,:,:) ! D_theta g_alpha,theta
2125 real(dp), pointer :: g23(:,:,:) ! g_psi,theta
2126 real(dp), pointer :: d3g23(:,:,:) ! D_theta g_psi,theta
2127 real(dp), pointer :: g33(:,:,:) ! g_theta,theta
2128 real(dp), pointer :: d3g33(:,:,:) ! D_theta g_theta,theta
2129 ! upper metric factors
2130 real(dp), pointer :: h12(:,:,:) ! h^alpha,psi
2131 real(dp), pointer :: d3h12(:,:,:) ! D_theta h^alpha,psi
2132 real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2133 real(dp), pointer :: d1h22(:,:,:) ! D_alpha h^psi,psi
2134 real(dp), pointer :: d3h22(:,:,:) ! D_theta h^psi,psi
2135 real(dp), pointer :: d13h22(:,:,:) ! D_alpha,theta h^psi,psi
2136 real(dp), pointer :: d33h22(:,:,:) ! D_theta,theta h^psi,psi
2137 real(dp), pointer :: h23(:,:,:) ! h^psi,theta
2138 real(dp), pointer :: d1h23(:,:,:) ! D_alpha h^psi,theta
2139 real(dp), pointer :: d3h23(:,:,:) ! D_theta h^psi,theta
2140 real(dp), pointer :: d13h23(:,:,:) ! D_alpha,theta h^psi,theta
2141 real(dp), pointer :: d33h23(:,:,:) ! D_theta,theta h^psi,theta
2142 ! helper variables
2143 real(dp), pointer :: ang_par_f(:,:,:) ! equilibrium parallel angle in flux coordinates
2144 real(dp), pointer :: ang_geo_f(:,:,:) ! equilibrium geodesical angle in flux coordinates
2145 real(dp), pointer :: q_saf(:), rot_t(:) ! safety factor and rotational transform in X grid
2146 real(dp), allocatable :: djq(:) ! either q' (pol. flux) or -iota' (tor. flux)
2147 real(dp), allocatable :: theta_3(:,:,:), d1theta_3(:,:,:), &
2148 &D3Theta_3(:,:,:) ! Theta^theta and derivatives
2149 real(dp), allocatable :: d13theta_3(:,:,:), d33theta_3(:,:,:) ! Theta^theta and derivatives
2150 real(dp), allocatable :: n_frac(:,:) ! nq-m (pol. flux) or n-iotam (tor. flux) for all modes
2151 real(dp), allocatable :: nm(:,:) ! either n*A_0 (pol. flux) or m (tor.flux)
2152 complex(dp), allocatable :: u_corr(:,:,:,:) ! correction to U_0 and U_1 for a certain (n,m)
2153 complex(dp), allocatable :: d3u_corr(:,:,:,:) ! D_theta U_corr
2154 ! U factors
2155 real(dp), allocatable, target :: t1(:,:,:,:) ! B_alpha/B_theta and par. deriv.
2156 real(dp), allocatable, target :: t2(:,:,:,:) ! Theta^alpha + q' theta and par. deriv.
2157 real(dp), allocatable, target :: t3(:,:,:,:) ! B_alpha/B_theta q' + J^2/B_theta mu_0 p' and par. deriv.
2158 real(dp), allocatable, target :: t4(:,:,:,:) ! B_alpha/B_theta q' theta - B_psi/B_theta and par. deriv.
2159 real(dp), allocatable, target :: t5(:,:,:,:) ! B_alpha/B_theta D3Theta^theta - D1Theta^theta and par. deriv.
2160 real(dp), allocatable, target :: t6(:,:,:,:) ! B_alpha/B_theta Theta^theta and par. deriv.
2161 real(dp), pointer :: t1_x(:,:,:,:) ! T1 in X grid and par. deriv.
2162 real(dp), pointer :: t2_x(:,:,:,:) ! T2 in X grid and par. deriv.
2163 real(dp), pointer :: t3_x(:,:,:,:) ! T3 in X grid and par. deriv.
2164 real(dp), pointer :: t4_x(:,:,:,:) ! T4 in X grid and par. deriv.
2165 real(dp), pointer :: t5_x(:,:,:,:) ! T5 in X grid and par. deriv.
2166 real(dp), pointer :: t6_x(:,:,:,:) ! T6 in X grid and par. deriv.
2167
2168 ! initialize ierr
2169 ierr = 0
2170
2171 ! allocate variables
2172 ! helper variables
2173 allocate(nm(grid_x%loc_n_r,x%n_mod))
2174 allocate(q_saf(grid_x%loc_n_r))
2175 allocate(rot_t(grid_x%loc_n_r))
2176 allocate(djq(grid_eq%loc_n_r))
2177 allocate(theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2178 allocate(d1theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2179 allocate(d3theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2180 allocate(d13theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2181 allocate(d33theta_3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2182 allocate(n_frac(grid_x%loc_n_r,x%n_mod))
2183 allocate(u_corr(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,2))
2184 allocate(d3u_corr(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,2))
2185 ! U factors
2186 select case (eq_style)
2187 case (1) ! VMEC
2188 t_size = 2
2189 case (2) ! HELENA
2190 t_size = 1
2191 end select
2192 allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2193 allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2194 allocate(t3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2195 allocate(t4(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2196 allocate(t5(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2197 allocate(t6(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r,t_size))
2198 select case (x_grid_style)
2199 case (1) ! equilibrium
2200 ! do nothing
2201 case (2,3) ! solution or enriched
2202 allocate(q_saf(grid_x%loc_n_r))
2203 allocate(rot_t(grid_x%loc_n_r))
2204 allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2205 allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2206 allocate(t3_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2207 allocate(t4_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2208 allocate(t5_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2209 allocate(t6_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r,t_size))
2210 end select
2211
2212 ! set pointers
2213 if (use_pol_flux_f) then
2214 ang_par_f => grid_eq%theta_F
2215 ang_geo_f => grid_eq%zeta_F
2216 else
2217 ang_par_f => grid_eq%zeta_F
2218 ang_geo_f => grid_eq%theta_F
2219 end if
2220 j => eq_2%jac_FD(:,:,:,0,0,0)
2221 d3j => eq_2%jac_FD(:,:,:,0,0,1)
2222 g13 => eq_2%g_FD(:,:,:,c([1,3],.true.),0,0,0)
2223 d3g13 => eq_2%g_FD(:,:,:,c([1,3],.true.),0,0,1)
2224 g23 => eq_2%g_FD(:,:,:,c([2,3],.true.),0,0,0)
2225 d3g23 => eq_2%g_FD(:,:,:,c([2,3],.true.),0,0,1)
2226 g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2227 d3g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,1)
2228 h12 => eq_2%h_FD(:,:,:,c([1,2],.true.),0,0,0)
2229 d3h12 => eq_2%h_FD(:,:,:,c([1,2],.true.),0,0,1)
2230 h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2231 d1h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),1,0,0)
2232 d3h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,1)
2233 d13h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),1,0,1)
2234 d33h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,2)
2235 h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,0)
2236 d1h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),1,0,0)
2237 d3h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,1)
2238 d13h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),1,0,1)
2239 d33h23 => eq_2%h_FD(:,:,:,c([2,3],.true.),0,0,2)
2240
2241 ! set up helper variables in eq grid
2242 if (use_pol_flux_f) then
2243 nm = x%n
2244 djq = eq_1%q_saf_FD(:,1)
2245 else
2246 djq = -eq_1%rot_t_FD(:,1)
2247 nm = x%m
2248 end if
2249 theta_3 = h23/h22
2250 select case (eq_style)
2251 case (1) ! VMEC
2252 d1theta_3 = d1h23/h22 - h23*d1h22/(h22**2)
2253 d3theta_3 = d3h23/h22 - h23*d3h22/(h22**2)
2254 d13theta_3 = d13h23/h22 - (d3h23*d1h22+d1h23*d3h22)/(h22**2) &
2255 &- h23*d13h22/(h22**2) + 2*h23*d3h22*d1h22/(h22**3)
2256 d33theta_3 = d33h23/h22 - 2*d3h23*d3h22/(h22**2) &
2257 &- h23*d33h22/(h22**2) + 2*h23*d3h22**2/(h22**3)
2258 case (2) ! HELENA
2259 ! geodesic derivative
2260 if (grid_eq%n(2).gt.norm_disc_prec_x+3) then ! need enough terms
2261 do kd = 1,grid_eq%loc_n_r
2262 do id = 1,grid_eq%n(1)
2263 ierr = spline(ang_geo_f(id,:,kd),theta_3(id,:,kd),&
2264 &ang_geo_f(id,:,kd),d1theta_3(id,:,kd),&
2265 &ord=norm_disc_prec_x,deriv=1)
2266 chckerr('')
2267 end do
2268 end do
2269 else
2270 d1theta_3 = 0._dp
2271 end if
2272 ! parallel derivative
2273 if (grid_eq%n(1).gt.norm_disc_prec_x+3) then ! need enough terms
2274 do kd = 1,grid_eq%loc_n_r
2275 do jd = 1,grid_eq%n(2)
2276 ierr = spline(ang_par_f(:,jd,kd),theta_3(:,jd,kd),&
2277 &ang_par_f(:,jd,kd),d3theta_3(:,jd,kd),&
2278 &ord=norm_disc_prec_x,deriv=1)
2279 chckerr('')
2280 end do
2281 end do
2282 else
2283 d3theta_3 = 0._dp
2284 end if
2285 end select
2286
2287 ! set up U factors in eq grid
2288 t1(:,:,:,1) = g13/g33
2289 do kd = 1,grid_eq%loc_n_r
2290 t2(:,:,kd,1) = h12(:,:,kd)/h22(:,:,kd) + djq(kd)*ang_par_f(:,:,kd)
2291 t3(:,:,kd,1) = t1(:,:,kd,1)*djq(kd) + &
2292 &j(:,:,kd)**2*vac_perm*eq_1%pres_FD(kd,1)/g33(:,:,kd)
2293 t4(:,:,kd,1) = t1(:,:,kd,1)*djq(kd)*ang_par_f(:,:,kd) - &
2294 &g23(:,:,kd)/g33(:,:,kd)
2295 end do
2296 t5(:,:,:,1) = t1(:,:,:,1)*d3theta_3 - d1theta_3
2297 t6(:,:,:,1) = t1(:,:,:,1)*theta_3
2298
2299 ! set up parallel derivatives of U in eq grid if VMEC is used
2300 if (eq_style.eq.1) then
2301 t1(:,:,:,2) = d3g13/g33 - g13*d3g33/g33**2
2302 do kd = 1,grid_eq%loc_n_r
2303 t2(:,:,kd,2) = d3h12(:,:,kd)/h22(:,:,kd) - &
2304 &h12(:,:,kd)*d3h22(:,:,kd)/h22(:,:,kd)**2 + djq(kd)
2305 t3(:,:,kd,2) = t1(:,:,kd,2)*djq(kd) + &
2306 &j(:,:,kd)*vac_perm*eq_1%pres_FD(kd,1)/g33(:,:,kd) * &
2307 &(2*d3j(:,:,kd)-d3g33(:,:,kd)*j(:,:,kd)/g33(:,:,kd))
2308 t4(:,:,kd,2) = (t1(:,:,kd,2)*ang_par_f(:,:,kd)+t1(:,:,kd,1))*&
2309 &djq(kd) - d3g23(:,:,kd)/g33(:,:,kd) + &
2310 &g23(:,:,kd)*d3g33(:,:,kd)/g33(:,:,kd)**2
2311 end do
2312 t5(:,:,:,2) = t1(:,:,:,2)*d3theta_3 + t1(:,:,:,1)*d33theta_3 - &
2313 &d13theta_3
2314 t6(:,:,:,2) = t1(:,:,:,2)*theta_3 + t1(:,:,:,1)*d3theta_3
2315 end if
2316
2317 select case (x_grid_style)
2318 case (1) ! equilibrium
2319 ! point
2320 q_saf => eq_1%q_saf_FD(:,0)
2321 rot_t => eq_1%rot_t_FD(:,0)
2322 t1_x => t1
2323 t2_x => t2
2324 t3_x => t3
2325 t4_x => t4
2326 t5_x => t5
2327 t6_x => t6
2328 case (2,3) ! solution or enriched
2329 ! interpolate
2330 ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
2331 &grid_x%loc_r_F,q_saf,ord=norm_disc_prec_x)
2332 chckerr('')
2333 ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
2334 &grid_x%loc_r_F,rot_t,ord=norm_disc_prec_x)
2335 chckerr('')
2336 do ld = 1,t_size
2337 do jd = 1,grid_x%n(2)
2338 do id = 1,grid_x%n(1)
2339 ierr = spline(grid_eq%loc_r_F,t1(id,jd,:,ld),&
2340 &grid_x%loc_r_F,t1_x(id,jd,:,ld),&
2341 &ord=norm_disc_prec_x)
2342 chckerr('')
2343 ierr = spline(grid_eq%loc_r_F,t2(id,jd,:,ld),&
2344 &grid_x%loc_r_F,t2_x(id,jd,:,ld),&
2345 &ord=norm_disc_prec_x)
2346 chckerr('')
2347 ierr = spline(grid_eq%loc_r_F,t3(id,jd,:,ld),&
2348 &grid_x%loc_r_F,t3_x(id,jd,:,ld),&
2349 &ord=norm_disc_prec_x)
2350 chckerr('')
2351 ierr = spline(grid_eq%loc_r_F,t4(id,jd,:,ld),&
2352 &grid_x%loc_r_F,t4_x(id,jd,:,ld),&
2353 &ord=norm_disc_prec_x)
2354 chckerr('')
2355 ierr = spline(grid_eq%loc_r_F,t5(id,jd,:,ld),&
2356 &grid_x%loc_r_F,t5_x(id,jd,:,ld),&
2357 &ord=norm_disc_prec_x)
2358 chckerr('')
2359 ierr = spline(grid_eq%loc_r_F,t6(id,jd,:,ld),&
2360 &grid_x%loc_r_F,t6_x(id,jd,:,ld),&
2361 &ord=norm_disc_prec_x)
2362 chckerr('')
2363 end do
2364 end do
2365 end do
2366 end select
2367
2368 ! set up n_frac
2369 do kd = 1,grid_x%loc_n_r
2370 if (use_pol_flux_f) then
2371 n_frac(kd,:) = x%n(kd,:)*q_saf(kd) - x%m(kd,:)
2372 else
2373 n_frac(kd,:) = -x%m(kd,:)*rot_t(kd) + x%n(kd,:)
2374 end if
2375 end do
2376
2377 ! calculate U_0 and U_1 for all modes
2378 do ld = 1,x%n_mod
2379 ! calculate order 1 of U_0 and U_1
2380 if (u_style.ge.1) then
2381 x%U_0(:,:,:,ld) = -t2_x(:,:,:,1)
2382 do kd = 1,grid_x%loc_n_r
2383 x%U_1(:,:,kd,ld) = iu/nm(kd,ld)
2384 end do
2385 end if
2386 ! include order 2
2387 if (u_style.ge.2) then
2388 do kd = 1,grid_x%loc_n_r
2389 u_corr(:,:,kd,1) = t3_x(:,:,kd,1) + &
2390 &iu*n_frac(kd,ld)*t4_x(:,:,kd,1)
2391 u_corr(:,:,kd,2) = n_frac(kd,ld)/nm(kd,ld)*t1_x(:,:,kd,1)
2392 x%U_0(:,:,kd,ld) = x%U_0(:,:,kd,ld) + iu/nm(kd,ld)*&
2393 &u_corr(:,:,kd,1)
2394 x%U_1(:,:,kd,ld) = x%U_1(:,:,kd,ld) + iu/nm(kd,ld)*&
2395 &u_corr(:,:,kd,2)
2396 end do
2397 end if
2398 ! include order 3
2399 if (u_style.ge.3) then
2400 do kd = 1,grid_x%loc_n_r
2401 u_corr(:,:,kd,1) = iu*n_frac(kd,ld)*&
2402 &(-t5_x(:,:,kd,1)-iu*n_frac(kd,ld)*t6_x(:,:,kd,1))
2403 x%U_0(:,:,kd,ld) = x%U_0(:,:,kd,ld) + (iu/nm(kd,ld))**2*&
2404 &u_corr(:,:,kd,1)
2405 end do
2406 end if
2407 if (u_style.ge.4) then
2408 call writo('The geodesic perturbation U is implemented up to &
2409 &order 3',warning=.true.)
2410 end if
2411 end do
2412
2413 ! calculate DU_0 and DU_1, starting with first part
2414 select case (eq_style)
2415 case (1) ! VMEC
2416 do ld = 1,x%n_mod
2417 ! calculate order 1 of DU_0 and DU_1
2418 if (u_style.ge.1) then
2419 x%DU_0(:,:,:,ld) = -t2_x(:,:,:,2)
2420 x%DU_1(:,:,:,ld) = 0._dp
2421 end if
2422 ! include order 2
2423 if (u_style.ge.2) then
2424 do kd = 1,grid_x%loc_n_r
2425 u_corr(:,:,kd,1) = t3_x(:,:,kd,2) + &
2426 &iu*n_frac(kd,ld)*t4_x(:,:,kd,2)
2427 u_corr(:,:,kd,2) = n_frac(kd,ld)/nm(kd,ld)*&
2428 &t1_x(:,:,kd,2)
2429 x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2430 &iu/nm(kd,ld)*u_corr(:,:,kd,1)
2431 x%DU_1(:,:,kd,ld) = x%DU_1(:,:,kd,ld) + &
2432 &iu/nm(kd,ld)*u_corr(:,:,kd,2)
2433 end do
2434 end if
2435 ! include order 3
2436 if (u_style.ge.3) then
2437 do kd = 1,grid_x%loc_n_r
2438 u_corr(:,:,kd,1) = iu*n_frac(kd,ld)*&
2439 &(-t5_x(:,:,kd,2)-iu*n_frac(kd,ld)*&
2440 &t6_x(:,:,kd,2))
2441 x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2442 &(iu/nm(kd,ld))**2*u_corr(:,:,kd,1)
2443 end do
2444 end if
2445 if (u_style.ge.4) then
2446 call writo('The geodesic perturbation U is implemented &
2447 &only up to order 3',warning=.true.)
2448 end if
2449 end do
2450#if ldebug
2451 if (ltest) then
2452 call writo('Test calculation of DU')
2453 if(get_log(.false.,ind=.true.)) then
2454 ierr = test_du()
2455 chckerr('')
2456 call pause_prog(ind=.true.)
2457 end if
2458 end if
2459#endif
2460 case (2) ! HELENA
2461 ! reset pointers
2462 if (use_pol_flux_f) then
2463 ang_par_f => grid_x%theta_F
2464 else
2465 ang_par_f => grid_x%zeta_F
2466 end if
2467 ! set up helper variable for derivative
2468 ! numerically derive U_0 and U_1
2469 do ld = 1,x%n_mod
2470 do kd = 1,grid_x%loc_n_r
2471 do jd = 1,grid_x%n(2)
2472 ierr = spline(ang_par_f(:,jd,kd),x%U_0(:,jd,kd,ld),&
2473 &ang_par_f(:,jd,kd),x%DU_0(:,jd,kd,ld),&
2474 &ord=norm_disc_prec_x,deriv=1)
2475 chckerr('')
2476 ierr = spline(ang_par_f(:,jd,kd),x%U_1(:,jd,kd,ld),&
2477 &ang_par_f(:,jd,kd),x%DU_1(:,jd,kd,ld),&
2478 &ord=norm_disc_prec_x,deriv=1)
2479 chckerr('')
2480 end do
2481 end do
2482 end do
2483 end select
2484
2485 ! add second part of derivatives ~ n_frac
2486 do kd = 1,grid_x%loc_n_r
2487 do ld = 1,x%n_mod
2488 x%DU_0(:,:,kd,ld) = x%DU_0(:,:,kd,ld) + &
2489 &iu*n_frac(kd,ld)*x%U_0(:,:,kd,ld)
2490 x%DU_1(:,:,kd,ld) = x%DU_1(:,:,kd,ld) + &
2491 &iu*n_frac(kd,ld)*x%U_1(:,:,kd,ld)
2492 end do
2493 end do
2494
2495 ! clean up
2496 nullify(ang_par_f,ang_geo_f)
2497 nullify(j,d3j)
2498 nullify(g13,d3g13)
2499 nullify(g23,d3g23)
2500 nullify(g33,d3g33)
2501 nullify(h12,d3h12)
2502 nullify(h22,d1h22,d3h22,d13h22,d33h22)
2503 nullify(h23,d1h23,d3h23,d13h23,d33h23)
2504 select case (x_grid_style)
2505 case (1) ! equilibrium
2506 ! do nothing
2507 case (2,3) ! solution or enriched
2508 deallocate(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x,t6_x)
2509 end select
2510 nullify(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x,t6_x)
2511#if ldebug
2512 contains
2513 ! Test calculation of DU by deriving U numerically.
2514 integer function test_du() result(ierr)
2515 use num_vars, only: norm_disc_prec_x
2516
2517 character(*), parameter :: rout_name = 'test_DU'
2518
2519 ! local variables
2520 integer :: jd, kd, ld ! counters
2521 complex(dp), allocatable :: du_0(:,:,:) ! alternative calculation for DU_0
2522 complex(dp), allocatable :: du_1(:,:,:) ! alternative calculation for DU_1
2523 character(len=max_str_ln) :: file_name ! name of plot file
2524 character(len=max_str_ln) :: description ! description of plot
2525
2526 ! initialize ierr
2527 ierr = 0
2528
2529 ! warning
2530 call writo('This test is representable only if there are enough &
2531 &parallel points in the grid!',warning=.true.)
2532
2533 ! output
2534 call writo('Going to test whether DU is consistent with U')
2535 call lvl_ud(1)
2536
2537 ! set up DU_0 and DU_1
2538 allocate(du_0(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2539 allocate(du_1(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2540
2541 ! reset pointers
2542 if (use_pol_flux_f) then
2543 ang_par_f => grid_x%theta_F
2544 else
2545 ang_par_f => grid_x%zeta_F
2546 end if
2547
2548 ! loop over all modes
2549 do ld = 1,x%n_mod
2550 ! derivate numerically
2551 do kd = 1,grid_x%loc_n_r
2552 do jd = 1,grid_x%n(2)
2553 ierr = spline(ang_par_f(:,jd,kd),x%U_0(:,jd,kd,ld),&
2554 &ang_par_f(:,jd,kd),du_0(:,jd,kd),&
2555 &ord=norm_disc_prec_x,deriv=1)
2556 chckerr('')
2557 ierr = spline(ang_par_f(:,jd,kd),x%U_1(:,jd,kd,ld),&
2558 &ang_par_f(:,jd,kd),du_1(:,jd,kd),&
2559 &ord=norm_disc_prec_x,deriv=1)
2560 chckerr('')
2561 end do
2562 end do
2563
2564 ! set some variables
2565 file_name = 'TEST_RE_DU_0_'//trim(i2str(ld))
2566 description = 'Verifying DU_0 by deriving U_0 numerically for &
2567 &mode '//trim(i2str(ld)//', real part')
2568
2569 ! plot difference for RE DU_0
2570 call plot_diff_hdf5(rp(du_0),rp(x%DU_0(:,:,:,ld)),&
2571 &file_name,descr=description,output_message=.true.)
2572
2573 ! set some variables
2574 file_name = 'TEST_IM_DU_0_'//trim(i2str(ld))
2575 description = 'Verifying DU_0 by deriving U_0 numerically for &
2576 &mode '//trim(i2str(ld)//', imaginary part')
2577
2578 ! plot difference for IM DU_0
2579 call plot_diff_hdf5(ip(du_0),ip(x%DU_0(:,:,:,ld)),&
2580 &file_name,descr=description,output_message=.true.)
2581
2582 ! set some variables
2583 file_name = 'TEST_RE_DU_1_'//trim(i2str(ld))
2584 description = 'Verifying DU_1 by deriving U_1 numerically for &
2585 &mode '//trim(i2str(ld)//', real part')
2586
2587 ! plot difference for RE DU_1
2588 call plot_diff_hdf5(rp(du_1),rp(x%DU_1(:,:,:,ld)),&
2589 &file_name,descr=description,output_message=.true.)
2590
2591 ! set some variables
2592 file_name = 'TEST_IM_DU_1_'//trim(i2str(ld))
2593 description = 'Verifying DU_1 by deriving U_1 numerically for &
2594 &mode '//trim(i2str(ld)//', imaginary part')
2595
2596 ! plot difference for IM DU_1
2597 call plot_diff_hdf5(ip(du_1),ip(x%DU_1(:,:,:,ld)),&
2598 &file_name,descr=description,output_message=.true.)
2599 end do
2600
2601 ! user output
2602 call lvl_ud(-1)
2603 call writo('Test complete')
2604 end function test_du
2605#endif
2606 end function calc_u
2607
2608 !> calculate \f$\widetilde{PV}_{k,m}^i\f$ (pol. flux) or
2609 !! \f$\widetilde{PV}_{l,n}^i\f$ (tor. flux) at all equilibrium \c loc_n_r
2610 !! values.
2611 !!
2612 !! Like in calc_U(), use is made of variables \c Ti that are set up in the
2613 !! equilibrium grid and are afterwards converted to \c Ti_X in the
2614 !! perturbation grid, which needs to have the same angular coordinates as
2615 !! the equilibrium grid:
2616 !! \f[\begin{aligned}
2617 !! T_1 &= \frac{h^{22}}{\mu_0} g_{33} \\
2618 !! T_2 &= \mathcal{J}S +
2619 !! \mu_0 \sigma \frac{g_{33}}{\mathcal{J}} h^{22} \\
2620 !! T_3 &= \frac{\sigma}{\mathcal{J}} T_2 \\
2621 !! T_4 &= \frac{1}{\mu_0} \mathcal{J}^2 h^{22} \\
2622 !! T_5 &= 2 p' \kappa_n \ ,
2623 !! \end{aligned}\f]
2624 !! with \f$T_i\f$ = \c Ti.
2625 !!
2626 !! The interpolated Ti_X are then used to calculate
2627 !! \f$\widetilde{PV}_{k,m}^i\f$:
2628 !! \f[\begin{aligned}
2629 !! \widetilde{PV}_{k,m}^0 &= T_1(DU_k^{0*} - T_2)(DU_m^0 - T_2) -
2630 !! T_3 + (nq-m)(nq-k)T_4 - T_5 \\
2631 !! \widetilde{PV}_{k,m}^1 &= T_1(DU_k^{0*} - T_2) DU_m^1 \\
2632 !! \widetilde{PV}_{k,m}^2 &= T_1 DU_k^{1*} DU_m^1 \ ,
2633 !! \end{aligned}\f]
2634 !! where
2635 !! \f[DU_k^i = i (nq-m) U_k^i + \frac{\partial U_k^i}{\partial\theta} . \f]
2636 !!
2637 !! This is valid for poloidal Flux coordinates and where \f$n\f$ is to be
2638 !! replaced by \f$m\f$ and \f$(nq-m)\f$ by \f$(n-\iota m)\f$ for toroidal
2639 !! Flux coordinates.
2640 !!
2641 !! \see See \cite weyens2014theory .
2642 !!
2643 !! \return ierr
2644 integer function calc_pv(grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,lim_sec_X) &
2645 &result(ierr)
2646 use num_vars, only: use_pol_flux_f, norm_disc_prec_x, x_grid_style
2647 use eq_vars, only: vac_perm
2648 use num_utilities, only: c, spline
2649 use x_utilities, only: is_necessary_x
2650 use x_vars, only: n_mod_x
2651
2652 character(*), parameter :: rout_name = 'calc_PV'
2653
2654 ! use input / output
2655 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
2656 type(grid_type), intent(in) :: grid_x !< perturbation grid
2657 type(eq_1_type), intent(in), target :: eq_1 !< flux equilibrium
2658 type(eq_2_type), intent(in), target :: eq_2 !< metric equilibrium
2659 type(x_1_type), intent(in) :: x_a !< vectorial perturbation variables of dimension 2
2660 type(x_1_type), intent(in) :: x_b !< vectorial perturbation variables of dimension 2
2661 type(x_2_type), intent(inout) :: x !< tensorial perturbation variables
2662 integer, intent(in), optional :: lim_sec_x(2,2) !< limits of \c m_X (pol flux) or \c n_X (tor flux) for both dimensions
2663
2664 ! local variables
2665 integer :: m, k ! counters
2666 integer :: id, jd, kd ! counters
2667 integer :: c_loc(2) ! local c for symmetric and asymmetric variables
2668 logical :: calc_this(2) ! whether this combination needs to be calculated
2669
2670 ! jacobian
2671 real(dp), pointer :: j(:,:,:) ! jac
2672 ! lower metric factors
2673 real(dp), pointer :: g33(:,:,:) ! h_theta,theta or h_zeta,zeta
2674 ! upper metric factors
2675 real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2676 ! helper variables
2677 real(dp), pointer :: q_saf(:), rot_t(:) ! safety factor and rotational transform in X grid
2678 real(dp), allocatable :: fac_n(:), fac_m(:) ! multiplicative factors for n and m
2679 ! PV factors
2680 real(dp), allocatable, target :: t1(:,:,:) ! h22/mu_0 g33
2681 real(dp), allocatable, target :: t2(:,:,:) ! JS + mu_0 sigma g33/J h22
2682 real(dp), allocatable, target :: t3(:,:,:) ! sigma/J T2
2683 real(dp), allocatable, target :: t4(:,:,:) ! 1/mu_0 J^2 h22
2684 real(dp), allocatable, target :: t5(:,:,:) ! 2 p' kappa_n
2685 real(dp), pointer :: t1_x(:,:,:) ! T1 in X grid and par. deriv.
2686 real(dp), pointer :: t2_x(:,:,:) ! T2 in X grid and par. deriv.
2687 real(dp), pointer :: t3_x(:,:,:) ! T3 in X grid and par. deriv.
2688 real(dp), pointer :: t4_x(:,:,:) ! T4 in X grid and par. deriv.
2689 real(dp), pointer :: t5_x(:,:,:) ! T5 in X grid and par. deriv.
2690
2691 ! initialize ierr
2692 ierr = 0
2693
2694 ! allocate variables
2695 ! helper variables
2696 allocate(q_saf(grid_x%loc_n_r))
2697 allocate(rot_t(grid_x%loc_n_r))
2698 allocate(fac_n(grid_x%loc_n_r),fac_m(grid_x%loc_n_r))
2699 ! PV factors
2700 allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2701 allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2702 allocate(t3(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2703 allocate(t4(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2704 allocate(t5(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2705 select case (x_grid_style)
2706 case (1) ! equilibrium
2707 ! do nothing
2708 case (2,3) ! solution or enriched
2709 allocate(q_saf(grid_x%loc_n_r))
2710 allocate(rot_t(grid_x%loc_n_r))
2711 allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2712 allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2713 allocate(t3_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2714 allocate(t4_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2715 allocate(t5_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2716 end select
2717
2718 ! set pointers
2719 j => eq_2%jac_FD(:,:,:,0,0,0)
2720 g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2721 h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2722
2723 ! set up PV factors in eq grid
2724 t1 = h22/(vac_perm*g33)
2725 t2 = j*eq_2%S + vac_perm*eq_2%sigma*g33/(j*h22)
2726 t3 = eq_2%sigma/j*t2
2727 t4 = 1._dp/(vac_perm*j**2*h22)
2728 do kd = 1,grid_eq%loc_n_r
2729 t5(:,:,kd) = 2*eq_1%pres_FD(kd,1)*eq_2%kappa_n(:,:,kd)
2730 end do
2731
2732 select case (x_grid_style)
2733 case (1) ! equilibrium
2734 ! point
2735 q_saf => eq_1%q_saf_FD(:,0)
2736 rot_t => eq_1%rot_t_FD(:,0)
2737 t1_x => t1
2738 t2_x => t2
2739 t3_x => t3
2740 t4_x => t4
2741 t5_x => t5
2742 case (2,3) ! solution or enriched
2743 ! interpolate
2744 ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
2745 &grid_x%loc_r_F,q_saf,ord=norm_disc_prec_x)
2746 chckerr('')
2747 ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
2748 &grid_x%loc_r_F,rot_t,ord=norm_disc_prec_x)
2749 chckerr('')
2750 do jd = 1,grid_x%n(2)
2751 do id = 1,grid_x%n(1)
2752 ierr = spline(grid_eq%loc_r_F,t1(id,jd,:),&
2753 &grid_x%loc_r_F,t1_x(id,jd,:),&
2754 &ord=norm_disc_prec_x)
2755 chckerr('')
2756 ierr = spline(grid_eq%loc_r_F,t2(id,jd,:),&
2757 &grid_x%loc_r_F,t2_x(id,jd,:),&
2758 &ord=norm_disc_prec_x)
2759 chckerr('')
2760 ierr = spline(grid_eq%loc_r_F,t3(id,jd,:),&
2761 &grid_x%loc_r_F,t3_x(id,jd,:),&
2762 &ord=norm_disc_prec_x)
2763 chckerr('')
2764 ierr = spline(grid_eq%loc_r_F,t4(id,jd,:),&
2765 &grid_x%loc_r_F,t4_x(id,jd,:),&
2766 &ord=norm_disc_prec_x)
2767 chckerr('')
2768 ierr = spline(grid_eq%loc_r_F,t5(id,jd,:),&
2769 &grid_x%loc_r_F,t5_x(id,jd,:),&
2770 &ord=norm_disc_prec_x)
2771 chckerr('')
2772 end do
2773 end do
2774 end select
2775
2776 ! set up fac_n and fac_m
2777 if (use_pol_flux_f) then
2778 fac_n = q_saf
2779 fac_m = 1.0_dp
2780 else
2781 fac_n = 1.0_dp
2782 fac_m = rot_t
2783 end if
2784
2785 ! loop over all modes
2786 do m = 1,x_b%n_mod
2787 do k = 1,x_a%n_mod
2788 ! check whether modes are correct
2789 do kd = 1,grid_x%loc_n_r
2790 if (x%n_1(kd,k).ne.x_a%n(kd,k) .or. &
2791 &x%n_2(kd,m).ne.x_b%n(kd,m) .or. &
2792 &x%m_1(kd,k).ne.x_a%m(kd,k) .or. &
2793 &x%m_2(kd,m).ne.x_b%m(kd,m)) then
2794 ierr = 1
2795 chckerr('Modes do not match')
2796 end if
2797 end do
2798
2799 ! check whether mode combination needs to be calculated
2800 calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
2801 calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
2802
2803 ! set up c_loc
2804 c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
2805 c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
2806
2807 ! calculate PV_0
2808 if (calc_this(1)) then
2809 do kd = 1,grid_x%loc_n_r
2810 x%PV_0(:,:,kd,c_loc(1)) = t1_x(:,:,kd)*&
2811 &(x_b%DU_0(:,:,kd,m) - t2_x(:,:,kd)) * &
2812 &(conjg(x_a%DU_0(:,:,kd,k)) - t2_x(:,:,kd)) - &
2813 &t3_x(:,:,kd) - t5_x(:,:,kd) + t4_x(:,:,kd)*&
2814 &(x_b%n(kd,m)*fac_n(kd)-x_b%m(kd,m)*fac_m(kd))*&
2815 &(x_a%n(kd,k)*fac_n(kd)-x_a%m(kd,k)*fac_m(kd))
2816 end do
2817 end if
2818
2819 ! calculate PV_1
2820 if (calc_this(2)) then
2821 x%PV_1(:,:,:,c_loc(2)) = t1_x * x_b%DU_1(:,:,:,m) * &
2822 &(conjg(x_a%DU_0(:,:,:,k))-t2_x)
2823 end if
2824
2825 ! calculate PV_2
2826 if (calc_this(1)) then
2827 x%PV_2(:,:,:,c_loc(1)) = t1_x * x_b%DU_1(:,:,:,m) * &
2828 &conjg(x_a%DU_1(:,:,:,k))
2829 end if
2830 end do
2831 end do
2832
2833 ! clean up
2834 nullify(j)
2835 nullify(g33)
2836 nullify(h22)
2837 select case (x_grid_style)
2838 case (1) ! equilibrium
2839 ! do nothing
2840 case (2,3) ! solution or enriched
2841 deallocate(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x)
2842 end select
2843 nullify(q_saf,rot_t,t1_x,t2_x,t3_x,t4_x,t5_x)
2844 end function calc_pv
2845
2846 !> calculate \f$\widetilde{KV}_{k,m}^i\f$ (pol. flux) or
2847 !! \f$\widetilde{KV}_{l,n}^i\f$ (tor. flux) at all equilibrium \c loc_n_r
2848 !! values.
2849 !!
2850 !! Like in calc_U(), use is made of variables \c Ti that are set up in the
2851 !! equilibrium grid and are afterwards converted to \c Ti_X in the
2852 !! perturbation grid, which needs to have the same angular coordinates as
2853 !! the equilibrium grid:
2854 !! \f[\begin{aligned}
2855 !! T_1 &= \rho \mathcal{J}^2 \frac{h^{22}}{\mu_0} g_{33} \\
2856 !! T_2 &= \rho \frac{1}{h^{22}} \ ,
2857 !! \end{aligned}\f]
2858 !! with \f$T_i\f$ = \c Ti.
2859 !!
2860 !! The interpolated Ti_X are then used to calculate
2861 !! \f$\widetilde{KV}_{k,m}^i\f$:
2862 !! \f[\begin{aligned}
2863 !! \widetilde{KV}_{k,m}^0 &= T_1 U_k^{0*} U_m^0 + T_2 \\
2864 !! \widetilde{KV}_{k,m}^1 &= T_1 U_k^{0*} U_m^1 \\
2865 !! \widetilde{KV}_{k,m}^2 &= T_1 U_k^{1*} U_m^1 \ ,
2866 !! \end{aligned}\f]
2867 !! where
2868 !! \f[DU_k^i = i (nq-m) U_k^i + \frac{\partial U_k^i}{\partial\theta} . \f]
2869 !!
2870 !! This is valid for poloidal Flux coordinates and where \f$n\f$ is to be
2871 !! replaced by \f$m\f$ and \f$(nq-m)\f$ by \f$(n-\iota m)\f$ for toroidal
2872 !! Flux coordinates.
2873 !!
2874 !! \see See \cite weyens2014theory .
2875 !!
2876 !! \return ierr
2877 integer function calc_kv(grid_eq,grid_X,eq_1,eq_2,X_a,X_b,X,lim_sec_X) &
2878 &result(ierr)
2879 use num_vars, only: k_style, norm_disc_prec_x, x_grid_style
2880 use num_utilities, only: c, spline
2881 use x_utilities, only: is_necessary_x
2882 use x_vars, only: n_mod_x
2883
2884 character(*), parameter :: rout_name = 'calc_KV'
2885
2886 ! use input / output
2887 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
2888 type(grid_type), intent(in) :: grid_x !< perturbation grid
2889 type(eq_1_type), intent(in) :: eq_1 !< flux equilibrium
2890 type(eq_2_type), intent(in), target :: eq_2 !< metric equilibrium
2891 type(x_1_type), intent(in) :: x_a !< vectorial perturbation variables of dimension 2
2892 type(x_1_type), intent(in) :: x_b !< vectorial perturbation variables of dimension 2
2893 type(x_2_type), intent(inout) :: x !< tensorial perturbation variables
2894 integer, intent(in), optional :: lim_sec_x(2,2) !< limits of \c m_X (pol flux) or \c n_X (tor flux) for both dimensions
2895
2896 ! local variables
2897 integer :: m, k ! counters
2898 integer :: id, jd, kd ! counters
2899 integer :: c_loc(2) ! local c for symmetric and asymmetric variables
2900 logical :: calc_this(2) ! whether this combination needs to be calculated
2901
2902 ! jacobian
2903 real(dp), pointer :: j(:,:,:) ! jac
2904 ! lower metric factors
2905 real(dp), pointer :: g33(:,:,:) ! h_theta,theta or h_zeta,zeta
2906 ! upper metric factors
2907 real(dp), pointer :: h22(:,:,:) ! h^psi,psi
2908 ! KV factors
2909 real(dp), allocatable, target :: t1(:,:,:) ! h22/mu_0 J g33
2910 real(dp), allocatable, target :: t2(:,:,:) ! JS + mu_0 sigma g33/h22
2911 real(dp), pointer :: t1_x(:,:,:) ! T1 in X grid and par. deriv.
2912 real(dp), pointer :: t2_x(:,:,:) ! T2 in X grid and par. deriv.
2913
2914 ! initialize ierr
2915 ierr = 0
2916
2917 ! allocate variables
2918 ! KV factors
2919 allocate(t1(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2920 allocate(t2(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
2921 select case (x_grid_style)
2922 case (1) ! equilibrium
2923 ! do nothing
2924 case (2,3) ! solution or enriched
2925 allocate(t1_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2926 allocate(t2_x(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
2927 end select
2928
2929 ! set pointers
2930 j => eq_2%jac_FD(:,:,:,0,0,0)
2931 g33 => eq_2%g_FD(:,:,:,c([3,3],.true.),0,0,0)
2932 h22 => eq_2%h_FD(:,:,:,c([2,2],.true.),0,0,0)
2933
2934 ! set up KV factors in eq grid
2935 do kd = 1,grid_eq%loc_n_r
2936 t1(:,:,kd) = eq_1%rho(kd)*j(:,:,kd)**2*h22(:,:,kd)/g33(:,:,kd)
2937 t2(:,:,kd) = eq_1%rho(kd)/h22(:,:,kd)
2938 end do
2939
2940 select case (x_grid_style)
2941 case (1) ! equilibrium
2942 ! point
2943 t1_x => t1
2944 t2_x => t2
2945 case (2,3) ! solution or enriched
2946 ! interpolate
2947 do jd = 1,grid_x%n(2)
2948 do id = 1,grid_x%n(1)
2949 ierr = spline(grid_eq%loc_r_F,t1(id,jd,:),&
2950 &grid_x%loc_r_F,t1_x(id,jd,:),&
2951 &ord=norm_disc_prec_x)
2952 chckerr('')
2953 ierr = spline(grid_eq%loc_r_F,t2(id,jd,:),&
2954 &grid_x%loc_r_F,t2_x(id,jd,:),&
2955 &ord=norm_disc_prec_x)
2956 chckerr('')
2957 end do
2958 end do
2959 end select
2960
2961 ! loop over all modes
2962 do m = 1,x_b%n_mod
2963 do k = 1,x_a%n_mod
2964 ! check whether modes are correct
2965 do kd = 1,grid_x%loc_n_r
2966 if (x%n_1(kd,k).ne.x_a%n(kd,k) .or. &
2967 &x%n_2(kd,m).ne.x_b%n(kd,m) .or. &
2968 &x%m_1(kd,k).ne.x_a%m(kd,k) .or. &
2969 &x%m_2(kd,m).ne.x_b%m(kd,m)) then
2970 ierr = 1
2971 chckerr('Modes do not match')
2972 end if
2973 end do
2974
2975 ! check whether mode combination needs to be calculated
2976 calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
2977 calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
2978
2979 ! set up c_loc
2980 c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
2981 c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
2982
2983 select case (k_style)
2984 case (1) ! normalization of full perpendicular component
2985 ! calculate KV_0
2986 if (calc_this(1)) then
2987 x%KV_0(:,:,:,c_loc(1)) = t2_x + t1_x * &
2988 &x_b%U_0(:,:,:,m) * conjg(x_a%U_0(:,:,:,k))
2989 end if
2990
2991 ! calculate KV_1
2992 if (calc_this(2)) then
2993 x%KV_1(:,:,:,c_loc(2)) = t1_x * &
2994 &x_b%U_1(:,:,:,m) * conjg(x_a%U_0(:,:,:,k))
2995 end if
2996
2997 ! calculate KV_2
2998 if (calc_this(1)) then
2999 x%KV_2(:,:,:,c_loc(1)) = t1_x * &
3000 &x_b%U_1(:,:,:,m) * conjg(x_a%U_1(:,:,:,k))
3001 end if
3002 case (2) ! normalization of only normal component
3003 ! calculate KV_0
3004 if (calc_this(1)) then
3005 x%KV_0(:,:,:,c_loc(1)) = t2_x
3006 end if
3007
3008 ! calculate KV_1
3009 if (calc_this(2)) then
3010 x%KV_1(:,:,:,c_loc(2)) = 0._dp
3011 end if
3012
3013 ! calculate KV_2
3014 if (calc_this(1)) then
3015 x%KV_2(:,:,:,c_loc(1)) = 0._dp
3016 end if
3017 end select
3018 end do
3019 end do
3020
3021 ! clean up
3022 nullify(j)
3023 nullify(g33)
3024 nullify(h22)
3025 select case (x_grid_style)
3026 case (1) ! equilibrium
3027 ! do nothing
3028 case (2,3) ! solution or enriched
3029 deallocate(t1_x,t2_x)
3030 end select
3031 nullify(t1_x,t2_x)
3032 end function calc_kv
3033
3034 !> Calculate the magnetic integrals from \f$\widetilde{PV}_{k,m}^i\f$ and
3035 !! \f$\widetilde{KV}_{k,m}^i\f$ in an equidistant grid where the step size
3036 !! can vary depending on the normal coordinate.
3037 !!
3038 !! All the variables should thus be field-line oriented. The result is saved
3039 !! in the first index of the X variables, the other can be ignored.
3040 !!
3041 !! Using \c prev_style, the results of this calculation on X can be combined
3042 !! with the ones already present in \c X_int.
3043 !! - \c prev_style=0: Overwrite [def].
3044 !! - \c prev_style=1: Add to current integral.
3045 !! - \c prev_style=2: Divide by 2 and add to current integral. Also modify the
3046 !! indices of the current integral:
3047 !! - for \c magn_int_style = 1: <tt>1 2 2 2 2 2 1 -> 2 2 2 2 2 2</tt>,
3048 !! - for \c magn_int_style = 2: <tt>1 3 3 2 3 3 1 -> 3 2 3 3 2 3</tt>.
3049 !! - \c prev_style=3: Add to current integral. Also modify the indices of
3050 !! the current integral as in \c prev_style = 2.
3051 !! - else: ignore previous magnetic integrals.
3052 !!
3053 !! Therefore, it is to be used in order. For example:
3054 !! -# \f$\phantom{\rightarrow}\f$ (R=1,E=1): style 0,
3055 !! -# \f$\rightarrow\f$ (R=1,E=2): style 1,
3056 !! -# \f$\rightarrow\f$ (R=2,E=1): style 2,
3057 !! -# \f$\rightarrow\f$ (R=2,E=2): style 3.
3058 !!
3059 !! Afterwards, it would cycle through 2 and 3, as style 0 and 1 are only for
3060 !! the first Richardson level.
3061 !!
3062 !! \note
3063 !! -# The variable type for the integrated tensorial perturbation variables
3064 !! is also \c X_2, but it is assumed that the first index has dimension 1,
3065 !! the rest is ignored.
3066 !! -# Simpson's 3/8 rule converges faster than the trapezoidal rule, but
3067 !! normally needs a better starting point (i.e. higher \c min_n_par_X)
3068 !! -# For alpha style 1, Weyl's lemma is used to convert the surface
3069 !! integral into a line integral along the magnetic field. This means
3070 !! that the resulting line integral still needs to be multiplied by
3071 !! \f$\frac{2\pi}{M}\f$, where \f$M\f$ is the number of times the poloidal
3072 !! or toroidal angle completes a full \f$2\pi\f$ tour, depending on \c
3073 !! use_pol_flux_F \cite Helander2014.
3074 !! -# For alpha style 2, this factor is just the step size in alpha.
3075 !!
3076 !! \see See x_vars.x_2_type.
3077 !!
3078 !! \return ierr
3079 integer function calc_magn_ints(grid_eq,grid_X,eq,X,X_int,prev_style,&
3080 &lim_sec_X) result(ierr)
3081 use num_vars, only: use_pol_flux_f, norm_disc_prec_x, magn_int_style, &
3082 &alpha_style, x_grid_style
3083 use x_utilities, only: is_necessary_x
3084 use x_vars, only: n_mod_x
3085 use num_utilities, only: c, spline
3086 use grid_vars, only: min_par_x, max_par_x, n_alpha
3087 use rich_vars, only: rich_lvl
3088
3089 character(*), parameter :: rout_name = 'calc_magn_ints'
3090
3091 ! input / output
3092 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
3093 type(grid_type), intent(in) :: grid_x !< perturbation grid
3094 type(eq_2_type), intent(in), target :: eq !< metric equilibrium
3095 type(x_2_type), intent(in) :: x !< tensorial perturbation variables
3096 type(x_2_type), intent(inout) :: x_int !< interpolated tensorial perturbation variables, containing all mode combinations
3097 integer, intent(in), optional :: prev_style !< style to treat X_prev
3098 integer, intent(in), optional :: lim_sec_x(2,2) !< limits of \c m_X (pol flux) or \c n_X (tor flux) for both dimensions
3099
3100 ! local variables, not used in subroutines
3101 logical :: calc_this(2) ! whether this combination needs to be calculated
3102 integer :: k, m ! counters
3103 integer :: id, jd, kd ! counters
3104 integer :: c_loc(2) ! local c for symmetric and asymmetric variables
3105 integer :: c_tot(2) ! total c for symmetric and asymmetric variables
3106 integer :: prev_style_loc ! local prev_style
3107 real(dp) :: alpha_fac ! factor for alpha (style 1: Weyl factor, style 2: step size)
3108 real(dp) :: prev_mult_fac ! multiplicative factor for previous style
3109 real(dp), allocatable :: step_size(:,:) ! step size for every field line and normal point
3110 real(dp), pointer :: ang_par_f(:,:,:) => null() ! parallel angle in flux coordinates
3111 complex(dp), allocatable :: v_int_work(:,:) ! work V_int
3112
3113 ! local variables, also used in subroutines
3114 integer :: nr_int_regions ! nr. of integration regions
3115 integer, allocatable :: int_dims(:,:) ! dimension information of integration regions
3116 real(dp), allocatable :: int_facs(:) ! integration factors for each of the regions
3117 complex(dp), allocatable :: j_exp_ang(:,:,:) ! J * exponential of Flux parallel angle
3118
3119 ! Makes use of nr_int_regions, int_dims, int_facs and J_exp_ang
3120
3121 ! jacobian
3122 real(dp), pointer :: j(:,:,:) ! jac
3123
3124 ! initialize ierr
3125 ierr = 0
3126
3127 ! user output
3128 call writo('Calculate field-line averages')
3129 call lvl_ud(1)
3130
3131 ! set up local prev_style
3132 prev_style_loc = 0
3133 if (present(prev_style)) prev_style_loc = prev_style
3134
3135 ! allocate variables
3136 allocate(j_exp_ang(grid_x%n(1),grid_x%n(2),grid_x%loc_n_r))
3137 select case (x_grid_style)
3138 case (1) ! equilibrium
3139 ! point
3140 j => eq%jac_FD(:,:,:,0,0,0)
3141 case (2,3) ! solution or enriched
3142 ! allocate
3143 allocate(j(grid_eq%n(1),grid_eq%n(2),grid_x%loc_n_r))
3144
3145 ! interpolate
3146 do jd = 1,grid_x%n(2)
3147 do id = 1,grid_x%n(1)
3148 ierr = spline(grid_eq%loc_r_F,eq%jac_FD(id,jd,:,0,0,0),&
3149 &grid_x%loc_r_F,j(id,jd,:),ord=norm_disc_prec_x)
3150 chckerr('')
3151 end do
3152 end do
3153 end select
3154
3155 ! set up parallel angle in flux coordinates
3156 if (use_pol_flux_f) then
3157 ang_par_f => grid_x%theta_F
3158 else
3159 ang_par_f => grid_x%zeta_F
3160 end if
3161
3162 ! set up alpha factor
3163 select case (alpha_style)
3164 case (1) ! single field line, multiple turns
3165 alpha_fac = (2*pi)**2/((max_par_x-min_par_x)*pi) ! Weyl factor
3166 case (2) ! multiple field lines, single turns
3167 alpha_fac = (2*pi)/n_alpha ! grid size
3168 end select
3169
3170 ! set up integration variables
3171
3172 ! set nr_int_regions
3173 select case (magn_int_style)
3174 case (1) ! Trapezoidal rule
3175 select case (prev_style_loc)
3176 case (2,3) ! change indices
3177 nr_int_regions = 1
3178 case default ! don't change indices
3179 nr_int_regions = 2
3180 end select
3181 case (2) ! Simpson's 3/8 rule
3182 select case (prev_style_loc)
3183 case (2,3) ! change indices
3184 nr_int_regions = 3
3185 case default ! don't change indices
3186 nr_int_regions = 4
3187 end select
3188 end select
3189
3190 ! set up integration factors and dimensions
3191 allocate(int_facs(nr_int_regions))
3192 allocate(int_dims(nr_int_regions,3))
3193 select case (magn_int_style)
3194 case (1) ! Trapezoidal rule
3195 select case (prev_style_loc)
3196 case (2,3) ! change indices
3197 int_facs = 0.5_dp*[2]
3198 int_dims(1,:) = [1,grid_x%n(1),1] ! region 1: all points
3199 case default ! don't change indices
3200 int_facs = 0.5_dp*[1,2]
3201 int_dims(1,:) = [1,grid_x%n(1),grid_x%n(1)-1] ! region 1: first and last points
3202 int_dims(2,:) = [2,grid_x%n(1)-1,1] ! region 2: intermediate points
3203 end select
3204 case (2) ! Simpson's 3/8 rule
3205 select case (prev_style_loc)
3206 case (2,3) ! change indices
3207 int_facs = 0.375_dp*[3,2,3]
3208 int_dims(1,:) = [1,grid_x%n(1)-2,3] ! region 1: intermediate points ~ 3
3209 int_dims(2,:) = [2,grid_x%n(1)-1,3] ! region 2: intermediate points ~ 2
3210 int_dims(3,:) = [3,grid_x%n(1),3] ! region 3: intermediate points ~ 3
3211 case default ! don't change indices
3212 int_facs = 0.375_dp*[1,3,3,2]
3213 int_dims(1,:) = [1,grid_x%n(1),grid_x%n(1)-1] ! region 1: first and last points
3214 int_dims(2,:) = [2,grid_x%n(1)-2,3] ! region 2: intermediate points ~ 3
3215 int_dims(3,:) = [3,grid_x%n(1)-1,3] ! region 3: intermediate points ~ 3
3216 int_dims(4,:) = [4,grid_x%n(1)-3,3] ! region 4: intermediate points ~ 2
3217 end select
3218 end select
3219
3220 ! set up multiplicative factor for previous level
3221 select case (prev_style_loc)
3222 case (1,3) ! add to integral of previous equilibrium job
3223 prev_mult_fac = 1.0_dp
3224 case (2) ! add to integral of previous Richardson level / 2
3225 prev_mult_fac = 0.5_dp
3226 case default ! don't add, overwrite
3227 prev_mult_fac = 0.0_dp
3228 end select
3229
3230 ! set up step size
3231 allocate(step_size(grid_x%n(2),grid_x%loc_n_r))
3232 step_size = ang_par_f(2,:,:)-ang_par_f(1,:,:)
3233 if (rich_lvl.gt.1) step_size = step_size*0.5_dp ! only half of points present: step size is actually half
3234
3235 ! initialize local V_int
3236 allocate(v_int_work(grid_x%n(2),grid_x%loc_n_r))
3237
3238 ! loop over all modes
3239 do m = 1,x%n_mod(2)
3240 do k = 1,x%n_mod(1)
3241 ! check whether mode combination needs to be calculated
3242 calc_this(1) = is_necessary_x(.true.,[k,m],lim_sec_x)
3243 calc_this(2) = is_necessary_x(.false.,[k,m],lim_sec_x)
3244
3245 ! set up local and total indices
3246 c_loc(1) = c([k,m],.true.,n_mod_x,lim_sec_x)
3247 c_loc(2) = c([k,m],.false.,n_mod_x,lim_sec_x)
3248 c_tot(1) = c([lim_sec_x(1,1)-1+k,lim_sec_x(1,2)-1+m],.true.,&
3249 &n_mod_x)
3250 c_tot(2) = c([lim_sec_x(1,1)-1+k,lim_sec_x(1,2)-1+m],.false.,&
3251 &n_mod_x)
3252
3253 ! calculate J_exp_ang
3254 do kd = 1,grid_x%loc_n_r
3255 if (use_pol_flux_f) then
3256 j_exp_ang(:,:,kd) = j(:,:,kd)*&
3257 &exp(iu*(x%m_1(kd,k)-x%m_2(kd,m))*&
3258 &ang_par_f(:,:,kd))
3259 else
3260 j_exp_ang(:,:,kd) = j(:,:,kd)*&
3261 &exp(-iu*(x%n_1(kd,k)-x%n_2(kd,m))*&
3262 &ang_par_f(:,:,kd))
3263 end if
3264 end do
3265
3266 ! integrate PV_0 and KV_0
3267 if (calc_this(1)) then
3268 x_int%PV_0(1,:,:,c_tot(1)) = &
3269 &prev_mult_fac*x_int%PV_0(1,:,:,c_tot(1))
3270 x_int%KV_0(1,:,:,c_tot(1)) = &
3271 &prev_mult_fac*x_int%KV_0(1,:,:,c_tot(1))
3272 call calc_magn_int_loc(x%PV_0(:,:,:,c_loc(1))*alpha_fac,&
3273 &x_int%PV_0(1,:,:,c_tot(1)),v_int_work,step_size)
3274 call calc_magn_int_loc(x%KV_0(:,:,:,c_loc(1))*alpha_fac,&
3275 &x_int%KV_0(1,:,:,c_tot(1)),v_int_work,step_size)
3276 end if
3277
3278 ! integrate PV_1 and KV_1
3279 if (calc_this(2)) then
3280 x_int%PV_1(1,:,:,c_tot(2)) = &
3281 &prev_mult_fac*x_int%PV_1(1,:,:,c_tot(2))
3282 x_int%KV_1(1,:,:,c_tot(2)) = &
3283 &prev_mult_fac*x_int%KV_1(1,:,:,c_tot(2))
3284 call calc_magn_int_loc(x%PV_1(:,:,:,c_loc(2))*alpha_fac,&
3285 &x_int%PV_1(1,:,:,c_tot(2)),v_int_work,step_size)
3286 call calc_magn_int_loc(x%KV_1(:,:,:,c_loc(2))*alpha_fac,&
3287 &x_int%KV_1(1,:,:,c_tot(2)),v_int_work,step_size)
3288 end if
3289
3290 ! integrate PV_2 and KV_2
3291 if (calc_this(1)) then
3292 x_int%PV_2(1,:,:,c_tot(1)) = &
3293 &prev_mult_fac*x_int%PV_2(1,:,:,c_tot(1))
3294 x_int%KV_2(1,:,:,c_tot(1)) = &
3295 &prev_mult_fac*x_int%KV_2(1,:,:,c_tot(1))
3296 call calc_magn_int_loc(x%PV_2(:,:,:,c_loc(1))*alpha_fac,&
3297 &x_int%PV_2(1,:,:,c_tot(1)),v_int_work,step_size)
3298 call calc_magn_int_loc(x%KV_2(:,:,:,c_loc(1))*alpha_fac,&
3299 &x_int%KV_2(1,:,:,c_tot(1)),v_int_work,step_size)
3300 end if
3301 end do
3302 end do
3303
3304 ! clean up
3305 nullify(ang_par_f)
3306 select case (x_grid_style)
3307 case (1) ! equilibrium
3308 ! do nothing
3309 case (2,3) ! solution or enriched
3310 deallocate(j)
3311 end select
3312 nullify(j)
3313
3314 ! user output
3315 call lvl_ud(-1)
3316 contains
3317 ! Integrate local magnetic integral, adding to the previous result.
3318 !
3319 ! Makes use of nr_int_regions, int_dims, int_facs and J_exp_ang.
3320 !> \private
3321 subroutine calc_magn_int_loc(V,V_int,V_int_work,step_size)
3322 ! input / output
3323 complex(dp), intent(in) :: v(:,:,:) ! V
3324 complex(dp), intent(inout) :: v_int(:,:) ! integrated V
3325 complex(dp), intent(inout) :: v_int_work(:,:) ! workint variable
3326 real(dp) :: step_size(:,:) ! step size for every normal point
3327
3328 ! local variables
3329 integer :: id, kd, ld ! counters
3330
3331 ! loop over regions
3332 do ld = 1,nr_int_regions
3333 ! initialize
3334 v_int_work = 0
3335
3336 ! integrate
3337 do kd = 1,size(step_size,2)
3338 do id = int_dims(ld,1),int_dims(ld,2),int_dims(ld,3)
3339 v_int_work(:,kd) = v_int_work(:,kd) + &
3340 &j_exp_ang(id,:,kd)*v(id,:,kd)*step_size(:,kd)
3341 end do
3342 end do
3343
3344 ! scale by integration factor and update V_int
3345 v_int = v_int + v_int_work*int_facs(ld)
3346 end do
3347 end subroutine calc_magn_int_loc
3348 end function calc_magn_ints
3349
3350 !> Divides the perturbation jobs.
3351 !!
3352 !! This concerns the calculation of the magnetic integrals of blocks of
3353 !! tensorial perturbation variables. These are set up using the equilibrium
3354 !! and vectorial perturbation variables that are stored in memory, but only
3355 !! the integrated tensorial result is stored in memory, with negligible
3356 !! variables size.
3357 !!
3358 !! The size of the (k,m) pairs to be calculated is determined by looking at
3359 !! what fits in memory when the equilibrium and vectorial perturbation
3360 !! variables are stored in memory first.
3361 !!
3362 !! \return ierr
3363 integer function divide_x_jobs(arr_size) result(ierr)
3364 use num_vars, only: max_x_mem, n_procs, x_jobs_lims
3365 use x_vars, only: n_mod_x
3366 use x_utilities, only: calc_memory_x
3367
3368 character(*), parameter :: rout_name = 'divide_X_jobs'
3369
3370 ! input / output
3371 integer, intent(in) :: arr_size !< array size (using loc_n_r)
3372
3373 ! local variables
3374 integer :: n_div ! factor by which to divide the total size
3375 real(dp) :: mem_size ! approximation of memory required for X variables
3376 integer :: n_mod_block ! nr. of modes in block
3377 integer, allocatable :: n_mod_loc(:) ! number of modes per block
3378 character(len=max_str_ln) :: block_message ! message about how many blocks
3379 character(len=max_str_ln) :: err_msg ! error message
3380
3381 ! initialize ierr
3382 ierr = 0
3383
3384 ! user output
3385 call writo('Dividing the perturbation jobs')
3386 call lvl_ud(1)
3387
3388 ! calculate largest possible block of (k,m) values
3389 n_div = 0
3390 mem_size = huge(1._dp)
3391 do while (mem_size.gt.max_x_mem)
3392 n_div = n_div + 1
3393 n_mod_block = ceiling(n_mod_x*1._dp/n_div)
3394 ierr = calc_memory_x(2,arr_size,n_mod_block,mem_size)
3395 chckerr('')
3396 if (n_div.gt.n_mod_x) then
3397 ierr = 1
3398 err_msg = 'The memory limit is too low'
3399 chckerr(err_msg)
3400 end if
3401 end do
3402 if (n_div.gt.1) then
3403 block_message = 'The '//trim(i2str(n_mod_x))//&
3404 &' Fourier modes are split into '//trim(i2str(n_div))//&
3405 &' and '//trim(i2str(n_div**2))//' jobs are done separately'
3406 if (n_procs.lt.n_div**2) then
3407 block_message = trim(block_message)//', '//&
3408 &trim(i2str(n_procs))//' at a time'
3409 else
3410 block_message = trim(block_message)//', but simultaneously'
3411 end if
3412 else
3413 block_message = 'The '//trim(i2str(n_mod_x))//' Fourier modes &
3414 &can be done without splitting them'
3415 end if
3416 call writo(block_message)
3417 call writo('The total memory for all processes together is estimated &
3418 &to be about '//trim(i2str(ceiling(mem_size)))//'MB')
3419 call lvl_ud(1)
3420 call writo('(maximum: '//trim(i2str(ceiling(max_x_mem)))//&
3421 &'MB left over from equilibrium variables)')
3422 call lvl_ud(-1)
3423
3424 ! set up jobs data as illustrated below for 3 divisions, order 1:
3425 ! [1,2,3]
3426 ! or order 2:
3427 ! [1,4,7]
3428 ! [2,5,8]
3429 ! [3,6,9]
3430 ! etc.
3431 ! Also initialize the jobs taken to false.
3432 allocate(n_mod_loc(n_div))
3433 n_mod_loc = n_mod_x/n_div ! number of radial points on this processor
3434 n_mod_loc(1:mod(n_mod_x,n_div)) = n_mod_loc(1:mod(n_mod_x,n_div)) + 1 ! add a mode to if there is a remainder
3435 call calc_x_jobs_lims(x_jobs_lims,n_mod_loc,2)
3436
3437 ! user output
3438 call lvl_ud(-1)
3439 call writo('Perturbation jobs divided')
3440 contains
3441 ! Calculate X_jobs_lims.
3442 !
3443 ! These are stored as follows: The job index is the second one, while
3444 ! the first index ranges from 1..2. For every job in the range of job
3445 ! indices, thee first index shows the limits of
3446 ! the different orders sequentially. For example:
3447 ! - for order 1: [3,4],
3448 ! which corresponds to the 1-D range 3..4,
3449 ! - for order 2: [3,4,1,5],
3450 ! which corresponds to the 2-D range 3..4 x 1..5,
3451 !
3452 ! etc.
3453 !> \private
3454 recursive subroutine calc_x_jobs_lims(res,n_mod,ord)
3455 ! input / output
3456 integer, intent(inout), allocatable :: res(:,:) ! result
3457 integer, intent(in) :: n_mod(:) ! X jobs data
3458 integer, intent(in) :: ord ! order of data
3459
3460 ! local variables
3461 integer, allocatable :: res_loc(:,:) ! local result
3462 integer :: n_div ! nr. of divisions of modes
3463 integer :: id ! counter
3464
3465 ! set up nr. of divisions
3466 n_div = size(n_mod)
3467
3468 ! (re)allocate result
3469 if (allocated(res)) deallocate(res)
3470 allocate(res(2*ord,n_div**ord))
3471
3472 ! loop over divisions
3473 do id = 1,n_div
3474 if (ord.gt.1) then ! call lower order
3475 call calc_x_jobs_lims(res_loc,n_mod,ord-1)
3476 res(1:2*ord-2,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3477 &res_loc
3478 end if ! set for own order
3479 res(2*ord-1,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3480 &sum(n_mod(1:id-1))+1
3481 res(2*ord,(id-1)*n_div**(ord-1)+1:id*n_div**(ord-1)) = &
3482 &sum(n_mod(1:id))
3483 end do
3484 end subroutine calc_x_jobs_lims
3485 end function divide_x_jobs
3486
3487#if ldebug
3488 !> Prints debug information for X_1 driver
3489 integer function print_debug_x_1(mds,grid_X,X_1) result(ierr)
3490 use num_vars, only: use_pol_flux_f
3491 use grid_utilities, only: trim_grid
3492 use grid_vars, only: alpha, n_alpha
3493 use eq_vars, only: max_flux_f
3494 use rich_vars, only: rich_lvl
3495
3496 character(*), parameter :: rout_name = 'print_debug_X_1'
3497
3498 ! input / output
3499 type(modes_type), intent(in) :: mds !< general modes variables
3500 type(grid_type), intent(in) :: grid_x !< perturbation grid
3501 type(x_1_type), intent(in) :: x_1 !< vectorial X variables
3502
3503 ! local variabbles
3504 type(grid_type) :: grid_trim ! trimmed grid
3505 character(len=max_str_ln), allocatable :: var_names(:) ! names of variables
3506 character(len=max_str_ln) :: file_name ! name of file
3507 integer :: k ! local row index in flux surface
3508 integer :: ld, kd, jd, rd, id ! counters
3509 integer :: min_nm_x ! minimal n (tor. flux) or m (pol. flux)
3510 integer :: n_mod_tot ! total number of modes
3511 integer :: kdl_tot(2) ! limits on total normal index for a mode
3512 integer :: rdl ! rd in local tables
3513 integer :: kd_loc ! local kd
3514 integer :: kd_loc_trim ! local kd on trimmed grid
3515 integer :: plot_dim(4) ! dimensions of plot
3516 integer :: plot_offset(4) ! local offset of plot
3517 complex(dp), allocatable :: u(:,:,:,:,:) ! U_i for all modes
3518 complex(dp), allocatable :: du(:,:,:,:,:) ! DU_i for all modes
3519 real(dp), allocatable :: x_plot(:,:,:,:) ! X of plot
3520 real(dp), allocatable :: y_plot(:,:,:,:) ! Y of plot
3521 real(dp), allocatable :: z_plot(:,:,:,:) ! Y of plot
3522
3523 ! initiaiize ierr
3524 ierr = 0
3525
3526 ! trim grid
3527 ierr = trim_grid(grid_x,grid_trim)
3528 chckerr('')
3529
3530 ! set local n_mod and allocate integrated quantities
3531 min_nm_x = minval(mds%sec(:,1),1)
3532 n_mod_tot = size(mds%sec,1)
3533 allocate(u(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,grid_trim%n(2),&
3534 &0:1))
3535 allocate(du(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,grid_trim%n(2),&
3536 &0:1))
3537 allocate(x_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3538 &grid_trim%n(2)))
3539 allocate(y_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3540 &grid_trim%n(2)))
3541 allocate(z_plot(n_mod_tot,grid_trim%n(1),grid_trim%loc_n_r,&
3542 &grid_trim%n(2)))
3543 u = 0._dp
3544 du = 0._dp
3545
3546 ! setup X, Y and Z of plot
3547 do kd = 1,grid_trim%loc_n_r
3548 z_plot(:,:,kd,1) = 1000*grid_trim%loc_r_F(kd)/max_flux_f*2*pi
3549 do rd = 1,n_mod_tot
3550 do id = 1,grid_trim%n(1)
3551 x_plot(rd,id,kd,:) = min_nm_x+rd-1
3552 do jd = 1,grid_trim%n(2)
3553 if (use_pol_flux_f) then
3554 y_plot(rd,id,kd,jd) = grid_trim%theta_F(id,jd,kd)*10
3555 else
3556 y_plot(rd,id,kd,jd) = grid_trim%zeta_F(id,jd,kd)*10
3557 end if
3558 end do
3559 end do
3560 end do
3561 end do
3562
3563 ! select all modes combinations
3564 do k = 1,size(mds%sec,1)
3565 ! set normal limits for mode k
3566 kdl_tot(1) = mds%sec(k,2)
3567 kdl_tot(2) = mds%sec(k,3)
3568
3569 ! limit to trimmed grid range
3570 kdl_tot(1) = max(kdl_tot(1),grid_trim%i_min)
3571 kdl_tot(2) = min(kdl_tot(2),grid_trim%i_max)
3572
3573 ! skip if out of normal range
3574 if (kdl_tot(1).gt.kdl_tot(2)) cycle
3575
3576 ! total mode index (starting from 1)
3577 rd = mds%sec(k,1)-min_nm_x+1
3578
3579 ! indices in local tables
3580 rdl = mds%sec(k,4)
3581
3582 ! loop over all normal points
3583 do kd = kdl_tot(1),kdl_tot(2)
3584 ! set indices
3585 kd_loc = kd - grid_x%i_min + 1
3586 kd_loc_trim = kd - grid_trim%i_min + 1
3587
3588 u(rd,:,kd_loc_trim,:,0) = x_1%U_0(:,:,kd_loc,rdl)
3589 u(rd,:,kd_loc_trim,:,1) = x_1%U_1(:,:,kd_loc,rdl)
3590 du(rd,:,kd_loc_trim,:,0) = x_1%DU_0(:,:,kd_loc,rdl)
3591 du(rd,:,kd_loc_trim,:,1) = x_1%DU_1(:,:,kd_loc,rdl)
3592 end do
3593 end do
3594
3595 ! plot
3596 allocate(var_names(n_alpha))
3597 do ld = 1,size(var_names)
3598 var_names(ld) = 'alpha = '//trim(r2strt(alpha(ld)))
3599 end do
3600 plot_dim = [n_mod_tot,grid_trim%n(1),grid_trim%n(3),n_alpha]
3601 plot_offset = [0,0,grid_trim%i_min-1,0]
3602 do id = 0,1
3603 file_name = 'U_'//trim(i2str(id))//'_R'//&
3604 &trim(i2str(rich_lvl))
3605 call plot_hdf5(var_names,'RE_'//trim(file_name),&
3606 &rp(u(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3607 &tot_dim=plot_dim, loc_offset=plot_offset,&
3608 &col_id=4,col=1)
3609 call plot_hdf5(var_names,'IM_'//trim(file_name),&
3610 &ip(u(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3611 &tot_dim=plot_dim, loc_offset=plot_offset,&
3612 &col_id=4,col=1)
3613
3614 file_name = 'DU_'//trim(i2str(id))//'_R'//&
3615 &trim(i2str(rich_lvl))
3616 call plot_hdf5(var_names,'RE_'//trim(file_name),&
3617 &rp(du(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3618 &tot_dim=plot_dim, loc_offset=plot_offset,&
3619 &col_id=4,col=1)
3620 call plot_hdf5(var_names,'IM_'//trim(file_name),&
3621 &ip(du(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3622 &tot_dim=plot_dim, loc_offset=plot_offset,&
3623 &col_id=4,col=1)
3624 end do
3625 deallocate(var_names)
3626
3627 ! clean up
3628 call grid_trim%dealloc()
3629 end function print_debug_x_1
3630
3631 !> Prints debug information for X_2 driver
3632 integer function print_debug_x_2(mds,grid_X,X_2_int) result(ierr)
3633 use num_utilities, only: c, con
3634 use x_vars, only: n_mod_x
3635 use grid_vars, only: alpha, n_alpha
3636 use grid_utilities, only: trim_grid
3637 use eq_vars, only: max_flux_f
3638 use rich_vars, only: rich_lvl
3639
3640 character(*), parameter :: rout_name = 'print_debug_X_2'
3641
3642 ! input / output
3643 type(modes_type), intent(in) :: mds !< general modes variables
3644 type(grid_type), intent(in) :: grid_x !< perturbation grid
3645 type(x_2_type), intent(in) :: x_2_int !< tensorial X variables
3646
3647 ! local variables
3648 type(grid_type) :: grid_trim ! trimmed grid
3649 character(len=max_str_ln), allocatable :: var_names(:) ! names of variables
3650 character(len=max_str_ln) :: file_name ! name of file
3651 integer :: k, m ! local row and column index in flux surface
3652 integer :: id, ld, kd, rd, cd ! counters
3653 integer :: min_nm_x ! minimal n (tor. flux) or m (pol. flux)
3654 integer :: n_mod_tot ! total number of modes
3655 integer :: kdl_tot(2) ! limits on total normal index for a mode combination
3656 integer :: kd_loc ! local kd
3657 integer :: kd_loc_trim ! local kd on trimmed grid
3658 integer :: plot_dim(4) ! dimensions of plot
3659 integer :: plot_offset(4) ! local offset of plot
3660 integer :: c_loc(2) ! table index for symmetric and asymmetric quantities
3661 integer :: rdl, cdl ! rd and cd in local tables
3662 complex(dp), allocatable :: pv_int(:,:,:,:,:) ! integrated PV_i for all mode combinations
3663 complex(dp), allocatable :: kv_int(:,:,:,:,:) ! integrated KV_i for all mode combinations
3664 real(dp), allocatable :: x_plot(:,:,:,:) ! X of plot
3665 real(dp), allocatable :: y_plot(:,:,:,:) ! Y of plot
3666 real(dp), allocatable :: z_plot(:,:,:,:) ! Y of plot
3667
3668 ! initiaiize ierr
3669 ierr = 0
3670
3671 ! trim grid
3672 ierr = trim_grid(grid_x,grid_trim)
3673 chckerr('')
3674
3675 ! set local n_mod and allocate integrated quantities
3676 min_nm_x = minval(mds%sec(:,1),1)
3677 n_mod_tot = size(mds%sec,1)
3678 allocate(pv_int(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,&
3679 &grid_trim%n(2),0:2))
3680 allocate(kv_int(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,&
3681 &grid_trim%n(2),0:2))
3682 allocate(x_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3683 allocate(y_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3684 allocate(z_plot(n_mod_tot,n_mod_tot,grid_trim%loc_n_r,1))
3685 pv_int = 0._dp
3686 kv_int = 0._dp
3687
3688 ! setup X, Y and Z of plot
3689 do kd = 1,grid_trim%loc_n_r
3690 z_plot(:,:,kd,1) = 1000*grid_trim%loc_r_F(kd)/max_flux_f*2*pi
3691 do cd = 1,n_mod_tot
3692 do rd = 1,n_mod_tot
3693 x_plot(rd,cd,kd,1) = min_nm_x+rd-1
3694 y_plot(rd,cd,kd,1) = min_nm_x+cd-1
3695 end do
3696 end do
3697 end do
3698
3699 ! select all mode combinations
3700 do m = 1,size(mds%sec,1)
3701 do k = 1,size(mds%sec,1)
3702 ! set normal limits for mode pair (k,m)
3703 kdl_tot(1) = max(mds%sec(k,2),mds%sec(m,2))
3704 kdl_tot(2) = min(mds%sec(k,3),mds%sec(m,3))
3705
3706 ! limit to trimmed grid range
3707 kdl_tot(1) = max(kdl_tot(1),grid_trim%i_min)
3708 kdl_tot(2) = min(kdl_tot(2),grid_trim%i_max)
3709
3710 ! skip if out of normal range
3711 if (kdl_tot(1).gt.kdl_tot(2)) cycle
3712
3713 ! total mode combinations indices (starting from 1)
3714 rd = mds%sec(k,1)-min_nm_x+1
3715 cd = mds%sec(m,1)-min_nm_x+1
3716
3717 ! indices in local tables
3718 rdl = mds%sec(k,4)
3719 cdl = mds%sec(m,4)
3720
3721 ! index in tables
3722 c_loc(1) = c([rdl,cdl],.true.,n_mod_x)
3723 c_loc(2) = c([rdl,cdl],.false.,n_mod_x)
3724
3725 ! loop over all normal points
3726 do kd = kdl_tot(1),kdl_tot(2)
3727 ! set indices
3728 kd_loc = kd - grid_x%i_min + 1
3729 kd_loc_trim = kd - grid_trim%i_min + 1
3730
3731 ! symmetric quantities
3732 pv_int(rd,cd,kd_loc_trim,:,0) = &
3733 &con(x_2_int%PV_0(1,:,kd_loc,c_loc(1)),&
3734 &[rdl,cdl],.true.,[n_alpha])
3735 pv_int(rd,cd,kd_loc_trim,:,2) = &
3736 &con(x_2_int%PV_2(1,:,kd_loc,c_loc(1)),&
3737 &[rdl,cdl],.true.,[n_alpha])
3738 kv_int(rd,cd,kd_loc_trim,:,0) = &
3739 &con(x_2_int%KV_0(1,:,kd_loc,c_loc(1)),&
3740 &[rdl,cdl],.true.,[n_alpha])
3741 kv_int(rd,cd,kd_loc_trim,:,2) = &
3742 &con(x_2_int%KV_2(1,:,kd_loc,c_loc(1)),&
3743 &[rdl,cdl],.true.,[n_alpha])
3744
3745 ! asymmetric quantities
3746 pv_int(rd,cd,kd_loc_trim,:,1) = &
3747 &x_2_int%PV_1(1,:,kd_loc,c_loc(2))
3748 kv_int(rd,cd,kd_loc_trim,:,1) = &
3749 &x_2_int%KV_1(1,:,kd_loc,c_loc(2))
3750 end do
3751 end do
3752 end do
3753
3754 ! plot
3755 allocate(var_names(n_alpha))
3756 do ld = 1,size(var_names)
3757 var_names(ld) = 'alpha = '//trim(r2strt(alpha(ld)))
3758 end do
3759 plot_dim = [n_mod_tot,n_mod_tot,grid_trim%n(3),n_alpha]
3760 plot_offset = [0,0,grid_trim%i_min-1,0]
3761 do id = 0,2
3762 file_name = 'PV_'//trim(i2str(id))//'_int_R'//&
3763 &trim(i2str(rich_lvl))
3764 call plot_hdf5(var_names,'RE_'//trim(file_name),&
3765 &rp(pv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3766 &tot_dim=plot_dim, loc_offset=plot_offset,&
3767 &col_id=4,col=1)
3768 call plot_hdf5(var_names,'IM_'//trim(file_name),&
3769 &ip(pv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3770 &tot_dim=plot_dim, loc_offset=plot_offset,&
3771 &col_id=4,col=1)
3772
3773 file_name = 'KV_'//trim(i2str(id))//'_int_R'//&
3774 &trim(i2str(rich_lvl))
3775 call plot_hdf5(var_names,'RE_'//trim(file_name),&
3776 &rp(kv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3777 &tot_dim=plot_dim, loc_offset=plot_offset,&
3778 &col_id=4,col=1)
3779 call plot_hdf5(var_names,'IM_'//trim(file_name),&
3780 &ip(kv_int(:,:,:,:,id)),x=x_plot,y=y_plot,z=z_plot,&
3781 &tot_dim=plot_dim, loc_offset=plot_offset,&
3782 &col_id=4,col=1)
3783 end do
3784 deallocate(var_names)
3785
3786 ! clean up
3787 call grid_trim%dealloc()
3788 end function print_debug_x_2
3789#endif
3790end module
test_du
integer function test_du()
Definition X_ops.f90:2515
grid_utilities::calc_eqd_grid
Calculate grid of equidistant points, where optionally the last point can be excluded.
Definition grid_utilities.f90:75
grid_utilities::coord_f2e
Converts Flux coordinates to Equilibrium coordinates .
Definition grid_utilities.f90:47
hdf5_vars::dealloc_var_1d
Deallocates 1D variable.
Definition HDF5_vars.f90:68
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
num_utilities::spline
Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the m...
Definition num_utilities.f90:276
output_ops::print_ex_2d
Print 2-D output on a file.
Definition output_ops.f90:47
x_ops::calc_x
Calculates either vectorial or tensorial perturbation variables.
Definition X_ops.f90:39
x_ops::print_output_x
Print either vectorial or tensorial perturbation quantities of a certain order to an output file.
Definition X_ops.f90:85
x_ops::redistribute_output_x
Redistribute the perturbation variables.
Definition X_ops.f90:51
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition eq_vars.f90:27
eq_vars::max_flux_f
real(dp), public max_flux_f
max. flux in Flux coordinates, set in calc_norm_range_PB3D_in
Definition eq_vars.f90:50
eq_vars::vac_perm
real(dp), public vac_perm
either usual mu_0 (default) or normalized
Definition eq_vars.f90:48
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_utilities::calc_xyz_grid
integer function, public calc_xyz_grid(grid_eq, grid_xyz, x, y, z, l, r)
Calculates , and on a grid grid_XYZ, determined through its E(quilibrium) coordinates.
Definition grid_utilities.f90:799
grid_vars
Variables pertaining to the different grids used.
Definition grid_vars.f90:4
grid_vars::alpha
real(dp), dimension(:), allocatable, public alpha
field line label alpha
Definition grid_vars.f90:28
grid_vars::min_par_x
real(dp), public min_par_x
min. of parallel coordinate [ ] in field-aligned grid
Definition grid_vars.f90:24
grid_vars::n_alpha
integer, public n_alpha
nr. of field-lines
Definition grid_vars.f90:23
grid_vars::max_par_x
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
Definition grid_vars.f90:25
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition HDF5_ops.f90:27
hdf5_ops::print_hdf5_arrs
integer function, public print_hdf5_arrs(vars, pb3d_name, head_name, rich_lvl, disp_info, ind_print, remove_previous_arrs)
Prints a series of arrays, in the form of an array of pointers, to an HDF5 file.
Definition HDF5_ops.f90:1132
hdf5_vars
Variables pertaining to HDF5 and XDMF.
Definition HDF5_vars.f90:4
hdf5_vars::max_dim_var_1d
integer, parameter, public max_dim_var_1d
maximum dimension of var_1D
Definition HDF5_vars.f90:21
input_utilities
Numerical utilities related to input.
Definition input_utilities.f90:4
input_utilities::pause_prog
subroutine, public pause_prog(ind)
Pauses the running of the program.
Definition input_utilities.f90:226
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::redistribute_var
integer function, public redistribute_var(var, dis_var, lims, lims_dis)
Redistribute variables according to new limits.
Definition MPI_utilities.f90:330
num_ops
Numerical operations.
Definition num_ops.f90:4
num_ops::calc_zero_zhang
character(len=max_str_ln) function, public calc_zero_zhang(zero, fun, x_int_in)
Finds the zero of a function using Zhang's method, which is simpler than Brent's method.
Definition num_ops.f90:462
num_utilities
Numerical utilities.
Definition num_utilities.f90:4
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::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::n_theta_plot
integer, public n_theta_plot
nr. of poloidal points in plot
Definition num_vars.f90:156
num_vars::max_name_ln
integer, parameter, public max_name_ln
maximum length of filenames
Definition num_vars.f90:51
num_vars::min_theta_plot
real(dp), public min_theta_plot
min. of theta_plot
Definition num_vars.f90:159
num_vars::pi
real(dp), parameter, public pi
Definition num_vars.f90:83
num_vars::max_zeta_plot
real(dp), public max_zeta_plot
max. of zeta_plot
Definition num_vars.f90:162
num_vars::tol_norm
real(dp), public tol_norm
tolerance for normal range (normalized to 0..1)
Definition num_vars.f90:134
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition num_vars.f90:85
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
num_vars::x_jobs_lims
integer, dimension(:,:), allocatable, public x_jobs_lims
data about X jobs: [ , , , ] for all jobs
Definition num_vars.f90:76
num_vars::u_style
integer, public u_style
style for calculation of U (1: ord.2, 2: ord.1, 1: ord.0)
Definition num_vars.f90:91
num_vars::n_zeta_plot
integer, public n_zeta_plot
nr. of toroidal points in plot
Definition num_vars.f90:157
num_vars::eq_jobs_lims
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
Definition num_vars.f90:77
num_vars::x_grid_style
integer, public x_grid_style
style for normal component of X grid (1: eq, 2: sol, 3: enriched)
Definition num_vars.f90:99
num_vars::min_zeta_plot
real(dp), public min_zeta_plot
min. of zeta_plot
Definition num_vars.f90:161
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition num_vars.f90:89
num_vars::pb3d_name
character(len=max_str_ln), public pb3d_name
name of PB3D output file
Definition num_vars.f90:139
num_vars::norm_disc_prec_eq
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
Definition num_vars.f90:120
num_vars::alpha_style
integer, public alpha_style
style for alpha (1: one field line, many turns, 2: many field lines, one turn)
Definition num_vars.f90:100
num_vars::rank
integer, public rank
MPI rank.
Definition num_vars.f90:68
num_vars::k_style
integer, public k_style
style for kinetic energy
Definition num_vars.f90:93
num_vars::norm_disc_prec_x
integer, public norm_disc_prec_x
precision for normal discretization for perturbation
Definition num_vars.f90:121
num_vars::x_style
integer, public x_style
style for secondary mode numbers (1: prescribed, 2: fast)
Definition num_vars.f90:95
num_vars::max_theta_plot
real(dp), public max_theta_plot
max. of theta_plot
Definition num_vars.f90:160
num_vars::eq_job_nr
integer, public eq_job_nr
nr. of eq job
Definition num_vars.f90:79
num_vars::use_pol_flux_f
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
Definition num_vars.f90:114
num_vars::magn_int_style
integer, public magn_int_style
style for magnetic integrals (1: trapezoidal, 2: Simpson 3/8)
Definition num_vars.f90:124
num_vars::max_x_mem
real(dp), public max_x_mem
maximum memory for perturbation calculations for all processes [MB]
Definition num_vars.f90:75
num_vars::no_plots
logical, public no_plots
no plots made
Definition num_vars.f90:140
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
output_ops::draw_ex
subroutine, public draw_ex(var_names, draw_name, nplt, draw_dim, plot_on_screen, ex_plot_style, data_name, draw_ops, extra_ops, is_animated, ranges, delay, persistent)
Use external program to draw a plot.
Definition output_ops.f90:1079
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
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::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition str_utilities.f90:42
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition str_utilities.f90:54
vmec_utilities
Numerical utilities related to the output of VMEC.
Definition VMEC_utilities.f90:4
vmec_utilities::calc_trigon_factors
integer function, public calc_trigon_factors(theta, zeta, trigon_factors)
Calculate the trigonometric cosine and sine factors.
Definition VMEC_utilities.f90:275
x_ops
Operations considering perturbation quantities.
Definition X_ops.f90:4
x_ops::calc_res_surf
integer function, public calc_res_surf(mds, grid_eq, eq, res_surf, info, jq)
Calculates resonating flux surfaces for the perturbation modes.
Definition X_ops.f90:1638
x_ops::debug_check_x_modes_2
logical, public debug_check_x_modes_2
plot debug information for check_x_modes_2()
Definition X_ops.f90:26
x_ops::setup_modes
integer function, public setup_modes(mds, grid_eq, grid, plot_name)
Sets up some variables concerning the mode numbers.
Definition X_ops.f90:1060
x_ops::divide_x_jobs
integer function, public divide_x_jobs(arr_size)
Divides the perturbation jobs.
Definition X_ops.f90:3364
x_ops::init_modes
integer function, public init_modes(grid_eq, eq)
Initializes some variables concerning the mode numbers.
Definition X_ops.f90:919
x_ops::print_debug_x_1
integer function, public print_debug_x_1(mds, grid_x, x_1)
Prints debug information for X_1 driver.
Definition X_ops.f90:3490
x_ops::check_x_modes
integer function, public check_x_modes(grid_eq, eq)
Checks whether the high-n approximation is valid:
Definition X_ops.f90:1350
x_ops::resonance_plot
integer function, public resonance_plot(mds, grid_eq, eq)
plot -profile or -profile in flux coordinates with or indicated if requested.
Definition X_ops.f90:1814
x_ops::print_debug_x_2
integer function, public print_debug_x_2(mds, grid_x, x_2_int)
Prints debug information for X_2 driver.
Definition X_ops.f90:3633
x_ops::calc_magn_ints
integer function, public calc_magn_ints(grid_eq, grid_x, eq, x, x_int, prev_style, lim_sec_x)
Calculate the magnetic integrals from and in an equidistant grid where the step size can vary depen...
Definition X_ops.f90:3081
x_utilities
Numerical utilities related to perturbation operations.
Definition X_utilities.f90:4
x_utilities::calc_memory_x
integer function, public calc_memory_x(ord, arr_size, n_mod, mem_size)
Calculate memory in MB necessary for X variables.
Definition X_utilities.f90:236
x_utilities::get_sec_x_range
subroutine, public get_sec_x_range(sec_x_range_loc, sec_x_range_tot, m, sym, lim_sec_x)
Gets one of the the local ranges of contiguous tensorial perturbation variables to be printed or read...
Definition X_utilities.f90:128
x_utilities::is_necessary_x
logical function, public is_necessary_x(sym, sec_x_id, lim_sec_x)
Determines whether a variable needs to be considered:
Definition X_utilities.f90:197
x_vars
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
x_vars::min_n_x
integer, dimension(:), allocatable, public min_n_x
lowest poloidal mode number m_X, in total eq grid
Definition X_vars.f90:131
x_vars::min_sec_x
integer, public min_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for X style 1)
Definition X_vars.f90:127
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
x_vars::max_n_x
integer, dimension(:), allocatable, public max_n_x
highest poloidal mode number m_X, in total eq grid
Definition X_vars.f90:132
x_vars::set_nn_mod
integer function, public set_nn_mod(sym, lim_sec_x)
Sets number of entries for tensorial perturbation variables.
Definition X_vars.f90:383
x_vars::max_m_x
integer, dimension(:), allocatable, public max_m_x
highest poloidal mode number m_X, in total eq grid
Definition X_vars.f90:134
x_vars::min_m_x
integer, dimension(:), allocatable, public min_m_x
lowest poloidal mode number m_X, in total eq grid
Definition X_vars.f90:133
x_vars::min_nm_x
integer, public min_nm_x
minimum for the high-n theory (debable)
Definition X_vars.f90:130
x_vars::max_sec_x
integer, public max_sec_x
m_X (pol. flux) or n_X (tor. flux) (only for\ c X style 1)
Definition X_vars.f90:128
x_vars::prim_x
integer, public prim_x
n_X (pol. flux) or m_X (tor. flux)
Definition X_vars.f90:126
eq_vars::eq_1_type
flux equilibrium type
Definition eq_vars.f90:63
eq_vars::eq_2_type
metric equilibrium type
Definition eq_vars.f90:114
grid_vars::grid_type
Type for grids.
Definition grid_vars.f90:59
hdf5_vars::var_1d_type
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition HDF5_vars.f90:48
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