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