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

Last change on this file since 1213 was 1213, checked in by raasch, 8 years ago

last commit documented

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