Changeset 8853

Timestamp:

Apr 27, 2013, 5:43:35 PM (12 years ago)

Author:

steve

Message:

Added some unit tests for rate operators using spatial rates

Location:

trunk/anuga_core/source/anuga/operators

Files:

• rate_operators.py (modified) (4 diffs)
• run_rate_spatial_operator.py (modified) (4 diffs)
• test_rate_operators.py (modified) (9 diffs)

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

@@ -63 +63 @@
         self.set_default_rate(default_rate)
 
+
         self.default_rate_invoked = False    # Flag
+
+        self.set_areas()
+        self.set_full_indices()
 
     def __call__(self):
@@ -88 +92 @@
                 y = self.coord_c[:,1]
             else:
-                x = self.coord_c[indices,0]
-                y = self.coord_c[indices,1]
+                x = self.coord_c[self.indices,0]
+                y = self.coord_c[self.indices,1]
             rate = self.get_spatial_rate(x,y,t)
         else:
@@ -171 +175 @@
 
         if x is None:
-            x = self.coord_c[:,0]
-
-        if y is None:
-            y = self.coord_c[:,1]
-
+            assert y is None
+            if self.indices is None:
+                x = self.coord_c[:,0]
+                y = self.coord_c[:,1]
+            else:
+                x = self.coord_c[self.indices,0]
+                y = self.coord_c[self.indices,1]
+
+        assert x is not None
+        assert y is not None
         #print x.shape,y.shape
         assert isinstance(t, (int, float))
@@ -235 +244 @@
         #print self.rate
         #print self.rate_spatial , self.rate_callable
+
+
+    def set_areas(self):
+
+        if self.indices is None:
+            self.areas = self.domain.areas
+            return
+
+        if self.indices is []:
+            self.areas = []
+            return
+
+        self.areas = self.domain.areas[self.indices]
+
+    def set_full_indices(self):
+
+        if self.indices is None:
+            self.full_indices = num.where(self.domain.tri_full_flag ==1)[0]
+            return
+
+        if self.indices is []:
+            self.full_indices = []
+            return
+
+        self.full_indices = num.where(self.domain.tri_full_flag[self.indices] == 1)[0]
+
+    def get_Q(self, full_only=True):
+        """ Calculate current overall discharge
+        """
+
+        if full_only:
+            if self.rate_spatial:
+                rate = self.get_spatial_rate() # rate is an array
+                fid = self.full_indices
+                return num.sum(self.areas[fid]*rate[fid])*self.factor
+            else:
+                rate = self.get_non_spatial_rate() # rate is a scalar
+                return num.sum(self.areas[fid]*rate)*self.factor
+        else:
+            if self.rate_spatial:
+                rate = self.get_spatial_rate() # rate is an array
+                return num.sum(self.areas*rate)*self.factor
+            else:
+                rate = self.get_non_spatial_rate() # rate is a scalar
+                return num.sum(self.areas*rate)*self.factor
 
     def set_default_rate(self, default_rate):

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

@@ -7 +7 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
+import numpy
 from anuga import rectangular_cross
 from anuga import Domain
@@ -26 +27 @@
                                                len1=length, len2=width)
 domain = Domain(points, vertices, boundary)
-domain.set_name('rate_polygon') # Output name
+domain.set_name('output_rate_spatial_operator') # Output name
 print domain.statistics()
 
@@ -123 +124 @@
 
 op1 = Polygonal_rate_operator(domain, rate=10.0, polygon=polygon2)
+
+area1 = numpy.sum(domain.areas[op1.indices])
+Q1 = 10.0*area1
+print 'op1 Q ',Q1
+
 op2 = Circular_rate_operator(domain, rate=10.0, radius=0.5, center=(10.0, 3.0))
-op3 = Rate_operator(domain, rate = lambda x,y,t: 0.01*(x+y))
 
-for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
+area2 = numpy.sum(domain.areas[op2.indices])
+Q2 = 10.0*area2
+print 'op2 Q ',Q2
+
+
+def rain(x,y,t):
+    """Function to calculate "rain"
+    input x,y should be considered to be numpy arrays
+    abd t a scalar
+    """
+    if t<10:
+        return (x+y)
+    else:
+        return 0*x
+
+
+#op3 = Rate_operator(domain, rate = rain, factor=1e-3)
+area3 = numpy.sum(domain.areas)
+Q3 = numpy.sum(op3.get_rate(t)*area3)
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+accum = 0.0
+yieldstep = 0.1
+finaltime = 40.0
+for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
     domain.print_timestepping_statistics()
     domain.print_operator_timestepping_statistics()
@@ -137 +168 @@
 
 
+    print indent + 'Exact accumultion = ', accum
+    accum += (Q1+Q2)*yieldstep
 
 

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

@@ -5 +5 @@
 import unittest, os
 import anuga
+import numpy
 from anuga import Domain
 from anuga import Reflective_boundary
@@ -403 +404 @@
         domain = Domain(points, vertices)
 
+
+        area = numpy.sum(domain.areas)
+
         #Flat surface with 1m of water
         domain.set_quantity('elevation', 0.0)
@@ -438 +442 @@
         operator()
 
+
         t = operator.get_time()
+        Q = operator.get_Q()
         x = operator.coord_c[:,0]
         y = operator.coord_c[:,1]
-        d = operator.get_timestep()*main_spatial_rate(x,y,t)*factor + 1
+        rate = main_spatial_rate(x,y,t)*factor
+        Q_ex = num.sum(domain.areas*rate)
+        d = operator.get_timestep()*rate + 1
 
         #print "d"
         #print d
+        #print area, Q, Q_ex
         stage_ex = num.array([ 1.0,  1.0,   1.0,  1.0])
         stage_ex[:] = d
@@ -457 +466 @@
         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):
+        assert num.allclose(Q_ex, Q)
+
+    def test_rate_operator_functions_spatial_with_ghost(self):
         from anuga.config import rho_a, rho_w, eta_w
         from math import pi, cos, sin
@@ -475 +485 @@
         domain = Domain(points, vertices)
 
+
+        area = numpy.sum(domain.areas)
+
         #Flat surface with 1m of water
         domain.set_quantity('elevation', 0.0)
@@ -492 +505 @@
 
         # Apply operator to these triangles
-        indices = [0,1,3]
         factor = 10.0
 
@@ -502 +514 @@
         default_rate = 0.0
 
+        # kludge to make a ghost cell
+        domain.tri_full_flag[1] = 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()
+        Q_all = operator.get_Q(full_only=False)
+        Q_full = operator.get_Q()
+        x = operator.coord_c[:,0]
+        y = operator.coord_c[:,1]
+        rate = main_spatial_rate(x,y,t)*factor
+        Q_ex_all = num.sum(domain.areas*rate)
+        Q_ex_full = num.sum(num.where(domain.tri_full_flag==1,domain.areas*rate,0.0))
+        d = operator.get_timestep()*rate + 1
+
+        #print "d"
+        #print d
+        #print Q_ex_full, Q_ex_all
+        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)
+        assert num.allclose(Q_ex_all, Q_all)
+        assert num.allclose(Q_ex_full, Q_full)
+
+    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
+
+
+        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, \
@@ -512 +608 @@
 
         t = operator.get_time()
+        Q = operator.get_Q()
         x = operator.coord_c[indices,0]
         y = operator.coord_c[indices,1]
-        d = operator.get_timestep()*main_spatial_rate(x,y,t)*factor + 1
+        rate = main_spatial_rate(x,y,t)*factor
+        Q_ex = num.sum(domain.areas[indices]*rate)
+        d = operator.get_timestep()*rate + 1
+
 
         #print "d"
@@ -530 +630 @@
         assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
         assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
-            
+        assert num.allclose(Q_ex, Q)
+
+    def test_rate_operator_functions_empty_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 = []
+        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
+
+        domain.tri_full_flag[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()
+        Q = operator.get_Q()
+        x = operator.coord_c[indices,0]
+        y = operator.coord_c[indices,1]
+        rate = main_spatial_rate(x,y,t)*factor
+        Q_ex = num.sum(domain.areas[indices]*rate)
+        d = operator.get_timestep()*rate + 1
+
+        #print Q_ex, Q
+        #print indices
+        #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)
+        assert num.allclose(Q_ex, Q)
+
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_rate_operators, 'test')