Changeset 9673 for trunk/anuga_core/anuga

Timestamp:

Feb 13, 2015, 4:58:07 PM (10 years ago)

Author:

davies

Message:

New pumping station function + tests, and a slight update to a bridge validation report

Location:

trunk/anuga_core/anuga/structures

Files:

• internal_boundary_functions.py (modified) (2 diffs)
• tests/test_internal_boundary_functions.py (modified) (5 diffs)

trunk/anuga_core/anuga/structures/internal_boundary_functions.py

@@ -401 +401 @@
     """
 
-    def __init__(self, pump_capacity, hw_to_start_pumping, verbose=True):
-        """
-            @param pump_capacity m^3/s
-
-        """
+    def __init__(self, domain, pump_capacity, hw_to_start_pumping, hw_to_stop_pumping,
+                 initial_pump_rate=0., pump_rate_of_increase = 1.0e+100, 
+                 pump_rate_of_decrease = 1.0e+100, verbose=True):
+        """
+            @param domain ANUGA domain
+            @param pump_capacity (m^3/s)
+            @param hw_to_start_pumping Turn pumps on if hw exceeds this (m)
+            @param hw_to_stop_pumping Turn pumps off if hw below this (m)
+            @param initial_pump_rate rate of pumps at start of simulation  (m^3/s)
+            @param pump_rate_of_increase Accelleration of pump rate when turning on (m^3/s/s)
+            @param pump_rate_of_decrease Decelleration of pump rate when turning off (m^3/s/s)
+            @param verbose 
+        """
+
         self.pump_capacity = pump_capacity
         self.hw_to_start_pumping = hw_to_start_pumping
+        self.hw_to_stop_pumping = hw_to_stop_pumping
+
+        self.pump_rate_of_increase = pump_rate_of_increase
+        self.pump_rate_of_decrease = pump_rate_of_decrease
+
+        self.domain=domain
+        self.last_time_called = domain.get_time()
+        self.time = domain.get_time()
+
+        if hw_to_start_pumping < hw_to_stop_pumping:
+            raise Exception('hw_to_start_pumping should be >= hw_to_stop_pumping')
+
+        if initial_pump_rate > pump_capacity:
+            raise Exception('Initial pump rate is > pump capacity')
+        
+        if ((self.pump_rate_of_increase < 0.) | (self.pump_rate_of_decrease < 0.)):
+            raise Exception('Pump rates of increase / decrease MUST be non-negative')
+
+        if ( (pump_capacity < 0.) | (initial_pump_rate < 0.) ):
+            raise Exception('Pump rates cannot be negative')
+
+        self.pump_rate = initial_pump_rate
 
         if verbose:
@@ -419 +450 @@
     def __call__(self, hw_in, tw_in):
         """
-            @param hw_in The stage (or energy) at the headwater site
-            @param tw_in The stage (or energy) at the tailwater site
-        """
-
-        if(hw_in > self.hw_to_start_pumping):
-            Q = self.pump_capacity
+            @param hw_in The stage (or energy) at the headwater site (m)
+            @param tw_in The stage (or energy) at the tailwater site (m)
+        """
+
+        # Compute the time since last called, so we can increase / decrease the pump rate if needed
+        self.time = self.domain.get_time()
+        if self.time > self.last_time_called:
+            dt = self.time - self.last_time_called
+            self.last_time_called = self.time
         else:
-            Q = 0.
-
-        return Q
-
+            dt = 0.
+            if self.time != self.last_time_called:
+                print self.time
+                print self.last_time_called
+                print self.time - self.last_time_called
+                raise Exception('Impossible timestepping, ask Gareth')
+          
+        # Increase / decrease the pump rate if needed 
+        if hw_in < self.hw_to_stop_pumping:
+            self.pump_rate = max(0., self.pump_rate - dt*self.pump_rate_of_decrease)
+        elif hw_in > self.hw_to_start_pumping:
+            self.pump_rate = min(self.pump_capacity, self.pump_rate + dt*self.pump_rate_of_increase)
+
+        return self.pump_rate
+

trunk/anuga_core/anuga/structures/tests/test_internal_boundary_functions.py

@@ -3 +3 @@
 
 import numpy
+import anuga
 from anuga.structures import internal_boundary_functions
 from anuga.structures.internal_boundary_functions import hecras_internal_boundary_function
+from anuga.structures.internal_boundary_functions import pumping_station_function
 from anuga.utilities.system_tools import get_pathname_from_package
 
+## Data for the domain used to test the pumping station function
+boundaryPolygon = [ [0., 0.], [0., 100.], [100.0, 100.0], [100.0, 0.0]]
+wallLoc = 50.
+# The boundary polygon + riverwall breaks the mesh into multiple regions
+# Must define the resolution in these areas with an xy point + maximum area
+# Otherwise triangle.c gets confused
+regionPtAreas = [ [99., 99., 20.0*20.0*0.5],
+                  [1., 1., 20.0*20.0*0.5] ]
 
 class Test_internal_boundary_functions(unittest.TestCase):
@@ -13 +23 @@
 
         path = get_pathname_from_package('anuga.structures')
-        self.input_file = filename1 = os.path.join(path, 'tests', 'data', 'hecras_bridge_table.csv')
-        self.hb = hecras_internal_boundary_function(self.input_file, verbose=False)
+        self.input_hecras_file = os.path.join(path, 'tests', 'data', 'hecras_bridge_table.csv')
 
         return
@@ -23 +32 @@
         return
 
+    def create_domain(self, wallHeight, InitialOceanStage, InitialLandStage, 
+                      riverWall=None, riverWall_Par=None):
+        # Riverwall = list of lists, each with a set of x,y,z (and optional QFactor) values
+
+        if(riverWall is None):
+            riverWall = { 'centralWall':
+                           [ [wallLoc, 0.0, wallHeight],
+                             [wallLoc, 100.0, wallHeight]] 
+                        }
+        if(riverWall_Par is None):
+            riverWall_Par = {'centralWall':{'Qfactor':1.0}}
+        # Make the domain
+        anuga.create_mesh_from_regions(boundaryPolygon, 
+                                 boundary_tags={'left': [0],
+                                                'top': [1],
+                                                'right': [2],
+                                                'bottom': [3]},
+                                   maximum_triangle_area = 200.,
+                                   minimum_triangle_angle = 28.0,
+                                   filename = 'testRiverwall.msh',
+                                   interior_regions =[ ], #[ [higherResPolygon, 1.*1.*0.5],
+                                                          #  [midResPolygon, 3.0*3.0*0.5]],
+                                   breaklines=riverWall.values(),
+                                   use_cache=False,
+                                   verbose=False,
+                                   regionPtArea=regionPtAreas)
+
+        domain = anuga.create_domain_from_file('testRiverwall.msh')
+
+        # 05/05/2014 -- riverwalls only work with DE0 and DE1
+        domain.set_flow_algorithm('DE0')
+        domain.set_name('test_riverwall')
+
+        def topography(x,y):
+            return -x/150. 
+
+        def stagefun(x,y):
+            stg = InitialOceanStage*(x>=wallLoc) + InitialLandStage*(x<wallLoc)
+            return stg 
+
+        # NOTE: Setting quantities at centroids is important for exactness of tests
+        domain.set_quantity('elevation',topography,location='centroids')     
+        domain.set_quantity('stage', stagefun,location='centroids')            
+        
+        domain.riverwallData.create_riverwalls(riverWall,riverWall_Par,verbose=False) 
+
+        # Boundary conditions
+        Br = anuga.Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+
+        return domain
+
     def test_basic_properties(self):
 
-        assert self.hb.name == self.input_file
+        self.hb = hecras_internal_boundary_function(self.input_hecras_file, verbose=False)
+        assert self.hb.name == self.input_hecras_file
         assert len(self.hb.hw_max_given_tw) == len(self.hb.nonfree_flow_tw)
         assert len(self.hb.nonfree_flow_curves) == len(self.hb.nonfree_flow_tw)
@@ -31 +93 @@
         return
 
-    def test_function_values(self):
+    def test_hecras_internal_boundary_function_values(self):
+
+        self.hb = hecras_internal_boundary_function(self.input_hecras_file, verbose=False)
 
         # Stationary states #
@@ -90 +154 @@
         return
 
+    def test_pumping_station_function(self):
+
+        domain = self.create_domain(wallHeight=10., 
+            InitialOceanStage=4., 
+            InitialLandStage=0.)
+
+        ps_function = pumping_station_function(
+            domain=domain,
+            pump_capacity=1.0,
+            hw_to_start_pumping=0.0,
+            hw_to_stop_pumping=-1.0,
+            initial_pump_rate=0., 
+            pump_rate_of_increase = 1.0, 
+            pump_rate_of_decrease = 1.0, 
+            verbose=False)
+
+
+        # Pump should be starting at zero       
+        assert(ps_function(1., 1.) == 0.)
+    
+        # As time increases, so will the pump rate 
+        domain.time = 0.5
+        assert(numpy.allclose(ps_function(1., 1.), 0.5))
+        assert(numpy.allclose(ps_function(1., 1.), 0.5))
+        domain.time = 0.75
+        assert(numpy.allclose(ps_function(1., 1.), 0.75))
+
+        # Pump rate maximum = 1.0
+        domain.time = 1.5
+        assert(numpy.allclose(ps_function(1., 1.), 1.0))
+
+        # Force the pump rate to zero, and it should start increasing as time
+        # increases
+        ps_function.pump_rate = 0.
+        assert (ps_function(1., 1.) == 0.)
+        domain.time = 1.75
+        assert (ps_function(1., 1.) == 0.25)
+        domain.time = 2.5
+        assert (ps_function(1., 1.) == 1.0)
+        # Can't increase any more
+        domain.time = 3.0 
+        assert (ps_function(1., 1.) == 1.0)
+       
+        # Let's try to decrease the pump
+        assert (ps_function(-1.5, 1.) == 1.0)
+        domain.time = 3.5
+        assert (ps_function(-1.5, 1.) == 0.5)
+        domain.time = 3.9
+
+        assert (numpy.allclose(ps_function(-1.5, 1.), 0.1))
+        domain.time = 4.0 
+        assert (numpy.allclose(ps_function(-1.5, 1.), 0.0))
+        domain.time = 5.0 
+        assert (numpy.allclose(ps_function(-1.5, 1.), 0.0))
+
+        
+        return 
+
 
 if __name__ == "__main__":