source: palm/trunk/SOURCE/poisfft_hybrid.f90 @ 392

Last change on this file since 392 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

  • Property svn:keywords set to Id
File size: 31.2 KB
Line 
1 MODULE poisfft_hybrid_mod
2!------------------------------------------------------------------------------
3!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: poisfft_hybrid.f90 392 2009-09-24 10:39:14Z raasch $
11!
12! 274 2009-03-26 15:11:21Z heinze
13! Output of messages replaced by message handling routine.
14!
15! Feb. 2007
16! RCS Log replace by Id keyword, revision history cleaned up
17!
18! Revision 1.11  2004/04/30 12:43:14  raasch
19! Renaming of fft routines, additional argument in calls of fft_y_m
20!
21! Revision 1.2  2002/12/19 16:08:31  raasch
22! Preprocessor directive KKMP introduced (OMP does NOT work),
23! array tri will be a shared array in OpenMP loop, to get better cache
24! utilization, the i index  (x-direction) will be executed in stride
25! "istride" as outer loop and in a shorter inner loop,
26! overlapping of computation and communication realized by new routine
27! poisfft_hybrid_nodes, name of old routine poisfft_hybrid changed to
28! poisfft_hybrid_omp, STOP statement replaced by call of subroutine local_stop
29!
30!
31! Description:
32! ------------
33! Solution of the Poisson equation with a 2D spectral method. 
34! Hybrid version for parallel computers using a 1D domain decomposition,
35! realized with MPI, along x and parallelization with OPEN-MP along y
36! (routine poisfft_hybrid_omp). In a second version (poisfft_hybrid_nodes),
37! optimization is realized by overlapping of computation and communication
38! and by simultaneously executing as many communication calls as switches
39! per logical partition (LPAR) are available. This version comes into
40! effect if more than one node is used and if the environment variable
41! tasks_per_node is set in a way that it can be devided by switch_per_lpar
42! without any rest.
43!
44! WARNING: In case of OpenMP, there are problems with allocating large
45! arrays in parallel regions.
46!
47! Copyright   Klaus Ketelsen / Siegfried Raasch   May 2002
48!------------------------------------------------------------------------------!
49
50    USE fft_xy
51    USE indices
52    USE pegrid
53
54    IMPLICIT NONE
55
56    PRIVATE
57    PUBLIC poisfft_hybrid, poisfft_hybrid_ini
58
59    INTEGER, PARAMETER ::  switch_per_lpar = 2
60
61    INTEGER, SAVE ::  nxl_a, nxr_a, &      ! total   x dimension
62                      nxl_p, nxr_p, &      ! partial x dimension
63                      nys_a, nyn_a, &      ! total   y dimension
64                      nys_p, nyn_p, &      ! partial y dimension
65
66                      npe_s,        &      ! total number of PEs for solver
67                      nwords,       &      ! number of points to be exchanged
68                                           ! with MPI_ALLTOALL
69                      n_omp_threads        ! number of OpenMP threads
70
71!
72!-- Variables for multi node version (cluster version) using routine
73!-- poisfft_hybrid_nodes
74    INTEGER, SAVE :: comm_nodes,          & ! communicater nodes
75                     comm_node_all,       & ! communicater all PEs node version
76                     comm_tasks,          & ! communicater tasks
77                     me, me_node, me_task,& ! identity of this PE
78                     nodes,               & ! number of nodes
79                     tasks_per_logical_node = -1    ! default no cluster
80                     
81
82!
83!-- Public interfaces
84    INTERFACE poisfft_hybrid_ini
85       MODULE PROCEDURE poisfft_hybrid_ini
86    END INTERFACE poisfft_hybrid_ini
87
88    INTERFACE poisfft_hybrid
89       MODULE PROCEDURE poisfft_hybrid
90    END INTERFACE poisfft_hybrid
91
92!
93!-- Private interfaces
94    INTERFACE poisfft_hybrid_omp
95       MODULE PROCEDURE poisfft_hybrid_omp
96    END INTERFACE poisfft_hybrid_omp
97
98    INTERFACE poisfft_hybrid_omp_vec
99       MODULE PROCEDURE poisfft_hybrid_omp_vec
100    END INTERFACE poisfft_hybrid_omp_vec
101
102    INTERFACE poisfft_hybrid_nodes
103       MODULE PROCEDURE poisfft_hybrid_nodes
104    END INTERFACE poisfft_hybrid_nodes
105
106    INTERFACE tridia_hybrid
107       MODULE PROCEDURE tridia_hybrid
108    END INTERFACE tridia_hybrid
109
110    INTERFACE cascade
111       MODULE PROCEDURE cascade
112    END INTERFACE cascade
113
114 CONTAINS
115 
116 
117    SUBROUTINE poisfft_hybrid_ini
118
119       USE control_parameters
120       USE pegrid
121
122       IMPLICIT NONE
123
124       CHARACTER(LEN=8)      ::  cdummy
125       INTEGER               ::  idummy, istat
126       INTEGER, DIMENSION(2) ::  coords, dims
127
128       LOGICAL, DIMENSION(2) ::  period = .false., re_dims
129
130
131!
132!--    Set the internal index values for the hybrid solver
133#if defined( __parallel )
134       npe_s = pdims(1)
135#else
136       npe_s = 1
137#endif
138       nxl_a = 0
139       nxr_a = nx
140       nxl_p = 0
141       nxr_p = ( ( nx+1 ) / npe_s ) - 1
142       nys_a = nys 
143       nyn_a = nyn
144       nys_p = 0
145       nyn_p = ( ( ny+1 ) / npe_s ) - 1
146
147       nwords = ( nxr_p-nxl_p+1 ) * nz * ( nyn_p-nys_p+1 )
148
149#if defined( __KKMP )
150       CALL LOCAL_GETENV( 'OMP_NUM_THREADS', 15, cdummy, idummy )
151       READ ( cdummy, '(I8)' )  n_omp_threads
152       IF ( n_omp_threads > 1 )  THEN
153          WRITE( message_string, * ) 'Number of OpenMP threads = ', &
154                                     n_omp_threads
155          CALL message( 'poisfft_hybrid_ini', 'PA0280', 0, 0, 0, 6, 0 ) 
156       ENDIF
157#else
158       n_omp_threads = 1
159#endif
160!
161!--    Initialize the one-dimensional FFT routines
162       CALL fft_init
163
164!
165!--    Setup for multi node version (poisfft_hybrid_nodes)
166       IF ( n_omp_threads == 1  .AND.  &
167            ( host(1:4) == 'ibmh' .OR. host(1:4) == 'ibmb' ) )  THEN
168
169          IF ( tasks_per_node /= -9999 )  THEN
170!
171!--          Multi node version requires that the available number of
172!--          switches per logical partition must be an integral divisor
173!--          of the chosen number of tasks per node
174             IF ( MOD( tasks_per_node, switch_per_lpar ) == 0 )  THEN
175!
176!--             Set the switch which decides about usage of the multi node
177!--             version
178                IF ( tasks_per_node / switch_per_lpar > 1  .AND. &
179                     numprocs > tasks_per_node )  THEN
180                   tasks_per_logical_node = tasks_per_node / switch_per_lpar
181                ENDIF
182
183                IF ( tasks_per_logical_node > -1 )  THEN
184
185                   WRITE( message_string, * ) 'running optimized ',         &
186                                              'multinode version',          &
187                                              '&switch_per_lpar        = ', &
188                                              switch_per_lpar,              &
189                                              '&tasks_per_lpar         = ', &
190                                              tasks_per_node,               &
191                                              'tasks_per_logical_node = ',  &
192                                              tasks_per_logical_node
193                   CALL message( 'poisfft_hybrid_ini', 'PA0281', 0, 0, 0, 6, 0 )
194
195                ENDIF
196
197             ENDIF
198          ENDIF
199       ENDIF
200
201!
202!--    Determine sub-topologies for multi node version
203       IF ( tasks_per_logical_node >= 2 )  THEN
204
205#if defined( __parallel )
206          nodes   = ( numprocs + tasks_per_logical_node - 1 ) / &
207                    tasks_per_logical_node
208          dims(1) = nodes
209          dims(2) = tasks_per_logical_node
210
211          CALL MPI_CART_CREATE( comm2d, 2, dims, period, .FALSE., &
212                                comm_node_all, istat )
213          CALL MPI_COMM_RANK( comm_node_all, me, istat )
214
215          re_dims(1) = .TRUE.
216          re_dims(2) = .FALSE.
217          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_nodes, istat )
218          CALL MPI_COMM_RANK( comm_nodes, me_node, istat )
219
220          re_dims(1) = .FALSE.
221          re_dims(2) = .TRUE.
222          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_tasks, istat )
223          CALL MPI_COMM_RANK( comm_tasks, me_task, istat )
224
225!          write(0,*) 'who am i',myid,me,me_node,me_task,nodes,&
226!                     tasks_per_logical_node
227#else
228          message_string = 'parallel environment (MPI) required'
229          CALL message( 'poisfft_hybrid_ini', 'PA0282', 1, 2, 0, 6, 0 )
230#endif
231       ENDIF
232
233    END SUBROUTINE poisfft_hybrid_ini
234
235
236    SUBROUTINE poisfft_hybrid( ar )
237
238       USE control_parameters
239       USE interfaces
240
241       IMPLICIT NONE
242
243       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
244
245       IF ( host(1:3) == 'nec' )  THEN
246          CALL poisfft_hybrid_omp_vec( ar )
247       ELSE
248          IF ( tasks_per_logical_node == -1 )  THEN
249             CALL poisfft_hybrid_omp( ar )
250          ELSE
251             CALL poisfft_hybrid_nodes( ar )
252          ENDIF
253       ENDIF
254
255    END SUBROUTINE poisfft_hybrid
256
257
258    SUBROUTINE poisfft_hybrid_omp ( ar )
259
260       USE cpulog
261       USE interfaces
262
263       IMPLICIT NONE
264
265       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
266       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
267
268       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
269
270       REAL, DIMENSION(0:nx)                 ::  fftx_ar
271       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
272 
273       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
274
275       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
276#if defined( __KKMP )
277       INTEGER ::  omp_get_thread_num
278       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
279       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
280#else
281       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
282#endif
283
284
285       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'start' )
286
287       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
288
289!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
290!$OMP  DO
291!
292!--    Store grid points to be transformed on a 1d-array, do the fft
293!--    and sample the results on a 4d-array
294       DO  iouter = nxl_p, nxr_p, istride     ! stride loop, better cache
295          iei = MIN( iouter+istride-1, nxr_p )
296          DO  k = 1, nz
297
298             DO  i = iouter, iei
299                ii = nxl + i
300                ir = i - iouter + 1
301
302                DO  j = nys_a, nyn_a
303                   ffty_ar(j,ir) = ar(k,j,ii)
304                ENDDO
305   
306                CALL fft_y( ffty_ar(:,ir), 'forward' )
307             ENDDO
308
309             m = nys_a
310             DO  n = 1, npe_s
311                DO  j = nys_p, nyn_p
312                   DO  i = iouter, iei
313                      ir = i - iouter + 1
314                      work1(i,k,j,n) = ffty_ar(m,ir)
315                   ENDDO
316                   m = m+1
317                ENDDO
318             ENDDO
319   
320          ENDDO
321       ENDDO
322!$OMP  END PARALLEL
323
324       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
325
326#if defined( __parallel )
327       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
328
329       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
330                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
331                          comm2d, istat )
332
333       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
334#else
335       work2 = work1
336#endif
337
338       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
339
340#if defined( __KKMP )
341!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,fftx_ar,tri_ar,jthread)
342!$OMP  DO
343       DO  j = nys_p, nyn_p
344          jthread = omp_get_thread_num() + 1
345#else
346       DO  j = nys_p, nyn_p
347          jthread = 1
348#endif
349          DO  k = 1, nz
350
351             m = nxl_a
352             DO  n = 1, npe_s
353                DO  i = nxl_p, nxr_p
354                   fftx_ar(m) = work2(i,k,j,n)
355                   m = m+1
356                ENDDO
357             ENDDO
358
359             CALL fft_x( fftx_ar, 'forward' )
360
361             DO  i = nxl_a, nxr_a
362                tri_ar(i,k) = fftx_ar(i)
363             ENDDO
364
365          ENDDO
366 
367          jj = myid * (nyn_p-nys_p+1) + j
368          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
369
370          DO  k = 1, nz
371             DO  i = nxl_a, nxr_a
372                fftx_ar(i) = tri_ar (i,k)
373             ENDDO
374
375             CALL fft_x( fftx_ar, 'backward' )
376
377             m = nxl_a
378             DO  n = 1, npe_s
379                DO  i = nxl_p, nxr_p
380                   work2(i,k,j,n) = fftx_ar(m)
381                   m = m+1
382                ENDDO
383             ENDDO
384
385          ENDDO
386       ENDDO
387#if defined( __KKMP )
388!$OMP  END PARALLEL
389#endif
390
391       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
392
393#if defined( __parallel )
394       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
395       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
396
397       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
398                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
399                          comm2d, istat )
400
401       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
402#else
403       work1 = work2
404#endif
405
406       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
407
408!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
409!$OMP  DO
410       DO  iouter = nxl_p, nxr_p, istride
411          iei = MIN( iouter+istride-1, nxr_p )
412          DO  k = 1, nz
413
414             m = nys_a
415             DO  n = 1, npe_s
416                DO  j = nys_p, nyn_p
417                   DO  i = iouter, iei
418                      ir = i - iouter + 1
419                      ffty_ar(m,ir) = work1 (i,k,j,n)
420                   ENDDO
421                   m = m+1
422                ENDDO
423             ENDDO
424
425             DO  i = iouter, iei
426                ii = nxl + i
427                ir = i - iouter + 1
428                CALL fft_y( ffty_ar(:,ir), 'backward' )
429
430                DO  j = nys_a, nyn_a
431                   ar(k,j,ii) = ffty_ar(j,ir)
432                ENDDO
433             ENDDO
434
435          ENDDO
436       ENDDO
437!$OMP  END PARALLEL
438
439       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
440
441       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'stop' )
442
443#if defined( __KKMP )
444       DEALLOCATE( tri )
445#endif
446
447    END SUBROUTINE poisfft_hybrid_omp
448
449
450    SUBROUTINE poisfft_hybrid_omp_vec ( ar )
451
452       USE cpulog
453       USE interfaces
454
455       IMPLICIT NONE
456
457       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
458       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
459
460       REAL, DIMENSION(0:nx,nz)               ::  tri_ar
461
462       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr)  ::  ar
463
464       REAL, DIMENSION(0:ny+3,nz,nxl_p:nxr_p) ::  ffty_ar3
465       REAL, DIMENSION(0:nx+3,nz,nys_p:nyn_p) ::  fftx_ar3
466 
467       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
468#if defined( __KKMP )
469       INTEGER ::  omp_get_thread_num
470       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
471       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
472#else
473       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
474#endif
475
476
477       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'start' )
478
479       CALL cpu_log( log_point_s(7), 'fft_y_m', 'start' )
480
481!$OMP  PARALLEL PRIVATE (i,j,k,m,n)
482!$OMP  DO
483!
484!--    Store grid points to be transformed on a 1d-array, do the fft
485!--    and sample the results on a 4d-array
486       DO  i = nxl_p, nxr_p
487
488          DO  j = nys_a, nyn_a
489             DO  k = 1, nz
490                ffty_ar3(j,k,i) = ar(k,j,i+nxl)
491             ENDDO
492          ENDDO
493
494          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'forward' )
495       ENDDO
496
497!$OMP  DO
498       DO  k = 1, nz
499          m = nys_a
500          DO  n = 1, npe_s
501             DO  j = nys_p, nyn_p
502                 DO  i = nxl_p, nxr_p
503                     work1(i,k,j,n) = ffty_ar3(m,k,i)
504                 ENDDO
505                 m = m+1
506             ENDDO
507          ENDDO
508       ENDDO
509!$OMP  END PARALLEL
510
511       CALL cpu_log( log_point_s(7), 'fft_y_m', 'pause' )
512
513#if defined( __parallel )
514       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
515       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
516                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
517                          comm2d, istat )
518       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
519#else
520       work2 = work1
521#endif
522
523       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'start' )
524
525#if defined( __KKMP )
526!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,tri_ar,jthread)
527!$OMP  DO
528       DO  j = nys_p, nyn_p
529          jthread = omp_get_thread_num() + 1
530#else
531       DO  j = nys_p, nyn_p
532          jthread = 1
533#endif
534          DO  k = 1, nz
535
536             m = nxl_a
537             DO  n = 1, npe_s
538                DO  i = nxl_p, nxr_p
539                   fftx_ar3(m,k,j) = work2(i,k,j,n)
540                   m = m+1
541                ENDDO
542             ENDDO
543          ENDDO
544
545          CALL fft_x_m( fftx_ar3(:,:,j), 'forward' )
546
547          DO  k = 1, nz
548             DO  i = nxl_a, nxr_a
549                tri_ar(i,k) = fftx_ar3(i,k,j)
550             ENDDO
551          ENDDO
552 
553          jj = myid * (nyn_p-nys_p+1) + j
554          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
555
556          DO  k = 1, nz
557             DO  i = nxl_a, nxr_a
558                fftx_ar3(i,k,j) = tri_ar (i,k)
559             ENDDO
560          ENDDO
561
562          CALL fft_x_m( fftx_ar3(:,:,j), 'backward' )
563
564          DO  k = 1, nz
565             m = nxl_a
566             DO  n = 1, npe_s
567                DO  i = nxl_p, nxr_p
568                   work2(i,k,j,n) = fftx_ar3(m,k,j)
569                   m = m+1
570                ENDDO
571             ENDDO
572          ENDDO
573
574       ENDDO
575#if defined( __KKMP )
576!$OMP  END PARALLEL
577#endif
578
579       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'stop' )
580
581#if defined( __parallel )
582       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
583       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
584       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
585                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
586                          comm2d, istat )
587       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
588#else
589       work1 = work2
590#endif
591
592       CALL cpu_log( log_point_s(7), 'fft_y_m', 'continue' )
593
594!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n)
595!$OMP  DO
596       DO  k = 1, nz
597          m = nys_a
598          DO  n = 1, npe_s
599             DO  j = nys_p, nyn_p
600                 DO  i = nxl_p, nxr_p
601                     ffty_ar3(m,k,i) = work1(i,k,j,n)
602                 ENDDO
603                 m = m+1
604             ENDDO
605          ENDDO
606       ENDDO 
607
608!$OMP  DO
609       DO  i = nxl_p, nxr_p
610          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'backward' )
611          DO  j = nys_a, nyn_a
612             DO  k = 1, nz
613                ar(k,j,i+nxl) = ffty_ar3(j,k,i)
614             ENDDO
615          ENDDO
616       ENDDO
617!$OMP  END PARALLEL
618
619       CALL cpu_log( log_point_s(7), 'fft_y_m', 'stop' )
620
621       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'stop' )
622
623#if defined( __KKMP )
624       DEALLOCATE( tri )
625#endif
626
627    END SUBROUTINE poisfft_hybrid_omp_vec
628
629
630    SUBROUTINE poisfft_hybrid_nodes ( ar )
631
632       USE cpulog
633       USE interfaces
634
635       IMPLICIT NONE
636
637       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
638       INTEGER            ::  i, iei, ii, iouter, ir, istat, j, jj, k, m, &
639                              n, nn, nt, nw1, nw2
640
641       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
642
643       REAL, DIMENSION(0:nx)                 ::  fftx_ar
644       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
645 
646       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
647
648       REAL, DIMENSION(nxl_p:nxr_p,nz,tasks_per_logical_node, &
649                          nodes,nys_p:nyn_p) ::  work1,work2
650       REAL, DIMENSION(5,0:nx,0:nz-1)        ::  tri
651
652
653       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'start' )
654
655       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
656
657!
658!--    Store grid points to be transformed on a 1d-array, do the fft
659!--    and sample the results on a 4d-array
660       DO  iouter = nxl_p, nxr_p, istride  ! stride loop, better cache
661          iei = MIN( iouter+istride-1, nxr_p )
662          DO  k = 1, nz
663
664             DO  i = iouter, iei
665                ii = nxl + i
666                ir = i - iouter + 1
667
668                DO  j = nys_a, nyn_a
669                   ffty_ar(j,ir) = ar(k,j,ii)
670                ENDDO
671   
672                CALL fft_y( ffty_ar(:,ir), 'forward' )
673             ENDDO
674
675             m = nys_a
676             DO  nn = 1, nodes
677                DO  nt = 1, tasks_per_logical_node
678                   DO  j = nys_p, nyn_p
679                      DO  i = iouter, iei
680                         ir = i - iouter + 1
681                         work1(i,k,nt,nn,j) = ffty_ar(m,ir)
682                      ENDDO
683                      m = m+1
684                   ENDDO
685                ENDDO
686             ENDDO
687   
688          ENDDO
689       ENDDO
690
691       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
692
693       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' )
694       nw1 = SIZE( work1, 1 ) * SIZE( work1, 2 )
695       DO  nn = 1, nodes
696           DO  j = nys_p, nyn_p
697#if defined( __parallel )
698              CALL MPI_ALLTOALL( work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
699                                 work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
700                                                     comm_tasks, istat )
701#endif
702           ENDDO
703       ENDDO
704       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
705
706
707       DO  j = nys_p, nyn_p
708       
709          CALL cascade( 1, j, nys_p, nyn_p )
710          nw2 = nw1 * SIZE( work1, 3 )
711          CALL cpu_log( log_point_s(37), 'alltoall_node', 'start' )
712#if defined( __parallel )
713          CALL MPI_ALLTOALL( work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
714                             work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
715                                                  comm_nodes, istat )
716#endif
717          CALL cpu_log( log_point_s(37), 'alltoall_node', 'pause' )
718          CALL cascade( 2, j, nys_p, nyn_p )
719
720          CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
721          DO  k = 1, nz
722
723             m = nxl_a
724             DO  nn = 1, nodes
725                DO  nt = 1, tasks_per_logical_node
726                   DO  i = nxl_p, nxr_p
727                      fftx_ar(m) = work1(i,k,nt,nn,j)
728                      m = m+1
729                   ENDDO
730                ENDDO
731             ENDDO
732
733             CALL fft_x( fftx_ar, 'forward' )
734
735             DO  i = nxl_a, nxr_a
736                tri_ar(i,k) = fftx_ar(i)
737             ENDDO
738
739          ENDDO
740 
741          jj = myid * (nyn_p-nys_p+1) + j
742          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:) )
743
744          DO  k = 1, nz
745             DO  i = nxl_a, nxr_a
746                fftx_ar(i) = tri_ar(i,k)
747             ENDDO
748
749             CALL fft_x( fftx_ar, 'backward' )
750
751             m = nxl_a
752             DO  nn = 1, nodes
753                DO  nt = 1, tasks_per_logical_node
754                   DO  i = nxl_p, nxr_p
755                      work1(i,k,nt,nn,j) = fftx_ar(m)
756                      m = m+1
757                   ENDDO
758                ENDDO
759             ENDDO
760          ENDDO
761
762          CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
763          nw2 = nw1 * SIZE( work1, 3 )
764          CALL cpu_log( log_point_s(37), 'alltoall_node', 'continue' )
765#if defined( __parallel )
766          CALL MPI_ALLTOALL( work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
767                             work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
768                                                  comm_nodes, istat )
769#endif
770          CALL cpu_log( log_point_s(37), 'alltoall_node', 'stop' )
771
772       ENDDO
773
774       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' ) 
775       DO  nn = 1, nodes
776           DO  j = nys_p, nyn_p
777#if defined( __parallel )
778              CALL MPI_ALLTOALL( work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
779                                 work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
780                                                     comm_tasks, istat )
781#endif
782           ENDDO
783       ENDDO
784       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
785
786       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
787
788       DO  iouter = nxl_p, nxr_p, istride
789          iei = MIN( iouter+istride-1, nxr_p )
790          DO  k = 1, nz
791
792             m = nys_a
793             DO  nn = 1, nodes
794                DO  nt = 1, tasks_per_logical_node
795                   DO  j = nys_p, nyn_p
796                      DO  i = iouter, iei
797                         ir = i - iouter + 1
798                         ffty_ar(m,ir) = work1(i,k,nt,nn,j)
799                      ENDDO
800                      m = m+1
801                   ENDDO
802                ENDDO
803             ENDDO
804
805             DO  i = iouter, iei
806                ii = nxl + i
807                ir = i - iouter + 1
808                CALL fft_y( ffty_ar(:,ir), 'backward' )
809
810                DO  j = nys_a, nyn_a
811                   ar(k,j,ii) = ffty_ar(j,ir)
812                ENDDO
813             ENDDO
814
815          ENDDO
816       ENDDO
817
818       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
819
820       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'stop' )
821
822    END SUBROUTINE poisfft_hybrid_nodes
823
824
825
826    SUBROUTINE tridia_hybrid( j, ar, tri )
827
828       USE arrays_3d
829       USE control_parameters
830       USE grid_variables
831
832       IMPLICIT NONE
833
834       INTEGER ::  i, j, k, nnyh
835       REAL, DIMENSION(0:nx,nz)       ::  ar
836       REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
837       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
838
839       nnyh = (ny+1) / 2
840
841       tri = 0.0
842!
843!--    Define constant elements of the tridiagonal matrix.
844       DO  k = 0, nz-1
845          DO  i = 0,nx
846             tri(2,i,k) = ddzu(k+1) * ddzw(k+1)
847             tri(3,i,k) = ddzu(k+2) * ddzw(k+1)
848          ENDDO
849       ENDDO
850
851       IF ( j <= nnyh )  THEN
852          CALL maketri_hybrid( j )
853       ELSE
854          CALL maketri_hybrid( ny+1-j)
855       ENDIF
856       CALL zerleg_hybrid
857       CALL substi_hybrid( ar, tri )
858
859    CONTAINS
860
861       SUBROUTINE maketri_hybrid( j )
862
863!----------------------------------------------------------------------!
864!                         maketri                                      !
865!                                                                      !
866!       computes the i- and j-dependent component of the matrix        !
867!----------------------------------------------------------------------!
868
869          USE constants
870
871          IMPLICIT NONE
872
873          INTEGER ::  i, j, k, nnxh
874          REAL    ::  a, c
875
876          REAL, DIMENSION(0:nx) ::  l
877
878
879          nnxh = (nx+1) / 2
880!
881!--       Provide the tridiagonal matrix for solution of the Poisson equation
882!--       in Fourier space. The coefficients are computed following the method
883!--       of Schmidt et al. (DFVLR-Mitteilung 84-15) --> departs from Stephan
884!--       Siano's original version.
885          DO  i = 0,nx
886             IF ( i >= 0 .AND. i < nnxh ) THEN
887                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
888                                        FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
889                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
890                                        FLOAT( ny+1 ) ) ) / ( dy * dy )
891             ELSEIF ( i == nnxh ) THEN
892                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
893                                         FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
894                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
895                                         FLOAT(ny+1) ) ) / ( dy * dy )
896             ELSE
897                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
898                                         FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
899                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
900                                         FLOAT( ny+1 ) ) ) / ( dy * dy )
901             ENDIF
902          ENDDO
903
904          DO  k = 0,nz-1
905             DO  i = 0, nx
906                a = -1.0 * ddzu(k+2) * ddzw(k+1)
907                c = -1.0 * ddzu(k+1) * ddzw(k+1)
908                tri(1,i,k) = a + c - l(i)
909             ENDDO
910          ENDDO
911          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
912             DO  i = 0,nx
913                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
914             ENDDO
915          ENDIF
916          IF ( ibc_p_t == 1 )  THEN
917             DO  i = 0,nx
918                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
919             ENDDO
920          ENDIF
921
922       END SUBROUTINE maketri_hybrid
923
924
925       SUBROUTINE zerleg_hybrid
926
927!----------------------------------------------------------------------!
928!                          zerleg                                      !
929!                                                                      !
930!       Splitting of the tridiagonal matrix (Thomas algorithm)         !
931!----------------------------------------------------------------------!
932
933          USE indices
934
935          IMPLICIT NONE
936
937          INTEGER ::  i, k
938
939!
940!--       Splitting
941          DO  i = 0, nx
942             tri(4,i,0) = tri(1,i,0)
943          ENDDO
944          DO  k = 1, nz-1
945             DO  i = 0,nx
946                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
947                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
948             ENDDO
949          ENDDO
950
951       END SUBROUTINE zerleg_hybrid
952
953       SUBROUTINE substi_hybrid( ar, tri )
954
955!----------------------------------------------------------------------!
956!                          substi                                      !
957!                                                                      !
958!       Substitution (Forward and Backward) (Thomas algorithm)         !
959!----------------------------------------------------------------------!
960
961          IMPLICIT NONE
962
963          INTEGER ::  i, j, k
964          REAL, DIMENSION(0:nx,nz)       ::  ar
965          REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
966          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
967
968!
969!--       Forward substitution
970          DO  i = 0, nx
971             ar1(i,0) = ar(i,1)
972          ENDDO
973          DO  k = 1, nz - 1
974             DO  i = 0,nx
975                ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
976             ENDDO
977          ENDDO
978
979!
980!--       Backward substitution
981          DO  i = 0,nx
982             ar(i,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
983          ENDDO
984          DO  k = nz-2, 0, -1
985             DO  i = 0,nx
986                ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
987                            / tri(4,i,k)
988             ENDDO
989          ENDDO
990
991       END SUBROUTINE substi_hybrid
992
993    END SUBROUTINE tridia_hybrid
994
995
996    SUBROUTINE cascade( loca, j, nys_p, nyn_p )
997
998       USE cpulog
999
1000       IMPLICIT NONE
1001
1002       INTEGER ::  ier, j, loca, nyn_p, nys_p, req, reqa(1)
1003       INTEGER, SAVE ::  tag = 10
1004#if defined( __parallel )
1005       INTEGER, DIMENSION(MPI_STATUS_SIZE) :: stat
1006#endif
1007
1008       REAL ::  buf, buf1
1009
1010
1011       buf  = 1.0
1012       buf1 = 1.1
1013       IF ( me_node == 0 )  THEN  ! first node only
1014
1015          SELECT CASE ( loca )
1016
1017             CASE ( 1 )  ! before alltoall
1018
1019                IF( me_task > 0 )  THEN  ! first task does not wait
1020#if defined( __parallel )
1021                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task-1, 0, &
1022                                      buf1, 1, MPI_REAL, me_task-1, 0, &
1023                                                 comm_tasks, stat,ierr )
1024#endif
1025                ELSEIF ( j > nys_p )  THEN
1026                   req = 0
1027                   tag = MOD( tag-10, 10 ) + 10
1028#if defined( __parallel )
1029                   CALL MPI_IRECV( buf, 1, MPI_REAL, tasks_per_logical_node-1,&
1030                                   tag, comm_tasks, req, ierr )
1031                   reqa = req
1032                   CALL MPI_WAITALL( 1, reqa, stat, ierr )
1033#endif
1034                ENDIF
1035
1036             CASE ( 2 )  ! after alltoall
1037
1038                IF ( me_task < tasks_per_logical_node-1 )  THEN  ! last task
1039#if defined( __parallel )
1040                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task+1, 0, &
1041                                      buf1, 1, MPI_REAL, me_task+1, 0, &
1042                                                 comm_tasks, stat, ierr)
1043#endif
1044                ELSEIF ( j < nyn_p )  THEN
1045                   req = 0
1046                   tag = MOD( tag-10, 10 ) + 10
1047#if defined( __parallel )
1048                   CALL MPI_ISEND( buf, 1, MPI_REAL, 0, tag, comm_tasks, req, &
1049                                   ierr )
1050#endif
1051                ENDIF
1052
1053          END SELECT
1054
1055       ENDIF
1056
1057    END SUBROUTINE cascade
1058
1059 END MODULE poisfft_hybrid_mod
Note: See TracBrowser for help on using the repository browser.