Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |