V3FIT
fourier.f90
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
12 !*******************************************************************************
13  MODULE fourier
14  USE stel_kinds
15  USE stel_constants
16  USE timer_mod
17  USE descriptor_mod, ONLY: inhessian, diagonaldone
18 
19  IMPLICIT NONE
20 
21 !*******************************************************************************
22 ! fourier module parameters
23 !*******************************************************************************
25  INTEGER, PARAMETER :: f_cos = 0
27  INTEGER, PARAMETER :: f_sin = 1
28 
29 ! Bit locations for fourier control flags.
32  INTEGER, PARAMETER :: f_none = 0 !0000
34  INTEGER, PARAMETER :: f_ds = 1 !0001
36  INTEGER, PARAMETER :: f_du = 2 !0010
38  INTEGER, PARAMETER :: f_dv = 4 !0100
40  INTEGER, PARAMETER :: f_dudv = ior(f_du, f_dv) !0110
43  INTEGER, PARAMETER :: f_sum = 8 !1000
44 
45 ! Test bit positions for Derivative flags.
47  INTEGER, PARAMETER :: b_ds = 0
49  INTEGER, PARAMETER :: b_du = 1
51  INTEGER, PARAMETER :: b_dv = 2
53  INTEGER, PARAMETER :: b_sum = 3
54 
55 ! Special Fourier modes.
56 ! Poloidal modes.
58  INTEGER, PARAMETER :: m0 = 0
60  INTEGER, PARAMETER :: m1 = 1
62  INTEGER, PARAMETER :: m2 = 2
63 ! Toroidal modes.
65  INTEGER, PARAMETER :: n0 = 0
67  INTEGER, PARAMETER :: n1 = 1
68 
69 !*******************************************************************************
70 ! DERIVED-TYPE DECLARATIONS
71 ! 1) fourier base class
72 !
73 !*******************************************************************************
74 !-------------------------------------------------------------------------------
76 !-------------------------------------------------------------------------------
77  TYPE :: fourier_class
79  REAL (dp), DIMENSION(:,:), POINTER :: orthonorm => null()
81  REAL (dp), DIMENSION(:,:), POINTER :: cosmu => null()
83  REAL (dp), DIMENSION(:,:), POINTER :: cosmum => null()
85  REAL (dp), DIMENSION(:,:), POINTER :: cosmui => null()
87  REAL (dp), DIMENSION(:,:), POINTER :: cosnv => null()
89  REAL (dp), DIMENSION(:,:), POINTER :: cosnvn => null()
91  REAL (dp), DIMENSION(:,:), POINTER :: sinmu => null()
93  REAL (dp), DIMENSION(:,:), POINTER :: sinmum => null()
95  REAL (dp), DIMENSION(:,:), POINTER :: sinmui => null()
97  REAL (dp), DIMENSION(:,:), POINTER :: sinnv => null()
99  REAL (dp), DIMENSION(:,:), POINTER :: sinnvn => null()
101  REAL (dp), DIMENSION(:,:), POINTER :: workmj1 => null()
103  REAL (dp), DIMENSION(:,:), POINTER :: workmj2 => null()
105  REAL (dp), DIMENSION(:,:), POINTER :: workin1 => null()
107  REAL (dp), DIMENSION(:,:), POINTER :: workin2 => null()
108  CONTAINS
109  PROCEDURE, pass :: toijsp_2d => fourier_toijsp_2d
110  PROCEDURE, pass :: toijsp_3d => fourier_toijsp_3d
111  generic :: toijsp => toijsp_2d, toijsp_3d
112  PROCEDURE, pass :: tomnsp_3d => fourier_tomnsp_3d
113  PROCEDURE, pass :: tomnsp_2d => fourier_tomnsp_2d
114  PROCEDURE, pass :: tomnsp_2d_u => fourier_tomnsp_2d_u
115  PROCEDURE, pass :: tomnsp_2d_v => fourier_tomnsp_2d_v
116  PROCEDURE, pass :: tomnsp_3d_pest => fourier_tomnsp_3d_pest
117  PROCEDURE, pass :: tomnsp_2d_pest => fourier_tomnsp_2d_pest
118  PROCEDURE, pass :: tomnsp_2d_u_pest => fourier_tomnsp_2d_u_pest
119  generic :: tomnsp => tomnsp_2d, tomnsp_3d, tomnsp_2d_u, &
120  tomnsp_2d_v, tomnsp_3d_pest, &
121  tomnsp_2d_pest, tomnsp_2d_u_pest
122  PROCEDURE :: get_origin => fourier_get_origin
123  final :: fourier_destruct
124  END TYPE
125 
126 !*******************************************************************************
127 ! INTERFACE BLOCKS
128 !*******************************************************************************
129 !-------------------------------------------------------------------------------
131 !-------------------------------------------------------------------------------
132  INTERFACE fourier_class
133  MODULE PROCEDURE fourier_construct
134  END INTERFACE
135 
136 !-------------------------------------------------------------------------------
138 !-------------------------------------------------------------------------------
139  INTERFACE check
140  MODULE PROCEDURE check_all, &
141  & check_mn
142  END INTERFACE
143 
144  PRIVATE :: check, check_all, check_mn
145 
146  CONTAINS
147 
148 !*******************************************************************************
149 ! CONSTRUCTION SUBROUTINES
150 !*******************************************************************************
151 !-------------------------------------------------------------------------------
178 !-------------------------------------------------------------------------------
179  FUNCTION fourier_construct(mpol, ntor, ntheta, nzeta, nfp, sym)
180 
181  IMPLICIT NONE
182 
183 ! Declare Arguments
184  class(fourier_class), POINTER :: fourier_construct
185  INTEGER, INTENT(in) :: mpol
186  INTEGER, INTENT(in) :: ntor
187  INTEGER, INTENT(in) :: ntheta
188  INTEGER, INTENT(in) :: nzeta
189  INTEGER, INTENT(in) :: nfp
190  LOGICAL, INTENT(in) :: sym
191 
192 ! Local Variables
193  REAL (dp) :: pi_norm
194  INTEGER :: i
195  INTEGER :: j
196  REAL (dp) :: arg
197  REAL (dp) :: dnorm
198 
199 ! Local Parameters
200  REAL (dp), PARAMETER :: nfactor = 1.41421356237310_dp
201 
202 ! Start of executable code.
203  ALLOCATE(fourier_construct)
204 
205  IF (sym) THEN
206  dnorm = one/(ntheta*nzeta)
207  pi_norm = twopi/ntheta
208  ELSE
209  dnorm = one/((ntheta - 1)*nzeta)
210  pi_norm = pi/(ntheta - 1)
211  END IF
212 
213  ALLOCATE(fourier_construct%cosmu(ntheta,0:mpol))
214  ALLOCATE(fourier_construct%cosmum(ntheta,0:mpol))
215  ALLOCATE(fourier_construct%cosmui(ntheta,0:mpol))
216  ALLOCATE(fourier_construct%sinmu(ntheta,0:mpol))
217  ALLOCATE(fourier_construct%sinmum(ntheta,0:mpol))
218  ALLOCATE(fourier_construct%sinmui(ntheta,0:mpol))
219 
220  DO i = 1, ntheta
221 ! Poloidal collocation points*m.
222  arg = pi_norm*(i - 1)
223 
224  DO j = 0, mpol
225  fourier_construct%cosmu(i, j) = cos(j*arg)
226  fourier_construct%sinmu(i, j) = sin(j*arg)
227 
228  CALL fourier_round_trig(fourier_construct%cosmu(i, j), &
229  & fourier_construct%sinmu(i, j))
230 
231  fourier_construct%cosmum(i, j) = j*fourier_construct%cosmu(i, j)
232  fourier_construct%sinmum(i, j) = j*fourier_construct%sinmu(i, j)
233 
234  fourier_construct%cosmui(i, j) = dnorm*fourier_construct%cosmu(i, j)
235  fourier_construct%sinmui(i, j) = dnorm*fourier_construct%sinmu(i, j)
236 
237  IF (.not.sym .and. (i .eq. 1 .or. i .eq. ntheta)) THEN
238  fourier_construct%cosmui(i, j) = 0.5 &
239  & * fourier_construct%cosmui(i, j)
240  fourier_construct%sinmui(i, j) = 0.5 &
241  & * fourier_construct%sinmui(i, j)
242  END IF
243  END DO
244  END DO
245 
246  ALLOCATE(fourier_construct%cosnv(nzeta,0:ntor))
247  ALLOCATE(fourier_construct%cosnvn(nzeta,0:ntor))
248  ALLOCATE(fourier_construct%sinnv(nzeta,0:ntor))
249  ALLOCATE(fourier_construct%sinnvn(nzeta,0:ntor))
250 
251 ! For now we are in one field period.
252  DO j = 0, ntor
253  DO i = 1, nzeta
254  arg = (twopi*(i - 1))/nzeta
255 
256  fourier_construct%cosnv(i,j) = cos(j*arg)
257  fourier_construct%sinnv(i,j) = sin(j*arg)
258 
259  CALL fourier_round_trig(fourier_construct%cosnv(i,j), &
260  & fourier_construct%sinnv(i,j))
261  END DO
262 
263  fourier_construct%cosnvn(:,j) = j*nfp*fourier_construct%cosnv(:,j)
264  fourier_construct%sinnvn(:,j) = j*nfp*fourier_construct%sinnv(:,j)
265  END DO
266 
267 ! Compute orthonorm factor for cos(mu + nv), sin(mu + nv) basis
268 ! (Note cos(mu)cos(nv) basis!) so that
269 !
270 ! <cos(mu + nv)*cos(m'u + n'v)>*orthonorm(m,n) = 1 (for m=m',n=n')
271 !
272 ! Orthonormal factor is Sqrt(2) for all modes except (0,0) mode. The (0,0) mode
273 ! is 1.
274  ALLOCATE(fourier_construct%orthonorm(0:mpol,-ntor:ntor))
275  fourier_construct%orthonorm = nfactor
276  fourier_construct%orthonorm(m0,n0) = 1
277  IF (ntheta .eq. mpol + 1) THEN
278  fourier_construct%orthonorm(mpol,:) = 1
279  END IF
280 
281 ! Allocate working buffers.
282  ALLOCATE(fourier_construct%workmj1(0:mpol,nzeta))
283  ALLOCATE(fourier_construct%workmj2(0:mpol,nzeta))
284  ALLOCATE(fourier_construct%workin1(ntheta,-ntor:ntor))
285  ALLOCATE(fourier_construct%workin2(ntheta,-ntor:ntor))
286 
287  END FUNCTION
288 
289 !*******************************************************************************
290 ! DESTRUCTION SUBROUTINES
291 !*******************************************************************************
292 !-------------------------------------------------------------------------------
298 !-------------------------------------------------------------------------------
299  SUBROUTINE fourier_destruct(this)
300 
301  IMPLICIT NONE
302 
303 ! Declare Arguments
304  TYPE (fourier_class), INTENT(inout) :: this
305 
306 ! Start of executable code.
307  IF (ASSOCIATED(this%orthonorm)) THEN
308  DEALLOCATE(this%orthonorm)
309  this%orthonorm => null()
310  END IF
311 
312  IF (ASSOCIATED(this%cosmu)) THEN
313  DEALLOCATE(this%cosmu)
314  this%cosmu => null()
315  END IF
316 
317  IF (ASSOCIATED(this%cosmum)) THEN
318  DEALLOCATE(this%cosmum)
319  this%cosmum => null()
320  END IF
321 
322  IF (ASSOCIATED(this%cosmui)) THEN
323  DEALLOCATE(this%cosmui)
324  this%cosmui => null()
325  END IF
326 
327  IF (ASSOCIATED(this%cosnv)) THEN
328  DEALLOCATE(this%cosnv)
329  this%cosnv => null()
330  END IF
331 
332  IF (ASSOCIATED(this%cosnvn)) THEN
333  DEALLOCATE(this%cosnvn)
334  this%cosnvn => null()
335  END IF
336 
337  IF (ASSOCIATED(this%sinmu)) THEN
338  DEALLOCATE(this%sinmu)
339  this%sinmu => null()
340  END IF
341 
342  IF (ASSOCIATED(this%sinmum)) THEN
343  DEALLOCATE(this%sinmum)
344  this%sinmum => null()
345  END IF
346 
347  IF (ASSOCIATED(this%sinmui)) THEN
348  DEALLOCATE(this%sinmui)
349  this%sinmui => null()
350  END IF
351 
352  IF (ASSOCIATED(this%sinnv)) THEN
353  DEALLOCATE(this%sinnv)
354  this%sinnv => null()
355  END IF
356 
357  IF (ASSOCIATED(this%sinnvn)) THEN
358  DEALLOCATE(this%sinnvn)
359  this%sinnvn => null()
360  END IF
361 
362  IF (ASSOCIATED(this%workmj1)) THEN
363  DEALLOCATE(this%workmj1)
364  this%workmj1 => null()
365  END IF
366 
367  IF (ASSOCIATED(this%workmj2)) THEN
368  DEALLOCATE(this%workmj2)
369  this%workmj2 => null()
370  END IF
371 
372  IF (ASSOCIATED(this%workin1)) THEN
373  DEALLOCATE(this%workin1)
374  this%workin1 => null()
375  END IF
376 
377  IF (ASSOCIATED(this%workin2)) THEN
378  DEALLOCATE(this%workin2)
379  this%workin2 => null()
380  END IF
381 
382  END SUBROUTINE
383 
384 !*******************************************************************************
385 ! UTILITY SUBROUTINES
386 !*******************************************************************************
387 !-------------------------------------------------------------------------------
397 !-------------------------------------------------------------------------------
398  SUBROUTINE fourier_tomnsp_3d(this, xuv, xmn, parity)
399 
400  IMPLICIT NONE
401 
402 ! Declare Arguments
403  class(fourier_class), INTENT(inout) :: this
404  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: xuv
405  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: xmn
406  INTEGER, INTENT(in) :: parity
407 
408 ! Local Variables
409  INTEGER :: js
410  REAL (dp) :: skston
411  REAL (dp) :: skstoff
412 
413 ! Start of executable code.
414  CALL second0(skston)
415 
416  DO js = 1, min(SIZE(xmn,3), SIZE(xuv,3))
417  CALL this%tomnsp(xuv(:,:,js), xmn(:,:,js), parity)
418  END DO
419 
420  CALL second0(skstoff)
421  IF (diagonaldone .and. inhessian) THEN
422  time_tomnsp = time_tomnsp + (skstoff - skston)
423  END IF
424 
425  time_tomnsp = time_tomnsp + (skstoff - skston)
426 
427  END SUBROUTINE
428 
429 !-------------------------------------------------------------------------------
442 !-------------------------------------------------------------------------------
443  SUBROUTINE fourier_tomnsp_3d_pest(this, xuv, xmn, lmns, lmnc, &
444  parity, asym)
445 
446  IMPLICIT NONE
447 
448 ! Declare Arguments
449  class(fourier_class), INTENT(inout) :: this
450  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: xuv
451  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: xmn
452  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: lmns
453  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: lmnc
454  INTEGER, INTENT(in) :: parity
455  LOGICAL, INTENT(in) :: asym
456 
457 ! Local Variables
458  INTEGER :: js
459  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: lambda
460  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: dupdu
461  INTEGER :: ns
462  REAL (dp) :: ton
463  REAL (dp) :: toff
464 
465 ! Start of executable code
466  CALL second0(ton)
467 
468  ns = min(SIZE(xmn,3), SIZE(xuv,3), SIZE(lmns, 3))
469 
470 ! Compute lambda.
471  ALLOCATE(lambda(SIZE(xuv,1),SIZE(xuv,2),ns))
472  ALLOCATE(dupdu(SIZE(xuv,1),SIZE(xuv,2),ns))
473 
474  CALL this%toijsp(lmns, lambda, f_none, f_sin)
475  CALL this%toijsp(lmns, dupdu, f_du, f_sin)
476  IF (asym) THEN
477  CALL this%toijsp(lmnc, lambda, f_sum, f_cos)
478  CALL this%toijsp(lmnc, dupdu, ior(f_sum,f_du), f_cos)
479  END IF
480 
481 ! Transform from up to u
482  dupdu = 1 + dupdu
483 
484 ! Interpolate to the full grid.
485  lambda(:,:,1) = lambda(:,:,2)
486  dupdu(:,:,1) = dupdu(:,:,2)
487  DO js = 2, ns - 1
488  lambda(:,:,js) = 0.5*(lambda(:,:,js) + lambda(:,:,js + 1))
489  lambda(:,:,js) = 0.5*(lambda(:,:,js) + lambda(:,:,js + 1))
490  END DO
491  lambda(:,:,ns) = 2*lambda(:,:,ns) - lambda(:,:,ns - 1)
492  dupdu(:,:,ns) = 2*dupdu(:,:,ns) - dupdu(:,:,ns - 1)
493 
494  DO js = 1, ns
495  CALL this%tomnsp(xuv(:,:,js), xmn(:,:,js), lambda(:,:,js), &
496  dupdu(:,:,js), parity, asym)
497  END DO
498 
499  DEALLOCATE(lambda)
500  DEALLOCATE(dupdu)
501 
502  CALL second0(toff)
503  IF (diagonaldone.AND.inhessian) THEN
504  tomnsp_time = tomnsp_time + (toff - ton)
505  END IF
506 
507  time_tomnsp = time_tomnsp + (toff - ton)
508 
509  END SUBROUTINE
510 
511 !-------------------------------------------------------------------------------
521 !-------------------------------------------------------------------------------
522  PURE SUBROUTINE fourier_tomnsp_2d(this, xuv, xmn, parity)
523 
524  IMPLICIT NONE
525 
526 ! Declare Arguments
527  class(fourier_class), INTENT(inout) :: this
528  REAL (dp), DIMENSION(:,:), INTENT(in) :: xuv
529  REAL (dp), DIMENSION(:,:), INTENT(out) :: xmn
530  INTEGER, INTENT(in) :: parity
531 
532 ! Start of executable code.
533 ! First, add over poloidal collocation angles. Calls fourier_tomnsp_2d_u.
534  CALL this%tomnsp(xuv)
535 ! Then add over toroidal collocation angles. Calls fourier_tomnsp_2d_v.
536  CALL this%tomnsp(xmn, parity)
537 
538  END SUBROUTINE
539 
540 !-------------------------------------------------------------------------------
553 !-------------------------------------------------------------------------------
554  PURE SUBROUTINE fourier_tomnsp_2d_pest(this, xuv, xmn, lambda, dupdu, &
555  parity, asym)
556 
557  IMPLICIT NONE
558 
559 ! Declare Arguments
560  class(fourier_class), INTENT(inout) :: this
561  REAL (dp), DIMENSION(:,:), INTENT(in) :: xuv
562  REAL (dp), DIMENSION(:,:), INTENT(out) :: xmn
563  REAL (dp), DIMENSION(:,:), INTENT(in) :: lambda
564  REAL (dp), DIMENSION(:,:), INTENT(in) :: dupdu
565  INTEGER, INTENT(in) :: parity
566  LOGICAL, INTENT(in) :: asym
567 
568 ! Start of executable code.
569 ! First, add over poloidal collocation angles. Calls fourier_tomnsp_2d_u.
570  CALL this%tomnsp(xuv, lambda, dupdu, asym)
571 ! Then add over toroidal collocation angles. Calls fourier_tomnsp_2d_v.
572  CALL this%tomnsp(xmn, parity)
573 
574  END SUBROUTINE
575 
576 !-------------------------------------------------------------------------------
585 !-------------------------------------------------------------------------------
586  PURE SUBROUTINE fourier_tomnsp_2d_u(this, xuv)
587 
588  IMPLICIT NONE
589 
590 ! Declare Arguments
591  class(fourier_class), INTENT(inout) :: this
592  REAL (dp), DIMENSION(:,:), INTENT(in) :: xuv
593 
594 ! Local Variables
595  INTEGER :: j
596  INTEGER :: m
597 
598 ! Start of executable code.
599 
600 ! First, add over poloidal collocation angles. Note use of cosmui/sinmui
601 ! include normalization factors.
602  DO j = 1, SIZE(xuv, 2) ! nzeta
603  DO m = m0, ubound(this%cosmu, 2) ! mpol
604  this%workmj1(m,j) = dot_product(xuv(:,j), this%cosmui(:,m))
605  this%workmj2(m,j) = dot_product(xuv(:,j), this%sinmui(:,m))
606  END DO
607  END DO
608 
609  END SUBROUTINE
610 
611 !-------------------------------------------------------------------------------
623 !-------------------------------------------------------------------------------
624  PURE SUBROUTINE fourier_tomnsp_2d_u_pest(this, xuv, lambda, dupdu, asym)
625 
626  IMPLICIT NONE
627 
628 ! Declare Arguments
629  class(fourier_class), INTENT(inout) :: this
630  REAL (dp), DIMENSION(:,:), INTENT(in) :: xuv
631  REAL (dp), DIMENSION(:,:), INTENT(in) :: lambda
632  REAL (dp), DIMENSION(:,:), INTENT(in) :: dupdu
633  LOGICAL, INTENT(in) :: asym
634 
635 ! Local Variables
636  INTEGER :: i
637  INTEGER :: j
638  INTEGER :: m
639  REAL (dp) :: cosml
640  REAL (dp) :: sinml
641  REAL (dp) :: templ
642  REAL (dp) :: cosmuip
643  REAL (dp) :: sinmuip
644 
645 ! Start of executable code.
646 
647 ! Use recurrence relations below for efficiency.
648 !
649 ! cosmuip => cos(m*u') = cosmiu(u,m)*cos(m*lambda) - sinmui(u,m)*sin(m*lambda)
650 ! sinmuip => sin(m*u') = sinmui(u,m)*cos(m*lambda) + cosmui(u,m)*sin(m*lambda)
651 !
652 ! Let cos(m*lambda) == cosml(u,v,m)
653 ! sin(m*lambda) == sinml(u,v,m)
654 ! use recurrence formulae to compute these efficiently
655 ! cosml(m=0) = 1
656 ! sinml(m=0) = 0
657 ! compute
658 ! cosml(m=1)
659 ! sinml(m=1)
660 ! then
661 ! cosml(m+1) = cosml(m)*cosml(1) - sinml(m)*sinml(1)
662 ! sinml(m+1) = sinml(m)*cosml(1) + cosml(m)*sinml(1)
663 
664 ! First, add over poloidal collocation angles. Note use of cosmui/sinmui
665 ! include normalization factors.
666  DO j = 1, SIZE(xuv, 2) ! nzeta
667  DO m = m0, ubound(this%cosmu, 2) ! mpol
668  this%workmj1(m,j) = 0
669  this%workmj2(m,j) = 0
670  DO i = 1, SIZE(xuv, 1) ! ntheta
671  IF (m .eq. m0) THEN
672  cosml = dupdu(i,j)
673  sinml = 0
674  ELSE IF (m .eq. m1) THEN
675  cosml = cos(lambda(i,j))*dupdu(i,j)
676  sinml = sin(lambda(i,j))*dupdu(i,j)
677  ELSE
678  templ = cosml
679  cosml = cosml*cos(lambda(i,j)) - sinml*sin(lambda(i,j))
680  sinml = sinml*cos(lambda(i,j)) + templ*sin(lambda(i,j))
681  END IF
682 
683  cosmuip = this%cosmui(i,m)*cosml - this%sinmui(i,m)*sinml
684  sinmuip = this%sinmui(i,m)*cosml + this%cosmui(i,m)*sinml
685 
686  this%workmj1(m,j) = this%workmj1(m,j) + xuv(i,j)*cosmuip
687  this%workmj2(m,j) = this%workmj2(m,j) + xuv(i,j)*sinmuip
688  END DO
689  END DO
690  END DO
691 
692  END SUBROUTINE
693 
694 !-------------------------------------------------------------------------------
705 !-------------------------------------------------------------------------------
706  PURE SUBROUTINE fourier_tomnsp_2d_v(this, xmn, parity)
707 
708  IMPLICIT NONE
709 
710 ! Declare Arguments
711  class(fourier_class), INTENT(inout) :: this
712  REAL (dp), DIMENSION(:,:), INTENT(out) :: xmn
713  INTEGER, INTENT(in) :: parity
714 
715 ! Local Variables
716  REAL (dp) :: workcos
717  REAL (dp) :: worksin
718  INTEGER :: j
719  INTEGER :: m
720  INTEGER :: n
721  INTEGER :: moff
722  INTEGER :: noff
723 
724 ! Start of executable code.
725  moff = lbound(xmn,1)
726  noff = lbound(xmn,2) + ubound(this%cosnv,2)
727 
728 ! Then add over toroidal collocation angles.
729  xmn(:,:) = 0
730  DO n = n0, ubound(this%cosnv, 2) ! ntor
731  DO j = 1, SIZE(this%cosnv, 1) ! nzeta
732  DO m = m0, ubound(this%cosmu, 2) ! mpol
733  IF (parity .eq. f_cos) THEN
734  workcos = this%workmj1(m,j)*this%cosnv(j,n)
735  worksin = this%workmj2(m,j)*this%sinnv(j,n)
736 
737  xmn(m + moff,n + noff) = xmn(m + moff,n + noff) &
738  + workcos - worksin
739 
740  IF (n .ne. n0 .and. m .ne. m0) THEN
741  xmn(m + moff,-n + noff) = xmn(m + moff,-n + noff) &
742  + workcos + worksin
743  END IF
744  ELSE
745  workcos = this%workmj2(m,j)*this%cosnv(j,n)
746  worksin = this%workmj1(m,j)*this%sinnv(j,n)
747  xmn(m + moff,n + noff) = xmn(m + moff,n + noff) &
748  + workcos + worksin
749  IF (n .ne. n0) THEN
750  xmn(m + moff,-n + noff) = xmn(m + moff,-n + noff) &
751  + workcos - worksin
752  END IF
753  END IF
754  END DO
755  END DO
756  END DO
757 
758  IF (parity .eq. f_cos) THEN
759  xmn(m0 + moff,:-n1 + noff) = zero ! Redundant: To avoid counting (0,-n) for non-zero n.
760  ELSE
761  xmn(m0 + moff,:n0 + noff) = zero
762  END IF
763 
764  xmn = this%orthonorm*xmn
765 
766  END SUBROUTINE
767 
768 !-------------------------------------------------------------------------------
781 !-------------------------------------------------------------------------------
782  SUBROUTINE fourier_toijsp_3d(this, xmn, xuv, dflag, parity)
783 
784  IMPLICIT NONE
785 
786 ! Declare Arguments
787  class(fourier_class), INTENT(inout) :: this
788  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: xmn
789  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: xuv
790  INTEGER, INTENT(in) :: dflag
791  INTEGER, INTENT(in) :: parity
792 
793 ! Local Variables
794  INTEGER :: js
795  REAL (dp) :: skston
796  REAL (dp) :: skstoff
797 
798 ! Start of executable code.
799  CALL second0(skston)
800 
801  DO js = 1, min(SIZE(xmn,3), SIZE(xuv,3))
802  CALL this%toijsp(xmn(:,:,js), xuv(:,:,js), dflag, parity)
803  END DO
804 
805  CALL second0(skstoff)
806  IF (diagonaldone .AND. inhessian) THEN
807  toijsp_time = toijsp_time + (skstoff - skston)
808  END IF
809 
810  time_toijsp = time_toijsp + (skstoff - skston)
811 
812  END SUBROUTINE
813 
814 !-------------------------------------------------------------------------------
824 !-------------------------------------------------------------------------------
825  PURE SUBROUTINE fourier_toijsp_2d(this, xmn, xuv, dflag, parity)
826 
827  IMPLICIT NONE
828 
829 ! Declare Arguments
830  class(fourier_class), INTENT(inout) :: this
831  REAL (dp), DIMENSION(:,:), INTENT(in) :: xmn
832  REAL (dp), DIMENSION(:,:), INTENT(inout) :: xuv
833  INTEGER, INTENT(in) :: parity
834  INTEGER, INTENT(in) :: dflag
835 
836 ! Local Variables
837  INTEGER :: j
838  INTEGER :: m
839  INTEGER :: n
840  INTEGER :: moff
841  INTEGER :: noff
842  REAL (dp) :: workmn
843  REAL (dp), DIMENSION(:,:), POINTER :: anglem1
844  REAL (dp), DIMENSION(:,:), POINTER :: anglem2
845  INTEGER :: msign
846  REAL (dp), DIMENSION(:,:), POINTER :: anglen1
847  REAL (dp), DIMENSION(:,:), POINTER :: anglen2
848  INTEGER :: nsign1
849  INTEGER :: nsign2
850  INTEGER :: nsign3
851 
852 ! Start of executable code.
853  moff = lbound(xmn,1)
854  noff = lbound(xmn,2) + ubound(this%cosnv,2)
855 
856 ! Poloidal derivative requested use cosmum, sinmum.
857  IF (btest(dflag, b_du)) THEN
858  anglem1 => this%sinmum
859  anglem2 => this%cosmum
860  msign = -1
861  ELSE
862  anglem1 => this%cosmu
863  anglem2 => this%sinmu
864  msign = 1
865  END IF
866 
867  this%workin1 = 0
868  this%workin2 = 0
869 
870 ! Sum over poloidal angles.
871  DO n = lbound(this%workin1, 2), ubound(this%workin1, 2) ! -ntor:ntor
872  DO m = 0, ubound(this%cosmu, 2) ! mpol
873  workmn = this%orthonorm(m,n)*xmn(m + moff,n + noff)
874 
875  this%workin1(:,n) = this%workin1(:,n) + msign*workmn*anglem1(:,m)
876  this%workin2(:,n) = this%workin2(:,n) + workmn*anglem2(:,m)
877  END DO
878  END DO
879 
880  IF (parity .eq. f_cos) THEN
881  nsign3 = 1
882  IF (btest(dflag, b_dv)) THEN
883  anglen1 => this%sinnvn
884  anglen2 => this%cosnvn
885  nsign1 = -1
886  nsign2 = -1
887  ELSE
888  anglen1 => this%cosnv
889  anglen2 => this%sinnv
890  nsign1 = 1
891  nsign2 = -1
892  END IF
893  ELSE
894  nsign3 = -1
895  IF (btest(dflag, b_dv)) THEN
896  anglen1 => this%cosnvn
897  anglen2 => this%sinnvn
898  nsign1 = 1
899  nsign2 = -1
900  ELSE
901  anglen1 => this%sinnv
902  anglen2 => this%cosnv
903  nsign1 = 1
904  nsign2 = 1
905  END IF
906  END IF
907 
908  IF (.not.btest(dflag, b_sum)) THEN
909  xuv(:,:) = 0
910  END IF
911 
912 ! First do the n = 0 mode.
913  IF (.not.btest(dflag, b_dv)) THEN
914  IF (parity .eq. f_cos) THEN
915  DO j = 1, SIZE(xuv, 2) ! nzeta
916  xuv(:,j) = this%workin1(:,n0) + xuv(:,j)
917  END DO
918  ELSE
919  DO j = 1, SIZE(xuv, 2) ! nzeta
920  xuv(:,j) = this%workin2(:,n0) + xuv(:,j)
921  END DO
922  END IF
923  END IF
924 
925 ! Sum over remaining toroidal modes.
926  DO n = n1, ubound(this%cosnv, 2) ! ntor
927  DO j = 1, SIZE(xuv, 2) ! nzeta
928  xuv(:,j) = xuv(:,j) &
929  + nsign1*(this%workin1(:,n) + &
930  nsign3*this%workin1(:,-n))*anglen1(j,n) &
931  + nsign2*(this%workin2(:,n) - &
932  nsign3*this%workin2(:,-n))*anglen2(j,n)
933  END DO
934  END DO
935 
936  END SUBROUTINE
937 
938 !-------------------------------------------------------------------------------
947 !-------------------------------------------------------------------------------
948  SUBROUTINE fourier_get_origin(this, xuv, mode, asym)
949 
950  IMPLICIT NONE
951 
952 ! Declare Arguments
953  class(fourier_class), INTENT(inout) :: this
954  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: xuv
955  INTEGER, INTENT(in) :: mode
956  LOGICAL, INTENT(in) :: asym
957 
958 ! Local Variables
959  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xmn_cos
960  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xmn_sin
961 
962 ! Local Parameters
963  INTEGER, PARAMETER :: js1 = 1
964  integer, PARAMETER :: js2 = 2
965 
966 ! Start of executable code.
967 
968 ! This routine is only called once from init_quantities so just allocate this
969 ! buffer in the call and not the fourier_class.
970  ALLOCATE(xmn_cos(0:ubound(this%cosmu,2), &
971  -ubound(this%cosnv,2):ubound(this%cosnv,2)))
972  CALL this%tomnsp(xuv(:,:,js2), xmn_cos(:,:), f_cos)
973 
974  IF (mode .gt. 0) THEN
975  xmn_cos(0:mode - 1,:) = 0
976  END IF
977  IF (mode .lt. ubound(this%cosmu,2)) THEN
978  xmn_cos(mode + 1,:) = 0
979  END IF
980 
981  IF (asym) THEN
982  ALLOCATE(xmn_sin(0:ubound(this%cosmu,2), &
983  -ubound(this%cosnv,2):ubound(this%cosnv,2)))
984  CALL this%tomnsp(xuv(:,:,js2), xmn_sin(:,:), f_sin)
985 
986  IF (mode .gt. 0) THEN
987  xmn_sin(0:mode - 1,:) = 0
988  END IF
989  IF (mode .lt. ubound(this%cosmu,2)) THEN
990  xmn_sin(mode + 1,:) = 0
991  END IF
992  END IF
993 
994  CALL this%toijsp(xmn_cos(:,:), xuv(:,:,js1), f_cos, f_none)
995  DEALLOCATE(xmn_cos)
996 
997  IF (asym) THEN
998  CALL this%toijsp(xmn_sin(:,:), xuv(:,:,js1), f_sin, f_sum)
999  DEALLOCATE(xmn_sin)
1000  END IF
1001 
1002  END SUBROUTINE
1003 
1004 !-------------------------------------------------------------------------------
1009 !-------------------------------------------------------------------------------
1010  SUBROUTINE fourier_round_trig(cosi, sini)
1011 
1012  IMPLICIT NONE
1013 
1014 ! Declare Arguments
1015  REAL (dp), INTENT(inout) :: cosi
1016  REAL (dp), INTENT(inout) :: sini
1017 
1018 ! Local variables
1019  REAL(dp) :: eps
1020 
1021 ! Start of executable code
1022  eps = 10*epsilon(eps)
1023 
1024  IF (abs(cosi - 1) .le. eps) THEN
1025  cosi = 1
1026  sini = 0
1027  ELSE IF (abs(cosi + 1) .le. eps) THEN
1028  cosi = -1
1029  sini = 0
1030  ELSE IF (abs(sini - 1) .le. eps) THEN
1031  cosi = 0
1032  sini = 1
1033  ELSE IF (abs(sini + 1) .le. eps) THEN
1034  cosi = 0
1035  sini = -1
1036  END IF
1037 
1038  END SUBROUTINE
1039 
1040 !*******************************************************************************
1041 ! UNIT TESTS
1042 !*******************************************************************************
1043 !-------------------------------------------------------------------------------
1048 !-------------------------------------------------------------------------------
1049  FUNCTION test_fourier()
1050  USE timer_mod, ONLY: test_fourier_time
1051 
1052  IMPLICIT NONE
1053 
1054 ! Declare Arguments
1055  LOGICAL :: test_fourier
1056 
1057 ! Local variables
1058  INTEGER :: m
1059  INTEGER :: n
1060  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: testij
1061  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: testmn
1062  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: result
1063  class(fourier_class), POINTER :: four => null()
1064 
1065 ! Local PARAMETERS
1066  INTEGER, PARAMETER :: pol = 10
1067  integer, PARAMETER :: tor = 10
1068  INTEGER, PARAMETER :: theta = 33
1069  INTEGER, PARAMETER :: zeta = 34
1070  INTEGER, PARAMETER :: nfp = 5
1071 
1072 ! Start of executable code.
1073  ALLOCATE(testij(theta,zeta))
1074  ALLOCATE(testmn(0:pol,-tor:tor))
1075  ALLOCATE(result(0:pol,-tor:tor))
1076 
1077  four => fourier_class(pol, tor, theta, zeta, nfp, .false.)
1078  test_fourier = .true.
1079 
1080  DO n = -tor, tor
1081  DO m = 0, pol
1082  IF (m .eq. m0 .and. n .lt. n0) THEN
1083  cycle
1084  END IF
1085 
1086 ! Test cosine parity transform.
1087  testmn = 0
1088  testmn(m,n) = 1
1089  CALL four%toijsp(testmn, testij, f_none, f_cos)
1090  CALL four%tomnsp(testij, result, f_cos)
1091  test_fourier = test_fourier .and. &
1092  check(testmn, result, pol, tor, &
1093  'cosine_parity_test')
1094 
1095 ! Test sine parity transform.
1096  IF (m .eq. m0 .and. n .eq. n0) THEN
1097  cycle
1098  END IF
1099  testmn = 0
1100  testmn(m,n) = 1
1101  CALL four%toijsp(testmn, testij, f_none, f_sin)
1102  CALL four%tomnsp(testij, result, f_sin)
1103  test_fourier = test_fourier .and. &
1104  check(testmn, result, pol, tor, &
1105  'sine_parity_test')
1106 
1107 ! Test u derivatives of cosine parity.
1108  testmn = 0
1109  testmn(m,n) = 1
1110  CALL four%toijsp(testmn, testij, f_du, f_cos)
1111  CALL four%tomnsp(testij, result, f_sin)
1112  testmn(m,n) = -m
1113  test_fourier = test_fourier .and. &
1114  check(testmn, result, pol, tor, &
1115  'du_cosine_parity_test')
1116 
1117 ! Test v derivatives of cosine parity.
1118  testmn = 0
1119  testmn(m,n) = 1
1120  CALL four%toijsp(testmn, testij, f_dv, f_cos)
1121  CALL four%tomnsp(testij, result, f_sin)
1122  testmn(m,n) = -n*nfp
1123  test_fourier = test_fourier .and. &
1124  check(testmn, result, pol, tor, &
1125  'dv_cosine_parity_test')
1126 
1127 ! Test u derivatives of sine parity.
1128  testmn = 0
1129  testmn(m,n) = 1
1130  CALL four%toijsp(testmn, testij, f_du, f_sin)
1131  CALL four%tomnsp(testij, result, f_cos)
1132  testmn(m,n) = m
1133  test_fourier = test_fourier .and. &
1134  check(testmn, result, pol, tor, &
1135  'du_sine_parity_test')
1136 
1137 ! Test v derivatives of sine parity.
1138  testmn = 0
1139  testmn(m,n) = 1
1140  CALL four%toijsp(testmn, testij, f_dv, f_sin)
1141  CALL four%tomnsp(testij, result, f_cos)
1142  testmn(m,n) = n*nfp
1143  test_fourier = test_fourier .and. &
1144  check(testmn, result, pol, tor, &
1145  'dv_sine_parity_test')
1146  END DO
1147  END DO
1148 
1149  DEALLOCATE(four)
1150 
1151  DEALLOCATE(testij)
1152  DEALLOCATE(testmn)
1153  DEALLOCATE(result)
1154 
1155  END FUNCTION
1156 
1157 !*******************************************************************************
1158 ! CHECK FUNCTIONS
1159 !*******************************************************************************
1160 !-------------------------------------------------------------------------------
1168 !-------------------------------------------------------------------------------
1169  FUNCTION check_mn(expected, received, m, n, name)
1170 
1171  IMPLICIT NONE
1172 
1173 ! Declare Arguments
1174  LOGICAL :: check_mn
1175  REAL (rprec), INTENT(in) :: expected
1176  REAL (rprec), INTENT(in) :: received
1177  INTEGER, INTENT(in) :: m
1178  INTEGER, INTENT(in) :: n
1179  CHARACTER (len=*), INTENT(in) :: name
1180 
1181 ! Local Parameters
1182  REAL (rprec), PARAMETER :: eps = 1.e-12_dp
1183 
1184 ! Start of executable code.
1185  check_mn = abs(expected - received) .lt. eps
1186  IF (.not.check_mn) THEN
1187  WRITE (*,1000) trim(name), m, n, expected, received
1188  END IF
1189 
1190 1000 FORMAT('fourier.f90: 'a,' m = ',i3,' n = ',i3,' test failed.',/, &
1191  'Expected ',e12.3,' Recieved ',e12.3)
1192 
1193  END FUNCTION
1194 
1195 !-------------------------------------------------------------------------------
1201 !-------------------------------------------------------------------------------
1202  FUNCTION check_all(expected, received, mpol, ntor, name)
1203 
1204  IMPLICIT NONE
1205 
1206 ! Declare Arguments
1207  LOGICAL :: check_all
1208  REAL (rprec), DIMENSION(:,:), INTENT(in) :: expected
1209  REAL (rprec), DIMENSION(:,:), INTENT(in) :: received
1210  INTEGER, INTENT(in) :: mpol
1211  INTEGER, INTENT(in) :: ntor
1212  CHARACTER (len=*), INTENT(in) :: name
1213 
1214 ! Local Variables
1215  INTEGER :: m
1216  INTEGER :: n
1217  INTEGER :: moff
1218  INTEGER :: noff
1219 
1220 ! Start of executable code.
1221  moff = lbound(expected, 1)
1222  noff = lbound(expected, 2) + ntor
1223 
1224  check_all = .true.
1225 
1226  DO m = 0, mpol
1227  DO n = -ntor, ntor
1228  check_all = check_all .and. check(expected(m + moff, n + noff), &
1229  received(m + moff, n + noff), &
1230  m, n, name)
1231  END DO
1232  END DO
1233 
1234  END FUNCTION
1235 
1236  END MODULE
fourier::fourier_construct
class(fourier_class) function, pointer fourier_construct(mpol, ntor, ntheta, nzeta, nfp, sym)
Construct a fourier_class object..
Definition: fourier.f90:180
fourier::n1
integer, parameter n1
n = 1 mode.
Definition: fourier.f90:67
fourier::b_dv
integer, parameter b_dv
Bit position of the f_dv flag.
Definition: fourier.f90:51
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
fourier::fourier_round_trig
subroutine fourier_round_trig(cosi, sini)
Round trig values to whole number for near numbers.
Definition: fourier.f90:1011
fourier::m0
integer, parameter m0
m = 0 mode.
Definition: fourier.f90:58
fourier::fourier_get_origin
subroutine fourier_get_origin(this, xuv, mode, asym)
Computes the origin value of a half-mesh quantity.
Definition: fourier.f90:949
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
fourier::check_all
logical function, private check_all(expected, received, mpol, ntor, name)
Check all test values.
Definition: fourier.f90:1203
fourier::fourier_tomnsp_2d_v
pure subroutine fourier_tomnsp_2d_v(this, xmn, parity)
Sum over toroidal angle.
Definition: fourier.f90:707
fourier::fourier_toijsp_2d
pure subroutine fourier_toijsp_2d(this, xmn, xuv, dflag, parity)
Fourier transform a 2D quantity to real space.
Definition: fourier.f90:826
fourier::b_sum
integer, parameter b_sum
Bit position of the f_sum flag.
Definition: fourier.f90:53
fourier::b_du
integer, parameter b_du
Bit position of the f_du flag.
Definition: fourier.f90:49
fourier::b_ds
integer, parameter b_ds
Bit position of the f_ds flag.
Definition: fourier.f90:47
fourier::n0
integer, parameter n0
n = 0 mode.
Definition: fourier.f90:65
fourier::f_ds
integer, parameter f_ds
Radial derivative flag.
Definition: fourier.f90:34
fourier::m2
integer, parameter m2
m = 2 mode.
Definition: fourier.f90:62
fourier::f_du
integer, parameter f_du
Poloidal derivative flag.
Definition: fourier.f90:36
fourier::fourier_destruct
subroutine fourier_destruct(this)
Deconstruct a fourier_class object.
Definition: fourier.f90:300
fourier::fourier_tomnsp_3d
subroutine fourier_tomnsp_3d(this, xuv, xmn, parity)
Convert a quantity from real to fourier space in parallel.
Definition: fourier.f90:399
fourier::fourier_tomnsp_3d_pest
subroutine fourier_tomnsp_3d_pest(this, xuv, xmn, lmns, lmnc,
Convert a quantity from real to fourier space in pest coordinates.
Definition: fourier.f90:444
fourier::fourier_tomnsp_2d_u
pure subroutine fourier_tomnsp_2d_u(this, xuv)
Sum over the poloidal part.
Definition: fourier.f90:587
fourier::fourier_tomnsp_2d_pest
pure subroutine fourier_tomnsp_2d_pest(this, xuv, xmn, lambda, dup
Fourier transform a 2D quantity to fourier space.
Definition: fourier.f90:555
fourier::f_dudv
integer, parameter f_dudv
Combined toroidal and poloidal derivatives.
Definition: fourier.f90:40
fourier::f_sum
integer, parameter f_sum
Sum fouier real space transformed quantity. This is used when a non-stellarator symmetric parity is b...
Definition: fourier.f90:43
fourier::m1
integer, parameter m1
m = 1 mode.
Definition: fourier.f90:60
fourier::fourier_toijsp_3d
subroutine fourier_toijsp_3d(this, xmn, xuv, dflag, parity)
Convert a quantity from fourier to real space in parallel.
Definition: fourier.f90:783
fourier::check_mn
logical function, private check_mn(expected, received, m, n, name)
Check test values.
Definition: fourier.f90:1170
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
fourier::f_dv
integer, parameter f_dv
Toroidal derivative flag.
Definition: fourier.f90:38
fourier::fourier_class
Base class containing fourier memory.
Definition: fourier.f90:77
fourier::f_none
integer, parameter f_none
Do not sum fouier real space transformed quantity. This is used when a stellarator symmetric parity i...
Definition: fourier.f90:32
fourier::fourier_tomnsp_2d_u_pest
pure subroutine fourier_tomnsp_2d_u_pest(this, xuv, lambda, dupdu,
Sum over the poloidal part.
Definition: fourier.f90:625
fourier::test_fourier
logical function test_fourier()
Test fourier routines.
Definition: fourier.f90:1050
fourier::check
Interface for checking results of the unit tests.
Definition: fourier.f90:139
fourier::fourier_tomnsp_2d
pure subroutine fourier_tomnsp_2d(this, xuv, xmn, parity)
Fourier transform a 2D quantity to fourier space.
Definition: fourier.f90:523
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17