source: anuga_core/source/anuga_parallel/distribute_mesh.py @ 7448

Last change on this file since 7448 was 7448, checked in by steve, 11 years ago

Reorganising parallel files

File size: 38.4 KB
Line 
1#########################################################
2#
3#
4#  Read in a data file and subdivide the triangle list
5#
6#
7#  The final routine, pmesh_divide_metis, does automatic
8# grid partitioning. Once testing has finished on this
9# routine the others should be removed.
10#
11#  Authors: Linda Stals and Matthew Hardy, June 2005
12#  Modified: Linda Stals, Nov 2005
13#            Jack Kelly, Nov 2005
14#            Steve Roberts, Aug 2009 (updating to numpy)
15#
16#
17#########################################################
18
19import sys
20from os import sep
21from sys import path
22from math import floor
23
24import numpy as num
25
26from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
27
28#########################################################
29#
30# If the triangles list is reordered, the quantities
31# assigned to the triangles must also be reorded.
32#
33# *) quantities contain the quantites in the old ordering
34# *) proc_sum[i] contains the number of triangles in
35# processor i
36# *) tri_index is a map from the old triangle ordering to
37# the new ordering, where the new number for triangle
38# i is proc_sum[tri_index[i][0]]+tri_index[i][1]
39#
40# -------------------------------------------------------
41#
42# *) The quantaties are returned in the new ordering
43#
44#########################################################
45
46def reorder(quantities, tri_index, proc_sum):
47
48    # Find the number triangles
49
50    N = len(tri_index)
51
52    # Temporary storage area
53
54    index = num.zeros(N, num.int)
55    q_reord = {}
56
57    # Find the new ordering of the triangles
58
59    for i in range(N):
60        bin = tri_index[i][0]
61        bin_off_set = tri_index[i][1]
62        index[i] = proc_sum[bin]+bin_off_set
63
64    # Reorder each quantity according to the new ordering
65
66    for k in quantities:
67        q_reord[k] = num.zeros((N, 3), num.float)
68        for i in range(N):
69            q_reord[k][index[i]]=quantities[k].vertex_values[i]
70    del index
71
72    return q_reord
73
74
75#########################################################
76#
77# Divide the mesh using a call to metis, through pymetis.
78#
79# -------------------------------------------------------
80#
81# *)  The nodes, triangles, boundary, and quantities are
82# returned. triangles_per_proc defines the subdivision.
83# The first triangles_per_proc[0] triangles are assigned
84# to processor 0, the next triangles_per_proc[1] are
85# assigned to processor 1 etc. The boundary and quantites
86# are ordered the same way as the triangles
87#
88#########################################################
89
90#path.append('..' + sep + 'pymetis')
91
92try:
93    from pymetis.metis import partMeshNodal
94except ImportError:
95    print "***************************************************"
96    print "         Metis is probably not compiled."
97    print "         Read \anuga_core\source\pymetis\README"
98    print "***************************************************"
99    raise ImportError
100
101def pmesh_divide_metis(domain, n_procs):
102   
103    # Initialise the lists
104    # List, indexed by processor of # triangles.
105   
106    triangles_per_proc = []
107   
108    # List of lists, indexed by processor of vertex numbers
109   
110    tri_list = []
111   
112    # List indexed by processor of cumulative total of triangles allocated.
113   
114    proc_sum = []
115    for i in range(n_procs):
116        tri_list.append([])
117        triangles_per_proc.append(0)
118        proc_sum.append([])
119
120    # Prepare variables for the metis call
121   
122    n_tri = len(domain.triangles)
123    if n_procs != 1: #Because metis chokes on it...
124        n_vert = domain.get_number_of_nodes()
125        t_list = domain.triangles.copy()
126        t_list = num.reshape(t_list, (-1,))
127   
128        # The 1 here is for triangular mesh elements.
129        edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
130        # print edgecut
131        # print npart
132        # print epart
133        del edgecut
134        del npart
135
136        # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero
137        # dimensional arrays. Correct this.
138        if type(epart[0]) == num.ndarray:
139            epart_new = num.zeros(len(epart), num.int)
140            for i in range(len(epart)):
141                epart_new[i] = epart[i][0]
142            epart = epart_new
143            del epart_new
144        # Assign triangles to processes, according to what metis told us.
145       
146        # tri_index maps triangle number -> processor, new triangle number
147        # (local to the processor)
148       
149        tri_index = {}
150        triangles = []       
151        for i in range(n_tri):
152            triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
153            tri_list[epart[i]].append(domain.triangles[i])
154            tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
155       
156        # Order the triangle list so that all of the triangles belonging
157        # to processor i are listed before those belonging to processor
158        # i+1
159
160        for i in range(n_procs):
161            for t in tri_list[i]:
162                triangles.append(t)
163           
164        # The boundary labels have to changed in accoradance with the
165        # new triangle ordering, proc_sum and tri_index help with this
166
167        proc_sum[0] = 0
168        for i in range(n_procs - 1):
169            proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]
170
171        # Relabel the boundary elements to fit in with the new triangle
172        # ordering
173
174        boundary = {}
175        for b in domain.boundary:
176            t =  tri_index[b[0]]
177            boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]
178
179        quantities = reorder(domain.quantities, tri_index, proc_sum)
180    else:
181        boundary = domain.boundary.copy()
182        triangles_per_proc[0] = n_tri
183        triangles = domain.triangles.copy()
184       
185        # This is essentially the same as a chunk of code from reorder.
186       
187        quantities = {}
188        for k in domain.quantities:
189            quantities[k] = num.zeros((n_tri, 3), num.float)
190            for i in range(n_tri):
191                quantities[k][i] = domain.quantities[k].vertex_values[i]
192       
193    # Extract the node list
194   
195    nodes = domain.get_nodes().copy()
196   
197    # Convert the triangle datastructure to be an array type,
198    # this helps with the communication
199
200    ttriangles = num.zeros((len(triangles), 3), num.int)
201    for i in range(len(triangles)):
202        ttriangles[i] = triangles[i]
203   
204    return nodes, ttriangles, boundary, triangles_per_proc, quantities
205
206
207#########################################################
208#
209# Subdivide the domain. This module is primarily
210# responsible for building the ghost layer and
211# communication pattern
212#
213#
214#  Author: Linda Stals, June 2005
215#  Modified: Linda Stals, Nov 2005 (optimise python code)
216#            Steve Roberts, Aug 2009 (convert to numpy)
217#
218#
219#########################################################
220
221
222
223#########################################################
224#
225# Subdivide the triangles into non-overlapping domains.
226#
227#  *)  The subdivision is controlled by triangles_per_proc.
228# The first triangles_per_proc[0] triangles are assigned
229# to the first processor, the second triangles_per_proc[1]
230# are assigned to the second processor etc.
231#
232#  *) nodes, triangles and boundary contains all of the
233# nodes, triangles and boundary tag information for the
234# whole domain. The triangles should be orientated in the
235# correct way and the nodes number consecutively from 0.
236#
237# -------------------------------------------------------
238#
239#  *) A dictionary containing the full_nodes, full_triangles
240# and full_boundary information for each processor is
241# returned. The node information consists of
242# [global_id, x_coord, y_coord].
243#
244#########################################################
245
246def submesh_full(nodes, triangles, boundary, triangles_per_proc):
247
248    # Initialise
249
250    tlower = 0
251    nproc = len(triangles_per_proc)
252    nnodes = len(nodes)
253    node_list = []
254    triangle_list = []
255    boundary_list = []
256    submesh = {}
257    node_range = num.reshape(num.arange(nnodes),(nnodes,1))
258
259    #print node_range
260    tsubnodes = num.concatenate((node_range, nodes), 1)
261
262
263    # Loop over processors
264
265    for p in range(nproc):
266
267        # Find triangles on processor p
268
269        tupper = triangles_per_proc[p]+tlower
270        subtriangles = triangles[tlower:tupper]
271        triangle_list.append(subtriangles)
272
273        # Find the boundary edges on processor p
274
275        subboundary = {}
276        for k in boundary:
277            if (k[0] >=tlower and k[0] < tupper):
278                subboundary[k]=boundary[k]
279        boundary_list.append(subboundary)
280
281        # Find nodes in processor p
282
283        nodemap = num.zeros(nnodes, 'i')
284        for t in subtriangles:
285            nodemap[t[0]]=1
286            nodemap[t[1]]=1
287            nodemap[t[2]]=1
288
289       
290        node_list.append(tsubnodes.take(num.flatnonzero(nodemap),axis=0))
291
292        # Move to the next processor
293
294        tlower = tupper
295
296    # Put the results in a dictionary
297
298    submesh["full_nodes"] = node_list
299    submesh["full_triangles"] = triangle_list
300    submesh["full_boundary"] = boundary_list
301
302    # Clean up before exiting
303
304    del (nodemap)
305
306    return submesh
307
308
309#########################################################
310#
311# Build the ghost layer of triangles
312#
313#  *) Given the triangle subpartion for the processor
314# build a ghost layer of triangles. The ghost layer
315# consists of two layers of neighbouring triangles.
316#
317#  *) The vertices in the ghost triangles must also
318# be added to the node list for the current processor
319#
320#
321# -------------------------------------------------------
322#
323#  *) The extra triangles and nodes are returned.
324#
325#  *)  The node information consists of
326# [global_id, x_coord, y_coord].
327#
328#  *) The triangle information consists of
329# [triangle number, t], where t = [v1, v2, v3].
330#
331#########################################################
332
333def ghost_layer(submesh, mesh, p, tupper, tlower):
334
335    ncoord = mesh.number_of_nodes
336    ntriangles = mesh.number_of_triangles
337
338    # Find the first layer of boundary triangles
339
340    trianglemap = num.zeros(ntriangles, 'i')
341    for t in range(tlower, tupper):
342       
343        n = mesh.neighbours[t, 0]
344
345        if n >= 0:
346            if n < tlower or n >= tupper:
347                trianglemap[n] = 1
348        n = mesh.neighbours[t, 1]
349        if n >= 0:
350            if n < tlower or n >= tupper:
351                trianglemap[n] = 1
352        n = mesh.neighbours[t, 2]
353        if n >= 0:
354            if n < tlower or n >= tupper:
355                trianglemap[n] = 1
356
357    # Find the second layer of boundary triangles
358
359    for t in range(len(trianglemap)):
360        if trianglemap[t]==1:
361            n = mesh.neighbours[t, 0]
362            if n >= 0:
363                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
364                    trianglemap[n] = 2
365            n = mesh.neighbours[t, 1]
366            if n >= 0:
367                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
368                    trianglemap[n] = 2
369            n = mesh.neighbours[t, 2]
370            if n >= 0:
371                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
372                    trianglemap[n] = 2
373
374    # Build the triangle list and make note of the vertices
375
376    nodemap = num.zeros(ncoord, 'i')
377    fullnodes = submesh["full_nodes"][p]
378
379    subtriangles = []
380    for i in range(len(trianglemap)):
381        if trianglemap[i] != 0:
382            t = list(mesh.triangles[i])
383            nodemap[t[0]] = 1
384            nodemap[t[1]] = 1
385            nodemap[t[2]] = 1
386
387    trilist = num.reshape(num.arange(ntriangles),(ntriangles,1))
388    tsubtriangles = num.concatenate((trilist, mesh.triangles), 1)
389    subtriangles = tsubtriangles.take(num.flatnonzero(trianglemap),axis=0)
390
391   
392    # Keep a record of the triangle vertices, if they are not already there
393
394    subnodes = []
395    for n in fullnodes:
396        nodemap[int(n[0])] = 0
397
398    nodelist = num.reshape(num.arange(ncoord),(ncoord,1))
399    tsubnodes = num.concatenate((nodelist, mesh.get_nodes()), 1)
400    subnodes = tsubnodes.take(num.flatnonzero(nodemap),axis=0)
401
402    # Clean up before exiting
403
404    del (nodelist)
405    del (trilist)
406    del (tsubnodes)
407    del (nodemap)
408    del (trianglemap)
409
410    # Return the triangles and vertices sitting on the boundary layer
411
412    return subnodes, subtriangles
413
414#########################################################
415#
416# Find the edges of the ghost trianlges that do not
417# have a neighbour in the current cell. These are
418# treated as a special type of boundary edge.
419#
420#  *) Given the ghost triangles in a particular
421# triangle, use the mesh to find its neigbours. If
422# the neighbour is not in the processor set it to
423# be a boundary edge
424#
425#  *) The vertices in the ghost triangles must also
426# be added to the node list for the current processor
427#
428#  *) The boundary edges for the ghost triangles are
429# ignored.
430#
431# -------------------------------------------------------
432#
433#  *) The type assigned to the ghost boundary edges is 'ghost'
434#
435#  *)  The boundary information is returned as a directorier
436# with the key = (triangle id, edge no) and the values
437# assigned to the key is 'ghost'
438#
439#
440#########################################################
441def is_in_processor(ghost_list, tlower, tupper, n):
442
443    return num.equal(ghost_list,n).any() or (tlower <= n and tupper > n)
444
445
446def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
447
448    ghost_list = []
449    subboundary = {}
450
451
452    for t in ghosttri:
453        ghost_list.append(t[0])
454   
455    for t in ghosttri:
456
457        n = mesh.neighbours[t[0], 0]
458        if not is_in_processor(ghost_list, tlower, tupper, n):
459            subboundary[t[0], 0] = 'ghost'
460
461        n = mesh.neighbours[t[0], 1]
462        if not is_in_processor(ghost_list, tlower, tupper, n):
463            subboundary[t[0], 1] = 'ghost'
464
465        n = mesh.neighbours[t[0], 2]
466        if not is_in_processor(ghost_list, tlower, tupper, n):
467            subboundary[t[0], 2] = 'ghost'
468           
469    return subboundary
470
471#########################################################
472#
473# The ghost triangles on the current processor will need
474# to get updated information from the neighbouring
475# processor containing the corresponding full triangles.
476#
477#  *) The tri_per_proc is used to determine which
478# processor contains the full node copy.
479#
480# -------------------------------------------------------
481#
482#  *) The ghost communication pattern consists of
483# [global node number, neighbour processor number].
484#
485#########################################################
486
487def ghost_commun_pattern(subtri, p, tri_per_proc):
488
489    # Loop over the ghost triangles
490
491    ghost_commun = num.zeros((len(subtri), 2), num.int)
492
493    for i in range(len(subtri)):
494        global_no = subtri[i][0]
495
496        # Find which processor contains the full triangle
497
498        nproc = len(tri_per_proc)
499        neigh = nproc-1
500        sum = 0
501        for q in range(nproc-1):
502            if (global_no < sum+tri_per_proc[q]):
503                neigh = q
504                break
505            sum = sum+tri_per_proc[q]
506
507        # Keep a copy of the neighbour processor number
508
509        ghost_commun[i] = [global_no, neigh]
510
511    return ghost_commun
512
513#########################################################
514#
515# The full triangles in this processor must communicate
516# updated information to neighbouring processor that
517# contain ghost triangles
518#
519#  *) The ghost communication pattern for all of the
520# processor must be built before calling this processor.
521#
522#  *) The full communication pattern is found by looping
523# through the ghost communication pattern for all of the
524# processors. Recall that this information is stored in
525# the form [global node number, neighbour processor number].
526# The full communication for the neighbour processor is
527# then updated.
528#
529# -------------------------------------------------------
530#
531#  *) The full communication pattern consists of
532# [global id, [p1, p2, ...]], where p1, p2 etc contain
533# a ghost node copy of the triangle global id.
534#
535#########################################################
536
537def full_commun_pattern(submesh, tri_per_proc):
538    tlower = 0
539    nproc = len(tri_per_proc)
540    full_commun = []
541
542    # Loop over the processor
543
544    for p in range(nproc):
545
546        # Loop over the full triangles in the current processor
547        # and build an empty dictionary
548
549        fcommun = {}
550        tupper = tri_per_proc[p]+tlower
551        for i in range(tlower, tupper):
552            fcommun[i] = []
553        full_commun.append(fcommun)
554        tlower = tupper
555
556    # Loop over the processor again
557
558    for p in range(nproc):
559
560        # Loop over the ghost triangles in the current processor,
561        # find which processor contains the corresponding full copy
562        # and note that the processor must send updates to this
563        # processor
564
565        for g in submesh["ghost_commun"][p]:
566            neigh = g[1]
567            full_commun[neigh][g[0]].append(p)
568
569    return full_commun
570
571
572#########################################################
573#
574# Given the non-overlapping grid partition, an extra layer
575# of triangles are included to help with the computations.
576# The triangles in this extra layer are not updated by
577# the processor, their updated values must be sent by the
578# processor containing the original, full, copy of the
579# triangle. The communication pattern that controls these
580# updates must also be built.
581#
582#  *) Assumes that full triangles, nodes etc have already
583# been found and stored in submesh
584#
585#  *) See the documentation for ghost_layer,
586# ghost_commun_pattern and full_commun_pattern
587#
588# -------------------------------------------------------
589#
590#  *) The additional information is added to the submesh
591# dictionary. See the documentation for ghost_layer,
592# ghost_commun_pattern and full_commun_pattern
593#
594#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
595# ghost_commun and full_commun is added to submesh
596#########################################################
597
598def submesh_ghost(submesh, mesh, triangles_per_proc):
599
600    nproc = len(triangles_per_proc)
601    tlower = 0
602    ghost_triangles = []
603    ghost_nodes = []
604    ghost_commun = []
605    ghost_bnd = []
606
607    # Loop over the processors
608
609    for p in range(nproc):
610
611        # Find the full triangles in this processor
612
613        tupper = triangles_per_proc[p]+tlower
614
615        # Build the ghost boundary layer
616
617        [subnodes, subtri] = \
618                   ghost_layer(submesh, mesh, p, tupper, tlower)
619        ghost_triangles.append(subtri)
620        ghost_nodes.append(subnodes)
621
622
623        # Find the boundary layer formed by the ghost triangles
624       
625        subbnd = ghost_bnd_layer(subtri, tlower, tupper, mesh, p)
626        ghost_bnd.append(subbnd)
627       
628        # Build the communication pattern for the ghost nodes
629
630        gcommun = \
631                ghost_commun_pattern(subtri, p, triangles_per_proc)
632        ghost_commun.append(gcommun)
633
634        # Move to the next processor
635
636        tlower = tupper
637
638
639    # Record the ghost layer and communication pattern
640
641    submesh["ghost_nodes"] = ghost_nodes
642    submesh["ghost_triangles"] = ghost_triangles
643    submesh["ghost_commun"] = ghost_commun
644    submesh["ghost_boundary"] = ghost_bnd
645   
646    # Build the communication pattern for the full triangles
647
648    full_commun = full_commun_pattern(submesh, triangles_per_proc)
649    submesh["full_commun"] = full_commun
650
651    # Return the submesh
652
653    return submesh
654
655
656#########################################################
657#
658# Certain quantities may be assigned to the triangles,
659# these quantities must be subdivided in the same way
660# as the triangles
661#
662#  *) The quantities are ordered in the same way as the
663# triangles
664#
665# -------------------------------------------------------
666#
667#  *) The quantites attached to the full triangles are
668# stored in full_quan
669#
670#  *) The quantities attached to the ghost triangles are
671# stored in ghost_quan
672#########################################################
673
674def submesh_quantities(submesh, quantities, triangles_per_proc):
675
676    nproc = len(triangles_per_proc)
677
678    lower = 0
679
680    # Build an empty dictionary to hold the quantites
681
682    submesh["full_quan"] = {}
683    submesh["ghost_quan"] = {}
684    for k in quantities:
685        submesh["full_quan"][k] = []
686        submesh["ghost_quan"][k] = []
687
688    # Loop trough the subdomains
689
690    for p in range(nproc):
691        upper =   lower+triangles_per_proc[p]
692
693        # Find the global ID of the ghost triangles
694
695        global_id = []
696        M = len(submesh["ghost_triangles"][p])
697        for j in range(M):
698            global_id.append(submesh["ghost_triangles"][p][j][0])
699
700        # Use the global ID to extract the quantites information from
701        # the full domain
702
703        for k in quantities:
704            submesh["full_quan"][k].append(quantities[k][lower:upper])
705            submesh["ghost_quan"][k].append(num.zeros( (M,3) , num.float))
706            for j in range(M):
707                submesh["ghost_quan"][k][p][j] = \
708                                               quantities[k][global_id[j]]
709
710        lower = upper
711
712    return submesh
713
714#########################################################
715#
716# Build the grid partition on the host.
717#
718#  *) See the documentation for submesh_ghost and
719# submesh_full
720#
721# -------------------------------------------------------
722#
723#  *) A dictionary containing the full_triangles,
724# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
725# ghost_boundary, ghost_commun and full_commun and true boundary polygon is returned.
726#
727#########################################################
728
729def build_submesh(nodes, triangles, edges, quantities,
730                  triangles_per_proc):
731
732    # Temporarily build the mesh to find the neighbouring
733    # triangles and true boundary polygon
734
735    mesh = Mesh(nodes, triangles)
736    boundary_polygon = mesh.get_boundary_polygon()
737   
738
739    # Subdivide into non-overlapping partitions
740
741    submeshf = submesh_full(nodes, triangles, edges, \
742                            triangles_per_proc)
743   
744    # Add any extra ghost boundary layer information
745
746    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
747
748    # Order the quantities information to be the same as the triangle
749    # information
750
751    submesh = submesh_quantities(submeshg, quantities, \
752                                 triangles_per_proc)
753
754    submesh["boundary_polygon"] = boundary_polygon
755    return submesh
756
757#########################################################
758#
759#  Given the subdivision of the grid assigned to the
760# current processor convert it into a form that is
761# appropriate for the GA datastructure.
762#
763#  The main function of these modules is to change the
764# node numbering. The GA datastructure assumes they
765# are numbered consecutively from 0.
766#
767#  The module also changes the communication pattern
768# datastructure into a form needed by parallel_advection
769#
770#  Authors: Linda Stals and Matthew Hardy, June 2005
771#  Modified: Linda Stals, Nov 2005 (optimise python code)
772#            Steve Roberts, Aug 2009 (updating to numpy)
773#
774#
775#########################################################
776
777 
778#########################################################
779# Convert the format of the data to that used by ANUGA
780#
781#
782# *) Change the nodes global ID's to an integer value,
783#starting from 0.
784#
785# *) The triangles and boundary edges must also be
786# updated accordingly.
787#
788# -------------------------------------------------------
789#
790# *) The nodes, triangles and boundary edges defined by
791# the new numbering scheme are returned
792#
793#########################################################
794
795def build_local_GA(nodes, triangles, boundaries, tri_index):
796
797    Nnodes =len(nodes)
798    Ntriangles = len(triangles)
799   
800    # Extract the nodes (using the local ID)
801   
802    GAnodes = num.take(nodes, (1, 2), 1)
803
804    # Build a global ID to local ID mapping
805
806    NGlobal = 0
807    for i in range(Nnodes):
808        if nodes[i][0] > NGlobal:
809            NGlobal = nodes[i][0]
810    index = num.zeros(int(NGlobal)+1, num.int)
811    num.put(index, num.take(nodes, (0,), 1).astype(num.int), \
812        num.arange(Nnodes))
813       
814    # Change the global IDs in the triangles to the local IDs
815
816    GAtriangles = num.zeros((Ntriangles, 3), num.int)
817    GAtriangles[:,0] = num.take(index, triangles[:,0])
818    GAtriangles[:,1] = num.take(index, triangles[:,1])
819    GAtriangles[:,2] = num.take(index, triangles[:,2])
820
821    # Change the triangle numbering in the boundaries
822
823    GAboundaries = {}
824    for b in boundaries:
825        GAboundaries[tri_index[b[0]], b[1]] = boundaries[b]
826       
827    del (index)
828   
829    return GAnodes, GAtriangles, GAboundaries
830
831
832#########################################################
833# Change the communication format to that needed by the
834# parallel advection file.
835#
836# *) The index contains [global triangle no,
837# local triangle no.]
838#
839# -------------------------------------------------------
840#
841# *) The ghost_recv and full_send dictionaries are
842# returned.
843#
844# *) ghost_recv dictionary is local id, global id, value
845#
846# *) full_recv dictionary is local id, global id, value
847#
848# *) The information is ordered by the global id. This
849# means that the communication order is predetermined and
850# local and global id do not need to be
851# compared when the information is sent/received.
852#
853#########################################################
854
855def build_local_commun(index, ghostc, fullc, nproc):
856
857    # Initialise
858
859    full_send = {}
860    ghost_recv = {}
861
862    # Build the ghost_recv dictionary (sort the
863    # information by the global numbering)
864   
865    ghostc = num.sort(ghostc, 0)
866   
867    for c in range(nproc):
868        s = ghostc[:,0]
869        d = num.compress(num.equal(ghostc[:,1],c), s)
870        if len(d) > 0:
871            ghost_recv[c] = [0, 0]
872            ghost_recv[c][1] = d
873            ghost_recv[c][0] = num.take(index, d)
874           
875    # Build a temporary copy of the full_send dictionary
876    # (this version allows the information to be stored
877    # by the global numbering)
878
879    tmp_send = {}
880    for global_id in fullc:
881        for i in range(len(fullc[global_id])):
882            neigh = fullc[global_id][i]
883            if not tmp_send.has_key(neigh):
884                tmp_send[neigh] = []
885            tmp_send[neigh].append([global_id, \
886                                    index[global_id]])
887
888    # Extract the full send information and put it in the form
889    # required for the full_send dictionary
890
891    for neigh in tmp_send:
892        neigh_commun = num.sort(tmp_send[neigh], 0)
893        full_send[neigh] = [0, 0]
894        full_send[neigh][0] = neigh_commun[:,1]
895        full_send[neigh][1] = neigh_commun[:,0]
896
897    return ghost_recv, full_send
898
899
900#########################################################
901# Convert the format of the data to that used by ANUGA
902#
903#
904# *) Change the nodes global ID's to an integer value,
905# starting from 0. The node numbering in the triangles
906# must also be updated to take this into account.
907#
908# *) The triangle number will also change, which affects
909# the boundary tag information and the communication
910# pattern.
911#
912# -------------------------------------------------------
913#
914# *) The nodes, triangles, boundary edges and communication
915# pattern defined by the new numbering scheme are returned
916#
917#########################################################
918
919def build_local_mesh(submesh, lower_t, upper_t, nproc):
920
921    # Combine the full nodes and ghost nodes
922
923    nodes = num.concatenate((submesh["full_nodes"], \
924                         submesh["ghost_nodes"]))
925   
926    # Combine the full triangles and ghost triangles
927
928    gtri =  num.take(submesh["ghost_triangles"],(1, 2, 3),1)
929    triangles = num.concatenate((submesh["full_triangles"], gtri))
930
931    # Combine the full boundaries and ghost boundaries
932
933    boundaries = submesh["full_boundary"]
934    for b in submesh["ghost_boundary"]:
935        boundaries[b]=submesh["ghost_boundary"][b]
936
937    # Make note of the new triangle numbers, including the ghost
938    # triangles
939
940    NGlobal = upper_t
941    for i in range(len(submesh["ghost_triangles"])):
942        id = submesh["ghost_triangles"][i][0]
943        if id > NGlobal:
944            NGlobal = id
945    index = num.zeros(int(NGlobal)+1, num.int)
946    index[lower_t:upper_t]=num.arange(upper_t-lower_t)
947    for i in range(len(submesh["ghost_triangles"])):
948        index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
949   
950    # Change the node numbering (and update the numbering in the
951    # triangles)
952
953    [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index)
954
955    # Extract the local quantities
956   
957    quantities ={}
958    for k in submesh["full_quan"]:
959        Nf = len(submesh["full_quan"][k])
960        Ng = len(submesh["ghost_quan"][k])
961        quantities[k] = num.zeros((Nf+Ng, 3), num.float)
962        quantities[k][0:Nf] = submesh["full_quan"][k] 
963        quantities[k][Nf:Nf+Ng] = submesh["ghost_quan"][k]
964                             
965    # Change the communication pattern into a form needed by
966    # the parallel_adv
967
968    gcommun = submesh["ghost_commun"]
969    fcommun = submesh["full_commun"]
970    [ghost_rec, full_send] = \
971                build_local_commun(index, gcommun, fcommun, nproc)
972
973    # Clean up before exiting
974
975    del(index)
976
977    return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \
978           full_send
979
980
981#########################################################
982#
983# Handle the communication between the host machine
984# (processor 0) and the processors. The host machine is
985# responsible for the doing the initial grid partitioning.
986#
987# The routines given below should be moved to the
988# build_submesh.py and build_local.py file to allow
989# overlapping of  communication and computation.
990# This should be done after more debugging.
991#
992#
993#  Author: Linda Stals, June 2005
994#  Modified: Linda Stals, Nov 2005 (optimise python code)
995#            Steve Roberts, Aug 2009 (update to numpy)
996#
997#
998#########################################################
999
1000
1001#########################################################
1002#
1003# Send the submesh to processor p.
1004#
1005# *) The order and form is strongly coupled with
1006# rec_submesh.
1007#
1008# -------------------------------------------------------
1009#
1010# *) All of the information has been sent to processor p.
1011#
1012#########################################################
1013
1014def send_submesh(submesh, triangles_per_proc, p, verbose=True):
1015
1016    import pypar
1017   
1018    myid = pypar.rank()
1019   
1020    if verbose: print 'process %d sending submesh to process %d' %(myid, p)
1021   
1022    # build and send the tagmap for the boundary conditions
1023   
1024    tagmap = {}
1025    counter = 1
1026    for b in submesh["full_boundary"][p]:
1027         bkey = submesh["full_boundary"][p][b]
1028         if not tagmap.has_key(bkey):
1029             tagmap[bkey] = counter
1030             counter = counter+1
1031    for b in submesh["ghost_boundary"][p]:
1032         bkey = submesh["ghost_boundary"][p][b]
1033         if not tagmap.has_key(bkey):
1034             tagmap[bkey] = counter
1035             counter = counter+1
1036
1037    pypar.send(tagmap, p)
1038
1039    # send the quantities key information
1040
1041    pypar.send(submesh["full_quan"].keys(), p)
1042
1043    # send the number of triangles per processor
1044
1045    pypar.send(triangles_per_proc, p)
1046
1047    # compress full_commun
1048
1049    flat_full_commun = []
1050
1051    for c in submesh["full_commun"][p]:
1052        for i in range(len(submesh["full_commun"][p][c])):
1053            flat_full_commun.append([c,submesh["full_commun"][p][c][i]])
1054
1055    # send the array sizes so memory can be allocated
1056
1057    setup_array = num.zeros((9,),num.int)
1058    setup_array[0] = len(submesh["full_nodes"][p])
1059    setup_array[1] = len(submesh["ghost_nodes"][p])
1060    setup_array[2] = len(submesh["full_triangles"][p])
1061    setup_array[3] = len(submesh["ghost_triangles"][p])
1062    setup_array[4] = len(submesh["full_boundary"][p])
1063    setup_array[5] = len(submesh["ghost_boundary"][p])
1064    setup_array[6] = len(submesh["ghost_commun"][p])
1065    setup_array[7] = len(flat_full_commun)
1066    setup_array[8] = len(submesh["full_quan"])
1067   
1068    pypar.send(num.array(setup_array, num.int), p)
1069   
1070    # send the nodes
1071
1072    pypar.send(num.array(submesh["full_nodes"][p], num.float), p)
1073    pypar.send(num.array(submesh["ghost_nodes"][p], num.float),p)
1074
1075    # send the triangles
1076
1077    pypar.send(num.array(submesh["full_triangles"][p],  num.int), p)
1078    pypar.send(num.array(submesh["ghost_triangles"][p], num.int), p)
1079
1080    # send the boundary
1081
1082    bc = []
1083    for b in submesh["full_boundary"][p]:
1084        bc.append([b[0], b[1], tagmap[submesh["full_boundary"][p][b]]])
1085
1086
1087    pypar.send(num.array(bc, num.int), p)
1088
1089    bc = []
1090    for b in submesh["ghost_boundary"][p]:
1091        bc.append([b[0], b[1], tagmap[submesh["ghost_boundary"][p][b]]])
1092
1093    pypar.send(num.array(bc, num.int), p)
1094
1095    # send the communication pattern
1096
1097    pypar.send(submesh["ghost_commun"][p], p)
1098
1099    pypar.send(num.array(flat_full_commun, num.int), p)
1100
1101    # send the quantities
1102   
1103    for k in submesh["full_quan"]:
1104        pypar.send(num.array(submesh["full_quan"][k][p], num.float), p)
1105       
1106    for k in submesh["ghost_quan"]:
1107        pypar.send(num.array(submesh["ghost_quan"][k][p], num.float),p)
1108       
1109
1110#########################################################
1111#
1112# Receive the submesh from processor p.
1113#
1114# *) The order and form is strongly coupled with
1115# send_submesh.
1116#
1117# -------------------------------------------------------
1118#
1119# *) All of the information has been received by the
1120# processor p and passed into build_local.
1121#
1122# *) The information is returned in a form needed by the
1123# GA datastructure.
1124#
1125#########################################################
1126
1127def rec_submesh_flat(p, verbose=True):
1128
1129    import pypar
1130   
1131    numproc = pypar.size()
1132    myid = pypar.rank()
1133
1134    submesh_cell = {}
1135   
1136    if verbose: print 'process %d receiving submesh from process %d' %(myid, p)
1137
1138    # receive the tagmap for the boundary conditions
1139   
1140    tagmap = pypar.receive(p)
1141
1142    itagmap = {}
1143    for t in tagmap:
1144        itagmap[tagmap[t]]=t
1145
1146    # receive the quantities key information
1147
1148    qkeys = pypar.receive(p)
1149
1150    # receive the number of triangles per processor
1151
1152    triangles_per_proc = pypar.receive(p)
1153
1154    # recieve information about the array sizes
1155
1156    setup_array = pypar.receive(p)
1157
1158    no_full_nodes      = setup_array[0]
1159    no_ghost_nodes     = setup_array[1]
1160    no_full_triangles  = setup_array[2]
1161    no_ghost_triangles = setup_array[3]
1162    no_full_boundary   = setup_array[4]
1163    no_ghost_boundary  = setup_array[5]
1164    no_ghost_commun    = setup_array[6]
1165    no_full_commun     = setup_array[7]
1166    no_quantities      = setup_array[8]
1167   
1168    # receive the full nodes
1169
1170    submesh_cell["full_nodes"] = pypar.receive(p)
1171
1172    # receive the ghost nodes
1173
1174    submesh_cell["ghost_nodes"] = pypar.receive(p)
1175   
1176    # receive the full triangles
1177
1178    submesh_cell["full_triangles"] = pypar.receive(p)
1179   
1180    # receive the ghost triangles
1181
1182    submesh_cell["ghost_triangles"] = pypar.receive(p)
1183
1184    # receive the full boundary
1185
1186    bnd_c = pypar.receive(p)
1187
1188    submesh_cell["full_boundary"] = {}
1189    for b in bnd_c:
1190        submesh_cell["full_boundary"][b[0],b[1]]=itagmap[b[2]]
1191
1192    # receive the ghost boundary
1193
1194    bnd_c = pypar.receive(p)
1195
1196    submesh_cell["ghost_boundary"] = {}
1197    for b in bnd_c:
1198        submesh_cell["ghost_boundary"][b[0],b[1]]=itagmap[b[2]]
1199
1200    # receive the ghost communication pattern
1201
1202    submesh_cell["ghost_commun"] = pypar.receive(p)
1203   
1204    # receive the full communication pattern
1205
1206    full_commun = pypar.receive(p)
1207
1208    submesh_cell["full_commun"] = {}
1209    for c in full_commun:
1210        submesh_cell["full_commun"][c[0]] = []
1211    for c in full_commun:
1212        submesh_cell["full_commun"][c[0]].append(c[1])
1213
1214    # receive the quantities
1215
1216    submesh_cell["full_quan"]={}
1217   
1218    for i in range(no_quantities):
1219        tmp = pypar.receive(p)
1220        submesh_cell["full_quan"][qkeys[i]]=num.zeros((no_full_triangles,3), num.float)
1221        submesh_cell["full_quan"][qkeys[i]][:] = tmp[:]
1222
1223    submesh_cell["ghost_quan"]={}
1224    for i in range(no_quantities):
1225        tmp = pypar.receive(p)
1226        submesh_cell["ghost_quan"][qkeys[i]]= num.zeros((no_ghost_triangles,3), num.float)
1227        submesh_cell["ghost_quan"][qkeys[i]][:] = tmp[:]
1228   
1229    return submesh_cell, triangles_per_proc,\
1230           no_full_nodes, no_full_triangles
1231
1232
1233
1234#########################################################
1235#
1236# Receive the submesh from processor p.
1237#
1238# *) The order and form is strongly coupled with
1239# send_submesh.
1240#
1241# -------------------------------------------------------
1242#
1243# *) All of the information has been received by the
1244# processor p and passed into build_local.
1245#
1246# *) The information is returned in a form needed by the
1247# GA datastructure.
1248#
1249#########################################################
1250
1251def rec_submesh(p, verbose=True):
1252
1253    import pypar
1254   
1255    numproc = pypar.size()
1256    myid = pypar.rank()
1257
1258    [submesh_cell, triangles_per_proc,\
1259     number_of_full_nodes, number_of_full_triangles] = rec_submesh_flat(p,verbose)
1260   
1261    # find the full triangles assigned to this processor
1262
1263    lower_t = 0
1264    for i in range(myid):
1265        lower_t = lower_t+triangles_per_proc[i]
1266    upper_t = lower_t+triangles_per_proc[myid]
1267
1268    # convert the information into a form needed by the GA
1269    # datastructure
1270
1271    [GAnodes, GAtriangles, boundary, quantities, ghost_rec, full_send] = \
1272              build_local_mesh(submesh_cell, lower_t, upper_t, \
1273                               numproc)
1274   
1275    return GAnodes, GAtriangles, boundary, quantities,\
1276           ghost_rec, full_send,\
1277           number_of_full_nodes, number_of_full_triangles
1278
1279
1280#########################################################
1281#
1282# Extract the submesh that will belong to the
1283# "host processor" (i.e. processor zero)
1284#
1285#  *) See the documentation for build_submesh
1286#
1287# -------------------------------------------------------
1288#
1289#  *) A dictionary containing the full_triangles,
1290# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
1291# ghost_boundary, ghost_commun and full_commun belonging
1292# to processor zero are returned.
1293#
1294#########################################################
1295def extract_hostmesh(submesh, triangles_per_proc):
1296
1297   
1298    submesh_cell = {}
1299    submesh_cell["full_nodes"] = submesh["full_nodes"][0]
1300    submesh_cell["ghost_nodes"] = submesh["ghost_nodes"][0]
1301    submesh_cell["full_triangles"] = submesh["full_triangles"][0]
1302    submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0]
1303    submesh_cell["full_boundary"] = submesh["full_boundary"][0]
1304    submesh_cell["ghost_boundary"] = submesh["ghost_boundary"][0]
1305    submesh_cell["ghost_commun"] = submesh["ghost_commun"][0]
1306    submesh_cell["full_commun"] = submesh["full_commun"][0]
1307    submesh_cell["full_quan"] ={}
1308    submesh_cell["ghost_quan"]={}
1309    for k in submesh["full_quan"]:
1310        submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
1311        submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
1312
1313    numprocs = len(triangles_per_proc)
1314    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
1315            build_local_mesh(submesh_cell, 0, triangles_per_proc[0], numprocs)
1316
1317
1318    return  points, vertices, boundary, quantities, ghost_recv_dict, \
1319           full_send_dict
1320
1321
1322
1323
Note: See TracBrowser for help on using the repository browser.