Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |