Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |