source: branches/Numeric_anuga_source/anuga_parallel/parallel_meshes.py @ 7952

Last change on this file since 7952 was 6721, checked in by steve, 16 years ago

Debugging parallel test function

File size: 10.1 KB
Line 
1"""parallel-meshes -
22D triangular domains for parallel finite-volume computations of
3the advection equation, with extra structures to define the
4sending and receiving communications define in dictionaries
5full_send_dict and ghost_recv_dict
6
7
8Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
9Geoscience Australia, 2005
10
11Modified by Linda Stals, March 2006, to include ghost boundaries
12
13"""
14
15
16import sys
17from Numeric import array, zeros, Float, Int, ones, sum
18
19import pypar
20
21from anuga.config import epsilon
22
23
24
25
26
27def parallel_rectangle(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
28
29
30    """Setup a rectangular grid of triangles
31    with m+1 by n+1 grid points
32    and side lengths len1, len2. If side lengths are omitted
33    the mesh defaults to the unit square, divided between all the
34    processors
35
36    len1: x direction (left to right)
37    len2: y direction (bottom to top)
38
39    """
40
41    processor = pypar.rank()
42    numproc   = pypar.size()
43
44    print 'numproc',numproc
45    print 'processor ',processor
46
47    m_low, m_high = pypar.balance(m_g, numproc, processor)
48   
49    n = n_g
50    m_low  = m_low-1
51    m_high = m_high+1
52
53    print 'm_low, m_high', m_low, m_high
54    m = m_high - m_low
55
56    delta1 = float(len1_g)/m_g
57    delta2 = float(len2_g)/n_g
58
59    len1 = len1_g*float(m)/float(m_g)
60    len2 = len2_g
61    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
62
63    #Calculate number of points
64    Np = (m+1)*(n+1)
65
66    class VIndex:
67
68        def __init__(self, n,m):
69            self.n = n
70            self.m = m
71
72        def __call__(self, i,j):
73            return j+i*(self.n+1)
74
75    class EIndex:
76
77        def __init__(self, n,m):
78            self.n = n
79            self.m = m
80
81        def __call__(self, i,j):
82            return 2*(j+i*self.n)
83
84
85    I = VIndex(n,m)
86    E = EIndex(n,m)
87
88    points = zeros( (Np,2), Float)
89
90    for i in range(m+1):
91        for j in range(n+1):
92
93            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
94
95    #Construct 2 triangles per rectangular element and assign tags to boundary
96    #Calculate number of triangles
97    Nt = 2*m*n
98
99
100    elements = zeros( (Nt,3), Int)
101    boundary = {}
102    Idgl = []
103    Idfl = []
104    Idgr = []
105    Idfr = []
106
107    full_send_dict = {}
108    ghost_recv_dict = {}
109    nt = -1
110    for i in range(m):
111        for j in range(n):
112
113            i1 = I(i,j+1)
114            i2 = I(i,j)
115            i3 = I(i+1,j+1)
116            i4 = I(i+1,j)
117
118            #Lower Element
119            nt = E(i,j)
120            if i == 0:
121                Idgl.append(nt)
122
123            if i == 1:
124                Idfl.append(nt)
125
126            if i == m-2:
127                Idfr.append(nt)
128
129            if i == m-1:
130                Idgr.append(nt)
131
132            if i == m-1:
133                if processor == numproc-1:
134                    boundary[nt, 2] = 'right'
135                else:
136                    boundary[nt, 2] = 'ghost'
137       
138            if j == 0:
139                boundary[nt, 1] = 'bottom'
140            elements[nt,:] = [i4,i3,i2]
141
142            #Upper Element
143            nt = E(i,j)+1
144            if i == 0:
145                Idgl.append(nt)
146
147            if i == 1:
148                Idfl.append(nt)
149
150            if i == m-2:
151                Idfr.append(nt)
152
153            if i == m-1:
154                Idgr.append(nt)
155
156            if i == 0:
157                if processor == 0:
158                    boundary[nt, 2] = 'left'
159                else:
160                    boundary[nt, 2] = 'ghost'
161            if j == n-1:
162                boundary[nt, 1] = 'top'
163            elements[nt,:] = [i1,i2,i3]
164
165    if numproc==1:
166        Idfl.extend(Idfr)
167        Idgr.extend(Idgl)
168
169        print Idfl
170        print Idgr
171       
172        Idfl = array(Idfl,Int)
173        Idgr = array(Idgr,Int)
174
175        print Idfl
176        print Idgr
177       
178        full_send_dict[processor]  = [Idfl, Idfl]
179        ghost_recv_dict[processor] = [Idgr, Idgr]
180
181        print  full_send_dict[processor]
182        print ghost_recv_dict[processor]
183    elif numproc == 2:
184        Idfl.extend(Idfr)
185        Idgr.extend(Idgl)
186        Idfl = array(Idfl,Int)
187        Idgr = array(Idgr,Int)
188        full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl]
189        ghost_recv_dict[(processor-1)%numproc] = [Idgr, Idgr]
190    else:
191        Idfl = array(Idfl,Int)
192        Idgl = array(Idgl,Int)
193
194        Idfr = array(Idfr,Int)
195        Idgr = array(Idgr,Int)
196
197        full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl]
198        ghost_recv_dict[(processor-1)%numproc] = [Idgl, Idgl]
199        full_send_dict[(processor+1)%numproc]  = [Idfr, Idfr]
200        ghost_recv_dict[(processor+1)%numproc] = [Idgr, Idgr]
201
202
203
204    print full_send_dict
205    print ghost_recv_dict       
206   
207    return  points, elements, boundary, full_send_dict, ghost_recv_dict
208
209
210
211def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
212
213
214    """Setup a rectangular grid of triangles
215    with m+1 by n+1 grid points
216    and side lengths len1, len2. If side lengths are omitted
217    the mesh defaults to the unit square.
218
219    len1: x direction (left to right)
220    len2: y direction (bottom to top)
221
222    Return to lists: points and elements suitable for creating a Mesh or
223    FVMesh object, e.g. Mesh(points, elements)
224    """
225
226    delta1 = float(len1)/m
227    delta2 = float(len2)/n
228
229    #Calculate number of points
230    Np = (m+1)*(n+1)
231
232    class VIndex:
233
234        def __init__(self, n,m):
235            self.n = n
236            self.m = m
237
238        def __call__(self, i,j):
239            return j+i*(self.n+1)
240
241    class EIndex:
242
243        def __init__(self, n,m):
244            self.n = n
245            self.m = m
246
247        def __call__(self, i,j):
248            return 2*(j+i*self.n)
249
250
251    I = VIndex(n,m)
252    E = EIndex(n,m)
253
254    points = zeros( (Np,2), Float)
255
256    for i in range(m+1):
257        for j in range(n+1):
258
259            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
260
261    #Construct 2 triangles per rectangular element and assign tags to boundary
262    #Calculate number of triangles
263    Nt = 2*m*n
264
265
266    elements = zeros( (Nt,3), Int)
267    boundary = {}
268    ghosts = {}
269    nt = -1
270    for i in range(m):
271        for j in range(n):
272
273            i1 = I(i,j+1)
274            i2 = I(i,j)
275            i3 = I(i+1,j+1)
276            i4 = I(i+1,j)
277
278            #Lower Element
279            nt = E(i,j)
280            if i == m-1:
281                ghosts[nt] = E(1,j)
282            if i == 0:
283                ghosts[nt] = E(m-2,j)
284
285            if j == n-1:
286                ghosts[nt] = E(i,1)
287
288            if j == 0:
289                ghosts[nt] = E(i,n-2)
290
291            if i == m-1:
292                if processor == numproc-1:
293                    boundary[nt, 2] = 'right'
294                else:
295                    boundary[nt, 2] = 'ghost'
296            if j == 0:
297                boundary[nt, 1] = 'bottom'
298            elements[nt,:] = [i4,i3,i2]
299
300            #Upper Element
301            nt = E(i,j)+1
302            if i == m-1:
303                ghosts[nt] = E(1,j)+1
304            if i == 0:
305                ghosts[nt] = E(m-2,j)+1
306
307            if j == n-1:
308                ghosts[nt] = E(i,1)+1
309
310            if j == 0:
311                ghosts[nt] = E(i,n-2)+1
312
313            if i == 0:
314                if processor == 0:
315                    boundary[nt, 2] = 'left'
316                else:
317                    boundary[nt, 2] = 'ghost'
318            if j == n-1:
319                boundary[nt, 1] = 'top'
320            elements[nt,:] = [i1,i2,i3]
321
322    #bottom left
323    nt = E(0,0)
324    nf = E(m-2,n-2)
325    ghosts[nt]   = nf
326    ghosts[nt+1] = nf+1
327
328    #bottom right
329    nt = E(m-1,0)
330    nf = E(1,n-2)
331    ghosts[nt]   = nf
332    ghosts[nt+1] = nf+1
333
334    #top left
335    nt = E(0,n-1)
336    nf = E(m-2,1)
337    ghosts[nt]   = nf
338    ghosts[nt+1] = nf+1
339
340    #top right
341    nt = E(m-1,n-1)
342    nf = E(1,1)
343    ghosts[nt]   = nf
344    ghosts[nt+1] = nf+1
345
346    return points, elements, boundary, ghosts
347
348def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
349
350
351    """Setup a rectangular grid of triangles
352    with m+1 by n+1 grid points
353    and side lengths len1, len2. If side lengths are omitted
354    the mesh defaults to the unit square.
355
356    len1: x direction (left to right)
357    len2: y direction (bottom to top)
358
359    Return to lists: points and elements suitable for creating a Mesh or
360    Domain object, e.g. Mesh(points, elements)
361    """
362
363    delta1 = float(len1)/m
364    delta2 = float(len2)/n
365
366    #Calculate number of points
367    Np = (m+1)*(n+1)
368
369    class VIndex:
370
371        def __init__(self, n,m):
372            self.n = n
373            self.m = m
374
375        def __call__(self, i,j):
376            return j+i*(self.n+1)
377
378    class EIndex:
379
380        def __init__(self, n,m):
381            self.n = n
382            self.m = m
383
384        def __call__(self, i,j):
385            return 2*(j+i*self.n)
386
387
388    I = VIndex(n,m)
389    E = EIndex(n,m)
390
391    points = zeros( (Np,2), Float)
392
393    for i in range(m+1):
394        for j in range(n+1):
395
396            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
397
398    #Construct 2 triangles per rectangular element and assign tags to boundary
399    #Calculate number of triangles
400    Nt = 2*m*n
401
402
403    elements = zeros( (Nt,3), Int)
404    boundary = {}
405    ghosts = {}
406    nt = -1
407    for i in range(m):
408        for j in range(n):
409
410            i1 = I(i,j+1)
411            i2 = I(i,j)
412            i3 = I(i+1,j+1)
413            i4 = I(i+1,j)
414
415            #Lower Element
416            nt = E(i,j)
417            if i == m-1:
418                ghosts[nt] = E(1,j)
419            if i == 0:
420                ghosts[nt] = E(m-2,j)
421
422            if i == m-1:
423                if processor == numproc-1:
424                    boundary[nt, 2] = 'right'
425                else:
426                    boundary[nt, 2] = 'ghost'
427            if j == 0:
428                boundary[nt, 1] = 'bottom'
429            elements[nt,:] = [i4,i3,i2]
430
431            #Upper Element
432            nt = E(i,j)+1
433            if i == m-1:
434                ghosts[nt] = E(1,j)+1
435            if i == 0:
436                ghosts[nt] = E(m-2,j)+1
437
438            if i == 0:
439                if processor == 0:
440                    boundary[nt, 2] = 'left'
441                else:
442                    boundary[nt, 2] = 'ghost'
443            if j == n-1:
444                boundary[nt, 1] = 'top'
445            elements[nt,:] = [i1,i2,i3]
446
447
448    return points, elements, boundary, ghosts
Note: See TracBrowser for help on using the repository browser.