Changeset 8850

Timestamp:

Apr 26, 2013, 8:33:25 AM (12 years ago)

Author:

steve

Message:

Added in possibility for spatially dependent rate in rate_operators.py

Location:

trunk/anuga_core/source/anuga/operators

Files:

• base_operator.py (modified) (1 diff)
• rate_operators.py (modified) (6 diffs)
• test_rate_operators.py (modified) (3 diffs)

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

@@ -27 +27 @@
         self.ymom_c  = self.domain.quantities['ymomentum'].centroid_values
         self.elev_c  = self.domain.quantities['elevation'].centroid_values
+        self.coord_c = self.domain.centroid_coordinates
         self.areas = self.domain.areas
 

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

@@ -57 +57 @@
         self.factor = factor
 
+        self.rate_callable = False
+        self.rate_spatial = False
+        
         self.set_rate(rate)
         self.set_default_rate(default_rate)
@@ -79 +82 @@
         indices = self.indices
 
-        rate = self.get_rate(t)
+
+        if self.rate_spatial:
+            if self.indices is None:
+                x = self.coord_c[:,0]
+                y = self.coord_c[:,1]
+            else:
+                x = self.coord_c[indices,0]
+                y = self.coord_c[indices,1]
+            rate = self.get_spatial_rate(x,y,t)
+        else:
+            rate = self.get_rate(t)
 
         if self.verbose is True:
@@ -85 +98 @@
                          % (self.quantity_name, domain.get_time(), rate))
 
-        if rate >= 0.0:
+        if num.all(rate >= 0.0):
             if self.indices is None:
                 self.stage_c[:] = self.stage_c[:]  \
@@ -107 +120 @@
             t = self.get_time()
 
-            
-        if callable(self.rate):
+
+        if  self.rate_callable:
             try:
                 rate = self.rate(t)
@@ -146 +159 @@
         return rate
 
+    def get_spatial_rate(self, x=None, y=None, t=None):
+        """Provide a rate to calculate added volume
+        only call if self.rate_spatial = True
+        """
+
+        assert self.rate_spatial
+
+        if t is None:
+            t = self.get_time()
+
+        if x is None:
+            x = self.coord_c[:,0]
+
+        if y is None:
+            y = self.coord_c[:,1]
+
+        #print xy
+        #print t
+        rate = self.rate(x,y,t)
+
+        return rate
+
 
     def set_rate(self, rate):
@@ -158 +193 @@
         assert (isinstance(rate, (int, float)) or
                 callable(rate)), msg
+
+        self.rate_callable = False
+        self.rate_spatial = False
+
         self.rate = rate
 
+        if callable(rate):
+            self.rate_callable = True
+
+            x = num.array([0.0, 1.0])
+            y = num.array([0.0, 2.0])
+            t =0.0
+
+            try:
+                self.rate(x,y,t)
+            except:
+                #print 'Problem calling with two arguments'
+                self.rate_spatial = False
+                self.rate_callable = False
+            else:
+                self.rate_spatial = True
+                self.rate_callable = True
+
+                #print self.rate_spatial , self.rate_callable
+                return
+
+            try:
+                self.rate(t)
+            except:
+                self.rate_callable = False
+                #print 'Problem calling 1 argument'
+            else:
+                self.rate_callable = True
+                self.rate_spatial = False
+
+        #print self.rate
+        #print self.rate_spatial , self.rate_callable
 
     def set_default_rate(self, default_rate):

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

@@ -1 +1 @@
 """  Test environmental forcing - rain, wind, etc.
 """
+import operator
 
 import unittest, os
@@ -188 +189 @@
                                         'Attribute2'])
 
+
         #Now try interpolation
         for i in range(20):
@@ -383 +385 @@
         assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
 
+
+    def test_rate_operator_functions_spatial(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0.0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0.0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        verbose = False
+
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+        factor = 10.0
+
+
+        def main_spatial_rate(x,y,t):
+            # x and y should be an n by 1 array
+            return x + y
+
+        default_rate = 0.0
+
+
+        operator = Rate_operator(domain, rate=main_spatial_rate, factor=factor, \
+                      indices=indices, default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        t = operator.get_time()
+        x = operator.coord_c[indices,0]
+        y = operator.coord_c[indices,1]
+        d = operator.get_timestep()*main_spatial_rate(x,y,t)*factor + 1
+
+        #print "d"
+        #print d
+        stage_ex = num.array([ 1.0,  1.0,   1.0,  1.0])
+        stage_ex[indices] = d
+
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
             
 if __name__ == "__main__":