SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
two_dim_metis_distribution.f90
Go to the documentation of this file.
3 PRIVATE
4  REAL(KIND=8) :: epsilon = 1.d-10
5 !!$ Dummy for metis...
6  INTEGER :: METIS_NOPTIONS=40, METIS_OPTION_NUMBERING=17
7 !!$ Dummy for metis...
8 #include "petsc/finclude/petsc.h"
9 CONTAINS
10  SUBROUTINE reorder_mesh(communicator,nb_proc,mesh,mesh_loc,list_of_interfaces)
11  USE def_type_mesh
12  USE my_util
13  USE sub_plot
14  IMPLICIT NONE
15  TYPE(mesh_type) :: mesh, mesh_loc
16  INTEGER, DIMENSION(2) :: me_loc, mes_loc, np_loc
17  INTEGER, DIMENSION(:), POINTER, OPTIONAL :: list_of_interfaces
18  INTEGER, DIMENSION(mesh%me+1) :: xind
19  INTEGER, DIMENSION(mesh%me) :: vwgt, part, old_m_to_new
20  LOGICAL, DIMENSION(mesh%np) :: virgin
21  INTEGER, DIMENSION(mesh%np) :: old_j_to_new
22  INTEGER, DIMENSION(mesh%mes) :: old_ms_to_new, parts
23  INTEGER, DIMENSION(2,mesh%mes) :: inter_news
24  INTEGER, DIMENSION(SIZE(mesh%jjs,1)) :: i_loc
25  LOGICAL, DIMENSION(mesh%mes) :: virgins
26  REAL(KIND=8), DIMENSION(mesh%np) :: r_old_j_to_new
27  INTEGER, DIMENSION(5) :: opts
28  INTEGER, DIMENSION(:), ALLOCATABLE :: xadj, adjwgt, nblmt_per_proc, start, displ
29  INTEGER, DIMENSION(:,:), ALLOCATABLE :: inter
30  INTEGER :: dim, nb_neigh, edge, m, ms, n, nb, numflag, p, wgtflag, new_dof, j, &
31  news, ns, nws, msop, nsop, proc, iop
32  LOGICAL :: test=.false.
33  REAL(KIND=8), DIMENSION(mesh%me) :: b
34  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: tpwgts
35  INTEGER, DIMENSION(METIS_NOPTIONS) :: metis_opt
36  REAL(KIND=8), DIMENSION(1) :: ubvec
37 !#include "petsc/finclude/petsc.h"
38  petscerrorcode :: ierr
39  petscmpiint :: rank , nb_proc
40  mpi_comm :: communicator
41  CALL mpi_comm_rank(communicator,rank,ierr)
42  !CALL MPI_Comm_size(PETSC_COMM_WORLD,nb_proc,ierr)
43 
44  IF (nb_proc==1) THEN
45  me_loc = (/ 1, mesh%me /)
46  mes_loc = (/ 1, mesh%mes /)
47  np_loc = (/ 1, mesh%np /)
48  CALL create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc)
49  RETURN
50  END IF
51 
52  ! Create the connectivity array based on neigh
53  nb_neigh = SIZE(mesh%neigh,1)
54  xind(1) = 1
55  DO m = 1, mesh%me
56  nb = 0
57  DO n = 1, nb_neigh
58  IF (mesh%neigh(n,m)==0) cycle
59  nb = nb + 1
60  END DO
61  xind(m+1) = xind(m) + nb
62  END DO
63  ALLOCATE(xadj(xind(mesh%me+1)-1))
64  p = 0
65  DO m = 1, mesh%me
66  DO n = 1, nb_neigh
67  IF (mesh%neigh(n,m)==0) cycle
68  p = p + 1
69  xadj(p) = mesh%neigh(n,m)
70  END DO
71  END DO
72  IF (p/=xind(mesh%me+1)-1) THEN
73  CALL error_petsc('BUG, p/=xind(mesh%me+1)-1')
74  END IF
75  ! End create the connectivity array based on neigh
76 
77  ALLOCATE(adjwgt(SIZE(xadj)))
78  opts = 0
79  adjwgt = 1
80  numflag = 1 ! Fortran numbering of processors
81  wgtflag = 2
82  vwgt = 1
83  ALLOCATE(tpwgts(nb_proc))
84  tpwgts=1.d0/nb_proc
85  CALL metis_setdefaultoptions(metis_opt)
86  metis_opt(metis_option_numbering)=1
87  ubvec=1.01
88  CALL metis_partgraphrecursive(mesh%me, 1, xind,xadj,vwgt, vwgt, adjwgt, &
89  nb_proc,tpwgts , ubvec, metis_opt, edge, part)
90  IF (rank==0) THEN
91  CALL plot_const_p1_label(mesh%jj, mesh%rr, 1.d0*part, 'dd.plt')
92  END IF
93 
94  ! Count elements on processors
95  ALLOCATE(nblmt_per_proc(nb_proc),start(nb_proc),displ(nb_proc))
96  nblmt_per_proc = 0
97  DO m = 1, mesh%me
98  nblmt_per_proc(part(m)) = nblmt_per_proc(part(m)) + 1
99  END DO
100  start(1) = 0
101  DO n = 2, nb_proc
102  start(n) = start(n-1) + nblmt_per_proc(n-1)
103  END DO
104  me_loc(1) = start(rank+1) + 1
105  me_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
106  displ = start
107  ! End count elements on processors
108 
109  ! Re-order elements
110  DO m = 1, mesh%me
111  start(part(m)) = start(part(m)) + 1
112  old_m_to_new(m) = start(part(m))
113  END DO
114  ! Re-order elements
115 
116 
117  !==Search on the boundary whether ms is on a cut.
118  nws = size( mesh%jjs,1)
119  news = 0
120  inter_news = 0
121  parts = part(mesh%neighs)
122  IF (present(list_of_interfaces)) THEN !==There is an interface
123  IF (SIZE(list_of_interfaces)/=0) THEN
124  virgins = .true.
125  news = 0
126  DO ms = 1, mesh%mes
127  IF (.NOT.virgins(ms)) cycle
128  IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle !==ms not on a cut
129  i_loc = mesh%jjs(:,ms)
130  DO msop = 1, mesh%mes
131  IF (msop==ms .OR. .NOT.virgins(msop)) cycle
132  IF (minval(abs(mesh%sides(msop)-list_of_interfaces))/=0) cycle !==msop not on a cut
133  DO ns = 1, nws
134  test = .false.
135  DO nsop = 1, nws
136  iop = mesh%jjs(nsop,msop)
137  IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon) THEN
138  test = .true.
139  EXIT
140  END IF
141  END DO
142  IF (.NOT.test) THEN
143  EXIT !==This msop does not coincide with ms
144  END IF
145  END DO
146  IF (test) EXIT
147  END DO
148  IF (.NOT.test) THEN
149  CALL error_petsc(.NOT.'BUG in create_local_mesh, test ')
150  END IF
151  IF (part(mesh%neighs(ms)) == part(mesh%neighs(msop))) cycle !==ms is an internal cut
152  proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
153  parts(ms) = proc
154  parts(msop) = proc
155  virgins(ms) = .false.
156  virgins(msop) = .false.
157  IF (proc /= rank+1) cycle !==ms and msop do not touch the current proc
158  news = news + 1
159  inter_news(1,news)=ms
160  inter_news(2,news)=msop
161  END DO
162  END IF
163  END IF
164 
165 
166  ! Re-order jj
167  ALLOCATE(inter(SIZE(mesh%jj,1),mesh%me))
168  DO m = 1, mesh%me
169  inter(:,old_m_to_new(m)) = mesh%jj(:,m)
170  END DO
171  mesh%jj = inter
172 
173  virgin = .true.
174  new_dof = 0
175  DO m = 1, mesh%me
176  DO n = 1, SIZE(mesh%jj,1)
177  j = mesh%jj(n,m)
178  IF (.NOT.virgin(j)) cycle
179  new_dof = new_dof + 1
180  virgin(j) = .false.
181  old_j_to_new(j) = new_dof
182  END DO
183  END DO
184  DO m = 1, mesh%me
185  inter(:,m) = old_j_to_new(mesh%jj(:,m))
186  END DO
187  mesh%jj = inter
188  DEALLOCATE(inter)
189 
190  IF (rank == 0) THEN
191  np_loc(1) = 1
192  ELSE
193  np_loc(1) = maxval(mesh%jj(:,displ(rank)+1:displ(rank)+nblmt_per_proc(rank))) + 1
194  END IF
195  np_loc(2) = maxval(mesh%jj(:,me_loc(1):me_loc(2)))
196  ! End re-order jj
197 
198  ! Re-order rr
199  DO n = 1, SIZE(mesh%rr,1)
200  r_old_j_to_new(old_j_to_new) = mesh%rr(n,:)
201  mesh%rr(n,:) = r_old_j_to_new(:)
202  END DO
203  ! Re-order rr
204 
205  ! Re-order neigh
206  ALLOCATE(inter(SIZE(mesh%neigh,1),mesh%me))
207  dim = SIZE(mesh%rr,1)
208  DO m = 1, mesh%me
209  DO n = 1, dim+1
210  IF (mesh%neigh(n,m) /=0) THEN
211  inter(n,old_m_to_new(m)) = old_m_to_new(mesh%neigh(n,m))
212  ELSE
213  inter(n,old_m_to_new(m)) = 0
214  END IF
215  END DO
216  END DO
217  mesh%neigh = inter
218  DEALLOCATE(inter)
219  ! End re-order neigh
220 
221  ! Re-order i_d
222  DEALLOCATE(xadj); ALLOCATE(xadj(mesh%me))
223  xadj(old_m_to_new) = mesh%i_d
224  mesh%i_d=xadj
225  ! End Re-order i_d
226 
227  ! Re-order neighs
228  DEALLOCATE(xadj); ALLOCATE(xadj(mesh%mes))
229 
230 
231  nblmt_per_proc=0
232  DO ms = 1, mesh%mes
233  n = parts(ms)
234  nblmt_per_proc(n) = nblmt_per_proc(n) + 1
235  END DO
236  start(1) = 0
237  DO n = 2, nb_proc
238  start(n) = start(n-1) + nblmt_per_proc(n-1)
239  END DO
240  mes_loc(1) = start(rank+1) + 1
241  mes_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
242 
243  DO ms = 1, mesh%mes
244  n = parts(ms)
245  start(n) = start(n) + 1
246  old_ms_to_new(ms) = start(n)
247  END DO
248  xadj(old_ms_to_new) = mesh%neighs
249  mesh%neighs=xadj
250  xadj = old_m_to_new(mesh%neighs)
251  mesh%neighs=xadj
252  ! End re-order neighs
253 
254  ! Re-order inter_news
255  xadj(1:news) = old_ms_to_new(inter_news(1,1:news))
256  inter_news(1,1:news) = xadj(1:news)
257  xadj(1:news) = old_ms_to_new(inter_news(2,1:news))
258  inter_news(2,1:news) = xadj(1:news)
259  ! End re-order inter_news
260 
261  ! Re-order sides
262  xadj(old_ms_to_new) = mesh%sides
263  mesh%sides=xadj
264  ! End re-order sides
265 
266  ! Re-order jjs
267  DO n = 1, SIZE(mesh%jjs,1)
268  xadj(old_ms_to_new) = old_j_to_new(mesh%jjs(n,:))
269  mesh%jjs(n,:)=xadj
270  END DO
271  ! End re-order jjs
272 
273  !==We create the local mesh now
274  CALL create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc, news, inter_news(:,1:news))
275  CALL mpi_comm_rank(mpi_comm_world,rank,ierr)
276 
277  DEALLOCATE(xadj,adjwgt,nblmt_per_proc,start,displ)
278  END SUBROUTINE reorder_mesh
279 
280 
281  SUBROUTINE create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc, news, inter_news)
282  USE def_type_mesh
283  USE my_util
284  IMPLICIT NONE
285  TYPE(mesh_type) :: mesh, mesh_loc
286  INTEGER, DIMENSION(2), INTENT(IN) :: me_loc, mes_loc, np_loc
287  INTEGER, DIMENSION(:,:), INTENT(IN), OPTIONAL :: inter_news
288  INTEGER, OPTIONAL :: news
289  INTEGER, DIMENSION(mesh%me) :: m_glob_to_loc,m_loc_to_glob
290  INTEGER, DIMENSION(mesh%np) :: glob_to_loc,loc_to_glob
291  LOGICAL, DIMENSION(mesh%np) :: virgin
292  INTEGER :: dim, nws, nw, m, ms, mop, msop, ns, msup, minf, dof, &
293  dom_me, nwc, dom_mes, dom_np, n, i
294 
295  dim = SIZE(mesh%rr,1)
296  nws = SIZE(mesh%jjs,1)
297  nw = SIZE(mesh%jj,1)
298  nwc = SIZE(mesh%neigh,1)
299  !==Test if one proc only
300  IF (me_loc(2) - me_loc(1) + 1==mesh%me) THEN
301  mesh_loc%me = mesh%me
302  mesh_loc%np = mesh%np
303  mesh_loc%mes = mesh%mes
304  mesh_loc%dom_me = mesh%me
305  mesh_loc%dom_np = mesh%np
306  mesh_loc%dom_mes = mesh%mes
307  ALLOCATE(mesh_loc%jj(nw,mesh%me))
308  mesh_loc%jj = mesh%jj
309  ALLOCATE(mesh_loc%loc_to_glob(mesh%np))
310  DO n = 1, mesh%np
311  mesh_loc%loc_to_glob(n) = n
312  END DO
313  ALLOCATE(mesh_loc%rr(dim,mesh%np))
314  mesh_loc%rr=mesh%rr
315  ALLOCATE(mesh_loc%neigh(nwc,mesh%me))
316  mesh_loc%neigh = mesh%neigh
317  ALLOCATE(mesh_loc%i_d(mesh%me))
318  mesh_loc%i_d = mesh%i_d
319  ALLOCATE(mesh_loc%neighs(mesh_loc%mes))
320  mesh_loc%neighs = mesh%neighs
321  ALLOCATE(mesh_loc%sides(mesh_loc%mes))
322  mesh_loc%sides = mesh%sides
323  ALLOCATE(mesh_loc%jjs(nws,mesh_loc%mes))
324  mesh_loc%jjs = mesh%jjs
325  RETURN
326  END IF
327  !==End test if one proc only
328 
329  IF (.NOT.present(news) .OR. .NOT.present( inter_news)) THEN
330  CALL error_petsc(.NOT..OR..NOT.'BUG in create_local_mesh present(news) present( inter_news)')
331  END IF
332 
333  !==Create the new mesh
334  dom_me = me_loc(2) - me_loc(1) + 1
335  dom_mes = mes_loc(2) - mes_loc(1) + 1
336  dom_np = np_loc(2) - np_loc(1) + 1
337  mesh_loc%me = dom_me + news
338  mesh_loc%mes = dom_mes
339  mesh_loc%dom_me = dom_me
340  mesh_loc%dom_np = dom_np
341  mesh_loc%dom_mes = dom_mes
342 
343  !==Re-order jj
344  virgin = .true.
345  dof = 0
346  DO m = me_loc(1), me_loc(2)
347  DO n = 1, nw
348  i = mesh%jj(n,m)
349  IF(.NOT.virgin(i) .OR. i.GE.np_loc(1)) cycle
350  virgin(i) = .false.
351  dof = dof + 1
352  END DO
353  END DO
354  ALLOCATE(mesh_loc%jj(nw,mesh_loc%me))
355 
356  m_glob_to_loc = 0
357  virgin = .true.
358  dof = dom_np
359  DO m = me_loc(1), me_loc(2)
360  DO n = 1, nw
361  i = mesh%jj(n,m)
362  IF(virgin(i)) THEN
363  virgin(i) = .false.
364  IF (i .LT. np_loc(1)) THEN
365  dof = dof + 1
366  glob_to_loc(i) = dof
367  loc_to_glob(dof) = i
368  ELSE
369  glob_to_loc(i) = i-np_loc(1) + 1
370  loc_to_glob(i-np_loc(1) + 1) = i
371  END IF
372  END IF
373  END DO
374  m_loc_to_glob(m-me_loc(1)+1) = m
375  m_glob_to_loc(m) = m-me_loc(1)+1
376  END DO
377 
378  DO n = 1, nw
379  mesh_loc%jj(n,1:dom_me) = glob_to_loc(mesh%jj(n,me_loc(1):me_loc(2)))
380  END DO
381  DO ns = 1, news
382  ms=inter_news(1,ns)
383  msop=inter_news(2,ns)
384  IF (mesh%neighs(ms) < me_loc(1) .OR. mesh%neighs(ms) > me_loc(2)) THEN
385  m = mesh%neighs(ms)
386  ELSE
387  m = mesh%neighs(msop)
388  END IF
389 
390  DO n = 1, nw
391  i = mesh%jj(n,m)
392  IF(virgin(i)) THEN
393  virgin(i) = .false.
394  IF (i.GE.np_loc(1) .AND. i.LE.np_loc(2)) THEN
395  CALL error_petsc(.GE..AND..LE.'BUG in create_local_mesh, inp_loc(1) inp_loc(2)')
396  END IF
397  dof = dof + 1
398  glob_to_loc(i) = dof
399  loc_to_glob(dof) = i
400  END IF
401  END DO
402  mesh_loc%jj(:,dom_me+ns) = glob_to_loc(mesh%jj(:,m))
403  m_loc_to_glob(dom_me+ns) = m
404  m_glob_to_loc(m) = dom_me+ns
405  END DO
406  !==End re-order jj
407 
408  !==Create mesh%loc_to_glob
409  IF (maxval(mesh_loc%jj)/=dof) THEN
410  CALL error_petsc('BUG in create_local_mesh, mesh_loc%jj)/=dof')
411  END IF
412  mesh_loc%np = dof
413  ALLOCATE(mesh_loc%loc_to_glob(mesh_loc%np))
414  mesh_loc%loc_to_glob = loc_to_glob(1:mesh_loc%np)
415  !==End create mesh%loc_to_glob
416 
417  !==Re-order rr
418  ALLOCATE(mesh_loc%rr(dim,mesh_loc%np))
419  DO n = 1, mesh_loc%np
420  mesh_loc%rr(:,n) = mesh%rr(:,mesh_loc%loc_to_glob(n))
421  END DO
422  !==End re-order rr
423 
424  !==Re-order neigh
425  ALLOCATE(mesh_loc%neigh(nwc,mesh_loc%me))
426  msup = maxval(m_loc_to_glob)
427  minf = minval(m_loc_to_glob)
428  DO m = 1, mesh_loc%me
429  DO n = 1, nwc
430  mop = mesh%neigh(n,m_loc_to_glob(m))
431  IF (mop==0) THEN
432  mesh_loc%neigh(n,m) = 0
433  ELSE IF(mop<minf .OR. mop>msup) THEN
434  mesh_loc%neigh(n,m) = -mop !JLG Aug 13 2015
435  ELSE
436  mesh_loc%neigh(n,m) = m_glob_to_loc(mop)
437  END IF
438  END DO
439  END DO
440  !==End re-order neigh
441 
442  !==Re-order i_d
443  ALLOCATE(mesh_loc%i_d(mesh_loc%me))
444  mesh_loc%i_d = mesh%i_d(m_loc_to_glob(1:mesh_loc%me))
445  !==End re-order i_d
446 
447  !==Re-order neighs
448  ALLOCATE(mesh_loc%neighs(mesh_loc%mes))
449  mesh_loc%neighs = m_glob_to_loc(mesh%neighs(mes_loc(1):mes_loc(2)))
450  !==End re-order neighs
451 
452 
453  !==Re-order sides
454  ALLOCATE(mesh_loc%sides(mesh_loc%mes))
455  mesh_loc%sides = mesh%sides(mes_loc(1):mes_loc(2))
456  !==End re-order sides
457 
458  !==Re-order jjs
459  ALLOCATE(mesh_loc%jjs(nws,mesh_loc%mes))
460  DO ns = 1, nws
461  mesh_loc%jjs(ns,:) = glob_to_loc(mesh%jjs(ns,mes_loc(1):mes_loc(2)))
462  END DO
463  !==End re-order jjs
464 
465  END SUBROUTINE create_local_mesh
466 
467  SUBROUTINE free_mesh_after(mesh)
468  USE def_type_mesh
469  USE my_util
470  IMPLICIT NONE
471  TYPE(mesh_type) :: mesh
472 
473  DEALLOCATE(mesh%jj)
474  DEALLOCATE(mesh%jjs)
475  DEALLOCATE(mesh%rr)
476  DEALLOCATE(mesh%neigh)
477  DEALLOCATE(mesh%sides)
478  DEALLOCATE(mesh%neighs)
479  DEALLOCATE(mesh%i_d)
480 
481  nullify(mesh%loc_to_glob)
482  nullify(mesh%disp)
483  nullify(mesh%domnp)
484  nullify(mesh%j_s)
485 
486  IF (mesh%edge_stab) THEN
487  DEALLOCATE(mesh%iis)
488  nullify(mesh%jji)
489  DEALLOCATE(mesh%jjsi)
490  DEALLOCATE(mesh%neighi)
491  END IF
492  mesh%dom_me = 0
493  mesh%dom_np = 0
494  mesh%dom_mes = 0
495  mesh%me = 0
496  mesh%mes = 0
497  mesh%np = 0
498  mesh%nps = 0
499  mesh%mi =0
500  mesh%edge_stab = .false.
501 
502  END SUBROUTINE free_mesh_after
503 END MODULE metis_reorder_elements
subroutine loc_to_glob(mesh_loc, mesh_glob, l_t_g)
subroutine create_local_mesh(mesh, mesh_loc, me_loc, mes_loc, np_loc, news, inter_news)
subroutine, public reorder_mesh(communicator, nb_proc, mesh, mesh_loc, list_of_interfaces)
subroutine plot_const_p1_label(jj, rr, uu, file_name)
Definition: sub_plot.f90:264
subroutine, public free_mesh_after(mesh)
subroutine error_petsc(string)
Definition: my_util.f90:15