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

Last change on this file since 668 was 668, checked in by suehring, 14 years ago

last commit documented

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