Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver.f90 @ 1220

Last change on this file since 1220 was 1217, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.9 KB

Line
1	MODULE tridia_solver
2
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	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: tridia_solver.f90 1217 2013-08-26 11:30:05Z heinze $
27	!
28	! 1216 2013-08-26 09:31:42Z raasch
29	! +tridia_substi_overlap for handling overlapping fft / transposition
30	!
31	! 1212 2013-08-15 08:46:27Z raasch
32	! Initial revision.
33	! Routines have been moved to seperate module from former file poisfft to here.
34	! The tridiagonal matrix coefficients of array tri are calculated only once at
35	! the beginning, i.e. routine split is called within tridia_init.
36	!
37	!
38	! Description:
39	! ------------
40	! solves the linear system of equations:
41	!
42	! -(4 pi^2(i^2/(dx^2nnx^2)+j^2/(dy^2nny^2))+
43	! 1/(dzu(k)dzw(k))+1/(dzu(k-1)dzw(k)))*p(i,j,k)+
44	! 1/(dzu(k)dzw(k))p(i,j,k+1)+1/(dzu(k-1)dzw(k))p(i,j,k-1)=d(i,j,k)
45	!
46	! by using the Thomas algorithm
47	!------------------------------------------------------------------------------!
48
49	USE indices
50	USE transpose_indices
51
52	IMPLICIT NONE
53
54	REAL, DIMENSION(:,:), ALLOCATABLE :: ddzuw
55
56	PRIVATE
57
58	INTERFACE tridia_substi
59	MODULE PROCEDURE tridia_substi
60	END INTERFACE tridia_substi
61
62	INTERFACE tridia_substi_overlap
63	MODULE PROCEDURE tridia_substi_overlap
64	END INTERFACE tridia_substi_overlap
65
66	PUBLIC tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
67
68	CONTAINS
69
70
71	SUBROUTINE tridia_init
72
73	USE arrays_3d, ONLY: ddzu_pres, ddzw
74
75	IMPLICIT NONE
76
77	INTEGER :: k
78
79	ALLOCATE( ddzuw(0:nz-1,3) )
80
81	DO k = 0, nz-1
82	ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
83	ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
84	ddzuw(k,3) = -1.0 * &
85	( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
86	ENDDO
87	!
88	!-- Calculate constant coefficients of the tridiagonal matrix
89	#if ! defined ( __check )
90	CALL maketri
91	CALL split
92	#endif
93
94	END SUBROUTINE tridia_init
95
96
97	SUBROUTINE maketri
98
99	!------------------------------------------------------------------------------!
100	! Computes the i- and j-dependent component of the matrix
101	!------------------------------------------------------------------------------!
102
103	USE arrays_3d, ONLY: tric
104	USE constants
105	USE control_parameters
106	USE grid_variables
107
108	IMPLICIT NONE
109
110	INTEGER :: i, j, k, nnxh, nnyh
111
112	!$acc declare create( ll )
113	REAL :: ll(nxl_z:nxr_z,nys_z:nyn_z)
114
115
116	nnxh = ( nx + 1 ) / 2
117	nnyh = ( ny + 1 ) / 2
118
119	!
120	!-- Provide the constant coefficients of the tridiagonal matrix for solution
121	!-- of the Poisson equation in Fourier space.
122	!-- The coefficients are computed following the method of
123	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
124	!-- Siano's original version by discretizing the Poisson equation,
125	!-- before it is Fourier-transformed.
126
127	!$acc kernels present( tric )
128	!$acc loop vector( 32 )
129	DO j = nys_z, nyn_z
130	DO i = nxl_z, nxr_z
131	IF ( j >= 0 .AND. j <= nnyh ) THEN
132	IF ( i >= 0 .AND. i <= nnxh ) THEN
133	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
134	REAL( nx+1 ) ) ) / ( dx * dx ) + &
135	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
136	REAL( ny+1 ) ) ) / ( dy * dy )
137	ELSE
138	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
139	REAL( nx+1 ) ) ) / ( dx * dx ) + &
140	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
141	REAL( ny+1 ) ) ) / ( dy * dy )
142	ENDIF
143	ELSE
144	IF ( i >= 0 .AND. i <= nnxh ) THEN
145	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
146	REAL( nx+1 ) ) ) / ( dx * dx ) + &
147	2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
148	REAL( ny+1 ) ) ) / ( dy * dy )
149	ELSE
150	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
151	REAL( nx+1 ) ) ) / ( dx * dx ) + &
152	2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
153	REAL( ny+1 ) ) ) / ( dy * dy )
154	ENDIF
155	ENDIF
156	ENDDO
157	ENDDO
158
159	!$acc loop
160	DO k = 0, nz-1
161	DO j = nys_z, nyn_z
162	!$acc loop vector( 32 )
163	DO i = nxl_z, nxr_z
164	tric(i,j,k) = ddzuw(k,3) - ll(i,j)
165	ENDDO
166	ENDDO
167	ENDDO
168	!$acc end kernels
169
170	IF ( ibc_p_b == 1 ) THEN
171	!$acc kernels present( tric )
172	!$acc loop
173	DO j = nys_z, nyn_z
174	DO i = nxl_z, nxr_z
175	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
176	ENDDO
177	ENDDO
178	!$acc end kernels
179	ENDIF
180	IF ( ibc_p_t == 1 ) THEN
181	!$acc kernels present( tric )
182	!$acc loop
183	DO j = nys_z, nyn_z
184	DO i = nxl_z, nxr_z
185	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
186	ENDDO
187	ENDDO
188	!$acc end kernels
189	ENDIF
190
191	END SUBROUTINE maketri
192
193
194	SUBROUTINE tridia_substi( ar )
195
196	!------------------------------------------------------------------------------!
197	! Substitution (Forward and Backward) (Thomas algorithm)
198	!------------------------------------------------------------------------------!
199
200	USE arrays_3d, ONLY: tri
201	USE control_parameters
202
203	IMPLICIT NONE
204
205	INTEGER :: i, j, k
206
207	REAL :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
208
209	!$acc declare create( ar1 )
210	REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1
211
212	!
213	!-- Forward substitution
214	DO k = 0, nz - 1
215	!$acc kernels present( ar, tri )
216	!$acc loop
217	DO j = nys_z, nyn_z
218	DO i = nxl_z, nxr_z
219
220	IF ( k == 0 ) THEN
221	ar1(i,j,k) = ar(i,j,k+1)
222	ELSE
223	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
224	ENDIF
225
226	ENDDO
227	ENDDO
228	!$acc end kernels
229	ENDDO
230
231	!
232	!-- Backward substitution
233	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
234	!-- by zero appearing if the pressure bc is set to neumann at the top of
235	!-- the model domain.
236	DO k = nz-1, 0, -1
237	!$acc kernels present( ar, tri )
238	!$acc loop
239	DO j = nys_z, nyn_z
240	DO i = nxl_z, nxr_z
241
242	IF ( k == nz-1 ) THEN
243	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
244	ELSE
245	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
246	/ tri(i,j,k,1)
247	ENDIF
248	ENDDO
249	ENDDO
250	!$acc end kernels
251	ENDDO
252
253	!
254	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
255	!-- The respective values of ar should be zero at all k-levels if
256	!-- acceleration of horizontally averaged vertical velocity is zero.
257	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
258	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
259	!$acc kernels loop present( ar )
260	DO k = 1, nz
261	ar(nxl_z,nys_z,k) = 0.0
262	ENDDO
263	ENDIF
264	ENDIF
265
266	END SUBROUTINE tridia_substi
267
268
269	SUBROUTINE tridia_substi_overlap( ar, jj )
270
271	!------------------------------------------------------------------------------!
272	! Substitution (Forward and Backward) (Thomas algorithm)
273	!------------------------------------------------------------------------------!
274
275	USE arrays_3d, ONLY: tri
276	USE control_parameters
277
278	IMPLICIT NONE
279
280	INTEGER :: i, j, jj, k
281
282	REAL :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
283
284	!$acc declare create( ar1 )
285	REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1
286
287	!
288	!-- Forward substitution
289	DO k = 0, nz - 1
290	!$acc kernels present( ar, tri )
291	!$acc loop
292	DO j = nys_z, nyn_z
293	DO i = nxl_z, nxr_z
294
295	IF ( k == 0 ) THEN
296	ar1(i,j,k) = ar(i,j,k+1)
297	ELSE
298	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
299	ENDIF
300
301	ENDDO
302	ENDDO
303	!$acc end kernels
304	ENDDO
305
306	!
307	!-- Backward substitution
308	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
309	!-- by zero appearing if the pressure bc is set to neumann at the top of
310	!-- the model domain.
311	DO k = nz-1, 0, -1
312	!$acc kernels present( ar, tri )
313	!$acc loop
314	DO j = nys_z, nyn_z
315	DO i = nxl_z, nxr_z
316
317	IF ( k == nz-1 ) THEN
318	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
319	ELSE
320	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
321	/ tri(i,jj,k,1)
322	ENDIF
323	ENDDO
324	ENDDO
325	!$acc end kernels
326	ENDDO
327
328	!
329	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
330	!-- The respective values of ar should be zero at all k-levels if
331	!-- acceleration of horizontally averaged vertical velocity is zero.
332	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
333	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
334	!$acc kernels loop present( ar )
335	DO k = 1, nz
336	ar(nxl_z,nys_z,k) = 0.0
337	ENDDO
338	ENDIF
339	ENDIF
340
341	END SUBROUTINE tridia_substi_overlap
342
343
344	SUBROUTINE split
345
346	!------------------------------------------------------------------------------!
347	! Splitting of the tridiagonal matrix (Thomas algorithm)
348	!------------------------------------------------------------------------------!
349
350	USE arrays_3d, ONLY: tri, tric
351
352	IMPLICIT NONE
353
354	INTEGER :: i, j, k
355
356	!
357	!-- Splitting
358	!$acc kernels present( tri, tric )
359	!$acc loop
360	DO j = nys_z, nyn_z
361	!$acc loop vector( 32 )
362	DO i = nxl_z, nxr_z
363	tri(i,j,0,1) = tric(i,j,0)
364	ENDDO
365	ENDDO
366	!$acc end kernels
367
368	DO k = 1, nz-1
369	!$acc kernels present( tri, tric )
370	!$acc loop
371	DO j = nys_z, nyn_z
372	!$acc loop vector( 32 )
373	DO i = nxl_z, nxr_z
374	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
375	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
376	ENDDO
377	ENDDO
378	!$acc end kernels
379	ENDDO
380
381	END SUBROUTINE split
382
383
384	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
385
386	!------------------------------------------------------------------------------!
387	! Solves the linear system of equations for a 1d-decomposition along x (see
388	! tridia)
389	!
390	! Attention: when using the intel compilers older than 12.0, array tri must
391	! be passed as an argument to the contained subroutines. Otherwise
392	! addres faults will occur. This feature can be activated with
393	! cpp-switch __intel11
394	! On NEC, tri should not be passed (except for routine substi_1dd)
395	! because this causes very bad performance.
396	!------------------------------------------------------------------------------!
397
398	USE arrays_3d
399	USE control_parameters
400
401	USE pegrid
402
403	IMPLICIT NONE
404
405	INTEGER :: i, j, k, nnyh, nx, ny, omp_get_thread_num, tn
406
407	REAL :: ddx2, ddy2
408
409	REAL, DIMENSION(0:nx,1:nz) :: ar
410	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
411
412
413	nnyh = ( ny + 1 ) / 2
414
415	!
416	!-- Define constant elements of the tridiagonal matrix.
417	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
418	!-- the exchanged loops create bank conflicts. The following directive
419	!-- prohibits loop exchange and the loops perform much better.
420	! tn = omp_get_thread_num()
421	! WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
422	! CALL local_flush( 120+tn )
423	!CDIR NOLOOPCHG
424	DO k = 0, nz-1
425	DO i = 0,nx
426	tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
427	tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
428	ENDDO
429	ENDDO
430	! WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop thread=', omp_get_thread_num()
431	! CALL local_flush( 120+tn )
432
433	IF ( j <= nnyh ) THEN
434	#if defined( __intel11 )
435	CALL maketri_1dd( j, tri )
436	#else
437	CALL maketri_1dd( j )
438	#endif
439	ELSE
440	#if defined( __intel11 )
441	CALL maketri_1dd( ny+1-j, tri )
442	#else
443	CALL maketri_1dd( ny+1-j )
444	#endif
445	ENDIF
446	#if defined( __intel11 )
447	CALL split_1dd( tri )
448	#else
449	CALL split_1dd
450	#endif
451	CALL substi_1dd( ar, tri )
452
453	CONTAINS
454
455	#if defined( __intel11 )
456	SUBROUTINE maketri_1dd( j, tri )
457	#else
458	SUBROUTINE maketri_1dd( j )
459	#endif
460
461	!------------------------------------------------------------------------------!
462	! computes the i- and j-dependent component of the matrix
463	!------------------------------------------------------------------------------!
464
465	USE constants
466
467	IMPLICIT NONE
468
469	INTEGER :: i, j, k, nnxh
470	REAL :: a, c
471
472	REAL, DIMENSION(0:nx) :: l
473
474	#if defined( __intel11 )
475	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
476	#endif
477
478
479	nnxh = ( nx + 1 ) / 2
480	!
481	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
482	!-- Fourier space. The coefficients are computed following the method of
483	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
484	!-- Siano's original version by discretizing the Poisson equation,
485	!-- before it is Fourier-transformed
486	DO i = 0, nx
487	IF ( i >= 0 .AND. i <= nnxh ) THEN
488	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
489	REAL( nx+1 ) ) ) * ddx2 + &
490	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
491	REAL( ny+1 ) ) ) * ddy2
492	ELSE
493	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
494	REAL( nx+1 ) ) ) * ddx2 + &
495	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
496	REAL( ny+1 ) ) ) * ddy2
497	ENDIF
498	ENDDO
499
500	DO k = 0, nz-1
501	DO i = 0, nx
502	a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
503	c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
504	tri(1,i,k) = a + c - l(i)
505	ENDDO
506	ENDDO
507	IF ( ibc_p_b == 1 ) THEN
508	DO i = 0, nx
509	tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
510	ENDDO
511	ENDIF
512	IF ( ibc_p_t == 1 ) THEN
513	DO i = 0, nx
514	tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
515	ENDDO
516	ENDIF
517
518	END SUBROUTINE maketri_1dd
519
520
521	#if defined( __intel11 )
522	SUBROUTINE split_1dd( tri )
523	#else
524	SUBROUTINE split_1dd
525	#endif
526
527	!------------------------------------------------------------------------------!
528	! Splitting of the tridiagonal matrix (Thomas algorithm)
529	!------------------------------------------------------------------------------!
530
531	IMPLICIT NONE
532
533	INTEGER :: i, k
534
535	#if defined( __intel11 )
536	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
537	#endif
538
539
540	!
541	!-- Splitting
542	DO i = 0, nx
543	tri(4,i,0) = tri(1,i,0)
544	ENDDO
545	DO k = 1, nz-1
546	DO i = 0, nx
547	tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
548	tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
549	ENDDO
550	ENDDO
551
552	END SUBROUTINE split_1dd
553
554
555	SUBROUTINE substi_1dd( ar, tri )
556
557	!------------------------------------------------------------------------------!
558	! Substitution (Forward and Backward) (Thomas algorithm)
559	!------------------------------------------------------------------------------!
560
561	IMPLICIT NONE
562
563	INTEGER :: i, k
564
565	REAL, DIMENSION(0:nx,nz) :: ar
566	REAL, DIMENSION(0:nx,0:nz-1) :: ar1
567	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
568
569	!
570	!-- Forward substitution
571	DO i = 0, nx
572	ar1(i,0) = ar(i,1)
573	ENDDO
574	DO k = 1, nz-1
575	DO i = 0, nx
576	ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
577	ENDDO
578	ENDDO
579
580	!
581	!-- Backward substitution
582	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
583	!-- by zero appearing if the pressure bc is set to neumann at the top of
584	!-- the model domain.
585	DO i = 0, nx
586	ar(i,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
587	ENDDO
588	DO k = nz-2, 0, -1
589	DO i = 0, nx
590	ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
591	/ tri(4,i,k)
592	ENDDO
593	ENDDO
594
595	!
596	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
597	!-- The respective values of ar should be zero at all k-levels if
598	!-- acceleration of horizontally averaged vertical velocity is zero.
599	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
600	IF ( j == 0 ) THEN
601	DO k = 1, nz
602	ar(0,k) = 0.0
603	ENDDO
604	ENDIF
605	ENDIF
606
607	END SUBROUTINE substi_1dd
608
609	END SUBROUTINE tridia_1dd
610
611
612	END MODULE tridia_solver

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |