PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
sol_utilities.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module sol_utilities
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
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 
16  implicit none
17  private
18  public calc_xuq, calc_tot_sol_vec, calc_loc_sol_vec, interp_v_spline
19 #if ldebug
20  public debug_calc_xuq_arr, debug_interp_v_spline
21 #endif
22 
23  ! global variables
24 #if ldebug
25  logical :: debug_calc_xuq_arr = .false.
26  logical :: debug_interp_v_spline = .false.
27 #endif
28 
29  ! interfaces
30 
56  interface calc_xuq
58  module procedure calc_xuq_ind
60  module procedure calc_xuq_arr
61  end interface
62 
63 contains
65  integer function calc_xuq_arr(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,&
66  &eq_2,X,sol,X_id,XUQ_style,time,XUQ,deriv) result(ierr)
67 
68  use num_vars, only: use_pol_flux_f, norm_disc_prec_sol, x_grid_style
69  use num_utilities, only: con2dis, spline
70  use x_utilities, only: trim_modes
71 #if ldebug
72  use eq_vars, only: max_flux_f
73  use num_utilities, only: calc_int
74 #endif
75 
76  character(*), parameter :: rout_name = 'calc_XUQ_arr'
77 
78  ! input / output
79  type(modes_type), intent(in), target :: mds_x
80  type(modes_type), intent(in), target :: mds_sol
81  type(grid_type), intent(in) :: grid_eq
82  type(grid_type), intent(in) :: grid_x
83  type(grid_type), intent(in) :: grid_sol
84  type(eq_1_type), intent(in) :: eq_1
85  type(eq_2_type), intent(in) :: eq_2
86  type(x_1_type), intent(in), target :: x
87  type(sol_type), intent(in) :: sol
88  integer, intent(in) :: x_id
89  integer, intent(in) :: xuq_style
90  real(dp), intent(in) :: time(:)
91  complex(dp), intent(inout) :: xuq(:,:,:,:)
92  logical, intent(in), optional :: deriv
93 
94  ! local variables
95  integer :: m ! secondary mode number
96  integer :: kd_tot ! kd in total grid
97  integer :: ind_x ! index of mode in X tables
98  integer :: ind_sol ! index of mode in sol tables
99  integer :: n_mod_tot ! total number of modes
100  integer :: n_t ! number of time points
101  integer :: n_loc, m_loc ! local mode numbers
102  integer :: id, jd, kd ! counters
103  integer :: xuq_dims(4) ! dimensions of XUQ
104  integer :: kdl_x(2), kdl_sol(2) ! limits on normal index for a mode in X and sol grids
105  integer :: id_lim_x(2), id_lim_sol(2) ! limits on total modes
106  integer, pointer :: sec_x_loc(:,:), sec_sol_loc(:,:) ! pointers to secondary mode variables
107  real(dp), allocatable :: expon(:,:,:) ! exponent in solution grid
108  real(dp), allocatable :: par_fac(:) ! multiplicative factor due to parallel derivative, without iu, in perturbation grid
109  real(dp), allocatable :: jq(:) ! iota or q, interpolated at perturbation grid
110  real(dp), allocatable :: s(:,:,:), inv_j(:,:,:) ! S and 1/J, interpolated at perturbation grid
111  real(dp), pointer :: r_x_loc(:), r_sol_loc(:) ! local r_X and r_sol for a mode
112  real(dp), pointer :: theta_f(:,:,:), zeta_f(:,:,:) ! angles in solution normal grid
113  complex(dp) :: sqrt_sol_val_norm ! normalized sqrt(sol_val)
114  complex(dp), allocatable :: fac_0(:,:,:), fac_1(:,:,:) ! factor to be multiplied with X and DX in perturbation grid
115  complex(dp), allocatable :: xuq_loc(:,:) ! complex XUQ without time at a normal point
116  complex(dp), allocatable :: sol_vec(:,:) ! total solution vector in solution grid
117  complex(dp), allocatable :: sol_vec_tot(:,:,:) ! total solution vector in solution grid and derivatives
118  complex(dp), allocatable :: dsol_vec(:,:) ! normal derivative of sol_vec in solution grid
119  complex(dp), pointer :: u0(:,:,:), u1(:,:,:) ! U0 and U1 or parallel derivative
120  logical :: extrap = .true. ! whether extrapolation is used
121  logical :: deriv_loc ! local copy of deriv
122  character(len=max_str_ln) :: err_msg ! error message
123 
124 #if ldebug
125  real(dp), allocatable :: sol_vec_alt(:) ! sol_vec calculated from Dsol_vec in solution grid
126 #endif
127 
128  ! initialize ierr
129  ierr = 0
130 
131  ! set up dimensions
132  n_t = size(time)
133  xuq_dims = shape(xuq)
134 
135  ! set up local copy of deriv
136  deriv_loc = .false.
137  if (present(deriv)) deriv_loc = deriv
138 
139  ! tests
140  if (grid_eq%n(1).ne.grid_x%n(1) .or. grid_eq%n(2).ne.grid_x%n(2)) then
141  ierr = 1
142  err_msg = 'Grids need to be compatible'
143  chckerr(err_msg)
144  end if
145  if (grid_eq%n(1).ne.xuq_dims(1) .or. grid_eq%n(2).ne.xuq_dims(2) .or. &
146  &grid_sol%loc_n_r.ne.xuq_dims(3) .or. n_t.ne.xuq_dims(4)) then
147  ierr = 1
148  err_msg = 'XUQ needs to have the correct dimensions'
149  chckerr(err_msg)
150  end if
151 
152  ! limit input modes to output modes
153  ! (X grid should comprise sol grid)
154  ierr = trim_modes(mds_x,mds_sol,id_lim_x,id_lim_sol)
155  chckerr('')
156  sec_x_loc => mds_x%sec(id_lim_x(1):id_lim_x(2),:)
157  sec_sol_loc => mds_sol%sec(id_lim_sol(1):id_lim_sol(2),:)
158 
159  n_mod_tot = size(sec_x_loc,1)
160 
161  if (size(sec_x_loc,1).ne.size(sec_sol_loc,1)) then
162  ierr = 1
163  err_msg = 'modes X and sol have inconsistent sizes'
164  chckerr(err_msg)
165  end if
166 
167  ! set up normalized sqrt(sol_val)
168  sqrt_sol_val_norm = sqrt(sol%val(x_id))
169  if (ip(sqrt_sol_val_norm).gt.0) sqrt_sol_val_norm = - sqrt_sol_val_norm ! exploding solution, not the decaying one
170  sqrt_sol_val_norm = sqrt_sol_val_norm / abs(sqrt_sol_val_norm) ! normalize
171 
172  ! allocate helper variables
173  allocate(jq(xuq_dims(3)))
174  allocate(s(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
175  if ((xuq_style.eq.3 .or. xuq_style.eq.4)) &
176  &allocate(inv_j(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
177  allocate(xuq_loc(xuq_dims(1),xuq_dims(2)))
178  allocate(dsol_vec(sol%n_mod,xuq_dims(3)))
179 
180  ! interpolate eq->sol
181  if (use_pol_flux_f) then
182  ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
183  &grid_sol%loc_r_F,jq,ord=norm_disc_prec_sol)
184  chckerr('')
185  else
186  ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
187  &grid_sol%loc_r_F,jq,ord=norm_disc_prec_sol)
188  chckerr('')
189  end if
190  do jd = 1,xuq_dims(2)
191  do id = 1,xuq_dims(1)
192  ierr = spline(grid_eq%loc_r_F,eq_2%S(id,jd,:),&
193  &grid_sol%loc_r_F,s(id,jd,:),ord=norm_disc_prec_sol)
194  chckerr('')
195  if ((xuq_style.eq.3 .or. xuq_style.eq.4)) then
196  ierr = spline(grid_eq%loc_r_F,&
197  &1._dp/eq_2%jac_FD(id,jd,:,0,0,0),grid_sol%loc_r_F,&
198  &inv_j(id,jd,:),ord=norm_disc_prec_sol)
199  chckerr('')
200  end if
201  end do
202  end do
203 
204  ! initialize XUQ
205  xuq = 0._dp
206 
207  ! set up normal derivative of X vec
208  ! Note: In X_style 2 (fast), the mode indices vary as a function as the
209  ! normal coordinate. All of these contiguous ranges are tabulated in
210  ! variable 'sec'.
211  if (xuq_style.eq.2 .or. xuq_style.eq.4) then
212  ! set up variables
213  allocate(sol_vec_tot(n_mod_tot,grid_sol%loc_n_r,0:1))
214  sol_vec_tot = 0._dp
215 
216  ! convert to total solution vector and derivate
217  do id = 0,1
218  ierr = calc_tot_sol_vec(mds_sol,grid_sol,sol%vec(:,:,x_id),&
219  &sol_vec,deriv=id)
220  chckerr('')
221  sol_vec_tot(:,:,id) = sol_vec
222  end do
223 
224 #if ldebug
225  if (debug_calc_xuq_arr) then
226  do m = 1,n_mod_tot
227  call writo('For mode '//trim(i2str(m))//', testing &
228  &whether Dsol_vec is correct by comparing its integral &
229  &with original sol_vec')
230  kdl_sol = sec_sol_loc(m,2:3)
231  allocate(sol_vec_alt(1:kdl_sol(2)-kdl_sol(1)+1))
232  sol_vec_alt = 0._dp
233  ierr = calc_int(rp(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),1)),&
234  &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),sol_vec_alt)
235  chckerr('')
236  call print_ex_2d(['TEST_RE_Dsol_vec'],'',reshape(&
237  &[rp(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),0)),&
238  &sol_vec_alt],[kdl_sol(2)-kdl_sol(1)+1,2]),x=reshape(&
239  &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),&
240  &[kdl_sol(2)-kdl_sol(1)+1,1])/max_flux_f*2*pi)
241  ierr = calc_int(ip(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),1)),&
242  &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),sol_vec_alt)
243  chckerr('')
244  call print_ex_2d(['TEST_IM_Dsol_vec'],'',reshape(&
245  &[ip(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),0)),&
246  &sol_vec_alt],[kdl_sol(2)-kdl_sol(1)+1,2]),x=reshape(&
247  &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),&
248  &[kdl_sol(2)-kdl_sol(1)+1,1])/max_flux_f*2*pi)
249  deallocate(sol_vec_alt)
250  end do
251  end if
252 #endif
253 
254  ! convert back to local solution vector
255  ierr = calc_loc_sol_vec(mds_sol,grid_sol%i_min,sol_vec_tot(:,:,1),&
256  &dsol_vec)
257  chckerr('')
258 
259  ! clean up
260  deallocate(sol_vec,sol_vec_tot)
261  else
262  dsol_vec = 0._dp
263  end if
264 
265  ! interpolate or point angles
266  select case (x_grid_style)
267  case (1,3) ! equilibrium or enriched
268  allocate(theta_f(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
269  allocate(zeta_f(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
270  do jd = 1,xuq_dims(2)
271  do id = 1,xuq_dims(1)
272  ierr = spline(grid_eq%loc_r_F,grid_eq%theta_F(id,jd,:),&
273  &grid_sol%loc_r_F,theta_f(id,jd,:))
274  chckerr('')
275  ierr = spline(grid_eq%loc_r_F,grid_eq%zeta_F(id,jd,:),&
276  &grid_sol%loc_r_F,zeta_f(id,jd,:))
277  chckerr('')
278  end do
279  end do
280  case (2) ! solution
281  theta_f => grid_x%theta_F
282  zeta_f => grid_x%zeta_F
283  end select
284 
285  ! iterate over all secondary modes
286  do m = 1,n_mod_tot
287  ! test
288  if (sec_x_loc(m,1).ne.sec_sol_loc(m,1)) then
289  ierr = 1
290  err_msg = 'For mode '//trim(i2str(m))//&
291  &', no consistency for modes in grid_X and grid_sol'
292  chckerr(err_msg)
293  end if
294 
295  ! set up and normal limits
296  kdl_x = sec_x_loc(m,2:3)
297  kdl_sol = sec_sol_loc(m,2:3)
298 
299  ! limit to grid ranges
300  kdl_x(1) = max(kdl_x(1),grid_x%i_min)
301  kdl_x(2) = min(kdl_x(2),grid_x%i_max)
302  kdl_sol(1) = max(kdl_sol(1),grid_sol%i_min)
303  kdl_sol(2) = min(kdl_sol(2),grid_sol%i_max)
304 
305  ! convert limits to local
306  kdl_x = kdl_x - grid_x%i_min + 1
307  kdl_sol = kdl_sol - grid_sol%i_min + 1
308 
309  ! cycle if mode does not exist anywhere
310  if (kdl_x(1).gt.kdl_x(2) .or. kdl_sol(1).gt.kdl_sol(2)) cycle
311 
312  ! get perturbation and solution normal ranges for this mode number
313  r_x_loc => grid_x%loc_r_F(kdl_x(1):kdl_x(2))
314  r_sol_loc => grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2))
315 
316  ! set up table indices
317  ind_x = sec_x_loc(m,4)
318  ind_sol = sec_sol_loc(m,4)
319 
320  ! Set up multiplicative factors fac_0 and fac_1 which are used in
321  ! XUQ = fac_0*X + fac_1*DX with fac_0 and fac_1:
322  ! 1. fac_0 = 1 (i(nq-m) or i(n-iota m) for deriv)
323  ! fac_1 = 0
324  ! 2. fac_0 = U_0 (DU_0 for deriv)
325  ! fac_1 = U_1 (DU_1 for deriv)
326  ! 3. fac_0 = i(nq-m)/J or i(n-iota m)/J (no deriv)
327  ! fac_1 = 0
328  ! 4. fac_0 = DU_0/J - S (no deriv)
329  ! fac_1 = DU_1/J
330  ! where par_fac = nq-m or n-iota m.
331  ! These quantities are then multiplied by the exponent for each mode
332  ! number.
333  allocate(expon(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2))) ! normal subset of grid_sol%loc_r_F
334  allocate(fac_0(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2))) ! normal subset of grid_sol%loc_r_F
335  allocate(fac_1(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2))) ! normal subset of grid_sol%loc_r_F
336  allocate(par_fac(kdl_sol(1):kdl_sol(2))) ! normal subset of grid_sol%loc_r_F
337 
338  ! calculate parallel multiplicative factor par_fac and exponent
339  do kd = kdl_sol(1),kdl_sol(2)
340  ! set up local mode numbers
341  kd_tot = kd + grid_sol%i_min - 1
342  n_loc = mds_sol%n(kd_tot,ind_sol)
343  m_loc = mds_sol%m(kd_tot,ind_sol)
344 #if ldebug
345  if (use_pol_flux_f .and. (m_loc.ne.sec_x_loc(m,1)) .or. &
346  &(.not.use_pol_flux_f .and. (n_loc.ne.sec_x_loc(m,1)))) then
347  ierr = 1
348  err_msg = 'mode numbers not consistent'
349  chckerr(err_msg)
350  end if
351 #endif
352 
353  if (use_pol_flux_f) then
354  par_fac(kd) = n_loc*jq(kd) - m_loc
355  else
356  par_fac(kd) = n_loc - m_loc*jq(kd)
357  end if
358  expon(:,:,kd) = n_loc*zeta_f(:,:,kd) - &
359  &m_loc*theta_f(:,:,kd)
360  end do
361 
362  ! set up multiplicative factors in solution grid
363  select case (xuq_style)
364  case (1) ! calculating X
365  if (deriv_loc) then ! parallel derivative
366  do kd = kdl_sol(1),kdl_sol(2)
367  fac_0(:,:,kd) = iu*par_fac(kd)
368  end do
369  else
370  fac_0 = 1._dp
371  end if
372  fac_1 = 0._dp
373  case (2) ! calculating U
374  ! interpolate X-> sol or point
375  select case (x_grid_style)
376  case (1,3) ! equilibrium or enriched
377  allocate(u0(xuq_dims(1),xuq_dims(2),&
378  &kdl_sol(1):kdl_sol(2)))
379  allocate(u1(xuq_dims(1),xuq_dims(2),&
380  &kdl_sol(1):kdl_sol(2)))
381  if (deriv_loc) then ! parallel derivative
382  ierr = interp_v_spline(&
383  &x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
384  &r_x_loc,r_sol_loc,extrap)
385  chckerr('')
386  ierr = interp_v_spline(&
387  &x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
388  &r_x_loc,r_sol_loc,extrap)
389  chckerr('')
390  else
391  ierr = interp_v_spline(&
392  &x%U_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
393  &r_x_loc,r_sol_loc,extrap)
394  chckerr('')
395  ierr = interp_v_spline(&
396  &x%U_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
397  &r_x_loc,r_sol_loc,extrap)
398  chckerr('')
399  end if
400  case (2) ! solution
401  if (deriv_loc) then ! parallel derivative
402  u0 => x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x)
403  u1 => x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x)
404  else
405  u0 => x%U_0(:,:,kdl_x(1):kdl_x(2),ind_x)
406  u1 => x%U_1(:,:,kdl_x(1):kdl_x(2),ind_x)
407  end if
408  end select
409 
410  ! assign
411  fac_0 = u0
412  fac_1 = u1
413 
414  ! clean up
415  select case (x_grid_style)
416  case (1,3) ! equilibrium or enriched
417  deallocate(u0,u1)
418  case (2) ! solution
419  ! do nothing
420  end select
421  case (3) ! calculating Qn
422  if (deriv_loc) then ! parallel derivative
423  ierr = 1
424  err_msg = 'For XUQ_style = '//trim(i2str(xuq_style))//&
425  &' calculation of parallel derivative not supported'
426  chckerr('')
427  else
428  do kd = kdl_sol(1),kdl_sol(2)
429  fac_0(:,:,kd) = iu*par_fac(kd)*inv_j(:,:,kd)
430  end do
431  end if
432  fac_1 = 0._dp
433  case (4) ! calculating Qg
434  ! interpolate X-> sol or point
435  err_msg = 'For XUQ_style = '//trim(i2str(xuq_style))//&
436  &' calculation of parallel derivative not supported'
437  chckerr('')
438  select case (x_grid_style)
439  case (1,3) ! equilibrium or enriched
440  allocate(u0(xuq_dims(1),xuq_dims(2),&
441  &kdl_sol(1):kdl_sol(2)))
442  allocate(u1(xuq_dims(1),xuq_dims(2),&
443  &kdl_sol(1):kdl_sol(2)))
444  if (deriv_loc) then ! parallel derivative
445  ierr = 1
446  chckerr(err_msg)
447  else
448  ierr = interp_v_spline(&
449  &x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
450  &r_x_loc,r_sol_loc,extrap)
451  chckerr('')
452  ierr = interp_v_spline(&
453  &x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
454  &r_x_loc,r_sol_loc,extrap)
455  chckerr('')
456  end if
457  case (2) ! solution
458  if (deriv_loc) then ! parallel derivative
459  ierr = 1
460  chckerr(err_msg)
461  else
462  u0 => x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x)
463  u1 => x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x)
464  end if
465  end select
466 
467  ! assign
468  fac_0 = u0*inv_j(:,:,kdl_sol(1):kdl_sol(2)) - &
469  &s(:,:,kdl_sol(1):kdl_sol(2))
470  fac_1 = u1*inv_j(:,:,kdl_sol(1):kdl_sol(2))
471 
472  ! clean up
473  select case (x_grid_style)
474  case (1,3) ! equilibrium or enriched
475  deallocate(u0,u1)
476  case (2) ! solution
477  ! do nothing
478  end select
479  case default
480  err_msg = 'No XUQ style associated with '//&
481  &trim(i2str(xuq_style))
482  ierr = 1
483  chckerr(err_msg)
484  end select
485 
486  ! iterate over all normal points in solution modes variables
487  do kd = kdl_sol(1),kdl_sol(2)
488  ! set up loc complex XUQ without time at this normal point
489  xuq_loc = exp(iu*expon(:,:,kd)) * &
490  &(sol%vec(ind_sol,kd,x_id)*fac_0(:,:,kd) + &
491  &dsol_vec(ind_sol,kd)*fac_1(:,:,kd))
492 
493  ! iterate over time steps
494  do id = 1,n_t
495  ! add current mode
496  xuq(:,:,kd,id) = xuq(:,:,kd,id) + xuq_loc * &
497  &exp(iu*sqrt_sol_val_norm*time(id)*2*pi)
498  end do
499  end do
500 
501  ! clean up
502  deallocate(fac_0,fac_1,expon,par_fac)
503  end do
504 
505  ! clean up
506  select case (x_grid_style)
507  case (1,3) ! equilibrium or enriched
508  deallocate(theta_f,zeta_f)
509  case (2) ! solution
510  ! do nothing
511  end select
512  nullify(r_x_loc,r_sol_loc)
513  nullify(theta_f,zeta_f)
514  nullify(u0,u1)
515  nullify(sec_x_loc,sec_sol_loc)
516  end function calc_xuq_arr
518  integer function calc_xuq_ind(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,&
519  &eq_2,X,sol,X_id,XUQ_style,time,XUQ,deriv) result(ierr)
520 
521  character(*), parameter :: rout_name = 'calc_XUQ_ind'
522 
523  ! input / output
524  type(modes_type), intent(in) :: mds_x
525  type(modes_type), intent(in) :: mds_sol
526  type(grid_type), intent(in) :: grid_eq
527  type(grid_type), intent(in) :: grid_x
528  type(grid_type), intent(in) :: grid_sol
529  type(eq_1_type), intent(in) :: eq_1
530  type(eq_2_type), intent(in) :: eq_2
531  type(x_1_type), intent(in) :: x
532  type(sol_type), intent(in) :: sol
533  integer, intent(in) :: x_id
534  integer, intent(in) :: xuq_style
535  real(dp), intent(in) :: time
536  complex(dp), intent(inout) :: xuq(:,:,:)
537  logical, intent(in), optional :: deriv
538 
539  ! local variables
540  complex(dp), allocatable :: xuq_arr(:,:,:,:) ! dummy variable to use array version
541 
542  ! allocate array version of XUQ
543  allocate(xuq_arr(size(xuq,1),size(xuq,2),size(xuq,3),1))
544 
545  ! call array version
546  ierr = calc_xuq_arr(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,eq_2,x,&
547  &sol,x_id,xuq_style,[time],xuq_arr,deriv=deriv)
548  chckerr('')
549 
550  ! copy array to individual XUQ
551  xuq = xuq_arr(:,:,:,1)
552 
553  ! deallocate array version of XUQ
554  deallocate(xuq_arr)
555  end function calc_xuq_ind
556 
579  integer function calc_tot_sol_vec(mds,grid_sol,sol_vec_loc,sol_vec_tot,&
580  &deriv) result(ierr)
581  use num_vars, only: x_style, norm_disc_prec_sol
582  use x_vars, only: n_mod_x
583  use num_utilities, only: spline
584 
585  character(*), parameter :: rout_name = 'calc_tot_sol_vec'
586 
587  ! input / output
588  type(modes_type), intent(in) :: mds
589  type(grid_type), intent(in) :: grid_sol
590  complex(dp), intent(in) :: sol_vec_loc(:,:)
591  complex(dp), intent(inout), allocatable :: sol_vec_tot(:,:)
592  integer, intent(in), optional :: deriv
593 
594  ! local variables
595  character(len=max_str_ln) :: err_msg ! error message
596  integer :: n_mod_loc ! local number of modes that are used in the calculations
597  integer :: n_mod_tot ! total number of modes
598  integer :: ld ! counter
599  integer :: kdl(2) ! kd limits
600  integer :: ld_loc ! local ld
601  integer :: n_r ! normal size of variables
602  integer :: deriv_loc ! local derivative
603  integer :: max_deriv ! maximum derivative
604  integer :: min_range ! minimum range for derivative
605  complex(dp), allocatable :: dsol_vec_loc(:) ! local Dsol_vec_tot in solution grid
606 
607  ! initialize ierr
608  ierr = 0
609 
610  ! set n_r
611  n_r = size(sol_vec_loc,2)
612 
613  ! set local derivative
614  deriv_loc = 0
615  if (present(deriv)) deriv_loc = deriv
616 
617  ! test
618  select case (norm_disc_prec_sol)
619  case (1) ! linear
620  max_deriv = 1
621  min_range = 2
622  case (2) ! Akita hermite
623  max_deriv = 1
624  min_range = 4
625  case (3) ! spline
626  max_deriv = 3
627  min_range = 4
628  case default
629  ierr = 1
630  err_msg = 'normal discretization precision for solution &
631  &cannot be higher than 3 or lower than 0'
632  chckerr(err_msg)
633  end select
634  if (deriv_loc.lt.0 .or. deriv_loc.gt.max_deriv) then
635  ierr = 1
636  err_msg = 'deriv must be between 0 and '//trim(i2str(max_deriv))//&
637  &' for normal discretization '//trim(i2str(norm_disc_prec_sol))
638  chckerr(err_msg)
639  endif
640 
641  ! operations depending on X style
642  select case (x_style)
643  case (1) ! prescribed
644  ! test allocation
645  if (allocated(sol_vec_tot)) then
646  if (size(sol_vec_tot,1).ne.size(sol_vec_loc,1) .or. &
647  &size(sol_vec_tot,2).ne.size(sol_vec_loc,2)) then
648  ierr = 1
649  err_msg = 'Local and total solution vectors need to &
650  &be of the same length'
651  chckerr(err_msg)
652  end if
653  else
654  allocate(sol_vec_tot(size(sol_vec_loc,1),n_r))
655  end if
656 
657  ! copy or derivate
658  if (deriv_loc.eq.0) then ! copy
659  sol_vec_tot = sol_vec_loc
660  else ! derivate
661  do ld = 1,size(sol_vec_tot,1)
662  ierr = spline(grid_sol%r_F,sol_vec_loc(ld,:),&
663  &grid_sol%r_F,sol_vec_tot(ld,:),&
664  &ord=norm_disc_prec_sol,deriv=deriv_loc)
665  chckerr('')
666  end do
667  end if
668  case (2) ! fast
669  ! set local and total nr. of modes
670  n_mod_loc = n_mod_x
671  n_mod_tot = size(mds%sec,1)
672 
673  ! test allocation
674  if (allocated(sol_vec_tot)) then
675  if (size(sol_vec_tot,1).ne.n_mod_tot .or. &
676  &size(sol_vec_tot,2).ne.size(sol_vec_loc,2)) then
677  ierr = 1
678  err_msg = 'Total solution vectors needs to be of &
679  &correct length '//trim(i2str(n_mod_tot))
680  chckerr(err_msg)
681  end if
682  else
683  allocate(sol_vec_tot(n_mod_tot,n_r))
684  end if
685 
686  ! track all possible modes and copy or derivate
687  sol_vec_tot = 0._dp
688  do ld = 1,size(mds%sec,1)
689  ld_loc = mds%sec(ld,1)-minval(mds%sec(:,1),1)+1
690  kdl = [mds%sec(ld,2),mds%sec(ld,3)]-grid_sol%i_min+1 ! normal range in solution vector
691  kdl = max(1,min(kdl,n_r)) ! limit to sol_vec_loc range
692  if (deriv_loc.eq.0) then
693  sol_vec_tot(ld_loc,kdl(1):kdl(2)) = &
694  &sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2))
695  else
696  allocate(dsol_vec_loc(kdl(2)-kdl(1)+1))
697  dsol_vec_loc = 0._dp
698 
699  if (kdl(2)-kdl(1)+1.ge.min_range) then ! only calculate if enough points
700  ierr = spline(grid_sol%r_F(kdl(1):kdl(2)),&
701  &sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2)),&
702  &grid_sol%r_F(kdl(1):kdl(2)),dsol_vec_loc,&
703  &ord=norm_disc_prec_sol,deriv=deriv_loc)
704  chckerr('')
705  end if
706  sol_vec_tot(ld_loc,kdl(1):kdl(2)) = dsol_vec_loc
707 
708  deallocate(dsol_vec_loc)
709  end if
710  end do
711  end select
712  end function calc_tot_sol_vec
713 
720  integer function calc_loc_sol_vec(mds,i_min,sol_vec_tot,sol_vec_loc) &
721  &result(ierr)
722  use num_vars, only: x_style
723  use x_vars, only: n_mod_x
724 
725  character(*), parameter :: rout_name = 'calc_loc_sol_vec'
726 
727  ! input / output
728  type(modes_type), intent(in) :: mds
729  integer, intent(in) :: i_min
730  complex(dp), intent(in) :: sol_vec_tot(:,:)
731  complex(dp), intent(inout), allocatable :: sol_vec_loc(:,:)
732 
733  ! local variables
734  character(len=max_str_ln) :: err_msg ! error message
735  integer :: n_mod_loc ! local number of modes that are used in the calculations
736  integer :: n_mod_tot ! total number of modes
737  integer :: ld ! counter
738  integer :: kdl(2) ! kd limits
739  integer :: ld_loc ! local ld
740  integer :: n_r ! normal size of variables
741 
742  ! set n_r
743  n_r = size(sol_vec_tot,2)
744 
745  ! initialize ierr
746  ierr = 0
747 
748  ! operations depending on X style
749  select case (x_style)
750  case (1) ! prescribed
751  ! test allocation
752  if (allocated(sol_vec_loc)) then
753  if (size(sol_vec_loc,1).ne.size(sol_vec_tot,1) .or. &
754  &size(sol_vec_loc,2).ne.size(sol_vec_tot,2)) then
755  ierr = 1
756  err_msg = 'Local and total solution vectors need to &
757  &be of the same length'
758  chckerr(err_msg)
759  end if
760  else
761  allocate(sol_vec_loc(size(sol_vec_tot,1),n_r))
762  end if
763  sol_vec_loc = sol_vec_tot
764  case (2) ! fast
765  ! set local and total nr. of modes
766  n_mod_loc = n_mod_x
767  n_mod_tot = size(mds%sec,1)
768 
769  ! test allocation
770  if (allocated(sol_vec_loc)) then
771  if (size(sol_vec_loc,1).ne.n_mod_loc .or. &
772  &size(sol_vec_loc,2).ne.size(sol_vec_tot,2)) then
773  ierr = 1
774  err_msg = 'Total solution vectors needs to be of &
775  &correct length '//trim(i2str(n_mod_loc))
776  chckerr(err_msg)
777  end if
778  else
779  allocate(sol_vec_loc(n_mod_loc,n_r))
780  end if
781 
782  ! track all possible modes and copy
783  sol_vec_loc = 0._dp
784  do ld = 1,size(mds%sec,1)
785  ld_loc = mds%sec(ld,1)-minval(mds%sec(:,1),1)+1
786  kdl = [mds%sec(ld,2),mds%sec(ld,3)]-i_min+1 ! range in sol_vec_loc
787  kdl = max(1,min(kdl,n_r)) ! limit to sol_vec_loc range
788  sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2)) = &
789  &sol_vec_tot(ld_loc,kdl(1):kdl(2))
790  end do
791  end select
792  end function calc_loc_sol_vec
793 
799  integer function interp_v_spline(V_i,V_o,r_i,r_o,extrap,ivs_stat) &
800  &result(ierr)
801 
802  use num_utilities, only: spline
803 
804  character(*), parameter :: rout_name = 'interp_V_spline'
805 
806  ! input / output
807  complex(dp), intent(in) :: v_i(:,:,:) ! input variable
808  complex(dp), intent(out) :: v_o(:,:,:) ! output variable
809  real(dp), intent(in) :: r_i(:) ! input r
810  real(dp), intent(in) :: r_o(:) ! output r
811  logical, intent(in) :: extrap ! extrapolation possible
812  integer, intent(out), optional :: ivs_stat ! which method was used (1: copy, 2: linear, 3: spline)
813 
814  ! local variables
815  integer :: id, jd, kd ! counters
816  character(len=max_str_ln) :: err_msg ! error message
817 
818  ! initialize ierr
819  ierr = 0
820 
821 #if ldebug
822  ! test sizes
823  if (size(r_i).ne.size(v_i,3)) then
824  ierr = 1
825  err_msg = 'r_i and V_i are not compatible'
826  chckerr(err_msg)
827  end if
828  if (size(r_o).ne.size(v_o,3)) then
829  ierr = 1
830  err_msg = 'r_o and V_o are not compatible'
831  chckerr(err_msg)
832  end if
833  if (size(v_i,1).ne.size(v_o,1) .or. size(v_i,2).ne.size(v_o,2)) then
834  ierr = 1
835  err_msg = 'V_i and V_o are not compatible'
836  chckerr(err_msg)
837  end if
838 #endif
839 
840  ! interpolate depending on normal size
841  select case (size(r_i))
842  case (1) ! copy directly
843  do kd = 1,size(v_o,3)
844  v_o(:,:,kd) = v_i(:,:,1)
845  end do
846  if (present(ivs_stat)) ivs_stat = 1
847  case (2:3) ! linear
848  do jd = 1,size(v_i,2)
849  do id = 1,size(v_i,1)
850  ierr = spline(r_i,v_i(id,jd,:),r_o,v_o(id,jd,:),&
851  &ord=1,extrap=extrap)
852  chckerr('')
853  end do
854  end do
855  if (present(ivs_stat)) ivs_stat = 2
856  case (4:) ! spline
857  do jd = 1,size(v_i,2)
858  do id = 1,size(v_i,1)
859  ierr = spline(r_i,v_i(id,jd,:),r_o,v_o(id,jd,:),&
860  &ord=3,extrap=extrap)
861  chckerr('')
862  end do
863  end do
864  if (present(ivs_stat)) ivs_stat = 3
865 #if ldebug
866  case default
867  ierr = 1
868  chckerr('Need more than 0 points')
869 #endif
870  end select
871 
872 #if ldebug
873  ! visualize
874  do jd = 1,size(v_i,2)
875  do id = 1,size(v_i,1)
876  if (debug_interp_v_spline) then
877  call plot_comp('real part',r_i,rp(v_i(id,jd,:)),&
878  &r_o,rp(v_o(id,jd,:)))
879  call plot_comp('imaginary part',r_i,ip(v_i(id,jd,:)),&
880  &r_o,ip(v_o(id,jd,:)))
881  end if
882  end do
883  end do
884  contains
886  subroutine plot_comp(plot_name,x_i,y_i,x_o,y_o)
887  ! input / output
888  character(len=*), intent(in) :: plot_name ! name of plot
889  real(dp), intent(in) :: x_i(:) ! x input
890  real(dp), intent(in) :: y_i(:) ! y input
891  real(dp), intent(in) :: x_o(:) ! x output
892  real(dp), intent(in) :: y_o(:) ! y output
893 
894  ! local variables
895  real(dp), allocatable :: x_plot(:,:) ! x for plot
896  real(dp), allocatable :: y_plot(:,:) ! y for plot
897  character(len=max_str_ln) :: var_names(2) ! variable names
898 
899  allocate(x_plot(max(size(x_i),size(x_o)),2))
900  allocate(y_plot(max(size(x_i),size(x_o)),2))
901 
902  x_plot(1:size(x_i),1) = x_i
903  x_plot(size(x_i):size(x_plot,1),1) = x_i(size(x_i))
904  x_plot(1:size(x_o),2) = x_o
905  x_plot(size(x_o):size(x_plot,1),2) = x_o(size(x_o))
906 
907  y_plot(1:size(x_i),1) = y_i(:)
908  y_plot(size(x_i):size(y_plot,1),1) = y_i(size(x_i))
909  y_plot(1:size(x_o),2) = y_o(:)
910  y_plot(size(x_o):size(y_plot,1),2) = y_o(size(x_o))
911 
912  var_names(1) = 'input '//plot_name
913  var_names(2) = 'output '//plot_name
914  call print_ex_2d(var_names,'',y_plot,x=x_plot)
915  end subroutine plot_comp
916 #endif
917  end function interp_v_spline
918 end module sol_utilities
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_utilities::calc_int
Integrates a function using the trapezoidal rule.
Definition: num_utilities.f90:160
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
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_utilities::con2dis
Convert between points from a continuous grid to a discrete grid.
Definition: num_utilities.f90:188
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
num_vars::norm_disc_prec_sol
integer, public norm_disc_prec_sol
precision for normal discretization for solution
Definition: num_vars.f90:122
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
num_vars::x_style
integer, public x_style
style for secondary mode numbers (1: prescribed, 2: fast)
Definition: num_vars.f90:95
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
grid_vars::grid_type
Type for grids.
Definition: grid_vars.f90:59
eq_vars::eq_1_type
flux equilibrium type
Definition: eq_vars.f90:63
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
sol_utilities::debug_calc_xuq_arr
logical, public debug_calc_xuq_arr
plot debug information for calc_XUQ_arr
Definition: sol_utilities.f90:25
sol_utilities::calc_loc_sol_vec
integer function, public calc_loc_sol_vec(mds, i_min, sol_vec_tot, sol_vec_loc)
Calculate the local version of the solution vector from the total version.
Definition: sol_utilities.f90:722
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
x_vars
Variables pertaining to the perturbation quantities.
Definition: X_vars.f90:4
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
x_utilities
Numerical utilities related to perturbation operations.
Definition: X_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
x_vars::modes_type
mode number type
Definition: X_vars.f90:36
sol_utilities::interp_v_spline
integer function, public interp_v_spline(V_i, V_o, r_i, r_o, extrap, ivs_stat)
Interpolation for a quantity V using splines.
Definition: sol_utilities.f90:801
sol_utilities::calc_xuq
Calculates the normal ( ) or geodesic ( ) components of the plasma perturbation or the magnetic pert...
Definition: sol_utilities.f90:56
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
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_vars
Variables pertaining to the different grids used.
Definition: grid_vars.f90:4
x_vars::n_mod_x
integer, public n_mod_x
size of m_X (pol. flux) or n_X (tor. flux)
Definition: X_vars.f90:129
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
sol_utilities::debug_interp_v_spline
logical, public debug_interp_v_spline
debug information for interp_v_spline
Definition: sol_utilities.f90:26
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
x_utilities::trim_modes
integer function, public trim_modes(mds_i, mds_o, id_lim_i, id_lim_o)
Limit input mode range to output mode range.
Definition: X_utilities.f90:293
eq_vars::eq_2_type
metric equilibrium type
Definition: eq_vars.f90:114
this documentation was created with Doxygen 1.8.17 by Toon Weyens