Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |