source: codes/icosagcm/trunk/src/domain.f90 @ 196

Last change on this file since 196 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 13.1 KB
Line 
1MODULE domain_mod
2  USE domain_param
3
4  TYPE t_domain
5    INTEGER :: face
6    INTEGER :: iim
7    INTEGER :: jjm
8    INTEGER :: ii_begin
9    INTEGER :: jj_begin
10    INTEGER :: ii_end
11    INTEGER :: jj_end
12    INTEGER :: ii_nb
13    INTEGER :: jj_nb
14    INTEGER :: ii_begin_glo
15    INTEGER :: jj_begin_glo
16    INTEGER :: ii_end_glo
17    INTEGER :: jj_end_glo
18    INTEGER,POINTER  :: assign_domain(:,:)
19    INTEGER,POINTER  :: assign_i(:,:)
20    INTEGER,POINTER  :: assign_j(:,:)
21    INTEGER,POINTER  :: edge_assign_domain(:,:,:)
22    INTEGER,POINTER  :: edge_assign_i(:,:,:)
23    INTEGER,POINTER  :: edge_assign_j(:,:,:)
24    INTEGER,POINTER  :: edge_assign_pos(:,:,:)
25    INTEGER,POINTER  :: edge_assign_sign(:,:,:)
26    REAL,POINTER     :: xyz(:,:,:)
27    REAL,POINTER     :: neighbour(:,:,:,:)
28    INTEGER,POINTER  :: delta(:,:)
29    REAL,POINTER     :: vertex(:,:,:,:)
30    LOGICAL,POINTER  :: own(:,:)
31    INTEGER,POINTER  :: ne(:,:,:)
32   
33    INTEGER :: t_right
34    INTEGER :: t_rup
35    INTEGER :: t_lup
36    INTEGER :: t_left
37    INTEGER :: t_ldown
38    INTEGER :: t_rdown
39
40    INTEGER :: u_right
41    INTEGER :: u_rup
42    INTEGER :: u_lup
43    INTEGER :: u_left
44    INTEGER :: u_ldown
45    INTEGER :: u_rdown
46
47    INTEGER :: z_rup
48    INTEGER :: z_up
49    INTEGER :: z_lup
50    INTEGER :: z_ldown
51    INTEGER :: z_down
52    INTEGER :: z_rdown
53   
54    INTEGER :: t_pos(6)
55    INTEGER :: u_pos(6)
56    INTEGER :: z_pos(6)
57     
58  END TYPE t_domain
59
60  INTEGER,SAVE :: ndomain
61  INTEGER,SAVE :: ndomain_glo
62
63  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:)
64  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:)
65
66  INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:)  ! size : ndomain_glo : mpi rank assigned to a domain
67  INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) ! size : ndomain_glo : corresponding local indice for a global domain indice
68  INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) ! size : ndomain : corresponding global indice for a local domain indice
69
70  LOGICAL, ALLOCATABLE, SAVE :: assigned_domain(:) ! size : ndomain : is an omp task is assigned to this domain
71!$OMP THREADPRIVATE(assigned_domain)
72   
73CONTAINS
74
75  SUBROUTINE create_domain
76  USE grid_param
77  USE mpipara
78  USE IOIPSL
79  IMPLICIT NONE
80  INTEGER :: ind,nf,ni,nj,i,j
81  INTEGER :: quotient, rest
82  INTEGER :: halo_i,halo_j
83 
84  TYPE(t_domain),POINTER :: d
85 
86    ndomain_glo=nsplit_i*nsplit_j*nb_face
87    ALLOCATE(domain_glo(ndomain_glo))
88    ALLOCATE(domglo_rank(ndomain_glo))
89    ALLOCATE(domglo_loc_ind(ndomain_glo))
90   
91    halo_i=0
92    CALL getin("halo_i",halo_i)
93    halo_j=1
94    CALL getin("halo_j",halo_j)
95
96    ind=0
97    DO nf=1,nb_face
98      DO nj=1,nsplit_j
99        DO ni=1,nsplit_i
100          ind=ind+1
101          d=>domain_glo(ind)
102          quotient=(iim_glo/nsplit_i)
103          rest=MOD(iim_glo,nsplit_i)
104          IF (ni-1 < rest) THEN
105            d%ii_nb=quotient+1
106            d%ii_begin_glo=(quotient+1)*(ni-1)+1
107          ELSE
108            d%ii_nb=quotient
109            d%ii_begin_glo=(quotient+1)*rest+(ni-1-rest) * quotient+1
110          ENDIF
111          d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1
112 
113          IF (ni/=1) THEN
114            d%ii_nb=d%ii_nb+1
115            d%ii_begin_glo=d%ii_begin_glo-1
116          ENDIF
117         
118       
119          quotient=(jjm_glo/nsplit_j)
120          rest=MOD(jjm_glo,nsplit_j)
121          IF (nj-1 < rest) THEN
122            d%jj_nb=quotient+1
123            d%jj_begin_glo=(quotient+1)*(nj-1)+1
124          ELSE
125            d%jj_nb=quotient
126            d%jj_begin_glo=(quotient+1)*rest+(nj-1-rest) * quotient+1
127          ENDIF
128          d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1
129
130          IF (nj/=1) THEN
131            d%jj_nb=d%jj_nb+1
132            d%jj_begin_glo=d%jj_begin_glo-1
133          ENDIF
134
135
136          d%iim=d%ii_nb+2*halo  + halo_i*2
137          d%jjm=d%jj_nb+2*halo  + halo_j*2
138          d%ii_begin=halo+1  + halo_i
139          d%jj_begin=halo+1  + halo_j
140          d%ii_end=d%ii_begin+d%ii_nb-1
141          d%jj_end=d%jj_begin+d%jj_nb-1
142          d%face=nf       
143          ALLOCATE(d%assign_domain(d%iim,d%jjm))
144          ALLOCATE(d%assign_i(d%iim,d%jjm))
145          ALLOCATE(d%assign_j(d%iim,d%jjm))
146          ALLOCATE(d%edge_assign_domain(0:5,d%iim,d%jjm))
147          ALLOCATE(d%edge_assign_i(0:5,d%iim,d%jjm))
148          ALLOCATE(d%edge_assign_j(0:5,d%iim,d%jjm))
149          ALLOCATE(d%edge_assign_pos(0:5,d%iim,d%jjm))
150          ALLOCATE(d%edge_assign_sign(0:5,d%iim,d%jjm))
151          ALLOCATE(d%delta(d%iim,d%jjm))
152          ALLOCATE(d%xyz(3,d%iim,d%jjm))
153          ALLOCATE(d%vertex(3,0:5,d%iim,d%jjm))
154          ALLOCATE(d%neighbour(3,0:5,d%iim,d%jjm))
155          ALLOCATE(d%own(d%iim,d%jjm))
156          ALLOCATE(d%ne(0:5,d%iim,d%jjm))
157         
158          IF (is_mpi_root) PRINT *,"Domain ",ind," : size ",d%ii_nb,"x",d%jj_nb
159         
160        END DO
161      END DO
162    END DO
163   
164  END SUBROUTINE create_domain
165 
166  SUBROUTINE copy_domain(d1,d2)
167  IMPLICIT NONE
168  INTEGER :: face
169  TYPE(t_domain),TARGET,INTENT(IN) :: d1
170  TYPE(t_domain), INTENT(OUT) :: d2
171 
172    d2%iim = d1%iim
173    d2%jjm = d1%jjm
174    d2%ii_begin = d1%ii_begin
175    d2%jj_begin = d1%jj_begin
176    d2%ii_end = d1%ii_end
177    d2%jj_end = d1%jj_end
178    d2%ii_nb =  d1%ii_nb
179    d2%jj_nb = d1%jj_nb
180    d2%ii_begin_glo = d1%ii_begin_glo
181    d2%jj_begin_glo = d1%jj_begin_glo
182    d2%ii_end_glo = d1%ii_end_glo
183    d2%jj_end_glo = d1%jj_end_glo
184    d2%assign_domain => d1%assign_domain
185    d2%assign_i => d1%assign_i
186    d2%assign_j => d1%assign_j
187    d2%edge_assign_domain => d1%edge_assign_domain
188    d2%edge_assign_i => d1%edge_assign_i
189    d2%edge_assign_j => d1%edge_assign_j
190    d2%edge_assign_pos => d1%edge_assign_pos
191    d2%edge_assign_sign => d1%edge_assign_sign
192    d2%xyz => d1%xyz
193    d2%neighbour => d1%neighbour
194    d2%delta => d1%delta
195    d2%vertex => d1%vertex
196    d2%own => d1%own
197    d2%ne => d1%ne
198   
199    d2%t_right = d1%t_right
200    d2%t_rup = d1%t_rup
201    d2%t_lup = d1%t_lup
202    d2%t_left = d1%t_left
203    d2%t_ldown = d1%t_ldown
204    d2%t_rdown = d1%t_rdown
205
206    d2%u_right = d1%u_right
207    d2%u_rup = d1%u_rup
208    d2%u_lup = d1%u_lup
209    d2%u_left = d1%u_left
210    d2%u_ldown = d1%u_ldown
211    d2%u_rdown = d1%u_rdown
212
213    d2%z_rup = d1%z_rup
214    d2%z_up = d1%z_up
215    d2%z_lup = d1%z_lup
216    d2%z_ldown = d1%z_ldown
217    d2%z_down = d1%z_down
218    d2%z_rdown = d1%z_rdown
219   
220    d2%t_pos = d1%t_pos
221    d2%u_pos = d1%u_pos
222    d2%z_pos = d1%z_pos
223   
224  END SUBROUTINE copy_domain
225 
226 
227  SUBROUTINE assign_cell
228  USE metric
229  IMPLICIT NONE
230    INTEGER :: ind_d,ind,ind2,e
231    INTEGER :: nf,nf2
232    INTEGER :: i,j,k,ii,jj
233    TYPE(t_domain),POINTER :: d
234    INTEGER :: delta
235     
236   
237    DO ind_d=1,ndomain_glo
238      d=>domain_glo(ind_d)
239      nf=d%face
240      DO j=d%jj_begin,d%jj_end
241        DO i=d%ii_begin,d%ii_end
242          ii=d%ii_begin_glo-d%ii_begin+i
243          jj=d%jj_begin_glo-d%jj_begin+j
244          ind=vertex_glo(ii,jj,nf)%ind
245          delta=vertex_glo(ii,jj,nf)%delta
246          IF (cell_glo(ind)%assign_face==nf) THEN
247            cell_glo(ind)%assign_domain=ind_d
248            cell_glo(ind)%assign_i=i
249            cell_glo(ind)%assign_j=j
250          ENDIF
251          IF ( i+1 <= d%ii_end ) CALL assign_edge(ind_d,ind,i,j,delta,0)
252          IF ( j+1 <= d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,1)
253          IF ( i-1 >= d%ii_begin .AND. j+1<=d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,2)
254          IF ( i-1 >= d%ii_begin ) CALL assign_edge(ind_d,ind,i,j,delta,3)
255          IF ( j-1 >= d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,4)
256          IF ( i+1 <= d%ii_end .AND. j-1 >=d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,5)
257        ENDDO
258      ENDDO
259    ENDDO
260   
261   
262    DO ind_d=1,ndomain_glo
263      d=>domain_glo(ind_d)
264      nf=d%face
265      DO j=d%jj_begin-1,d%jj_end+1
266        DO i=d%ii_begin-1,d%ii_end+1
267          ii=d%ii_begin_glo-d%ii_begin+i
268          jj=d%jj_begin_glo-d%jj_begin+j
269          ind=vertex_glo(ii,jj,nf)%ind
270          d%assign_domain(i,j)=cell_glo(ind)%assign_domain
271          d%assign_i(i,j)=cell_glo(ind)%assign_i
272          d%assign_j(i,j)=cell_glo(ind)%assign_j
273          delta=vertex_glo(ii,jj,nf)%delta
274          d%delta(i,j)=vertex_glo(ii,jj,nf)%delta
275          DO k=0,5
276            ind2=vertex_glo(ii,jj,nf)%neighbour(k)
277            d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:)
278
279            d%ne(k,i,j)=1-2*MOD(k,2)
280
281            e=cell_glo(ind)%edge(MOD(k+delta+6,6))
282            d%edge_assign_domain(k,i,j)=edge_glo(e)%assign_domain
283            d%edge_assign_i(k,i,j)=edge_glo(e)%assign_i
284            d%edge_assign_j(k,i,j)=edge_glo(e)%assign_j
285            d%edge_assign_pos(k,i,j)=edge_glo(e)%assign_pos
286            nf2=domain_glo(edge_glo(e)%assign_domain)%face
287            d%edge_assign_sign(k,i,j)=1-2*MOD(12+tab_index(nf,nf2,0),2)
288            IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) & 
289                /= MOD(edge_glo(e)%assign_pos+3,6)) THEN
290              d%edge_assign_sign(k,i,j)=-d%edge_assign_sign(k,i,j)
291            ENDIF
292             
293          ENDDO
294          d%xyz(:,i,j)=cell_glo(ind)%xyz(:)
295          IF (d%assign_domain(i,j)==ind_d) THEN
296           d%own(i,j)=.TRUE.
297          ELSE
298           d%own(i,j)=.FALSE.
299          ENDIF
300        ENDDO
301      ENDDO
302    ENDDO
303
304  CONTAINS
305
306    SUBROUTINE assign_edge(ind_d,ind,i,j,delta,k)
307      INTEGER :: ind_d,ind,i,j,delta,k
308      INTEGER :: e
309      e=cell_glo(ind)%edge(MOD(k+delta+6,6))
310      edge_glo(e)%assign_domain=ind_d
311      edge_glo(e)%assign_i=i
312      edge_glo(e)%assign_j=j
313      edge_glo(e)%assign_pos=k
314      edge_glo(e)%assign_delta=delta
315
316     END  SUBROUTINE assign_edge
317         
318  END SUBROUTINE assign_cell
319
320  SUBROUTINE compute_boundary
321  USE spherical_geom_mod
322  IMPLICIT NONE
323    INTEGER :: ind_d
324    INTEGER :: i,j,k
325    TYPE(t_domain),POINTER :: d 
326    REAL(rstd) :: ng1(3),ng2(3) 
327
328    DO ind_d=1,ndomain_glo
329      d=>domain_glo(ind_d)
330      DO j=d%jj_begin-1,d%jj_end+1
331        DO i=d%ii_begin-1,d%ii_end+1
332          DO k=0,5
333            ng1=d%neighbour(:,MOD(k,6),i,j)
334            ng2=d%neighbour(:,MOD(k+1,6),i,j)
335            IF (sqrt(sum((ng1-ng2)**2))<1e-16) ng2=d%neighbour(:,MOD(k+2,6),i,j)
336            CALL circumcenter(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j))
337          ENDDO
338        ENDDO
339      ENDDO
340    ENDDO       
341  END SUBROUTINE compute_boundary
342
343  SUBROUTINE set_neighbour_indice
344  USE metric
345  IMPLICIT NONE
346    INTEGER :: ind_d
347    TYPE(t_domain),POINTER :: d 
348   
349    DO ind_d=1,ndomain_glo
350      d=>domain_glo(ind_d)
351      d%t_right=1
352      d%t_left=-1
353      d%t_rup=d%iim
354      d%t_lup=d%iim-1
355      d%t_ldown=-d%iim
356      d%t_rdown=-d%iim+1
357     
358      d%u_right=0
359      d%u_lup=d%iim*d%jjm
360      d%u_ldown=2*d%iim*d%jjm
361     
362      d%u_rup=d%t_rup+d%u_ldown
363      d%u_left=d%t_left+d%u_right
364      d%u_rdown=d%t_rdown+d%u_lup
365     
366      d%z_up=0
367      d%z_down=d%iim*d%jjm
368      d%z_rup=d%t_rup+d%z_down
369      d%z_lup=d%t_lup+d%z_down
370      d%z_ldown=d%t_ldown+d%z_up
371      d%z_rdown=d%t_rdown+d%z_up
372     
373      d%t_pos(right)=d%t_right
374      d%t_pos(rup)=d%t_rup
375      d%t_pos(lup)=d%t_lup
376      d%t_pos(left)=d%t_left
377      d%t_pos(ldown)=d%t_ldown
378      d%t_pos(rdown)=d%t_rdown
379
380      d%u_pos(right)=d%u_right
381      d%u_pos(rup)=d%u_rup
382      d%u_pos(lup)=d%u_lup
383      d%u_pos(left)=d%u_left
384      d%u_pos(ldown)=d%u_ldown
385      d%u_pos(rdown)=d%u_rdown
386     
387      d%z_pos(vrup)=d%z_rup
388      d%z_pos(vup)=d%z_up
389      d%z_pos(vlup)=d%z_lup
390      d%z_pos(vldown)=d%z_ldown
391      d%z_pos(vdown)=d%z_down
392      d%z_pos(vrdown)=d%z_rdown
393     
394    ENDDO 
395 
396  END SUBROUTINE set_neighbour_indice
397     
398  SUBROUTINE assign_domain
399  USE mpipara
400  IMPLICIT NONE
401    INTEGER :: nb_domain(0:mpi_size-1)
402    INTEGER :: rank, ind,ind_glo
403   
404    DO rank=0,mpi_size-1
405      nb_domain(rank)=ndomain_glo/mpi_size
406      IF ( rank < MOD(ndomain_glo,mpi_size) ) nb_domain(rank)=nb_domain(rank)+1
407    ENDDO
408   
409    ndomain=nb_domain(mpi_rank)
410    ALLOCATE(domain(ndomain))
411    ALLOCATE(domloc_glo_ind(ndomain))
412   
413    rank=0
414    ind=0
415    DO ind_glo=1,ndomain_glo
416      ind=ind+1
417      domglo_rank(ind_glo)=rank
418      domglo_loc_ind(ind_glo)=ind
419      IF (rank==mpi_rank) THEN
420        CALL copy_domain(domain_glo(ind_glo),domain(ind))
421        domloc_glo_ind(ind)=ind_glo
422      ENDIF
423     
424      IF (ind==nb_domain(rank)) THEN
425        rank=rank+1
426        ind=0
427      ENDIF
428    ENDDO
429
430!$OMP PARALLEL
431    CALL assign_domain_omp
432!$OMP END PARALLEL
433       
434   
435  END SUBROUTINE  assign_domain     
436   
437  SUBROUTINE assign_domain_omp
438  USE omp_para
439  USE mpipara
440  IMPLICIT NONE
441    INTEGER :: nb_domain
442    INTEGER :: rank, ind, i
443
444    ALLOCATE(assigned_domain(ndomain))
445   
446    ind=0
447    DO rank=0,omp_size-1
448      nb_domain=ndomain/omp_size
449      IF ( rank < MOD(ndomain,omp_size) ) nb_domain=nb_domain+1
450   
451      DO i=1,nb_domain
452       ind=ind+1
453       IF (rank==omp_rank) THEN
454         assigned_domain(ind)=.TRUE.
455         PRINT *,"Rank ",mpi_rank," task ",rank," local domain : ",ind," global domain ",domloc_glo_ind(ind)
456       ELSE
457         assigned_domain(ind)=.FALSE.
458       ENDIF 
459      ENDDO
460   
461    ENDDO
462   
463  END SUBROUTINE assign_domain_omp
464   
465
466         
467  SUBROUTINE compute_domain
468  IMPLICIT NONE
469    CALL init_domain_param
470    CALL create_domain
471    CALL assign_cell
472    CALL compute_boundary
473    CALL set_neighbour_indice
474    CALL assign_domain
475     
476  END SUBROUTINE compute_domain
477         
478END MODULE domain_mod 
479 
Note: See TracBrowser for help on using the repository browser.