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

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

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:


openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:


bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

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