Changeset 5294

Timestamp:

May 8, 2008, 10:58:15 PM (16 years ago)

Author:

ole

Message:

Work on Inflow and Rain forcing terms.
Also looked at a few backwards compatibility issues with a view to
commit new limiter options

Location:

anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (9 diffs)
• test_shallow_water_domain.py (modified) (14 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -96 +96 @@
      import Transmissive_boundary
 
-from anuga.utilities.numerical_tools import gradient, mean
+from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
 from anuga.config import minimum_storable_height
 from anuga.config import minimum_allowed_height, maximum_allowed_speed
@@ -106 +106 @@
 from anuga.config import use_edge_limiter
 from anuga.config import use_centroid_velocities
+
+
+from anuga.utilities.polygon import inside_polygon, polygon_area        
+
 
 
@@ -1479 +1483 @@
 
 
-class Rainfall:
+
+
+class General_forcing:
+    """Class General_forcing - general explicit forcing term for update of quantity
+    
+    This is used by Inflow and Rainfall for instance
+    
+
+    General_forcing(quantity_name, rate, center, radius, polygon)
+
+    domain:     ANUGA computational domain
+    quantity_name: Name of quantity to update. It must be a known conserved quantity.
+    rate [?/s]: Total rate of change over the specified area.  
+                This parameter can be either a constant or a
+                function of time. Positive values indicate increases, 
+                negative values indicate decreases.
+                Rate can be None at initialisation but must be specified
+                before forcting term is applied (i.e. simulation has started).
+
+    center [m]: Coordinates at center of flow point
+    radius [m]: Size of circular area
+    polygon:    Arbitrary polygon.
+
+
+    Either center, radius or polygon can be specified but not both.
+    If neither are specified the entire domain gets updated.
+    
+    See Inflow or Rainfall for examples of use
+    """
+
+
+    # FIXME (AnyOne) : Add various methods to allow spatial variations
+
+    def __init__(self,
+                 domain,
+                 quantity_name,
+                 rate=None,
+                 center=None, radius=None,
+                 polygon=None,
+                 verbose=False):
+                     
+
+        from math import pi
+
+        self.domain = domain
+        self.quantity_name = quantity_name
+        self.rate = rate
+        self.center = ensure_numeric(center)
+        self.radius = radius
+        self.polygon = polygon        
+        self.verbose = verbose
+
+        # Update area if applicable
+        self.area = None        
+        if center is not None and radius is not None:
+            assert len(center) == 2
+            msg = 'Polygon cannot be specified when center and radius are'
+            assert polygon is None, msg
+
+            self.area = radius**2*pi
+        
+        if polygon is not None:
+            self.area = polygon_area(self.polygon)
+
+
+        # Pointer to update vector
+        self.update = domain.quantities[self.quantity_name].explicit_update            
+
+        # Determine indices in flow area
+        N = len(domain)    
+        points = domain.get_centroid_coordinates(absolute=True)
+
+        self.indices = None
+        if self.center is not None and self.radius is not None:
+            # Inlet is circular
+            
+            self.indices = []
+            for k in range(N):
+                x, y = points[k,:] # Centroid
+                if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
+                    self.indices.append(k)
+                    
+        if self.polygon is not None:                    
+            # Inlet is polygon
+            self.indices = inside_polygon(points, self.polygon)
+            
+
+            
+
+
+    def __call__(self, domain):
+        """Apply inflow function at time specified in domain and update stage
+        """
+
+        # Call virtual method allowing local modifications
+        rate = self.update_rate(domain.get_time())
+        if rate is None:
+            msg = 'Attribute rate must be specified in General_forcing'
+            msg += ' or its descendants before attempting to call it'
+            raise Exception, msg
+        
+
+        # Now rate is a number
+        if self.verbose is True:
+            print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
+                                                     domain.get_time(),
+                                                     rate)
+
+
+        if self.indices is None:
+            self.update[:] += rate
+        else:
+            # Brute force assignment of restricted rate
+            for k in self.indices:
+                self.update[k] += rate
+
+
+    def update_rate(self, t):
+        """Virtual method allowing local modifications by writing an
+        overriding version in descendant
+        
+        """
+        if callable(self.rate):
+            rate = self.rate(t)
+        else:
+            rate = self.rate
+
+        return rate
+
+
+    def get_quantity_values(self):
+        """Return values for specified quantity restricted to opening 
+        """
+        return self.domain.quantities[self.quantity_name].get_values(indices=self.indices)
+    
+
+    def set_quantity_values(self, val):
+        """Set values for specified quantity restricted to opening 
+        """
+        self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices)    
+
+
+
+class Rainfall(General_forcing):
     """Class Rainfall - general 'rain over entire domain' forcing term.
     
@@ -1490 +1637 @@
     
     Rainfall(rain)
-        
+
+    domain    
     rain [mm/s]:  Total rain rate over the specified domain.  
                   NOTE: Raingauge Data needs to reflect the time step.
@@ -1504 +1652 @@
                   The specified flow will be divided by the area of
                   the inflow region and then applied to update the
-                  quantity in question. 
+                  stage quantity.
 
     polygon: Specifies a polygon to restrict the rainfall.
@@ -1524 +1672 @@
     """
 
