source: palm/trunk/SOURCE/poisfft_hybrid.f90 @ 1036

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

  • Property svn:keywords set to Id
File size: 32.9 KB
Line 
1 MODULE poisfft_hybrid_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poisfft_hybrid.f90 1036 2012-10-22 13:43:42Z raasch $
27!
28! 1013 2012-09-21 07:03:55Z raasch
29! FLOAT type conversion replaced by REAL
30!
31! 809 2012-01-30 13:32:58Z maronga
32! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
33!
34! 807 2012-01-25 11:53:51Z maronga
35! New cpp directive "__check" implemented which is used by check_namelist_files
36! (most of the code is unneeded by check_namelist_files).
37!
38! 667 2010-12-23 12:06:00Z suehring/gryschka
39! ddzu replaced by ddzu_pres due to changes in zu(0)
40!
41! 415 2009-12-15 10:26:23Z raasch
42! Dimension of array stat in cascade change to prevent type problems with___
43! mpi2 libraries
44!
45! 274 2009-03-26 15:11:21Z heinze
46! Output of messages replaced by message handling routine.
47!
48! Feb. 2007
49! RCS Log replace by Id keyword, revision history cleaned up
50!
51! Revision 1.11  2004/04/30 12:43:14  raasch
52! Renaming of fft routines, additional argument in calls of fft_y_m
53!
54! Revision 1.2  2002/12/19 16:08:31  raasch
55! Preprocessor directive KKMP introduced (OMP does NOT work),
56! array tri will be a shared array in OpenMP loop, to get better cache
57! utilization, the i index  (x-direction) will be executed in stride
58! "istride" as outer loop and in a shorter inner loop,
59! overlapping of computation and communication realized by new routine
60! poisfft_hybrid_nodes, name of old routine poisfft_hybrid changed to
61! poisfft_hybrid_omp, STOP statement replaced by call of subroutine local_stop
62!
63!
64! Description:
65! ------------
66! Solution of the Poisson equation with a 2D spectral method. 
67! Hybrid version for parallel computers using a 1D domain decomposition,
68! realized with MPI, along x and parallelization with OPEN-MP along y
69! (routine poisfft_hybrid_omp). In a second version (poisfft_hybrid_nodes),
70! optimization is realized by overlapping of computation and communication
71! and by simultaneously executing as many communication calls as switches
72! per logical partition (LPAR) are available. This version comes into
73! effect if more than one node is used and if the environment variable
74! tasks_per_node is set in a way that it can be devided by switch_per_lpar
75! without any rest.
76!
77! WARNING: In case of OpenMP, there are problems with allocating large
78! arrays in parallel regions.
79!
80! Copyright   Klaus Ketelsen / Siegfried Raasch   May 2002
81!------------------------------------------------------------------------------!
82
83    USE fft_xy
84    USE indices
85    USE pegrid
86
87    IMPLICIT NONE
88
89    INTEGER, PARAMETER ::  switch_per_lpar = 2
90
91    INTEGER, SAVE ::  nxl_a, nxr_a, &      ! total   x dimension
92                      nxl_p, nxr_p, &      ! partial x dimension
93                      nys_a, nyn_a, &      ! total   y dimension
94                      nys_p, nyn_p, &      ! partial y dimension
95
96                      npe_s,        &      ! total number of PEs for solver
97                      nwords,       &      ! number of points to be exchanged
98                                           ! with MPI_ALLTOALL
99                      n_omp_threads        ! number of OpenMP threads
100
101!
102!-- Variables for multi node version (cluster version) using routine
103!-- poisfft_hybrid_nodes
104    INTEGER, SAVE :: comm_nodes,          & ! communicater nodes
105                     comm_node_all,       & ! communicater all PEs node version
106                     comm_tasks,          & ! communicater tasks
107                     me, me_node, me_task,& ! identity of this PE
108                     nodes,               & ! number of nodes
109                     tasks_per_logical_node = -1    ! default no cluster
110
111
112    PRIVATE
113
114
115#if ! defined ( __check )
116    PUBLIC poisfft_hybrid, poisfft_hybrid_ini
117
118
119!
120!-- Public interfaces
121    INTERFACE poisfft_hybrid_ini
122       MODULE PROCEDURE poisfft_hybrid_ini
123    END INTERFACE poisfft_hybrid_ini
124
125    INTERFACE poisfft_hybrid
126       MODULE PROCEDURE poisfft_hybrid
127    END INTERFACE poisfft_hybrid
128
129!
130!-- Private interfaces
131    INTERFACE poisfft_hybrid_omp
132       MODULE PROCEDURE poisfft_hybrid_omp
133    END INTERFACE poisfft_hybrid_omp
134
135    INTERFACE poisfft_hybrid_omp_vec
136       MODULE PROCEDURE poisfft_hybrid_omp_vec
137    END INTERFACE poisfft_hybrid_omp_vec
138
139    INTERFACE poisfft_hybrid_nodes
140       MODULE PROCEDURE poisfft_hybrid_nodes
141    END INTERFACE poisfft_hybrid_nodes
142
143    INTERFACE tridia_hybrid
144       MODULE PROCEDURE tridia_hybrid
145    END INTERFACE tridia_hybrid
146
147    INTERFACE cascade
148       MODULE PROCEDURE cascade
149    END INTERFACE cascade
150#else
151    PUBLIC poisfft_hybrid_ini
152
153!
154!-- Public interfaces
155    INTERFACE poisfft_hybrid_ini
156       MODULE PROCEDURE poisfft_hybrid_ini
157    END INTERFACE poisfft_hybrid_ini
158#endif
159
160 CONTAINS
161 
162
163    SUBROUTINE poisfft_hybrid_ini
164
165       USE control_parameters
166       USE pegrid
167
168       IMPLICIT NONE
169
170       CHARACTER(LEN=8)      ::  cdummy
171       INTEGER               ::  idummy, istat
172       INTEGER, DIMENSION(2) ::  coords, dims
173
174       LOGICAL, DIMENSION(2) ::  period = .false., re_dims
175
176
177!
178!--    Set the internal index values for the hybrid solver
179#if defined( __parallel )
180       npe_s = pdims(1)
181#else
182       npe_s = 1
183#endif
184       nxl_a = 0
185       nxr_a = nx
186       nxl_p = 0
187       nxr_p = ( ( nx+1 ) / npe_s ) - 1
188       nys_a = nys 
189       nyn_a = nyn
190       nys_p = 0
191       nyn_p = ( ( ny+1 ) / npe_s ) - 1
192
193       nwords = ( nxr_p-nxl_p+1 ) * nz * ( nyn_p-nys_p+1 )
194
195#if defined( __KKMP ) && ! defined ( __check )
196       CALL LOCAL_GETENV( 'OMP_NUM_THREADS', 15, cdummy, idummy )
197       READ ( cdummy, '(I8)' )  n_omp_threads
198       IF ( n_omp_threads > 1 )  THEN
199          WRITE( message_string, * ) 'Number of OpenMP threads = ', &
200                                     n_omp_threads
201          CALL message( 'poisfft_hybrid_ini', 'PA0280', 0, 0, 0, 6, 0 ) 
202       ENDIF
203#else
204       n_omp_threads = 1
205#endif
206!
207!--    Initialize the one-dimensional FFT routines
208       CALL fft_init
209
210!
211!--    Setup for multi node version (poisfft_hybrid_nodes)
212       IF ( n_omp_threads == 1  .AND.  &
213            ( host(1:4) == 'ibmh' .OR. host(1:4) == 'ibmb' ) )  THEN
214
215          IF ( tasks_per_node /= -9999 )  THEN
216!
217!--          Multi node version requires that the available number of
218!--          switches per logical partition must be an integral divisor
219!--          of the chosen number of tasks per node
220             IF ( MOD( tasks_per_node, switch_per_lpar ) == 0 )  THEN
221!
222!--             Set the switch which decides about usage of the multi node
223!--             version
224                IF ( tasks_per_node / switch_per_lpar > 1  .AND. &
225                     numprocs > tasks_per_node )  THEN
226                   tasks_per_logical_node = tasks_per_node / switch_per_lpar
227                ENDIF
228
229                IF ( tasks_per_logical_node > -1 )  THEN
230
231                   WRITE( message_string, * ) 'running optimized ',         &
232                                              'multinode version',          &
233                                              '&switch_per_lpar        = ', &
234                                              switch_per_lpar,              &
235                                              '&tasks_per_lpar         = ', &
236                                              tasks_per_node,               &
237                                              'tasks_per_logical_node = ',  &
238                                              tasks_per_logical_node
239                   CALL message( 'poisfft_hybrid_ini', 'PA0281', 0, 0, 0, 6, 0 )
240
241                ENDIF
242
243             ENDIF
244          ENDIF
245       ENDIF
246
247!
248!--    Determine sub-topologies for multi node version
249       IF ( tasks_per_logical_node >= 2 )  THEN
250
251#if defined( __parallel ) && ! defined ( __check )
252          nodes   = ( numprocs + tasks_per_logical_node - 1 ) / &
253                    tasks_per_logical_node
254          dims(1) = nodes
255          dims(2) = tasks_per_logical_node
256
257          CALL MPI_CART_CREATE( comm2d, 2, dims, period, .FALSE., &
258                                comm_node_all, istat )
259          CALL MPI_COMM_RANK( comm_node_all, me, istat )
260
261          re_dims(1) = .TRUE.
262          re_dims(2) = .FALSE.
263          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_nodes, istat )
264          CALL MPI_COMM_RANK( comm_nodes, me_node, istat )
265
266          re_dims(1) = .FALSE.
267          re_dims(2) = .TRUE.
268          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_tasks, istat )
269          CALL MPI_COMM_RANK( comm_tasks, me_task, istat )
270
271!          write(0,*) 'who am i',myid,me,me_node,me_task,nodes,&
272!                     tasks_per_logical_node
273#elif ! defined( __parallel )
274          message_string = 'parallel environment (MPI) required'
275          CALL message( 'poisfft_hybrid_ini', 'PA0282', 1, 2, 0, 6, 0 )
276#endif
277       ENDIF
278
279    END SUBROUTINE poisfft_hybrid_ini
280
281#if ! defined ( __check )
282    SUBROUTINE poisfft_hybrid( ar )
283
284       USE control_parameters
285       USE interfaces
286
287       IMPLICIT NONE
288
289       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
290
291       IF ( host(1:3) == 'nec' )  THEN
292          CALL poisfft_hybrid_omp_vec( ar )
293       ELSE
294          IF ( tasks_per_logical_node == -1 )  THEN
295             CALL poisfft_hybrid_omp( ar )
296          ELSE
297             CALL poisfft_hybrid_nodes( ar )
298          ENDIF
299       ENDIF
300
301    END SUBROUTINE poisfft_hybrid
302
303
304    SUBROUTINE poisfft_hybrid_omp ( ar )
305
306       USE cpulog
307       USE interfaces
308
309       IMPLICIT NONE
310
311       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
312       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
313
314       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
315
316       REAL, DIMENSION(0:nx)                 ::  fftx_ar
317       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
318 
319       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
320
321       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
322#if defined( __KKMP )
323       INTEGER ::  omp_get_thread_num
324       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
325       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
326#else
327       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
328#endif
329
330
331       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'start' )
332
333       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
334
335!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
336!$OMP  DO
337!
338!--    Store grid points to be transformed on a 1d-array, do the fft
339!--    and sample the results on a 4d-array
340       DO  iouter = nxl_p, nxr_p, istride     ! stride loop, better cache
341          iei = MIN( iouter+istride-1, nxr_p )
342          DO  k = 1, nz
343
344             DO  i = iouter, iei
345                ii = nxl + i
346                ir = i - iouter + 1
347
348                DO  j = nys_a, nyn_a
349                   ffty_ar(j,ir) = ar(k,j,ii)
350                ENDDO
351   
352                CALL fft_y( ffty_ar(:,ir), 'forward' )
353             ENDDO
354
355             m = nys_a
356             DO  n = 1, npe_s
357                DO  j = nys_p, nyn_p
358                   DO  i = iouter, iei
359                      ir = i - iouter + 1
360                      work1(i,k,j,n) = ffty_ar(m,ir)
361                   ENDDO
362                   m = m+1
363                ENDDO
364             ENDDO
365   
366          ENDDO
367       ENDDO
368!$OMP  END PARALLEL
369
370       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
371
372#if defined( __parallel )
373       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
374
375       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
376                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
377                          comm2d, istat )
378
379       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
380#else
381       work2 = work1
382#endif
383
384       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
385
386#if defined( __KKMP )
387!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,fftx_ar,tri_ar,jthread)
388!$OMP  DO
389       DO  j = nys_p, nyn_p
390          jthread = omp_get_thread_num() + 1
391#else
392       DO  j = nys_p, nyn_p
393          jthread = 1
394#endif
395          DO  k = 1, nz
396
397             m = nxl_a
398             DO  n = 1, npe_s
399                DO  i = nxl_p, nxr_p
400                   fftx_ar(m) = work2(i,k,j,n)
401                   m = m+1
402                ENDDO
403             ENDDO
404
405             CALL fft_x( fftx_ar, 'forward' )
406
407             DO  i = nxl_a, nxr_a
408                tri_ar(i,k) = fftx_ar(i)
409             ENDDO
410
411          ENDDO
412 
413          jj = myid * (nyn_p-nys_p+1) + j
414          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
415
416          DO  k = 1, nz
417             DO  i = nxl_a, nxr_a
418                fftx_ar(i) = tri_ar (i,k)
419             ENDDO
420
421             CALL fft_x( fftx_ar, 'backward' )
422
423             m = nxl_a
424             DO  n = 1, npe_s
425                DO  i = nxl_p, nxr_p
426                   work2(i,k,j,n) = fftx_ar(m)
427                   m = m+1
428                ENDDO
429             ENDDO
430
431          ENDDO
432       ENDDO
433#if defined( __KKMP )
434!$OMP  END PARALLEL
435#endif
436
437       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
438
439#if defined( __parallel )
440       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
441       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
442
443       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
444                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
445                          comm2d, istat )
446
447       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
448#else
449       work1 = work2
450#endif
451
452       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
453
454!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
455!$OMP  DO
456       DO  iouter = nxl_p, nxr_p, istride
457          iei = MIN( iouter+istride-1, nxr_p )
458          DO  k = 1, nz
459
460             m = nys_a
461             DO  n = 1, npe_s
462                DO  j = nys_p, nyn_p
463                   DO  i = iouter, iei
464                      ir = i - iouter + 1
465                      ffty_ar(m,ir) = work1 (i,k,j,n)
466                   ENDDO
467                   m = m+1
468                ENDDO
469             ENDDO
470
471             DO  i = iouter, iei
472                ii = nxl + i
473                ir = i - iouter + 1
474                CALL fft_y( ffty_ar(:,ir), 'backward' )
475
476                DO  j = nys_a, nyn_a
477                   ar(k,j,ii) = ffty_ar(j,ir)
478                ENDDO
479             ENDDO
480
481          ENDDO
482       ENDDO
483!$OMP  END PARALLEL
484
485       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
486
487       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'stop' )
488
489#if defined( __KKMP )
490       DEALLOCATE( tri )
491#endif
492
493    END SUBROUTINE poisfft_hybrid_omp
494
495
496    SUBROUTINE poisfft_hybrid_omp_vec ( ar )
497
498       USE cpulog
499       USE interfaces
500
501       IMPLICIT NONE
502
503       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
504       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
505
506       REAL, DIMENSION(0:nx,nz)               ::  tri_ar
507
508       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr)  ::  ar
509
510       REAL, DIMENSION(0:ny+3,nz,nxl_p:nxr_p) ::  ffty_ar3
511       REAL, DIMENSION(0:nx+3,nz,nys_p:nyn_p) ::  fftx_ar3
512 
513       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
514#if defined( __KKMP )
515       INTEGER ::  omp_get_thread_num
516       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
517       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
518#else
519       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
520#endif
521
522
523       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'start' )
524
525       CALL cpu_log( log_point_s(7), 'fft_y_m', 'start' )
526
527!$OMP  PARALLEL PRIVATE (i,j,k,m,n)
528!$OMP  DO
529!
530!--    Store grid points to be transformed on a 1d-array, do the fft
531!--    and sample the results on a 4d-array
532       DO  i = nxl_p, nxr_p
533
534          DO  j = nys_a, nyn_a
535             DO  k = 1, nz
536                ffty_ar3(j,k,i) = ar(k,j,i+nxl)
537             ENDDO
538          ENDDO
539
540          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'forward' )
541       ENDDO
542
543!$OMP  DO
544       DO  k = 1, nz
545          m = nys_a
546          DO  n = 1, npe_s
547             DO  j = nys_p, nyn_p
548                 DO  i = nxl_p, nxr_p
549                     work1(i,k,j,n) = ffty_ar3(m,k,i)
550                 ENDDO
551                 m = m+1
552             ENDDO
553          ENDDO
554       ENDDO
555!$OMP  END PARALLEL
556
557       CALL cpu_log( log_point_s(7), 'fft_y_m', 'pause' )
558
559#if defined( __parallel )
560       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
561       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
562                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
563                          comm2d, istat )
564       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
565#else
566       work2 = work1
567#endif
568
569       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'start' )
570
571#if defined( __KKMP )
572!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,tri_ar,jthread)
573!$OMP  DO
574       DO  j = nys_p, nyn_p
575          jthread = omp_get_thread_num() + 1
576#else
577       DO  j = nys_p, nyn_p
578          jthread = 1
579#endif
580          DO  k = 1, nz
581
582             m = nxl_a
583             DO  n = 1, npe_s
584                DO  i = nxl_p, nxr_p
585                   fftx_ar3(m,k,j) = work2(i,k,j,n)
586                   m = m+1
587                ENDDO
588             ENDDO
589          ENDDO
590
591          CALL fft_x_m( fftx_ar3(:,:,j), 'forward' )
592
593          DO  k = 1, nz
594             DO  i = nxl_a, nxr_a
595                tri_ar(i,k) = fftx_ar3(i,k,j)
596             ENDDO
597          ENDDO
598 
599          jj = myid * (nyn_p-nys_p+1) + j
600          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
601
602          DO  k = 1, nz
603             DO  i = nxl_a, nxr_a
604                fftx_ar3(i,k,j) = tri_ar (i,k)
605             ENDDO
606          ENDDO
607
608          CALL fft_x_m( fftx_ar3(:,:,j), 'backward' )
609
610          DO  k = 1, nz
611             m = nxl_a
612             DO  n = 1, npe_s
613                DO  i = nxl_p, nxr_p
614                   work2(i,k,j,n) = fftx_ar3(m,k,j)
615                   m = m+1
616                ENDDO
617             ENDDO
618          ENDDO
619
620       ENDDO
621#if defined( __KKMP )
622!$OMP  END PARALLEL
623#endif
624
625       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'stop' )
626
627#if defined( __parallel )
628       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
629       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
630       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
631                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
632                          comm2d, istat )
633       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
634#else
635       work1 = work2
636#endif
637
638       CALL cpu_log( log_point_s(7), 'fft_y_m', 'continue' )
639
640!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n)
641!$OMP  DO
642       DO  k = 1, nz
643          m = nys_a
644          DO  n = 1, npe_s
645             DO  j = nys_p, nyn_p
646                 DO  i = nxl_p, nxr_p
647                     ffty_ar3(m,k,i) = work1(i,k,j,n)
648                 ENDDO
649                 m = m+1
650             ENDDO
651          ENDDO
652       ENDDO 
653
654!$OMP  DO
655       DO  i = nxl_p, nxr_p
656          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'backward' )
657          DO  j = nys_a, nyn_a
658             DO  k = 1, nz
659                ar(k,j,i+nxl) = ffty_ar3(j,k,i)
660             ENDDO
661          ENDDO
662       ENDDO
663!$OMP  END PARALLEL
664
665       CALL cpu_log( log_point_s(7), 'fft_y_m', 'stop' )
666
667       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'stop' )
668
669#if defined( __KKMP )
670       DEALLOCATE( tri )
671#endif
672
673    END SUBROUTINE poisfft_hybrid_omp_vec
674
675
676    SUBROUTINE poisfft_hybrid_nodes ( ar )
677
678       USE cpulog
679       USE interfaces
680
681       IMPLICIT NONE
682
683       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
684       INTEGER            ::  i, iei, ii, iouter, ir, istat, j, jj, k, m, &
685                              n, nn, nt, nw1, nw2
686
687       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
688
689       REAL, DIMENSION(0:nx)                 ::  fftx_ar
690       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
691 
692       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
693
694       REAL, DIMENSION(nxl_p:nxr_p,nz,tasks_per_logical_node, &
695                          nodes,nys_p:nyn_p) ::  work1,work2
696       REAL, DIMENSION(5,0:nx,0:nz-1)        ::  tri
697
698
699       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'start' )
700
701       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
702
703!
704!--    Store grid points to be transformed on a 1d-array, do the fft
705!--    and sample the results on a 4d-array
706       DO  iouter = nxl_p, nxr_p, istride  ! stride loop, better cache
707          iei = MIN( iouter+istride-1, nxr_p )
708          DO  k = 1, nz
709
710             DO  i = iouter, iei
711                ii = nxl + i
712                ir = i - iouter + 1
713
714                DO  j = nys_a, nyn_a
715                   ffty_ar(j,ir) = ar(k,j,ii)
716                ENDDO
717   
718                CALL fft_y( ffty_ar(:,ir), 'forward' )
719             ENDDO
720
721             m = nys_a
722             DO  nn = 1, nodes
723                DO  nt = 1, tasks_per_logical_node
724                   DO  j = nys_p, nyn_p
725                      DO  i = iouter, iei
726                         ir = i - iouter + 1
727                         work1(i,k,nt,nn,j) = ffty_ar(m,ir)
728                      ENDDO
729                      m = m+1
730                   ENDDO
731                ENDDO
732             ENDDO
733   
734          ENDDO
735       ENDDO
736
737       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
738
739       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' )
740       nw1 = SIZE( work1, 1 ) * SIZE( work1, 2 )
741       DO  nn = 1, nodes
742           DO  j = nys_p, nyn_p
743#if defined( __parallel )
744              CALL MPI_ALLTOALL( work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
745                                 work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
746                                                     comm_tasks, istat )
747#endif
748           ENDDO
749       ENDDO
750       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
751
752
753       DO  j = nys_p, nyn_p
754       
755          CALL cascade( 1, j, nys_p, nyn_p )
756          nw2 = nw1 * SIZE( work1, 3 )
757          CALL cpu_log( log_point_s(37), 'alltoall_node', 'start' )
758#if defined( __parallel )
759          CALL MPI_ALLTOALL( work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
760                             work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
761                                                  comm_nodes, istat )
762#endif
763          CALL cpu_log( log_point_s(37), 'alltoall_node', 'pause' )
764          CALL cascade( 2, j, nys_p, nyn_p )
765
766          CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
767          DO  k = 1, nz
768
769             m = nxl_a
770             DO  nn = 1, nodes
771                DO  nt = 1, tasks_per_logical_node
772                   DO  i = nxl_p, nxr_p
773                      fftx_ar(m) = work1(i,k,nt,nn,j)
774                      m = m+1
775                   ENDDO
776                ENDDO
777             ENDDO
778
779             CALL fft_x( fftx_ar, 'forward' )
780
781             DO  i = nxl_a, nxr_a
782                tri_ar(i,k) = fftx_ar(i)
783             ENDDO
784
785          ENDDO
786 
787          jj = myid * (nyn_p-nys_p+1) + j
788          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:) )
789
790          DO  k = 1, nz
791             DO  i = nxl_a, nxr_a
792                fftx_ar(i) = tri_ar(i,k)
793             ENDDO
794
795             CALL fft_x( fftx_ar, 'backward' )
796
797             m = nxl_a
798             DO  nn = 1, nodes
799                DO  nt = 1, tasks_per_logical_node
800                   DO  i = nxl_p, nxr_p
801                      work1(i,k,nt,nn,j) = fftx_ar(m)
802                      m = m+1
803                   ENDDO
804                ENDDO
805             ENDDO
806          ENDDO
807
808          CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
809          nw2 = nw1 * SIZE( work1, 3 )
810          CALL cpu_log( log_point_s(37), 'alltoall_node', 'continue' )
811#if defined( __parallel )
812          CALL MPI_ALLTOALL( work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
813                             work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
814                                                  comm_nodes, istat )
815#endif
816          CALL cpu_log( log_point_s(37), 'alltoall_node', 'stop' )
817
818       ENDDO
819
820       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' ) 
821       DO  nn = 1, nodes
822           DO  j = nys_p, nyn_p
823#if defined( __parallel )
824              CALL MPI_ALLTOALL( work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
825                                 work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
826                                                     comm_tasks, istat )
827#endif
828           ENDDO
829       ENDDO
830       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
831
832       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
833
834       DO  iouter = nxl_p, nxr_p, istride
835          iei = MIN( iouter+istride-1, nxr_p )
836          DO  k = 1, nz
837
838             m = nys_a
839             DO  nn = 1, nodes
840                DO  nt = 1, tasks_per_logical_node
841                   DO  j = nys_p, nyn_p
842                      DO  i = iouter, iei
843                         ir = i - iouter + 1
844                         ffty_ar(m,ir) = work1(i,k,nt,nn,j)
845                      ENDDO
846                      m = m+1
847                   ENDDO
848                ENDDO
849             ENDDO
850
851             DO  i = iouter, iei
852                ii = nxl + i
853                ir = i - iouter + 1
854                CALL fft_y( ffty_ar(:,ir), 'backward' )
855
856                DO  j = nys_a, nyn_a
857                   ar(k,j,ii) = ffty_ar(j,ir)
858                ENDDO
859             ENDDO
860
861          ENDDO
862       ENDDO
863
864       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
865
866       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'stop' )
867
868    END SUBROUTINE poisfft_hybrid_nodes
869
870
871
872    SUBROUTINE tridia_hybrid( j, ar, tri )
873
874       USE arrays_3d
875       USE control_parameters
876       USE grid_variables
877
878       IMPLICIT NONE
879
880       INTEGER ::  i, j, k, nnyh
881       REAL, DIMENSION(0:nx,nz)       ::  ar
882       REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
883       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
884
885       nnyh = (ny+1) / 2
886
887       tri = 0.0
888!
889!--    Define constant elements of the tridiagonal matrix.
890       DO  k = 0, nz-1
891          DO  i = 0,nx
892             tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
893             tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
894          ENDDO
895       ENDDO
896
897       IF ( j <= nnyh )  THEN
898          CALL maketri_hybrid( j )
899       ELSE
900          CALL maketri_hybrid( ny+1-j)
901       ENDIF
902       CALL zerleg_hybrid
903       CALL substi_hybrid( ar, tri )
904
905    CONTAINS
906
907       SUBROUTINE maketri_hybrid( j )
908
909!----------------------------------------------------------------------!
910!                         maketri                                      !
911!                                                                      !
912!       computes the i- and j-dependent component of the matrix        !
913!----------------------------------------------------------------------!
914
915          USE constants
916
917          IMPLICIT NONE
918
919          INTEGER ::  i, j, k, nnxh
920          REAL    ::  a, c
921
922          REAL, DIMENSION(0:nx) ::  l
923
924
925          nnxh = (nx+1) / 2
926!
927!--       Provide the tridiagonal matrix for solution of the Poisson equation
928!--       in Fourier space. The coefficients are computed following the method
929!--       of Schmidt et al. (DFVLR-Mitteilung 84-15) --> departs from Stephan
930!--       Siano's original version.
931          DO  i = 0,nx
932             IF ( i >= 0 .AND. i < nnxh ) THEN
933                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
934                                        REAL( nx+1 ) ) ) / ( dx * dx ) + &
935                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
936                                        REAL( ny+1 ) ) ) / ( dy * dy )
937             ELSEIF ( i == nnxh ) THEN
938                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
939                                         REAL( nx+1 ) ) ) / ( dx * dx ) + &
940                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
941                                         REAL(ny+1) ) ) / ( dy * dy )
942             ELSE
943                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
944                                         REAL( nx+1 ) ) ) / ( dx * dx ) + &
945                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
946                                         REAL( ny+1 ) ) ) / ( dy * dy )
947             ENDIF
948          ENDDO
949
950          DO  k = 0,nz-1
951             DO  i = 0, nx
952                a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
953                c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
954                tri(1,i,k) = a + c - l(i)
955             ENDDO
956          ENDDO
957          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
958             DO  i = 0,nx
959                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
960             ENDDO
961          ENDIF
962          IF ( ibc_p_t == 1 )  THEN
963             DO  i = 0,nx
964                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
965             ENDDO
966          ENDIF
967
968       END SUBROUTINE maketri_hybrid
969
970
971       SUBROUTINE zerleg_hybrid
972
973!----------------------------------------------------------------------!
974!                          zerleg                                      !
975!                                                                      !
976!       Splitting of the tridiagonal matrix (Thomas algorithm)         !
977!----------------------------------------------------------------------!
978
979          USE indices
980
981          IMPLICIT NONE
982
983          INTEGER ::  i, k
984
985!
986!--       Splitting
987          DO  i = 0, nx
988             tri(4,i,0) = tri(1,i,0)
989          ENDDO
990          DO  k = 1, nz-1
991             DO  i = 0,nx
992                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
993                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
994             ENDDO
995          ENDDO
996
997       END SUBROUTINE zerleg_hybrid
998
999       SUBROUTINE substi_hybrid( ar, tri )
1000
1001!----------------------------------------------------------------------!
1002!                          substi                                      !
1003!                                                                      !
1004!       Substitution (Forward and Backward) (Thomas algorithm)         !
1005!----------------------------------------------------------------------!
1006
1007          IMPLICIT NONE
1008
1009          INTEGER ::  i, j, k
1010          REAL, DIMENSION(0:nx,nz)       ::  ar
1011          REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
1012          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1013
1014!
1015!--       Forward substitution
1016          DO  i = 0, nx
1017             ar1(i,0) = ar(i,1)
1018          ENDDO
1019          DO  k = 1, nz - 1
1020             DO  i = 0,nx
1021                ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
1022             ENDDO
1023          ENDDO
1024
1025!
1026!--       Backward substitution
1027          DO  i = 0,nx
1028             ar(i,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
1029          ENDDO
1030          DO  k = nz-2, 0, -1
1031             DO  i = 0,nx
1032                ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
1033                            / tri(4,i,k)
1034             ENDDO
1035          ENDDO
1036
1037       END SUBROUTINE substi_hybrid
1038
1039    END SUBROUTINE tridia_hybrid
1040
1041
1042    SUBROUTINE cascade( loca, j, nys_p, nyn_p )
1043
1044       USE cpulog
1045
1046       IMPLICIT NONE
1047
1048       INTEGER ::  ier, j, loca, nyn_p, nys_p, req, reqa(1)
1049       INTEGER, SAVE ::  tag = 10
1050#if defined( __parallel )
1051       INTEGER, DIMENSION(MPI_STATUS_SIZE)   :: stat
1052       INTEGER, DIMENSION(MPI_STATUS_SIZE,1) :: stata
1053#endif
1054
1055       REAL ::  buf, buf1
1056
1057
1058       buf  = 1.0
1059       buf1 = 1.1
1060       IF ( me_node == 0 )  THEN  ! first node only
1061
1062          SELECT CASE ( loca )
1063
1064             CASE ( 1 )  ! before alltoall
1065
1066                IF( me_task > 0 )  THEN  ! first task does not wait
1067#if defined( __parallel )
1068                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task-1, 0, &
1069                                      buf1, 1, MPI_REAL, me_task-1, 0, &
1070                                                 comm_tasks, stat, ierr )
1071#endif
1072                ELSEIF ( j > nys_p )  THEN
1073                   req = 0
1074                   tag = MOD( tag-10, 10 ) + 10
1075#if defined( __parallel )
1076                   CALL MPI_IRECV( buf, 1, MPI_REAL, tasks_per_logical_node-1,&
1077                                   tag, comm_tasks, req, ierr )
1078                   reqa = req
1079                   CALL MPI_WAITALL( 1, reqa, stata, ierr )
1080#endif
1081                ENDIF
1082
1083             CASE ( 2 )  ! after alltoall
1084
1085                IF ( me_task < tasks_per_logical_node-1 )  THEN  ! last task
1086#if defined( __parallel )
1087                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task+1, 0, &
1088                                      buf1, 1, MPI_REAL, me_task+1, 0, &
1089                                                 comm_tasks, stat, ierr)
1090#endif
1091                ELSEIF ( j < nyn_p )  THEN
1092                   req = 0
1093                   tag = MOD( tag-10, 10 ) + 10
1094#if defined( __parallel )
1095                   CALL MPI_ISEND( buf, 1, MPI_REAL, 0, tag, comm_tasks, req, &
1096                                   ierr )
1097#endif
1098                ENDIF
1099
1100          END SELECT
1101
1102       ENDIF
1103
1104    END SUBROUTINE cascade
1105#endif
1106 END MODULE poisfft_hybrid_mod
Note: See TracBrowser for help on using the repository browser.