PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
sol_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module 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.
28  logical :: debug_plot_harmonics = .false.
29  logical :: debug_calc_e = .false.
30  logical :: debug_x_norm = .false.
31  logical :: debug_du = .false.
32 #endif
33 
34 contains
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
41  integer, intent(in) :: last_unstable_id
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 
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
134  type(modes_type), intent(in) :: mds_sol
135  type(grid_type), intent(in) :: grid_eq
136  type(grid_type), intent(in) :: grid_x
137  type(grid_type), intent(in) :: grid_sol
138  type(eq_1_type), intent(in) :: eq_1
139  type(eq_2_type), intent(in) :: eq_2
140  type(x_1_type), intent(in) :: x
141  type(sol_type), intent(in) :: sol
142  real(dp), intent(in) :: xyz(:,:,:,:)
143  integer, intent(in) :: x_id
144  logical, intent(in) :: plot_var(2)
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 
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
659  type(grid_type), intent(in) :: grid_sol
660  type(sol_type), intent(in) :: sol
661  integer, intent(in) :: x_id
662  real(dp), intent(in) :: res_surf(:,:)
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 
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
999  type(modes_type), intent(in) :: mds_sol
1000  type(grid_type), intent(in) :: grid_eq
1001  type(grid_type), intent(in) :: grid_x
1002  type(grid_type), intent(in) :: grid_sol
1003  type(eq_1_type), intent(in) :: eq_1
1004  type(eq_2_type), intent(in) :: eq_2
1005  type(x_1_type), intent(in) :: x
1006  type(sol_type), intent(in) :: sol
1007  type(vac_type), intent(in) :: vac
1008  integer, intent(in) :: x_id
1009  logical, intent(in) :: b_aligned
1010  real(dp), intent(in), optional :: xyz(:,:,:,:)
1011  complex(dp), intent(inout), optional :: e_pot_int(7)
1012  complex(dp), intent(inout), optional :: e_kin_int(2)
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 
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
1220  type(modes_type), intent(in) :: mds_sol
1221  type(grid_type), intent(in) :: grid_eq
1222  type(grid_type), intent(in) :: grid_x
1223  type(grid_type), intent(in) :: grid_sol
1224  type(eq_1_type), intent(in) :: eq_1
1225  type(eq_2_type), intent(in) :: eq_2
1226  type(x_1_type), intent(in) :: x
1227  type(sol_type), intent(in) :: sol
1228  type(vac_type), intent(in) :: vac
1229  logical, intent(in) :: b_aligned
1230  integer, intent(in) :: x_id
1231  complex(dp), intent(inout), allocatable :: e_pot(:,:,:,:)
1232  complex(dp), intent(inout), allocatable :: e_kin(:,:,:,:)
1233  complex(dp), intent(inout) :: e_pot_int(7)
1234  complex(dp), intent(inout) :: e_kin_int(2)
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 
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
1577  type(sol_type), intent(in) :: sol
1578  character(len=*), intent(in) :: data_name
1579  integer, intent(in), optional :: rich_lvl
1580  logical, intent(in), optional :: remove_previous_arrs
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
1677 end module sol_ops
num_utilities::c
integer function, public c(ij, sym, n, lim_n)
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
Definition: num_utilities.f90:2556
mpi_utilities::get_ser_var
Gather parallel variable in serial version on group master.
Definition: MPI_utilities.f90:55
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
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
hdf5_vars::dealloc_var_1d
Deallocates 1D variable.
Definition: HDF5_vars.f90:68
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
eq_vars
Variables that have to do with equilibrium quantities and the grid used in the calculations:
Definition: eq_vars.f90:27
num_vars::use_normalization
logical, public use_normalization
whether to use normalization or not
Definition: num_vars.f90:115
eq_vars::vac_perm
real(dp), public vac_perm
either usual mu_0 (default) or normalized
Definition: eq_vars.f90:48
mpi_utilities
Numerical utilities related to MPI.
Definition: MPI_utilities.f90:20
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
num_vars::ex_max_size
integer, parameter, public ex_max_size
maximum size of matrices for external plot
Definition: num_vars.f90:137
sol_ops::debug_du
logical, public debug_du
plot debug information for calculation of DU
Definition: sol_ops.f90:31
num_vars::norm_disc_prec_eq
integer, public norm_disc_prec_eq
precision for normal discretization for equilibrium
Definition: num_vars.f90:120
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
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
hdf5_ops
Operations on HDF5 and XDMF variables.
Definition: HDF5_ops.f90:27
num_vars::eq_job_nr
integer, public eq_job_nr
nr. of eq job
Definition: num_vars.f90:79
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
num_vars::n_procs
integer, public n_procs
nr. of MPI processes
Definition: num_vars.f90:69
num_vars::norm_disc_prec_sol
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition: num_vars.f90:122
hdf5_vars
Variables pertaining to HDF5 and XDMF.
Definition: HDF5_vars.f90:4
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
sol_ops::calc_e
integer function calc_e(mds_X, mds_sol, grid_eq, grid_X, grid_sol, eq_1, eq_2, X, sol, vac, B_aligned, X_id, E_pot, E_kin, E_pot_int, E_kin_int)
Calculate the energy terms in the energy decomposition.
Definition: sol_ops.f90:1205
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
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
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition: str_utilities.f90:54
sol_ops::debug_calc_e
logical, public debug_calc_e
plot debug information for calc_E()
Definition: sol_ops.f90:29
eq_vars::r_0
real(dp), public r_0
independent normalization constant for nondimensionalization
Definition: eq_vars.f90:42
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
vac_vars
Variables pertaining to the vacuum quantities.
Definition: vac_vars.f90:4
grid_vars::alpha
real(dp), dimension(:), allocatable, public alpha
field line label alpha
Definition: grid_vars.f90:28
sol_vars::sol_type
solution type
Definition: sol_vars.f90:30
sol_utilities
Numerical utilities related to the solution vectors.
Definition: sol_utilities.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
hdf5_vars::var_1d_type
1D equivalent of multidimensional variables, used for internal HDF5 storage.
Definition: HDF5_vars.f90:48
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::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
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
num_vars::eq_style
integer, public eq_style
either 1 (VMEC) or 2 (HELENA)
Definition: num_vars.f90:89
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
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
sol_ops::debug_x_norm
logical, public debug_x_norm
plot debug information X_norm
Definition: sol_ops.f90:30
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
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
x_vars::modes_type
mode number type
Definition: X_vars.f90:36
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
sol_utilities::calc_xuq
Calculates the normal ( ) or geodesic ( ) components of the plasma perturbation or the magnetic pert...
Definition: sol_utilities.f90:56
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
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
num_vars::k_style
integer, public k_style
style for kinetic energy
Definition: num_vars.f90:93
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
num_vars::pi
real(dp), parameter, public pi
Definition: num_vars.f90:83
grid_utilities
Numerical utilities related to the grids and different coordinate systems.
Definition: grid_utilities.f90:4
num_vars::pb3d_name
character(len=max_str_ln), public pb3d_name
name of PB3D output file
Definition: num_vars.f90:139
vac_vars::vac_type
vacuum type
Definition: vac_vars.f90:46
sol_ops::plot_sol_vals
subroutine, public plot_sol_vals(sol, last_unstable_id)
Plots Eigenvalues.
Definition: sol_ops.f90:37
sol_ops
Operations on the solution vectors such as decompososing the energy, etc...
Definition: sol_ops.f90:4
grid_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
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
messages::lvl_ud
subroutine, public lvl_ud(inc)
Increases/decreases lvl of output.
Definition: messages.f90:254
num_vars::no_plots
logical, public no_plots
no plots made
Definition: num_vars.f90:140
sol_ops::debug_plot_sol_vec
logical, public debug_plot_sol_vec
plot debug information for plot_sol_vec()
Definition: sol_ops.f90:27
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::use_pol_flux_f
logical, public use_pol_flux_f
whether poloidal flux is used in F coords.
Definition: num_vars.f90:114
x_vars::x_1_type
vectorial perturbation type
Definition: X_vars.f90:51
sol_vars
Variables pertaining to the solution quantities.
Definition: sol_vars.f90:4
num_vars::eq_jobs_lims
integer, dimension(:,:), allocatable, public eq_jobs_lims
data about eq jobs: [ , ] for all jobs
Definition: num_vars.f90:77
grid_vars::n_alpha
integer, public n_alpha
nr. of field-lines
Definition: grid_vars.f90:23
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
num_vars::sol_n_procs
integer, public sol_n_procs
nr. of MPI processes for solution with SLEPC
Definition: num_vars.f90:70
num_vars::rank
integer, public rank
MPI rank.
Definition: num_vars.f90:68
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
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
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::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
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
eq_vars::b_0
real(dp), public b_0
derived normalization constant for nondimensionalization
Definition: eq_vars.f90:45
num_vars::norm_disc_prec_x
integer, public norm_disc_prec_x
precision for normal discretization for perturbation
Definition: num_vars.f90:121
this documentation was created with Doxygen 1.8.17 by Toon Weyens