PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
sol_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Operations on the solution vectors such as decompososing the energy, etc...
3!------------------------------------------------------------------------------!
4module sol_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, pi, rank
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, modes_type
14 use sol_vars, only: sol_type
15 use vac_vars, only: vac_type
16
17 implicit none
18 private
19 public plot_harmonics, plot_sol_vec, decompose_energy, print_output_sol, &
20 &plot_sol_vals
21#if ldebug
22 public debug_plot_sol_vec, debug_calc_e, debug_du, debug_x_norm
23#endif
24
25 ! global variables
26#if ldebug
27 logical :: debug_plot_sol_vec = .false. !< plot debug information for plot_sol_vec() \ldebug
28 logical :: debug_plot_harmonics = .false. !< plot debug information for plot_harmonics() \ldebug
29 logical :: debug_calc_e = .false. !< plot debug information for calc_E() \ldebug
30 logical :: debug_x_norm = .false. !< plot debug information \c X_norm \ldebug
31 logical :: debug_du = .false. !< plot debug information for calculation of \c DU \ldebug
32#endif
33
34contains
35 !> Plots Eigenvalues.
36 subroutine plot_sol_vals(sol,last_unstable_id)
37 use num_vars, only: rank
38
39 ! input / output
40 type(sol_type), intent(in) :: sol !< solution variables
41 integer, intent(in) :: last_unstable_id !< index of last unstable EV
42
43 ! local variables
44 character(len=max_str_ln) :: plot_title ! title for plots
45 character(len=max_str_ln) :: plot_name ! file name for plots
46 integer :: n_sol_found ! how many solutions found and saved
47 integer :: id ! counter
48
49 ! set local variables
50 n_sol_found = size(sol%val)
51
52 ! only let master plot
53 if (rank.eq.0) then
54 ! Last Eigenvalues
55 plot_title = 'final Eigenvalues omega^2 [log]'
56 plot_name = 'Eigenvalues'
57 call print_ex_2d(plot_title,plot_name,&
58 &log10(abs(rp(sol%val(1:n_sol_found)))),draw=.false.)
59 call draw_ex([plot_title],plot_name,1,1,.false.,ex_plot_style=1)
60
61 ! Last Eigenvalues: unstable range
62 if (last_unstable_id.gt.0) then
63 plot_title = 'final unstable Eigenvalues omega^2'
64 plot_name = 'Eigenvalues_unstable'
65 call print_ex_2d(plot_title,plot_name,&
66 &rp(sol%val(1:last_unstable_id)),&
67 &x=[(id*1._dp,id=1,last_unstable_id)],draw=.false.)
68 call draw_ex([plot_title],plot_name,1,1,.false.,ex_plot_style=1)
69 end if
70
71 ! Last Eigenvalues: stable range
72 if (last_unstable_id.lt.n_sol_found) then
73 plot_title = 'final stable Eigenvalues omega^2'
74 plot_name = 'Eigenvalues_stable'
75 call print_ex_2d(plot_title,plot_name,&
76 &rp(sol%val(last_unstable_id+1:n_sol_found)),&
77 &x=[(id*1._dp,id=last_unstable_id+1,n_sol_found)],&
78 &draw=.false.)
79 call draw_ex([plot_title],plot_name,1,1,.false.,ex_plot_style=1)
80 end if
81 end if
82 end subroutine plot_sol_vals
83
84 !> Plots Eigenvectors.
85 !!
86 !! This is done using the angular part of the the provided equilibrium grid
87 !! and the normal part of the provided solution grid.
88 !!
89 !! The perturbation grid is assumed to have the same angular coordinates as
90 !! the equilibrium grid, and the normal coordinates correspond to either the
91 !! equilibrium grid (\c X_grid_style 1), the solution grid (\c X_grid_style
92 !! 2) or the enriched equilibrium grid (\c X_grid_style 3).
93 !!
94 !! The output grid, furthermore, has the angular part corresponding to the
95 !! equilibrium grid, and the normal part to the solution grid.
96 !!
97 !! \note The normalization factors are taken into account and the output is
98 !! transformed back to unnormalized values:
99 !! \f[\begin{aligned}
100 !! \vec{\xi} &\sim \frac{X_0}{R_0 B_0} \ , \\
101 !! \vec{Q} &\sim \frac{X_0}{R_0^2} \ , \\
102 !! \end{aligned}\f]
103 !!
104 !! which translates to
105 !! \f[\begin{aligned}
106 !! X &\sim X_0 \ , \\
107 !! U &\sim \frac{X_0}{R_0^2 B_0} \ , \\
108 !! Qn &\sim \frac{X_0 B_0}{R_0} \ , \\
109 !! Qg &\sim \frac{X_0}{R_0^3} \ ,
110 !! \end{aligned}\f]
111 !! where \f$X_0\f$ is not determined but is common to all factors.
112 !!
113 !! \return ierr
114 integer function plot_sol_vec(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,&
115 &eq_2,X,sol,XYZ,X_id,plot_var) result(ierr)
116
117 use num_vars, only: no_plots, tol_zero, pert_mult_factor_post, &
118 &eq_job_nr, eq_jobs_lims, eq_job_nr, norm_disc_prec_eq, &
119 &norm_disc_prec_x, norm_disc_prec_sol, use_normalization, &
120 &x_grid_style, eq_style
121 use grid_utilities, only: trim_grid, calc_vec_comp
122 use sol_utilities, only: calc_xuq
123 use eq_vars, only: r_0, b_0
124 use num_utilities, only: c, spline
125 use mpi_utilities, only: get_ser_var
126#if ldebug
127 use num_vars, only: use_pol_flux_f
128#endif
129
130 character(*), parameter :: rout_name = 'plot_sol_vec'
131
132 ! input / output
133 type(modes_type), intent(in) :: mds_x !< general modes variables in perturbation grid
134 type(modes_type), intent(in) :: mds_sol !< general modes variables in solution grid
135 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
136 type(grid_type), intent(in) :: grid_x !< perturbation grid
137 type(grid_type), intent(in) :: grid_sol !< solution grid
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(in) :: x !< perturbation variables
141 type(sol_type), intent(in) :: sol !< solution variables
142 real(dp), intent(in) :: xyz(:,:,:,:) !< X, Y and Z of extended eq_grid
143 integer, intent(in) :: x_id !< nr. of Eigenvalue (for output name)
144 logical, intent(in) :: plot_var(2) !< whether variables are plotted (1: plasma perturbation, 2: magnetic perturbation)
145
146 ! local variables
147 type(grid_type) :: grid_out ! output grid (see description)
148 type(grid_type) :: grid_out_trim ! trimmed output grid
149 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
150 integer :: id, jd, jd2, kd, td ! counters
151 integer :: n_t(2) ! nr. of time steps in quarter period, nr. of quarter periods
152 integer :: plot_dim(4) ! dimensions of plot
153 integer :: plot_offset(4) ! local offset of plot
154 integer :: col ! collection type for HDF5 plot
155 logical :: cont_plot ! continued plot
156 logical :: plot_var_loc(2) ! local plot_var
157 logical :: xyz_plot_setup ! XYZ_plot_setup has been set up
158 real(dp) :: x_0 ! X at theta = zeta = 0
159 real(dp), allocatable :: x_0_ser(:) ! X_0 for all processes
160 real(dp), allocatable :: time(:) ! fraction of Alfvén time
161 real(dp), allocatable :: xyz_plot(:,:,:,:,:) ! copies of XYZ
162 real(dp), allocatable :: xyz_vec(:,:,:,:,:) ! copies of XYZ for vector plot
163 real(dp), allocatable :: h12(:,:,:) ! interpolated h_FD(1,2) on solution grid
164 real(dp), allocatable :: h22(:,:,:) ! interpolated h_FD(2,2) on solution grid
165 real(dp), allocatable :: h23(:,:,:) ! interpolated h_FD(3,2) on solution grid
166 real(dp), allocatable :: g13(:,:,:) ! interpolated g_FD(1,3) on solution grid
167 real(dp), allocatable :: g33(:,:,:) ! interpolated g_FD(3,3) on solution grid
168 real(dp), allocatable :: jac_fd_int(:,:,:) ! interpolated jac_FD on solution grid
169 real(dp), allocatable :: f_plot_phase(:,:,:,:) ! phase of f_plot
170 real(dp), allocatable :: ccomp(:,:,:,:,:) ! Cart. components of perturbation
171 complex(dp) :: omega ! sqrt of Eigenvalue
172 complex(dp), allocatable :: f_plot(:,:,:,:,:) ! the function to plot
173 character(len=max_str_ln) :: err_msg ! error message
174 character(len=max_str_ln) :: var_name(2) ! name of variable that is plot
175 character(len=max_str_ln) :: file_name(2) ! name of file
176 character(len=max_str_ln) :: description(2) ! description
177 character(len=2) :: sol_name ! name of solution vector ('xi' or 'Q')
178#if ldebug
179 integer :: nm ! n (pol. flux) or m (tor. flux)
180 integer :: ld ! counter
181 real(dp), allocatable :: dq_saf(:) ! norm. deriv. of q_saf interpolated on solution grid
182 real(dp), allocatable :: drot_t(:) ! norm. deriv. of rot_t interpolated on solution grid
183 complex(dp), allocatable :: u_inf(:,:,:,:) ! ideal ballooning U
184 complex(dp), allocatable :: u_inf_prop(:,:,:) ! proportional part of U_inf
185#endif
186
187 ! initialize ierr
188 ierr = 0
189
190 ! bypass plots if no_plots
191 if (no_plots) return
192
193 ! return if no variables are plot
194 if (count(plot_var).eq.0) return
195 plot_var_loc = plot_var
196 if (abs(pert_mult_factor_post).ge.tol_zero) then
197 ! set up X_0
198 x_0 = 0._dp
199 do kd = 1,grid_sol%loc_n_r
200 x_0 = max(x_0,abs(sum(sol%vec(:,kd,x_id))))
201 end do
202 ierr = get_ser_var([x_0],x_0_ser,scatter=.true.)
203 chckerr('')
204 x_0 = maxval(x_0_ser)
205
206 ! user output
207 call writo('Perturbing the position vector (X,Y,Z) with &
208 &multiplicative factor '//trim(r2strt(pert_mult_factor_post)))
209 if (.not.plot_var(1)) then
210 call writo('For nonzero "pert_mult_factor_POST", need to also &
211 &plot xi',warning=.true.)
212 plot_var_loc(1) = .true.
213 end if
214 end if
215
216 ! user output
217 call writo('Plot the solution vector')
218 call lvl_ud(1)
219
220#if ldebug
221 ! tests
222 if (size(xyz,4).ne.3) then
223 ierr = 1
224 err_msg = 'X, Y and Z needed'
225 chckerr(err_msg)
226 end if
227 if (grid_eq%n(1).ne.size(xyz,1) .or. &
228 &grid_eq%n(2).ne.size(xyz,2) .or. &
229 &grid_sol%loc_n_r.ne.size(xyz,3)) then
230 ierr = 1
231 err_msg = 'XYZ needs to have the correct dimensions'
232 chckerr(err_msg)
233 end if
234#endif
235
236 ! set up n_t
237 ! if the Eigenvalue is negative, the Eigenfunction explodes, so limit
238 ! n_t(2) to 1. If it is positive, the Eigenfunction oscilates, so choose
239 ! n_t(2) = 8 for 2 whole periods
240 if (rp(sol%val(x_id)).lt.0) then ! exploding, unstable
241 n_t(1) = 1 ! 1 point per quarter period
242 n_t(2) = 1 ! 1 quarter period
243 else ! oscillating, stable
244 n_t(1) = 2 ! 2 points per quarter period
245 n_t(2) = 4 ! 4 quarter periods
246 end if
247
248 ! set collection type
249 if (product(n_t).gt.1) then
250 col = 2 ! temporal collection
251 else
252 col = 1 ! no collection
253 end if
254
255 ! set up output grid
256 ierr = grid_out%init([grid_eq%n(1:2),grid_sol%n(3)],&
257 &i_lim=[grid_sol%i_min,grid_sol%i_max],divided=grid_sol%divided)
258 chckerr('')
259 grid_out%r_F = grid_sol%r_F
260 grid_out%r_E = grid_sol%r_E
261 grid_out%loc_r_F = grid_sol%loc_r_F
262 grid_out%loc_r_E = grid_sol%loc_r_E
263 if (eq_style.eq.1) allocate(grid_out%trigon_factors(&
264 &size(grid_x%trigon_factors,1),size(grid_x%trigon_factors,2),&
265 &size(grid_x%trigon_factors,3),grid_out%loc_n_r,&
266 &size(grid_x%trigon_factors,5))) ! VMEC
267 select case (x_grid_style)
268 case (1,3) ! equilibrium or enriched
269 ! interpolate X->sol
270 do jd = 1,grid_x%n(2)
271 do id = 1,grid_x%n(1)
272 ierr = spline(grid_x%loc_r_F,grid_x%theta_F(id,jd,:),&
273 &grid_sol%loc_r_F,grid_out%theta_F(id,jd,:),&
274 &ord=norm_disc_prec_x)
275 chckerr('')
276 ierr = spline(grid_x%loc_r_F,grid_x%theta_E(id,jd,:),&
277 &grid_sol%loc_r_F,grid_out%theta_E(id,jd,:),&
278 &ord=norm_disc_prec_x)
279 chckerr('')
280 ierr = spline(grid_x%loc_r_F,grid_x%zeta_F(id,jd,:),&
281 &grid_sol%loc_r_F,grid_out%zeta_F(id,jd,:),&
282 &ord=norm_disc_prec_x)
283 chckerr('')
284 ierr = spline(grid_x%loc_r_F,grid_x%zeta_E(id,jd,:),&
285 &grid_sol%loc_r_F,grid_out%zeta_E(id,jd,:),&
286 &ord=norm_disc_prec_x)
287 chckerr('')
288 if (eq_style.eq.1) then ! VMEC
289 do td = 1,size(grid_x%trigon_factors,5)
290 do kd = 1,size(grid_x%trigon_factors,1)
291 ierr = spline(grid_x%loc_r_F,&
292 &grid_x%trigon_factors(kd,id,jd,:,td),&
293 &grid_sol%loc_r_F,&
294 &grid_out%trigon_factors&
295 &(kd,id,jd,:,td),ord=norm_disc_prec_x)
296 chckerr('')
297 end do
298 end do
299 end if
300 end do
301 end do
302 case (2) ! solution
303 ! copy
304 grid_out%theta_F = grid_x%theta_F
305 grid_out%theta_E = grid_x%theta_E
306 grid_out%zeta_F = grid_x%zeta_F
307 grid_out%zeta_E = grid_x%zeta_E
308 end select
309
310 ! trim output grid
311 ierr = trim_grid(grid_out,grid_out_trim,norm_id)
312 chckerr('')
313
314 ! set up plot dimensions and offset
315 plot_dim = [grid_out_trim%n,product(n_t)]
316 plot_offset = [0,0,grid_out_trim%i_min-1,0]
317
318 ! possibly modify if multiple equilibrium parallel jobs
319 if (size(eq_jobs_lims,2).gt.1) then
320 plot_dim(1) = eq_jobs_lims(2,size(eq_jobs_lims,2)) - &
321 &eq_jobs_lims(1,1) + 1
322 plot_offset(1) = eq_jobs_lims(1,eq_job_nr) - 1
323 end if
324
325 ! set up continued plot
326 cont_plot = eq_job_nr.gt.1
327
328 ! For each time step, calculate the time (as fraction of Alfvén time)
329 allocate(time(product(n_t)))
330 do td = 1,product(n_t)
331 if (n_t(1).eq.1) then
332 time(td) = (td-1) * 0.25
333 else
334 time(td) = (mod(td-1,n_t(1))*1._dp/n_t(1) + &
335 &(td-1)/n_t(1)) * 0.25
336 end if
337 end do
338
339 ! set up interpolated metric factors (not trimmed)
340 allocate(h12(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
341 allocate(h22(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
342 allocate(h23(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
343 allocate(g13(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
344 allocate(g33(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
345 allocate(jac_fd_int(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r))
346#if ldebug
347 allocate(dq_saf(grid_out%loc_n_r))
348 allocate(drot_t(grid_out%loc_n_r))
349#endif
350
351 ! interpolate eq->sol
352 do jd = 1,grid_eq%n(2)
353 do id = 1,grid_eq%n(1)
354 ierr = spline(grid_eq%loc_r_F,&
355 &eq_2%h_FD(id,jd,:,c([1,2],.true.),0,0,0),grid_sol%loc_r_F,&
356 &h12(id,jd,:),ord=norm_disc_prec_eq)
357 chckerr('')
358 ierr = spline(grid_eq%loc_r_F,&
359 &eq_2%h_FD(id,jd,:,c([2,2],.true.),0,0,0),grid_sol%loc_r_F,&
360 &h22(id,jd,:),ord=norm_disc_prec_eq)
361 chckerr('')
362 ierr = spline(grid_eq%loc_r_F,&
363 &eq_2%h_FD(id,jd,:,c([2,3],.true.),0,0,0),grid_sol%loc_r_F,&
364 &h23(id,jd,:),ord=norm_disc_prec_eq)
365 chckerr('')
366 ierr = spline(grid_eq%loc_r_F,&
367 &eq_2%g_FD(id,jd,:,c([1,3],.true.),0,0,0),grid_sol%loc_r_F,&
368 &g13(id,jd,:),ord=norm_disc_prec_eq)
369 chckerr('')
370 ierr = spline(grid_eq%loc_r_F,&
371 &eq_2%g_FD(id,jd,:,c([3,3],.true.),0,0,0),grid_sol%loc_r_F,&
372 &g33(id,jd,:),ord=norm_disc_prec_eq)
373 chckerr('')
374 ierr = spline(grid_eq%loc_r_F,&
375 &eq_2%jac_FD(id,jd,:,0,0,0),grid_sol%loc_r_F,&
376 &jac_fd_int(id,jd,:),ord=norm_disc_prec_eq)
377 chckerr('')
378 end do
379 end do
380#if ldebug
381 ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,1),grid_sol%loc_r_F,&
382 &dq_saf,ord=norm_disc_prec_eq)
383 chckerr('')
384 ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,1),grid_sol%loc_r_F,&
385 &drot_t,ord=norm_disc_prec_eq)
386 chckerr('')
387#endif
388
389 ! calculate omega = sqrt(sol_val) and make sure to select the decaying
390 ! solution
391 omega = sqrt(sol%val(x_id))
392 if (ip(omega).gt.0) omega = - omega ! exploding solution, not the decaying one
393
394 ! set up f_plot (not trimmed) and XYZ_plot (trimmed)
395 ! set up:
396 ! * f_plot (not trimmed): X, U, Q
397 ! * ccomp (not trimmed): Cartesian components of perturbation
398 ! * XYZ_plot (trimmed): copies of XYZ for plot for different times
399 ! * XYZ_vec (trimmed): copies of XYZ for vector plot
400 allocate(f_plot(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r,&
401 &product(n_t),2))
402 allocate(ccomp(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r,3,2))
403 allocate(xyz_plot(grid_out%n(1),grid_out%n(2),grid_out_trim%loc_n_r,&
404 &size(time),3))
405 allocate(xyz_vec(grid_out%n(1),grid_out%n(2),grid_out_trim%loc_n_r,3,&
406 &3))
407 xyz_plot_setup = .false.
408
409 ! operations for plasma perturbation and magnetic perturbation
410 do jd = 1,2 ! first plasma perturbation, then magnetic perturbation
411 if (.not.plot_var_loc(jd)) cycle
412
413 ! calculate vector components and set solution name, variable name,
414 ! file name and description
415 select case (jd)
416 case (1) ! plasma perturbation
417 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,&
418 &eq_2,x,sol,x_id,1,time,f_plot(:,:,:,:,1))
419 chckerr('')
420 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,&
421 &eq_2,x,sol,x_id,2,time,f_plot(:,:,:,:,2))
422 chckerr('')
423 sol_name = 'xi'
424 file_name(1) = trim(i2str(x_id))//'_sol_X'
425 file_name(2) = trim(i2str(x_id))//'_sol_U'
426 case (2) ! magnetic perturbation
427 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,&
428 &eq_2,x,sol,x_id,3,time,f_plot(:,:,:,:,1))
429 chckerr('')
430 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,&
431 &eq_2,x,sol,x_id,4,time,f_plot(:,:,:,:,2))
432 chckerr('')
433 sol_name = 'Q'
434 file_name(1) = trim(i2str(x_id))//'_sol_Qn'
435 file_name(2) = trim(i2str(x_id))//'_sol_Qg'
436 end select
437 var_name(1) = 'Normal comp. of '//trim(sol_name)
438 var_name(2) = 'Geodesic comp. of '//trim(sol_name)
439 description(1) = 'Normal component of vector '//trim(sol_name)//&
440 &' for Eigenvalue '//trim(i2str(x_id))//' with omega = '//&
441 &trim(r2str(rp(omega)))
442 description(2) = 'Geodesic component of vector '//trim(sol_name)//&
443 &' for Eigenvalue '//trim(i2str(x_id))//' with omega = '//&
444 &trim(r2str(rp(omega)))
445
446 ! calculate vector components at every time point
447 do td = 1,product(n_t)
448 ! covariant components:
449 ! xi . e_alpha = U J^2 h22/g33
450 ! xi . e_psi = X 1/h22 - U J^2 h12/g33
451 ! xi . e_theta = 0
452 ! and similar for Q
453 ccomp(:,:,:,1,1) = rp(f_plot(:,:,:,td,2)) * &
454 &jac_fd_int**2*h22/g33
455 ccomp(:,:,:,2,1) = rp(f_plot(:,:,:,td,1)) * &
456 &1._dp/h22 - &
457 &rp(f_plot(:,:,:,td,2)) * &
458 &jac_fd_int**2*h12/g33
459 ccomp(:,:,:,3,1) = 0._dp
460
461 ! contravariant components:
462 ! xi . nabla alpha = X h12/h22 + U
463 ! xi . nabla psi = X
464 ! xi . nabla theta = X h23/h22 - U g13/g33
465 ! and similar for Q
466 ccomp(:,:,:,1,2) = rp(f_plot(:,:,:,td,1)) * &
467 &h12/h22 + rp(f_plot(:,:,:,td,2))
468 ccomp(:,:,:,2,2) = rp(f_plot(:,:,:,td,1))
469 ccomp(:,:,:,3,2) = rp(f_plot(:,:,:,td,1)) * &
470 &h23/h22 - rp(f_plot(:,:,:,td,2)) * g13/g33
471
472 ! multiplication factor if xi
473 if (abs(pert_mult_factor_post).ge.tol_zero .and. jd.eq.1) &
474 &ccomp = ccomp*pert_mult_factor_post/x_0
475
476 ! transform normalized quantities to unnormalized
477 if (use_normalization) then
478 select case (jd)
479 case (1) ! plasma perturbation
480 ccomp = ccomp / (r_0*b_0) ! xi ~ X_0/(R_0*B_0)
481 case (2) ! magnetic perturbation
482 ccomp = ccomp / (r_0**2) ! Q ~ X_0/(R_0^2)
483 end select
484 end if
485
486 ! transform to Cartesian components
487 ierr = calc_vec_comp(grid_out,grid_eq,eq_1,eq_2,ccomp,&
488 &norm_disc_prec_sol)
489 chckerr('')
490
491 ! vector plot
492 do id = 1,3
493 xyz_vec(:,:,:,id,:) = xyz(:,:,norm_id(1):norm_id(2),:)
494 end do
495 call plot_hdf5([trim(sol_name)],trim(sol_name)//'_T_'//&
496 &trim(i2str(td)),ccomp(:,:,norm_id(1):norm_id(2),:,1),&
497 &tot_dim=[plot_dim(1:3),3],&
498 &loc_offset=[plot_offset(1:3),0],&
499 &x=xyz_vec(:,:,:,:,1),y=xyz_vec(:,:,:,:,2),&
500 &z=xyz_vec(:,:,:,:,3),col=4,cont_plot=cont_plot,&
501 &descr='perturbation for time t = '//&
502 &trim(r2strt(time(td))))
503
504 ! perturb position vector if xi; it does not matter whether
505 ! we take covariant or contravariant components
506 if (.not.xyz_plot_setup) then
507 xyz_plot(:,:,:,td,:) = xyz(:,:,norm_id(1):norm_id(2),:)
508 if (abs(pert_mult_factor_post).ge.tol_zero .and. jd.eq.1) &
509 &xyz_plot(:,:,:,td,:) = xyz_plot(:,:,:,td,:) + &
510 &ccomp(:,:,norm_id(1):norm_id(2),:,1)
511 end if
512 end do
513
514 ! XYZ_plot has been set up if we get here
515 xyz_plot_setup = .true.
516
517#if ldebug
518 if (debug_plot_sol_vec .and. jd.eq.1) then
519 call writo('Checking how well U approximates the pure &
520 &ballooning result')
521 call lvl_ud(1)
522
523 ! set up ideal U
524 allocate(u_inf(grid_out%n(1),grid_out%n(2),grid_out%loc_n_r,&
525 &product(n_t)))
526
527 ! derive the X vector
528 do ld = 1,product(n_t)
529 do jd2 = 1,grid_out%n(2)
530 do id = 1,grid_out%n(1)
531 ierr = spline(grid_out%loc_r_F,&
532 &f_plot(id,jd2,:,ld,1),grid_out%loc_r_F,&
533 &u_inf(id,jd2,:,ld),ord=norm_disc_prec_eq,&
534 &deriv=1)
535 chckerr('')
536 end do
537 end do
538 end do
539
540 ! set up dummy variable Theta^alpha + q' theta and nm
541 allocate(u_inf_prop(grid_out%n(1),grid_out%n(2),&
542 grid_out%loc_n_r))
543 if (use_pol_flux_f) then
544 do kd = 1,grid_out%loc_n_r
545 u_inf_prop(:,:,kd) = h12(:,:,kd)/h22(:,:,kd) + &
546 &dq_saf(kd)*grid_out%theta_F(:,:,kd)
547 end do
548 nm = sol%n(1,1)
549 else
550 do kd = 1,grid_out%loc_n_r
551 u_inf_prop(:,:,kd) = h12(:,:,kd)/h22(:,:,kd) - &
552 &drot_t(kd)*grid_out%zeta_F(:,:,kd)
553 end do
554 nm = sol%m(1,1)
555 end if
556
557 ! multiply by i/n or i/m and add the proportional part
558 do ld = 1,product(n_t)
559 u_inf(:,:,:,ld) = iu/nm * u_inf(:,:,:,ld) - &
560 &f_plot(:,:,:,ld,1) * u_inf_prop
561 end do
562 deallocate(u_inf_prop)
563
564 ! plotting real part
565 call plot_hdf5('RE X','TEST_RE_X_'//trim(i2str(x_id)),&
566 &rp(f_plot(:,:,norm_id(1):norm_id(2),1,1)),&
567 &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3))
568 call plot_diff_hdf5(rp(f_plot(:,:,norm_id(1):norm_id(2),1,2)),&
569 &rp(u_inf(:,:,norm_id(1):norm_id(2),1)),&
570 &'TEST_RE_U_inf_'//trim(i2str(x_id)),tot_dim=plot_dim(1:3),&
571 &loc_offset=plot_offset(1:3),descr='To test whether U &
572 &approximates the ideal ballooning result',&
573 &output_message=.true.)
574
575 ! plotting imaginary part
576 call plot_hdf5('IM X','TEST_IM_X_'//trim(i2str(x_id)),&
577 &ip(f_plot(:,:,norm_id(1):norm_id(2),1,1)),&
578 &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3))
579 call plot_diff_hdf5(ip(f_plot(:,:,norm_id(1):norm_id(2),1,2)),&
580 &ip(u_inf(:,:,norm_id(1):norm_id(2),1)),&
581 &'TEST_IM_U_inf_'//trim(i2str(x_id)),tot_dim=plot_dim(1:3),&
582 &loc_offset=plot_offset(1:3),descr='To test whether U &
583 &approximates the ideal ballooning result',&
584 &output_message=.true.)
585
586 call writo('Checking done')
587 call lvl_ud(-1)
588 end if
589#endif
590
591 ! set up temporary variable for phase
592 allocate(f_plot_phase(grid_out%n(1),grid_out%n(2),&
593 &grid_out_trim%loc_n_r,product(n_t)))
594
595 ! transform normalized quantities to unnormalized
596 if (use_normalization) then
597 select case (jd)
598 case (1) ! plasma perturbation
599 !f_plot(:,:,:,:,1) = f_plot(:,:,:,:,1) ! X
600 f_plot(:,:,:,:,2) = f_plot(:,:,:,:,2) / (r_0**2*b_0) ! U
601 case (2) ! magnetic perturbation
602 f_plot(:,:,:,:,1) = f_plot(:,:,:,:,1) * b_0/r_0 ! Qn
603 f_plot(:,:,:,:,2) = f_plot(:,:,:,:,2) / (r_0**3) ! Qg
604 end select
605 end if
606
607 do kd = 1,2
608 call plot_hdf5([var_name(kd)],trim(file_name(kd))//'_RE',&
609 &rp(f_plot(:,:,norm_id(1):norm_id(2),:,kd)),&
610 &tot_dim=plot_dim,loc_offset=plot_offset,&
611 &x=xyz_plot(:,:,:,:,1),y=xyz_plot(:,:,:,:,2),&
612 &z=xyz_plot(:,:,:,:,3),col=col,cont_plot=cont_plot,&
613 &descr=description(kd))
614 call plot_hdf5([var_name(kd)],trim(file_name(kd))//'_IM',&
615 &ip(f_plot(:,:,norm_id(1):norm_id(2),:,kd)),&
616 &tot_dim=plot_dim,loc_offset=plot_offset,&
617 &x=xyz_plot(:,:,:,:,1),y=xyz_plot(:,:,:,:,2),&
618 &z=xyz_plot(:,:,:,:,3),col=col,cont_plot=cont_plot,&
619 &descr=description(kd))
620 f_plot_phase = atan2(&
621 &ip(f_plot(:,:,norm_id(1):norm_id(2),:,kd)),&
622 &rp(f_plot(:,:,norm_id(1):norm_id(2),:,kd)))
623 where (f_plot_phase.lt.0) f_plot_phase = f_plot_phase + 2*pi
624 call plot_hdf5([var_name(kd)],trim(file_name(kd))//'_PH',&
625 &f_plot_phase,tot_dim=plot_dim,loc_offset=plot_offset,&
626 &x=xyz_plot(:,:,:,:,1),y=xyz_plot(:,:,:,:,2),&
627 &z=xyz_plot(:,:,:,:,3),col=col,cont_plot=cont_plot,&
628 &descr=description(kd))
629 end do
630
631 ! clean up
632 deallocate(f_plot_phase)
633 end do
634
635 ! clean up
636 call grid_out%dealloc()
637 call grid_out_trim%dealloc()
638
639 call lvl_ud(-1)
640 end function plot_sol_vec
641
642 !> Plots the harmonics and their maximum in 2-D.
643 !!
644 !! \return ierr
645 integer function plot_harmonics(mds,grid_sol,sol,X_id,res_surf) result(ierr)
646 use mpi_utilities, only: get_ser_var
647 use num_vars, only: ex_max_size, rank, no_plots, use_pol_flux_f
648 use eq_vars, only: max_flux_f
649 use sol_utilities, only: calc_tot_sol_vec
650 use grid_utilities, only: trim_grid
651#if ldebug
652 use num_vars, only: norm_disc_prec_sol
653#endif
654
655 character(*), parameter :: rout_name = 'plot_harmonics'
656
657 ! input / output
658 type(modes_type), intent(in) :: mds !< general modes variables in solution grid
659 type(grid_type), intent(in) :: grid_sol !< solution grid
660 type(sol_type), intent(in) :: sol !< solution variables
661 integer, intent(in) :: x_id !< nr. of Eigenvalue (for output name)
662 real(dp), intent(in) :: res_surf(:,:) !< resonant surfaces
663
664 ! local variables
665 type(grid_type) :: grid_sol_trim ! trimmed sol grid
666 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
667 integer :: id, ld ! counters
668 integer :: ld_loc ! local ld
669 integer :: max_deriv ! maximum derivative for sol_vec_ser_tot
670 integer :: n_mod_tot ! total number of modes that can resonate
671 integer :: n_mod_loc ! local number of modes that are used in the calculations
672 integer :: min_nm_x ! minimum n or m in total modes
673 character(len=max_str_ln) :: file_name ! name of file of plots of this proc.
674 character(len=max_str_ln), allocatable :: plot_title(:) ! title for plots
675 real(dp) :: norm_factor ! conversion factor max_flux/2pi from flux to normal coordinate
676 real(dp), allocatable :: x_plot(:,:) ! x values of plot
677 real(dp), allocatable :: y_plot(:,:) ! y values of plot
678 complex(dp), allocatable :: sol_vec_ser(:,:) ! serial MPI Eigenvector for local number of modes
679 complex(dp), allocatable :: sol_vec_dum(:,:) ! dummy vector for sol_vec calculations
680 complex(dp), allocatable :: sol_vec_ser_tot(:,:,:) ! serial MPI Eigenvector for total number of modes and derivatives
681 complex(dp), allocatable :: sol_vec_ser_loc(:) ! local sol_vec_ser
682 real(dp), allocatable :: sol_vec_phase(:,:) ! phase of sol_vec
683 character :: nm_x ! n or m
684 character(len=max_str_ln) :: err_msg ! error message
685#if ldebug
686 real(dp), allocatable :: sol_vec_sum(:,:) ! for test output
687#endif
688
689 ! initialize ierr
690 ierr = 0
691
692 ! bypass plots if no_plots
693 if (no_plots) return
694
695 ! trim sol grid
696 ierr = trim_grid(grid_sol,grid_sol_trim,norm_id)
697 chckerr('')
698
699 ! set up local and total nr. of modes, which can be different for X
700 ! style 2 (see discussion in sol_utilities)
701 n_mod_loc = sol%n_mod
702 n_mod_tot = size(mds%sec,1)
703 min_nm_x = minval(mds%sec(:,1),1)
704
705 ! set up serial sol_vec on master
706 if (rank.eq.0) &
707 &allocate(sol_vec_ser(1:n_mod_loc,1:grid_sol_trim%n(3)))
708 do ld = 1,n_mod_loc
709 ierr = get_ser_var(sol%vec(ld,norm_id(1):norm_id(2),x_id),&
710 &sol_vec_ser_loc,scatter=.true.)
711 chckerr('')
712 if (rank.eq.0) sol_vec_ser(ld,:) = sol_vec_ser_loc
713 deallocate(sol_vec_ser_loc)
714 end do
715
716 ! convert to full mode on master and set up normal factor
717 if (rank.eq.0) then
718#if ldebug
719 select case (norm_disc_prec_sol)
720 case (1:2)
721 max_deriv = 1
722 case (3)
723 max_deriv = 2
724 end select
725#else
726 max_deriv = 0
727#endif
728 allocate(sol_vec_ser_tot(1:n_mod_tot,1:grid_sol_trim%n(3),&
729 &0:max_deriv))
730 do id = 0,max_deriv
731 ierr = calc_tot_sol_vec(mds,grid_sol,sol_vec_ser,&
732 &sol_vec_dum,deriv=id) ! Note: need to use untrimmed grid
733 chckerr('')
734 sol_vec_ser_tot(:,:,id) = sol_vec_dum
735 end do
736 deallocate(sol_vec_dum)
737
738 norm_factor = max_flux_f/(2*pi)
739
740 ! set up x_plot, y_plot and sol_vec_phase
741 allocate(x_plot(grid_sol_trim%n(3),n_mod_tot))
742 allocate(y_plot(grid_sol_trim%n(3),n_mod_tot))
743 allocate(sol_vec_phase(grid_sol_trim%n(3),n_mod_tot))
744 if (use_pol_flux_f) then
745 nm_x = 'm'
746 else
747 nm_x = 'n'
748 end if
749 do ld = 1,n_mod_tot
750 x_plot(:,ld) = grid_sol_trim%r_F
751 y_plot(:,ld) = min_nm_x+ld-1
752 end do
753
754 ! scale x_plot by normal factor
755 x_plot = x_plot/norm_factor
756
757 do id = 0,max_deriv
758 ! 1. prepare 2-D plots
759 allocate(plot_title(n_mod_tot))
760 do ld = 1,n_mod_tot
761 plot_title(ld) = 'EV - midplane ('//nm_x//' = '//&
762 &trim(i2str(min_nm_x+ld-1))//')'
763 end do
764
765 ! 2. plot real part at midplane
766
767 ! set up file name of this rank and plot title
768 file_name = trim(i2str(x_id))//'_EV_midplane_RE'
769 if (id.gt.0) file_name = trim(file_name)//'_D'//trim(i2str(id))
770
771 ! print real amplitude of harmonics of eigenvector at midplane
772 call print_ex_2d(plot_title,file_name,&
773 &rp(transpose(sol_vec_ser_tot(:,:,id))),x=x_plot,&
774 &draw=.false.)
775 call draw_ex(plot_title,file_name,n_mod_tot,1,&
776 &.false.,draw_ops=['with lines'],ex_plot_style=1)
777 call draw_ex(plot_title,file_name,n_mod_tot,1,&
778 &.false.,ex_plot_style=2)
779
780 ! plot in file using decoupled 3D if not too big
781 if (n_mod_tot*grid_sol_trim%n(3).le.ex_max_size) then
782 call draw_ex(plot_title,trim(file_name)//'_3D',&
783 &n_mod_tot,3,.false.,draw_ops=['with lines'],&
784 &data_name=file_name,ex_plot_style=1)
785 end if
786
787 ! 3. plot imaginary part at midplane
788
789 ! set up file name of this rank and plot title
790 file_name = trim(i2str(x_id))//'_EV_midplane_IM'
791 if (id.gt.0) file_name = trim(file_name)//'_D'//trim(i2str(id))
792
793 ! print imag amplitude of harmonics of eigenvector at midplane
794 call print_ex_2d(plot_title,file_name,&
795 &ip(transpose(sol_vec_ser_tot(:,:,id))),x=x_plot,&
796 &draw=.false.)
797
798 ! plot in file
799 call draw_ex(plot_title,file_name,n_mod_tot,1,.false.,&
800 &draw_ops=['with lines'],ex_plot_style=1)
801
802 ! plot in file using decoupled 3D if not too big
803 if (n_mod_tot*grid_sol_trim%n(3).le.ex_max_size) then
804 call draw_ex(plot_title,trim(file_name)//'_3D',&
805 &n_mod_tot,3,.false.,draw_ops=['with lines'],&
806 &data_name=file_name,ex_plot_style=1)
807 end if
808
809 ! 4. plot phase at midplane
810 ! set up file name of this rank and plot title
811 file_name = trim(i2str(x_id))//'_EV_midplane_PH'
812 if (id.gt.0) file_name = trim(file_name)//'_D'//trim(i2str(id))
813
814 ! set up phase
815 sol_vec_phase = atan2(ip(transpose(sol_vec_ser_tot(:,:,id))),&
816 &rp(transpose(sol_vec_ser_tot(:,:,id))))
817 where (sol_vec_phase.lt.0) sol_vec_phase = sol_vec_phase + 2*pi
818
819 ! print imag amplitude of harmonics of eigenvector at midplane
820 call print_ex_2d(plot_title,file_name,sol_vec_phase,&
821 &x=x_plot,draw=.false.)
822
823 ! plot in file
824 call draw_ex(plot_title,file_name,n_mod_tot,1,&
825 &.false.,draw_ops=['with lines'],ex_plot_style=1)
826
827 ! plot in file using decoupled 3D if not too big
828 if (n_mod_tot*grid_sol_trim%n(3).le.ex_max_size) then
829 call draw_ex(plot_title,trim(file_name)//'_3D',&
830 &n_mod_tot,3,.false.,draw_ops=['with lines'],&
831 &data_name=file_name,ex_plot_style=1)
832 end if
833
834 ! deallocate variables
835 deallocate(plot_title)
836
837 ! 5. prepare 3_D plots
838 allocate(plot_title(3))
839 plot_title(1) = 'real part'
840 plot_title(2) = 'imaginary part'
841 plot_title(3) = 'phase'
842 file_name = trim(i2str(x_id))//'_EV_midplane'
843 if (id.gt.0) file_name = trim(file_name)//'_D'//trim(i2str(id))
844
845 ! 6. plot above in HDf5
846
847 ! plot
848 call plot_hdf5(plot_title,trim(file_name),&
849 &reshape([rp(transpose(sol_vec_ser_tot(:,:,id))),&
850 &ip(transpose(sol_vec_ser_tot(:,:,id))),&
851 &sol_vec_phase],[grid_sol_trim%n(3),n_mod_tot,1,3]),&
852 &x=reshape(x_plot,[grid_sol_trim%n(3),n_mod_tot,1,1]),&
853 &y=reshape(y_plot,[grid_sol_trim%n(3),n_mod_tot,1,1]))
854
855 ! deallocate variables
856 deallocate(plot_title)
857 end do
858
859#if ldebug
860 if (debug_plot_harmonics) then
861 allocate(plot_title(4))
862 allocate(sol_vec_sum(1:grid_sol_trim%n(3),4))
863 plot_title(1) = 'EV - outboard midplane - REAL'
864 plot_title(2) = 'EV - outboard midplane - IMAG'
865 plot_title(3) = 'EV - inboard midplane - REAL'
866 plot_title(4) = 'EV - inboard midplane - IMAG'
867 do id = 0,max_deriv
868 file_name = 'TEST_harmonics'
869 if (id.gt.0) file_name = trim(file_name)//'_D'//&
870 &trim(i2str(id))
871 sol_vec_sum(:,1) = rp(sum(sol_vec_ser_tot(:,:,id),1))
872 sol_vec_sum(:,2) = ip(sum(sol_vec_ser_tot(:,:,id),1))
873 sol_vec_sum(:,3) = rp(&
874 &sum(sol_vec_ser_tot(1:n_mod_tot:2,:,id),1) - &
875 &sum(sol_vec_ser_tot(2:n_mod_tot:2,:,id),1))
876 sol_vec_sum(:,4) = ip(&
877 &sum(sol_vec_ser_tot(1:n_mod_tot:2,:,id),1) - &
878 &sum(sol_vec_ser_tot(2:n_mod_tot:2,:,id),1))
879 call print_ex_2d(plot_title,file_name,sol_vec_sum,&
880 &x=x_plot(:,1:1),draw=.false.)
881 call draw_ex(plot_title,file_name,4,1,&
882 &.false.,draw_ops=['with lines'],ex_plot_style=1)
883 call draw_ex(plot_title,file_name,4,1,&
884 &.false.,ex_plot_style=2)
885 end do
886 deallocate(plot_title)
887 deallocate(sol_vec_sum)
888 end if
889#endif
890
891 ! deallocate variables
892 deallocate(x_plot)
893 deallocate(y_plot)
894 deallocate(sol_vec_phase)
895
896 ! only plot maximum if resonant surfaces found
897 if (size(res_surf,1).gt.0) then
898 ! set up file name of this rank and plot title
899 file_name = trim(i2str(x_id))//'_EV_max'
900 allocate(plot_title(2))
901 plot_title(1) = 'mode maximum'
902 plot_title(2) = 'resonating surface'
903
904 ! set up plot variables
905 allocate(x_plot(n_mod_tot,2))
906 allocate(y_plot(n_mod_tot,2))
907
908 ! set up maximum of each mode at midplane
909 ld_loc = 0
910 do ld = 1,n_mod_tot
911 if (maxval(abs(rp(sol_vec_ser_tot(ld,:,0)))).gt.0) then ! mode resonates somewhere
912 ld_loc = ld_loc + 1
913 x_plot(ld_loc,1) = grid_sol_trim%r_F(&
914 &maxloc(abs(rp(sol_vec_ser_tot(ld,:,0))),1))
915 y_plot(ld_loc,1) = min_nm_x+ld-1
916 end if
917 end do
918 if (ld_loc.eq.0) then
919 ierr = 1
920 err_msg = 'None of the modes resonates'
921 chckerr(err_msg)
922 end if
923 x_plot(ld_loc+1:n_mod_tot,1) = x_plot(ld_loc,1) ! identical copies of last point for numerical reasons
924 y_plot(ld_loc+1:n_mod_tot,1) = y_plot(ld_loc,1) ! identical copies of last point for numerical reasons
925
926 ! set up resonating surface for each mode
927 x_plot(1:size(res_surf,1),2) = res_surf(:,2) ! normal position of resonating mode
928 x_plot(size(res_surf,1)+1:n_mod_tot,2) = &
929 &res_surf(size(res_surf,1),2) ! identical copies of last point for numerical reasons
930 y_plot(1:size(res_surf,1),2) = res_surf(:,1) ! total mode index of resonating mode
931 y_plot(size(res_surf,1)+1:n_mod_tot,2) = &
932 &res_surf(size(res_surf,1),1) ! identical copies of last point for numerical reasons
933
934 ! scale x_plot by normal factor
935 x_plot = x_plot/norm_factor
936
937 ! plot the maximum at midplane
938 call print_ex_2d(plot_title,file_name,y_plot,x=x_plot,&
939 &draw=.false.)
940
941 ! draw plot in file
942 call draw_ex(plot_title,file_name,2,1,.false.,extra_ops=&
943 &'set xrange ['//&
944 &trim(r2str(grid_sol_trim%r_F(1)/norm_factor))//':'//&
945 &trim(r2str(grid_sol_trim%r_F(grid_sol_trim%n(3))/&
946 &norm_factor))//']',ex_plot_style=1)
947 end if
948 end if
949
950 ! clean up
951 call grid_sol_trim%dealloc()
952 end function plot_harmonics
953
954 !> Decomposes the plasma potential and kinetic energy in its individual
955 !! terms for an individual Eigenvalue.
956 !!
957 !! Use is made of variables representing the potential and kinetic energy
958 !! \cite Weyens3D.
959 !! - \f$E_\text{pot}\f$:
960 !! - normal line bending term:
961 !! \f$\frac{1}{\mu_0} \frac{Q_n^2}{h^{22}}\f$,
962 !! - geodesic line bending term:
963 !! \f$\frac{1}{\mu_0} \mathcal{J}^2 \frac{h^{22}}{g_{33}} Q_g^2\f$,
964 !! - normal ballooning term:
965 !! \f$-2 p' X^2 \kappa_n\f$,
966 !! - geodesic ballooning term:
967 !! \f$-2 p' X U^* \kappa_g\f$,
968 !! - normal kink term:
969 !! \f$-\sigma X^*Q_g\f$,
970 !! - geodesic kink term:
971 !! \f$\sigma U^*Q_n\f$,
972 !! - \f$E_\text{kin}\f$:
973 !! - normal kinetic term:
974 !! \f$\rho \frac{X^2}{h^{22}}\f$,
975 !! - geodesic kinetic term:
976 !! \f$\rho \mathcal{J}^2 \frac{h^{22}}{g_{33}} U^2\f$.
977 !!
978 !! The energy terms are calculated normally on the sol grid, interpolating
979 !! the quantities defined on the eq grid, and angularly in the eq grid.
980 !!
981 !! Optionally, the results can be plotted by providing X, Y and Z. By
982 !! default, they are instead written to an output file.
983 !!
984 !! Also, the fraction between potential and kinetic energy can be returned,
985 !! compared with the eigenvalue.
986 !!
987 !! \return ierr
988 integer function decompose_energy(mds_X,mds_sol,grid_eq,grid_X,grid_sol,&
989 &eq_1,eq_2,X,sol,vac,X_id,B_aligned,XYZ,E_pot_int,E_kin_int) &
990 &result(ierr)
991
992 use grid_utilities, only: trim_grid
993 use num_vars, only: no_plots, eq_job_nr, eq_jobs_lims, eq_job_nr
994
995 character(*), parameter :: rout_name = 'decompose_energy'
996
997 ! input / output
998 type(modes_type), intent(in) :: mds_x !< general modes variables in pertubation grid
999 type(modes_type), intent(in) :: mds_sol !< general modes variables in solution grid
1000 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
1001 type(grid_type), intent(in) :: grid_x !< perturbation grid
1002 type(grid_type), intent(in) :: grid_sol !< solution grid
1003 type(eq_1_type), intent(in) :: eq_1 !< flux equilibrium
1004 type(eq_2_type), intent(in) :: eq_2 !< metric equilibrium
1005 type(x_1_type), intent(in) :: x !< perturbation variables
1006 type(sol_type), intent(in) :: sol !< solution variables
1007 type(vac_type), intent(in) :: vac !< vacuum variables
1008 integer, intent(in) :: x_id !< nr. of Eigenvalue
1009 logical, intent(in) :: b_aligned !< whether grid is field-aligned
1010 real(dp), intent(in), optional :: xyz(:,:,:,:) !< X, Y and Z for plotting
1011 complex(dp), intent(inout), optional :: e_pot_int(7) !< integrated potential energy
1012 complex(dp), intent(inout), optional :: e_kin_int(2) !< integrated kinetic energy
1013
1014 ! local variables
1015 type(grid_type) :: grid_sol_trim ! trimmed sol grid
1016 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
1017 integer :: kd ! counter
1018 integer :: loc_dim(3) ! local dimension
1019 integer :: plot_dim(3) ! total dimensions for plot
1020 integer :: plot_offset(3) ! local offsets for plot
1021 logical :: cont_plot ! continued plot
1022 complex(dp), allocatable, target :: e_pot(:,:,:,:) ! potential energy
1023 complex(dp), allocatable, target :: e_kin(:,:,:,:) ! kinetic energy
1024 complex(dp) :: e_pot_int_loc(7) ! integrated potential energy for this parallel job
1025 complex(dp) :: e_kin_int_loc(2) ! integrated kinetic energy for this parallel job
1026 complex(dp), pointer :: e_pot_trim(:,:,:,:) => null() ! trimmed part of E_pot
1027 complex(dp), pointer :: e_kin_trim(:,:,:,:) => null() ! trimmed part of E_kin
1028 real(dp), allocatable, target :: x_tot(:,:,:,:), y_tot(:,:,:,:), &
1029 &Z_tot(:,:,:,:) ! multiple copies of X, Y and Z for collection plot
1030 real(dp), pointer :: x_tot_trim(:,:,:,:) => null() ! trimmed part of X
1031 real(dp), pointer :: y_tot_trim(:,:,:,:) => null() ! trimmed part of Y
1032 real(dp), pointer :: z_tot_trim(:,:,:,:) => null() ! trimmed part of Z
1033 character(len=max_str_ln), allocatable :: var_names_pot(:) ! name of potential energy variables
1034 character(len=max_str_ln), allocatable :: var_names_kin(:) ! name of kinetic energy variables
1035 character(len=max_str_ln), allocatable :: var_names(:) ! name of other variables that are plot
1036 character(len=max_str_ln) :: file_name ! name of file
1037 character(len=max_str_ln) :: description ! description
1038
1039 ! initialize ierr
1040 ierr = 0
1041
1042 ! user output
1043 call writo('Calculate energy terms')
1044 call lvl_ud(1)
1045
1046 ! calculate for this parallel job
1047 ierr = calc_e(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,eq_2,x,sol,&
1048 &vac,b_aligned,x_id,e_pot,e_kin,e_pot_int_loc,e_kin_int_loc)
1049 chckerr('')
1050
1051 ! add to totals if requested
1052 if (present(e_pot_int)) e_pot_int = e_pot_int + e_pot_int_loc
1053 if (present(e_kin_int)) e_kin_int = e_kin_int + e_kin_int_loc
1054
1055 call lvl_ud(-1)
1056
1057 ! plot if wanted
1058 if (present(xyz)) then
1059 ! bypass plots if no_plots
1060 if (no_plots) return
1061
1062 ! user output
1063 call writo('Preparing plots')
1064 call lvl_ud(1)
1065
1066 ! set up continued plot
1067 cont_plot = eq_job_nr.gt.1
1068
1069 ! set up potential energy variable names
1070 allocate(var_names_pot(6))
1071 var_names_pot(1) = '1. normal line bending term ~ Q_n^2'
1072 var_names_pot(2) = '2. geodesic line bending term ~ Q_g^2'
1073 var_names_pot(3) = '3. normal ballooning term ~ -p'' X^2 kappa_n'
1074 var_names_pot(4) = '4. geodesic ballooning term ~ -p'' X U* kappa_g'
1075 var_names_pot(5) = '5. normal kink term ~ -sigma X*Q_g'
1076 var_names_pot(6) = '6. geodesic kink term ~ sigma U*Q_n'
1077
1078 ! set up kinetic energy variable names
1079 allocate(var_names_kin(2))
1080 var_names_kin(1) = '1. normal kinetic term ~ rho X^2'
1081 var_names_kin(2) = '2. geodesic kinetic term ~ rho U^2'
1082
1083 ! trim sol grid
1084 ierr = trim_grid(grid_sol,grid_sol_trim,norm_id)
1085 chckerr('')
1086
1087 ! set plot variables
1088 loc_dim = [grid_eq%n(1),grid_eq%n(2),grid_sol_trim%loc_n_r]
1089 plot_dim = [grid_eq%n(1),grid_eq%n(2),grid_sol_trim%n(3)]
1090 plot_offset = [0,0,grid_sol_trim%i_min-1]
1091 allocate(x_tot(loc_dim(1),loc_dim(2),loc_dim(3),6))
1092 allocate(y_tot(loc_dim(1),loc_dim(2),loc_dim(3),6))
1093 allocate(z_tot(loc_dim(1),loc_dim(2),loc_dim(3),6))
1094 do kd = 1,6
1095 x_tot(:,:,:,kd) = xyz(:,:,norm_id(1):norm_id(2),1)
1096 y_tot(:,:,:,kd) = xyz(:,:,norm_id(1):norm_id(2),2)
1097 z_tot(:,:,:,kd) = xyz(:,:,norm_id(1):norm_id(2),3)
1098 end do
1099 allocate(var_names(2))
1100
1101 ! possibly modify if multiple equilibrium parallel jobs
1102 if (size(eq_jobs_lims,2).gt.1) then
1103 plot_dim(1) = eq_jobs_lims(2,size(eq_jobs_lims,2)) - &
1104 &eq_jobs_lims(1,1) + 1
1105 plot_offset(1) = eq_jobs_lims(1,eq_job_nr) - 1
1106 end if
1107
1108 ! point to the trimmed versions of energy
1109 e_pot_trim => e_pot(:,:,norm_id(1):norm_id(2),:)
1110 e_kin_trim => e_kin(:,:,norm_id(1):norm_id(2),:)
1111
1112 call lvl_ud(-1)
1113
1114 ! user output
1115 call writo('Plot energy terms')
1116 call lvl_ud(1)
1117
1118 ! E_kin
1119 file_name = trim(i2str(x_id))//'_E_kin'
1120 description = 'Kinetic energy constituents'
1121 x_tot_trim => x_tot(:,:,:,1:2)
1122 y_tot_trim => y_tot(:,:,:,1:2)
1123 z_tot_trim => z_tot(:,:,:,1:2)
1124 call plot_hdf5(var_names_kin,trim(file_name)//'_RE',rp(e_kin_trim),&
1125 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1126 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1127 &descr=description)
1128 call plot_hdf5(var_names_kin,trim(file_name)//'_IM',ip(e_kin_trim),&
1129 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1130 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1131 &descr=description)
1132 nullify(x_tot_trim,y_tot_trim,z_tot_trim)
1133
1134 ! E_pot
1135 file_name = trim(i2str(x_id))//'_E_pot'
1136 description = 'Potential energy constituents'
1137 x_tot_trim => x_tot(:,:,:,1:6)
1138 y_tot_trim => y_tot(:,:,:,1:6)
1139 z_tot_trim => z_tot(:,:,:,1:6)
1140 call plot_hdf5(var_names_pot,trim(file_name)//'_RE',rp(e_pot_trim),&
1141 &tot_dim=[plot_dim,6],loc_offset=[plot_offset,0],&
1142 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1143 &descr=description)
1144 call plot_hdf5(var_names_pot,trim(file_name)//'_IM',ip(e_pot_trim),&
1145 &tot_dim=[plot_dim,6],loc_offset=[plot_offset,0],&
1146 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1147 &descr=description)
1148 nullify(x_tot_trim,y_tot_trim,z_tot_trim)
1149
1150 ! E_stab
1151 var_names(1) = '1. stabilizing term'
1152 var_names(2) = '2. destabilizing term'
1153 file_name = trim(i2str(x_id))//'_E_stab'
1154 description = '(de)stabilizing energy'
1155 x_tot_trim => x_tot(:,:,:,1:2)
1156 y_tot_trim => y_tot(:,:,:,1:2)
1157 z_tot_trim => z_tot(:,:,:,1:2)
1158 call plot_hdf5(var_names,trim(file_name)//'_RE',rp(reshape(&
1159 &[sum(e_pot_trim(:,:,:,1:2),4),sum(e_pot_trim(:,:,:,3:6),4)],&
1160 &[loc_dim(1),loc_dim(2),loc_dim(3),2])),&
1161 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1162 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,&
1163 &cont_plot=cont_plot,descr=description)
1164 call plot_hdf5(var_names,trim(file_name)//'_IM',ip(reshape(&
1165 &[sum(e_pot_trim(:,:,:,1:2),4),sum(e_pot_trim(:,:,:,3:6),4)],&
1166 &[loc_dim(1),loc_dim(2),loc_dim(3),2])),&
1167 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1168 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,&
1169 &cont_plot=cont_plot,descr=description)
1170
1171 ! E
1172 var_names(1) = '1. potential energy'
1173 var_names(2) = '2. kinetic energy'
1174 file_name = trim(i2str(x_id))//'_E'
1175 description = 'total potential and kinetic energy'
1176 call plot_hdf5(var_names,trim(file_name)//'_RE',rp(reshape(&
1177 &[sum(e_pot_trim,4),sum(e_kin_trim,4)],&
1178 &[loc_dim(1),loc_dim(2),loc_dim(3),2])),&
1179 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1180 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1181 &descr=description)
1182 call plot_hdf5(var_names,trim(file_name)//'_IM',ip(reshape(&
1183 &[sum(e_pot_trim,4),sum(e_kin_trim,4)],&
1184 &[loc_dim(1),loc_dim(2),loc_dim(3),2])),&
1185 &tot_dim=[plot_dim,2],loc_offset=[plot_offset,0],&
1186 &x=x_tot_trim,y=y_tot_trim,z=z_tot_trim,cont_plot=cont_plot,&
1187 &descr=description)
1188
1189 call lvl_ud(-1)
1190
1191 ! clean up
1192 nullify(e_kin_trim,e_pot_trim)
1193 nullify(x_tot_trim,y_tot_trim,z_tot_trim)
1194 call grid_sol_trim%dealloc()
1195 end if
1196 end function decompose_energy
1197
1198 !> Calculate the energy terms in the energy decomposition.
1199 !!
1200 !! \note see explanation of routine in decompose_energy().
1201 !!
1202 !! \return ierr
1203 integer function calc_e(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,eq_2,X,&
1204 &sol,vac,B_aligned,X_id,E_pot,E_kin,E_pot_int,E_kin_int) result(ierr)
1205
1206 use num_vars, only: use_pol_flux_f, n_procs, k_style, &
1207 &norm_disc_prec_eq, norm_disc_prec_x, norm_disc_prec_sol, rank, &
1208 &eq_job_nr, eq_jobs_lims, x_grid_style
1209 use eq_vars, only: vac_perm
1210 use num_utilities, only: c, spline
1211 use grid_utilities, only: calc_int_vol, trim_grid, untrim_grid
1212 use grid_vars, only: alpha, n_alpha
1213 use mpi_utilities, only: get_ser_var
1214 use sol_utilities, only: calc_xuq
1215
1216 character(*), parameter :: rout_name = 'calc_E'
1217
1218 ! input / output
1219 type(modes_type), intent(in) :: mds_x !< general modes variables in perturbation grid
1220 type(modes_type), intent(in) :: mds_sol !< general modes variables in solution grid
1221 type(grid_type), intent(in) :: grid_eq !< equilibrium grid
1222 type(grid_type), intent(in) :: grid_x !< perturbation grid
1223 type(grid_type), intent(in) :: grid_sol !< and solution grid
1224 type(eq_1_type), intent(in) :: eq_1 !< flux equilibrium
1225 type(eq_2_type), intent(in) :: eq_2 !< metric equilibrium
1226 type(x_1_type), intent(in) :: x !< perturbation variables
1227 type(sol_type), intent(in) :: sol !< solution variables
1228 type(vac_type), intent(in) :: vac !< vacuum variables
1229 logical, intent(in) :: b_aligned !< whether grid is field-aligned
1230 integer, intent(in) :: x_id !< nr. of Eigenvalue
1231 complex(dp), intent(inout), allocatable :: e_pot(:,:,:,:) !< potential energy
1232 complex(dp), intent(inout), allocatable :: e_kin(:,:,:,:) !< kinetic energy
1233 complex(dp), intent(inout) :: e_pot_int(7) !< integrated potential energy
1234 complex(dp), intent(inout) :: e_kin_int(2) !< integrated kinetic energy
1235
1236 ! local variables
1237 integer :: norm_id(2) ! untrimmed normal indices for trimmed grids
1238 integer :: id, jd, kd ! counters
1239 integer :: loc_dim(3) ! local dimension
1240 type(grid_type) :: grid_sol_trim ! trimmed sol grid
1241 real(dp), allocatable :: h22(:,:,:), g33(:,:,:), j(:,:,:) ! interpolated h_FD(2,2), g_FD(3,3) and J_FD
1242 real(dp), allocatable :: kappa_n(:,:,:), kappa_g(:,:,:) ! interpolated kappa_n and kappa_g
1243 real(dp), allocatable :: sigma(:,:,:) ! interpolated sigma
1244 real(dp), allocatable :: d2p(:,:,:), rho(:,:,:) ! interpolated Dpres_FD, rho
1245 real(dp), allocatable :: ang_1(:,:,:), ang_2(:,:,:), norm(:) ! coordinates of grid in which to calculate total energy
1246 complex(dp), allocatable :: xuq(:,:,:,:) ! X, U, Q_n and Q_g
1247 complex(dp), allocatable :: e_int_tot(:) ! integrated potential or kinetic energy of all processes
1248 character(len=max_str_ln) :: err_msg ! error message
1249#if ldebug
1250 integer :: plot_dim(3) ! dimensions of plot
1251 integer :: plot_offset(3) ! local offset of plot
1252 real(dp), allocatable :: s(:,:,:) ! interpolated S
1253 complex(dp), allocatable :: du(:,:,:) ! D_par U
1254 complex(dp), allocatable :: du_alt(:,:,:) ! DU calculated from U
1255#endif
1256
1257 ! set loc_dim
1258 loc_dim = [grid_eq%n(1:2),grid_sol%loc_n_r] ! includes ghost regions of width norm_disc_prec_sol
1259
1260#if ldebug
1261 ! set up plot dimensions and offset
1262 plot_dim = [grid_eq%n(1:2),grid_sol%n(3)]
1263 plot_offset = [0,0,grid_sol%i_min-1]
1264#endif
1265
1266 ! allocate local variables
1267 allocate(h22(loc_dim(1),loc_dim(2),loc_dim(3)))
1268 allocate(g33(loc_dim(1),loc_dim(2),loc_dim(3)))
1269 allocate(j(loc_dim(1),loc_dim(2),loc_dim(3)))
1270 allocate(kappa_n(loc_dim(1),loc_dim(2),loc_dim(3)))
1271 allocate(kappa_g(loc_dim(1),loc_dim(2),loc_dim(3)))
1272 allocate(sigma(loc_dim(1),loc_dim(2),loc_dim(3)))
1273 allocate(d2p(loc_dim(1),loc_dim(2),loc_dim(3)))
1274 allocate(rho(loc_dim(1),loc_dim(2),loc_dim(3)))
1275 allocate(ang_1(loc_dim(1),loc_dim(2),loc_dim(3)))
1276 allocate(ang_2(loc_dim(1),loc_dim(2),loc_dim(3)))
1277 allocate(norm(loc_dim(3)))
1278 allocate(xuq(loc_dim(1),loc_dim(2),loc_dim(3),4))
1279 allocate(e_pot(loc_dim(1),loc_dim(2),loc_dim(3),6))
1280 allocate(e_kin(loc_dim(1),loc_dim(2),loc_dim(3),2))
1281#if ldebug
1282 if (debug_calc_e .or. debug_du) then
1283 allocate(du(loc_dim(1),loc_dim(2),loc_dim(3)))
1284 allocate(s(loc_dim(1),loc_dim(2),loc_dim(3)))
1285
1286 ! calculate D_par U
1287 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,eq_2,x,&
1288 &sol,x_id,2,0._dp,du,deriv=.true.)
1289 chckerr('')
1290 end if
1291#endif
1292
1293#if lwith_intel
1294 ! set J to zero to avoid INTEL compiler bug
1295 j = 0._dp
1296#endif
1297
1298 ! interpolate eq->sol
1299 do jd = 1,grid_eq%n(2)
1300 do id = 1,grid_eq%n(1)
1301 ierr = spline(grid_eq%loc_r_F,&
1302 &eq_2%h_FD(id,jd,:,c([2,2],.true.),0,0,0),grid_sol%loc_r_F,&
1303 &h22(id,jd,:),ord=norm_disc_prec_eq)
1304 chckerr('')
1305 ierr = spline(grid_eq%loc_r_F,&
1306 &eq_2%g_FD(id,jd,:,c([3,3],.true.),0,0,0),grid_sol%loc_r_F,&
1307 &g33(id,jd,:),ord=norm_disc_prec_eq)
1308 chckerr('')
1309 ierr = spline(grid_eq%loc_r_F,&
1310 &eq_2%jac_FD(id,jd,:,0,0,0),grid_sol%loc_r_F,&
1311 &j(id,jd,:),ord=norm_disc_prec_eq)
1312 chckerr('')
1313 ierr = spline(grid_eq%loc_r_F,&
1314 &eq_2%kappa_n(id,jd,:),grid_sol%loc_r_F,&
1315 &kappa_n(id,jd,:),ord=norm_disc_prec_eq)
1316 chckerr('')
1317 ierr = spline(grid_eq%loc_r_F,&
1318 &eq_2%kappa_g(id,jd,:),grid_sol%loc_r_F,&
1319 &kappa_g(id,jd,:),ord=norm_disc_prec_eq)
1320 chckerr('')
1321 ierr = spline(grid_eq%loc_r_F,&
1322 &eq_2%sigma(id,jd,:),grid_sol%loc_r_F,&
1323 &sigma(id,jd,:),ord=norm_disc_prec_eq)
1324 chckerr('')
1325#if ldebug
1326 if (debug_calc_e) then
1327 ierr = spline(grid_eq%loc_r_F,&
1328 &eq_2%S(id,jd,:),grid_sol%loc_r_F,&
1329 &s(id,jd,:),ord=norm_disc_prec_eq)
1330 chckerr('')
1331 end if
1332#endif
1333 end do
1334 end do
1335 ierr = spline(grid_eq%loc_r_F,eq_1%pres_FD(:,1),grid_sol%loc_r_F,&
1336 &d2p(1,1,:),ord=norm_disc_prec_eq)
1337 chckerr('')
1338 ierr = spline(grid_eq%loc_r_F,eq_1%rho,grid_sol%loc_r_F,&
1339 &rho(1,1,:),ord=norm_disc_prec_eq)
1340 chckerr('')
1341 do kd = 1,loc_dim(3)
1342 d2p(:,:,kd) = d2p(1,1,kd)
1343 rho(:,:,kd) = rho(1,1,kd)
1344 end do
1345
1346 ! set angles and norm
1347 select case (x_grid_style)
1348 case (1,3) ! equilibrium or enriched
1349 ! interpolate X->sol
1350 if (b_aligned) then
1351 if (use_pol_flux_f) then
1352 do jd = 1,grid_eq%n(2)
1353 do id = 1,grid_eq%n(1)
1354 ierr = spline(grid_x%loc_r_F,&
1355 &grid_x%theta_F(id,jd,:),grid_sol%loc_r_F,&
1356 &ang_1(id,jd,:),ord=norm_disc_prec_x)
1357 chckerr('')
1358 end do
1359 end do
1360 else
1361 do jd = 1,grid_eq%n(2)
1362 do id = 1,grid_eq%n(1)
1363 ierr = spline(grid_x%loc_r_F,&
1364 &grid_x%zeta_F(id,jd,:),grid_sol%loc_r_F,&
1365 &ang_1(id,jd,:),ord=norm_disc_prec_x)
1366 chckerr('')
1367 end do
1368 end do
1369 end if
1370 do jd = 1,n_alpha
1371 ang_2(:,jd,:) = alpha(jd)
1372 end do
1373 else
1374 do jd = 1,grid_eq%n(2)
1375 do id = 1,grid_eq%n(1)
1376 ierr = spline(grid_x%loc_r_F,&
1377 &grid_x%theta_F(id,jd,:),grid_sol%loc_r_F,&
1378 &ang_1(id,jd,:),ord=norm_disc_prec_x)
1379 chckerr('')
1380 ierr = spline(grid_x%loc_r_F,&
1381 &grid_x%zeta_F(id,jd,:),grid_sol%loc_r_F,&
1382 &ang_2(id,jd,:),ord=norm_disc_prec_x)
1383 chckerr('')
1384 end do
1385 end do
1386 end if
1387 case (2) ! solution
1388 ! copy
1389 if (b_aligned) then
1390 if (use_pol_flux_f) then
1391 ang_1 = grid_x%theta_F
1392 else
1393 ang_1 = grid_x%zeta_F
1394 end if
1395 do jd = 1,n_alpha
1396 ang_2(:,jd,:) = alpha(jd)
1397 end do
1398 else
1399 ang_1 = grid_x%theta_F
1400 ang_2 = grid_x%zeta_F
1401 end if
1402 end select
1403 norm = grid_sol%loc_r_F
1404
1405 ! calculate X, U, Q_n and Q_g
1406 do kd = 1,4
1407 ierr = calc_xuq(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,eq_2,x,&
1408 &sol,x_id,kd,0._dp,xuq(:,:,:,kd))
1409 chckerr('')
1410 end do
1411
1412 ! calc kinetic energy
1413 e_kin(:,:,:,1) = rho/h22*xuq(:,:,:,1)*conjg(xuq(:,:,:,1))
1414 select case (k_style)
1415 case (1) ! normalization of full perpendicular component
1416 e_kin(:,:,:,2) = &
1417 &rho*h22*j**2/g33*xuq(:,:,:,2)*conjg(xuq(:,:,:,2))
1418 case (2) ! normalization of only normal component
1419 e_kin(:,:,:,2) = 0._dp
1420 case default
1421 err_msg = 'No normalization style associated with '//&
1422 &trim(i2str(k_style))
1423 ierr = 1
1424 chckerr(err_msg)
1425 end select
1426
1427 ! calc potential energy
1428 e_pot(:,:,:,1) = 1._dp/vac_perm*1._dp/h22*xuq(:,:,:,3)*&
1429 &conjg(xuq(:,:,:,3))
1430 e_pot(:,:,:,2) = 1._dp/vac_perm*h22*j**2/g33*xuq(:,:,:,4)*&
1431 &conjg(xuq(:,:,:,4))
1432 e_pot(:,:,:,3) = -2*d2p*kappa_n*xuq(:,:,:,1)*conjg(xuq(:,:,:,1))
1433 e_pot(:,:,:,4) = -2*d2p*kappa_g*xuq(:,:,:,1)*conjg(xuq(:,:,:,2))
1434 e_pot(:,:,:,5) = -sigma*xuq(:,:,:,4)*conjg(xuq(:,:,:,1))
1435 e_pot(:,:,:,6) = sigma*xuq(:,:,:,3)*conjg(xuq(:,:,:,2))
1436
1437#if ldebug
1438 if (debug_calc_e) then
1439 call writo('Testing whether the total potential energy is &
1440 &equal to the alternative form given in [ADD REF]')
1441 call lvl_ud(1)
1442
1443 ! alternative formulation for E_pot, always real
1444 e_pot(:,:,:,3) = (-2*d2p*kappa_n+sigma*s)*&
1445 &xuq(:,:,:,1)*conjg(xuq(:,:,:,1))
1446 e_pot(:,:,:,4) = -sigma/j*2*rp(xuq(:,:,:,1)*conjg(du))
1447 e_pot(:,:,:,5:6) = 0._dp
1448
1449 call lvl_ud(-1)
1450 end if
1451
1452 if (debug_x_norm) then
1453 call writo('Plotting the squared amplitude of the normal &
1454 &component')
1455 call lvl_ud(1)
1456
1457 call plot_hdf5('X_norm', 'TEST_X_norm_POST_'//trim(i2str(x_id)),&
1458 &rp(xuq(:,:,:,1)*conjg(xuq(:,:,:,1))),&
1459 &tot_dim=plot_dim,loc_offset=plot_offset)
1460
1461 call lvl_ud(-1)
1462 end if
1463
1464 if (debug_du .and. b_aligned) then
1465 call writo('Testing whether DU is indeed the parallel &
1466 &derivative of U')
1467 call lvl_ud(1)
1468
1469 allocate(du_alt(loc_dim(1),loc_dim(2),loc_dim(3)))
1470
1471 ! derive
1472 do kd = 1,loc_dim(3)
1473 do jd = 1,loc_dim(2)
1474 ierr = spline(ang_1(:,jd,kd),xuq(:,jd,kd,2),ang_1(:,jd,kd),&
1475 &du_alt(:,jd,kd),ord=norm_disc_prec_sol,deriv=1)
1476 chckerr('')
1477 end do
1478 end do
1479
1480 ! plot real part
1481 call plot_hdf5('RE U','TEST_RE_U_'//&
1482 &trim(i2str(x_id)),rp(xuq(:,:,:,2)),tot_dim=plot_dim,&
1483 &loc_offset=plot_offset)
1484 call plot_diff_hdf5(rp(du),rp(du_alt),&
1485 &'TEST_RE_DU_'//trim(i2str(x_id)),plot_dim,&
1486 &plot_offset,descr='To test whether DU is &
1487 &parallel derivative of U',output_message=.true.)
1488
1489 ! plot imaginary part
1490 call plot_hdf5('IM U','TEST_IM_U_'//&
1491 &trim(i2str(x_id)),ip(xuq(:,:,:,2)),tot_dim=plot_dim,&
1492 &loc_offset=plot_offset)
1493 call plot_diff_hdf5(ip(du),ip(du_alt),&
1494 &'TEST_IM_DU_'//trim(i2str(x_id)),plot_dim,&
1495 &plot_offset,descr='To test whether DU is &
1496 &parallel derivative of U',output_message=.true.)
1497
1498 deallocate(du_alt)
1499
1500 call lvl_ud(-1)
1501 end if
1502#endif
1503
1504 ! trim sol grid
1505 ierr = trim_grid(grid_sol,grid_sol_trim,norm_id)
1506 chckerr('')
1507
1508 ! add ghost region of width one to the right of the interval
1509 if (rank.lt.n_procs-1) norm_id(2) = norm_id(2)+1 ! adjust norm_id as well
1510
1511 ! integrate energy using ghosted variables
1512 ierr = calc_int_vol(ang_1(:,:,norm_id(1):norm_id(2)),&
1513 &ang_2(:,:,norm_id(1):norm_id(2)),norm(norm_id(1):norm_id(2)),&
1514 &j(:,:,norm_id(1):norm_id(2)),&
1515 &e_kin(:,:,norm_id(1):norm_id(2),:),e_kin_int)
1516 chckerr('')
1517 ierr = calc_int_vol(ang_1(:,:,norm_id(1):norm_id(2)),&
1518 &ang_2(:,:,norm_id(1):norm_id(2)),norm(norm_id(1):norm_id(2)),&
1519 &j(:,:,norm_id(1):norm_id(2)),&
1520 &e_pot(:,:,norm_id(1):norm_id(2),:),e_pot_int(1:6))
1521 chckerr('')
1522
1523 ! add vacuum energy if last process and first equilibrium job
1524 e_pot_int(7) = 0._dp
1525 write(*,*) '¡¡¡¡¡ NO VACUUM !!!!!'
1526 if (1.eq.2 .and. rank.eq.n_procs-1 .and. eq_job_nr.eq.size(eq_jobs_lims,2)) then
1527 do jd = 1,sol%n_mod
1528 do kd = 1,sol%n_mod
1529 e_pot_int(7) = e_pot_int(7) + &
1530 &conjg(sol%vec(kd,grid_sol%loc_n_r,x_id)) * &
1531 &vac%res(kd,jd) * sol%vec(jd,grid_sol%loc_n_r,x_id)
1532 end do
1533 end do
1534 end if
1535
1536 ! bundle all processes
1537 do kd = 1,2
1538 ierr = get_ser_var([e_kin_int(kd)],e_int_tot,scatter=.true.)
1539 chckerr('')
1540 e_kin_int(kd) = sum(e_int_tot)
1541 deallocate(e_int_tot)
1542 end do
1543 do kd = 1,7
1544 ierr = get_ser_var([e_pot_int(kd)],e_int_tot,scatter=.true.)
1545 chckerr('')
1546 e_pot_int(kd) = sum(e_int_tot)
1547 deallocate(e_int_tot)
1548 end do
1549
1550 ! clean up
1551 call grid_sol_trim%dealloc()
1552 end function calc_e
1553
1554 !> Print solution quantities to an output file:
1555 !!
1556 !! - sol:
1557 !! - \c val
1558 !! - \c vec
1559 !!
1560 !! If \c rich_lvl is provided, \c "_[R_rich_lvl]" is appended to the data
1561 !! name if it is \c >0.
1562 !!
1563 !! \return ierr
1564 integer function print_output_sol(grid,sol,data_name,rich_lvl,&
1565 &remove_previous_arrs) result(ierr)
1566
1567 use num_vars, only: pb3d_name, sol_n_procs, rank
1568 use hdf5_ops, only: print_hdf5_arrs
1569 use hdf5_vars, only: dealloc_var_1d, var_1d_type, &
1570 &max_dim_var_1d
1571 use grid_utilities, only: trim_grid
1572
1573 character(*), parameter :: rout_name = 'print_output_sol'
1574
1575 ! input / output
1576 type(grid_type), intent(in) :: grid !< solution grid variables
1577 type(sol_type), intent(in) :: sol !< solution variables
1578 character(len=*), intent(in) :: data_name !< name under which to store
1579 integer, intent(in), optional :: rich_lvl !< Richardson level to print
1580 logical, intent(in), optional :: remove_previous_arrs !< remove previous variables if present
1581
1582 ! local variables
1583 type(var_1d_type), allocatable, target :: sol_1d(:) ! 1D equivalent of eq. variables
1584 type(var_1d_type), pointer :: sol_1d_loc => null() ! local element in sol_1D
1585 type(grid_type) :: grid_trim ! trimmed solution grid
1586 integer :: id ! counter
1587
1588 ! initialize ierr
1589 ierr = 0
1590
1591 ! only write to output file if at least one solution
1592 if (size(sol%val).gt.0) then
1593 ! user output
1594 call writo('Write solution variables to output file')
1595 call lvl_ud(1)
1596
1597 ! trim grids
1598 ierr = trim_grid(grid,grid_trim)
1599 chckerr('')
1600
1601 ! Set up the 1D equivalents of the solution variables
1602 allocate(sol_1d(max_dim_var_1d))
1603
1604 id = 1
1605
1606 ! RE_sol_val
1607 sol_1d_loc => sol_1d(id); id = id+1
1608 sol_1d_loc%var_name = 'RE_sol_val'
1609 allocate(sol_1d_loc%tot_i_min(1),sol_1d_loc%tot_i_max(1))
1610 allocate(sol_1d_loc%loc_i_min(1),sol_1d_loc%loc_i_max(1))
1611 sol_1d_loc%tot_i_min = [1]
1612 sol_1d_loc%tot_i_max = [size(sol%val)]
1613 sol_1d_loc%loc_i_min = [1]
1614 sol_1d_loc%loc_i_max = [size(sol%val)]
1615 allocate(sol_1d_loc%p(size(sol%val)))
1616 sol_1d_loc%p = rp(sol%val)
1617
1618 ! IM_sol_val
1619 sol_1d_loc => sol_1d(id); id = id+1
1620 sol_1d_loc%var_name = 'IM_sol_val'
1621 allocate(sol_1d_loc%tot_i_min(1),sol_1d_loc%tot_i_max(1))
1622 allocate(sol_1d_loc%loc_i_min(1),sol_1d_loc%loc_i_max(1))
1623 sol_1d_loc%tot_i_min = [1]
1624 sol_1d_loc%tot_i_max = [size(sol%val)]
1625 sol_1d_loc%loc_i_min = [1]
1626 sol_1d_loc%loc_i_max = [size(sol%val)]
1627 allocate(sol_1d_loc%p(size(sol%val)))
1628 sol_1d_loc%p = ip(sol%val)
1629
1630 ! RE_sol_vec
1631 sol_1d_loc => sol_1d(id); id = id+1
1632 sol_1d_loc%var_name = 'RE_sol_vec'
1633 allocate(sol_1d_loc%tot_i_min(3),sol_1d_loc%tot_i_max(3))
1634 allocate(sol_1d_loc%loc_i_min(3),sol_1d_loc%loc_i_max(3))
1635 sol_1d_loc%loc_i_min = [1,grid_trim%i_min,1]
1636 sol_1d_loc%loc_i_max = [sol%n_mod,grid_trim%i_max,size(sol%vec,3)]
1637 sol_1d_loc%tot_i_min = [1,1,1]
1638 sol_1d_loc%tot_i_max = [sol%n_mod,grid_trim%n(3),size(sol%vec,3)]
1639 allocate(sol_1d_loc%p(size(sol%vec)))
1640 sol_1d_loc%p = reshape(rp(sol%vec),[size(sol%vec)])
1641
1642 ! IM_sol_vec
1643 sol_1d_loc => sol_1d(id); id = id+1
1644 sol_1d_loc%var_name = 'IM_sol_vec'
1645 allocate(sol_1d_loc%tot_i_min(3),sol_1d_loc%tot_i_max(3))
1646 allocate(sol_1d_loc%loc_i_min(3),sol_1d_loc%loc_i_max(3))
1647 sol_1d_loc%loc_i_min = [1,grid_trim%i_min,1]
1648 sol_1d_loc%loc_i_max = [sol%n_mod,grid_trim%i_max,size(sol%vec,3)]
1649 sol_1d_loc%tot_i_min = [1,1,1]
1650 sol_1d_loc%tot_i_max = [sol%n_mod,grid_trim%n(3),size(sol%vec,3)]
1651 allocate(sol_1d_loc%p(size(sol%vec)))
1652 sol_1d_loc%p = reshape(ip(sol%vec),[size(sol%vec)])
1653
1654 ! write
1655 ! Note: if processes that are not used in the solve do not have the
1656 ! solution vector and should therefore not be used to write the
1657 ! output.
1658 if (sol_n_procs.gt.1 .or. rank.eq.0) then
1659 ierr = print_hdf5_arrs(sol_1d(1:id-1),pb3d_name,&
1660 &trim(data_name),rich_lvl=rich_lvl,&
1661 &remove_previous_arrs=remove_previous_arrs)
1662 chckerr('')
1663 end if
1664
1665 ! clean up
1666 call grid_trim%dealloc()
1667 call dealloc_var_1d(sol_1d)
1668 nullify(sol_1d_loc)
1669
1670 ! user output
1671 call lvl_ud(-1)
1672 else
1673 ! user output
1674 call writo('No solutions to write')
1675 end if
1676 end function print_output_sol
1677end module sol_ops
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::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::plot_hdf5
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Definition output_ops.f90:137
output_ops::print_ex_2d
Print 2-D output on a file.
Definition output_ops.f90:47
sol_utilities::calc_xuq
Calculates the normal ( ) or geodesic ( ) components of the plasma perturbation or the magnetic pert...
Definition sol_utilities.f90:56
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition eq_vars.f90:27
eq_vars::r_0
real(dp), public r_0
independent normalization constant for nondimensionalization
Definition eq_vars.f90:42
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
eq_vars::b_0
real(dp), public b_0
derived normalization constant for nondimensionalization
Definition eq_vars.f90:45
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition grid_utilities.f90:4
grid_utilities::untrim_grid
integer function, public untrim_grid(grid_in, grid_out, size_ghost)
Untrims a trimmed grid by introducing an asymmetric ghost regions at the right.
Definition grid_utilities.f90:1764
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_vec_comp
integer function, public calc_vec_comp(grid, grid_eq, eq_1, eq_2, v_com, norm_disc_prec, v_mag, base_name, max_transf, v_flux_tor, v_flux_pol, xyz, compare_tor_pos)
Calculates contra- and covariant components of a vector in multiple coordinate systems.
Definition grid_utilities.f90:1859
grid_utilities::calc_int_vol
integer function, public calc_int_vol(ang_1, ang_2, norm, j, f, f_int)
Calculates volume integral on a 3D grid.
Definition grid_utilities.f90:1320
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::n_alpha
integer, public n_alpha
nr. of field-lines
Definition grid_vars.f90:23
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
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
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_vars
Numerical variables used by most other modules.
Definition num_vars.f90:4
num_vars::ex_max_size
integer, parameter, public ex_max_size
maximum size of matrices for external plot
Definition num_vars.f90:137
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::norm_disc_prec_sol
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition num_vars.f90:122
num_vars::pi
real(dp), parameter, public pi
Definition num_vars.f90:83
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::use_normalization
logical, public use_normalization
whether to use normalization or not
Definition num_vars.f90:115
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::max_deriv
integer, parameter, public max_deriv
highest derivatives for metric factors in Flux coords.
Definition num_vars.f90:52
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::pert_mult_factor_post
real(dp), public pert_mult_factor_post
factor with which to multiply perturbation strength for POST
Definition num_vars.f90:176
num_vars::ex_plot_style
integer, public ex_plot_style
external plot style (1: GNUPlot, 2: Bokeh for 2D, Mayavi for 3D)
Definition num_vars.f90:175
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::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
sol_ops
Operations on the solution vectors such as decompososing the energy, etc...
Definition sol_ops.f90:4
sol_ops::print_output_sol
integer function, public print_output_sol(grid, sol, data_name, rich_lvl, remove_previous_arrs)
Print solution quantities to an output file:
Definition sol_ops.f90:1566
sol_ops::plot_sol_vals
subroutine, public plot_sol_vals(sol, last_unstable_id)
Plots Eigenvalues.
Definition sol_ops.f90:37
sol_ops::decompose_energy
integer function, public decompose_energy(mds_x, mds_sol, grid_eq, grid_x, grid_sol, eq_1, eq_2, x, sol, vac, x_id, b_aligned, xyz, e_pot_int, e_kin_int)
Decomposes the plasma potential and kinetic energy in its individual terms for an individual Eigenval...
Definition sol_ops.f90:991
sol_ops::debug_calc_e
logical, public debug_calc_e
plot debug information for calc_E()
Definition sol_ops.f90:29
sol_ops::debug_x_norm
logical, public debug_x_norm
plot debug information X_norm
Definition sol_ops.f90:30
sol_ops::plot_sol_vec
integer function, public plot_sol_vec(mds_x, mds_sol, grid_eq, grid_x, grid_sol, eq_1, eq_2, x, sol, xyz, x_id, plot_var)
Plots Eigenvectors.
Definition sol_ops.f90:116
sol_ops::plot_harmonics
integer function, public plot_harmonics(mds, grid_sol, sol, x_id, res_surf)
Plots the harmonics and their maximum in 2-D.
Definition sol_ops.f90:646
sol_ops::debug_du
logical, public debug_du
plot debug information for calculation of DU
Definition sol_ops.f90:31
sol_ops::debug_plot_sol_vec
logical, public debug_plot_sol_vec
plot debug information for plot_sol_vec()
Definition sol_ops.f90:27
sol_utilities
Numerical utilities related to the solution vectors.
Definition sol_utilities.f90:4
sol_utilities::calc_tot_sol_vec
integer function, public calc_tot_sol_vec(mds, grid_sol, sol_vec_loc, sol_vec_tot, deriv)
Calculate the total version of the solution vector from the local version.
Definition sol_utilities.f90:581
sol_vars
Variables pertaining to the solution quantities.
Definition sol_vars.f90:4
str_utilities
Operations on strings.
Definition str_utilities.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition str_utilities.f90:18
str_utilities::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
vac_vars
Variables pertaining to the vacuum quantities.
Definition vac_vars.f90:4
x_vars
Variables pertaining to the perturbation quantities.
Definition X_vars.f90:4
x_vars::mds_x
type(modes_type), public mds_x
modes variables for perturbation grid
Definition X_vars.f90:124
x_vars::mds_sol
type(modes_type), public mds_sol
modes variables for solution grid
Definition X_vars.f90:125
x_vars::min_nm_x
integer, public min_nm_x
minimum for the high-n theory (debable)
Definition X_vars.f90:130
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
sol_vars::sol_type
solution type
Definition sol_vars.f90:30
vac_vars::vac_type
vacuum type
Definition vac_vars.f90:46
x_vars::modes_type
mode number type
Definition X_vars.f90:36
x_vars::x_1_type
vectorial perturbation type
Definition X_vars.f90:51
this documentation was created with Doxygen 1.16.1 by Toon Weyens