Home

Context Navigation

source: palm/trunk/SOURCE/init_pegrid.f90 @ 29

Last change on this file since 29 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 24.9 KB

Line
1	SUBROUTINE init_pegrid
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: init_pegrid.f90 4 2007-02-13 11:33:16Z raasch $
11	! RCS Log replace by Id keyword, revision history cleaned up
12	!
13	! Revision 1.28 2006/04/26 13:23:32 raasch
14	! lcmuk does not understand the !$ comment so a cpp-directive is required
15	!
16	! Revision 1.1 1997/07/24 11:15:09 raasch
17	! Initial revision
18	!
19	!
20	! Description:
21	! ------------
22	! Determination of the virtual processor topology (if not prescribed by the
23	! user)and computation of the grid point number and array bounds of the local
24	! domains.
25	!------------------------------------------------------------------------------!
26
27	USE control_parameters
28	USE fft_xy
29	USE indices
30	USE pegrid
31	USE poisfft_mod
32	USE poisfft_hybrid_mod
33	USE statistics
34	USE transpose_indices
35
36
37	IMPLICIT NONE
38
39	INTEGER :: gathered_size, i, ind(5), j, k, maximum_grid_level_l, &
40	mg_switch_to_pe0_level_l, mg_levels_x, mg_levels_y, &
41	mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, nnz_y, &
42	numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l, nzb_l, &
43	nzt_l, omp_get_num_threads, subdomain_size
44
45	INTEGER, DIMENSION(:), ALLOCATABLE :: ind_all, nxlf, nxrf, nynf, nysf
46
47	LOGICAL :: found
48
49	!
50	!-- Get the number of OpenMP threads
51	!$OMP PARALLEL
52	#if defined( __lcmuk )
53	threads_per_task = omp_get_num_threads()
54	#else
55	!$ threads_per_task = omp_get_num_threads()
56	#endif
57	!$OMP END PARALLEL
58
59
60	#if defined( __parallel )
61	!
62	!-- Determine the processor topology or check it, if prescribed by the user
63	IF ( npex == -1 .AND. npey == -1 ) THEN
64
65	!
66	!-- Automatic determination of the topology
67	!-- The default on SMP- and cluster-hosts is a 1d-decomposition along x
68	#if defined( __lcmuk )
69	host = 'lcmuk'
70	#endif
71	IF ( host(1:3) == 'ibm' .OR. host(1:3) == 'nec' .OR. &
72	host(1:2) == 'lc' .OR. host(1:3) == 'dec' ) THEN
73
74	pdims(1) = numprocs
75	pdims(2) = 1
76
77	ELSE
78
79	numproc_sqr = SQRT( REAL( numprocs ) )
80	pdims(1) = MAX( numproc_sqr , 1 )
81	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
82	pdims(1) = pdims(1) - 1
83	ENDDO
84	pdims(2) = numprocs / pdims(1)
85
86	ENDIF
87
88	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
89
90	!
91	!-- Prescribed by user. Number of processors on the prescribed topology
92	!-- must be equal to the number of PEs available to the job
93	IF ( ( npex * npey ) /= numprocs ) THEN
94	PRINT*, '+++ init_pegrid:'
95	PRINT, ' number of PEs of the prescribed topology (', npexnpey, &
96	') does not match the number of PEs available to the ', &
97	'job (', numprocs, ')'
98	CALL local_stop
99	ENDIF
100	pdims(1) = npex
101	pdims(2) = npey
102
103	ELSE
104	!
105	!-- If the processor topology is prescribed by the user, the number of
106	!-- PEs must be given in both directions
107	PRINT*, '+++ init_pegrid:'
108	PRINT*, ' if the processor topology is prescribed by the user, ', &
109	'both values of "npex" and "npey" must be given in the ', &
110	'NAMELIST-parameter file'
111	CALL local_stop
112
113	ENDIF
114
115	!
116	!-- The hybrid solver can only be used in case of a 1d-decomposition along x
117	IF ( pdims(2) /= 1 .AND. psolver == 'poisfft_hybrid' ) THEN
118	IF ( myid == 0 ) THEN
119	PRINT, '** init_pegrid: psolver = "poisfft_hybrid" can only be'
120	PRINT*, ' used in case of a 1d-decomposition along x'
121	ENDIF
122	ENDIF
123
124	!
125	!-- If necessary, set horizontal boundary conditions to non-cyclic
126	IF ( bc_lr /= 'cyclic' ) cyclic(1) = .FALSE.
127	IF ( bc_ns /= 'cyclic' ) cyclic(2) = .FALSE.
128
129	!
130	!-- Create the virtual processor grid
131	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
132	comm2d, ierr )
133	CALL MPI_COMM_RANK( comm2d, myid, ierr )
134	WRITE (myid_char,'(''_'',I4.4)') myid
135
136	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
137	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
138	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
139
140	!
141	!-- Determine sub-topologies for transpositions
142	!-- Transposition from z to x:
143	remain_dims(1) = .TRUE.
144	remain_dims(2) = .FALSE.
145	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
146	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
147	!
148	!-- Transposition from x to y
149	remain_dims(1) = .FALSE.
150	remain_dims(2) = .TRUE.
151	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
152	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
153
154
155	!
156	!-- Find a grid (used for array d) which will match the transposition demands
157	IF ( grid_matching == 'strict' ) THEN
158
159	nxa = nx; nya = ny; nza = nz
160
161	ELSE
162
163	found = .FALSE.
164	xn: DO nxa = nx, 2*nx
165	!
166	!-- Meet conditions for nx
167	IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
168	MOD( nxa+1, pdims(2) ) /= 0 ) CYCLE xn
169
170	yn: DO nya = ny, 2*ny
171	!
172	!-- Meet conditions for ny
173	IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
174	MOD( nya+1, pdims(1) ) /= 0 ) CYCLE yn
175
176
177	zn: DO nza = nz, 2*nz
178	!
179	!-- Meet conditions for nz
180	IF ( ( MOD( nza, pdims(1) ) /= 0 .AND. pdims(1) /= 1 .AND. &
181	pdims(2) /= 1 ) .OR. &
182	( MOD( nza, pdims(2) ) /= 0 .AND. dt_dosp /= 9999999.9 &
183	) ) THEN
184	CYCLE zn
185	ELSE
186	found = .TRUE.
187	EXIT xn
188	ENDIF
189
190	ENDDO zn
191
192	ENDDO yn
193
194	ENDDO xn
195
196	IF ( .NOT. found ) THEN
197	IF ( myid == 0 ) THEN
198	PRINT*,'+++ init_pegrid: no matching grid for transpositions found'
199	ENDIF
200	CALL local_stop
201	ENDIF
202
203	ENDIF
204
205	!
206	!-- Calculate array bounds in x-direction for every PE.
207	!-- The last PE along x may get less grid points than the others
208	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
209	nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
210
211	IF ( MOD( nxa+1 , pdims(1) ) /= 0 ) THEN
212	IF ( myid == 0 ) THEN
213	PRINT*,'+++ x-direction: gridpoint number (',nx+1,') is not an'
214	PRINT*,' integral divisor of the number of proces', &
215	&'sors (', pdims(1),')'
216	ENDIF
217	CALL local_stop
218	ELSE
219	nnx = ( nxa + 1 ) / pdims(1)
220	IF ( nnx*pdims(1) - ( nx + 1) > nnx ) THEN
221	IF ( myid == 0 ) THEN
222	PRINT*,'+++ x-direction: nx does not match the requirements ', &
223	'given by the number of PEs'
224	PRINT*,' used'
225	PRINT,' please use nx = ', nx - ( pdims(1) - ( nnxpdims(1) &
226	- ( nx + 1 ) ) ), ' instead of nx =', nx
227	ENDIF
228	CALL local_stop
229	ENDIF
230	ENDIF
231
232	!
233	!-- Left and right array bounds, number of gridpoints
234	DO i = 0, pdims(1)-1
235	nxlf(i) = i * nnx
236	nxrf(i) = ( i + 1 ) * nnx - 1
237	nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
238	ENDDO
239
240	!
241	!-- Calculate array bounds in y-direction for every PE.
242	IF ( MOD( nya+1 , pdims(2) ) /= 0 ) THEN
243	IF ( myid == 0 ) THEN
244	PRINT*,'+++ y-direction: gridpoint number (',ny+1,') is not an'
245	PRINT*,' integral divisor of the number of proces', &
246	&'sors (', pdims(2),')'
247	ENDIF
248	CALL local_stop
249	ELSE
250	nny = ( nya + 1 ) / pdims(2)
251	IF ( nny*pdims(2) - ( ny + 1) > nny ) THEN
252	IF ( myid == 0 ) THEN
253	PRINT*,'+++ x-direction: nx does not match the requirements ', &
254	'given by the number of PEs'
255	PRINT*,' used'
256	PRINT,' please use nx = ', nx - ( pdims(1) - ( nnxpdims(1) &
257	- ( nx + 1 ) ) ), ' instead of nx =', nx
258	ENDIF
259	CALL local_stop
260	ENDIF
261	ENDIF
262
263	!
264	!-- South and north array bounds
265	DO j = 0, pdims(2)-1
266	nysf(j) = j * nny
267	nynf(j) = ( j + 1 ) * nny - 1
268	nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
269	ENDDO
270
271	!
272	!-- Local array bounds of the respective PEs
273	nxl = nxlf(pcoord(1))
274	nxra = nxrf(pcoord(1))
275	nxr = MIN( nx, nxra )
276	nys = nysf(pcoord(2))
277	nyna = nynf(pcoord(2))
278	nyn = MIN( ny, nyna )
279	nzb = 0
280	nzta = nza
281	nzt = MIN( nz, nzta )
282	nnz = nza
283
284	!
285	!-- Calculate array bounds and gridpoint numbers for the transposed arrays
286	!-- (needed in the pressure solver)
287	!-- For the transposed arrays, cyclic boundaries as well as top and bottom
288	!-- boundaries are omitted, because they are obstructive to the transposition
289
290	!
291	!-- 1. transposition z --> x
292	!-- This transposition is not neccessary in case of a 1d-decomposition along x,
293	!-- except that the uptream-spline method is switched on
294	IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. &
295	scalar_advec == 'ups-scheme' ) THEN
296
297	IF ( pdims(2) == 1 .AND. ( momentum_advec == 'ups-scheme' .OR. &
298	scalar_advec == 'ups-scheme' ) ) THEN
299	IF ( myid == 0 ) THEN
300	PRINT*,'+++ WARNING: init_pegrid: 1d-decomposition along x ', &
301	&'chosen but nz restrictions may occur'
302	PRINT*,' since ups-scheme is activated'
303	ENDIF
304	ENDIF
305	nys_x = nys
306	nyn_xa = nyna
307	nyn_x = nyn
308	nny_x = nny
309	IF ( MOD( nza , pdims(1) ) /= 0 ) THEN
310	IF ( myid == 0 ) THEN
311	PRINT*,'+++ transposition z --> x:'
312	PRINT*,' nz=',nz,' is not an integral divisior of pdims(1)=', &
313	&pdims(1)
314	ENDIF
315	CALL local_stop
316	ENDIF
317	nnz_x = nza / pdims(1)
318	nzb_x = 1 + myidx * nnz_x
319	nzt_xa = ( myidx + 1 ) * nnz_x
320	nzt_x = MIN( nzt, nzt_xa )
321
322	sendrecvcount_zx = nnx * nny * nnz_x
323
324	ENDIF
325
326	!
327	!-- 2. transposition x --> y
328	nnz_y = nnz_x
329	nzb_y = nzb_x
330	nzt_ya = nzt_xa
331	nzt_y = nzt_x
332	IF ( MOD( nxa+1 , pdims(2) ) /= 0 ) THEN
333	IF ( myid == 0 ) THEN
334	PRINT*,'+++ transposition x --> y:'
335	PRINT*,' nx+1=',nx+1,' is not an integral divisor of ',&
336	&'pdims(2)=',pdims(2)
337	ENDIF
338	CALL local_stop
339	ENDIF
340	nnx_y = (nxa+1) / pdims(2)
341	nxl_y = myidy * nnx_y
342	nxr_ya = ( myidy + 1 ) * nnx_y - 1
343	nxr_y = MIN( nx, nxr_ya )
344
345	sendrecvcount_xy = nnx_y * nny_x * nnz_y
346
347	!
348	!-- 3. transposition y --> z (ELSE: x --> y in case of 1D-decomposition
349	!-- along x)
350	IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. &
351	scalar_advec == 'ups-scheme' ) THEN
352	!
353	!-- y --> z
354	!-- This transposition is not neccessary in case of a 1d-decomposition
355	!-- along x, except that the uptream-spline method is switched on
356	nnx_z = nnx_y
357	nxl_z = nxl_y
358	nxr_za = nxr_ya
359	nxr_z = nxr_y
360	IF ( MOD( nya+1 , pdims(1) ) /= 0 ) THEN
361	IF ( myid == 0 ) THEN
362	PRINT*,'+++ Transposition y --> z:'
363	PRINT*,' ny+1=',ny+1,' is not an integral divisor of ',&
364	&'pdims(1)=',pdims(1)
365	ENDIF
366	CALL local_stop
367	ENDIF
368	nny_z = (nya+1) / pdims(1)
369	nys_z = myidx * nny_z
370	nyn_za = ( myidx + 1 ) * nny_z - 1
371	nyn_z = MIN( ny, nyn_za )
372
373	sendrecvcount_yz = nnx_y * nny_z * nnz_y
374
375	ELSE
376	!
377	!-- x --> y. This condition must be fulfilled for a 1D-decomposition along x
378	IF ( MOD( nya+1 , pdims(1) ) /= 0 ) THEN
379	IF ( myid == 0 ) THEN
380	PRINT*,'+++ Transposition x --> y:'
381	PRINT*,' ny+1=',ny+1,' is not an integral divisor of ',&
382	&'pdims(1)=',pdims(1)
383	ENDIF
384	CALL local_stop
385	ENDIF
386
387	ENDIF
388
389	!
390	!-- Indices for direct transpositions z --> y (used for calculating spectra)
391	IF ( dt_dosp /= 9999999.9 ) THEN
392	IF ( MOD( nza, pdims(2) ) /= 0 ) THEN
393	IF ( myid == 0 ) THEN
394	PRINT*,'+++ Direct transposition z --> y (needed for spectra):'
395	PRINT*,' nz=',nz,' is not an integral divisor of ',&
396	&'pdims(2)=',pdims(2)
397	ENDIF
398	CALL local_stop
399	ELSE
400	nxl_yd = nxl
401	nxr_yda = nxra
402	nxr_yd = nxr
403	nzb_yd = 1 + myidy * ( nza / pdims(2) )
404	nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
405	nzt_yd = MIN( nzt, nzt_yda )
406
407	sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
408	ENDIF
409	ENDIF
410
411	!
412	!-- Indices for direct transpositions y --> x (they are only possible in case
413	!-- of a 1d-decomposition along x)
414	IF ( pdims(2) == 1 ) THEN
415	nny_x = nny / pdims(1)
416	nys_x = myid * nny_x
417	nyn_xa = ( myid + 1 ) * nny_x - 1
418	nyn_x = MIN( ny, nyn_xa )
419	nzb_x = 1
420	nzt_xa = nza
421	nzt_x = nz
422	sendrecvcount_xy = nnx * nny_x * nza
423	ENDIF
424
425	!
426	!-- Indices for direct transpositions x --> y (they are only possible in case
427	!-- of a 1d-decomposition along y)
428	IF ( pdims(1) == 1 ) THEN
429	nnx_y = nnx / pdims(2)
430	nxl_y = myid * nnx_y
431	nxr_ya = ( myid + 1 ) * nnx_y - 1
432	nxr_y = MIN( nx, nxr_ya )
433	nzb_y = 1
434	nzt_ya = nza
435	nzt_y = nz
436	sendrecvcount_xy = nnx_y * nny * nza
437	ENDIF
438
439	!
440	!-- Arrays for storing the array bounds are needed any more
441	DEALLOCATE( nxlf , nxrf , nynf , nysf )
442
443	#if defined( __print )
444	!
445	!-- Control output
446	IF ( myid == 0 ) THEN
447	PRINT, ' processor topology *'
448	PRINT*, ' '
449	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',&
450	&' nys: nyn'
451	PRINT*, '------------------------------------------------------------',&
452	&'-----------'
453	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
454	myidx, myidy, nxl, nxr, nys, nyn
455	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
456	2(2X,I4,':',I4))
457
458	!
459	!-- Recieve data from the other PEs
460	DO i = 1,numprocs-1
461	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
462	ierr )
463	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
464	ENDDO
465	ELSE
466
467	!
468	!-- Send data to PE0
469	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
470	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
471	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
472	ibuf(12) = nyn
473	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
474	ENDIF
475	#endif
476
477	#else
478
479	!
480	!-- Array bounds when running on a single PE (respectively a non-parallel
481	!-- machine)
482	nxl = 0
483	nxr = nx
484	nxra = nx
485	nnx = nxr - nxl + 1
486	nys = 0
487	nyn = ny
488	nyna = ny
489	nny = nyn - nys + 1
490	nzb = 0
491	nzt = nz
492	nzta = nz
493	nnz = nz
494
495	!
496	!-- For non-cyclic boundaries extend array u (v) by one gridpoint
497	IF ( bc_lr /= 'cyclic' ) uxrp = 1
498	IF ( bc_ns /= 'cyclic' ) vynp = 1
499
500	!
501	!-- Array bounds for the pressure solver (in the parallel code, these bounds
502	!-- are the ones for the transposed arrays)
503	nys_x = nys
504	nyn_x = nyn
505	nyn_xa = nyn
506	nzb_x = nzb + 1
507	nzt_x = nzt
508	nzt_xa = nzt
509
510	nxl_y = nxl
511	nxr_y = nxr
512	nxr_ya = nxr
513	nzb_y = nzb + 1
514	nzt_y = nzt
515	nzt_ya = nzt
516
517	nxl_z = nxl
518	nxr_z = nxr
519	nxr_za = nxr
520	nys_z = nys
521	nyn_z = nyn
522	nyn_za = nyn
523
524	#endif
525
526	!
527	!-- Calculate number of grid levels necessary for the multigrid poisson solver
528	!-- as well as the gridpoint indices on each level
529	IF ( psolver == 'multigrid' ) THEN
530
531	!
532	!-- First calculate number of possible grid levels for the subdomains
533	mg_levels_x = 1
534	mg_levels_y = 1
535	mg_levels_z = 1
536
537	i = nnx
538	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
539	i = i / 2
540	mg_levels_x = mg_levels_x + 1
541	ENDDO
542
543	j = nny
544	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
545	j = j / 2
546	mg_levels_y = mg_levels_y + 1
547	ENDDO
548
549	k = nnz
550	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
551	k = k / 2
552	mg_levels_z = mg_levels_z + 1
553	ENDDO
554
555	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
556
557	!
558	!-- Find out, if the total domain allows more levels. These additional
559	!-- levels are processed on PE0 only.
560	IF ( numprocs > 1 ) THEN
561	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
562	mg_switch_to_pe0_level_l = maximum_grid_level
563
564	mg_levels_x = 1
565	mg_levels_y = 1
566
567	i = nx+1
568	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
569	i = i / 2
570	mg_levels_x = mg_levels_x + 1
571	ENDDO
572
573	j = ny+1
574	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
575	j = j / 2
576	mg_levels_y = mg_levels_y + 1
577	ENDDO
578
579	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
580
581	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
582	mg_switch_to_pe0_level_l = maximum_grid_level_l - &
583	mg_switch_to_pe0_level_l + 1
584	ELSE
585	mg_switch_to_pe0_level_l = 0
586	ENDIF
587	ELSE
588	mg_switch_to_pe0_level_l = 0
589	maximum_grid_level_l = maximum_grid_level
590	ENDIF
591
592	!
593	!-- Use switch level calculated above only if it is not pre-defined
594	!-- by user
595	IF ( mg_switch_to_pe0_level == 0 ) THEN
596
597	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
598	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
599	maximum_grid_level = maximum_grid_level_l
600	ENDIF
601
602	ELSE
603	!
604	!-- Check pre-defined value and reset to default, if neccessary
605	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
606	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
607	IF ( myid == 0 ) THEN
608	PRINT*, '+++ WARNING init_pegrid: mg_switch_to_pe0_level ', &
609	'out of range and reset to default (=0)'
610	ENDIF
611	mg_switch_to_pe0_level = 0
612	ELSE
613	!
614	!-- Use the largest number of possible levels anyway and recalculate
615	!-- the switch level to this largest number of possible values
616	maximum_grid_level = maximum_grid_level_l
617
618	ENDIF
619	ENDIF
620
621	ENDIF
622
623	ALLOCATE( grid_level_count(maximum_grid_level), &
624	nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
625	nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
626	nzt_mg(maximum_grid_level) )
627
628	grid_level_count = 0
629	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
630
631	DO i = maximum_grid_level, 1 , -1
632
633	IF ( i == mg_switch_to_pe0_level ) THEN
634	#if defined( __parallel )
635	!
636	!-- Save the grid size of the subdomain at the switch level, because
637	!-- it is needed in poismg.
638	!-- Array bounds of the local subdomain grids are gathered on PE0
639	ind(1) = nxl_l; ind(2) = nxr_l
640	ind(3) = nys_l; ind(4) = nyn_l
641	ind(5) = nzt_l
642	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
643	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
644	MPI_INTEGER, comm2d, ierr )
645	DO j = 0, numprocs-1
646	DO k = 1, 5
647	mg_loc_ind(k,j) = ind_all(k+j*5)
648	ENDDO
649	ENDDO
650	DEALLOCATE( ind_all )
651	!
652	!-- Calculate the grid size of the total domain gathered on PE0
653	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
654	nxl_l = 0
655	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
656	nys_l = 0
657	!
658	!-- The size of this gathered array must not be larger than the
659	!-- array tend, which is used in the multigrid scheme as a temporary
660	!-- array
661	subdomain_size = ( nxr - nxl + 3 ) * ( nyn - nys + 3 ) * &
662	( nzt - nzb + 2 )
663	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
664	( nzt_l - nzb + 2 )
665
666	IF ( gathered_size > subdomain_size ) THEN
667	IF ( myid == 0 ) THEN
668	PRINT*, '+++ init_pegrid: not enough memory for storing ', &
669	'gathered multigrid data on PE0'
670	ENDIF
671	CALL local_stop
672	ENDIF
673	#else
674	PRINT*, '+++ init_pegrid: multigrid gather/scatter impossible ', &
675	'in non parallel mode'
676	CALL local_stop
677	#endif
678	ENDIF
679
680	nxl_mg(i) = nxl_l
681	nxr_mg(i) = nxr_l
682	nys_mg(i) = nys_l
683	nyn_mg(i) = nyn_l
684	nzt_mg(i) = nzt_l
685
686	nxl_l = nxl_l / 2
687	nxr_l = nxr_l / 2
688	nys_l = nys_l / 2
689	nyn_l = nyn_l / 2
690	nzt_l = nzt_l / 2
691	ENDDO
692
693	ELSE
694
695	maximum_grid_level = 1
696
697	ENDIF
698
699	grid_level = maximum_grid_level
700
701	#if defined( __parallel )
702	!
703	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
704	ngp_y = nyn - nys + 1
705
706	!
707	!-- Define a new MPI derived datatype for the exchange of ghost points in
708	!-- y-direction for 2D-arrays (line)
709	CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )
710	CALL MPI_TYPE_COMMIT( type_x, ierr )
711	CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_INTEGER, type_x_int, ierr )
712	CALL MPI_TYPE_COMMIT( type_x_int, ierr )
713
714	!
715	!-- Calculate gridpoint numbers for the exchange of ghost points along x
716	!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
717	!-- exchange of ghost points in y-direction (xz-plane).
718	!-- Do these calculations for the model grid and (if necessary) also
719	!-- for the coarser grid levels used in the multigrid method
720	ALLOCATE ( ngp_yz(maximum_grid_level), type_xz(maximum_grid_level) )
721
722	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
723
724	DO i = maximum_grid_level, 1 , -1
725	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
726
727	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
728	MPI_REAL, type_xz(i), ierr )
729	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
730
731	nxl_l = nxl_l / 2
732	nxr_l = nxr_l / 2
733	nys_l = nys_l / 2
734	nyn_l = nyn_l / 2
735	nzt_l = nzt_l / 2
736	ENDDO
737	#endif
738
739	#if defined( __parallel )
740	!
741	!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
742	!-- horizontal boundary conditions. Set variables for extending array u (v)
743	!-- by one gridpoint on the left/rightmost (northest/southest) processor
744	IF ( pleft == MPI_PROC_NULL ) THEN
745	IF ( bc_lr == 'dirichlet/neumann' ) THEN
746	inflow_l = .TRUE.
747	ELSEIF ( bc_lr == 'neumann/dirichlet' ) THEN
748	outflow_l = .TRUE.
749	ENDIF
750	ENDIF
751
752	IF ( pright == MPI_PROC_NULL ) THEN
753	IF ( bc_lr == 'dirichlet/neumann' ) THEN
754	outflow_r = .TRUE.
755	ELSEIF ( bc_lr == 'neumann/dirichlet' ) THEN
756	inflow_r = .TRUE.
757	ENDIF
758	uxrp = 1
759	ENDIF
760
761	IF ( psouth == MPI_PROC_NULL ) THEN
762	IF ( bc_ns == 'dirichlet/neumann' ) THEN
763	outflow_s = .TRUE.
764	ELSEIF ( bc_ns == 'neumann/dirichlet' ) THEN
765	inflow_s = .TRUE.
766	ENDIF
767	ENDIF
768
769	IF ( pnorth == MPI_PROC_NULL ) THEN
770	IF ( bc_ns == 'dirichlet/neumann' ) THEN
771	inflow_n = .TRUE.
772	ELSEIF ( bc_ns == 'neumann/dirichlet' ) THEN
773	outflow_n = .TRUE.
774	ENDIF
775	vynp = 1
776	ENDIF
777
778	!
779	!-- Additional MPI derived data type for the exchange of ghost points along x
780	!-- needed in case of non-cyclic boundary conditions along y on the northmost
781	!-- processors (for the exchange of the enlarged v array)
782	IF ( bc_ns /= 'cyclic' .AND. pnorth == MPI_PROC_NULL ) THEN
783	ngp_yz_p = ( nzt - nzb + 2 ) * ( nyn + vynp - nys + 3 )
784	CALL MPI_TYPE_VECTOR( nxr-nxl+3, nzt-nzb+2, ngp_yz_p, &
785	MPI_REAL, type_xz_p, ierr )
786	CALL MPI_TYPE_COMMIT( type_xz_p, ierr )
787	ENDIF
788	#else
789	IF ( bc_lr == 'dirichlet/neumann' ) THEN
790	inflow_l = .TRUE.
791	outflow_r = .TRUE.
792	uxrp = 1
793	ELSEIF ( bc_lr == 'neumann/dirichlet' ) THEN
794	outflow_l = .TRUE.
795	inflow_r = .TRUE.
796	ENDIF
797
798	IF ( bc_ns == 'dirichlet/neumann' ) THEN
799	inflow_n = .TRUE.
800	outflow_s = .TRUE.
801	ELSEIF ( bc_ns == 'neumann/dirichlet' ) THEN
802	outflow_n = .TRUE.
803	inflow_s = .TRUE.
804	vynp = 1
805	ENDIF
806	#endif
807
808	IF ( psolver == 'poisfft_hybrid' ) THEN
809	CALL poisfft_hybrid_ini
810	ELSE
811	CALL poisfft_init
812	ENDIF
813
814	END SUBROUTINE init_pegrid

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |