Changeset 9673

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

Files:

• anuga_core/anuga/structures/internal_boundary_functions.py (modified) (2 diffs)
• anuga_core/anuga/structures/tests/test_internal_boundary_functions.py (modified) (5 diffs)
• anuga_core/validation_tests/behaviour_only/bridge_hecras2/results.tex (modified) (1 diff)
• anuga_work/development/gareth/tests/ras_bridge/channel_floodplain1.py (modified) (1 diff)

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__":

trunk/anuga_core/validation_tests/behaviour_only/bridge_hecras2/results.tex

@@ -61 +61 @@
 diffusion causes additional drag for overbank flows. 
 
-Several additional factors will contribute to deviations between the models. ANUGA
-models cross-channel variations in the water surface elevation, which is
+Several additional factors will contribute to deviations between the models.
+ANUGA models cross-channel variations in the water surface elevation, which is
 assumed constant in HECRAS. The under-bridge flux in ANUGA is based on the
 water elevations in the central channel up-and-down stream of the bridge, thus
-the 2D representation will have an effect on the under-bridge flow. The models
+the 2D representation will have an effect on the under-bridge flow. Further,
+since the ANUGA model here is fairly coarsely resolved, the flow details near
+the bridge are noticably affected by the choice of flow algorithm. The models
 also flux the momentum under the bridge in different ways. ANUGA's method is to
 compute the average momentum in each direction along the upstream bridge

trunk/anuga_work/development/gareth/tests/ras_bridge/channel_floodplain1.py

@@ -187 +187 @@
     hecras_discharge_function,
     exchange_lines=[bridge_in, bridge_out],
-    enquiry_gap=5.0,
+    enquiry_gap=0.01,
     zero_outflow_momentum=False,
-    smoothing_timescale=0.0,
+    smoothing_timescale=30.0,
     logging=True)