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

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

last commit documented

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