source: palm/trunk/SOURCE/poisfft_mod.f90 @ 4436

Last change on this file since 4436 was 4429, checked in by raasch, 5 years ago

serial (non-MPI) test case added, several bugfixes for the serial mode

  • Property svn:keywords set to Id
File size: 44.0 KB
Line 
1!> @file poisfft_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poisfft_mod.f90 4429 2020-02-27 15:24:30Z oliver.maas $
27! Statements added to avoid compile errors due to unused dummy arguments in serial mode
28!
29! 4366 2020-01-09 08:12:43Z raasch
30! modification concerning NEC vectorizatio
31!
32! 4360 2020-01-07 11:25:50Z suehring
33! Corrected "Former revisions" section
34!
35! 3690 2019-01-22 22:56:42Z knoop
36! OpenACC port for SPEC
37!
38! Revision 1.1  1997/07/24 11:24:14  raasch
39! Initial revision
40!
41!
42! Description:
43! ------------
44!> Solves the Poisson equation with a 2D spectral method
45!>        d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s
46!>
47!> Input:
48!> real    ar   contains (nnz,nny,nnx) elements of the velocity divergence,
49!>              starting from (1,nys,nxl)
50!>
51!> Output:
52!> real    ar   contains the solution for perturbation pressure p
53!------------------------------------------------------------------------------!
54 MODULE poisfft_mod
55 
56
57    USE fft_xy,                                                                &
58        ONLY:  fft_init, fft_y, fft_y_1d, fft_y_m, fft_x, fft_x_1d, fft_x_m,   &
59               temperton_fft_vec
60
61    USE indices,                                                               &
62        ONLY:  nnx, nny, nx, nxl, nxr, ny, nys, nyn, nz
63
64    USE transpose_indices,                                                     &
65        ONLY:  nxl_y, nxl_z, nxr_y, nxr_z, nys_x, nys_z, nyn_x, nyn_z, nzb_x,  &
66               nzb_y, nzt_x, nzt_y
67
68    USE tridia_solver,                                                         &
69        ONLY:  tridia_1dd, tridia_init, tridia_substi, tridia_substi_overlap
70
71    IMPLICIT NONE
72
73    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
74
75    PRIVATE
76
77    PUBLIC  poisfft, poisfft_init
78
79    INTERFACE poisfft
80       MODULE PROCEDURE poisfft
81    END INTERFACE poisfft
82
83    INTERFACE poisfft_init
84       MODULE PROCEDURE poisfft_init
85    END INTERFACE poisfft_init
86
87
88 CONTAINS
89
90!------------------------------------------------------------------------------!
91! Description:
92! ------------
93!> Setup coefficients for FFT and the tridiagonal solver
94!------------------------------------------------------------------------------!
95    SUBROUTINE poisfft_init
96
97       IMPLICIT NONE
98
99
100       CALL fft_init
101
102       CALL tridia_init
103
104       poisfft_initialized = .TRUE.
105
106    END SUBROUTINE poisfft_init
107
108
109
110!------------------------------------------------------------------------------!
111! Description:
112! ------------
113!> Two-dimensional Fourier Transformation in x- and y-direction.
114!------------------------------------------------------------------------------!
115    SUBROUTINE poisfft( ar )
116
117       USE control_parameters,                                                 &
118           ONLY:  transpose_compute_overlap
119
120       USE cpulog,                                                             &
121           ONLY:  cpu_log, cpu_log_nowait, log_point_s
122
123       USE kinds
124
125       USE pegrid
126
127       IMPLICIT NONE
128
129       INTEGER(iwp) ::  ii           !<
130       INTEGER(iwp) ::  iind         !<
131       INTEGER(iwp) ::  inew         !<
132       INTEGER(iwp) ::  jj           !<
133       INTEGER(iwp) ::  jind         !<
134       INTEGER(iwp) ::  jnew         !<
135       INTEGER(iwp) ::  ki           !<
136       INTEGER(iwp) ::  kk           !<
137       INTEGER(iwp) ::  knew         !<
138       INTEGER(iwp) ::  n            !<
139       INTEGER(iwp) ::  nblk         !<
140       INTEGER(iwp) ::  nnx_y        !<
141       INTEGER(iwp) ::  nny_z        !<
142       INTEGER(iwp) ::  nnz_x        !<
143       INTEGER(iwp) ::  nxl_y_bound  !<
144       INTEGER(iwp) ::  nxr_y_bound  !<
145
146       INTEGER(iwp), DIMENSION(4) ::  isave  !<
147
148       REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar      !<
149       REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv  !<
150
151#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
152#if __acc_fft_device
153       !$ACC DECLARE CREATE(ar_inv)
154#endif
155
156       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  ar1      !<
157       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_in     !<
158       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_inv    !<
159       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_y  !<
160       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_z  !<
161
162
163       CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
164
165       IF ( .NOT. poisfft_initialized )  CALL poisfft_init
166
167#if !__acc_fft_device
168       !$ACC UPDATE HOST(ar)
169#endif
170
171#ifndef _OPENACC
172!
173!--    Two-dimensional Fourier Transformation in x- and y-direction.
174       IF ( pdims(2) == 1  .AND.  pdims(1) > 1 )  THEN
175
176!
177!--       1d-domain-decomposition along x:
178!--       FFT along y and transposition y --> x
179          CALL ffty_tr_yx( ar, ar )
180
181!
182!--       FFT along x, solving the tridiagonal system and backward FFT
183          CALL fftx_tri_fftx( ar )
184
185!
186!--       Transposition x --> y and backward FFT along y
187          CALL tr_xy_ffty( ar, ar )
188
189       ELSEIF ( pdims(1) == 1 .AND. pdims(2) > 1 )  THEN
190
191!
192!--       1d-domain-decomposition along y:
193!--       FFT along x and transposition x --> y
194          CALL fftx_tr_xy( ar, ar )
195
196!
197!--       FFT along y, solving the tridiagonal system and backward FFT
198          CALL ffty_tri_ffty( ar )
199
200!
201!--       Transposition y --> x and backward FFT along x
202          CALL tr_yx_fftx( ar, ar )
203
204       ELSEIF ( .NOT. transpose_compute_overlap )  THEN
205#endif
206
207!
208!--       2d-domain-decomposition or no decomposition (1 PE run)
209!--       Transposition z --> x
210          CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
211          CALL resort_for_zx( ar, ar_inv )
212          CALL transpose_zx( ar_inv, ar )
213          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
214
215          CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
216          IF ( temperton_fft_vec )  THEN
217!
218!--          Vector version outputs a transformed array ar_inv that does not require resorting
219!--          (which is done for ar further below)
220             CALL fft_x( ar, 'forward',  ar_inv=ar_inv)
221          ELSE
222             CALL fft_x( ar, 'forward')
223          ENDIF
224          CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
225
226!
227!--       Transposition x --> y
228          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
229          IF( .NOT. temperton_fft_vec )  CALL resort_for_xy( ar, ar_inv )
230          CALL transpose_xy( ar_inv, ar )
231          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
232
233          CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
234          IF ( temperton_fft_vec )  THEN
235!
236!--          Input array ar_inv from fft_x can be directly used here.
237!--          The output (also in array ar_inv) does not require resorting below.
238             CALL fft_y( ar, 'forward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
239                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
240          ELSE
241             CALL fft_y( ar, 'forward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,      &
242                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
243          ENDIF
244          CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
245
246!
247!--       Transposition y --> z
248          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
249          IF ( .NOT. temperton_fft_vec )  CALL resort_for_yz( ar, ar_inv )
250          CALL transpose_yz( ar_inv, ar )
251          CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
252
253!
254!--       Solve the tridiagonal equation system along z
255          CALL cpu_log( log_point_s(6), 'tridia', 'start' )
256          CALL tridia_substi( ar )
257          CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
258
259!
260!--       Inverse Fourier Transformation
261!--       Transposition z --> y
262          CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
263          CALL transpose_zy( ar, ar_inv )
264!
265!--       The fft_y below (vector branch) can directly process ar_inv (i.e. does not require a
266!--       resorting)
267          IF ( .NOT. temperton_fft_vec )  CALL resort_for_zy( ar_inv, ar )
268          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
269
270          CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
271          IF ( temperton_fft_vec )  THEN
272!
273!--          Output array ar_inv can be used as input to the below fft_x routine without resorting
274             CALL fft_y( ar, 'backward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,&
275                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
276          ELSE
277             CALL fft_y( ar, 'backward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,     &
278                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
279          ENDIF
280
281          CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
282
283!
284!--       Transposition y --> x
285          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
286          CALL transpose_yx( ar, ar_inv )
287          IF ( .NOT. temperton_fft_vec )  CALL resort_for_yx( ar_inv, ar )
288          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
289
290          CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
291          IF ( temperton_fft_vec )  THEN
292             CALL fft_x( ar, 'backward',  ar_inv=ar_inv )
293          ELSE
294             CALL fft_x( ar, 'backward' )
295          ENDIF
296          CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
297
298!
299!--       Transposition x --> z
300          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
301          CALL transpose_xz( ar, ar_inv )
302          CALL resort_for_xz( ar_inv, ar )
303          CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
304
305#ifndef _OPENACC
306       ELSE
307
308!
309!--       2d-domain-decomposition or no decomposition (1 PE run) with
310!--       overlapping transposition / fft
311!--       cputime logging must not use barriers, which would prevent overlapping
312          ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
313                    f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
314!
315!--       Transposition z --> x + subsequent fft along x
316          ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
317          CALL resort_for_zx( ar, f_inv )
318!
319!--       Save original indices and gridpoint counter
320          isave(1) = nz
321          isave(2) = nzb_x
322          isave(3) = nzt_x
323          isave(4) = sendrecvcount_zx
324!
325!--       Set new indices for transformation
326          nblk  = nz / pdims(1)
327          nz    = pdims(1)
328          nnz_x = 1
329          nzb_x = 1 + myidx * nnz_x
330          nzt_x = ( myidx + 1 ) * nnz_x
331          sendrecvcount_zx = nnx * nny * nnz_x
332
333          ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
334          ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
335
336          DO  kk = 1, nblk
337
338             IF ( kk == 1 )  THEN
339                CALL cpu_log( log_point_s(5), 'transpo forward', 'start', cpu_log_nowait )
340             ELSE
341                CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
342             ENDIF
343
344             DO  knew = 1, nz
345                ki = kk + nblk * ( knew - 1 )
346                f_in(:,:,knew) = f_inv(:,:,ki)
347             ENDDO
348
349             CALL transpose_zx( f_in, ar1(:,:,:))
350             CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
351
352             IF ( kk == 1 )  THEN
353                CALL cpu_log( log_point_s(4), 'fft_x', 'start', cpu_log_nowait )
354             ELSE
355                CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
356             ENDIF
357
358             n = isave(2) + kk - 1
359             CALL fft_x( ar1(:,:,:), 'forward',  ar_2d = f_out_z(:,:,n))
360             CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
361
362          ENDDO
363!
364!--       Restore original indices/counters
365          nz               = isave(1)
366          nzb_x            = isave(2)
367          nzt_x            = isave(3)
368          sendrecvcount_zx = isave(4)
369
370          DEALLOCATE( ar1, f_in, f_inv )
371
372!
373!--       Transposition x --> y + subsequent fft along y
374          ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
375          CALL resort_for_xy( f_out_z, f_inv )
376!
377!--       Save original indices and gridpoint counter
378          isave(1) = nx
379          isave(2) = nxl_y
380          isave(3) = nxr_y
381          isave(4) = sendrecvcount_xy
382!
383!--       Set new indices for transformation
384          nblk  = ( ( nx+1 ) / pdims(2) ) - 1
385          nx    = pdims(2)
386          nnx_y = 1
387          nxl_y = myidy * nnx_y
388          nxr_y = ( myidy + 1 ) * nnx_y - 1
389          sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
390
391          ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
392          ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
393
394          DO  ii = 0, nblk
395
396             CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
397
398             DO  inew = 0, nx-1
399                iind = ii + ( nblk + 1 ) * inew
400                f_in(:,:,inew) = f_inv(:,:,iind)
401             ENDDO
402
403             CALL transpose_xy( f_in, ar1(:,:,:) )
404
405             CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
406
407             IF ( ii == 1 )  THEN
408                CALL cpu_log( log_point_s(7), 'fft_y', 'start', cpu_log_nowait )
409             ELSE
410                CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
411             ENDIF
412
413             nxl_y_bound = isave(2)
414             nxr_y_bound = isave(3)
415             n           = isave(2) + ii
416             CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y,               &
417                         nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
418                         nxl_y_l = n, nxr_y_l = n )
419
420             CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
421
422          ENDDO
423!
424!--       Restore original indices/counters
425          nx               = isave(1)
426          nxl_y            = isave(2)
427          nxr_y            = isave(3)
428          sendrecvcount_xy = isave(4)
429
430          DEALLOCATE( ar1, f_in, f_inv )
431
432!
433!--       Transposition y --> z + subsequent tridia + resort for z --> y
434          ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
435          CALL resort_for_yz( f_out_y, f_inv )
436!
437!--       Save original indices and gridpoint counter
438          isave(1) = ny
439          isave(2) = nys_z
440          isave(3) = nyn_z
441          isave(4) = sendrecvcount_yz
442!
443!--       Set new indices for transformation
444          nblk             = ( ( ny+1 ) / pdims(1) ) - 1
445          ny               = pdims(1)
446          nny_z            = 1
447          nys_z            = myidx * nny_z
448          nyn_z            = ( myidx + 1 ) * nny_z - 1
449          sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
450
451          ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) )
452          ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
453
454          DO  jj = 0, nblk
455!
456!--          Forward Fourier Transformation
457!--          Transposition y --> z
458             CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
459
460             DO  jnew = 0, ny-1
461                jind = jj + ( nblk + 1 ) * jnew
462                f_in(:,:,jnew) = f_inv(:,:,jind)
463             ENDDO
464
465             CALL transpose_yz( f_in, ar1(:,:,:) )
466
467             IF ( jj == nblk )  THEN
468                CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
469             ELSE
470                CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
471             ENDIF
472
473!
474!--          Solve the tridiagonal equation system along z
475             CALL cpu_log( log_point_s(6), 'tridia', 'start', cpu_log_nowait )
476
477             n = isave(2) + jj
478             CALL tridia_substi_overlap( ar1(:,:,:), n )
479
480             CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
481
482!
483!--          Inverse Fourier Transformation
484!--          Transposition z --> y
485!--          Only one thread should call MPI routines, therefore forward and
486!--          backward tranpose are in the same section
487             IF ( jj == 0 )  THEN
488                CALL cpu_log( log_point_s(8), 'transpo invers', 'start', cpu_log_nowait )
489             ELSE
490                CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
491             ENDIF
492
493             CALL transpose_zy( ar1(:,:,:), f_in )
494
495             DO  jnew = 0, ny-1
496                jind = jj + ( nblk + 1 ) * jnew
497                f_inv(:,:,jind) = f_in(:,:,jnew)
498             ENDDO
499
500             CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
501
502          ENDDO
503!
504!--       Restore original indices/counters
505          ny               = isave(1)
506          nys_z            = isave(2)
507          nyn_z            = isave(3)
508          sendrecvcount_yz = isave(4)
509
510          CALL resort_for_zy( f_inv, f_out_y )
511
512          DEALLOCATE( ar1, f_in, f_inv )
513
514!
515!--       fft along y backward + subsequent transposition y --> x
516          ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
517!
518!--       Save original indices and gridpoint counter
519          isave(1) = nx
520          isave(2) = nxl_y
521          isave(3) = nxr_y
522          isave(4) = sendrecvcount_xy
523!
524!--       Set new indices for transformation
525          nblk             = (( nx+1 ) / pdims(2) ) - 1
526          nx               = pdims(2)
527          nnx_y            = 1
528          nxl_y            = myidy * nnx_y
529          nxr_y            = ( myidy + 1 ) * nnx_y - 1
530          sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
531
532          ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
533          ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
534
535          DO  ii = 0, nblk
536
537             CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
538
539             n = isave(2) + ii
540             nxl_y_bound = isave(2)
541             nxr_y_bound = isave(3)
542
543             CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y,              &
544                         nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
545                         nxl_y_l = n, nxr_y_l = n )
546
547             IF ( ii == nblk )  THEN
548                CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
549             ELSE
550                CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
551             ENDIF
552
553             CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
554
555             CALL transpose_yx( ar1(:,:,:), f_in )
556
557             DO  inew = 0, nx-1
558                iind = ii + (nblk+1) * inew
559                f_inv(:,:,iind) = f_in(:,:,inew)
560             ENDDO
561
562             CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
563
564          ENDDO
565!
566!--       Restore original indices/counters
567          nx               = isave(1)
568          nxl_y            = isave(2)
569          nxr_y            = isave(3)
570          sendrecvcount_xy = isave(4)
571
572          CALL resort_for_yx( f_inv, f_out_z )
573
574          DEALLOCATE( ar1, f_in, f_inv )
575
576!
577!--       fft along x backward + subsequent final transposition x --> z
578          ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
579!
580!--       Save original indices and gridpoint counter
581          isave(1) = nz
582          isave(2) = nzb_x
583          isave(3) = nzt_x
584          isave(4) = sendrecvcount_zx
585!
586!--       Set new indices for transformation
587          nblk             = nz / pdims(1)
588          nz               = pdims(1)
589          nnz_x            = 1
590          nzb_x            = 1 + myidx * nnz_x
591          nzt_x            = ( myidx + 1 ) * nnz_x
592          sendrecvcount_zx = nnx * nny * nnz_x
593
594          ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
595          ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
596
597          DO  kk = 1, nblk
598
599             CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
600
601             n = isave(2) + kk - 1
602             CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n))
603
604             IF ( kk == nblk )  THEN
605                CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
606             ELSE
607                CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
608             ENDIF
609
610             CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
611
612             CALL transpose_xz( ar1(:,:,:), f_in )
613
614             DO  knew = 1, nz
615                ki = kk + nblk * (knew-1)
616                f_inv(:,:,ki) = f_in(:,:,knew)
617             ENDDO
618
619             IF ( kk == nblk )  THEN
620                CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
621             ELSE
622                CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
623             ENDIF
624
625          ENDDO
626!
627!--       Restore original indices/counters
628          nz               = isave(1)
629          nzb_x            = isave(2)
630          nzt_x            = isave(3)
631          sendrecvcount_zx = isave(4)
632
633          CALL resort_for_xz( f_inv, ar )
634
635          DEALLOCATE( ar1, f_in, f_inv )
636
637       ENDIF
638#endif
639
640#if !__acc_fft_device
641       !$ACC UPDATE DEVICE(ar)
642#endif
643
644       CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
645
646    END SUBROUTINE poisfft
647
648
649!------------------------------------------------------------------------------!
650! Description:
651! ------------
652!> Fourier-transformation along y with subsequent transposition y --> x for
653!> a 1d-decomposition along x.
654!>
655!> @attention The performance of this routine is much faster on the NEC-SX6,
656!>            if the first index of work_ffty_vec is odd. Otherwise
657!>            memory bank conflicts may occur (especially if the index is a
658!>            multiple of 128). That's why work_ffty_vec is dimensioned as
659!>            0:ny+1.
660!>            Of course, this will not work if users are using an odd number
661!>            of gridpoints along y.
662!------------------------------------------------------------------------------!
663    SUBROUTINE ffty_tr_yx( f_in, f_out )
664
665       USE control_parameters,                                                 &
666           ONLY:  loop_optimization
667
668       USE cpulog,                                                             &
669           ONLY:  cpu_log, log_point_s
670
671       USE kinds
672
673       USE pegrid
674
675       IMPLICIT NONE
676
677       INTEGER(iwp)            ::  i            !<
678       INTEGER(iwp)            ::  iend         !<
679       INTEGER(iwp)            ::  iouter       !<
680       INTEGER(iwp)            ::  ir           !<
681       INTEGER(iwp)            ::  j            !<
682       INTEGER(iwp)            ::  k            !<
683
684       INTEGER(iwp), PARAMETER ::  stridex = 4  !<
685
686       REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_in   !<
687       REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_out  !<
688       REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
689
690       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  work_ffty      !<
691       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec  !<
692
693!
694!--    Carry out the FFT along y, where all data are present due to the
695!--    1d-decomposition along x. Resort the data in a way that x becomes
696!--    the first index.
697       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
698
699       IF ( loop_optimization == 'vector' )  THEN
700
701          ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
702!
703!--       Code optimized for vector processors
704          !$OMP PARALLEL PRIVATE ( i, j, k )
705          !$OMP DO
706          DO  i = nxl, nxr
707
708             DO  j = 0, ny
709                DO  k = 1, nz
710                   work_ffty_vec(j,k,i) = f_in(k,j,i)
711                ENDDO
712             ENDDO
713
714             CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
715
716          ENDDO
717
718          !$OMP DO
719          DO  k = 1, nz
720             DO  j = 0, ny
721                DO  i = nxl, nxr
722                   work(i,k,j) = work_ffty_vec(j,k,i)
723                ENDDO
724             ENDDO
725          ENDDO
726          !$OMP END PARALLEL
727
728          DEALLOCATE( work_ffty_vec )
729
730       ELSE
731!
732!--       Cache optimized code.
733          ALLOCATE( work_ffty(0:ny,stridex) )
734!
735!--       The i-(x-)direction is split into a strided outer loop and an inner
736!--       loop for better cache performance
737          !$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
738          !$OMP DO
739          DO  iouter = nxl, nxr, stridex
740
741             iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
742
743             DO  k = 1, nz
744
745                DO  i = iouter, iend
746
747                   ir = i-iouter+1  ! counter within a stride
748                   DO  j = 0, ny
749                      work_ffty(j,ir) = f_in(k,j,i)
750                   ENDDO
751!
752!--                FFT along y
753                   CALL fft_y_1d( work_ffty(:,ir), 'forward' )
754
755                ENDDO
756
757!
758!--             Resort
759                DO  j = 0, ny
760                   DO  i = iouter, iend
761                      work(i,k,j) = work_ffty(j,i-iouter+1)
762                   ENDDO
763                ENDDO
764
765             ENDDO
766
767          ENDDO
768          !$OMP END PARALLEL
769
770          DEALLOCATE( work_ffty )
771
772       ENDIF
773
774       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
775
776!
777!--    Transpose array
778#if defined( __parallel )
779       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
780       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
781       CALL MPI_ALLTOALL( work(nxl,1,0),      sendrecvcount_xy, MPI_REAL, &
782                          f_out(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
783                          comm1dx, ierr )
784       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
785#else
786!
787!--    Next line required to avoid compile error about unused dummy argument in serial mode
788       i = SIZE( f_out )
789#endif
790
791    END SUBROUTINE ffty_tr_yx
792
793
794!------------------------------------------------------------------------------!
795! Description:
796! ------------
797!>  Transposition x --> y with a subsequent backward Fourier transformation for
798!>  a 1d-decomposition along x
799!------------------------------------------------------------------------------!
800    SUBROUTINE tr_xy_ffty( f_in, f_out )
801
802       USE control_parameters,                                                 &
803           ONLY:  loop_optimization
804
805       USE cpulog,                                                             &
806           ONLY:  cpu_log, log_point_s
807
808       USE kinds
809
810       USE pegrid
811
812       IMPLICIT NONE
813
814       INTEGER(iwp)            ::  i            !<
815       INTEGER(iwp)            ::  iend         !<
816       INTEGER(iwp)            ::  iouter       !<
817       INTEGER(iwp)            ::  ir           !<
818       INTEGER(iwp)            ::  j            !<
819       INTEGER(iwp)            ::  k            !<
820
821       INTEGER(iwp), PARAMETER ::  stridex = 4  !<
822
823       REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_in   !<
824       REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_out  !<
825       REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
826
827       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  work_ffty         !<
828       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec     !<
829
830!
831!--    Transpose array
832#if defined( __parallel )
833       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
834       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
835       CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
836                          work(nxl,1,0),     sendrecvcount_xy, MPI_REAL, &
837                          comm1dx, ierr )
838       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
839#else
840!
841!--    Next line required to avoid compile error about unused dummy argument in serial mode
842       i = SIZE( f_in )
843#endif
844
845!
846!--    Resort the data in a way that y becomes the first index and carry out the
847!--    backward fft along y.
848       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
849
850       IF ( loop_optimization == 'vector' )  THEN
851
852          ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
853!
854!--       Code optimized for vector processors
855          !$OMP PARALLEL PRIVATE ( i, j, k )
856          !$OMP DO
857          DO  k = 1, nz
858             DO  j = 0, ny
859                DO  i = nxl, nxr
860                   work_ffty_vec(j,k,i) = work(i,k,j)
861                ENDDO
862             ENDDO
863          ENDDO
864
865          !$OMP DO
866          DO  i = nxl, nxr
867
868             CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
869
870             DO  j = 0, ny
871                DO  k = 1, nz
872                   f_out(k,j,i) = work_ffty_vec(j,k,i)
873                ENDDO
874             ENDDO
875
876          ENDDO
877          !$OMP END PARALLEL
878
879          DEALLOCATE( work_ffty_vec )
880
881       ELSE
882!
883!--       Cache optimized code.
884          ALLOCATE( work_ffty(0:ny,stridex) )
885!
886!--       The i-(x-)direction is split into a strided outer loop and an inner
887!--       loop for better cache performance
888          !$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
889          !$OMP DO
890          DO  iouter = nxl, nxr, stridex
891
892             iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
893
894             DO  k = 1, nz
895!
896!--             Resort
897                DO  j = 0, ny
898                   DO  i = iouter, iend
899                      work_ffty(j,i-iouter+1) = work(i,k,j)
900                   ENDDO
901                ENDDO
902
903                DO  i = iouter, iend
904
905!
906!--                FFT along y
907                   ir = i-iouter+1  ! counter within a stride
908                   CALL fft_y_1d( work_ffty(:,ir), 'backward' )
909
910                   DO  j = 0, ny
911                      f_out(k,j,i) = work_ffty(j,ir)
912                   ENDDO
913                ENDDO
914
915             ENDDO
916
917          ENDDO
918          !$OMP END PARALLEL
919
920          DEALLOCATE( work_ffty )
921
922       ENDIF
923
924       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
925
926    END SUBROUTINE tr_xy_ffty
927
928
929!------------------------------------------------------------------------------!
930! Description:
931! ------------
932!> FFT along x, solution of the tridiagonal system and backward FFT for
933!> a 1d-decomposition along x
934!>
935!> @warning this subroutine may still not work for hybrid parallelization
936!>          with OpenMP (for possible necessary changes see the original
937!>          routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
938!------------------------------------------------------------------------------!
939    SUBROUTINE fftx_tri_fftx( ar )
940
941       USE control_parameters,                                                 &
942           ONLY:  loop_optimization
943
944       USE cpulog,                                                             &
945           ONLY:  cpu_log, log_point_s
946
947       USE grid_variables,                                                     &
948           ONLY:  ddx2, ddy2
949
950       USE kinds
951
952       USE pegrid
953
954       IMPLICIT NONE
955
956       INTEGER(iwp) ::  i                   !<
957       INTEGER(iwp) ::  j                   !<
958       INTEGER(iwp) ::  k                   !<
959       INTEGER(iwp) ::  m                   !<
960       INTEGER(iwp) ::  n                   !<
961!$     INTEGER(iwp) ::  omp_get_thread_num  !<
962       INTEGER(iwp) ::  tn                  !<
963
964       REAL(wp), DIMENSION(0:nx)                          ::  work_fftx  !<
965       REAL(wp), DIMENSION(0:nx,1:nz)                     ::  work_trix  !<
966       REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  ar         !<
967       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
968
969
970       CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
971
972       ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
973
974       tn = 0              ! Default thread number in case of one thread
975!$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
976       DO  j = nys_x, nyn_x
977
978!$        tn = omp_get_thread_num()
979
980          IF ( loop_optimization == 'vector' )  THEN
981!
982!--          Code optimized for vector processors
983             DO  k = 1, nz
984
985                m = 0
986                DO  n = 1, pdims(1)
987                   DO  i = 1, nnx
988                      work_trix(m,k) = ar(i,k,j,n)
989                      m = m + 1
990                   ENDDO
991                ENDDO
992
993             ENDDO
994
995             CALL fft_x_m( work_trix, 'forward' )
996
997          ELSE
998!
999!--          Cache optimized code
1000             DO  k = 1, nz
1001
1002                m = 0
1003                DO  n = 1, pdims(1)
1004                   DO  i = 1, nnx
1005                      work_fftx(m) = ar(i,k,j,n)
1006                      m = m + 1
1007                   ENDDO
1008                ENDDO
1009
1010                CALL fft_x_1d( work_fftx, 'forward' )
1011
1012                DO  i = 0, nx
1013                   work_trix(i,k) = work_fftx(i)
1014                ENDDO
1015
1016             ENDDO
1017
1018          ENDIF
1019
1020!
1021!--       Solve the linear equation system
1022          CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
1023
1024          IF ( loop_optimization == 'vector' )  THEN
1025!
1026!--          Code optimized for vector processors
1027             CALL fft_x_m( work_trix, 'backward' )
1028
1029             DO  k = 1, nz
1030
1031                m = 0
1032                DO  n = 1, pdims(1)
1033                   DO  i = 1, nnx
1034                      ar(i,k,j,n) = work_trix(m,k)
1035                      m = m + 1
1036                   ENDDO
1037                ENDDO
1038
1039             ENDDO
1040
1041          ELSE
1042!
1043!--          Cache optimized code
1044             DO  k = 1, nz
1045
1046                DO  i = 0, nx
1047                   work_fftx(i) = work_trix(i,k)
1048                ENDDO
1049
1050                CALL fft_x_1d( work_fftx, 'backward' )
1051
1052                m = 0
1053                DO  n = 1, pdims(1)
1054                   DO  i = 1, nnx
1055                      ar(i,k,j,n) = work_fftx(m)
1056                      m = m + 1
1057                   ENDDO
1058                ENDDO
1059
1060             ENDDO
1061
1062          ENDIF
1063
1064       ENDDO
1065
1066       DEALLOCATE( tri )
1067
1068       CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
1069
1070    END SUBROUTINE fftx_tri_fftx
1071
1072
1073!------------------------------------------------------------------------------!
1074! Description:
1075! ------------
1076!> Fourier-transformation along x with subsequent transposition x --> y for
1077!> a 1d-decomposition along y.
1078!>
1079!> @attention NEC-branch of this routine may significantly profit from
1080!>            further optimizations. So far, performance is much worse than
1081!>            for routine ffty_tr_yx (more than three times slower).
1082!------------------------------------------------------------------------------!
1083    SUBROUTINE fftx_tr_xy( f_in, f_out )
1084
1085
1086       USE control_parameters,                                                 &
1087           ONLY:  loop_optimization
1088
1089       USE cpulog,                                                             &
1090           ONLY:  cpu_log, log_point_s
1091
1092       USE kinds
1093
1094       USE pegrid
1095
1096       IMPLICIT NONE
1097
1098       INTEGER(iwp) ::  i  !<
1099       INTEGER(iwp) ::  j  !<
1100       INTEGER(iwp) ::  k  !<
1101
1102       REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1103       REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_in       !<
1104       REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_out      !<
1105       REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1106
1107!
1108!--    Carry out the FFT along x, where all data are present due to the
1109!--    1d-decomposition along y. Resort the data in a way that y becomes
1110!--    the first index.
1111       CALL cpu_log( log_point_s(4), 'fft_x_1d', 'start' )
1112
1113       IF ( loop_optimization == 'vector' )  THEN
1114!
1115!--       Code for vector processors
1116!$OMP     PARALLEL PRIVATE ( i, j, k )
1117!$OMP     DO
1118          DO  i = 0, nx
1119
1120             DO  j = nys, nyn
1121                DO  k = 1, nz
1122                   work_fftx(i,k,j) = f_in(k,j,i)
1123                ENDDO
1124             ENDDO
1125
1126          ENDDO
1127
1128!$OMP     DO
1129          DO  j = nys, nyn
1130
1131             CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1132
1133             DO  k = 1, nz
1134                DO  i = 0, nx
1135                   work(j,k,i) = work_fftx(i,k,j)
1136                ENDDO
1137             ENDDO
1138
1139          ENDDO
1140!$OMP     END PARALLEL
1141
1142       ELSE
1143
1144!
1145!--       Cache optimized code (there might be still a potential for better
1146!--       optimization).
1147!$OMP     PARALLEL PRIVATE (i,j,k)
1148!$OMP     DO
1149          DO  i = 0, nx
1150
1151             DO  j = nys, nyn
1152                DO  k = 1, nz
1153                   work_fftx(i,k,j) = f_in(k,j,i)
1154                ENDDO
1155             ENDDO
1156
1157          ENDDO
1158
1159!$OMP     DO
1160          DO  j = nys, nyn
1161             DO  k = 1, nz
1162
1163                CALL fft_x_1d( work_fftx(0:nx,k,j), 'forward' )
1164
1165                DO  i = 0, nx
1166                   work(j,k,i) = work_fftx(i,k,j)
1167                ENDDO
1168             ENDDO
1169
1170          ENDDO
1171!$OMP     END PARALLEL
1172
1173       ENDIF
1174       CALL cpu_log( log_point_s(4), 'fft_x_1d', 'pause' )
1175
1176!
1177!--    Transpose array
1178#if defined( __parallel )
1179       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1180       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1181       CALL MPI_ALLTOALL( work(nys,1,0),      sendrecvcount_xy, MPI_REAL, &
1182                          f_out(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1183                          comm1dy, ierr )
1184       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1185#else
1186!
1187!--    Next line required to avoid compile error about unused dummy argument in serial mode
1188       i = SIZE( f_out )
1189#endif
1190
1191    END SUBROUTINE fftx_tr_xy
1192
1193
1194!------------------------------------------------------------------------------!
1195! Description:
1196! ------------
1197!> Transposition y --> x with a subsequent backward Fourier transformation for
1198!> a 1d-decomposition along x.
1199!------------------------------------------------------------------------------!
1200    SUBROUTINE tr_yx_fftx( f_in, f_out )
1201
1202
1203       USE control_parameters,                                                 &
1204           ONLY:  loop_optimization
1205
1206       USE cpulog,                                                             &
1207           ONLY:  cpu_log, log_point_s
1208
1209       USE kinds
1210
1211       USE pegrid
1212
1213       IMPLICIT NONE
1214
1215       INTEGER(iwp) ::  i  !<
1216       INTEGER(iwp) ::  j  !<
1217       INTEGER(iwp) ::  k  !<
1218
1219       REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1220       REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_in       !<
1221       REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_out      !<
1222       REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1223
1224
1225!
1226!--    Transpose array
1227#if defined( __parallel )
1228       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1229       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1230       CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1231                          work(nys,1,0),     sendrecvcount_xy, MPI_REAL, &
1232                          comm1dy, ierr )
1233       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1234#else
1235!
1236!--    Next line required to avoid compile error about unused dummy argument in serial mode
1237       i = SIZE( f_in )
1238#endif
1239
1240!
1241!--    Carry out the FFT along x, where all data are present due to the
1242!--    1d-decomposition along y. Resort the data in a way that y becomes
1243!--    the first index.
1244       CALL cpu_log( log_point_s(4), 'fft_x_1d', 'continue' )
1245
1246       IF ( loop_optimization == 'vector' )  THEN
1247!
1248!--       Code optimized for vector processors
1249!$OMP     PARALLEL PRIVATE ( i, j, k )
1250!$OMP     DO
1251          DO  j = nys, nyn
1252
1253             DO  k = 1, nz
1254                DO  i = 0, nx
1255                   work_fftx(i,k,j) = work(j,k,i)
1256                ENDDO
1257             ENDDO
1258
1259             CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1260
1261          ENDDO
1262
1263!$OMP     DO
1264          DO  i = 0, nx
1265             DO  j = nys, nyn
1266                DO  k = 1, nz
1267                   f_out(k,j,i) = work_fftx(i,k,j)
1268                ENDDO
1269             ENDDO
1270          ENDDO
1271!$OMP     END PARALLEL
1272
1273       ELSE
1274
1275!
1276!--       Cache optimized code (there might be still a potential for better
1277!--       optimization).
1278!$OMP     PARALLEL PRIVATE (i,j,k)
1279!$OMP     DO
1280          DO  j = nys, nyn
1281             DO  k = 1, nz
1282
1283                DO  i = 0, nx
1284                   work_fftx(i,k,j) = work(j,k,i)
1285                ENDDO
1286
1287                CALL fft_x_1d( work_fftx(0:nx,k,j), 'backward' )
1288
1289             ENDDO
1290          ENDDO
1291
1292!$OMP     DO
1293          DO  i = 0, nx
1294             DO  j = nys, nyn
1295                DO  k = 1, nz
1296                   f_out(k,j,i) = work_fftx(i,k,j)
1297                ENDDO
1298             ENDDO
1299          ENDDO
1300!$OMP     END PARALLEL
1301
1302       ENDIF
1303       CALL cpu_log( log_point_s(4), 'fft_x_1d', 'stop' )
1304
1305    END SUBROUTINE tr_yx_fftx
1306
1307
1308!------------------------------------------------------------------------------!
1309! Description:
1310! ------------
1311!> FFT along y, solution of the tridiagonal system and backward FFT for
1312!> a 1d-decomposition along y.
1313!>
1314!> @warning this subroutine may still not work for hybrid parallelization
1315!>          with OpenMP (for possible necessary changes see the original
1316!>          routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
1317!------------------------------------------------------------------------------!
1318    SUBROUTINE ffty_tri_ffty( ar )
1319
1320
1321       USE control_parameters,                                                 &
1322           ONLY:  loop_optimization
1323
1324       USE cpulog,                                                             &
1325           ONLY:  cpu_log, log_point_s
1326
1327       USE grid_variables,                                                     &
1328           ONLY:  ddx2, ddy2
1329
1330       USE kinds
1331
1332       USE pegrid
1333
1334       IMPLICIT NONE
1335
1336       INTEGER(iwp) ::  i                   !<
1337       INTEGER(iwp) ::  j                   !<
1338       INTEGER(iwp) ::  k                   !<
1339       INTEGER(iwp) ::  m                   !<
1340       INTEGER(iwp) ::  n                   !<
1341!$     INTEGER(iwp) ::  omp_get_thread_num  !<
1342       INTEGER(iwp) ::  tn                  !<
1343
1344       REAL(wp), DIMENSION(0:ny)                          ::  work_ffty  !<
1345       REAL(wp), DIMENSION(0:ny,1:nz)                     ::  work_triy  !<
1346       REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  ar         !<
1347       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
1348
1349
1350       CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'start' )
1351
1352       ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1353
1354       tn = 0           ! Default thread number in case of one thread
1355!$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1356       DO  i = nxl_y, nxr_y
1357
1358!$        tn = omp_get_thread_num()
1359
1360          IF ( loop_optimization == 'vector' )  THEN
1361!
1362!--          Code optimized for vector processors
1363             DO  k = 1, nz
1364
1365                m = 0
1366                DO  n = 1, pdims(2)
1367                   DO  j = 1, nny
1368                      work_triy(m,k) = ar(j,k,i,n)
1369                      m = m + 1
1370                   ENDDO
1371                ENDDO
1372
1373             ENDDO
1374
1375             CALL fft_y_m( work_triy, ny, 'forward' )
1376
1377          ELSE
1378!
1379!--          Cache optimized code
1380             DO  k = 1, nz
1381
1382                m = 0
1383                DO  n = 1, pdims(2)
1384                   DO  j = 1, nny
1385                      work_ffty(m) = ar(j,k,i,n)
1386                      m = m + 1
1387                   ENDDO
1388                ENDDO
1389
1390                CALL fft_y_1d( work_ffty, 'forward' )
1391
1392                DO  j = 0, ny
1393                   work_triy(j,k) = work_ffty(j)
1394                ENDDO
1395
1396             ENDDO
1397
1398          ENDIF
1399
1400!
1401!--       Solve the linear equation system
1402          CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1403
1404          IF ( loop_optimization == 'vector' )  THEN
1405!
1406!--          Code optimized for vector processors
1407             CALL fft_y_m( work_triy, ny, 'backward' )
1408
1409             DO  k = 1, nz
1410
1411                m = 0
1412                DO  n = 1, pdims(2)
1413                   DO  j = 1, nny
1414                      ar(j,k,i,n) = work_triy(m,k)
1415                      m = m + 1
1416                   ENDDO
1417                ENDDO
1418
1419             ENDDO
1420
1421          ELSE
1422!
1423!--          Cache optimized code
1424             DO  k = 1, nz
1425
1426                DO  j = 0, ny
1427                   work_ffty(j) = work_triy(j,k)
1428                ENDDO
1429
1430                CALL fft_y_1d( work_ffty, 'backward' )
1431
1432                m = 0
1433                DO  n = 1, pdims(2)
1434                   DO  j = 1, nny
1435                      ar(j,k,i,n) = work_ffty(m)
1436                      m = m + 1
1437                   ENDDO
1438                ENDDO
1439
1440             ENDDO
1441
1442          ENDIF
1443
1444       ENDDO
1445
1446       DEALLOCATE( tri )
1447
1448       CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'stop' )
1449
1450    END SUBROUTINE ffty_tri_ffty
1451
1452 END MODULE poisfft_mod
Note: See TracBrowser for help on using the repository browser.