PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
grid_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations that have to do with the grids and different coordinate systems.
3!------------------------------------------------------------------------------!
4module grid_ops
5#include <PB3D_macros.h>
6 use str_utilities
7 use output_ops
8 use messages
9 use num_vars, only: dp, pi, max_str_ln
10 use grid_vars, only: grid_type
11 use eq_vars, only: eq_1_type
12
13 implicit none
14 private
15 public calc_norm_range, calc_ang_grid_eq_b, magn_grid_plot, setup_grid_eq, &
16 &setup_grid_eq_b, setup_grid_x, setup_grid_sol, print_output_grid, &
17 &redistribute_output_grid
18#if ldebug
19 public debug_calc_ang_grid_eq_b
20#endif
21
22 ! global variables
23#if ldebug
24 !> \ldebug
25 logical :: debug_calc_ang_grid_eq_b = .false. !< plot debug information for calc_ang_grid_eq_b()
26#endif
27
28contains
29 !> Calculates normal range for the input grid, the equilibrium grid and/or
30 !! the solution grid.
31 !!
32 !! General workings, depending on \c X_grid_style:
33 !!
34 !! | 1 (equilibrium) | 2 (solution ) | 3 (enriched) |
35 !! | ------------------- | ------------------- | ---------------------- |
36 !! | calc_eq() | calc_eq() | calc_eq() |
37 !! | | redistribute to sol | add points to optimize |
38 !! | | interpolate to sol | interpolate to X |
39 !! | calc_x() | calc_x() | calc_x() |
40 !! | redistribute to sol | copy to sol | redistribute to sol |
41 !! | interpolate to sol | | interpolate to sol |
42 !!
43 !! \return ierr
44 integer function calc_norm_range(style,in_limits,eq_limits,X_limits,&
45 &sol_limits,r_F_eq,r_F_X,r_F_sol,jq) result(ierr)
46
47 character(*), parameter :: rout_name = 'calc_norm_range'
48
49 ! input / output
50 character(len=*), intent(in) :: style !< style of calculation (PB3D: in, eq, X or sol; POST)
51 integer, intent(inout), optional :: in_limits(2) !< min. and max. index of in grid
52 integer, intent(inout), optional :: eq_limits(2) !< min. and max. index of eq grid for this process
53 integer, intent(inout), optional :: x_limits(2) !< min. and max. index of X grid for this process
54 integer, intent(inout), optional :: sol_limits(2) !< min. and max. index of sol grid for this process
55 real(dp), intent(inout), optional :: r_f_eq(:) !< equilibrium r_F
56 real(dp), intent(inout), allocatable, optional :: r_f_x(:) !< perturbation r_F
57 real(dp), intent(inout), optional :: r_f_sol(:) !< solution r_F
58 real(dp), intent(in), optional :: jq(:) !< q_saf (pol. flux) or rot_t (tor. flux) in total Flux variables
59
60 ! local variables
61 character(len=max_str_ln) :: err_msg ! error message
62
63 ! initialize ierr
64 ierr = 0
65
66 ! select depending on style
67 select case (style)
68 case ('PB3D_in') ! PB3D: in
69 if (present(in_limits)) then
70 ierr = calc_norm_range_pb3d_in(in_limits)
71 chckerr('')
72 else
73 ierr = 1
74 err_msg = 'for PB3D: in, in_limits need to be provided'
75 chckerr(err_msg)
76 end if
77 case ('PB3D_eq') ! PB3D: eq
78 if (present(eq_limits)) then
79 call calc_norm_range_pb3d_eq(eq_limits)
80 else
81 ierr = 1
82 err_msg = 'for PB3D: eq, eq_limits need to be provided'
83 chckerr(err_msg)
84 end if
85 case ('PB3D_X') ! PB3D: X
86 if (present(eq_limits).and.present(x_limits).and.&
87 &present(r_f_eq).and.present(r_f_x).and.(present(jq))) then
88 ierr = calc_norm_range_pb3d_x(eq_limits,x_limits,&
89 &r_f_eq,r_f_x,jq)
90 chckerr('')
91 else
92 ierr = 1
93 err_msg = 'for PB3D: X, eq_limits, X_limits, r_F_eq, &
94 &r_F_X and jq need to be provided'
95 chckerr(err_msg)
96 end if
97 case ('PB3D_sol') ! PB3D: sol
98 if (present(sol_limits).and.present(r_f_sol)) then
99 ierr = calc_norm_range_pb3d_sol(sol_limits,r_f_sol)
100 chckerr('')
101 else
102 ierr = 1
103 err_msg = 'for PB3D: sol, sol_limits and r_F_sol need to &
104 &be provided'
105 chckerr(err_msg)
106 end if
107 case ('POST') ! POST
108 if (present(eq_limits) .and. present(x_limits) .and. &
109 &present(sol_limits) .and. present(r_f_eq) .and. &
110 &present(r_f_x) .and. present(r_f_sol)) then
111 call calc_norm_range_post(eq_limits,x_limits,sol_limits,&
112 &r_f_eq,r_f_x,r_f_sol)
113 else
114 ierr = 1
115 err_msg = 'for POST, eq_limits, X_limits, sol_limits, &
116 &r_F_eq, r_F_X and r_F_sol need to be provided'
117 chckerr(err_msg)
118 end if
119 case default
120 ierr = 1
121 err_msg = 'Incorrect style "'//trim(style)//'"'
122 chckerr(err_msg)
123 end select
124 contains
125 !> \public PB3D input version
126 !!
127 !! The normal range is calculated by finding the tightest range of the
128 !! input variables encompassing the entire solution range. Additionally,
129 !! the global max_flux variables are set as well..
130 !!
131 !! \note This is the global, undivided range. The division information
132 !! is calculated in 'eq_limits'. Furthermore, the input variables have
133 !! to be tabulated on the full grid provided by the equilibrium code to
134 !! be able to be used in PB3D and POST.
135 integer function calc_norm_range_pb3d_in(in_limits) &
136 &result(ierr) ! PB3D version for equilibrium grid
137 use num_vars, only: use_pol_flux_e, use_pol_flux_f, eq_style, &
138 &norm_disc_prec_eq
139 use num_utilities, only: con2dis, dis2con, round_with_tol, spline
140 use vmec_vars, only: flux_p_v, flux_t_v
141 use helena_vars, only: flux_p_h, flux_t_h
142 use x_vars, only: min_r_sol, max_r_sol
143 use eq_vars, only: max_flux_e, max_flux_f
144 use grid_vars, only: n_r_in
145
146 character(*), parameter :: rout_name = 'calc_norm_range_PB3D_in'
147
148 ! input / output
149 integer, intent(inout) :: in_limits(2) !< total min. and max. index of eq. grid for this process
150
151 ! local variables
152 real(dp), allocatable :: flux_f(:), flux_e(:) ! either pol. or tor. flux in F and E
153 real(dp) :: tot_min_r_in_f_con ! tot_min_r_in in continuous F coords.
154 real(dp) :: tot_max_r_in_f_con ! tot_max_r_in in continuous F coords.
155 real(dp) :: tot_min_r_in_e_con(1) ! tot_min_r_in in continuous E coords.
156 real(dp) :: tot_max_r_in_e_con(1) ! tot_max_r_in in continuous E coords.
157 real(dp) :: tot_min_r_in_e_dis ! tot_min_r_in in discrete E coords., unrounded
158 real(dp) :: tot_max_r_in_e_dis ! tot_max_r_in in discrete E coords., unrounded
159
160 ! initialize ierr
161 ierr = 0
162
163 ! set up flux_F and flux_E, depending on which equilibrium style is
164 ! being used:
165 ! 1: VMEC
166 ! 2: HELENA
167 allocate(flux_f(n_r_in),flux_e(n_r_in))
168 select case (eq_style)
169 case (1) ! VMEC
170 ! set up F flux
171 if (use_pol_flux_f) then
172 flux_f = flux_p_v(:,0)
173 else
174 flux_f = - flux_t_v(:,0) ! conversion VMEC LH -> RH coord. system
175 end if
176 ! set up E flux
177 if (use_pol_flux_e) then
178 flux_f = flux_p_v(:,0)
179 else
180 flux_e = flux_t_v(:,0)
181 end if
182 case (2) ! HELENA
183 ! set up F flux
184 if (use_pol_flux_f) then
185 flux_f = flux_p_h(:,0)
186 else
187 flux_f = flux_t_h(:,0)
188 end if
189 ! set up E flux
190 if (use_pol_flux_e) then
191 flux_e = flux_p_h(:,0)
192 else
193 flux_e = flux_t_h(:,0)
194 end if
195 end select
196
197 ! set max_flux
198 max_flux_e = flux_e(n_r_in)
199 max_flux_f = flux_f(n_r_in)
200
201 ! normalize flux_F and flux_E to (0..1) using max_flux
202 flux_e = flux_e/max_flux_e
203 flux_f = flux_f/max_flux_f
204
205 ! get lower and upper bound of total solution range
206 ! 1. start
207 tot_min_r_in_f_con = min_r_sol
208 tot_max_r_in_f_con = max_r_sol
209 ! 2. round with tolerance
210 ierr = round_with_tol(tot_min_r_in_f_con,0._dp,1._dp)
211 chckerr('')
212 ierr = round_with_tol(tot_max_r_in_f_con,0._dp,1._dp)
213 chckerr('')
214 ! 3. continuous E coords
215 ierr = spline(flux_f,flux_e,[tot_min_r_in_f_con],&
216 &tot_min_r_in_e_con,ord=norm_disc_prec_eq)
217 ierr = spline(flux_f,flux_e,[tot_max_r_in_f_con],&
218 &tot_max_r_in_e_con,ord=norm_disc_prec_eq)
219 ! 4. round with tolerance
220 ierr = round_with_tol(tot_min_r_in_e_con,minval(flux_e),&
221 &maxval(flux_e))
222 chckerr('')
223 ierr = round_with_tol(tot_max_r_in_e_con,minval(flux_e),&
224 &maxval(flux_e))
225 chckerr('')
226 ! 5. discrete E index, unrounded
227 ierr = con2dis(tot_min_r_in_e_con(1),tot_min_r_in_e_dis,flux_e)
228 chckerr('')
229 ierr = con2dis(tot_max_r_in_e_con(1),tot_max_r_in_e_dis,flux_e)
230 chckerr('')
231 ! 6. discrete E index, rounded
232 in_limits(1) = floor(tot_min_r_in_e_dis)
233 in_limits(2) = ceiling(tot_max_r_in_e_dis)
234 end function calc_norm_range_pb3d_in
235
236 !> \public PB3D equilibrium version
237 !!
238 !! The normal range is calculated by dividing the full equilibrium range
239 !! passed from the input phase between the processes.
240 !!
241 !! \note For X_style 2 (solution) or 3 (optimized), at the end of the
242 !! equilibrium phase, there will be an exchange of data from each
243 !! process to each process so that they all have the tightest possible
244 !! fit of data of the corresponding perturbation limits.
245 subroutine calc_norm_range_pb3d_eq(eq_limits) ! PB3D version for equilibrium grid
246 use num_vars, only: n_procs, rank, norm_disc_prec_eq
247 use grid_vars, only: n_r_eq
248
249 ! input / output
250 integer, intent(inout) :: eq_limits(2) !< min. and max. index of eq. grid for this process
251
252 ! set local equilibrium limits
253 eq_limits(1) = nint(1 + 1._dp*rank*(n_r_eq-1)/n_procs)
254 eq_limits(2) = nint(1 + (1._dp+rank)*(n_r_eq-1)/n_procs)
255
256 ! increase lower limits for processes that are not first
257 if (rank.gt.0) eq_limits(1) = eq_limits(1)+1
258
259 ! ghost regions of width 4*norm_disc_prec_eq
260 ! Note: When finite differences were used, a factor of 2 was enough,
261 ! but it looks like for splines it has to be higher. Even for 4,
262 ! there is some slight discrepancy... This should be looked into in
263 ! the future.
264 eq_limits(1) = max(eq_limits(1)-4*norm_disc_prec_eq,1)
265 eq_limits(2) = min(eq_limits(2)+4*norm_disc_prec_eq,n_r_eq)
266 end subroutine calc_norm_range_pb3d_eq
267
268 !> \public PB3D perturbation version
269 !!
270 !! The normal range is determined according to X_grid style:
271 !! 1. taken identical to the equilibrium grid, which means that after
272 !! the perturbation phase the integrated tensorial quantities have to
273 !! be interpolated in interp_v() in driver_sol() after having been
274 !! redistributed at the start of the solution driver.
275 !! 2. taken identical to the solution grid, which means that the
276 !! the equilibrium quantities have to be interpolated in calc_u(),
277 !! calc_kv() and calc_pv() after having been redistributed at the end
278 !! of the equilibrium driver.
279 !! 3. by considering an optimal interpolation extension of the
280 !! equilibrium range, adding intermediairy points between grid points
281 !! where the safety factor changes too quickly, according to the
282 !! variable 'max_njq_change'. This uses \c prim_X.
283 integer function calc_norm_range_pb3d_x(eq_limits,X_limits,r_F_eq,&
284 &r_F_X,jq) result(ierr) ! PB3D version for perturbation grid
285
286 use num_vars, only: x_grid_style, max_njq_change, rank
287 use grid_vars, only: n_r_eq, n_r_sol
288 use x_vars, only: prim_x
289 use grid_utilities, only: calc_eqd_grid
290 use eq_vars, only: max_flux_f
291
292 character(*), parameter :: rout_name = 'calc_norm_range_PB3D_X'
293
294 ! input / output
295 integer, intent(in) :: eq_limits(2) !< min. and max. index of eq grid for this process
296 integer, intent(inout) :: x_limits(2) !< min. and max. index of X grid for this process
297 real(dp), intent(in) :: r_f_eq(:) !< equilibrium r_F
298 real(dp), intent(inout), allocatable :: r_f_x(:) !< perturbation r_F
299 real(dp), intent(in) :: jq(:) !< q_saf (pol. flux) or rot_t (tor. flux) in total Flux variables
300
301 ! local variables
302 integer :: kd ! counter
303 integer, allocatable :: div(:) ! number of extra divisions for this grid interval
304 real(dp), allocatable :: r_f_plot(:,:) ! plot of r_F_eq and r_F_X
305
306 ! initialize ierr
307 ierr = 0
308
309 select case (x_grid_style)
310 case (1) ! equilibrium
311 ! copy equilibrium
312 allocate(r_f_x(n_r_eq))
313 r_f_x = r_f_eq
314 x_limits = eq_limits
315
316 ! output
317 call writo('Using the equilibrium grid')
318 case (2) ! solution
319 ! calculate solution
320 allocate(r_f_x(n_r_sol))
321 ierr = calc_norm_range_pb3d_sol(sol_limits=x_limits,&
322 &r_f_sol=r_f_x)
323 chckerr('')
324
325 ! output
326 call writo('Using the solution grid')
327 case (3) ! enriched
328 ! decide extra divisions for each interval
329 allocate(div(size(r_f_eq)-1))
330 do kd = 1,size(r_f_eq)-1
331 div(kd) = &
332 &floor(prim_x*abs(jq(kd+1)-jq(kd))/max_njq_change)
333 end do
334
335 ! set up r_F_X with divisions
336 allocate(r_f_x(size(r_f_eq)+sum(div)))
337 do kd = 1,size(r_f_eq)-1
338 ierr = calc_eqd_grid(&
339 &r_f_x(kd+sum(div(1:kd-1)):kd+sum(div(1:kd))),&
340 &r_f_eq(kd),r_f_eq(kd+1),excl_last=.true.)
341 chckerr('')
342 end do
343 r_f_x(size(r_f_x)) = r_f_eq(size(r_f_eq))
344
345 ! set up normal limits
346 x_limits(1) = eq_limits(1)+sum(div(1:eq_limits(1)-1))
347 x_limits(2) = eq_limits(2)+sum(div(1:eq_limits(2)-1))
348
349 ! plot
350 if (rank.eq.0) then
351 allocate(r_f_plot(size(r_f_x),4))
352 r_f_plot(1:size(r_f_eq),1) = r_f_eq
353 r_f_plot(1:size(r_f_eq),3) = [(kd,kd=1,size(r_f_eq))]
354 r_f_plot(size(r_f_eq)+1:size(r_f_x),1) = &
355 &r_f_eq(size(r_f_eq))
356 r_f_plot(size(r_f_eq)+1:size(r_f_x),3) = size(r_f_eq)
357 r_f_plot(:,2) = r_f_x
358 r_f_plot(:,4) = [(kd,kd=1,size(r_f_x))]
359 r_f_plot(:,1:2) = r_f_plot(:,1:2)/max_flux_f*2*pi
360 call print_ex_2d(['eq','X '],'r_F_X',&
361 &r_f_plot(:,3:4),x=r_f_plot(:,1:2),draw=.false.)
362 call draw_ex(['eq','X '],'r_F_X',2,1,&
363 &.false.)
364 end if
365
366 ! output
367 call writo('Enriching the equilibrium grid to have jq &
368 &change not more than '//trim(r2strt(max_njq_change)))
369 call lvl_ud(1)
370 if (any(div.gt.0)) then
371 call writo('Added '//trim(i2str(sum(div)))//' points')
372 else
373 call writo('But it was not necessary to add points')
374 end if
375 call lvl_ud(-1)
376 end select
377 end function calc_norm_range_pb3d_x
378
379 !> \public PB3D solution version
380 !!
381 !! The normal range is determined by simply dividing the solution range,
382 !! a ghost range is required.
383 !!
384 !! By default, this routine uses \c sol_n_procs processes, but this can
385 !! be overruled.
386 integer function calc_norm_range_pb3d_sol(sol_limits,r_F_sol,n_procs) &
387 &result(ierr) ! PB3D version for solution grid
388 use num_vars, only: sol_n_procs, rank, norm_disc_prec_sol
389 use eq_vars, only: max_flux_f
390 use x_vars, only: min_r_sol, max_r_sol
391 use num_utilities, only: round_with_tol
392 use grid_vars, only: n_r_sol
393
394 character(*), parameter :: rout_name = 'calc_norm_range_PB3D_sol'
395
396 ! input / output
397 integer, intent(inout) :: sol_limits(2) !< min. and max. index of sol grid for this process
398 real(dp), intent(inout) :: r_f_sol(:) !< solution r_F
399 integer, intent(in), optional :: n_procs !< how many processes used
400
401 ! local variables
402 integer :: id ! counter
403 integer, allocatable :: loc_n_r_sol(:) ! local nr. of normal points in solution grid
404 integer :: n_procs_loc ! local n_procs
405
406 ! initialize ierr
407 ierr = 0
408
409 ! set local n_procs
410 n_procs_loc = sol_n_procs
411 if (present(n_procs)) n_procs_loc = n_procs
412
413 ! set loc_n_r_sol
414 allocate(loc_n_r_sol(n_procs_loc))
415 loc_n_r_sol = n_r_sol/n_procs_loc ! number of radial points on this processor
416 loc_n_r_sol(1:mod(n_r_sol,n_procs_loc)) = &
417 &loc_n_r_sol(1:mod(n_r_sol,n_procs_loc)) + 1 ! add a mode to if there is a remainder
418
419 ! set sol_limits
420 if (rank.lt.n_procs_loc) then
421 sol_limits = &
422 &[sum(loc_n_r_sol(1:rank))+1,sum(loc_n_r_sol(1:rank+1))]
423
424 ! ghost regions of width 2*norm_disc_prec_sol
425 sol_limits(1) = max(sol_limits(1)-norm_disc_prec_sol,1)
426 sol_limits(2) = min(sol_limits(2)+norm_disc_prec_sol,n_r_sol)
427 else
428 sol_limits(1) = n_r_sol
429 sol_limits(2) = n_r_sol
430 end if
431
432 ! calculate r_F_sol in range from min_r_sol to max_r_sol
433 r_f_sol = [(min_r_sol + (id-1.0_dp)/(n_r_sol-1.0_dp)*&
434 &(max_r_sol-min_r_sol),id=1,n_r_sol)]
435
436 ! round with standard tolerance
437 ierr = round_with_tol(r_f_sol,0.0_dp,1.0_dp)
438 chckerr('')
439
440 ! translate to the real normal variable in range from
441 ! 0..flux/2pi
442 r_f_sol = r_f_sol*max_flux_f/(2*pi)
443 end function calc_norm_range_pb3d_sol
444
445 !> \public POST version
446 !!
447 !! The normal range is determined by simply dividing a possible subset
448 !! of the solution range, indicated by \c min_r_plot and \c max_r_sol,
449 !! including a ghost range and getting a bounding equilibrium and
450 !! perturbation range.
451 subroutine calc_norm_range_post(eq_limits,X_limits,sol_limits,r_F_eq,&
452 &r_F_X,r_F_sol) ! POST version
453
454 use num_vars, only: n_procs, rank, norm_disc_prec_sol, &
455 &min_r_plot, max_r_plot, x_grid_style
456 use eq_vars, only: max_flux_f
457 use grid_utilities, only: find_compr_range
458
459 ! input / output
460 integer, intent(inout) :: eq_limits(2) !< min. and max. index of eq grid for this process
461 integer, intent(inout) :: X_limits(2) !< min. and max. index of X grid for this process
462 integer, intent(inout) :: sol_limits(2) !< min. and max. index of sol grid for this process
463 real(dp), intent(in) :: r_F_eq(:) !< eq r_F
464 real(dp), intent(in) :: r_F_X(:) !< X r_F
465 real(dp), intent(in) :: r_F_sol(:) !< sol r_F
466
467 ! local variables
468 integer :: sol_limits_tot(2) ! total solution limits
469 integer :: n_r_sol ! total nr. of normal points in solution grid
470 integer, allocatable :: loc_n_r_sol(:) ! local nr. of normal points in solution grid
471 real(dp) :: min_sol, max_sol ! min. and max. of r_F_sol in range of this process
472
473 ! initialize ierr
474 ierr = 0
475
476 ! get min and max of solution range
477 min_sol = max(minval(r_f_sol),min_r_plot*max_flux_f/(2*pi))
478 max_sol = min(maxval(r_f_sol),max_r_plot*max_flux_f/(2*pi))
479
480 ! find the solution index that comprises this range
481 call find_compr_range(r_f_sol,[min_sol,max_sol],sol_limits_tot)
482
483 ! initialize n_r_sol
484 n_r_sol = sol_limits_tot(2)-sol_limits_tot(1)+1
485 allocate(loc_n_r_sol(n_procs))
486
487 ! divide the solution grid equally over all the processes
488 loc_n_r_sol = n_r_sol/n_procs ! number of radial points on this processor
489 loc_n_r_sol(1:mod(n_r_sol,n_procs)) = &
490 &loc_n_r_sol(1:mod(n_r_sol,n_procs)) + 1 ! add a mode to if there is a remainder
491
492 ! set sol_limits
493 sol_limits = [sum(loc_n_r_sol(1:rank))+1,sum(loc_n_r_sol(1:rank+1))]
494 sol_limits = sol_limits + sol_limits_tot(1) - 1
495 if (rank.gt.0) sol_limits(1) = sol_limits(1)-norm_disc_prec_sol ! ghost region for num. deriv.
496 if (rank.lt.n_procs-1) sol_limits(2) = &
497 &sol_limits(2)+norm_disc_prec_sol+1 ! ghost region for num. deriv. and overlap for int_vol
498 min_sol = minval(r_f_sol(sol_limits(1):sol_limits(2)))
499 max_sol = maxval(r_f_sol(sol_limits(1):sol_limits(2)))
500
501 ! determine eq_limits: smallest eq and X range comprising sol range
502 call find_compr_range(r_f_eq,[min_sol,max_sol],eq_limits)
503 select case (x_grid_style)
504 case (1,3) ! equilibrium or enriched
505 call find_compr_range(r_f_x,[min_sol,max_sol],x_limits)
506 case (2) ! solution
507 x_limits = sol_limits
508 end select
509 end subroutine calc_norm_range_post
510 end function calc_norm_range
511
512 !> Sets up the equilibrium grid.
513 !!
514 !! The following variables are calculated:
515 !! - equilibrium variables (eq)
516 !! - perturbation variables (X)
517 !!
518 !! For the implementation of the equilibrium grid the normal part of the
519 !! grid is always given by the output of the equilibrium code, but the
520 !! angular part depends on the style:
521 !! - VMEC: The output is analytic in the angular coordinates, so
522 !! field-aligned coordinates are used. Later, in calc_ang_grid_eq_b(), the
523 !! angular variables are calculated.
524 !! - HELENA: The output is given on a poloidal grid (no toroidal dependency
525 !! due to axisymmetry), which is used.
526 !!
527 !! \return ierr
528 integer function setup_grid_eq(grid_eq,eq_limits) result(ierr)
529 use num_vars, only: eq_style, eq_jobs_lims, eq_job_nr
530 use grid_vars, only: n_r_eq, n_alpha
531 use helena_vars, only: nchi, chi_h
532 use rich_vars, only: n_par_x
533
534 character(*), parameter :: rout_name = 'setup_grid_eq'
535
536 ! input / output
537 type(grid_type), intent(inout) :: grid_eq !< equilibrium grid
538 integer, intent(in) :: eq_limits(2) !< min. and max. index of eq grid of this process
539
540 ! local variables
541 integer :: id ! counter
542
543 ! initialize ierr
544 ierr = 0
545
546 ! set up general equilibrium grid:
547 ! choose which equilibrium style is being used:
548 ! 1: VMEC
549 ! 2: HELENA
550 select case (eq_style)
551 case (1) ! VMEC
552 ! user output
553 call writo('Field-aligned with '//trim(i2str(n_r_eq))//&
554 &' normal and '//trim(i2str(n_par_x))//' parallel points')
555
556 ! create grid with eq_jobs_lims
557 ierr = grid_eq%init([eq_jobs_lims(2,eq_job_nr)-&
558 &eq_jobs_lims(1,eq_job_nr)+1,n_alpha,n_r_eq],eq_limits)
559 chckerr('')
560 case (2) ! HELENA
561 ! user output
562 call writo('Identical to the equilibrium input grid')
563 call writo(trim(i2str(n_r_eq))//' normal and '//&
564 &trim(i2str(nchi))//' angular points')
565 call writo('Poloidal range '//trim(r2strt(chi_h(1)))//'..'//&
566 &trim(r2strt(chi_h(nchi))))
567 call writo('Will be interpolated to a field-aligned grid &
568 &later')
569
570 ! create grid
571 ierr = grid_eq%init([nchi,1,n_r_eq],eq_limits) ! axisymmetric equilibrium
572 chckerr('')
573
574 ! copy angular grid from HELENA
575 do id = 1,grid_eq%n(1)
576 grid_eq%theta_E(id,:,:) = chi_h(id)
577 end do
578 grid_eq%zeta_E = 0._dp
579
580 ! convert to Flux coordinates (trivial)
581 grid_eq%theta_F = grid_eq%theta_E
582 grid_eq%zeta_F = grid_eq%zeta_E
583 end select
584 end function setup_grid_eq
585
586 !> Sets up the field-aligned equilibrium grid.
587 !!
588 !! This serves as a bridge to the solution grid, as it contains the same
589 !! normal coordinate as the general grid, but the angular coordinates are
590 !! defined by the solution grid.
591 !!
592 !! Optionally, only half the grid can be calculated (i.e. only the even
593 !! points), which is used for Richardson levels greater than 1.
594 !!
595 !! \note In contrast to \c setup_grid_eq, the angular coordinates are also
596 !! calculated here.
597 !!
598 !! \return ierr
599 integer function setup_grid_eq_b(grid_eq,grid_eq_B,eq,only_half_grid) &
600 &result(ierr)
601 use num_vars, only: eq_style, eq_jobs_lims, eq_job_nr
602 use grid_vars, only: n_alpha
603
604 character(*), parameter :: rout_name = 'setup_grid_eq_B'
605
606 ! input / output
607 type(grid_type), intent(inout) :: grid_eq !< general equilibrium grid
608 type(grid_type), intent(inout) :: grid_eq_b !< field-aligned equilibrium grid
609 type(eq_1_type), intent(in) :: eq !< flux equilibrium variables
610 logical, intent(in), optional :: only_half_grid !< calculate only half grid with even points
611
612 ! local variables
613 character(len=max_str_ln) :: err_msg ! error message
614
615 ! initialize ierr
616 ierr = 0
617
618 ! choose which equilibrium style is being used:
619 ! 1: VMEC
620 ! 2: HELENA
621 select case (eq_style)
622 case (1) ! VMEC
623 ! the grid is already field aligned
624 ierr = 1
625 err_msg = 'The grid is already field-aligned for VMEC'
626 chckerr(err_msg)
627 case (2) ! HELENA
628 ! create grid with eq_jobs_lims
629 ierr = grid_eq_b%init([eq_jobs_lims(2,eq_job_nr)-&
630 &eq_jobs_lims(1,eq_job_nr)+1,n_alpha,grid_eq%n(3)],&
631 &[grid_eq%i_min,grid_eq%i_max]) ! only one field line
632 chckerr('')
633
634 ! copy the normal coords.
635 grid_eq_b%loc_r_E = grid_eq%loc_r_E
636 grid_eq_b%loc_r_F = grid_eq%loc_r_F
637 grid_eq_b%r_E = grid_eq%r_E
638 grid_eq_b%r_F = grid_eq%r_F
639
640 ! calculate the angular grid that follows the magnetic field
641 ierr = calc_ang_grid_eq_b(grid_eq_b,eq,only_half_grid)
642 chckerr('')
643 end select
644 end function setup_grid_eq_b
645
646 !> Sets up the general perturbation grid, in which the perturbation
647 !! variables are calculated.
648 !!
649 !! For \c X_grid_style 1, this grid is identical to the equilibrium grid,
650 !! and for \c X_grid_style 2,it has the same angular extent but with
651 !! different normal points, indicated by the variable \c r_F_X.
652 !!
653 !! \return ierr
654 integer function setup_grid_x(grid_eq,grid_X,r_F_X,X_limits) result(ierr)
655 use num_vars, only: norm_disc_prec_x, x_grid_style
656 use grid_utilities, only: coord_f2e
657 use num_utilities, only: spline
658
659 character(*), parameter :: rout_name = 'setup_grid_X'
660
661 ! input / output
662 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
663 type(grid_type), intent(inout) :: grid_x !< perturbation grid
664 real(dp), intent(in), allocatable :: r_f_x(:) !< points of perturbation grid
665 integer, intent(in) :: x_limits(2) !< min. and max. index of perturbation grid of this process
666
667 ! local variables
668 integer :: id, jd ! counters
669
670 ! initialize ierr
671 ierr = 0
672
673 select case (x_grid_style)
674 case (1) ! equilibrium
675 ! X grid identical to equilibrium grid
676 ierr = grid_eq%copy(grid_x)
677 chckerr('')
678 case (2,3) ! solution and enriched
679 ! create grid
680 ierr = grid_x%init([grid_eq%n(1:2),size(r_f_x)],x_limits)
681 chckerr('')
682
683 ! set Flux variables
684 grid_x%r_F = r_f_x
685 grid_x%loc_r_F = r_f_x(x_limits(1):x_limits(2))
686
687 ! convert to Equilibrium variables
688 ierr = coord_f2e(grid_eq,grid_x%r_F,grid_x%r_E,&
689 &r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
690 chckerr('')
691 ierr = coord_f2e(grid_eq,grid_x%loc_r_F,grid_x%loc_r_E,&
692 &r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
693 chckerr('')
694
695 ! interpolate
696 do jd = 1,grid_eq%n(2)
697 do id = 1,grid_eq%n(1)
698 ierr = spline(grid_eq%loc_r_F,grid_eq%theta_E(id,jd,:),&
699 &grid_x%loc_r_F,grid_x%theta_E(id,jd,:),&
700 &ord=norm_disc_prec_x)
701 chckerr('')
702 ierr = spline(grid_eq%loc_r_F,grid_eq%zeta_E(id,jd,:),&
703 &grid_x%loc_r_F,grid_x%zeta_E(id,jd,:),&
704 &ord=norm_disc_prec_x)
705 chckerr('')
706 ierr = spline(grid_eq%loc_r_F,grid_eq%theta_F(id,jd,:),&
707 &grid_x%loc_r_F,grid_x%theta_F(id,jd,:),&
708 &ord=norm_disc_prec_x)
709 chckerr('')
710 ierr = spline(grid_eq%loc_r_F,grid_eq%zeta_F(id,jd,:),&
711 &grid_x%loc_r_F,grid_x%zeta_F(id,jd,:),&
712 &ord=norm_disc_prec_x)
713 chckerr('')
714 end do
715 end do
716 end select
717 end function setup_grid_x
718
719 !> Sets up the general solution grid, in which the solution variables are
720 !! calculated.
721 !!
722 !! For the solution grid, only one parallel point is used, but possibly
723 !! multiple geodesic points, equal to the number of field lines, \c n_alpha.
724 !!
725 !! \return ierr
726 integer function setup_grid_sol(grid_eq,grid_X,grid_sol,r_F_sol,&
727 &sol_limits) result(ierr)
728
729 use num_vars, only: n_procs, x_grid_style
730 use grid_utilities, only: coord_f2e
731
732 character(*), parameter :: rout_name = 'setup_grid_sol'
733
734 ! input / output
735 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
736 type(grid_type), intent(in) :: grid_x !< perturbation grid
737 type(grid_type), intent(inout) :: grid_sol !< solution grid
738 real(dp), intent(in) :: r_f_sol(:) !< points of solution grid
739 integer, intent(in) :: sol_limits(2) !< min. and max. index of sol grid of this process
740
741 ! initialize ierr
742 ierr = 0
743
744 ! output
745 call writo(trim(i2str(size(r_f_sol)))//' normal points')
746
747 select case (x_grid_style)
748 case (1,3) ! equilibrium or enriched
749 ! create grid
750 ierr = grid_sol%init([0,0,size(r_f_sol)],sol_limits,&
751 &divided=n_procs.gt.1)
752 chckerr('')
753
754 ! set Flux variables
755 grid_sol%r_F = r_f_sol
756 grid_sol%loc_r_F = r_f_sol(sol_limits(1):sol_limits(2))
757
758 ! convert to Equilibrium variables
759 ierr = coord_f2e(grid_eq,grid_sol%r_F,grid_sol%r_E,&
760 &r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
761 chckerr('')
762 ierr = coord_f2e(grid_eq,grid_sol%loc_r_F,grid_sol%loc_r_E,&
763 &r_f_array=grid_eq%r_F,r_e_array=grid_eq%r_E)
764 chckerr('')
765 case (2) ! solution
766 ! solution grid identical to perturation grid but with only 1
767 ! parallel point.
768 ierr = grid_sol%init([0,0,grid_x%n(3)],&
769 &[grid_x%i_min,grid_x%i_max],grid_x%divided)
770 chckerr('')
771 grid_sol%r_F = grid_x%r_F
772 grid_sol%loc_r_F = grid_x%r_F(sol_limits(1):sol_limits(2))
773 grid_sol%r_E = grid_x%r_E
774 grid_sol%loc_r_E = grid_x%r_E(sol_limits(1):sol_limits(2))
775 end select
776 end function setup_grid_sol
777
778 !> Calculate equilibrium grid that follows magnetic field lines.
779 !!
780 !! This grid is different from the equilibrium grid from setup_grid_eq for
781 !! HELENA, as the latter is the output grid from HELENA, which is situated
782 !! in a single poloidal cross-section, as opposed to a really field-aligne
783 !! grid.
784 !!
785 !! For VMEC, this is not used as the grid is field-aligned from the start.
786 !!
787 !! \note
788 !! -# The end-points are included for the grids in the parallel direction.
789 !! This is to facilitate working with the trapezoidal rule or Simpson's 3/8
790 !! rule for integration. This is \b NOT valid in general!
791 !! -# by setting the flag \c only_half_grid, only the even points of the
792 !! parallel grid are calculated, which is useful for higher Richardson
793 !! levels with VMEC so that only new angular points are calculated and the
794 !! old ones reused.
795 !!
796 !! \return ierr
797 integer function calc_ang_grid_eq_b(grid_eq,eq,only_half_grid) result(ierr)
798 use num_vars, only: use_pol_flux_f, use_pol_flux_e, &
799 &eq_style, tol_zero, eq_job_nr, eq_jobs_lims
800 use grid_vars, only: min_par_x, max_par_x, alpha, n_alpha
801 use eq_vars, only: max_flux_e
802 use grid_utilities, only: coord_f2e, calc_eqd_grid, calc_n_par_x_rich
803 use eq_utilities, only: print_info_eq
804 use rich_vars, only: n_par_x
805 use x_vars, only: min_r_sol, max_r_sol
806#if ldebug
807 use grid_utilities, only: coord_e2f
808#endif
809
810 character(*), parameter :: rout_name = 'calc_ang_grid_eq_B'
811
812 ! input / output
813 type(grid_type), intent(inout) :: grid_eq !< equilibrium grid of which to calculate angular part
814 type(eq_1_type), intent(in), target :: eq !< flux equilibrium variables
815 logical, intent(in), optional :: only_half_grid !< calculate only half grid with even points
816
817 ! local variables
818 character(len=max_str_ln) :: err_msg ! error message
819 real(dp), allocatable :: r_e_loc(:) ! flux in Equilibrium coords.
820 real(dp), pointer :: flux_f(:) => null() ! flux that the F uses as normal coord.
821 real(dp), pointer :: flux_e(:) => null() ! flux that the E uses as normal coord.
822 real(dp) :: r_f_factor, r_e_factor ! mult. factors for r_F and r_E
823 integer :: pmone ! plus or minus one
824 integer :: id, jd, kd ! counters
825 integer :: n_par_x_loc ! local n_par_X
826 logical :: only_half_grid_loc ! local only_half_grid
827 real(dp), allocatable :: theta_f_loc(:,:,:) ! local theta_F
828 real(dp), allocatable :: zeta_f_loc(:,:,:) ! local zeta_F
829
830 ! initialize ierr
831 ierr = 0
832
833 call lvl_ud(1)
834
835 ! set local only_half_grid
836 only_half_grid_loc = .false.
837 if (present(only_half_grid)) only_half_grid_loc = only_half_grid
838
839 ! set local n_par_X_rich
840 ierr = calc_n_par_x_rich(n_par_x_loc,only_half_grid_loc)
841 chckerr('')
842
843 ! user output
844 call writo('for '//trim(i2str(grid_eq%n(3)))//&
845 &' values on normal range '//trim(r2strt(min_r_sol))//'..'//&
846 &trim(r2strt(max_r_sol)))
847 call writo('for '//trim(i2str(n_par_x))//' values on parallel &
848 &range '//trim(r2strt(min_par_x))//'pi..'//&
849 &trim(r2strt(max_par_x))//'pi')
850 call lvl_ud(1)
851 if (only_half_grid_loc) call writo('for this Richardson level, only &
852 &the even points are setup up')
853 call print_info_eq(n_par_x_loc)
854 call lvl_ud(-1)
855
856 ! set up flux_E and plus minus one
857 ! Note: this routine is similar to calc_loc_r, but that routine cannot
858 ! be used here because it is possible that the FD quantities are not yet
859 ! defined.
860 ! choose which equilibrium style is being used:
861 ! 1: VMEC
862 ! 2: HELENA
863 select case (eq_style)
864 case (1) ! VMEC
865 if (use_pol_flux_e) then ! normalized poloidal flux
866 flux_e => eq%flux_p_E(:,0)
867 else ! normalized toroidal flux
868 flux_e => eq%flux_t_E(:,0)
869 end if
870 r_e_factor = max_flux_e
871 pmone = -1 ! conversion VMEC LH -> RH coord. system
872 case (2) ! HELENA
873 flux_e => eq%flux_p_E(:,0)
874 r_e_factor = 2*pi
875 pmone = 1
876 end select
877
878 ! set up flux_F
879 if (use_pol_flux_f) then ! poloidal flux / 2pi
880 flux_f => eq%flux_p_E(:,0)
881 r_f_factor = 2*pi
882 else ! toroidal flux / 2pi
883 flux_f => eq%flux_t_E(:,0)
884 r_f_factor = pmone*2*pi ! possible conversion VMEC LH -> RH coord. system
885 end if
886
887 ! set up parallel angle in Flux coordinates on equidistant grid
888 ! and use this to calculate the other angle as well
889 allocate(theta_f_loc(n_par_x,n_alpha,grid_eq%loc_n_r))
890 allocate(zeta_f_loc(n_par_x,n_alpha,grid_eq%loc_n_r))
891 if (use_pol_flux_f) then ! parallel angle theta
892 ierr = calc_eqd_grid(theta_f_loc,min_par_x*pi,&
893 &max_par_x*pi,1,excl_last=.false.) ! first index corresponds to parallel angle
894 chckerr('')
895 do kd = 1,grid_eq%loc_n_r
896 zeta_f_loc(:,:,kd) = pmone*eq%q_saf_E(kd,0)*&
897 &theta_f_loc(:,:,kd)
898 end do
899 do jd = 1,n_alpha
900 zeta_f_loc(:,jd,:) = zeta_f_loc(:,jd,:) + alpha(jd)
901 end do
902 else ! parallel angle zeta
903 ierr = calc_eqd_grid(zeta_f_loc,min_par_x*pi,&
904 &max_par_x*pi,1,excl_last=.false.) ! first index corresponds to parallel angle
905 chckerr('')
906 do kd = 1,grid_eq%loc_n_r
907 theta_f_loc(:,:,kd) = pmone*eq%rot_t_E(kd,0)*&
908 &zeta_f_loc(:,:,kd)
909 end do
910 do jd = 1,n_alpha
911 theta_f_loc(:,jd,:) = theta_f_loc(:,jd,:) - alpha(jd)
912 end do
913 end if
914
915 ! set up grid_eq angular F coordinates
916 if (only_half_grid_loc) then
917 do id = eq_jobs_lims(1,eq_job_nr),eq_jobs_lims(2,eq_job_nr)
918 grid_eq%theta_F(id-eq_jobs_lims(1,eq_job_nr)+1,:,:) = &
919 &theta_f_loc(2*id,:,:)
920 grid_eq%zeta_F(id-eq_jobs_lims(1,eq_job_nr)+1,:,:) = &
921 &zeta_f_loc(2*id,:,:)
922 end do
923 else
924 grid_eq%theta_F = theta_f_loc(eq_jobs_lims(1,eq_job_nr):&
925 &eq_jobs_lims(2,eq_job_nr),:,:)
926 grid_eq%zeta_F = zeta_f_loc(eq_jobs_lims(1,eq_job_nr):&
927 &eq_jobs_lims(2,eq_job_nr),:,:)
928 end if
929
930 ! allocate local r_E
931 allocate(r_e_loc(size(flux_f)))
932
933 ! convert Flux coordinates to Equilibrium coordinates (use
934 ! custom flux_E and flux_F because the Flux quantities are not
935 ! yet calculated)
936 call writo('convert F to E coordinates')
937 call lvl_ud(1)
938 ierr = coord_f2e(grid_eq,grid_eq%loc_r_F,grid_eq%theta_F,&
939 &grid_eq%zeta_F,r_e_loc,grid_eq%theta_E,grid_eq%zeta_E,&
940 &r_f_array=flux_f/r_f_factor,r_e_array=flux_e/r_e_factor)
941 chckerr('')
942 call lvl_ud(-1)
943
944 ! test whether r_E_loc indeed corresponds to loc_r_E of eq grid
945 if (maxval(abs(grid_eq%loc_r_E-r_e_loc)).gt.10*tol_zero) then
946 ierr = 1
947 err_msg = 'loc_r_E of equilibrium grid is not recovered'
948 chckerr(err_msg)
949 end if
950
951#if ldebug
952 if (debug_calc_ang_grid_eq_b) then
953 ! test whether F variables recovered
954 deallocate(theta_f_loc,zeta_f_loc)
955 allocate(theta_f_loc(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
956 allocate(zeta_f_loc(grid_eq%n(1),grid_eq%n(2),grid_eq%loc_n_r))
957 ierr = coord_e2f(grid_eq,grid_eq%loc_r_E,grid_eq%theta_E,&
958 &grid_eq%zeta_E,grid_eq%loc_r_F,theta_f_loc,zeta_f_loc,&
959 &r_e_array=flux_e/r_e_factor,r_f_array=flux_f/r_f_factor)
960 chckerr('')
961
962 ! plot difference
963 call plot_diff_hdf5(grid_eq%theta_F,theta_f_loc,'TEST_theta_F',&
964 &grid_eq%n,[0,0,grid_eq%i_min-1],'test whether F variable are &
965 &recovered',output_message=.true.)
966 call plot_diff_hdf5(grid_eq%zeta_F,zeta_f_loc,'TEST_zeta_F',&
967 &grid_eq%n,[0,0,grid_eq%i_min-1],'test whether F variable are &
968 &recovered',output_message=.true.)
969 end if
970#endif
971
972 ! deallocate local variables
973 deallocate(r_e_loc)
974 nullify(flux_f,flux_e)
975
976 call lvl_ud(-1)
977 end function calc_ang_grid_eq_b
978
979 !> Redistribute a grid to match the normal distribution of solution grid.
980 !!
981 !! The routine first calculates the smallest eq range that comprises the sol
982 !! range. Then, it gets the lowest equilibrium limits able to setup an
983 !! output grid that starts at index 1. After determining the output grid, it
984 !! then sends the variables to their new processes using MPI.
985 !!
986 !! \note
987 !! -# Only the Flux variables are saved.
988 !! -# the redistributed grid has trimmed outer limits, i.e. it starts at 1
989 !! and ends at the upper limit of the last process. This can be turned off
990 !! optionally using \c no_outer_trim.
991 !!
992 !! \return ierr
993 integer function redistribute_output_grid(grid,grid_out,no_outer_trim) &
994 &result(ierr)
995 use grid_utilities, only: find_compr_range
996 use mpi_utilities, only: redistribute_var, get_ser_var
997 use grid_vars, only: n_r_sol
998 use num_vars, only: n_procs
999
1000 character(*), parameter :: rout_name = 'redistribute_output_grid'
1001
1002 ! input / output
1003 type(grid_type), intent(in) :: grid !< equilibrium grid variables
1004 type(grid_type), intent(inout) :: grid_out !< redistributed equilibrium grid variables
1005 logical, intent(in), optional :: no_outer_trim !< do not trim the outer limits
1006
1007 ! local variables
1008 integer :: id ! counter
1009 integer :: eq_limits(2) ! normal limits for equilibrium variables
1010 integer :: sol_limits(2) ! normal limits for perturbation variables
1011 integer :: n_out(3) ! n of grid_out
1012 integer :: i_lim_tot(2) ! total limits of grid
1013 integer :: i_lim_out(2) ! limits of grid_out
1014 integer :: lims(2), lims_dis(2) ! limits and distributed limits, taking into account the angular extent
1015 real(dp), allocatable :: r_f_sol(:) ! perturbation r_F
1016 real(dp), allocatable :: temp_var(:) ! temporary variable
1017 integer, allocatable :: temp_lim(:) ! temporary limit
1018 integer, allocatable :: eq_limits_tot(:,:) ! total equilibrium limits
1019 logical :: no_outer_trim_loc ! local no_outer_trim
1020
1021 ! initialize ierr
1022 ierr = 0
1023
1024 ! calculate normal range for solution and save in perturbation variables
1025 allocate(r_f_sol(n_r_sol))
1026 ierr = calc_norm_range('PB3D_sol',sol_limits=sol_limits,r_f_sol=r_f_sol)
1027 chckerr('')
1028
1029 ! determine eq_limits: smallest eq range comprising X range
1030 call find_compr_range(grid%r_F,r_f_sol(sol_limits),eq_limits)
1031
1032 ! get lowest equilibrium limits to be able to setup an output grid that
1033 ! starts at index 1
1034 allocate(eq_limits_tot(2,n_procs))
1035 do id = 1,2
1036 ierr = get_ser_var(eq_limits(id:id),temp_lim,scatter=.true.)
1037 chckerr('')
1038 eq_limits_tot(id,:) = temp_lim
1039 end do
1040
1041 ! set up redistributed grid
1042 no_outer_trim_loc = .false.
1043 if (present(no_outer_trim)) no_outer_trim_loc = no_outer_trim
1044 if (no_outer_trim_loc) then
1045 n_out = grid%n
1046 i_lim_tot = [1,n_out(3)]
1047 i_lim_out = eq_limits
1048 else
1049 n_out(1:2) = grid%n(1:2)
1050 n_out(3) = eq_limits_tot(2,n_procs)-eq_limits_tot(1,1)+1
1051 i_lim_tot = [eq_limits_tot(1,1),eq_limits_tot(2,n_procs)]
1052 i_lim_out = eq_limits-eq_limits_tot(1,1)+1
1053 end if
1054 ierr = grid_out%init(n_out,i_lim_out)
1055 chckerr('')
1056
1057 ! set up limits taking into account angular extent and temporary var
1058 lims(1) = product(grid%n(1:2))*(grid%i_min-1)+1
1059 lims(2) = product(grid%n(1:2))*grid%i_max
1060 lims_dis(1) = product(grid%n(1:2))*(eq_limits(1)-1)+1
1061 lims_dis(2) = product(grid%n(1:2))*eq_limits(2)
1062 allocate(temp_var(lims_dis(2)-lims_dis(1)+1))
1063
1064 ! copy total variables
1065 ! r_F
1066 grid_out%r_F = grid%r_F(i_lim_tot(1):i_lim_tot(2))
1067 ! r_E
1068 grid_out%r_E = grid%r_E(i_lim_tot(1):i_lim_tot(2))
1069
1070 ! distribute local variables
1071 ! theta_F
1072 ierr = redistribute_var(reshape(grid%theta_F,[size(grid%theta_F)]),&
1073 &temp_var,lims,lims_dis)
1074 chckerr('')
1075 grid_out%theta_F = reshape(temp_var,shape(grid_out%theta_F))
1076 ! zeta_F
1077 ierr = redistribute_var(reshape(grid%zeta_F,[size(grid%zeta_F)]),&
1078 &temp_var,lims,lims_dis)
1079 chckerr('')
1080 grid_out%zeta_F = reshape(temp_var,shape(grid_out%zeta_F))
1081 ! theta_E
1082 ierr = redistribute_var(reshape(grid%theta_E,[size(grid%theta_E)]),&
1083 &temp_var,lims,lims_dis)
1084 chckerr('')
1085 grid_out%theta_E = reshape(temp_var,shape(grid_out%theta_E))
1086 ! zeta_E
1087 ierr = redistribute_var(reshape(grid%zeta_E,[size(grid%zeta_E)]),&
1088 &temp_var,lims,lims_dis)
1089 chckerr('')
1090 grid_out%zeta_E = reshape(temp_var,shape(grid_out%zeta_E))
1091 ! loc_r_F and loc_R_E
1092 if (grid_out%divided) then
1093 grid_out%loc_r_F = grid_out%r_F(grid_out%i_min:grid_out%i_max)
1094 grid_out%loc_r_E = grid_out%r_E(grid_out%i_min:grid_out%i_max)
1095 end if
1096 end function redistribute_output_grid
1097
1098 !> Plots the grid in real 3-D space.
1099 !!
1100 !! This creates an animation that can be used by ParaView or VisIt.
1101 !!
1102 !! The equilibrium grid should contain the fieldline-oriented angles with \c
1103 !! ang_1 the parallel angle and \c ang_2 the field line label.
1104 !!
1105 !! \see See grid_vars.grid_type for a discussion on \c ang_1 and \c ang_2.
1106 !!
1107 !! \note
1108 !! -# This procedure does not use \c n_theta_plot and \c n_zeta_plot from
1109 !! num_vars, but instead temporarily overwrites them with its own, since it
1110 !! is suposed to be 3-D also in the axisymmetric case.
1111 !! -# The implementation is currently very slow.
1112 !!
1113 !! \return ierr
1114 integer function magn_grid_plot(grid) result(ierr)
1115 use num_vars, only: rank, no_plots, n_theta_plot, n_zeta_plot, &
1116 &eq_style, min_theta_plot, max_theta_plot, min_zeta_plot, &
1117 &max_zeta_plot
1118 use grid_utilities, only: trim_grid, extend_grid_f, calc_xyz_grid
1119 use vmec_utilities, only: calc_trigon_factors
1120
1121 character(*), parameter :: rout_name = 'magn_grid_plot'
1122
1123 ! input / output
1124 type(grid_type), intent(in) :: grid !< fieldline-oriented equilibrium grid
1125
1126 ! local variables
1127 real(dp), allocatable :: x_1(:,:,:), y_1(:,:,:), z_1(:,:,:) ! X, Y and Z of surface in Axisymmetric coordinates
1128 real(dp), allocatable :: x_2(:,:,:), y_2(:,:,:), z_2(:,:,:) ! X, Y and Z of magnetic field lines in Axisymmetric coordinates
1129 real(dp), pointer :: x_1_tot(:,:,:) => null() ! total X
1130 real(dp), pointer :: y_1_tot(:,:,:) => null() ! total Y
1131 real(dp), pointer :: z_1_tot(:,:,:) => null() ! total Z
1132 real(dp), pointer :: x_2_tot(:,:,:) => null() ! total X
1133 real(dp), pointer :: y_2_tot(:,:,:) => null() ! total Y
1134 real(dp), pointer :: z_2_tot(:,:,:) => null() ! total Z
1135 type(grid_type) :: grid_ext ! angularly extended grid
1136 type(grid_type) :: grid_plot ! grid for plotting
1137 integer :: n_theta_plot_old ! backup of n_theta_plot
1138 integer :: n_zeta_plot_old ! backup of n_zeta_plot
1139 real(dp) :: min_theta_plot_old, max_theta_plot_old ! backup of min and max_theta_plot
1140 real(dp) :: min_zeta_plot_old, max_zeta_plot_old ! backup of min and max_zeta_plot
1141 character(len=max_str_ln) :: anim_name ! name of animation
1142
1143 ! initialize ierr
1144 ierr = 0
1145
1146 ! bypass plots if no_plots
1147 if (no_plots) return
1148
1149 call writo('Plotting magnetic field and flux surfaces')
1150
1151 call lvl_ud(1)
1152
1153 ! 0. set up variables
1154
1155 ! save n_theta_plot and n_zeta_plot and change them
1156 n_theta_plot_old = n_theta_plot
1157 n_zeta_plot_old = n_zeta_plot
1158 min_theta_plot_old = min_theta_plot
1159 max_theta_plot_old = max_theta_plot
1160 min_zeta_plot_old = min_zeta_plot
1161 max_zeta_plot_old = max_zeta_plot
1162 n_theta_plot = 40
1163 n_zeta_plot = 160
1164 min_theta_plot = 1
1165 max_theta_plot = 3
1166 min_zeta_plot = 0
1167 max_zeta_plot = 2
1168
1169 ! extend grid
1170 ierr = extend_grid_f(grid,grid_ext,grid_eq=grid)
1171 chckerr('')
1172
1173 ! restore n_theta_plot and n_zeta_plot
1174 n_theta_plot = n_theta_plot_old
1175 n_zeta_plot = n_zeta_plot_old
1176 min_theta_plot = min_theta_plot_old
1177 max_theta_plot = max_theta_plot_old
1178 min_zeta_plot = min_zeta_plot_old
1179 max_zeta_plot = max_zeta_plot_old
1180
1181 ! trim extended grid into plot grid
1182 ierr = trim_grid(grid_ext,grid_plot)
1183 chckerr('')
1184
1185 ! if VMEC, calculate trigonometric factors of plot grid
1186 if (eq_style.eq.1) then
1187 ierr = calc_trigon_factors(grid_plot%theta_E,grid_plot%zeta_E,&
1188 &grid_plot%trigon_factors)
1189 chckerr('')
1190 end if
1191
1192 ! set animation name
1193 anim_name = 'Magnetic field in flux surfaces'
1194
1195 ! 1. plot flux surfaces
1196 call writo('writing flux surfaces')
1197
1198 ! calculate X_1,Y_1 and Z_1
1199 allocate(x_1(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1200 allocate(y_1(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1201 allocate(z_1(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1202 ierr = calc_xyz_grid(grid,grid_plot,x_1,y_1,z_1)
1203 chckerr('')
1204
1205 ! dealloc grid
1206 call grid_plot%dealloc()
1207
1208 ! get pointers to full X, Y and Z
1209 ! The reason for this is that the plot is not as simple as usual, so no
1210 ! divided plots are used, and also efficiency is not the biggest
1211 ! priority. Therefore, all the plotting of the local is handled by a
1212 ! single process, the master.
1213 call get_full_xyz(x_1,y_1,z_1,x_1_tot,y_1_tot,z_1_tot,'flux surfaces')
1214
1215 ! 2. plot field lines
1216 call writo('writing field lines')
1217
1218 ! trim grid into plot grid
1219 ierr = trim_grid(grid,grid_plot)
1220 chckerr('')
1221
1222 ! if VMEC, calculate trigonometric factors of plot grid
1223 if (eq_style.eq.1) then
1224 ierr = calc_trigon_factors(grid_plot%theta_E,grid_plot%zeta_E,&
1225 &grid_plot%trigon_factors)
1226 chckerr('')
1227 end if
1228
1229 ! calculate X_2,Y_2 and Z_2
1230 allocate(x_2(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1231 allocate(y_2(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1232 allocate(z_2(grid_plot%n(1),grid_plot%n(2),grid_plot%loc_n_r))
1233 ierr = calc_xyz_grid(grid,grid_plot,x_2,y_2,z_2)
1234 chckerr('')
1235
1236 ! dealloc grids
1237 call grid_plot%dealloc()
1238 call grid_ext%dealloc()
1239
1240 ! get pointers to full X, Y and Z
1241 ! The reason for this is that the plot is not as simple as usual, so no
1242 ! divided plots are used, and also efficiency is not the biggest
1243 ! priority. Therefore, all the plotting of the local is handled by a
1244 ! single process, the master.
1245 call get_full_xyz(x_2,y_2,z_2,x_2_tot,y_2_tot,z_2_tot,'field lines')
1246
1247 ierr = magn_grid_plot_hdf5(x_1_tot,x_2_tot,y_1_tot,y_2_tot,&
1248 &z_1_tot,z_2_tot,anim_name)
1249 chckerr('')
1250
1251 ! clean up
1252 deallocate(x_1,y_1,z_1)
1253 deallocate(x_2,y_2,z_2)
1254 nullify(x_1_tot,y_1_tot,z_1_tot)
1255 nullify(x_2_tot,y_2_tot,z_2_tot)
1256
1257 call lvl_ud(-1)
1258
1259 call writo('Done plotting magnetic field and flux surfaces')
1260 contains
1261 ! get pointer to full plotting variables X, Y and Z
1262 !> \private
1263 subroutine get_full_xyz(X,Y,Z,X_tot,Y_tot,Z_tot,merge_name)
1264 use mpi_utilities, only: get_ser_var
1265
1266 ! input / output
1267 real(dp), intent(in), target :: x(:,:,:), y(:,:,:), z(:,:,:) ! X, Y and Z of either flux surfaces or magnetic field lines
1268 real(dp), intent(inout), pointer :: x_tot(:,:,:) ! pointer to full X
1269 real(dp), intent(inout), pointer :: y_tot(:,:,:) ! pointer to full Y
1270 real(dp), intent(inout), pointer :: z_tot(:,:,:) ! pointer to full Z
1271 character(len=*) :: merge_name ! name of variable to be merged
1272
1273 ! local variables
1274 real(dp), allocatable :: ser_xyz_loc(:) ! serial copy of XYZ_loc
1275 integer, allocatable :: tot_dim(:) ! total dimensions for plot
1276
1277 ! merge plots for flux surfaces if more than one process
1278 if (grid%divided) then ! merge local plots
1279 ! user output
1280 call writo('Merging local plots for '//merge_name)
1281
1282 ! get total dimension
1283 ierr = get_ser_var([size(x,3)],tot_dim)
1284 chckerr('')
1285
1286 ! allocate total X, Y and Z (should have same sizes)
1287 if (rank.eq.0) then
1288 allocate(x_tot(size(x,1),size(x,2),sum(tot_dim)))
1289 allocate(y_tot(size(x,1),size(x,2),sum(tot_dim)))
1290 allocate(z_tot(size(x,1),size(x,2),sum(tot_dim)))
1291 end if
1292
1293 allocate(ser_xyz_loc(size(x(:,:,1))*sum(tot_dim)))
1294 ierr = get_ser_var(reshape(x,[size(x)]),ser_xyz_loc)
1295 chckerr('')
1296 if (rank.eq.0) x_tot = reshape(ser_xyz_loc,shape(x_tot))
1297 ierr = get_ser_var(reshape(y,[size(y)]),ser_xyz_loc)
1298 chckerr('')
1299 if (rank.eq.0) y_tot = reshape(ser_xyz_loc,shape(y_tot))
1300 ierr = get_ser_var(reshape(z,[size(z)]),ser_xyz_loc)
1301 chckerr('')
1302 if (rank.eq.0) z_tot = reshape(ser_xyz_loc,shape(z_tot))
1303 else ! just point
1304 x_tot => x
1305 y_tot => y
1306 z_tot => z
1307 end if
1308 end subroutine get_full_xyz
1309
1310 ! Plot with HDF5
1311 !> \private
1312 integer function magn_grid_plot_hdf5(X_1,X_2,Y_1,Y_2,Z_1,Z_2,&
1313 &anim_name) result(ierr)
1314 use hdf5_ops, only: open_hdf5_file, add_hdf5_item, &
1315 &print_hdf5_top, print_hdf5_geom, print_hdf5_3d_data_item, &
1316 &print_hdf5_grid, close_hdf5_file
1317 use hdf5_vars, only: dealloc_xml_str, &
1318 & xml_str_type, hdf5_file_type
1319 use rich_vars, only: rich_lvl
1320 use grid_vars, only: n_alpha
1321
1322 character(*), parameter :: rout_name = 'magn_grid_plot_HDF5'
1323
1324 ! input / output
1325 real(dp), intent(in), pointer :: x_1(:,:,:), y_1(:,:,:), z_1(:,:,:) ! X, Y and Z of surface in Axisymmetric coordinates
1326 real(dp), intent(in), pointer :: x_2(:,:,:), y_2(:,:,:), z_2(:,:,:) ! X, Y and Z of magnetic field lines in Axisymmetric coordinates
1327 character(len=*), intent(in) :: anim_name ! name of animation
1328
1329 ! local variables
1330 character(len=max_str_ln) :: plot_title(2) ! plot title for flux surface and for field lines
1331 character(len=max_str_ln) :: file_name ! name of file
1332 integer :: id, jd ! counter
1333 type(hdf5_file_type) :: file_info ! file info
1334 type(xml_str_type), allocatable :: grids(:) ! grid in spatial collection (flux surface, field line)
1335 type(xml_str_type), allocatable :: space_col_grids(:) ! grid with space collection at different times
1336 type(xml_str_type) :: time_col_grid ! grid with time collection
1337 type(xml_str_type) :: top ! topology
1338 type(xml_str_type) :: xyz(3) ! data items for geometry
1339 type(xml_str_type) :: geom ! geometry
1340 integer :: loc_dim(3,2) ! local dimensions (flux surfaces, field lines)
1341 integer :: n_r ! nr. of normal points
1342
1343 ! initialize ierr
1344 ierr = 0
1345
1346 ! only master
1347 if (rank.eq.0) then
1348 ! user output
1349 call writo('Drawing animation with HDF5')
1350 call lvl_ud(1)
1351
1352 ! set up loc_dim and n_r
1353 loc_dim(:,1) = [size(x_1,1),size(x_1,2),1]
1354 loc_dim(:,2) = [size(x_2,1),1,1]
1355 n_r = size(x_1,3) - 1 ! should be same for all other X_i, Y_i and Z_i
1356
1357 ! set up plot titles and file name
1358 plot_title(1) = 'Magnetic Flux Surfaces'
1359 plot_title(2) = 'Magnetic Field Lines'
1360 file_name = 'field_lines_in_flux_surfaces_R_'//&
1361 &trim(i2str(rich_lvl))
1362
1363 ! open HDF5 file
1364 ierr = open_hdf5_file(file_info,trim(file_name),&
1365 &descr=anim_name,ind_plot=.true.)
1366 chckerr('')
1367
1368 ! create grid for flux surface and all field lines and one for
1369 ! time collection
1370 allocate(grids(1+n_alpha))
1371 allocate(space_col_grids(n_r))
1372
1373 ! loop over all normal points
1374 do id = 1,n_r
1375 ! A. start with flux surface
1376
1377 ! print topology
1378 call print_hdf5_top(top,2,loc_dim(:,1))
1379
1380 ! print data item for X
1381 ierr = print_hdf5_3d_data_item(xyz(1),file_info,&
1382 &'X_surf_'//trim(i2str(id)),x_1(:,:,id:id),loc_dim(:,1))
1383 chckerr('')
1384
1385 ! print data item for Y
1386 ierr = print_hdf5_3d_data_item(xyz(2),file_info,&
1387 &'Y_surf_'//trim(i2str(id)),y_1(:,:,id:id),loc_dim(:,1))
1388 chckerr('')
1389
1390 ! print data item for Z
1391 ierr = print_hdf5_3d_data_item(xyz(3),file_info,&
1392 &'Z_surf_'//trim(i2str(id)),z_1(:,:,id:id),loc_dim(:,1))
1393 chckerr('')
1394
1395 ! print geometry with X, Y and Z data item
1396 call print_hdf5_geom(geom,2,xyz,.true.)
1397
1398 ! create a grid with the topology and the geometry
1399 ierr = print_hdf5_grid(grids(1),plot_title(1),1,&
1400 &grid_top=top,grid_geom=geom,reset=.true.)
1401 chckerr('')
1402
1403 ! B. end with magnetic field
1404
1405 do jd = 1,n_alpha ! iterate over all field lines
1406 ! print topology
1407 call print_hdf5_top(top,2,loc_dim(:,2))
1408
1409 ! print data item for X of this field line
1410 ierr = print_hdf5_3d_data_item(xyz(1),file_info,&
1411 &'X_field_'//trim(i2str(id))//'_'//trim(i2str(jd)),&
1412 &x_2(:,jd:jd,id+1:id+1),loc_dim(:,2))
1413 chckerr('')
1414
1415 ! print data item for Y of this field line
1416 ierr = print_hdf5_3d_data_item(xyz(2),file_info,&
1417 &'Y_field_'//trim(i2str(id))//'_'//trim(i2str(jd)),&
1418 &y_2(:,jd:jd,id+1:id+1),loc_dim(:,2))
1419 chckerr('')
1420
1421 ! print data item for Z of this field line
1422 ierr = print_hdf5_3d_data_item(xyz(3),file_info,&
1423 &'Z_field_'//trim(i2str(id))//'_'//trim(i2str(jd)),&
1424 &z_2(:,jd:jd,id+1:id+1),loc_dim(:,2))
1425 chckerr('')
1426
1427 ! print geometry with X, Y and Z data item
1428 call print_hdf5_geom(geom,2,xyz,.true.)
1429
1430 ! create a grid with the topology and the geometry
1431 ierr = print_hdf5_grid(grids(jd+1),plot_title(2),1,&
1432 &grid_top=top,grid_geom=geom,reset=.true.)
1433 chckerr('')
1434 end do
1435
1436 ! C. merge flux surface and magnetic field in to spatial
1437 ! collection
1438 ierr = print_hdf5_grid(space_col_grids(id),&
1439 &'spatial collection',3,grid_time=id*1._dp,&
1440 &grid_grids=grids,reset=.true.)
1441 chckerr('')
1442 end do
1443
1444 ! create grid collection from individual grids and reset them
1445 ierr = print_hdf5_grid(time_col_grid,'time collection',2,&
1446 &grid_grids=space_col_grids,reset=.true.)
1447 chckerr('')
1448
1449 ! add collection grid to HDF5 file and reset it
1450 ierr = add_hdf5_item(file_info,time_col_grid,reset=.true.)
1451 chckerr('')
1452
1453 ! close HDF5 file
1454 ierr = close_hdf5_file(file_info)
1455 chckerr('')
1456
1457 call lvl_ud(-1)
1458
1459 ! clean up
1460 do id = 1,2
1461 call dealloc_xml_str(grids(id))
1462 end do
1463 call dealloc_xml_str(space_col_grids)
1464 call dealloc_xml_str(time_col_grid)
1465 call dealloc_xml_str(top)
1466 do id = 1,3
1467 call dealloc_xml_str(xyz(id))
1468 end do
1469 call dealloc_xml_str(geom)
1470 end if
1471 end function magn_grid_plot_hdf5
1472 end function magn_grid_plot
1473
1474 !> Print grid variables to an output file.
1475 !!
1476 !! If \c rich_lvl is provided, <tt>_R_[rich_lvl]</tt> is appended to the
1477 !! data name if it is > 0.
1478 !!
1479 !! Optionally, it can be specified that this is a divided parallel grid,
1480 !! corresponding to the variable \c eq_jobs_lims with index \c eq_job_nr. In
1481 !! this case, the total grid size is adjusted to the one specified by \c
1482 !! eq_jobs_lims and the grid is written as a subset.
1483 !!
1484 !! \note <tt>grid_</tt> is added in front the data_name.
1485 !!
1486 !! \return ierr
1487 integer function print_output_grid(grid,grid_name,data_name,rich_lvl,&
1488 &par_div,remove_previous_arrs) result(ierr)
1489 use num_vars, only: pb3d_name, eq_jobs_lims, eq_job_nr
1490 use hdf5_ops, only: print_hdf5_arrs
1491 use hdf5_vars, only: var_1d_type, &
1492 &max_dim_var_1d
1493 use grid_utilities, only: trim_grid
1494
1495 character(*), parameter :: rout_name = 'print_output_grid'
1496
1497 ! input / output
1498 type(grid_type), intent(in) :: grid !< grid variables
1499 character(len=*), intent(in) :: grid_name !< name to display
1500 character(len=*), intent(in) :: data_name !< name under which to store
1501 integer, intent(in), optional :: rich_lvl !< Richardson level to reconstruct
1502 logical, intent(in), optional :: par_div !< is a parallely divided grid
1503 logical, intent(in), optional :: remove_previous_arrs !< remove previous variables if present
1504
1505 ! local variables
1506 integer :: n_tot(3) ! total n
1507 integer :: par_id(2) ! local parallel interval
1508 type(grid_type) :: grid_trim ! trimmed grid
1509 type(var_1d_type), allocatable, target :: grid_1d(:) ! 1D equivalent of grid variables
1510 type(var_1d_type), pointer :: grid_1d_loc => null() ! local element in grid_1D
1511 integer :: id ! counter
1512 logical :: par_div_loc ! local par_div
1513
1514 ! initialize ierr
1515 ierr = 0
1516
1517 ! user output
1518 call writo('Write '//trim(grid_name)//' grid variables to output &
1519 &file')
1520 call lvl_ud(1)
1521
1522 ! trim grid
1523 ierr = trim_grid(grid,grid_trim)
1524 chckerr('')
1525
1526 ! set local par_div
1527 par_div_loc = .false.
1528 if (present(par_div)) par_div_loc = par_div
1529
1530 ! set total n and parallel interval
1531 n_tot = grid_trim%n
1532 par_id = [1,n_tot(1)]
1533 if (grid_trim%n(1).gt.0 .and. par_div_loc) then ! total grid includes all equilibrium jobs
1534 n_tot(1) = maxval(eq_jobs_lims)-minval(eq_jobs_lims)+1
1535 par_id = eq_jobs_lims(:,eq_job_nr)
1536 end if
1537
1538 ! Set up the 1D equivalents of the perturbation variables
1539 allocate(grid_1d(max_dim_var_1d))
1540
1541 ! set up variables grid_1D
1542 id = 1
1543
1544 ! n
1545 grid_1d_loc => grid_1d(id); id = id+1
1546 grid_1d_loc%var_name = 'n'
1547 allocate(grid_1d_loc%tot_i_min(1),grid_1d_loc%tot_i_max(1))
1548 allocate(grid_1d_loc%loc_i_min(1),grid_1d_loc%loc_i_max(1))
1549 grid_1d_loc%tot_i_min = [1]
1550 grid_1d_loc%tot_i_max = [3]
1551 grid_1d_loc%loc_i_min = [1]
1552 grid_1d_loc%loc_i_max = [3]
1553 allocate(grid_1d_loc%p(3))
1554 grid_1d_loc%p = n_tot
1555
1556 ! r_F
1557 grid_1d_loc => grid_1d(id); id = id+1
1558 grid_1d_loc%var_name = 'r_F'
1559 allocate(grid_1d_loc%tot_i_min(1),grid_1d_loc%tot_i_max(1))
1560 allocate(grid_1d_loc%loc_i_min(1),grid_1d_loc%loc_i_max(1))
1561 grid_1d_loc%tot_i_min = [1]
1562 grid_1d_loc%tot_i_max = [n_tot(3)]
1563 grid_1d_loc%loc_i_min = [grid_trim%i_min]
1564 grid_1d_loc%loc_i_max = [grid_trim%i_max]
1565 allocate(grid_1d_loc%p(size(grid_trim%loc_r_F)))
1566 grid_1d_loc%p = grid_trim%loc_r_F
1567
1568 ! r_E
1569 grid_1d_loc => grid_1d(id); id = id+1
1570 grid_1d_loc%var_name = 'r_E'
1571 allocate(grid_1d_loc%tot_i_min(1),grid_1d_loc%tot_i_max(1))
1572 allocate(grid_1d_loc%loc_i_min(1),grid_1d_loc%loc_i_max(1))
1573 grid_1d_loc%tot_i_min = [1]
1574 grid_1d_loc%tot_i_max = [n_tot(3)]
1575 grid_1d_loc%loc_i_min = [grid_trim%i_min]
1576 grid_1d_loc%loc_i_max = [grid_trim%i_max]
1577 allocate(grid_1d_loc%p(size(grid_trim%loc_r_E)))
1578 grid_1d_loc%p = grid_trim%loc_r_E
1579
1580 ! only for 3D grids
1581 if (product(grid%n(1:2)).ne.0) then
1582 ! theta_F
1583 grid_1d_loc => grid_1d(id); id = id+1
1584 grid_1d_loc%var_name = 'theta_F'
1585 allocate(grid_1d_loc%tot_i_min(3),grid_1d_loc%tot_i_max(3))
1586 allocate(grid_1d_loc%loc_i_min(3),grid_1d_loc%loc_i_max(3))
1587 grid_1d_loc%tot_i_min = [1,1,1]
1588 grid_1d_loc%tot_i_max = n_tot
1589 grid_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
1590 grid_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
1591 allocate(grid_1d_loc%p(size(grid_trim%theta_F)))
1592 grid_1d_loc%p = reshape(grid_trim%theta_F,[size(grid_trim%theta_F)])
1593
1594 ! theta_E
1595 grid_1d_loc => grid_1d(id); id = id+1
1596 grid_1d_loc%var_name = 'theta_E'
1597 allocate(grid_1d_loc%tot_i_min(3),grid_1d_loc%tot_i_max(3))
1598 allocate(grid_1d_loc%loc_i_min(3),grid_1d_loc%loc_i_max(3))
1599 grid_1d_loc%tot_i_min = [1,1,1]
1600 grid_1d_loc%tot_i_max = n_tot
1601 grid_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
1602 grid_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
1603 allocate(grid_1d_loc%p(size(grid_trim%theta_E)))
1604 grid_1d_loc%p = reshape(grid_trim%theta_E,[size(grid_trim%theta_E)])
1605
1606 ! zeta_F
1607 grid_1d_loc => grid_1d(id); id = id+1
1608 grid_1d_loc%var_name = 'zeta_F'
1609 allocate(grid_1d_loc%tot_i_min(3),grid_1d_loc%tot_i_max(3))
1610 allocate(grid_1d_loc%loc_i_min(3),grid_1d_loc%loc_i_max(3))
1611 grid_1d_loc%tot_i_min = [1,1,1]
1612 grid_1d_loc%tot_i_max = n_tot
1613 grid_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
1614 grid_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
1615 allocate(grid_1d_loc%p(size(grid_trim%zeta_F)))
1616 grid_1d_loc%p = reshape(grid_trim%zeta_F,[size(grid_trim%zeta_F)])
1617
1618 ! zeta_E
1619 grid_1d_loc => grid_1d(id); id = id+1
1620 grid_1d_loc%var_name = 'zeta_E'
1621 allocate(grid_1d_loc%tot_i_min(3),grid_1d_loc%tot_i_max(3))
1622 allocate(grid_1d_loc%loc_i_min(3),grid_1d_loc%loc_i_max(3))
1623 grid_1d_loc%tot_i_min = [1,1,1]
1624 grid_1d_loc%tot_i_max = n_tot
1625 grid_1d_loc%loc_i_min = [par_id(1),1,grid_trim%i_min]
1626 grid_1d_loc%loc_i_max = [par_id(2),n_tot(2),grid_trim%i_max]
1627 allocate(grid_1d_loc%p(size(grid_trim%zeta_E)))
1628 grid_1d_loc%p = reshape(grid_trim%zeta_E,[size(grid_trim%zeta_E)])
1629 end if
1630
1631 ! write
1632 ierr = print_hdf5_arrs(grid_1d(1:id-1),pb3d_name,&
1633 &'grid_'//trim(data_name),rich_lvl=rich_lvl,&
1634 &ind_print=.not.grid_trim%divided,&
1635 &remove_previous_arrs=remove_previous_arrs)
1636 chckerr('')
1637
1638 ! clean up
1639 call grid_trim%dealloc()
1640
1641 ! clean up
1642 nullify(grid_1d_loc)
1643
1644 ! user output
1645 call lvl_ud(-1)
1646 end function print_output_grid
1647end module grid_ops
calc_norm_range_post
subroutine calc_norm_range_post(eq_limits, x_limits, sol_limits, r_f_eq, r_f_x, r_f_sol)
POST version.
Definition grid_ops.f90:453
calc_norm_range_pb3d_sol
integer function calc_norm_range_pb3d_sol(sol_limits, r_f_sol, n_procs)
PB3D solution version.
Definition grid_ops.f90:388
calc_norm_range_pb3d_eq
subroutine calc_norm_range_pb3d_eq(eq_limits)
PB3D equilibrium version.
Definition grid_ops.f90:246
calc_norm_range_pb3d_in
integer function calc_norm_range_pb3d_in(in_limits)
PB3D input version.
Definition grid_ops.f90:137
calc_norm_range_pb3d_x
integer function calc_norm_range_pb3d_x(eq_limits, x_limits, r_f_eq, r_f_x, jq)
PB3D perturbation version.
Definition grid_ops.f90:285
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_e2f
Converts Equilibrium coordinates . to Flux coordinates .
Definition grid_utilities.f90:66
grid_utilities::coord_f2e
Converts Flux coordinates to Equilibrium coordinates .
Definition grid_utilities.f90:47
hdf5_vars::dealloc_xml_str
Deallocates XML_str_type.
Definition HDF5_vars.f90:60
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition MPI_utilities.f90:55
num_utilities::con2dis
Convert between points from a continuous grid to a discrete grid.
Definition num_utilities.f90:188
num_utilities::dis2con
Convert between points from a discrete grid to a continuous grid.
Definition num_utilities.f90:206
num_utilities::round_with_tol
Rounds an arry of values to limits, with a tolerance that can optionally be modified.
Definition num_utilities.f90:169
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
eq_utilities
Numerical utilities related to equilibrium variables.
Definition eq_utilities.f90:4
eq_utilities::print_info_eq
subroutine, public print_info_eq(n_par_x_rich)
Prints information for equilibrium parallel job.
Definition eq_utilities.f90:986
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::max_flux_e
real(dp), public max_flux_e
max. flux in Equilibrium coordinates, set in calc_norm_range_PB3D_in
Definition eq_vars.f90:49
grid_ops
Operations that have to do with the grids and different coordinate systems.
Definition grid_ops.f90:4
grid_ops::calc_ang_grid_eq_b
integer function, public calc_ang_grid_eq_b(grid_eq, eq, only_half_grid)
Calculate equilibrium grid that follows magnetic field lines.
Definition grid_ops.f90:798
grid_ops::print_output_grid
integer function, public print_output_grid(grid, grid_name, data_name, rich_lvl, par_div, remove_previous_arrs)
Print grid variables to an output file.
Definition grid_ops.f90:1489
grid_ops::calc_norm_range
integer function, public calc_norm_range(style, in_limits, eq_limits, x_limits, sol_limits, r_f_eq, r_f_x, r_f_sol, jq)
Calculates normal range for the input grid, the equilibrium grid and/or the solution grid.
Definition grid_ops.f90:46
grid_ops::redistribute_output_grid
integer function, public redistribute_output_grid(grid, grid_out, no_outer_trim)
Redistribute a grid to match the normal distribution of solution grid.
Definition grid_ops.f90:995
grid_ops::setup_grid_sol
integer function, public setup_grid_sol(grid_eq, grid_x, grid_sol, r_f_sol, sol_limits)
Sets up the general solution grid, in which the solution variables are calculated.
Definition grid_ops.f90:728
grid_ops::debug_calc_ang_grid_eq_b
logical, public debug_calc_ang_grid_eq_b
plot debug information for calc_ang_grid_eq_b()
Definition grid_ops.f90:25
grid_ops::setup_grid_eq_b
integer function, public setup_grid_eq_b(grid_eq, grid_eq_b, eq, only_half_grid)
Sets up the field-aligned equilibrium grid.
Definition grid_ops.f90:601
grid_ops::magn_grid_plot
integer function, public magn_grid_plot(grid)
Plots the grid in real 3-D space.
Definition grid_ops.f90:1115
grid_ops::setup_grid_eq
integer function, public setup_grid_eq(grid_eq, eq_limits)
Sets up the equilibrium grid.
Definition grid_ops.f90:529
grid_ops::setup_grid_x
integer function, public setup_grid_x(grid_eq, grid_x, r_f_x, x_limits)
Sets up the general perturbation grid, in which the perturbation variables are calculated.
Definition grid_ops.f90:655
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition grid_utilities.f90:4
grid_utilities::find_compr_range
subroutine, public find_compr_range(r_f, lim_r, lim_id)
finds smallest range that comprises a minimum and maximum value.
Definition grid_utilities.f90:2995
grid_utilities::extend_grid_f
integer function, public extend_grid_f(grid_in, grid_ext, grid_eq, n_theta_plot, n_zeta_plot, lim_theta_plot, lim_zeta_plot)
Extend a grid angularly.
Definition grid_utilities.f90:1096
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_n_par_x_rich
integer function, public calc_n_par_x_rich(n_par_x_rich, only_half_grid)
Calculates the local number of parallel grid points for this Richardson level, taking into account th...
Definition grid_utilities.f90:2862
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::n_r_eq
integer, public n_r_eq
nr. of normal points in equilibrium grid
Definition grid_vars.f90:20
grid_vars::max_par_x
real(dp), public max_par_x
max. of parallel coordinate [ ] in field-aligned grid
Definition grid_vars.f90:25
grid_vars::n_r_in
integer, public n_r_in
nr. of normal points in input grid
Definition grid_vars.f90:19
grid_vars::n_r_sol
integer, public n_r_sol
nr. of normal points in solution grid
Definition grid_vars.f90:22
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition HDF5_ops.f90:27
hdf5_ops::print_hdf5_grid
integer function, public print_hdf5_grid(grid_id, grid_name, grid_type, grid_time, grid_top, grid_geom, grid_atts, grid_grids, reset, ind_plot)
Prints an HDF5 grid.
Definition HDF5_ops.f90:886
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_ops::add_hdf5_item
integer function, public add_hdf5_item(file_info, xdmf_item, reset, ind_plot)
Add an XDMF item to a HDF5 file.
Definition HDF5_ops.f90:340
hdf5_ops::open_hdf5_file
integer function, public open_hdf5_file(file_info, file_name, sym_type, descr, ind_plot, cont_plot)
Opens an HDF5 file and accompanying xmf file and returns the handles.
Definition HDF5_ops.f90:78
hdf5_ops::print_hdf5_geom
subroutine, public print_hdf5_geom(geom_id, geom_type, geom_dataitems, reset, ind_plot)
Prints an HDF5 geometry.
Definition HDF5_ops.f90:805
hdf5_ops::close_hdf5_file
integer function, public close_hdf5_file(file_info, ind_plot, cont_plot)
Closes an HDF5 file and writes the accompanying xmf file.
Definition HDF5_ops.f90:264
hdf5_ops::print_hdf5_3d_data_item
integer function, public print_hdf5_3d_data_item(dataitem_id, file_info, var_name, var, dim_tot, loc_dim, loc_offset, init_val, ind_plot, cont_plot)
Prints an HDF5 data item.
Definition HDF5_ops.f90:396
hdf5_ops::print_hdf5_top
subroutine, public print_hdf5_top(top_id, top_type, top_n_elem, ind_plot)
Prints an HDF5 topology.
Definition HDF5_ops.f90:755
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
helena_vars
Variables that have to do with HELENA quantities.
Definition HELENA_vars.f90:4
helena_vars::nchi
integer, public nchi
nr. of poloidal points
Definition HELENA_vars.f90:35
helena_vars::chi_h
real(dp), dimension(:), allocatable, public chi_h
poloidal angle
Definition HELENA_vars.f90:23
helena_vars::flux_p_h
real(dp), dimension(:,:), allocatable, public flux_p_h
poloidal flux
Definition HELENA_vars.f90:24
helena_vars::flux_t_h
real(dp), dimension(:,:), allocatable, public flux_t_h
toroidal flux
Definition HELENA_vars.f90:25
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_utilities
Numerical utilities.
Definition num_utilities.f90:4
num_vars
Numerical variables used by most other modules.
Definition num_vars.f90:4
num_vars::sol_n_procs
integer, public sol_n_procs
nr. of MPI processes for solution with SLEPC
Definition num_vars.f90:70
num_vars::dp
integer, parameter, public dp
double precision
Definition num_vars.f90:46
num_vars::n_theta_plot
integer, public n_theta_plot
nr. of poloidal points in plot
Definition num_vars.f90:156
num_vars::norm_disc_prec_sol
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition num_vars.f90:122
num_vars::max_r_plot
real(dp), public max_r_plot
max. of r_plot
Definition num_vars.f90:164
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::max_njq_change
real(dp), public max_njq_change
maximum change of prim. mode number times saf. fac. / rot. transf. when using X_style 2 (fast)
Definition num_vars.f90:119
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition num_vars.f90:69
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
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::min_r_plot
real(dp), public min_r_plot
min. of r_plot
Definition num_vars.f90:163
num_vars::rank
integer, public rank
MPI rank.
Definition num_vars.f90:68
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::max_theta_plot
real(dp), public max_theta_plot
max. of theta_plot
Definition num_vars.f90:160
num_vars::use_pol_flux_e
logical, public use_pol_flux_e
whether poloidal flux is used in E coords.
Definition num_vars.f90:113
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::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition num_vars.f90:133
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::plot_diff_hdf5
subroutine, public plot_diff_hdf5(a, b, file_name, tot_dim, loc_offset, descr, output_message)
Takes two input vectors and plots these as well as the relative and absolute difference in a HDF5 fil...
Definition output_ops.f90:1765
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
rich_vars::n_par_x
integer, public n_par_x
nr. of parallel points in field-aligned grid
Definition rich_vars.f90:20
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::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
vmec_vars
Variables that concern the output of VMEC.
Definition VMEC_vars.f90:4
vmec_vars::flux_t_v
real(dp), dimension(:,:), allocatable, public flux_t_v
toroidal flux
Definition VMEC_vars.f90:34
vmec_vars::flux_p_v
real(dp), dimension(:,:), allocatable, public flux_p_v
poloidal flux
Definition VMEC_vars.f90:35
x_vars
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
x_vars::max_r_sol
real(dp), public max_r_sol
max. normal range for pert.
Definition X_vars.f90:136
x_vars::min_r_sol
real(dp), public min_r_sol
min. normal range for pert.
Definition X_vars.f90:135
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
grid_vars::grid_type
Type for grids.
Definition grid_vars.f90:59
hdf5_vars::hdf5_file_type
HDF5 data type, containing the information about HDF5 files.
Definition HDF5_vars.f90:40
hdf5_vars::var_1d_type
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition HDF5_vars.f90:48
hdf5_vars::xml_str_type
XML strings used in XDMF.
Definition HDF5_vars.f90:33
this documentation was created with Doxygen 1.16.1 by Toon Weyens