Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |