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^2*nnx^2)+j^2/(dy^2*nny^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.