Home

Context Navigation

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

Last change on this file since 4651 was 4649, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 45.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |