Home

Context Navigation

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

Last change on this file since 1024 was 1014, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 32.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |