Changeset 8361

Timestamp:

Mar 18, 2012, 9:07:07 PM (13 years ago)

Author:

steve

Message:

Working on rate operators (rain) to work in parallel

Location:

trunk/anuga_core/source/anuga

Files:

• geometry/polygon.py (modified) (2 diffs)
• geometry/polygon_ext.c (modified) (1 diff)
• operators/base_operator.py (modified) (4 diffs)
• operators/rate_operators.py (modified) (9 diffs)
• structures/inlet_operator.py (modified) (8 diffs)
• structures/run_culvert.py (moved) (moved from trunk/anuga_core/source/anuga/structures/testing_culvert.py)
• structures/run_inlet_operator.py (moved) (moved from trunk/anuga_core/source/anuga/structures/testing_inlet_operator.py) (2 diffs)

trunk/anuga_core/source/anuga/geometry/polygon.py

@@ -229 +229 @@
 
 def line_intersect(triangles, line, verbose=False):
-    """Determine if a polygon and triangle overlap
+    """Determine which of a list of trianglee intersect a line
 
     """
@@ -242 +242 @@
 
     if verbose:
-        log.critical('Found %d triangles (out of %d) that overlap the polygon' % (count, M))
+        log.critical('Found %d triangles (out of %d) that intersect line' % (count, M))
 
     return indices[:count]

trunk/anuga_core/source/anuga/geometry/polygon_ext.c

@@ -514 +514 @@
     if (A >= 1 && B >= 1)
     {
-        return 1; //polygon sits completely inside a triangle
+        return 1; //line sits completely inside a triangle
     }
     

trunk/anuga_core/source/anuga/operators/base_operator.py

@@ -6 +6 @@
     """ 
 
-    def __init__(self, domain):
+    counter = 0
+
+    def __init__(self,
+                 domain,
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False):
         
         self.domain = domain
@@ -15 +22 @@
             assert self.__parallel_safe(), msg
 
-        self.logging = False
-        self.logging_file = ''
+        if description == None:
+            self.description = ' '
+        else:
+            self.description = description
+
+
+        if label == None:
+            self.label = "inlet_%g" % Operator.counter
+        else:
+            self.label = label + '_%g' % Operator.counter
+
+
+        self.verbose = verbose
+
+        # Keep count of inlet operator
+        Operator.counter += 1
+
+        self.set_logging(logging)
+
 
     def __call__(self):
@@ -37 +61 @@
         return message
     
+    def timestepping_statistics(self):
+
+        message  = 'You need to implement timestepping statistics for your operator'
+        return message
+
 
     def print_statistics(self):
 
         print self.statistics()
-
-
-    def timestepping_statistics(self):
-
-        message  = 'You need to implement timestepping statistics for your operator'
-        return message
 
     def print_timestepping_statistics(self):
@@ -58 +81 @@
         if self.logging:
             log_to_file(self.log_filename, self.timestepping_statistics())
-            
+
+
+
+    def set_logging(self, flag=True):
+
+        self.logging = flag
+
+        # If flag is true open file with mode = "w" to form a clean file for logging
+        if self.logging:
+            self.log_filename = self.label + '.log'
+            log_to_file(self.log_filename, self.statistics(), mode='w')
+            log_to_file(self.log_filename, 'time,Q')
+
+            #log_to_file(self.log_filename, self.culvert_type)
+
+

trunk/anuga_core/source/anuga/operators/rate_operators.py

@@ -14 +14 @@
 import numpy as num
 import anuga.utilities.log as log
+
+from anuga.geometry.polygon import inside_polygon
 
 from anuga.operators.base_operator import Operator
@@ -36 +38 @@
                  indices=None,
                  default_rate=None,
-                 verbose=False):
-
-        Operator.__init__(self,domain)
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False):
+
+
+        anuga.Operator.__init__(self, domain, description, label, logging, verbose)
 
         #------------------------------------------
@@ -81 +87 @@
 
     def __call__(self):
+        """
+        Apply rate to those triangles defined in indices
+
+        indices == [] don't apply anywhere
+        indices == None apply everywhere
+        otherwise applyfor the specific indices
+        """
 
         if self.indices is []:
@@ -103 +116 @@
 
     def update_rate(self, t):
-        """Virtual method allowing local modifications by writing an
-        overriding version in descendant
+        """Provide a rate to calculate added volume
         """
 
@@ -146 +158 @@
 
     def __parallel_safe(self):
-        """Operator is just applied independently to each cell and
+        """Operator is applied independently to each cell and
         so is parallel safe.
         """
@@ -166 +178 @@
 
 #===============================================================================
-# Specific Rate Operators for special regions.
+# Specific Rate Operators for circular region.
 #===============================================================================
 class Circular_rate_operator(Rate_operator):
@@ -179 +191 @@
     def __init__(self, domain,
                  rate=0.0,
-                 center = None,
-                 radius = None,
+                 center=None,
+                 radius=None,
                  default_rate=None,
                  verbose=False):
-
-        if verbose: log.critical(self.__name__+': Beginning Initialisation')
-
 
         assert center is not None
@@ -205 +214 @@
 
             if ((x-c[0])**2+(y-c[1])**2) < r**2:
-                self.indices.append(k)
+                indices.append(k)
 
 
@@ -224 +233 @@
 
 
+#===============================================================================
+# Specific Rate Operators for polygonal region.
+#===============================================================================
+class Polygonal_rate_operator(Rate_operator):
+    """
+    Add water at certain rate (ms^{-1} = vol/Area/sec) over a
+    polygonal region
+
+    rate can be a function of time.
+
+    """
+
+    def __init__(self, domain,
+                 rate=0.0,
+                 polygon=None,
+                 default_rate=None,
+                 verbose=False):
+
+        assert center is not None
+        assert radius is not None
+
+
+        # Determine indices in update region
+        N = domain.get_number_of_triangles()
+        points = domain.get_centroid_coordinates(absolute=True)
+
+
+        indices = inside_polygon(points, polygon)
+
+
+        # It is possible that circle doesn't intersect with mesh (as can happen
+        # for parallel runs
+
+
+        Rate_operator.__init__(self,
+                               domain,
+                               rate=rate,
+                               indices=indices,
+                               default_rate=default_rate,
+                               verbose=verbose)
+
+
+        self.polygon = polygon
+        
+
+
+
+

trunk/anuga_core/source/anuga/structures/inlet_operator.py

@@ -1 +1 @@
 import anuga
 import numpy
-import math
 import inlet
 
-from anuga.utilities.system_tools import log_to_file
 
-
-class Inlet_operator:
+class Inlet_operator(anuga.Operator):
     """Inlet Operator - add water to an inlet.
     Sets up the geometry of problem
@@ -17 +14 @@
     """ 
 
-    counter = 0
 
     def __init__(self,
@@ -27 +23 @@
                  logging = False,
                  verbose = False):
+
+
+        anuga.Operator.__init__(self, domain, description, label, logging, verbose)
+
         
-        self.domain = domain
-        self.domain.set_fractional_step_operator(self)
         self.line = numpy.array(line, dtype='d')
 
@@ -35 +33 @@
         self.Q = Q
 
-        if description == None:
-            self.description = ' '
-        else:
-            self.description = description
-        
-
-        if label == None:
-            self.label = "inlet_%g" % Inlet_operator.counter
-        else:
-            self.label = label + '_%g' % Inlet_operator.counter
-
-
-        self.verbose = verbose
-
-        # Keep count of inlet operator
-        Inlet_operator.counter += 1
 
         self.enquiry_point = 0.5*(self.line[0] + self.line[1])
@@ -56 +38 @@
         self.inlet = inlet.Inlet(self.domain, self.line, verbose= verbose)
 
-        self.set_logging(logging)
+        self.applied_Q = 0.0
 
     def __call__(self):
@@ -63 +45 @@
 
         t = self.domain.get_time()
+
         Q1 = self.update_Q(t)
         Q2 = self.update_Q(t + timestep)
+
+        Q = 0.5*(Q1+Q2)
+        volume = Q*timestep
         
-        volume = 0.5*(Q1+Q2)*timestep
-        
-        assert 0.5*(Q1+Q2) >= 0.0, 'Q < 0: Water to be removed from an inlet!'
-        
+        assert volume >= 0.0, 'Q < 0: Water to be removed from an inlet!'
+
+        # store last discharge
+        self.applied_Q = Q
+
         # Distribute volume so as to obtain flat surface
         self.inlet.set_stages_evenly(volume)
@@ -77 +64 @@
         
     def update_Q(self, t):
-        """Virtual method allowing local modifications by writing an
-        overriding version in descendant
+        """Allowing local modifications of Q
         """
         
@@ -121 +107 @@
 
 
-    def print_statistics(self):
-
-        print self.statistics()
-
-
-    def print_timestepping_statistics(self):
+    def timestepping_statistics(self):
 
         message = '---------------------------\n'
         message += 'Inlet report for %s:\n' % self.label
         message += '--------------------------\n'
-        message += 'Q [m^3/s]: %.2f\n' % self.Q
-        
-
-        print message
-
-
-    def set_logging(self, flag=True):
-
-        self.logging = flag
-
-        # If flag is true open file with mode = "w" to form a clean file for logging
-        if self.logging:
-            self.log_filename = self.label + '.log'
-            log_to_file(self.log_filename, self.statistics(), mode='w')
-            log_to_file(self.log_filename, 'time,Q')
-
-            #log_to_file(self.log_filename, self.culvert_type)
-
-
-    def timestepping_statistics(self):
-
-        message  = '%.5f, ' % self.domain.get_time()
-        message += '%.5f, ' % self.Q
+        message += 'Q [m^3/s]: %.2f\n' % self.applied_Q
 
         return message
-
-    def log_timestepping_statistics(self):
-
-         if self.logging:
-             log_to_file(self.log_filename, self.timestepping_statistics())
-
 
 

trunk/anuga_core/source/anuga/structures/run_inlet_operator.py

@@ -102 +102 @@
 
 
-line = [[0.0, 5.0], [0.0, 10.0]]
+line0 = [[0.0, 5.0], [0.0, 10.0]]
 
 Q = file_function('test_hydrograph.tms', quantities=['hydrograph'])
 
-Inlet_operator(domain, line, Q)
+inlet0 = Inlet_operator(domain, line0, Q, label='first inlet')
 
+
+line1 = [[1.0, 5.0], [2.0, 10.0]]
+inlet1 = Inlet_operator(domain, line1, 2.0)
 
 
@@ -150 +153 @@
 
     print domain.volumetric_balance_statistics()
-    
+
+    inlet0.print_timestepping_statistics()
+    print inlet1.timestepping_statistics()
     pass