Home

Context Navigation

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

Last change on this file since 4327 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 41.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |