Changeset 8852

Timestamp:

Apr 26, 2013, 9:04:35 PM (12 years ago)

Author:

steve

Message:

Adding in another script to test spatial rate operator

Location:

trunk/anuga_core/source/anuga/operators

Files:

• rate_operators.py (modified) (4 diffs)
• run_rate_spatial_operator.py (added)
• test_rate_operators.py (modified) (2 diffs)

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

@@ -92 +92 @@
             rate = self.get_spatial_rate(x,y,t)
         else:
-            rate = self.get_rate(t)
+            rate = self.get_non_spatial_rate(t)
 
         if self.verbose is True:
@@ -113 +113 @@
                        + factor*rate*timestep, self.elev_c[indices])
 
-    def get_rate(self, t=None):
+
+    def get_non_spatial_rate(self, t=None):
         """Provide a rate to calculate added volume
         """
@@ -175 +176 @@
             y = self.coord_c[:,1]
 
+        #print x.shape,y.shape
+        assert isinstance(t, (int, float))
+        assert len(x) == len(y)
+
         #print xy
         #print t
@@ -274 +279 @@
     def timestepping_statistics(self):
 
-        message  = indent + self.label + ': Rate = ' + str(self.get_rate())
+        if self.rate_spatial:
+            rate = self.get_spatial_rate()
+            min_rate = num.min(rate)
+            max_rate = num.max(rate)
+            message  = indent + self.label + ': Min rate = %g, Max rate = %g '% (min_rate,max_rate)
+        else:
+            message  = indent + self.label + ': Rate = ' + str(self.get_non_spatial_rate())
+
+
         return message
 

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

@@ -420 +420 @@
 
         # Apply operator to these triangles
+        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, \
+                      default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        t = operator.get_time()
+        x = operator.coord_c[:,0]
+        y = operator.coord_c[:,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[:] = 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)
+
+    def test_rate_operator_functions_spatial_indices(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
@@ -458 +530 @@
         assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
         assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
-
             
 if __name__ == "__main__":