-    # FIXME (OLE): Add a polygon as an alternative.
-    # FIXME (AnyOne) : Add various methods to allow spatial variations
-    # FIXME (OLE): Generalise to all quantities
-
     
     def __init__(self,
-                 rain=0.0,
-                 quantity_name='stage',
-                 polygon=None):
-
-        self.rain = rain
-        self.quantity_name = quantity_name
-        self.polygon = polygon
-    
-    def __call__(self, domain):
-        
-        # FIXME(Ole): Move this to top of file eventually
-        from anuga.utilities.polygon import inside_polygon        
-
-        # Update rainfall
-        if callable(self.rain):
-            rain = self.rain(domain.get_time())
-        else:
-            rain = self.rain
-
-        # Now rain is a number
-        # Converting mm/s to m/s to apply in ANUGA
-        # 1mm of rainfall is equivalent to 1 litre /m2 
-        # Flow is expressed as m3/s converted to a stage height in (m)
-        
-        # Note 1m3 = 1x10^9mm3 (mls)
-        # or is that m3 to Litres ??? Check this how is it applied !!!
-        rain_fall = rain/1000
-
-        # Find indices for coordinates inside polygon if specified
-        indices = None
-        if self.polygon is not None:
-            points = domain.get_centroid_coordinates(absolute=True)
-            indices = inside_polygon(points, self.polygon)
-            
-
-        # Assign rainfall value to explicit update vector to be used as
-        # a forcing term
-        quantity = domain.quantities[self.quantity_name].explicit_update
-        if indices is None:
-            quantity[:] += rain_fall
+                 domain,
+                 rate=0.0,
+                 center=None, radius=None,
+                 polygon=None,
+                 verbose=False):
+
+        # Converting mm/s to m/s to apply in ANUGA)
+        if callable(rate):
+            rain = lambda t: rate(t)/1000.0
         else:
-            # Brute force assignment of restricted rainfall
-            for i in indices:
-                quantity[i] += rain_fall
-
-        
-
-
-
-
-class Inflow:
+            rain = rate/1000.0            
+            
+        General_forcing.__init__(self,
+                                 domain,
+                                 'stage',
+                                 rate=rain,
+                                 center=center, radius=radius,
+                                 polygon=polygon,
+                                 verbose=verbose)
+
+        
+
+
+
+
+class Inflow(General_forcing):
     """Class Inflow - general 'rain and drain' forcing term.
     
     Useful for implementing flows in and out of the domain.
     
-    Inflow(center, radius, flow)
-    
-    center [m]: Coordinates at center of flow point
-    radius [m]: Size of circular area
+    Inflow(flow, center, radius, polygon)
+
+    domain
     flow [m^3/s]: Total flow rate over the specified area.  
                   This parameter can be either a constant or a
@@ -1595 +1713 @@
                   The specified flow will be divided by the area of
                   the inflow region and then applied to update the
-                  quantity in question. 
+                  quantity in question.     
+    center [m]: Coordinates at center of flow point
+    radius [m]: Size of circular area
+    polygon:    Arbitrary polygon.
+
+    Either center, radius or polygon must be specified
     
     Examples
@@ -1612 +1735 @@
     Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
 
+
     #--------------------------------------------------------------------------
     # Setup specialised forcing terms
@@ -1625 +1749 @@
     """
 
-    # FIXME (OLE): Add a polygon as an alternative.
-    # FIXME (OLE): Generalise to all quantities
 
     def __init__(self,
+                 domain,
+                 rate=0.0,
                  center=None, radius=None,
-                 flow=0.0,
-                 quantity_name='stage',
-                 verbose=False):
-
-        from math import pi
-
-        
-                 
-        if center is not None and radius is not None:
-            assert len(center) == 2
+                 polygon=None,
+                 verbose=False):                 
+
+
+        #msg = 'Class Inflow must have either center & radius or a polygon specified.'
+        #assert center is not None and radius is not None or\
+        #       polygon is not None, msg
+
+
+        # Create object first to make area is available
+        General_forcing.__init__(self,
+                                 domain,
+                                 'stage',
+                                 rate=rate,
+                                 center=center, radius=radius,
+                                 polygon=polygon,
+                                 verbose=verbose)
+
+    def update_rate(self, t):
+        """Virtual method allowing local modifications by writing an
+        overriding version in descendant
+
+        This one converts m^3/s to m/s which can be added directly to 'stage' in ANUGA
+        """
+
+        
+        
+        if callable(self.rate):
+            _rate = self.rate(t)/self.area
         else:
