source: codes/icosagcm/devel/src/parallel/domain.f90 @ 603

Last change on this file since 603 was 533, checked in by dubos, 7 years ago

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