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

Last change on this file since 1222 was 1222, checked in by raasch, 8 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 raasch $
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^2*nnx^2)+j^2/(dy^2*nny^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.