-            msg = 'Both center and radius must be specified'
-            raise Exception, msg
-    
-        self.center = center
-        self.radius = radius
-        self.area = radius**2*pi
-        self.flow = flow
-        self.quantity_name = quantity_name
-        self.verbose = verbose
-    
-    def __call__(self, domain):
-
-        # Determine indices in flow area
-        if not hasattr(self, 'indices'):
-            center = self.center
-            radius = self.radius
-            
-            N = len(domain)    
-            self.indices = []
-            coordinates = domain.get_centroid_coordinates()     
-            for k in range(N):
-                x, y = coordinates[k,:] # Centroid
-                if ((x-center[0])**2+(y-center[1])**2) < radius**2:
-                    self.indices.append(k)    
-
-        # Update inflow
-        if callable(self.flow):
-            flow = self.flow(domain.get_time())
-        else:
-            flow = self.flow
-
-        # Now flow is a number
-        if self.verbose is True:
-            print 'Inflow at time = %.2f = %f' %(domain.get_time(),
-                                                 flow)
-        
-        quantity = domain.quantities[self.quantity_name].explicit_update
-        for k in self.indices:
-            quantity[k] += flow/self.area                       
+            _rate = self.rate/self.area
+
+        return _rate
+
+
 
 

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -2194 +2194 @@
         # Setup only one forcing term, constant rainfall
         domain.forcing_terms = []
-        domain.forcing_terms.append( Rainfall(rain=2.0) )
+        domain.forcing_terms.append( Rainfall(domain, rate=2.0) )
 
         domain.compute_forcing_terms()
@@ -2231 +2231 @@
         # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        domain.forcing_terms.append( Rainfall(rain=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
+        domain.forcing_terms.append( Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
 
         domain.compute_forcing_terms()
@@ -2242 +2242 @@
         # FIXME: Do Time dependent rainfall with poly
 
+
+
+
+    def test_inflow_using_circle(self):
+        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)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) )
+
+        domain.compute_forcing_terms()
+        #print domain.quantities['stage'].explicit_update
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
+        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+
+
+    def test_inflow_using_circle_function(self):
+        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)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) )
+
+        domain.compute_forcing_terms()
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
+        assert allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+        
 
 
@@ -3279 +3354 @@
         from Numeric import array
 
-        #Create basic mesh
+        # Create basic mesh
         points, vertices, boundary = rectangular(6, 6)
 
-        #Create shallow water domain
+        # Create shallow water domain
         domain = Domain(points, vertices, boundary)
         domain.smooth = False
-        domain.default_order=2
-
-        #IC
+        domain.default_order = 2
+
+        # IC
         def x_slope(x, y):
             return x/3
@@ -3293 +3368 @@
         domain.set_quantity('elevation', x_slope)
         domain.set_quantity('friction', 0)
-        domain.set_quantity('stage', 0.4) #Steady
+        domain.set_quantity('stage', 0.4) # Steady
 
         # Boundary conditions (reflective everywhere)
@@ -3308 +3383 @@
 
 
-        #Evolution
+        # Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
             stage =  domain.quantities['stage'].get_integral()
@@ -3314 +3389 @@
             ymom = domain.quantities['ymomentum'].get_integral()
 
-            if allclose(t, 6):  #Steady state reached
+            if allclose(t, 6):  # Steady state reached
                 steady_xmom = domain.quantities['xmomentum'].get_integral()
                 steady_ymom = domain.quantities['ymomentum'].get_integral()
@@ -3321 +3396 @@
             if t > 6:
                 #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
-                assert allclose(xmom, steady_xmom)
+                msg = 'xmom=%.2f, steady_xmom=%.2f' %(xmom, steady_xmom)
+                assert allclose(xmom, steady_xmom), msg
                 assert allclose(ymom, steady_ymom)
                 assert allclose(stage, steady_stage)
@@ -3405 +3481 @@
         domain = Domain(points, vertices, boundary)
         domain.smooth = False
-        domain.default_order=2
+        domain.default_order = 2
         domain.beta_w      = 0.9
         domain.beta_w_dry  = 0.9
@@ -3421 +3497 @@
         domain.check_integrity()
 
-        #Evolution
+        # Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
             pass# domain.write_time()
 
-        #Data from earlier version of abstract_2d_finite_volumes
+        # Data from earlier version of abstract_2d_finite_volumes
         assert allclose(domain.min_timestep, 0.0396825396825)
         assert allclose(domain.max_timestep, 0.0396825396825)
@@ -3566 +3642 @@
         #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
         domain.H0 = 0 # Backwards compatibility (6/2/7)
+        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
         domain.set_maximum_allowed_speed(1.0)        
 
@@ -3575 +3652 @@
         domain.check_integrity()
 
-        #Evolution
+        # Evolution
         for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
             pass
 
-
-        assert allclose(domain.min_timestep, 0.0210448446782)
+        msg = 'min step was %f instead of %f' %(domain.min_timestep,
+                                                0.0210448446782) 
+
+        assert allclose(domain.min_timestep, 0.0210448446782), msg
         assert allclose(domain.max_timestep, 0.0210448446782)
 
@@ -3630 +3709 @@
         domain.beta_vh_dry = 0.9        
         domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
-        domain.H0 = 0 # Backwards compatibility (6/2/7)        
+        domain.H0 = 0 # Backwards compatibility (6/2/7)
+        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)        
 
         # Boundary conditions
@@ -3682 +3762 @@
         domain.beta_vh_dry = 0.9
         domain.H0 = 0 # Backwards compatibility (6/2/7)
+        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)        
         domain.set_maximum_allowed_speed(1.0)