source: inundation/parallel/parallel_meshes.py @ 3120

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

changed run_parallel_advection_prof.py to increase the size of the grid

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