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

Last change on this file since 1258 was 1258, 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 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^2*nnx^2)+j^2/(dy^2*nny^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.