source: inundation/parallel/pmesh_divide.py @ 2769

Last change on this file since 2769 was 2769, checked in by linda, 18 years ago

Removed the old pmesh_divide routines, the mesh is now divided using metis. Commented run_parallel_advection, run_parallel_merimbula_metis and run_parallel_sw_merimbula_metis. Add run_parallel_merimbula_tes and run_parallel_sw_merimbula_test as files that can be used to test new changes etc

File size: 5.7 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#
15#
16#########################################################
17
18from os import sep
19from sys import path
20from math import floor
21
22from Numeric import zeros, Float, Int, reshape, argsort
23
24
25#########################################################
26#
27# If the triangles list is reordered, the quantities
28# assigned to the triangles must also be reorded.
29#
30# *) quantities contain the quantites in the old ordering
31# *) proc_sum[i] contains the number of triangles in
32# processor i
33# *) tri_index is a map from the old triangle ordering to
34# the new ordering, where the new number for triangle
35# i is proc_sum[tri_index[i][0]]+tri_index[i][1]
36#
37# -------------------------------------------------------
38#
39# *) The quantaties are returned in the new ordering
40#
41#########################################################
42
43def reorder(quantities, tri_index, proc_sum):
44
45    # Find the number triangles
46
47    N = len(tri_index)
48
49    # Temporary storage area
50
51    index = zeros(N, Int)
52    q_reord = {}
53
54    # Find the new ordering of the triangles
55
56    for i in range(N):
57        bin = tri_index[i][0]
58        bin_off_set = tri_index[i][1]
59        index[i] = proc_sum[bin]+bin_off_set
60
61    # Reorder each quantity according to the new ordering
62
63    for k in quantities:
64        q_reord[k] = zeros((N, 3), Float)
65        for i in range(N):
66            q_reord[k][index[i]]=quantities[k].vertex_values[i]
67    del index
68
69    return q_reord
70
71
72#########################################################
73#
74# Divide the mesh using a call to metis, through pymetis.
75#
76# -------------------------------------------------------
77#
78# *)  The nodes, triangles, boundary, and quantities are
79# returned. triangles_per_proc defines the subdivision.
80# The first triangles_per_proc[0] triangles are assigned
81# to processor 0, the next triangles_per_proc[1] are
82# assigned to processor 1 etc. The boundary and quantites
83# are ordered the same way as the triangles
84#
85#########################################################
86path.append('..' + sep + 'pymetis')
87
88from metis import partMeshNodal
89
90def pmesh_divide_metis(domain, n_procs):
91   
92    # Initialise the lists
93    # List, indexed by processor of # triangles.
94   
95    triangles_per_proc = []
96   
97    # List of lists, indexed by processor of vertex numbers
98   
99    tri_list = []
100   
101    # List indexed by processor of cumulative total of triangles allocated.
102   
103    proc_sum = []
104    for i in range(n_procs):
105        tri_list.append([])
106        triangles_per_proc.append(0)
107        proc_sum.append([])
108
109    # Prepare variables for the metis call
110   
111    n_tri = len(domain.triangles)
112    if n_procs != 1: #Because metis chokes on it...
113        n_vert = len(domain.coordinates)
114        t_list = domain.triangles.copy()
115        t_list = reshape(t_list, (-1,))
116   
117        # The 1 here is for triangular mesh elements.
118       
119        edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
120        # print edgecut
121        # print npart
122        # print epart
123        del edgecut
124        del npart
125        # Assign triangles to processes, according to what metis told us.
126       
127        # tri_index maps triangle number -> processor, new triangle number
128        # (local to the processor)
129       
130        tri_index = {}
131        triangles = []       
132        for i in range(n_tri):
133            triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
134            tri_list[epart[i]].append(domain.triangles[i])
135            tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
136       
137        # Order the triangle list so that all of the triangles belonging
138        # to processor i are listed before those belonging to processor
139        # i+1
140
141        for i in range(n_procs):
142            for t in tri_list[i]:
143                triangles.append(t)
144           
145        # The boundary labels have to changed in accoradance with the
146        # new triangle ordering, proc_sum and tri_index help with this
147
148        proc_sum[0] = 0
149        for i in range(n_procs - 1):
150            proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]
151
152        # Relabel the boundary elements to fit in with the new triangle
153        # ordering
154
155        boundary = {}
156        for b in domain.boundary:
157            t =  tri_index[b[0]]
158            boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]
159
160        quantities = reorder(domain.quantities, tri_index, proc_sum)
161    else:
162        boundary = domain.boundary.copy()
163        triangles_per_proc[0] = n_tri
164        triangles = domain.triangles.copy()
165       
166        # This is essentially the same as a chunk of code from reorder.
167       
168        quantities = {}
169        for k in domain.quantities:
170            quantities[k] = zeros((n_tri, 3), Float)
171            for i in range(n_tri):
172                quantities[k][i] = domain.quantities[k].vertex_values[i]
173       
174    # Extract the node list
175   
176    nodes = domain.coordinates.copy()
177   
178    # Convert the triangle datastructure to be an array type,
179    # this helps with the communication
180
181    ttriangles = zeros((len(triangles), 3), Int)
182    for i in range(len(triangles)):
183        ttriangles[i] = triangles[i]
184       
185    return nodes, ttriangles, boundary, triangles_per_proc, quantities
Note: See TracBrowser for help on using the repository browser.