source: trunk/anuga_core/source/anuga_parallel/parallel_inlet.py @ 9426

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

Made changes to operators to take a region as well as poly and line

File size: 16.4 KB
Line 
1# To change this template, choose Tools | Templates
2# and open the template in the editor.
3
4import anuga.geometry.polygon
5from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
6from anuga.config import velocity_protection, g
7import math
8
9import numpy as num
10from anuga.structures.inlet import Inlet
11
12
13class Parallel_Inlet(Inlet):
14    """Contains information associated with each inlet
15    """
16
17    """
18    Parallel inlet:
19   
20    master_proc - coordinates all processors associated with this inlet
21    usually the processors with domains which contains parts of this inlet.
22   
23    procs - is the list of all processors associated with this inlet.
24
25    (We assume that the above arguments are determined correctly by the parallel_operator_factory)
26    """
27
28    def __init__(self, domain, poly, master_proc = 0, procs = None, verbose=False):
29
30        self.domain = domain
31        self.poly = num.asarray(poly, dtype=num.float64)
32        self.verbose = verbose
33
34        self.line = True
35        if len(self.poly) > 2:
36            self.line = False
37
38        self.master_proc = master_proc
39
40        if procs is None:
41            self.procs = [self.master_proc]
42        else:
43            self.procs = procs
44
45        import pypar
46        self.myid = pypar.rank()
47
48        self.compute_triangle_indices()
49        self.compute_area()
50        #self.compute_inlet_length()
51
52
53
54    def compute_triangle_indices(self):
55
56
57        domain_centroids = self.domain.get_centroid_coordinates(absolute=True)
58        vertex_coordinates = self.domain.get_full_vertex_coordinates(absolute=True)
59
60
61        if self.line: # poly is a line
62            self.triangle_indices = line_intersect(vertex_coordinates, self.poly)
63
64        else: # poly is a polygon
65
66            tris_0 = line_intersect(vertex_coordinates, [self.poly[0],self.poly[1]])
67            tris_1 = inside_polygon(domain_centroids, self.poly)
68            self.triangle_indices = num.union1d(tris_0, tris_1)
69            #print self.triangle_indices
70
71        for i in self.triangle_indices:
72            assert self.domain.tri_full_flag[i] == 1
73           
74
75    def compute_area(self):
76
77        # Compute inlet area as the sum of areas of triangles identified
78        # by line. Must be called after compute_inlet_triangle_indices().
79        if len(self.triangle_indices) == 0:
80            region = 'Inlet line=%s' % (self.line)
81            msg = 'No triangles have been identified in region '
82            print "WARNING: " + msg
83
84        self.area = 0.0
85        for j in self.triangle_indices:
86            self.area += self.domain.areas[j]
87
88        msg = 'Inlet exchange area has area = %f' % self.area
89        assert self.area >= 0.0
90
91    def compute_inlet_length(self):
92        """ Compute the length of the inlet within this domain (as
93        defined by the input line
94        """
95
96        point0 = self.line[0]
97        point1 = self.line[1]
98
99        self.inlet_length = anuga.geometry.polygon.line_length(self.line)
100
101
102    def get_inlet_length(self):
103        # LOCAL
104        msg = "Warning: compute_inlet_length not implemented"
105        warnings.warn(msg)
106        return self.inlet_length
107
108    def get_line(self):
109        return self.line
110
111    def get_area(self):
112        # LOCAL
113        return self.area
114
115    def get_global_area(self):
116        # GLOBAL: Master processor gathers area from all child processors, and returns value
117
118        # WARNING: requires synchronization, must be called by all procs associated
119        # with this inlet
120
121        import pypar
122        local_area = self.area
123        area = local_area
124
125        if self.myid == self.master_proc:
126
127            for i in self.procs:
128                if i == self.master_proc: continue
129
130                val = pypar.receive(i)
131                area = area + val
132        else:
133            pypar.send(area, self.master_proc)
134
135        return area
136
137
138    def get_areas(self):
139        # Must be called after compute_inlet_triangle_indices().
140        # LOCAL
141       
142        return self.domain.areas.take(self.triangle_indices)
143
144
145    def get_stages(self):
146        # LOCAL
147
148        return self.domain.quantities['stage'].centroid_values.take(self.triangle_indices)
149
150
151    def get_average_stage(self):
152        # LOCAL
153
154        return num.sum(self.get_stages()*self.get_areas())/self.area
155
156    def get_global_average_stage(self):
157        # GLOBAL: Master processor gathers stages from all child processors, and returns average
158
159        # WARNING: requires synchronization, must be called by all procs associated
160        # with this inlet
161
162        import pypar
163        local_stage = num.sum(self.get_stages()*self.get_areas())
164        global_area = self.get_global_area()
165
166
167
168        global_stage = local_stage
169
170        if self.myid == self.master_proc:
171            for i in self.procs:
172                if i == self.master_proc: continue
173
174                val = pypar.receive(i)
175                global_stage = global_stage + val
176        else:
177            pypar.send(local_stage, self.master_proc)
178
179
180        if global_area > 0.0:
181            return global_stage/global_area
182        else:
183            return 0.0
184
185    def get_elevations(self):
186        # LOCAL
187        return self.domain.quantities['elevation'].centroid_values.take(self.triangle_indices)
188
189    def get_average_elevation(self):
190        # LOCAL
191
192        if self.area > 0:
193            return num.sum(self.get_elevations()*self.get_areas())/self.area
194        else:
195            return 0.0
196
197
198    def get_xmoms(self):
199        # LOCAL
200        return self.domain.quantities['xmomentum'].centroid_values.take(self.triangle_indices)
201
202
203    def get_average_xmom(self):
204        # LOCAL
205
206        if self.area > 0:
207            return num.sum(self.get_xmoms()*self.get_areas())/self.area
208        else:
209            return 0.0
210
211    def get_global_average_xmom(self):
212        # GLOBAL: master proc gathers all xmom values and returns average
213        # WARNING: requires synchronization, must be called by all procs associated
214        # with this inlet
215
216        import pypar
217        global_area = self.get_global_area()
218        local_xmoms = num.sum(self.get_xmoms()*self.get_areas())
219        global_xmoms = local_xmoms
220
221        if self.myid == self.master_proc:
222            for i in self.procs:
223                if i == self.master_proc: continue
224
225                val = pypar.receive(i)
226                global_xmoms = global_xmoms + val
227        else:
228            pypar.send(local_xmoms, self.master_proc)
229
230
231        if global_area > 0.0:
232            return global_xmoms/global_area
233        else:
234            return 0.0
235
236
237    def get_ymoms(self):
238        # LOCAL
239        return self.domain.quantities['ymomentum'].centroid_values.take(self.triangle_indices)
240
241
242    def get_average_ymom(self):
243        # LOCAL
244        return num.sum(self.get_ymoms()*self.get_areas())/self.area
245
246    def get_global_average_ymom(self):
247        # GLOBAL: master proc gathers all ymom values and returns average
248        # WARNING: requires synchronization, must be called by all procs associated
249        # with this inlet
250
251        import pypar
252        global_area = self.get_global_area()
253        local_ymoms = num.sum(self.get_ymoms()*self.get_areas())
254        global_ymoms = local_ymoms
255
256        if self.myid == self.master_proc:
257            for i in self.procs:
258                if i == self.master_proc: continue
259
260                val = pypar.receive(i)
261                global_ymoms = global_ymoms + val
262        else:
263            pypar.send(local_ymoms, self.master_proc)
264
265
266        if global_area > 0.0:
267            return global_ymoms/global_area
268        else:
269            return 0.0
270
271    def get_depths(self):
272        # LOCAL
273        return self.get_stages() - self.get_elevations()
274
275
276    def get_total_water_volume(self):
277        # LOCAL
278       return num.sum(self.get_depths()*self.get_areas())
279
280    def get_global_total_water_volume(self):
281        # GLOBAL: master proc gathers total water volumes from each proc and returns average
282        # WARNING: requires synchronization, must be called by all procs associated
283        # with this inlet
284
285        import pypar
286        local_volume = num.sum(self.get_depths()*self.get_areas())
287        volume = local_volume
288
289        if self.myid == self.master_proc:
290
291            for i in self.procs:
292                if i == self.master_proc: continue
293
294                val = pypar.receive(i)
295                volume = volume + val
296        else:
297            pypar.send(volume, self.master_proc)
298
299        return volume
300
301    def get_average_depth(self):
302        # LOCAL
303
304        if self.area > 0.0:
305            return self.get_total_water_volume()/self.area
306        else:
307            return 0.0
308
309    def get_global_average_depth(self):
310        # GLOBAL: master proc gathers all depth values and returns average
311        # WARNING: requires synchronization, must be called by all procs associated
312        # with this inlet
313       
314        area = self.get_global_area()
315        total_water_volume = self.get_global_total_water_volume()
316
317
318        if area > 0.0:
319            return total_water_volume / area
320        else:
321            return 0.0
322
323
324    def get_velocities(self):
325        #LOCAL
326        depths = self.get_depths()
327        u = depths*self.get_xmoms()/(depths**2 + velocity_protection)
328        v = depths*self.get_ymoms()/(depths**2 + velocity_protection)
329
330        return u, v
331
332
333    def get_xvelocities(self):
334        #LOCAL
335        depths = self.get_depths()
336        return depth*self.get_xmoms()/(depths**2 + velocity_protection)
337
338    def get_yvelocities(self):
339        #LOCAL
340        depths = self.get_depths()
341        return depths*self.get_ymoms()/(depths**2 + velocity_protection)
342
343
344    def get_average_speed(self):
345        #LOCAL
346        u, v = self.get_velocities()
347
348        average_u = num.sum(u*self.get_areas())/self.area
349        average_v = num.sum(v*self.get_areas())/self.area
350
351        return math.sqrt(average_u**2 + average_v**2)
352
353
354    def get_average_velocity_head(self):
355        #LOCAL
356        return 0.5*self.get_average_speed()**2/g
357
358
359    def get_average_total_energy(self):
360        #LOCAL
361        return self.get_average_velocity_head() + self.get_average_stage()
362
363
364    def get_average_specific_energy(self):
365        #LOCAL
366        return self.get_average_velocity_head() + self.get_average_depth()
367
368
369# Set routines (ALL LOCAL)
370
371    def set_depths(self,depth):
372
373        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + depth)
374
375
376    def set_stages(self,stage):
377
378        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
379
380
381    def set_xmoms(self,xmom):
382
383        self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
384
385
386    def set_ymoms(self,ymom):
387
388        self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
389
390
391    def set_elevations(self,elevation):
392
393        self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation)
394
395
396    def set_stages_evenly(self,volume):
397        """ Distribute volume of water over
398        inlet exchange region so that stage is level
399        """
400        # WARNING: requires synchronization, must be called by all procs associated
401        # with this inlet
402
403        import pypar
404        centroid_coordinates = self.domain.get_full_centroid_coordinates(absolute=True)
405        areas = self.get_areas()
406        stages = self.get_stages()
407
408        stages_order = stages.argsort()
409
410        # PETE: send stages and areas, apply merging procedure
411
412        s_areas = {}
413        s_stages = {}
414        s_stages_order = {}
415        total_stages = len(stages)
416
417        if self.myid == self.master_proc:
418            s_areas[self.myid] = areas
419            s_stages[self.myid] = stages
420            s_stages_order[self.myid] = stages_order
421     
422            # Recieve areas, stages, and stages order
423            for i in self.procs:
424                if i != self.master_proc:
425                    s_areas[i] = pypar.receive(i)
426                    s_stages[i] = pypar.receive(i)
427                    s_stages_order[i] = pypar.receive(i)
428                    total_stages = total_stages + len(s_stages[i])
429
430        else:
431            # Send areas, stages, and stages order to master proc of inlet
432            pypar.send(areas, self.master_proc)
433            pypar.send(stages, self.master_proc)
434            pypar.send(stages_order, self.master_proc)
435
436        # merge sorted stage order
437        if self.myid == self.master_proc:
438            pos = {}
439            summed_volume = 0.
440            summed_areas = 0.
441            prev_stage = 0.
442            num_stages = 0.
443            first = True
444           
445            for i in self.procs:
446                pos[i] = 0
447
448            while num_stages < total_stages:
449                # Determine current minimum stage of all the processors in s_stages
450                num_stages = num_stages + 1
451                current_stage = num.finfo(num.float32).max
452                index = -1
453
454                for i in self.procs:
455                    if pos[i] >= len(s_stages[i]):
456                        continue
457                       
458                    if s_stages[i][s_stages_order[i][pos[i]]] < current_stage:
459                        current_stage = s_stages[i][s_stages_order[i][pos[i]]]
460                        index = i
461
462                # If first iteration, then only update summed_areas, position, and prev|current stage
463               
464                if first:
465                    first = False
466                    summed_areas = s_areas[index][s_stages_order[index][pos[index]]]
467                    pos[index] = pos[index] + 1
468                    prev_stage = current_stage
469                    continue
470
471                assert index >= 0, "Index out of bounds"
472
473                # Update summed volume and summed areas
474                tmp_volume = summed_volume + (summed_areas * (current_stage - prev_stage))
475
476                # Terminate if volume exceeded
477                if tmp_volume >= volume:
478                    break
479
480                summed_areas = summed_areas + s_areas[index][s_stages_order[index][pos[index]]]
481                pos[index] = pos[index] + 1
482                summed_volume = tmp_volume
483               
484                # Update position of index processor and current stage
485                prev_stage = current_stage
486
487            # Calculate new stage
488            new_stage = prev_stage + (volume - summed_volume) / summed_areas
489
490            # Send postion and new stage to all processors
491            for i in self.procs:
492                if i != self.master_proc:
493                    pypar.send(pos[i], i)
494                    pypar.send(new_stage, i)
495
496            # Update own depth
497            stages[stages_order[0:pos[self.myid]]] = new_stage
498        else:
499            pos = pypar.receive(self.master_proc)
500            new_stage = pypar.receive(self.master_proc)
501            stages[stages_order[0:pos]] = new_stage
502
503        self.set_stages(stages)
504
505        stages = self.get_stages()
506        stages_order = stages.argsort()
507
508    def set_depths_evenly(self,volume):
509        """ Distribute volume over all exchange
510        cells with equal depth of water
511        """
512        new_depth = self.get_average_depth() + (volume/self.get_area())
513        self.set_depths(new_depth)
514
515    def get_master_proc(self):
516        return self.master_proc
517
518    def parallel_safe(self):
519
520        return True
521
522    def statistics(self):
523        # WARNING: requires synchronization, must be called by all procs associated
524        # with this inlet
525
526        import pypar
527
528        message = ''
529
530        tri_indices = {}
531
532        if self.myid == self.master_proc:
533            tri_indices[self.myid] = self.triangle_indices
534           
535            for proc in self.procs:
536                if proc == self.master_proc: continue
537               
538                tri_indices[proc] = pypar.receive(proc)
539
540        else:
541            pypar.send(self.triangle_indices, self.master_proc)
542
543
544        if self.myid == self.master_proc:
545            message += '=====================================\n'
546            message +=  'Inlet\n' 
547            message += '=====================================\n'
548
549            for proc in self.procs:
550                message += '======> inlet triangle indices and centres at P%d\n' %(proc)
551                message += '%s' % tri_indices[proc]
552                message += '\n'
553               
554                message += '%s' % self.domain.get_centroid_coordinates()[tri_indices[proc]]
555                message += '\n'
556               
557            message += 'line\n'
558            message += '%s' % self.line
559            message += '\n'
560
561        return message
562
563__author__="pete"
564__date__ ="$16/08/2011 6:49:42 PM$"
565
566if __name__ == "__main__":
567    print "Hello World"
Note: See TracBrowser for help on using the repository browser.