4 REAL(KIND=8) :: epsilon = 1.d-10
6 INTEGER :: METIS_NOPTIONS=40, METIS_OPTION_NUMBERING=17
8 #include "petsc/finclude/petsc.h"
10 SUBROUTINE reorder_mesh(communicator,nb_proc,mesh,mesh_loc,list_of_interfaces)
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
38 petscerrorcode :: ierr
39 petscmpiint :: rank , nb_proc
40 mpi_comm :: communicator
41 CALL mpi_comm_rank(communicator,rank,ierr)
45 me_loc = (/ 1, mesh%me /)
46 mes_loc = (/ 1, mesh%mes /)
47 np_loc = (/ 1, mesh%np /)
53 nb_neigh =
SIZE(mesh%neigh,1)
58 IF (mesh%neigh(n,m)==0) cycle
61 xind(m+1) = xind(m) + nb
63 ALLOCATE(xadj(xind(mesh%me+1)-1))
67 IF (mesh%neigh(n,m)==0) cycle
69 xadj(p) = mesh%neigh(n,m)
72 IF (p/=xind(mesh%me+1)-1)
THEN
77 ALLOCATE(adjwgt(
SIZE(xadj)))
83 ALLOCATE(tpwgts(nb_proc))
85 CALL metis_setdefaultoptions(metis_opt)
86 metis_opt(metis_option_numbering)=1
88 CALL metis_partgraphrecursive(mesh%me, 1, xind,xadj,vwgt, vwgt, adjwgt, &
89 nb_proc,tpwgts , ubvec, metis_opt, edge, part)
95 ALLOCATE(nblmt_per_proc(nb_proc),start(nb_proc),displ(nb_proc))
98 nblmt_per_proc(part(m)) = nblmt_per_proc(part(m)) + 1
102 start(n) = start(n-1) + nblmt_per_proc(n-1)
104 me_loc(1) = start(rank+1) + 1
105 me_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
111 start(part(m)) = start(part(m)) + 1
112 old_m_to_new(m) = start(part(m))
118 nws =
size( mesh%jjs,1)
121 parts = part(mesh%neighs)
122 IF (present(list_of_interfaces))
THEN
123 IF (
SIZE(list_of_interfaces)/=0)
THEN
127 IF (.NOT.virgins(ms)) cycle
128 IF (minval(abs(mesh%sides(ms)-list_of_interfaces))/=0) cycle
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
136 iop = mesh%jjs(nsop,msop)
137 IF (maxval(abs(mesh%rr(:,i_loc(ns))-mesh%rr(:,iop))).LT.epsilon)
THEN
149 CALL
error_petsc(.NOT.
'BUG in create_local_mesh, test ')
151 IF (part(mesh%neighs(ms)) == part(mesh%neighs(msop))) cycle
152 proc = min(part(mesh%neighs(ms)),part(mesh%neighs(msop)))
155 virgins(ms) = .false.
156 virgins(msop) = .false.
157 IF (proc /= rank+1) cycle
159 inter_news(1,news)=ms
160 inter_news(2,news)=msop
167 ALLOCATE(inter(
SIZE(mesh%jj,1),mesh%me))
169 inter(:,old_m_to_new(m)) = mesh%jj(:,m)
176 DO n = 1,
SIZE(mesh%jj,1)
178 IF (.NOT.virgin(j)) cycle
179 new_dof = new_dof + 1
181 old_j_to_new(j) = new_dof
185 inter(:,m) = old_j_to_new(mesh%jj(:,m))
193 np_loc(1) = maxval(mesh%jj(:,displ(rank)+1:displ(rank)+nblmt_per_proc(rank))) + 1
195 np_loc(2) = maxval(mesh%jj(:,me_loc(1):me_loc(2)))
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(:)
206 ALLOCATE(inter(
SIZE(mesh%neigh,1),mesh%me))
207 dim =
SIZE(mesh%rr,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))
213 inter(n,old_m_to_new(m)) = 0
222 DEALLOCATE(xadj);
ALLOCATE(xadj(mesh%me))
223 xadj(old_m_to_new) = mesh%i_d
228 DEALLOCATE(xadj);
ALLOCATE(xadj(mesh%mes))
234 nblmt_per_proc(n) = nblmt_per_proc(n) + 1
238 start(n) = start(n-1) + nblmt_per_proc(n-1)
240 mes_loc(1) = start(rank+1) + 1
241 mes_loc(2) = start(rank+1) + nblmt_per_proc(rank+1)
245 start(n) = start(n) + 1
246 old_ms_to_new(ms) = start(n)
248 xadj(old_ms_to_new) = mesh%neighs
250 xadj = old_m_to_new(mesh%neighs)
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)
262 xadj(old_ms_to_new) = mesh%sides
267 DO n = 1,
SIZE(mesh%jjs,1)
268 xadj(old_ms_to_new) = old_j_to_new(mesh%jjs(n,:))
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)
277 DEALLOCATE(xadj,adjwgt,nblmt_per_proc,start,displ)
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
295 dim =
SIZE(mesh%rr,1)
296 nws =
SIZE(mesh%jjs,1)
298 nwc =
SIZE(mesh%neigh,1)
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))
311 mesh_loc%loc_to_glob(n) = n
313 ALLOCATE(mesh_loc%rr(dim,mesh%np))
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
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)')
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
346 DO m = me_loc(1), me_loc(2)
349 IF(.NOT.virgin(i) .OR. i.GE.np_loc(1)) cycle
354 ALLOCATE(mesh_loc%jj(nw,mesh_loc%me))
359 DO m = me_loc(1), me_loc(2)
364 IF (i .LT. np_loc(1))
THEN
369 glob_to_loc(i) = i-np_loc(1) + 1
374 m_loc_to_glob(m-me_loc(1)+1) = m
375 m_glob_to_loc(m) = m-me_loc(1)+1
379 mesh_loc%jj(n,1:dom_me) = glob_to_loc(mesh%jj(n,me_loc(1):me_loc(2)))
383 msop=inter_news(2,ns)
384 IF (mesh%neighs(ms) < me_loc(1) .OR. mesh%neighs(ms) > me_loc(2))
THEN
387 m = mesh%neighs(msop)
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)')
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
409 IF (maxval(mesh_loc%jj)/=dof)
THEN
410 CALL
error_petsc(
'BUG in create_local_mesh, mesh_loc%jj)/=dof')
413 ALLOCATE(mesh_loc%loc_to_glob(mesh_loc%np))
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))
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
430 mop = mesh%neigh(n,m_loc_to_glob(m))
432 mesh_loc%neigh(n,m) = 0
433 ELSE IF(mop<minf .OR. mop>msup)
THEN
434 mesh_loc%neigh(n,m) = -mop
436 mesh_loc%neigh(n,m) = m_glob_to_loc(mop)
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))
448 ALLOCATE(mesh_loc%neighs(mesh_loc%mes))
449 mesh_loc%neighs = m_glob_to_loc(mesh%neighs(mes_loc(1):mes_loc(2)))
454 ALLOCATE(mesh_loc%sides(mesh_loc%mes))
455 mesh_loc%sides = mesh%sides(mes_loc(1):mes_loc(2))
459 ALLOCATE(mesh_loc%jjs(nws,mesh_loc%mes))
461 mesh_loc%jjs(ns,:) = glob_to_loc(mesh%jjs(ns,mes_loc(1):mes_loc(2)))
476 DEALLOCATE(mesh%neigh)
477 DEALLOCATE(mesh%sides)
478 DEALLOCATE(mesh%neighs)
481 nullify(mesh%loc_to_glob)
486 IF (mesh%edge_stab)
THEN
489 DEALLOCATE(mesh%jjsi)
490 DEALLOCATE(mesh%neighi)
500 mesh%edge_stab = .false.
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)
subroutine, public free_mesh_after(mesh)
subroutine error_petsc(string)