Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 1212 for palm/trunk/SOURCE


Ignore:
Timestamp:
Aug 15, 2013 8:46:27 AM (11 years ago)
Author:
raasch
Message:

tridia-solver moved to seperate module; the tridiagonal matrix coefficients of array tri are calculated only once at the beginning

Location:
palm/trunk/SOURCE
Files:
1 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r1211 r1212  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# +tridia_solver, -poisfft_hybrid
    2323#
    2424# Former revisions:
    … …  
    157157        message.f90 microphysics.f90 modules.f90 netcdf.f90 package_parin.f90 \
    158158        palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 \
    159         poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
     159        poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
    160160        production_e.f90 prognostic_equations.f90 random_function.f90 \
    161161        random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
    … …  
    164164        surface_coupler.f90 swap_timelevel.f90 temperton_fft.f90 \
    165165        time_integration.f90 time_to_string.f90 timestep.f90 \
    166         timestep_scheme_steering.f90 transpose.f90 user_3d_data_averaging.f90 \
     166        timestep_scheme_steering.f90 transpose.f90 tridia_solver.f90 \
     167        user_3d_data_averaging.f90 \
    167168        user_actions.f90 user_additional_routines.f90 \
    168169        user_check_data_output.f90 user_check_data_output_pr.f90 \
    … …  
    307308parin.o: modules.o
    308309plant_canopy_model.o: modules.o
    309 poisfft.o: modules.o fft_xy.o
    310 poisfft_hybrid.o: modules.o fft_xy.o
     310poisfft.o: modules.o fft_xy.o tridia_solver.o
    311311poismg.o: modules.o
    312312prandtl_fluxes.o: modules.o
    313 pres.o: modules.o poisfft.o poisfft_hybrid.o
     313pres.o: modules.o poisfft.o
    314314print_1d.o: modules.o
    315315production_e.o: modules.o wall_fluxes.o
    … …  
    337337timestep_scheme_steering.o: modules.o
    338338transpose.o: modules.o
     339tridia_solver.o: modules.o
    339340user_3d_data_averaging.o: modules.o user_module.o
    340341user_actions.o: modules.o user_module.o

  • palm/trunk/SOURCE/Makefile_check

    r1112 r1212  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# object file list replaced by one line statement,
     23# -poisfft_hybrid
    2324#
    2425# Former revisions:
    … …  
    6162PROG = check_namelist_files.x
    6263
    63 RCS = check_open.f90 check_namelist_files.f90 check_parameters.f90 \
     64SOURCES = check_open.f90 check_namelist_files.f90 check_parameters.f90 \
    6465      close_file.f90 cpu_log.f90 cuda_fft_interfaces.f90 exchange_horiz.f90 \
    6566      exchange_horiz_2d.f90 fft_xy.f90 init_grid.f90 init_masks.f90 \
    6667      init_cloud_physics.f90 init_pegrid.f90 local_flush.f90 local_stop.f90 \
    6768      local_system.f90 message.f90 modules.f90 package_parin.f90 parin.f90 \
    68       poisfft.f90 poisfft_hybrid.f90 random_function.f90 singleton.f90 \
     69      poisfft.f90 random_function.f90 singleton.f90 \
    6970      subsidence.f90 temperton_fft.f90 \
    7071      user_3d_data_averaging.f90 user_actions.f90 \
    … …  
    8182
    8283
    83 
    84 
    85 OBJS = check_open.o check_namelist_files.o check_parameters.o close_file.o \
    86        cpu_log.o cuda_fft_interfaces.o exchange_horiz.o exchange_horiz_2d.o \
    87        fft_xy.o init_grid.o init_masks.o init_pegrid.o init_cloud_physics.o\
    88        local_flush.o local_stop.o local_system.o message.o \
    89        modules.o package_parin.o parin.o poisfft.o \
    90        poisfft_hybrid.o random_function.o singleton.o subsidence.o temperton_fft.o \
    91        user_3d_data_averaging.o user_actions.o user_additional_routines.o \
    92        user_check_data_output.o user_check_data_output_pr.o \
    93        user_check_parameters.o user_data_output_2d.o user_data_output_3d.o \
    94        user_data_output_mask.o user_data_output_dvrp.o \
    95        user_define_netcdf_grid.o user_dvrp_coltab.o user_header.o \
    96        user_init.o user_init_3d_model.o user_init_grid.o \
    97        user_init_plant_canopy.o user_last_actions.o user_lpm_advec.o \
    98        user_lpm_init.o user_lpm_set_attributes.o user_module.o user_parin.o \
    99        user_read_restart_data.o user_spectra.o user_statistics.o \
     84OBJS=$(SOURCES:.f90=.o)
    10085
    10186CC = cc
    … …  
    144129parin.o: modules.o
    145130poisfft.o: modules.o fft_xy.o
    146 poisfft_hybrid.o: modules.o fft_xy.o
    147131random_function.o: modules.o
    148132singleton.o: singleton.f90

  • palm/trunk/SOURCE/check_parameters.f90

    r1211 r1212  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! checks for poisfft_hybrid removed
    2323!
    2424! Former revisions:
    … …  
    674674!-- Check whether there are any illegal values
    675675!-- Pressure solver:
    676     IF ( psolver /= 'poisfft'  .AND.  psolver /= 'poisfft_hybrid'  .AND. &
     676    IF ( psolver /= 'poisfft'  .AND. &
    677677         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
    678678       message_string = 'unknown solver for perturbation pressure: psolver' // &
    … …  
    680680       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
    681681    ENDIF
    682 
    683 #if defined( __parallel )
    684     IF ( psolver == 'poisfft_hybrid'  .AND.  pdims(2) /= 1 )  THEN
    685        message_string = 'psolver = "' // TRIM( psolver ) // '" only works ' // &
    686                         'for a 1d domain-decomposition along x & please do' // &
    687                         ' not set npey/=1 in the parameter file'
    688        CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
    689     ENDIF
    690 #else
    691     IF ( psolver == 'poisfft_hybrid' )  THEN
    692        message_string = 'psolver = "' // TRIM( psolver ) // '" only works' // &
    693                         ' for a parallel environment'
    694        CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
    695     ENDIF
    696 #endif
    697682
    698683    IF ( psolver == 'multigrid' )  THEN

  • palm/trunk/SOURCE/header.f90

    r1182 r1212  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! output for poisfft_hybrid removed
    2323!
    2424! Former revisions:
    … …  
    329329    IF ( psolver(1:7) == 'poisfft' )  THEN
    330330       WRITE ( io, 111 )  TRIM( fft_method )
    331        IF ( psolver == 'poisfft_hybrid' )  WRITE ( io, 138 )
    332331    ELSEIF ( psolver == 'sor' )  THEN
    333332       WRITE ( io, 112 )  nsor_ini, nsor, omega_sor
    … …  
    16801679            '     gridpoints of coarsest domain (x,y,z):    (',I3,',',I3,',', &
    16811680                  I3,')')
    1682 138 FORMAT ('     Using hybrid version for 1d-domain-decomposition')
    16831681139 FORMAT (' --> Loop optimization method: ',A)
    16841682140 FORMAT ('     maximum residual allowed:                ',E10.3)

  • palm/trunk/SOURCE/init_3d_model.f90

    r1196 r1212  
    2323! Current revisions:
    2424! ------------------
    25 !
     25! array tri is allocated and included in data copy statement
    2626!
    2727! Former revisions:
    … …  
    364364!-- Array for storing constant coeffficients of the tridiagonal solver
    365365    IF ( psolver == 'poisfft' )  THEN
     366       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
    366367       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
    367368    ENDIF
    … …  
    518519    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence  .OR.  &
    519520         num_acc_per_node > 0 )  THEN
    520        print*, '*** allocating diss'
    521521       ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    522522    ENDIF
    … …  
    14961496       CALL disturb_field( nzb_v_inner, tend, v )
    14971497       n_sor = nsor_ini
    1498        !$acc data copy( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, p, tric, u, v, w, weight_pres, weight_substep, tend )
     1498       !$acc data copy( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, p, tri, tric, u, v, w, weight_pres, weight_substep, tend )
    14991499       CALL pres
    15001500       !$acc end data

  • palm/trunk/SOURCE/init_pegrid.f90

    r1160 r1212  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! error message for poisfft_hybrid removed
    2323!
    2424! Former revisions:
    … …  
    232232       CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
    233233
    234     ENDIF
    235 
    236 !
    237 !-- The hybrid solver can only be used in case of a 1d-decomposition along x
    238     IF ( pdims(2) /= 1  .AND.  psolver == 'poisfft_hybrid' )  THEN
    239        message_string = 'psolver = "poisfft_hybrid" can only be' // &
    240                         '& used in case of a 1d-decomposition along x'
    241        CALL message( 'init_pegrid', 'PA0223', 1, 2, 0, 6, 0 )
    242234    ENDIF
    243235

  • palm/trunk/SOURCE/modules.f90

    r1182 r1212  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +tri
    2323!
    2424! Former revisions:
    … …  
    463463#endif
    464464
    465     REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
     465    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, tri
    466466
    467467    REAL, DIMENSION(:,:,:), ALLOCATABLE :: var_x, var_y, var_z, gamma_x,       &

  • palm/trunk/SOURCE/palm.f90

    r1182 r1212  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! +tri in copyin statement
    2323!
    2424! Former revisions:
    … …  
    244244!-- host values
    245245    !$acc  data copyin( d, diss, e, e_p, kh, km, p, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
    246     !$acc       copyin( tric, dzu, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ref_state, ug, u_init, vg, v_init, zu, zw )   &
     246    !$acc       copyin( tri, tric, dzu, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ref_state, ug, u_init, vg, v_init, zu, zw )   &
    247247    !$acc       copyin( hom, qs, qsws, qswst, rif, rif_wall, shf, ts, tswst, us, usws, uswst, vsws, vswst, z0, z0h )      &
    248248    !$acc       copyin( fxm, fxp, fym, fyp, fwxm, fwxp, fwym, fwyp, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u )       &

  • palm/trunk/SOURCE/poisfft.f90

    r1209 r1212  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! tridia routines moved to seperate module tridia_solver
    2323!
    2424! Former revisions:
    … …  
    171171    USE indices
    172172    USE transpose_indices
     173    USE tridia_solver
    173174
    174175    IMPLICIT NONE
    175176
    176177    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
    177 
    178     REAL, DIMENSION(:,:), ALLOCATABLE ::  ddzuw
    179178
    180179    PRIVATE
    … …  
    211210       CALL fft_init
    212211
    213        ALLOCATE( ddzuw(0:nz-1,3) )
    214 
    215        DO  k = 0, nz-1
    216           ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
    217           ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
    218           ddzuw(k,3) = -1.0 * &
    219                        ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
    220        ENDDO
    221 !
    222 !--    Calculate constant coefficients of the tridiagonal matrix
    223 #if ! defined ( __check )
    224        CALL maketri
    225 #endif
     212       CALL tridia_init
    226213
    227214       poisfft_initialized = .TRUE.
    … …  
    323310!--       Solve the tridiagonal equation system along z
    324311          CALL cpu_log( log_point_s(6), 'tridia', 'start' )
    325           CALL tridia( ar )
     312          CALL tridia_substi( ar )
    326313          CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
    327314
    … …  
    10911078    END SUBROUTINE ffty_tri_ffty
    10921079
    1093 
    1094     SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
    1095 
    1096 !------------------------------------------------------------------------------!
    1097 ! Solves the linear system of equations for a 1d-decomposition along x (see
    1098 ! tridia)
    1099 !
    1100 ! Attention:  when using the intel compilers older than 12.0, array tri must
    1101 !             be passed as an argument to the contained subroutines. Otherwise
    1102 !             addres faults will occur. This feature can be activated with
    1103 !             cpp-switch __intel11
    1104 !             On NEC, tri should not be passed (except for routine substi_1dd)
    1105 !             because this causes very bad performance.
    1106 !------------------------------------------------------------------------------!
    1107 
    1108        USE arrays_3d
    1109        USE control_parameters
    1110 
    1111        USE pegrid
    1112 
    1113        IMPLICIT NONE
    1114 
    1115        INTEGER ::  i, j, k, nnyh, nx, ny, omp_get_thread_num, tn
    1116 
    1117        REAL    ::  ddx2, ddy2
    1118 
    1119        REAL, DIMENSION(0:nx,1:nz)     ::  ar
    1120        REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    1121 
    1122 
    1123        nnyh = ( ny + 1 ) / 2
    1124 
    1125 !
    1126 !--    Define constant elements of the tridiagonal matrix.
    1127 !--    The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
    1128 !--    the exchanged loops create bank conflicts. The following directive
    1129 !--    prohibits loop exchange and the loops perform much better.
    1130 !       tn = omp_get_thread_num()
    1131 !       WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
    1132 !       CALL local_flush( 120+tn )
    1133 !CDIR NOLOOPCHG
    1134        DO  k = 0, nz-1
    1135           DO  i = 0,nx
    1136              tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
    1137              tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
    1138           ENDDO
    1139        ENDDO
    1140 !       WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop   thread=', omp_get_thread_num()
    1141 !       CALL local_flush( 120+tn )
    1142 
    1143        IF ( j <= nnyh )  THEN
    1144 #if defined( __intel11 )
    1145           CALL maketri_1dd( j, tri )
    1146 #else
    1147           CALL maketri_1dd( j )
    11481080#endif
    1149        ELSE
    1150 #if defined( __intel11 )
    1151           CALL maketri_1dd( ny+1-j, tri )
    1152 #else
    1153           CALL maketri_1dd( ny+1-j )
    1154 #endif
    1155        ENDIF
    1156 #if defined( __intel11 )
    1157        CALL split_1dd( tri )
    1158 #else
    1159        CALL split_1dd
    1160 #endif
    1161        CALL substi_1dd( ar, tri )
    1162 
    1163     CONTAINS
    1164 
    1165 #if defined( __intel11 )
    1166        SUBROUTINE maketri_1dd( j, tri )
    1167 #else
    1168        SUBROUTINE maketri_1dd( j )
    1169 #endif
    1170 
    1171 !------------------------------------------------------------------------------!
    1172 ! computes the i- and j-dependent component of the matrix
    1173 !------------------------------------------------------------------------------!
    1174 
    1175           USE constants
    1176 
    1177           IMPLICIT NONE
    1178 
    1179           INTEGER ::  i, j, k, nnxh
    1180           REAL    ::  a, c
    1181 
    1182           REAL, DIMENSION(0:nx) ::  l
    1183 
    1184 #if defined( __intel11 )
    1185           REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    1186 #endif
    1187 
    1188 
    1189           nnxh = ( nx + 1 ) / 2
    1190 !
    1191 !--       Provide the tridiagonal matrix for solution of the Poisson equation in
    1192 !--       Fourier space. The coefficients are computed following the method of
    1193 !--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
    1194 !--       Siano's original version by discretizing the Poisson equation,
    1195 !--       before it is Fourier-transformed
    1196           DO  i = 0, nx
    1197              IF ( i >= 0 .AND. i <= nnxh ) THEN
    1198                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    1199                                          REAL( nx+1 ) ) ) * ddx2 + &
    1200                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    1201                                          REAL( ny+1 ) ) ) * ddy2
    1202              ELSE
    1203                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    1204                                          REAL( nx+1 ) ) ) * ddx2 + &
    1205                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    1206                                          REAL( ny+1 ) ) ) * ddy2
    1207              ENDIF
    1208           ENDDO
    1209 
    1210           DO  k = 0, nz-1
    1211              DO  i = 0, nx
    1212                 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
    1213                 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
    1214                 tri(1,i,k) = a + c - l(i)
    1215              ENDDO
    1216           ENDDO
    1217           IF ( ibc_p_b == 1 )  THEN
    1218              DO  i = 0, nx
    1219                 tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
    1220              ENDDO
    1221           ENDIF
    1222           IF ( ibc_p_t == 1 )  THEN
    1223              DO  i = 0, nx
    1224                 tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
    1225              ENDDO
    1226           ENDIF
    1227 
    1228        END SUBROUTINE maketri_1dd
    1229 
    1230 
    1231 #if defined( __intel11 )
    1232        SUBROUTINE split_1dd( tri )
    1233 #else
    1234        SUBROUTINE split_1dd
    1235 #endif
    1236 
    1237 !------------------------------------------------------------------------------!
    1238 ! Splitting of the tridiagonal matrix (Thomas algorithm)
    1239 !------------------------------------------------------------------------------!
    1240 
    1241           IMPLICIT NONE
    1242 
    1243           INTEGER ::  i, k
    1244 
    1245 #if defined( __intel11 )
    1246           REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    1247 #endif
    1248 
    1249 
    1250 !
    1251 !--       Splitting
    1252           DO  i = 0, nx
    1253              tri(4,i,0) = tri(1,i,0)
    1254           ENDDO
    1255           DO  k = 1, nz-1
    1256              DO  i = 0, nx
    1257                 tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
    1258                 tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
    1259              ENDDO
    1260           ENDDO
    1261 
    1262        END SUBROUTINE split_1dd
    1263 
    1264 
    1265        SUBROUTINE substi_1dd( ar, tri )
    1266 
    1267 !------------------------------------------------------------------------------!
    1268 ! Substitution (Forward and Backward) (Thomas algorithm)
    1269 !------------------------------------------------------------------------------!
    1270 
    1271           IMPLICIT NONE
    1272 
    1273           INTEGER ::  i, k
    1274 
    1275           REAL, DIMENSION(0:nx,nz)       ::  ar
    1276           REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
    1277           REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    1278 
    1279 !
    1280 !--       Forward substitution
    1281           DO  i = 0, nx
    1282              ar1(i,0) = ar(i,1)
    1283           ENDDO
    1284           DO  k = 1, nz-1
    1285              DO  i = 0, nx
    1286                 ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
    1287              ENDDO
    1288           ENDDO
    1289 
    1290 !
    1291 !--       Backward substitution
    1292 !--       Note, the add of 1.0E-20 in the denominator is due to avoid divisions
    1293 !--       by zero appearing if the pressure bc is set to neumann at the top of
    1294 !--       the model domain.
    1295           DO  i = 0, nx
    1296              ar(i,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
    1297           ENDDO
    1298           DO  k = nz-2, 0, -1
    1299              DO  i = 0, nx
    1300                 ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
    1301                             / tri(4,i,k)
    1302              ENDDO
    1303           ENDDO
    1304 
    1305 !
    1306 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    1307 !--       The respective values of ar should be zero at all k-levels if
    1308 !--       acceleration of horizontally averaged vertical velocity is zero.
    1309           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    1310              IF ( j == 0 )  THEN
    1311                 DO  k = 1, nz
    1312                    ar(0,k) = 0.0
    1313                 ENDDO
    1314              ENDIF
    1315           ENDIF
    1316 
    1317        END SUBROUTINE substi_1dd
    1318 
    1319     END SUBROUTINE tridia_1dd
    1320 
    1321 
    1322     SUBROUTINE tridia( ar )
    1323 
    1324 !------------------------------------------------------------------------------!
    1325 ! solves the linear system of equations:
    1326 !
    1327 ! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
    1328 !   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
    1329 ! 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)
    1330 !
    1331 ! by using the Thomas algorithm
    1332 !------------------------------------------------------------------------------!
    1333 
    1334        USE arrays_3d
    1335 
    1336        IMPLICIT NONE
    1337 
    1338        INTEGER ::  i, j, k
    1339 
    1340        !$acc declare create( tri )
    1341        REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
    1342 
    1343        REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
    1344 
    1345 
    1346        CALL split( tri )
    1347        CALL substi( ar, tri )
    1348 
    1349     END SUBROUTINE tridia
    1350 
    1351 
    1352     SUBROUTINE maketri
    1353 
    1354 !------------------------------------------------------------------------------!
    1355 ! Computes the i- and j-dependent component of the matrix
    1356 !------------------------------------------------------------------------------!
    1357 
    1358           USE arrays_3d,  ONLY: tric
    1359           USE constants
    1360           USE control_parameters
    1361           USE grid_variables
    1362 
    1363           IMPLICIT NONE
    1364 
    1365           INTEGER ::  i, j, k, nnxh, nnyh
    1366 
    1367           !$acc declare create( ll )
    1368           REAL    ::  ll(nxl_z:nxr_z,nys_z:nyn_z)
    1369 
    1370 
    1371           nnxh = ( nx + 1 ) / 2
    1372           nnyh = ( ny + 1 ) / 2
    1373 
    1374 !
    1375 !--       Provide the constant coefficients of the tridiagonal matrix for solution
    1376 !--       of the Poisson equation in Fourier space.
    1377 !--       The coefficients are computed following the method of
    1378 !--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
    1379 !--       Siano's original version by discretizing the Poisson equation,
    1380 !--       before it is Fourier-transformed.
    1381 
    1382           !$acc kernels present( tric )
    1383           !$acc loop vector( 32 )
    1384           DO  j = nys_z, nyn_z
    1385              DO  i = nxl_z, nxr_z
    1386                 IF ( j >= 0  .AND.  j <= nnyh )  THEN
    1387                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
    1388                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    1389                                                   REAL( nx+1 ) ) ) / ( dx * dx ) + &
    1390                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    1391                                                   REAL( ny+1 ) ) ) / ( dy * dy )
    1392                    ELSE
    1393                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    1394                                                   REAL( nx+1 ) ) ) / ( dx * dx ) + &
    1395                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    1396                                                   REAL( ny+1 ) ) ) / ( dy * dy )
    1397                    ENDIF
    1398                 ELSE
    1399                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
    1400                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    1401                                                   REAL( nx+1 ) ) ) / ( dx * dx ) + &
    1402                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
    1403                                                   REAL( ny+1 ) ) ) / ( dy * dy )
    1404                    ELSE
    1405                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    1406                                                   REAL( nx+1 ) ) ) / ( dx * dx ) + &
    1407                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
    1408                                                   REAL( ny+1 ) ) ) / ( dy * dy )
    1409                    ENDIF
    1410                 ENDIF
    1411              ENDDO
    1412           ENDDO
    1413 
    1414           !$acc loop
    1415           DO  k = 0, nz-1
    1416              DO  j = nys_z, nyn_z
    1417                 !$acc loop vector( 32 )
    1418                 DO  i = nxl_z, nxr_z
    1419                    tric(i,j,k) = ddzuw(k,3) - ll(i,j)
    1420                 ENDDO
    1421              ENDDO
    1422           ENDDO
    1423           !$acc end kernels
    1424 
    1425           IF ( ibc_p_b == 1 )  THEN
    1426              !$acc kernels present( tric )
    1427              !$acc loop
    1428              DO  j = nys_z, nyn_z
    1429                 DO  i = nxl_z, nxr_z
    1430                    tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
    1431                 ENDDO
    1432              ENDDO
    1433              !$acc end kernels
    1434           ENDIF
    1435           IF ( ibc_p_t == 1 )  THEN
    1436              !$acc kernels present( tric )
    1437              !$acc loop
    1438              DO  j = nys_z, nyn_z
    1439                 DO  i = nxl_z, nxr_z
    1440                    tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
    1441                 ENDDO
    1442              ENDDO
    1443              !$acc end kernels
    1444           ENDIF
    1445 
    1446     END SUBROUTINE maketri
    1447 
    1448 
    1449     SUBROUTINE substi( ar, tri )
    1450 
    1451 !------------------------------------------------------------------------------!
    1452 ! Substitution (Forward and Backward) (Thomas algorithm)
    1453 !------------------------------------------------------------------------------!
    1454 
    1455           USE control_parameters
    1456 
    1457           IMPLICIT NONE
    1458 
    1459           INTEGER ::  i, j, k
    1460 
    1461           REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
    1462           REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
    1463 
    1464           !$acc declare create( ar1 )
    1465           REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1
    1466 
    1467 !
    1468 !--       Forward substitution
    1469           DO  k = 0, nz - 1
    1470              !$acc kernels present( ar, tri )
    1471              !$acc loop
    1472              DO  j = nys_z, nyn_z
    1473                 DO  i = nxl_z, nxr_z
    1474 
    1475                    IF ( k == 0 )  THEN
    1476                       ar1(i,j,k) = ar(i,j,k+1)
    1477                    ELSE
    1478                       ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
    1479                    ENDIF
    1480 
    1481                 ENDDO
    1482              ENDDO
    1483              !$acc end kernels
    1484           ENDDO
    1485 
    1486 !
    1487 !--       Backward substitution
    1488 !--       Note, the 1.0E-20 in the denominator is due to avoid divisions
    1489 !--       by zero appearing if the pressure bc is set to neumann at the top of
    1490 !--       the model domain.
    1491           DO  k = nz-1, 0, -1
    1492              !$acc kernels present( ar, tri )
    1493              !$acc loop
    1494              DO  j = nys_z, nyn_z
    1495                 DO  i = nxl_z, nxr_z
    1496 
    1497                    IF ( k == nz-1 )  THEN
    1498                       ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
    1499                    ELSE
    1500                       ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
    1501                               / tri(i,j,k,1)
    1502                    ENDIF
    1503                 ENDDO
    1504              ENDDO
    1505              !$acc end kernels
    1506           ENDDO
    1507 
    1508 !
    1509 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    1510 !--       The respective values of ar should be zero at all k-levels if
    1511 !--       acceleration of horizontally averaged vertical velocity is zero.
    1512           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    1513              IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
    1514                 !$acc kernels loop present( ar )
    1515                 DO  k = 1, nz
    1516                    ar(nxl_z,nys_z,k) = 0.0
    1517                 ENDDO
    1518              ENDIF
    1519           ENDIF
    1520 
    1521     END SUBROUTINE substi
    1522 
    1523 
    1524     SUBROUTINE split( tri )
    1525 
    1526 !------------------------------------------------------------------------------!
    1527 ! Splitting of the tridiagonal matrix (Thomas algorithm)
    1528 !------------------------------------------------------------------------------!
    1529 
    1530           USE arrays_3d,  ONLY: tric
    1531 
    1532           IMPLICIT NONE
    1533 
    1534           INTEGER ::  i, j, k
    1535 
    1536           REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
    1537 
    1538 !
    1539 !--       Splitting
    1540           !$acc kernels present( tri, tric )
    1541           !$acc loop
    1542           DO  j = nys_z, nyn_z
    1543              !$acc loop vector( 32 )
    1544              DO  i = nxl_z, nxr_z
    1545                 tri(i,j,0,1) = tric(i,j,0)
    1546              ENDDO
    1547           ENDDO
    1548           !$acc end kernels
    1549 
    1550           DO  k = 1, nz-1
    1551              !$acc kernels present( tri, tric )
    1552              !$acc loop
    1553              DO  j = nys_z, nyn_z
    1554                 !$acc loop vector( 32 )
    1555                 DO  i = nxl_z, nxr_z
    1556                    tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
    1557                    tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
    1558                 ENDDO
    1559              ENDDO
    1560              !$acc end kernels
    1561           ENDDO
    1562 
    1563     END SUBROUTINE split
    1564 
    1565 #endif
    15661081
    15671082 END MODULE poisfft_mod

  • palm/trunk/SOURCE/pres.f90

    r1118 r1212  
    1919!
    2020! Current revisions:
    21 ! -----------------
     21! ------------------
     22! call of poisfft_hybrid removed
    2223!
    2324! Former revisions:
    … …  
    131132    USE pegrid
    132133    USE poisfft_mod
    133     USE poisfft_hybrid_mod
    134134    USE statistics
    135135
    … …  
    405405!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
    406406       IF ( psolver == 'poisfft' )  THEN
    407 !
    408 !--       Solver for 2d-decomposition
     407
    409408          CALL poisfft( d, tend )
    410 
    411        ELSEIF ( psolver == 'poisfft_hybrid' )  THEN
    412 !
    413 !--       Solver for 1d-decomposition (using MPI and OpenMP).
    414 !--       The old hybrid-solver is still included here, as long as there
    415 !--       are some optimization problems in poisfft
    416           CALL poisfft_hybrid( d )
    417409
    418410       ENDIF
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |