Home

Context Navigation

source: palm/trunk/SOURCE/pmc_particle_interface.f90 @ 2846

Last change on this file since 2846 was 2841, checked in by knoop, 7 years ago
Bugfix: wrong placement of include 'mpif.h' corrected
Property svn:keywords set to `Id`
File size: 42.9 KB

Line
1	MODULE pmc_particle_interface
2
3	!------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the
7	! terms of the GNU General Public License as published by the Free Software
8	! Foundation, either version 3 of the License, or (at your option) any later
9	! version.
10	!
11	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14	!
15	! You should have received a copy of the GNU General Public License along with
16	! PALM. If not, see <http://www.gnu.org/licenses/>.
17	!
18	! Copyright 1997-2018 Leibniz Universitaet Hannover
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------!
27	! $Id: pmc_particle_interface.f90 2841 2018-02-27 15:02:57Z raasch $
28	! Bugfix: wrong placement of include 'mpif.h' corrected
29	!
30	! 2809 2018-02-15 09:55:58Z schwenkel
31	! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE.
32	!
33	!
34	! Initial Version (by Klaus Ketelsen)
35	!
36	!
37	! Description:
38	! ------------
39	! Introduce particle transfer in nested models.
40	! Interface to palm lpm model to handel particle transfer between parent and
41	! child model.
42	!------------------------------------------------------------------------------!
43
44	USE, INTRINSIC :: ISO_C_BINDING
45
46	#if !defined( __mpifh )
47	USE MPI
48	#endif
49
50	USE kinds
51
52	USE pegrid, &
53	ONLY: myidx,myidy
54
55	USE indices, &
56	ONLY: nx, ny, nxl, nxr, nys, nyn, nzb, nzt, nxlg, nxrg, nysg, nyng, nbgp
57
58	USE grid_variables, &
59	ONLY: dx, dy
60
61	USE arrays_3d, &
62	ONLY: zw
63
64	USE control_parameters, &
65	ONLY: message_string
66
67	USE particle_attributes, &
68	ONLY: prt_count, particles, grid_particles, &
69	particle_type, number_of_particles, zero_particle, &
70	ibc_par_t, ibc_par_lr, ibc_par_ns, alloc_factor
71
72	USE lpm_pack_and_sort_mod
73
74	USE pmc_general, &
75	ONLY: pedef
76
77	USE pmc_parent, &
78	ONLY: children, pmc_s_fillbuffer, pmc_s_getdata_from_buffer, &
79	pmc_s_get_child_npes
80
81	USE pmc_child, &
82	ONLY: me, pmc_c_getbuffer, pmc_c_putbuffer
83
84	USE pmc_interface, &
85	ONLY: cpl_id, get_number_of_childs, nr_part, part_adr, nested_run, &
86	get_childid, get_child_edges, nr_partc, part_adrc, &
87	coarse_bound, coord_x, coord_y, cg, get_child_gridspacing, &
88	lower_left_coord_x, lower_left_coord_y
89
90	USE pmc_handle_communicator, &
91	ONLY: pmc_parent_for_child
92
93	USE pmc_mpi_wrapper, &
94	ONLY: pmc_send_to_parent, pmc_recv_from_child
95
96	USE lpm_exchange_horiz_mod, &
97	ONLY: realloc_particles_array
98
99	IMPLICIT NONE
100
101	#if defined( __mpifh )
102	INCLUDE "mpif.h"
103	#endif
104
105	PRIVATE
106	SAVE
107
108	TYPE coarse_particle_def
109	INTEGER(iwp) :: nr_particle
110
111	TYPE(particle_type),ALLOCATABLE,DIMENSION(:) :: parent_particles
112	END TYPE coarse_particle_def
113
114	INTEGER(iwp),PARAMETER :: Min_particles_per_column = 100 !<
115	INTEGER(iwp),PARAMETER :: max_nr_particle_in_rma_win = 100000 !<
116
117	INTEGER(iwp) :: nr_fine_in_coarse !< Number of fine grid cells in coarse grid (one direction)
118	INTEGER(iwp) :: particle_win_child !<
119
120	INTEGER(iwp),ALLOCATABLE,DIMENSION(:) :: particle_win_parent !<
121
122	TYPE(C_PTR), ALLOCATABLE,DIMENSION(:) :: buf_ptr !<
123
124	TYPE(particle_type), DIMENSION(:),POINTER :: particle_in_win !<
125
126	TYPE(coarse_particle_def),ALLOCATABLE,DIMENSION(:,:) :: coarse_particles !<
127
128
129
130	INTERFACE pmcp_g_init
131	MODULE PROCEDURE pmcp_g_init
132	END INTERFACE pmcp_g_init
133
134	INTERFACE pmcp_g_alloc_win
135	MODULE PROCEDURE pmcp_g_alloc_win
136	END INTERFACE pmcp_g_alloc_win
137
138	INTERFACE pmcp_c_get_particle_from_parent
139	MODULE PROCEDURE pmcp_c_get_particle_from_parent
140	END INTERFACE pmcp_c_get_particle_from_parent
141
142	INTERFACE pmcp_c_send_particle_to_parent
143	MODULE PROCEDURE pmcp_c_send_particle_to_parent
144	END INTERFACE pmcp_c_send_particle_to_parent
145
146	INTERFACE pmcp_p_fill_particle_win
147	MODULE PROCEDURE pmcp_p_fill_particle_win
148	END INTERFACE pmcp_p_fill_particle_win
149
150	INTERFACE pmcp_p_empty_particle_win
151	MODULE PROCEDURE pmcp_p_empty_particle_win
152	END INTERFACE pmcp_p_empty_particle_win
153
154	INTERFACE pmcp_g_print_number_of_particles
155	MODULE PROCEDURE pmcp_g_print_number_of_particles
156	END INTERFACE pmcp_g_print_number_of_particles
157
158	INTERFACE pmcp_p_delete_particles_in_fine_grid_area
159	MODULE PROCEDURE pmcp_p_delete_particles_in_fine_grid_area
160	END INTERFACE pmcp_p_delete_particles_in_fine_grid_area
161
162	PUBLIC pmcp_g_init, pmcp_g_alloc_win, pmcp_c_get_particle_from_parent
163	PUBLIC pmcp_c_send_particle_to_parent, pmcp_p_fill_particle_win, pmcp_g_print_number_of_particles
164	PUBLIC pmcp_p_empty_particle_win, pmcp_p_delete_particles_in_fine_grid_area
165
166	CONTAINS
167
168	!------------------------------------------------------------------------------!
169	! Description:
170	! ------------
171	!> general routine:
172	!> Initializing actions of the particle interface
173	!> check particle boundary conditions for the child models
174	!------------------------------------------------------------------------------!
175	SUBROUTINE pmcp_g_init
176
177	IMPLICIT NONE
178
179	INTEGER(iwp) :: nr_childs !< Number of child models of the current model
180
181	nr_childs = get_number_of_childs()
182	!
183	!-- Check if the current model has child models
184	IF ( nr_childs > 0 ) THEN
185	ALLOCATE( nr_part(nysg:nyng, nxlg:nxrg) )
186	ALLOCATE( part_adr(nysg:nyng, nxlg:nxrg) )
187	nr_part = 0
188	part_adr = 0
189	ENDIF
190
191	!
192	!-- Set the boundary conditions to nested for all non root (i.e child) models
193	IF ( cpl_id > 1 ) THEN
194
195	IF ( ibc_par_t /= 3 ) THEN
196	ibc_par_t = 3
197	message_string = 'In Child model: bc_par_t is automatically set to nested '
198	CALL message( 'pmcp_g_init ', 'PA0477', 0, 1, 0, 6, 0 ) ! PA number has to be adjusted
199	ENDIF
200
201	IF ( ibc_par_lr /= 3 ) THEN
202	ibc_par_lr = 3
203	message_string = 'In Child model: bc_par_lr is automatically set to nested '
204	CALL message( 'pmcp_g_init ', 'PA0478', 0, 1, 0, 6, 0 ) ! PA number has to be adjusted
205	ENDIF
206
207	IF ( ibc_par_ns /= 3 ) THEN
208	ibc_par_ns = 3
209	message_string = 'In Child model: bc_par_ns is automatically set to nested '
210	CALL message( 'pmcp_g_init ', 'PA0479', 0, 1, 0, 6, 0 ) ! PA number has to be adjusted
211	ENDIF
212
213	ENDIF
214
215	END SUBROUTINE pmcp_g_init
216	!------------------------------------------------------------------------------!
217	! Description:
218	! ------------
219	!> general routine:
220	!> allocate the MPI windows
221	!------------------------------------------------------------------------------!
222	SUBROUTINE pmcp_g_alloc_win
223
224	IMPLICIT NONE
225
226	INTEGER(iwp) :: m !< loop index
227	INTEGER(iwp) :: ierr !< error code
228	INTEGER(iwp) :: icl !< left boundary in coarse(parent) index space
229	INTEGER(iwp) :: icr !< right boundary in coarse(parent) index space
230	INTEGER(iwp) :: jcs !< south boundary in coarse(parent) index space
231	INTEGER(iwp) :: jcn !< north boundary in coarse(parent) index space
232	INTEGER(iwp) :: child_id !< Id of a child model
233	INTEGER(iwp) :: nr_childs !< Number of child models of the current model
234
235	INTEGER(KIND=MPI_ADDRESS_KIND) :: parsize !<
236	INTEGER(KIND=MPI_ADDRESS_KIND) :: winsize !<
237	TYPE(C_PTR), SAVE :: ptr !<
238
239	TYPE(particle_type),DIMENSION(:),POINTER :: win_buffer !<
240
241	INTEGER(iwp),DIMENSION(1) :: buf_shape !<
242
243	!
244	!-- If the model has a parent model prepare the structures for transfer
245	IF ( cpl_id > 1 ) THEN
246
247	parsize = STORAGE_SIZE(zero_particle)/8
248
249	CALL MPI_ALLOC_MEM( parsize , MPI_INFO_NULL, ptr, ierr )
250	buf_shape(1) = 1
251	CALL C_F_POINTER( ptr, win_buffer, buf_shape )
252	CALL MPI_WIN_CREATE( win_buffer, parsize, parsize, &
253	MPI_INFO_NULL, me%intra_comm, particle_win_child, &
254	ierr )
255
256	!
257	!-- Child domain boundaries in the parent index space
258	icl = coarse_bound(1)
259	icr = coarse_bound(2)
260	jcs = coarse_bound(3)
261	jcn = coarse_bound(4)
262
263	! ALLOCATE( coarse_particles(jcs:jcn, icl:icr) )
264
265	coarse_particles(:,:)%nr_particle = 0
266	ENDIF
267
268	!
269	!-- If the model has child models prepare the structures for transfer
270	nr_childs = get_number_of_childs()
271	IF ( nr_childs > 0 ) THEN
272	ALLOCATE( particle_win_parent(nr_childs) )
273	ALLOCATE( buf_ptr(nr_childs) )
274	DO m = 1, nr_childs
275	child_id = get_childid(m)
276	parsize = STORAGE_SIZE(zero_particle)/8
277
278	winsize = max_nr_particle_in_rma_win * parsize
279	CALL MPI_ALLOC_MEM( winsize , MPI_INFO_NULL, buf_ptr(m), ierr )
280	buf_shape(1) = max_nr_particle_in_rma_win
281	CALL C_F_POINTER( buf_ptr(m), win_buffer, buf_shape )
282	CALL MPI_WIN_CREATE( win_buffer, winsize, parsize, &
283	MPI_INFO_NULL, children(child_id)%intra_comm, &
284	particle_win_parent(m), ierr )
285	ENDDO
286	ENDIF
287
288	END SUBROUTINE pmcp_g_alloc_win
289	!------------------------------------------------------------------------------!
290	! Description:
291	! ------------
292	!> child routine:
293	!> Read/get particles out of the parent MPI window
294	!------------------------------------------------------------------------------!
295	SUBROUTINE pmcp_c_get_particle_from_parent
296
297	IMPLICIT NONE
298
299	INTEGER(iwp) :: i !< x grid index
300	INTEGER(iwp) :: icl !< left boundary in coarse(parent) index space
301	INTEGER(iwp) :: ierr !< error code
302	INTEGER(iwp) :: ij !< combined xy index for the buffer array
303	INTEGER(iwp) :: ip !< loop index (child PEs)
304	INTEGER(iwp) :: j !< y grid index
305	INTEGER(iwp) :: jcs !< south boundary in coarse(parent) index space
306	INTEGER(iwp) :: nr !< number of particles to receive from a parent box
307
308	INTEGER(KIND=MPI_ADDRESS_KIND) :: parsize !<
309	INTEGER(KIND=MPI_ADDRESS_KIND) :: target_disp !<
310
311	TYPE(pedef), POINTER :: ape !< TO_DO Klaus: give a description and better name of the variable
312
313	IF ( cpl_id > 1 ) THEN
314
315	CALL pmc_c_getbuffer( particle_transfer = .TRUE. ) !Get number of particle/column and offset in RMA window xx
316
317	!
318	!-- Wait for buffer to fill.
319	!
320	!-- The parent side (in pmc_s_fillbuffer) is filling the buffer in the MPI RMA window
321	!-- When the filling is complete, a MPI_BARRIER is called.
322	!-- The child is not allowd to access the parent-buffer before it is completely filled
323	!-- Synchronization is done implicitely in pmc_c_getbuffer and pmc_s_fillbuffer on the parent side
324
325	icl = coarse_bound(1)
326	jcs = coarse_bound(3)
327
328	DO ip = 1, me%inter_npes
329
330	ape => me%pes(ip)
331
332	DO ij = 1, ape%nrele
333	j = ape%locind(ij)%j + jcs - 1
334	i = ape%locind(ij)%i + icl - 1
335	nr = nr_partc(j,i)
336	IF ( nr > 0 ) THEN
337
338	CALL check_and_alloc_coarse_particle (i, j, nr)
339	parsize = STORAGE_SIZE(zero_particle)/8
340	target_disp = part_adrc(j,i) - 1
341	CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , ip - 1, 0, particle_win_child, ierr )
342	CALL MPI_GET( coarse_particles(j,i)%parent_particles, &
343	nr * parsize, MPI_BYTE, ip - 1, &
344	target_disp, nr * parsize, &
345	MPI_BYTE, particle_win_child, ierr )
346	CALL MPI_WIN_UNLOCK( ip - 1, particle_win_child, ierr )
347	ENDIF
348	coarse_particles(j,i)%nr_particle = nr
349	ENDDO
350	ENDDO
351
352	CALL c_copy_particle_to_child_grid
353	ENDIF
354
355	END SUBROUTINE pmcp_c_get_particle_from_parent
356	!------------------------------------------------------------------------------!
357	! Description:
358	! ------------
359	!> child routine:
360	!> Write/put particles into the parent MPI window
361	!------------------------------------------------------------------------------!
362	SUBROUTINE pmcp_c_send_particle_to_parent
363
364	IMPLICIT NONE
365
366	INTEGER(iwp) :: disp_offset !<
367	INTEGER(iwp) :: i !< x loop index
368	INTEGER(iwp) :: icl !< left boundary in coarse(parent) index space
369	INTEGER(iwp) :: icr !< right boundary in coarse(parent) index space
370	INTEGER(iwp) :: ierr !< error code
371	INTEGER(iwp) :: ij !< combined xy index for the buffer array
372	INTEGER(iwp) :: ip !< loop index (child PEs)
373	INTEGER(iwp) :: j !< y loop index
374	INTEGER(iwp) :: jcs !< south boundary in coarse(parent) index space
375	INTEGER(iwp) :: jcn !< north boundary in coarse(parent) index space
376	INTEGER(iwp) :: max_nr_particle_per_pe !< maximum number of particles per PE (depending on grid apect ratio)
377	INTEGER(iwp) :: n !< shorter variable name for nr_fine_in_coarse
378	INTEGER(iwp) :: nr !< shorter variabel name for nr_partc
379	INTEGER(iwp) :: pe_offset !< offset index of the current PE
380
381	INTEGER(KIND=MPI_ADDRESS_KIND) :: parsize !<
382	INTEGER(KIND=MPI_ADDRESS_KIND) :: target_disp !<
383
384	REAL(wp) :: eps=0.00001 !< used in calculations to avoid rounding errors
385	REAL(wp) :: xx !< number of fine grid cells inside a coarse grid cell in x-direction
386	REAL(wp) :: yy !< number of fine grid cells inside a coarse grid cell in y-direction
387
388	TYPE(pedef), POINTER :: ape !< TO_DO Klaus: give a description and better name of the variable
389
390	! TYPE(particle_type) :: dummy_part !< dummy particle (needed for size calculations)
391
392
393	IF ( cpl_id > 1 ) THEN
394	CALL c_copy_particle_to_coarse_grid
395
396	!
397	!-- Child domain boundaries in the parent index space
398
399	icl = coarse_bound(1)
400	icr = coarse_bound(2)
401	jcs = coarse_bound(3)
402	jcn = coarse_bound(4)
403
404	nr_partc = 0
405
406	DO i = icl, icr
407	DO j = jcs, jcn
408	nr_partc(j,i) = coarse_particles(j,i)%nr_particle
409	ENDDO
410	ENDDO
411	part_adrc = 0
412
413	!
414	!-- compute number of fine grid cells in coarse grid (one direction)
415	xx = (cg%dx + eps) / dx ! +eps to avoid rounding error
416	yy = (cg%dy + eps) / dy
417	nr_fine_in_coarse = MAX( INT(xx), INT(yy) )
418
419	IF ( MOD( coord_x(0), cg%dx ) /= 0.0 .OR. MOD( coord_y(0), cg%dy ) /= 0.0 ) THEN
420	nr_fine_in_coarse = nr_fine_in_coarse + 1
421	ENDIF
422
423	!
424	!-- Assign a number to my child PE to select different areas in the RMA window on server side
425	!-- With this number a square of child PEs is defined which share the same coarse grid cells
426
427	n = nr_fine_in_coarse ! local variable n to make folloing statements shorter
428	pe_offset = MOD( myidx, n ) * n + MOD( myidy, n)
429	max_nr_particle_per_pe = max_nr_particle_in_rma_win / ( n * n )
430	disp_offset = pe_offset * max_nr_particle_per_pe
431	parsize = STORAGE_SIZE(zero_particle)/8
432	DO ip = 1, me%inter_npes
433
434	ape => me%pes(ip)
435
436	target_disp = disp_offset
437	DO ij = 1, ape%nrele
438	j = ape%locind(ij)%j + jcs - 1
439	i = ape%locind(ij)%i + icl - 1
440	nr = nr_partc(j,i)
441	IF( nr > 0 ) THEN
442	IF ( target_disp + nr - disp_offset >= max_nr_particle_per_pe ) THEN
443	WRITE(9,*) 'RMA window too small on child ', &
444	target_disp + nr - disp_offset, &
445	max_nr_particle_per_pe, max_nr_particle_in_rma_win
446	message_string = 'RMA window too small on child'
447	CALL message( 'pmci_create_child_arrays', 'PA0480', 3, 2, 0, 6, 0 ) ! PA number has to be adjusted
448	ENDIF
449	CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , ip - 1, 0, particle_win_child, ierr )
450	CALL MPI_PUT( coarse_particles(j,i)%parent_particles, &
451	nr * parsize, MPI_BYTE, ip - 1, &
452	target_disp, nr * parsize , &
453	MPI_BYTE, particle_win_child, ierr )
454	CALL MPI_WIN_UNLOCK( ip - 1, particle_win_child, ierr )
455	part_adrc(j,i) = target_disp + 1
456	target_disp = target_disp + nr
457	ENDIF
458	ENDDO
459	ENDDO
460
461	CALL pmc_c_putbuffer ( particle_transfer = .TRUE. ) !Send new number of particle/column and offset to parent
462
463	ENDIF
464
465	END SUBROUTINE pmcp_c_send_particle_to_parent
466	!------------------------------------------------------------------------------!
467	! Description:
468	! ------------
469	!> parent routine:
470	!> write particles into the MPI window
471	!------------------------------------------------------------------------------!
472	SUBROUTINE pmcp_p_fill_particle_win
473
474	IMPLICIT NONE
475
476	LOGICAL :: active_particle !< Particles located in the fine/child grid area are marked as active (to be transferred)
477	LOGICAL,SAVE :: lfirst = .TRUE. !<
478
479	INTEGER(iwp) :: child_id !< id of the child model
480	INTEGER(iwp) :: i !< x grid box index
481	INTEGER(iwp) :: ierr !< error code
482	INTEGER(iwp) :: ij !< combined xy index for the buffer array
483	INTEGER(iwp) :: ip !< loop index (child PEs)
484	INTEGER(iwp) :: j !< y grid box index
485	INTEGER(iwp) :: k !< z grid box index
486	INTEGER(iwp) :: m !< loop index (number of childs)
487	INTEGER(iwp) :: n !< loop index (number of particles)
488	INTEGER(iwp) :: nr_part_col !< Number of particles to transfer per column
489	INTEGER(iwp) :: number_of_particles !<
490	INTEGER(iwp) :: pindex !<
491	INTEGER(iwp) :: tot_particle_count !< Total number of particles per child
492
493	REAL(wp) :: dx_child !< child grid spacing
494	REAL(wp) :: dy_child !< child grid spacing
495	REAL(wp) :: dz_child !< child grid spacing
496	REAL(wp) :: ny_coord !< north coordinate of child grid
497	REAL(wp) :: ny_coord_b !< north coordinate of child grid boundary
498	REAL(wp) :: lx_coord !< left coordinate of child grid
499	REAL(wp) :: lx_coord_b !< left coordinate of child grid boundary
500	REAL(wp) :: rx_coord !< right coordinate of child grid
501	REAL(wp) :: rx_coord_b !< right coordinate of child grid boundary
502	REAL(wp) :: sy_coord !< south coordinate of child grid
503	REAL(wp) :: sy_coord_b !< south coordinate of child grid boundary
504	REAL(wp) :: uz_coord !< upper coordinate of child grid
505	REAL(wp) :: uz_coord_b !< upper coordinate of child grid boundary
506	REAL(wp) :: x !< particle position
507	REAL(wp) :: xo !< origin of particle
508	REAL(wp) :: y !< particle position
509	REAL(wp) :: yo !< origin of particle
510	REAL(wp) :: z !< particle position
511
512	INTEGER(iwp),DIMENSION(1) :: buf_shape !<
513
514	TYPE(pedef), POINTER :: ape !< TO_DO Klaus: give a description and better name of the variable
515
516	DO m = 1, get_number_of_childs()
517
518	child_id = get_childid(m)
519
520	CALL get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b, &
521	sy_coord, sy_coord_b, ny_coord, ny_coord_b, &
522	uz_coord, uz_coord_b)
523
524	CALL get_child_gridspacing( m, dx_child, dy_child, dz_child )
525
526	IF ( lfirst ) THEN
527	WRITE(9,'(a,5f10.2)') 'edges ',lx_coord,rx_coord,sy_coord,ny_coord,uz_coord
528	WRITE(9,'(a,5f10.2)') 'edges boundary ',lx_coord_b,rx_coord_b,sy_coord_b,ny_coord_b,uz_coord_b
529	WRITE(9,'(a,5f10.2)') 'child spacing ',dx_child, dy_child, dz_child,lower_left_coord_x,lower_left_coord_y
530	ENDIF
531	!
532	!-- reset values for every child
533	tot_particle_count = 0
534	nr_part = 0
535	part_adr = 0
536	pindex = 1
537
538	buf_shape(1) = max_nr_particle_in_rma_win
539	CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
540
541	DO ip = 1, children(child_id)%inter_npes
542
543	ape => children(child_id)%pes(ip)
544
545	nr_part_col = 0
546
547	DO ij = 1, ape%nrele
548
549	!
550	!-- Inside the PMC adressing of 3d arrays starts with 1
551	i = ape%locind(ij)%i + nxl - nbgp - 1
552	j = ape%locind(ij)%j + nys - nbgp - 1
553	nr_part_col = 0 ! Number of particles to transfer per column
554	part_adr(j,i) = pindex
555
556	DO k = nzb + 1, nzt
557	number_of_particles = prt_count(k,j,i)
558
559	IF ( number_of_particles <= 0 ) CYCLE
560
561	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
562
563	! Select particles within boundary area
564
565	DO n = 1, number_of_particles
566	x = particles(n)%x
567	y = particles(n)%y
568	z = particles(n)%z
569	!
570	!-- check if the particle is located in the fine grid area
571	active_particle = ((x > lx_coord .AND. x < rx_coord) .AND. &
572	(y > sy_coord .AND. y < ny_coord) .AND. &
573	(z > 0.000000001 .AND. z < uz_coord))
574	IF ( active_particle .AND. particles(n)%particle_mask ) THEN
575
576	particle_in_win(pindex) = particles(n)
577	!
578	!-- Change particle positions and origin relative to global origin
579	particle_in_win(pindex)%x = particle_in_win(pindex)%x + lower_left_coord_x
580	particle_in_win(pindex)%y = particle_in_win(pindex)%y + lower_left_coord_y
581	particle_in_win(pindex)%origin_x = particle_in_win(pindex)%origin_x + lower_left_coord_x
582	particle_in_win(pindex)%origin_y = particle_in_win(pindex)%origin_y + lower_left_coord_y
583
584	tot_particle_count = tot_particle_count + 1
585	nr_part_col = nr_part_col + 1
586	pindex = pindex + 1
587	IF ( pindex > max_nr_particle_in_rma_win ) THEN
588	WRITE(9,*) 'RMA window too small on parent ',pindex, max_nr_particle_in_rma_win
589	message_string = 'RMA window too small on parent'
590	CALL message( 'pmci_create_child_arrays', 'PA0481', 3, 2, 0, 6, 0 ) ! PA number has to be adjusted
591	ENDIF
592	END IF
593	ENDDO
594	ENDDO
595	nr_part(j,i) = nr_part_col
596	ENDDO
597	ENDDO
598
599	CALL pmc_s_fillbuffer( child_id, particle_transfer = .TRUE. )
600	ENDDO
601
602	lfirst = .FALSE.
603
604	END SUBROUTINE pmcp_p_fill_particle_win
605
606	!------------------------------------------------------------------------------!
607	! Description:
608	! ------------
609	!> parent routine:
610	!> delete particles from the MPI window
611	!------------------------------------------------------------------------------!
612	SUBROUTINE pmcp_p_empty_particle_win
613	IMPLICIT NONE
614
615	INTEGER(iwp) :: child_id !< model id of the child
616	INTEGER(iwp) :: ip !< loop index (child PEs)
617	INTEGER(iwp) :: m !< loop index (number of childs)
618
619	INTEGER(iwp),DIMENSION(1) :: buf_shape !<
620
621	DO m = 1, get_number_of_childs()
622
623	child_id = get_childid(m)
624
625	buf_shape(1) = max_nr_particle_in_rma_win
626	CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
627
628	!
629	!-- In some cells of the coarse grid, there are contributions from more than one child process
630	!-- Therfore p_copy_particle_to_org_grid is done for one child process per call
631
632	DO ip = 1, pmc_s_get_child_npes( child_id )
633	nr_part = 0
634	part_adr = 0
635
636	CALL pmc_s_getdata_from_buffer( child_id, particle_transfer = .TRUE., child_process_nr = ip )
637
638	CALL p_copy_particle_to_org_grid( m, child_id )
639	ENDDO
640
641	ENDDO
642
643	END SUBROUTINE pmcp_p_empty_particle_win
644
645	!------------------------------------------------------------------------------!
646	! Description:
647	! ------------
648	!> parent routine:
649	!> After the transfer mark all parent particles that are still inside on of the
650	!> child areas for deletion.
651	!------------------------------------------------------------------------------!
652	SUBROUTINE pmcp_p_delete_particles_in_fine_grid_area
653
654	IMPLICIT NONE
655
656	LOGICAL :: to_delete !< particles outside of model domain are marked as to_delete
657
658	INTEGER(iwp) :: i !< loop index (x grid)
659	INTEGER(iwp) :: j !< loop index (y grid)
660	INTEGER(iwp) :: k !< loop index (z grid)
661	INTEGER(iwp) :: m !< loop index (number of particles)
662	INTEGER(iwp) :: n !< loop index (number of childs)
663
664	REAL(wp) :: dx_child !< child grid spacing
665	REAL(wp) :: dy_child !< child grid spacing
666	REAL(wp) :: dz_child !< child grid spacing
667	REAL(wp) :: ny_coord !< north coordinate of child grid
668	REAL(wp) :: ny_coord_b !< north coordinate of child grid boundary
669	REAL(wp) :: lx_coord !< left coordinate of child grid
670	REAL(wp) :: lx_coord_b !< left coordinate of child grid boundary
671	REAL(wp) :: rx_coord !< right coordinate of child grid
672	REAL(wp) :: rx_coord_b !< right coordinate of child grid boundary
673	REAL(wp) :: sy_coord !< south coordinate of child grid
674	REAL(wp) :: sy_coord_b !< south coordinate of child grid boundary
675	REAL(wp) :: uz_coord !< upper coordinate of child grid
676	REAL(wp) :: uz_coord_b !< upper coordinate of child grid boundary
677	REAL(wp) :: x !< particle position
678	REAL(wp) :: y !< particle position
679	REAL(wp) :: z !< particle position
680
681	DO m = 1, get_number_of_childs()
682	CALL get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b, &
683	sy_coord, sy_coord_b, ny_coord, ny_coord_b, &
684	uz_coord, uz_coord_b )
685
686
687	CALL get_child_gridspacing( m, dx_child, dy_child, dz_child )
688
689	DO i = nxl, nxr
690	DO j = nys, nyn
691	DO k = nzb, nzt
692	number_of_particles = prt_count(k,j,i)
693
694	IF ( number_of_particles == 0 ) CYCLE
695
696	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
697
698	DO n = 1, number_of_particles
699	x = particles(n)%x
700	y = particles(n)%y
701	z = particles(n)%z
702
703	to_delete = ((x > lx_coord .AND. x < rx_coord) .AND. &
704	(y > sy_coord .AND. y < ny_coord) .AND. &
705	(z > 0.000000001 .AND. z < uz_coord))
706
707	IF ( to_delete ) THEN
708	particles(n)%particle_mask = .FALSE.
709	ENDIF
710	ENDDO
711	ENDDO
712	ENDDO
713	ENDDO
714
715	ENDDO
716
717	END SUBROUTINE pmcp_p_delete_particles_in_fine_grid_area
718	!------------------------------------------------------------------------------!
719	! Description:
720	! ------------
721	!> general routine:
722	!> print the total number of of the current model and its child models into a file
723	!------------------------------------------------------------------------------!
724	SUBROUTINE pmcp_g_print_number_of_particles (my_time,local_nr_particles)
725
726	USE pegrid, &
727	ONLY: myid
728
729	IMPLICIT NONE
730
731	INTEGER(iwp),INTENT(IN) :: local_nr_particles !<
732
733	REAL(wp),INTENT(IN) :: my_time !<
734
735	LOGICAL, SAVE :: is_file_open=.FALSE. !<
736
737	INTEGER(iwp) :: child_id !<
738	INTEGER(iwp) :: child_nr_particles !< total number of particles in all child models
739	INTEGER(iwp) :: ierr !< error code
740	INTEGER(iwp) :: m !< loop index (number of childs)
741
742	CHARACTER(LEN=16) :: fname !< filename
743
744	INTEGER(iwp),DIMENSION(2) :: ivalr !< integer value to be received
745	INTEGER(iwp),DIMENSION(2) :: ivals !< integer value to be send
746
747
748	child_nr_particles = 0
749	IF ( myid == 0 ) THEN
750	IF ( cpl_id > 1 ) THEN
751	ivals(1) = local_nr_particles
752	CALL pmc_send_to_parent( ivals, 1, 0, 400, ierr )
753	ENDIF
754	DO m = 1, SIZE( pmc_parent_for_child ) - 1
755
756	child_id = pmc_parent_for_child(m)
757	CALL pmc_recv_from_child( child_id, ivalr, 1, 0, 400, ierr )
758	child_nr_particles = child_nr_particles + ivalr(1)
759	ENDDO
760
761	IF ( SIZE( pmc_parent_for_child ) > 1 ) THEN
762	IF ( .NOT. is_file_open ) THEN !kk muss noch auf file_open umgestellt werden
763	WRITE(fname,'(a,i2.2)') 'nr_particles_',cpl_id
764	OPEN (333,file = fname)
765	is_file_open = .true.
766	ENDIF
767	WRITE(333,'(f12.2,3i10)') my_time,local_nr_particles + child_nr_particles,local_nr_particles,child_nr_particles
768	ENDIF
769	ENDIF
770
771	END SUBROUTINE pmcp_g_print_number_of_particles
772
773	!------------------------------------------------------------------------------!
774	!------------------------------------------------------------------------------!
775	! Private subroutines
776	!------------------------------------------------------------------------------!
777	!------------------------------------------------------------------------------!
778
779	!------------------------------------------------------------------------------!
780	! Description:
781	! ------------
782	!> child routine
783	!> update the size of the local buffer (coarse_particles)
784	!------------------------------------------------------------------------------!
785	SUBROUTINE check_and_alloc_coarse_particle (ic,jc,nr,with_copy)
786
787	IMPLICIT NONE
788
789	INTEGER(iwp),INTENT(IN) :: ic !< coarse x grid index
790	INTEGER(iwp),INTENT(IN) :: jc !< coarse y grid index
791	INTEGER(iwp),INTENT(IN) :: nr !< number of particles to be transferred/stored in local buffer
792
793	LOGICAL,INTENT(IN),OPTIONAL :: with_copy !< copy particles in buffer? or reallocate empty buffer
794
795	LOGICAL :: with_copy_lo !< local variable of with copy
796
797	INTEGER(iwp) :: new_size !< new size of the local buffer
798	INTEGER(iwp) :: old_size !< old size of the local buffer
799
800	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
801
802	with_copy_lo = .FALSE.
803	IF ( PRESENT( with_copy ) ) with_copy_lo = with_copy
804
805	IF ( .NOT. ALLOCATED( coarse_particles(jc,ic)%parent_particles ) ) THEN
806	new_size = MAX( nr, Min_particles_per_column )
807	ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
808	ELSEIF ( nr > SIZE( coarse_particles(jc,ic)%parent_particles ) ) THEN
809
810	old_size = SIZE( coarse_particles(jc,ic)%parent_particles )
811	new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
812	new_size = MAX( nr, new_size, old_size + Min_particles_per_column )
813
814	!
815	!-- Copy existing particles to new particle buffer
816	IF ( with_copy_lo ) THEN
817	ALLOCATE( tmp_particles_d(old_size) )
818	tmp_particles_d(1:old_size) = coarse_particles(jc,ic)%parent_particles
819
820	DEALLOCATE( coarse_particles(jc,ic)%parent_particles )
821	ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
822
823	coarse_particles(jc,ic)%parent_particles(1:old_size) = tmp_particles_d(1:old_size)
824	coarse_particles(jc,ic)%parent_particles(old_size+1:new_size) = zero_particle
825
826	DEALLOCATE( tmp_particles_d )
827	!
828	!-- allocate or reallocate an empty buffer
829	ELSE
830	DEALLOCATE( coarse_particles(jc,ic)%parent_particles )
831	ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
832	ENDIF
833	ENDIF
834
835	END SUBROUTINE check_and_alloc_coarse_particle
836
837	!------------------------------------------------------------------------------!
838	! Description:
839	! ------------
840	!> child routine:
841	!> copy/sort particles out of the local buffer into the respective grid boxes
842	!------------------------------------------------------------------------------!
843	SUBROUTINE c_copy_particle_to_child_grid
844
845	IMPLICIT NONE
846
847	INTEGER(iwp) :: ic !< coarse x grid index
848	INTEGER(iwp) :: icl !< left boundary in coarse(parent) index space
849	INTEGER(iwp) :: icr !< right boundary in coarse(parent) index space
850	INTEGER(iwp) :: ip !< x grid index
851	INTEGER(iwp) :: jc !< coarse y grid index
852	INTEGER(iwp) :: jcn !< north boundary in coarse(parent) index space
853	INTEGER(iwp) :: jcs !< south boundary in coarse(parent) index space
854	INTEGER(iwp) :: jp !< y grid index
855	INTEGER(iwp) :: kp !< z grid index
856	INTEGER(iwp) :: n !< loop index (number of particles)
857	INTEGER(iwp) :: nr !< short variable for name number or particles
858
859	REAL(wp) :: xc !< child x coordinate
860	REAL(wp) :: xoc !< child x origin
861	REAL(wp) :: yc !< child y coordinate
862	REAL(wp) :: yoc !< child y origin
863	REAL(wp) :: zc !< child z coordinate
864
865	!
866	!-- Child domain boundaries in the parent index space
867	icl = coarse_bound(1)
868	icr = coarse_bound(2)
869	jcs = coarse_bound(3)
870	jcn = coarse_bound(4)
871
872	DO ic = icl, icr
873	DO jc = jcs, jcn
874	nr = coarse_particles(jc,ic)%nr_particle
875
876	IF ( nr > 0 ) THEN
877	DO n = 1, nr
878	xc = coarse_particles(jc,ic)%parent_particles(n)%x-lower_left_coord_x
879	yc = coarse_particles(jc,ic)%parent_particles(n)%y-lower_left_coord_y
880	zc = coarse_particles(jc,ic)%parent_particles(n)%z
881	xoc = coarse_particles(jc,ic)%parent_particles(n)%origin_x-lower_left_coord_x
882	yoc = coarse_particles(jc,ic)%parent_particles(n)%origin_y-lower_left_coord_y
883	ip = xc / dx
884	jp = yc / dy
885	kp = nzt
886	DO WHILE ( zw(kp-1) > zc .AND. kp > nzb + 1 ) ! kk search loop has to be optimzed !!!
887	kp = kp - 1
888	ENDDO
889
890	IF ( ip >= nxl .AND. ip <= nxr .AND. jp >= nys .AND. jp <= nyn ) THEN
891	prt_count(kp,jp,ip) = prt_count(kp,jp,ip) + 1
892	IF ( prt_count(kp,jp,ip) > SIZE( grid_particles(kp,jp,ip)%particles ) ) THEN
893	CALL realloc_particles_array( ip, jp, kp )
894	ENDIF
895	coarse_particles(jc,ic)%parent_particles(n)%x = xc ! Adjust coordinates to child grid
896	coarse_particles(jc,ic)%parent_particles(n)%y = yc
897	coarse_particles(jc,ic)%parent_particles(n)%origin_x = xoc ! Adjust origins to child grid
898	coarse_particles(jc,ic)%parent_particles(n)%origin_y = yoc
899	grid_particles(kp,jp,ip)%particles(prt_count(kp,jp,ip)) = coarse_particles(jc,ic)%parent_particles(n)
900	ENDIF
901	ENDDO
902	ENDIF
903	ENDDO
904	ENDDO
905
906	END SUBROUTINE c_copy_particle_to_child_grid
907	!------------------------------------------------------------------------------!
908	! Description:
909	! ------------
910	!> child routine:
911	!> copy particles which left the model area into the local buffer
912	!------------------------------------------------------------------------------!
913	SUBROUTINE c_copy_particle_to_coarse_grid
914
915	IMPLICIT NONE
916
917	LOGICAL :: boundary_particle !<
918
919	INTEGER(iwp) :: i !< loop index (x grid)
920	INTEGER(iwp) :: ic !< loop index (coarse x grid)
921	INTEGER(iwp) :: icl !< left boundary in coarse(parent) index space
922	INTEGER(iwp) :: icr !< left boundary in coarse(parent) index space
923	INTEGER(iwp) :: ierr !< error code
924	INTEGER(iwp) :: j !< loop index (y grid)
925	INTEGER(iwp) :: jc !< loop index (coarse y grid)
926	INTEGER(iwp) :: jcs !< south boundary in coarse(parent) index space
927	INTEGER(iwp) :: jcn !< north boundary in coarse(parent) index space
928	INTEGER(iwp) :: k !< loop index (z grid)
929	INTEGER(iwp) :: n !< loop index (number of particles)
930
931	REAL(iwp) :: x !< x coordinate
932	REAL(iwp) :: xo !< x origin
933	REAL(iwp) :: x_left !< absolute left boundary
934	REAL(iwp) :: x_right !< absolute right boundary
935	REAL(iwp) :: y !< left boundary in coarse(parent) index space
936	REAL(iwp) :: yo !< y origin
937	REAL(iwp) :: y_north !< absolute north boundary
938	REAL(iwp) :: y_south !< absoulte south boundary
939	REAL(iwp) :: z !< z coordinate
940
941	!
942	!-- Child domain boundaries in the parent index space
943
944	icl = coarse_bound(1)
945	icr = coarse_bound(2)
946	jcs = coarse_bound(3)
947	jcn = coarse_bound(4)
948
949	!
950	!-- absolute coordinates
951	x_left = coord_x(0)
952	x_right = coord_x(nx) + dx
953	y_south = coord_y(0)
954	y_north = coord_y(ny) + dy
955
956	! Clear Particle Buffer
957
958	DO ic = icl, icr
959	DO jc = jcs, jcn
960	coarse_particles(jc,ic)%nr_particle = 0
961	ENDDO
962	ENDDO
963
964	!
965	!-- Determine particles which leave the inner area in east or west dirextion
966	!-- Compute only first (nxl) and last (nxr) loop iterration.
967	DO i = nxl, nxr
968	DO j = nys, nyn
969	DO k = nzb + 1, nzt
970	number_of_particles = prt_count(k,j,i)
971	IF ( number_of_particles <= 0 ) CYCLE
972	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
973	DO n = 1, number_of_particles
974	IF ( particles(n)%particle_mask ) THEN
975	x = particles(n)%x+ lower_left_coord_x
976	y = particles(n)%y+ lower_left_coord_y
977	xo = particles(n)%origin_x + lower_left_coord_x
978	yo = particles(n)%origin_y + lower_left_coord_y
979	z = particles(n)%z
980
981	boundary_particle = .FALSE.
982	boundary_particle = boundary_particle .OR. (x < x_left .OR. x > x_right)
983	boundary_particle = boundary_particle .OR. (y < y_south .OR. y > y_north)
984	boundary_particle = boundary_particle .OR. (z > zw(nzt))
985
986	IF ( boundary_particle ) THEN
987	ic = x / cg%dx !TODO anpassen auf Mehrfachnesting
988	jc = y / cg%dy
989
990	IF ( ic >= icl .AND. ic <= icr .AND. jc >= jcs .AND. jc <= jcn ) THEN
991	coarse_particles(jc,ic)%nr_particle = coarse_particles(jc,ic)%nr_particle + 1
992	CALL check_and_alloc_coarse_particle( ic, jc, coarse_particles(jc,ic)%nr_particle, with_copy=.TRUE. )
993
994	coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle) = particles(n)
995	coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%x = x !adapt to global coordinates
996	coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%y = y
997	coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%origin_x = xo
998	coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%origin_y = yo
999	!-- Mark particle as deleted after copying it to the transfer buffer
1000	grid_particles(k,j,i)%particles(n)%particle_mask = .FALSE.
1001	ELSE
1002	WRITE(9,'(a,10i6)') 'This should not happen ',i,j,k,ic,jc,icl,icr,jcs,jcn
1003	CALL MPI_Abort( MPI_COMM_WORLD, 9999, ierr )
1004	ENDIF
1005	ENDIF
1006	ENDIF
1007	ENDDO
1008	ENDDO
1009	ENDDO
1010	ENDDO
1011
1012	!
1013	!- Pack particles (eliminate those marked for deletion),
1014	!- determine new number of particles
1015	CALL lpm_sort_in_subboxes
1016
1017	END SUBROUTINE c_copy_particle_to_coarse_grid
1018	!------------------------------------------------------------------------------!
1019	! Description:
1020	! ------------
1021	!> parent routine:
1022	!> copy/sort particles from the MPI window into the respective grid boxes
1023	!------------------------------------------------------------------------------!
1024	SUBROUTINE p_copy_particle_to_org_grid( m, child_id )
1025
1026	IMPLICIT NONE
1027
1028	INTEGER(iwp),INTENT(IN) :: child_id !<
1029	INTEGER(iwp),INTENT(IN) :: m !<
1030
1031	INTEGER(iwp) :: i !< loop index (x grid)
1032	INTEGER(iwp) :: j !< loop index (y grid)
1033	INTEGER(iwp) :: k !< loop index (z grid)
1034	INTEGER(iwp) :: n !< loop index (nr part)
1035	INTEGER(iwp) :: nr !< short variable name for nr_part
1036	INTEGER(iwp) :: pindex !< short variable name part_adr
1037
1038	REAL(iwp) :: x !< x coordinate
1039	REAL(iwp) :: xo !< x origin
1040	REAL(iwp) :: y !< y coordinate
1041	REAL(iwp) :: yo !< y origin
1042	REAL(iwp) :: z !< z coordinate
1043
1044	INTEGER(iwp),DIMENSION(1) :: buf_shape !<
1045
1046	buf_shape(1) = max_nr_particle_in_rma_win
1047	CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
1048
1049	DO i = nxl, nxr
1050	DO j = nys, nyn
1051	nr = nr_part(j,i)
1052	IF ( nr > 0 ) THEN
1053	pindex = part_adr(j,i)
1054	DO n = 1, nr
1055	x = particle_in_win(pindex)%x-lower_left_coord_x
1056	y = particle_in_win(pindex)%y-lower_left_coord_y
1057	z = particle_in_win(pindex)%z
1058	xo = particle_in_win(pindex)%origin_x-lower_left_coord_x
1059	yo = particle_in_win(pindex)%origin_y-lower_left_coord_y
1060	k = nzt + 1
1061	DO WHILE ( zw(k-1) > z .AND. k > nzb + 1 ) ! kk search loop has to be optimzed !!!
1062	k = k - 1
1063	END DO
1064
1065	prt_count(k,j,i) = prt_count(k,j,i) + 1
1066	IF ( prt_count(k,j,i) > SIZE( grid_particles(k,j,i)%particles ) ) THEN
1067	CALL realloc_particles_array( i, j, k )
1068	ENDIF
1069	grid_particles(k,j,i)%particles(prt_count(k,j,i)) = particle_in_win(pindex)
1070
1071	!
1072	!-- Update particle positions and origins relative to parent domain
1073	grid_particles(k,j,i)%particles(prt_count(k,j,i))%x = x
1074	grid_particles(k,j,i)%particles(prt_count(k,j,i))%y = y
1075	grid_particles(k,j,i)%particles(prt_count(k,j,i))%origin_x = xo
1076	grid_particles(k,j,i)%particles(prt_count(k,j,i))%origin_y = yo
1077	pindex = pindex + 1
1078	ENDDO
1079	ENDIF
1080	ENDDO
1081	ENDDO
1082
1083	END SUBROUTINE p_copy_particle_to_org_grid
1084
1085
1086	END MODULE pmc_particle_interface

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |