PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
VMEC_utilities.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module vmec_utilities
5 #include <PB3D_macros.h>
6  use str_utilities
7  use output_ops
8  use messages
9  use num_vars, only: &
10  &dp, max_str_ln, pi
11  use vmec_vars
12 
13  implicit none
14  private
15  public calc_trigon_factors, fourier2real
16 #if ldebug
17  public debug_calc_trigon_factors
18 #endif
19 
20  ! global variables
21 #if ldebug
22  logical :: debug_calc_trigon_factors = .false.
23 #endif
24 
25  ! interfaces
26 
56  interface fourier2real
58  module procedure fourier2real_1
60  module procedure fourier2real_2
61  end interface
62 
63 contains
65  integer function fourier2real_1(varf_c,varf_s,trigon_factors,varr,sym,&
66  &deriv) result(ierr)
67 
68  character(*), parameter :: rout_name = 'fourier2real_1'
69 
70  ! input / output
71  real(dp), intent(in) :: varf_c(:,:)
72  real(dp), intent(in) :: varf_s(:,:)
73  real(dp), intent(in) :: trigon_factors(:,:,:,:,:)
74  real(dp), intent(inout) :: varr(:,:,:)
75  logical, intent(in), optional :: sym(2)
76  integer, intent(in), optional :: deriv(2)
77 
78  ! local variables
79  character(len=max_str_ln) :: err_msg ! error message
80  integer :: dims(3) ! dimensions of varr
81  integer :: id, kd ! counters
82  integer :: deriv_loc(2) ! local derivative
83  logical :: sym_loc(2) ! local sym
84  real(dp) :: deriv_fac ! factur due to derivatives
85 
86  ! initialize ierr
87  ierr = 0
88 
89  ! set local deriv and sym
90  deriv_loc = [0,0]
91  sym_loc = [.true.,.true.]
92  if (present(deriv)) deriv_loc = deriv
93  if (present(sym)) sym_loc = sym
94 
95  ! set dimensions
96  dims = shape(varr)
97 
98  ! test
99  if (.not.sym_loc(1) .and. .not.sym_loc(2)) then
100  ierr = 1
101  err_msg = 'Need at least the cosine or the sine factor'
102  chckerr(err_msg)
103  end if
104  if (size(trigon_factors,1).ne.mnmax_v .and. &
105  &size(trigon_factors,2).ne.dims(1) .and. &
106  &size(trigon_factors,3).ne.dims(2) .and. &
107  &size(trigon_factors,4).ne.dims(3) .and. &
108  &size(trigon_factors,5).ne.2) then
109  ierr = 1
110  err_msg = 'Trigonometric factors are not set up'
111  chckerr(err_msg)
112  end if
113 
114  ! initialize
115  varr = 0.0_dp
116 
117  ! sum over modes
118  do id = 1,mnmax_v
119  ! setup derivative factor for varf_c
120  deriv_fac = mn_v(id,1)**deriv_loc(1)*(-mn_v(id,2))**deriv_loc(2)*&
121  &(-1)**((sum(deriv_loc)+1)/2)
122 
123  ! add terms ~ varf_c
124  if (sym_loc(1)) then
125  if (mod(sum(deriv_loc),2).eq.0) then ! even number of derivatives
126  do kd = 1,dims(3)
127  varr(:,:,kd) = varr(:,:,kd) + &
128  &varf_c(id,kd) * deriv_fac * &
129  &trigon_factors(id,:,:,kd,1)
130  end do
131  else ! odd number of derivatives
132  do kd = 1,dims(3)
133  varr(:,:,kd) = varr(:,:,kd) + &
134  &varf_c(id,kd) * deriv_fac * &
135  &trigon_factors(id,:,:,kd,2)
136  end do
137  end if
138  end if
139 
140  ! setup derivative factor for varf_s
141  deriv_fac = mn_v(id,1)**deriv_loc(1)*(-mn_v(id,2))**deriv_loc(2)*&
142  &(-1)**(sum(deriv_loc)/2)
143 
144  ! add terms ~ varf_s
145  if (sym_loc(2)) then
146  if (mod(sum(deriv_loc),2).eq.0) then ! even number of derivatives
147  do kd = 1,dims(3)
148  varr(:,:,kd) = varr(:,:,kd) + &
149  &varf_s(id,kd) * deriv_fac * &
150  &trigon_factors(id,:,:,kd,2)
151  end do
152  else ! odd number of derivatives
153  do kd = 1,dims(3)
154  varr(:,:,kd) = varr(:,:,kd) + &
155  &varf_s(id,kd) * deriv_fac * &
156  &trigon_factors(id,:,:,kd,1)
157  end do
158  end if
159  end if
160  end do
161  end function fourier2real_1
163  integer function fourier2real_2(varf_c,varf_s,theta,zeta,varr,sym,deriv) &
164  &result(ierr)
165 
166  character(*), parameter :: rout_name = 'fourier2real_2'
167 
168  ! input / output
169  real(dp), intent(in) :: varf_c(:,:)
170  real(dp), intent(in) :: varf_s(:,:)
171  real(dp), intent(in) :: theta(:,:,:)
172  real(dp), intent(in) :: zeta(:,:,:)
173  real(dp), intent(inout) :: varr(:,:,:)
174  logical, intent(in), optional :: sym(2)
175  integer, intent(in), optional :: deriv(2)
176 
177 
178  ! local variables
179  character(len=max_str_ln) :: err_msg ! error message
180  integer :: dims(3) ! dimensions of varr
181  integer :: id, kd ! counters
182  integer :: deriv_loc(2) ! local derivative
183  logical :: sym_loc(2) ! local sym
184  real(dp) :: deriv_fac ! factur due to derivatives
185  real(dp), allocatable :: ang(:,:,:) ! angle of trigonometric functions
186 
187  ! initialize ierr
188  ierr = 0
189 
190  ! set local deriv and sym
191  deriv_loc = [0,0]
192  sym_loc = [.true.,.true.]
193  if (present(deriv)) deriv_loc = deriv
194  if (present(sym)) sym_loc = sym
195 
196  ! test
197  if (.not.sym_loc(1) .and. .not.sym_loc(2)) then
198  ierr = 1
199  err_msg = 'Need at least the cosine or the sine factor'
200  chckerr(err_msg)
201  end if
202 
203  ! set dimensions
204  dims = shape(varr)
205 
206  ! initialize angle
207  allocate(ang(dims(1),dims(2),dims(3)))
208 
209  ! initialize output
210  varr = 0.0_dp
211 
212  ! sum over modes
213  do id = 1,mnmax_v
214  ! set angle
215  ang = mn_v(id,1)*theta - mn_v(id,2)*zeta
216 
217  ! setup derivative factor for varf_c
218  deriv_fac = mn_v(id,1)**deriv_loc(1)*(-mn_v(id,2))**deriv_loc(2)*&
219  &(-1)**((sum(deriv_loc)+1)/2)
220 
221  ! add terms ~ varf_c
222  if (sym_loc(1)) then
223  if (mod(sum(deriv_loc),2).eq.0) then ! even number of derivatives
224  do kd = 1,dims(3)
225  varr(:,:,kd) = varr(:,:,kd) + &
226  &varf_c(id,kd) * deriv_fac * cos(ang(:,:,kd))
227  end do
228  else ! odd number of derivatives
229  do kd = 1,dims(3)
230  varr(:,:,kd) = varr(:,:,kd) + &
231  &varf_c(id,kd) * deriv_fac * sin(ang(:,:,kd))
232  end do
233  end if
234  end if
235 
236  ! setup derivative factor for varf_s
237  deriv_fac = mn_v(id,1)**deriv_loc(1)*(-mn_v(id,2))**deriv_loc(2)*&
238  &(-1)**(sum(deriv_loc)/2)
239 
240  ! add terms ~ varf_s
241  if (sym_loc(2)) then
242  if (mod(sum(deriv_loc),2).eq.0) then ! even number of derivatives
243  do kd = 1,dims(3)
244  varr(:,:,kd) = varr(:,:,kd) + &
245  &varf_s(id,kd) * deriv_fac * sin(ang(:,:,kd))
246  end do
247  else ! odd number of derivatives
248  do kd = 1,dims(3)
249  varr(:,:,kd) = varr(:,:,kd) + &
250  &varf_s(id,kd) * deriv_fac * cos(ang(:,:,kd))
251  end do
252  end if
253  end if
254  end do
255  end function fourier2real_2
256 
273  integer function calc_trigon_factors(theta,zeta,trigon_factors) &
274  &result(ierr)
275 
276  character(*), parameter :: rout_name = 'calc_trigon_factors'
277 
278  ! input / output
279  real(dp), intent(in) :: theta(:,:,:)
280  real(dp), intent(in) :: zeta(:,:,:)
281  real(dp), intent(inout), allocatable :: trigon_factors(:,:,:,:,:)
282 
283  ! local variables
284  character(len=max_str_ln) :: err_msg ! error message
285  integer :: n_ang_1, n_ang_2, n_r ! sizes of 3D real output array
286  real(dp), allocatable :: cos_theta(:,:,:,:) ! cos(m theta) for all m
287  real(dp), allocatable :: sin_theta(:,:,:,:) ! sin(m theta) for all m
288  real(dp), allocatable :: cos_zeta(:,:,:,:) ! cos(n theta) for all n
289  real(dp), allocatable :: sin_zeta(:,:,:,:) ! sin(n theta) for all n
290  integer :: id, m, n ! counters
291 
292  ! initialize ierr
293  ierr = 0
294 
295  ! set n_ang_1 and n_r
296  n_ang_1 = size(theta,1)
297  n_ang_2 = size(theta,2)
298  n_r = size(theta,3)
299 
300 #if ldebug
301  if (debug_calc_trigon_factors) then
302  call writo('Calculate trigonometric factors for grid of size ['//&
303  &trim(i2str(n_ang_1))//','//trim(i2str(n_ang_2))//','//&
304  &trim(i2str(n_r))//']')
305  end if
306 #endif
307 
308  ! tests
309  if (size(zeta,1).ne.n_ang_1 .or. size(zeta,2).ne.n_ang_2 .or. &
310  &size(zeta,3).ne.n_r) then
311  ierr = 1
312  err_msg = 'theta and zeta need to have the same size'
313  chckerr(err_msg)
314  end if
315 
316  ! setup cos_theta, sin_theta, cos_zeta and sin_zeta
317  allocate(cos_theta(0:mpol_v-1,n_ang_1,n_ang_2,n_r))
318  allocate(sin_theta(0:mpol_v-1,n_ang_1,n_ang_2,n_r))
319  allocate(cos_zeta(-ntor_v:ntor_v,n_ang_1,n_ang_2,n_r))
320  allocate(sin_zeta(-ntor_v:ntor_v,n_ang_1,n_ang_2,n_r))
321 
322  do m = 0,mpol_v-1
323 #if ldebug
324  if (debug_calc_trigon_factors) then
325  call writo('Calculating for m = '//trim(i2str(m))//&
326  &' of [0..'//trim(i2str(mpol_v-1))//']')
327  end if
328 #endif
329  cos_theta(m,:,:,:) = cos(m*theta)
330  sin_theta(m,:,:,:) = sin(m*theta)
331  end do
332  do n = -ntor_v,ntor_v
333 #if ldebug
334  if (debug_calc_trigon_factors) then
335  call writo('Calculating for n = '//trim(i2str(n))//&
336  &' of ['//trim(i2str(-ntor_v))//'..'//&
337  &trim(i2str(ntor_v))//']')
338  end if
339 #endif
340  cos_zeta(n,:,:,:) = cos(n*nfp_v*zeta)
341  sin_zeta(n,:,:,:) = sin(n*nfp_v*zeta)
342  end do
343 
344  ! initialize trigon_factors
345  allocate(trigon_factors(mnmax_v,n_ang_1,n_ang_2,n_r,2),stat=ierr)
346  chckerr('Failed to allocate trigonometric factors')
347  trigon_factors = 0.0_dp
348 
349  ! calculate cos(m theta - n zeta) =
350  ! cos(m theta) cos(n zeta) + sin(m theta) sin(n zeta)
351  ! and sin(m theta - n zeta) =
352  ! sin(m theta) cos(n zeta) - cos(m theta) sin(n zeta)
353  ! Note: need to scale the indices mn_V(:,2) by nfp_V.
354  do id = 1,mnmax_v
355 #if ldebug
356  if (debug_calc_trigon_factors) then
357  call writo('Calculating for id = '//trim(i2str(id))//&
358  &' of [1..'//trim(i2str(mnmax_v))//']')
359  end if
360 #endif
361  trigon_factors(id,:,:,:,1) = &
362  &cos_theta(mn_v(id,1),:,:,:)*&
363  &cos_zeta(mn_v(id,2)/nfp_v,:,:,:) + &
364  &sin_theta(mn_v(id,1),:,:,:)*&
365  &sin_zeta(mn_v(id,2)/nfp_v,:,:,:)
366  trigon_factors(id,:,:,:,2) = &
367  &sin_theta(mn_v(id,1),:,:,:)*&
368  &cos_zeta(mn_v(id,2)/nfp_v,:,:,:) - &
369  &cos_theta(mn_v(id,1),:,:,:)*&
370  &sin_zeta(mn_v(id,2)/nfp_v,:,:,:)
371  end do
372 
373  ! deallocate variables
374  deallocate(cos_theta,sin_theta)
375  deallocate(cos_zeta,sin_zeta)
376  end function calc_trigon_factors
377 end module vmec_utilities
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
vmec_vars::mn_v
integer, dimension(:,:), allocatable, public mn_v
m and n of modes
Definition: VMEC_vars.f90:32
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
vmec_utilities::calc_trigon_factors
integer function, public calc_trigon_factors(theta, zeta, trigon_factors)
Calculate the trigonometric cosine and sine factors.
Definition: VMEC_utilities.f90:275
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
vmec_utilities::debug_calc_trigon_factors
logical, public debug_calc_trigon_factors
plot debug information for calc_trigon_factors()
Definition: VMEC_utilities.f90:22
vmec_utilities
Numerical utilities related to the output of VMEC.
Definition: VMEC_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
vmec_vars
Variables that concern the output of VMEC.
Definition: VMEC_vars.f90:4
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5
vmec_utilities::fourier2real
Inverse Fourier transformation, from VMEC.
Definition: VMEC_utilities.f90:56
this documentation was created with Doxygen 1.8.17 by Toon Weyens