source: trunk/anuga_core/source/anuga_parallel/parallel_meshes.py @ 8026

Last change on this file since 8026 was 8026, checked in by steve, 14 years ago

fixed up pypar_available

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