Changeset 1682 for palm/trunk/SOURCE/poisfft.f90
- Timestamp:
- Oct 7, 2015 11:56:08 PM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/poisfft.f90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/poisfft.f90
r1483 r1682 1 MODULE poisfft_mod 2 1 !> @file poisfft.f90 3 2 !--------------------------------------------------------------------------------! 4 3 ! This file is part of PALM. … … 20 19 ! Current revisions: 21 20 ! ----------------- 22 ! 21 ! Code annotations made doxygen readable 23 22 ! 24 23 ! Former revisions: … … 116 115 ! Description: 117 116 ! ------------ 118 ! Original version by Stephan Siano (pois3d), as of July 23, 1996 119 ! Adapted for 2D-domain-decomposition by Siegfried Raasch, July 3, 1997 120 ! 121 ! Solves the Poisson equation with a 2D spectral method 122 ! d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s 123 ! 124 ! Input: 125 ! real ar contains (nnz,nny,nnx) elements of the velocity divergence, 126 ! starting from (1,nys,nxl) 127 ! 128 ! Output: 129 ! real ar contains the solution for perturbation pressure p 130 !------------------------------------------------------------------------------! 117 !> Solves the Poisson equation with a 2D spectral method 118 !> d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s 119 !> 120 !> Input: 121 !> real ar contains (nnz,nny,nnx) elements of the velocity divergence, 122 !> starting from (1,nys,nxl) 123 !> 124 !> Output: 125 !> real ar contains the solution for perturbation pressure p 126 !------------------------------------------------------------------------------! 127 MODULE poisfft_mod 128 131 129 132 130 USE fft_xy, & … … 169 167 CONTAINS 170 168 169 !------------------------------------------------------------------------------! 170 ! Description: 171 ! ------------ 172 !> @todo Missing subroutine description. 173 !------------------------------------------------------------------------------! 171 174 SUBROUTINE poisfft_init 172 175 … … 178 181 IMPLICIT NONE 179 182 180 INTEGER(iwp) :: k ! :183 INTEGER(iwp) :: k !< 181 184 182 185 … … 191 194 192 195 #if ! defined ( __check ) 196 !------------------------------------------------------------------------------! 197 ! Description: 198 ! ------------ 199 !> Two-dimensional Fourier Transformation in x- and y-direction. 200 !------------------------------------------------------------------------------! 193 201 SUBROUTINE poisfft( ar ) 194 202 … … 205 213 IMPLICIT NONE 206 214 207 INTEGER(iwp) :: ii ! :208 INTEGER(iwp) :: iind ! :209 INTEGER(iwp) :: inew ! :210 INTEGER(iwp) :: jj ! :211 INTEGER(iwp) :: jind ! :212 INTEGER(iwp) :: jnew ! :213 INTEGER(iwp) :: ki ! :214 INTEGER(iwp) :: kk ! :215 INTEGER(iwp) :: knew ! :216 INTEGER(iwp) :: n ! :217 INTEGER(iwp) :: nblk ! :218 INTEGER(iwp) :: nnx_y ! :219 INTEGER(iwp) :: nny_z ! :220 INTEGER(iwp) :: nnz_t ! :221 INTEGER(iwp) :: nnz_x ! :222 INTEGER(iwp) :: nxl_y_bound ! :223 INTEGER(iwp) :: nxr_y_bound ! :224 225 INTEGER(iwp), DIMENSION(4) :: isave ! :226 227 REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) :: ar ! :228 REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv ! :215 INTEGER(iwp) :: ii !< 216 INTEGER(iwp) :: iind !< 217 INTEGER(iwp) :: inew !< 218 INTEGER(iwp) :: jj !< 219 INTEGER(iwp) :: jind !< 220 INTEGER(iwp) :: jnew !< 221 INTEGER(iwp) :: ki !< 222 INTEGER(iwp) :: kk !< 223 INTEGER(iwp) :: knew !< 224 INTEGER(iwp) :: n !< 225 INTEGER(iwp) :: nblk !< 226 INTEGER(iwp) :: nnx_y !< 227 INTEGER(iwp) :: nny_z !< 228 INTEGER(iwp) :: nnz_t !< 229 INTEGER(iwp) :: nnz_x !< 230 INTEGER(iwp) :: nxl_y_bound !< 231 INTEGER(iwp) :: nxr_y_bound !< 232 233 INTEGER(iwp), DIMENSION(4) :: isave !< 234 235 REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) :: ar !< 236 REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv !< 229 237 !$acc declare create( ar_inv ) 230 238 231 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ar1 ! :232 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_in ! :233 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_inv ! :234 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_y ! :235 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_z ! :239 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ar1 !< 240 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_in !< 241 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_inv !< 242 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_y !< 243 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_z !< 236 244 237 245 … … 708 716 709 717 710 718 !------------------------------------------------------------------------------! 719 ! Description: 720 ! ------------ 721 !> Fourier-transformation along y with subsequent transposition y --> x for 722 !> a 1d-decomposition along x. 723 !> 724 !> @attention The performance of this routine is much faster on the NEC-SX6, 725 !> if the first index of work_ffty_vec is odd. Otherwise 726 !> memory bank conflicts may occur (especially if the index is a 727 !> multiple of 128). That's why work_ffty_vec is dimensioned as 728 !> 0:ny+1. 729 !> Of course, this will not work if users are using an odd number 730 !> of gridpoints along y. 731 !------------------------------------------------------------------------------! 711 732 SUBROUTINE ffty_tr_yx( f_in, f_out ) 712 733 713 !------------------------------------------------------------------------------!714 ! Fourier-transformation along y with subsequent transposition y --> x for715 ! a 1d-decomposition along x716 !717 ! ATTENTION: The performance of this routine is much faster on the NEC-SX6,718 ! if the first index of work_ffty_vec is odd. Otherwise719 ! memory bank conflicts may occur (especially if the index is a720 ! multiple of 128). That's why work_ffty_vec is dimensioned as721 ! 0:ny+1.722 ! Of course, this will not work if users are using an odd number723 ! of gridpoints along y.724 !------------------------------------------------------------------------------!725 734 726 735 USE control_parameters, & … … 736 745 IMPLICIT NONE 737 746 738 INTEGER(iwp) :: i ! :739 INTEGER(iwp) :: iend ! :740 INTEGER(iwp) :: iouter ! :741 INTEGER(iwp) :: ir ! :742 INTEGER(iwp) :: j ! :743 INTEGER(iwp) :: k ! :744 745 INTEGER(iwp), PARAMETER :: stridex = 4 ! :746 747 REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty ! :747 INTEGER(iwp) :: i !< 748 INTEGER(iwp) :: iend !< 749 INTEGER(iwp) :: iouter !< 750 INTEGER(iwp) :: ir !< 751 INTEGER(iwp) :: j !< 752 INTEGER(iwp) :: k !< 753 754 INTEGER(iwp), PARAMETER :: stridex = 4 !< 755 756 REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !< 748 757 #if defined( __nec ) 749 REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec ! :758 REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !< 750 759 #endif 751 REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_in ! :752 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_out ! :753 REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work ! :760 REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_in !< 761 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_out !< 762 REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !< 754 763 755 764 ! … … 844 853 845 854 855 !------------------------------------------------------------------------------! 856 ! Description: 857 ! ------------ 858 !> Transposition x --> y with a subsequent backward Fourier transformation for 859 !> a 1d-decomposition along x 860 !------------------------------------------------------------------------------! 846 861 SUBROUTINE tr_xy_ffty( f_in, f_out ) 847 862 848 !------------------------------------------------------------------------------!849 ! Transposition x --> y with a subsequent backward Fourier transformation for850 ! a 1d-decomposition along x851 !------------------------------------------------------------------------------!852 863 853 864 USE control_parameters, & … … 863 874 IMPLICIT NONE 864 875 865 INTEGER(iwp) :: i ! :866 INTEGER(iwp) :: iend ! :867 INTEGER(iwp) :: iouter ! :868 INTEGER(iwp) :: ir ! :869 INTEGER(iwp) :: j ! :870 INTEGER(iwp) :: k ! :871 872 INTEGER(iwp), PARAMETER :: stridex = 4 ! :873 874 REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty ! :876 INTEGER(iwp) :: i !< 877 INTEGER(iwp) :: iend !< 878 INTEGER(iwp) :: iouter !< 879 INTEGER(iwp) :: ir !< 880 INTEGER(iwp) :: j !< 881 INTEGER(iwp) :: k !< 882 883 INTEGER(iwp), PARAMETER :: stridex = 4 !< 884 885 REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !< 875 886 #if defined( __nec ) 876 REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec ! :887 REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !< 877 888 #endif 878 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_in ! :879 REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_out ! :880 REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work ! :889 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_in !< 890 REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_out !< 891 REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !< 881 892 882 893 ! … … 970 981 971 982 983 !------------------------------------------------------------------------------! 984 ! Description: 985 ! ------------ 986 !> FFT along x, solution of the tridiagonal system and backward FFT for 987 !> a 1d-decomposition along x 988 !> 989 !> @warning this subroutine may still not work for hybrid parallelization 990 !> with OpenMP (for possible necessary changes see the original 991 !> routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002) 992 !------------------------------------------------------------------------------! 972 993 SUBROUTINE fftx_tri_fftx( ar ) 973 994 974 !------------------------------------------------------------------------------!975 ! FFT along x, solution of the tridiagonal system and backward FFT for976 ! a 1d-decomposition along x977 !978 ! WARNING: this subroutine may still not work for hybrid parallelization979 ! with OpenMP (for possible necessary changes see the original980 ! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)981 !------------------------------------------------------------------------------!982 995 983 996 USE control_parameters, & … … 996 1009 IMPLICIT NONE 997 1010 998 INTEGER(iwp) :: i ! :999 INTEGER(iwp) :: j ! :1000 INTEGER(iwp) :: k ! :1001 INTEGER(iwp) :: m ! :1002 INTEGER(iwp) :: n ! :1003 INTEGER(iwp) :: omp_get_thread_num ! :1004 INTEGER(iwp) :: tn ! :1005 1006 REAL(wp), DIMENSION(0:nx) :: work_fftx ! :1007 REAL(wp), DIMENSION(0:nx,1:nz) :: work_trix ! :1008 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: ar ! :1009 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri ! :1011 INTEGER(iwp) :: i !< 1012 INTEGER(iwp) :: j !< 1013 INTEGER(iwp) :: k !< 1014 INTEGER(iwp) :: m !< 1015 INTEGER(iwp) :: n !< 1016 INTEGER(iwp) :: omp_get_thread_num !< 1017 INTEGER(iwp) :: tn !< 1018 1019 REAL(wp), DIMENSION(0:nx) :: work_fftx !< 1020 REAL(wp), DIMENSION(0:nx,1:nz) :: work_trix !< 1021 REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: ar !< 1022 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !< 1010 1023 1011 1024 … … 1113 1126 1114 1127 1128 !------------------------------------------------------------------------------! 1129 ! Description: 1130 ! ------------ 1131 !> Fourier-transformation along x with subsequent transposition x --> y for 1132 !> a 1d-decomposition along y. 1133 !> 1134 !> @attention NEC-branch of this routine may significantly profit from 1135 !> further optimizations. So far, performance is much worse than 1136 !> for routine ffty_tr_yx (more than three times slower). 1137 !------------------------------------------------------------------------------! 1115 1138 SUBROUTINE fftx_tr_xy( f_in, f_out ) 1116 1139 1117 !------------------------------------------------------------------------------!1118 ! Fourier-transformation along x with subsequent transposition x --> y for1119 ! a 1d-decomposition along y1120 !1121 ! ATTENTION: The NEC-branch of this routine may significantly profit from1122 ! further optimizations. So far, performance is much worse than1123 ! for routine ffty_tr_yx (more than three times slower).1124 !------------------------------------------------------------------------------!1125 1140 1126 1141 USE control_parameters, & … … 1136 1151 IMPLICIT NONE 1137 1152 1138 INTEGER(iwp) :: i ! :1139 INTEGER(iwp) :: j ! :1140 INTEGER(iwp) :: k ! :1141 1142 REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx ! :1143 REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_in ! :1144 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_out ! :1145 REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work ! :1153 INTEGER(iwp) :: i !< 1154 INTEGER(iwp) :: j !< 1155 INTEGER(iwp) :: k !< 1156 1157 REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !< 1158 REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_in !< 1159 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_out !< 1160 REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !< 1146 1161 1147 1162 ! … … 1228 1243 1229 1244 1245 !------------------------------------------------------------------------------! 1246 ! Description: 1247 ! ------------ 1248 !> Transposition y --> x with a subsequent backward Fourier transformation for 1249 !> a 1d-decomposition along x. 1250 !------------------------------------------------------------------------------! 1230 1251 SUBROUTINE tr_yx_fftx( f_in, f_out ) 1231 1252 1232 !------------------------------------------------------------------------------!1233 ! Transposition y --> x with a subsequent backward Fourier transformation for1234 ! a 1d-decomposition along x1235 !------------------------------------------------------------------------------!1236 1253 1237 1254 USE control_parameters, & … … 1247 1264 IMPLICIT NONE 1248 1265 1249 INTEGER(iwp) :: i ! :1250 INTEGER(iwp) :: j ! :1251 INTEGER(iwp) :: k ! :1252 1253 REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx ! :1254 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_in ! :1255 REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_out ! :1256 REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work ! :1266 INTEGER(iwp) :: i !< 1267 INTEGER(iwp) :: j !< 1268 INTEGER(iwp) :: k !< 1269 1270 REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !< 1271 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_in !< 1272 REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_out !< 1273 REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !< 1257 1274 1258 1275 ! … … 1335 1352 1336 1353 1354 !------------------------------------------------------------------------------! 1355 ! Description: 1356 ! ------------ 1357 !> FFT along y, solution of the tridiagonal system and backward FFT for 1358 !> a 1d-decomposition along y. 1359 !> 1360 !> @warning this subroutine may still not work for hybrid parallelization 1361 !> with OpenMP (for possible necessary changes see the original 1362 !> routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002) 1363 !------------------------------------------------------------------------------! 1337 1364 SUBROUTINE ffty_tri_ffty( ar ) 1338 1365 1339 !------------------------------------------------------------------------------!1340 ! FFT along y, solution of the tridiagonal system and backward FFT for1341 ! a 1d-decomposition along y1342 !1343 ! WARNING: this subroutine may still not work for hybrid parallelization1344 ! with OpenMP (for possible necessary changes see the original1345 ! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)1346 !------------------------------------------------------------------------------!1347 1366 1348 1367 USE control_parameters, & … … 1361 1380 IMPLICIT NONE 1362 1381 1363 INTEGER(iwp) :: i ! :1364 INTEGER(iwp) :: j ! :1365 INTEGER(iwp) :: k ! :1366 INTEGER(iwp) :: m ! :1367 INTEGER(iwp) :: n ! :1368 INTEGER(iwp) :: omp_get_thread_num ! :1369 INTEGER(iwp) :: tn ! :1370 1371 REAL(wp), DIMENSION(0:ny) :: work_ffty ! :1372 REAL(wp), DIMENSION(0:ny,1:nz) :: work_triy ! :1373 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: ar ! :1374 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri ! :1382 INTEGER(iwp) :: i !< 1383 INTEGER(iwp) :: j !< 1384 INTEGER(iwp) :: k !< 1385 INTEGER(iwp) :: m !< 1386 INTEGER(iwp) :: n !< 1387 INTEGER(iwp) :: omp_get_thread_num !< 1388 INTEGER(iwp) :: tn !< 1389 1390 REAL(wp), DIMENSION(0:ny) :: work_ffty !< 1391 REAL(wp), DIMENSION(0:ny,1:nz) :: work_triy !< 1392 REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: ar !< 1393 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !< 1375 1394 1376 1395
Note: See TracChangeset
for help on using the changeset viewer.
