Changeset 242

Timestamp:

Aug 30, 2004, 4:55:48 PM (21 years ago)

Author:

ole

Message:

Subclassed class Quantity into Conserved_quantity

Location:

inundation/ga/storm_surge/pyvolution

Files:

• config.py (modified) (1 diff)
• domain.py (modified) (2 diffs)
• quantity.py (modified) (8 diffs)
• shallow_water.py (modified) (2 diffs)
• test_quantity.py (modified) (10 diffs)

inundation/ga/storm_surge/pyvolution/config.py

@@ -62 +62 @@
 
 use_psyco = True  #Use psyco optimisations
-#use_psyco = False  #Do not use psyco optimisations
+use_psyco = False  #Do not use psyco optimisations
 
 

inundation/ga/storm_surge/pyvolution/domain.py

@@ -19 +19 @@
 
         from Numeric import zeros, Float 
-        from quantity import Quantity
+        from quantity import Quantity, Conserved_quantity
 
         #List of quantity names entering
@@ -34 +34 @@
         #Build dictionary of Quantity instances keyed by quantity names
         self.quantities = {}
-        for name in self.conserved_quantities + other_quantities:
+        for name in self.conserved_quantities:
+            self.quantities[name] = Conserved_quantity(self)
+        for name in other_quantities:
             self.quantities[name] = Quantity(self)
 

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -15 +15 @@
 """
 
-#FIXME: Make Conserved_quantity a subclass of Quantity 
-
-from mesh import Mesh
 
 class Quantity:
@@ -23 +20 @@
     def __init__(self, domain, vertex_values=None):
 
+        from mesh import Mesh
         from Numeric import array, zeros, Float
 
@@ -53 +51 @@
         self.edge_values = zeros((N, 3), Float)        
 
-        #Allocate space for boundary values
-        L = len(domain.boundary)
-        self.boundary_values = zeros(L, Float)
-
-        #Allocate space for updates of conserved quantities by
-        #flux calculations and forcing functions
-        self.explicit_update = zeros(N, Float )
-        self.semi_implicit_update = zeros(N, Float )
-                
         #Intialise centroid and edge_values
         self.interpolate()
@@ -71 +60 @@
         """
 
-        #FIXME: Write using vector operations (or in C)
         N = self.vertex_values.shape[0]        
         for i in range(N):
@@ -79 +67 @@
             
             self.centroid_values[i] = (v0 + v1 + v2)/3
-            
-            #FIXME: This is a duplicate from 'interpolate_from_vertices_to_edges'
-            self.edge_values[i, 0] = 0.5*(v1 + v2)
-            self.edge_values[i, 1] = 0.5*(v0 + v2)            
-            self.edge_values[i, 2] = 0.5*(v0 + v1)
-
+
+        self.interpolate_from_vertices_to_edges()
+
+
+    def interpolate_from_vertices_to_edges(self):
+        for k in range(self.vertex_values.shape[0]):
+            q0 = self.vertex_values[k, 0]
+            q1 = self.vertex_values[k, 1]
+            q2 = self.vertex_values[k, 2]
+            
+            self.edge_values[k, 0] = 0.5*(q1+q2)
+            self.edge_values[k, 1] = 0.5*(q0+q2)  
+            self.edge_values[k, 2] = 0.5*(q0+q1)
+
+            
+            
+    def set_values(self, X, location='vertices'):
+        """Set values for quantity
+
+        X: Compatible list, Numeric array (see below), constant or function
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If values are described a function, it will be evaluated at specified points
+
+        If selected location is vertices, values for centroid and edges
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+        """
+
+        if location not in ['vertices', 'centroids', 'edges']:
+            msg = 'Invalid location: %s' %location
+            raise msg
+
+        if X is None:
+            msg = 'Given values are None'
+            raise msg            
+        
+        import types
+        
+        if callable(X):
+            #Use function specific method
+            self.set_function_values(X, location)            
+        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
+            if location == 'centroids':
+                self.centroid_values[:] = X
+            elif location == 'edges':
+                self.edge_values[:] = X
+            else:
+                self.vertex_values[:] = X                
+
+        else:
+            #Use array specific method
+            self.set_array_values(X, location)
+
+        if location == 'vertices':
+            #Intialise centroid and edge_values
+            self.interpolate()
+            
+
+
+    def set_function_values(self, f, location='vertices'):
+        """Set values for quantity using specified function
+
+        f: x, y -> z Function where x, y and z are arrays
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+                  Default is "vertices"
+        """
+
+        if location == 'centroids':
+            P = self.domain.centroids
+            self.set_values(f(P[:,0], P[:,1]), location)
+        elif location == 'edges':
+            raise 'Not implemented: %s' %location
+        else:
+            #Vertices
+            P = self.domain.get_vertex_coordinates()
+            for i in range(3):
+                self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
+
+            
+    def set_array_values(self, values, location='vertices'):
+        """Set values for quantity
+
+        values: Numeric array
+        location: Where values are to be stored.
+                  Permissible options are: vertices, edges, centroid
+                  Default is "vertices"
+
+        In case of location == 'centroid' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. Otherwise it must be of dimension Nx3
+
+        The values will be stored in elements following their
+        internal ordering.
+
+        If selected location is vertices, values for centroid and edges
+        will be assigned interpolated values.
+        In any other case, only values for the specified locations
+        will be assigned and the others will be left undefined.
+        """
+
+        from Numeric import array, Float
+
+        values = array(values).astype(Float)
+
+        N = self.centroid_values.shape[0]
+        
+        msg = 'Number of values must match number of elements'
+        assert values.shape[0] == N, msg
+        
+        if location == 'centroids':
+            assert len(values.shape) == 1, 'Values array must be 1d'
+            self.centroid_values = values                        
+        elif location == 'edges':
+            assert len(values.shape) == 2, 'Values array must be 2d'
+            msg = 'Array must be N x 3'            
+            assert values.shape[1] == 3, msg
+            
+            self.edge_values = values
+        else:
+            assert len(values.shape) == 2, 'Values array must be 2d'
+            msg = 'Array must be N x 3'            
+            assert values.shape[1] == 3, msg
+            
+            self.vertex_values = values
+
+
+
+class Conserved_quantity(Quantity):
+    """Class conserved quantity adds to Quantity:
+
+    boundary values, storage and method for updating, and
+    methods for extrapolation from centropid to vertices inluding
+    gradients and limiters
+    """
+
+    def __init__(self, domain, vertex_values=None):
+        Quantity.__init__(self, domain, vertex_values)
+        
+        from Numeric import zeros, Float
+
+        #Allocate space for boundary values
+        L = len(domain.boundary)
+        self.boundary_values = zeros(L, Float)
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+
+        N = domain.number_of_elements
+        self.explicit_update = zeros(N, Float )
+        self.semi_implicit_update = zeros(N, Float )
+                
 
     def update(self, timestep):
@@ -134 +279 @@
                 #We have zero neighbouring volumes -
                 #Fall back to first order scheme
-
                 pass
             
@@ -162 +306 @@
                     a[k] = (y1*q0 - y0*q1)/det
                     b[k] = (x0*q1 - x1*q0)/det
-
 
             else:
@@ -288 +431 @@
 
             
-
-    def interpolate_from_vertices_to_edges(self):
-        #FIXME: Write using vector operations (or in C)
-        for k in range(self.vertex_values.shape[0]):
-            q0 = self.vertex_values[k, 0]
-            q1 = self.vertex_values[k, 1]
-            q2 = self.vertex_values[k, 2]
-            
-            self.edge_values[k, 0] = 0.5*(q1+q2)
-            self.edge_values[k, 1] = 0.5*(q0+q2)  
-            self.edge_values[k, 2] = 0.5*(q0+q1)
-
-            
-            
-    def set_values(self, X, location='vertices'):
-        """Set values for quantity
-
-        X: Compatible list, Numeric array (see below), constant or function
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
-                  Default is "vertices"
-
-        In case of location == 'centroid' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements in the mesh. Otherwise it must be of dimension Nx3
-
-        The values will be stored in elements following their
-        internal ordering.
-
-        If values are described a function, it will be evaluated at specified points
-
-        If selected location is vertices, values for centroid and edges
-        will be assigned interpolated values.
-        In any other case, only values for the specified locations
-        will be assigned and the others will be left undefined.
-        """
-
-
-        #FIXME: Should do location error check here only
-        
-        if location not in ['vertices', 'centroids', 'edges']:
-            msg = 'Invalid location: %s' %location
-            raise msg
-
-        if X is None:
-            msg = 'Given values are None'
-            raise msg            
-        
-        import types
-        
-        if callable(X):
-            #Use function specific method
-            self.set_function_values(X, location)            
-        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
-            if location == 'centroids':
-                self.centroid_values[:] = X
-            elif location == 'edges':
-                self.edge_values[:] = X
-            else:
-                self.vertex_values[:] = X                
-
-        else:
-            #Use array specific method
-            self.set_array_values(X, location)
-
-        if location == 'vertices':
-            #Intialise centroid and edge_values
-            self.interpolate()
-            
-
-
-    def set_function_values(self, f, location='vertices'):
-        """Set values for quantity using specified function
-
-        f: x, y -> z Function where x, y and z are arrays
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
-                  Default is "vertices"
-        """
-
-        if location == 'centroids':
-            P = self.domain.centroids
-            self.set_values(f(P[:,0], P[:,1]), location)
-        elif location == 'edges':
-            raise 'Not yet implemented: %s' %location
-        else:
-            #Vertices
-            P = self.domain.get_vertex_coordinates()
-            for i in range(3):
-                self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
-
-            
-    def set_array_values(self, values, location='vertices'):
-        """Set values for quantity
-
-        values: Numeric array
-        location: Where values are to be stored.
-                  Permissible options are: vertices, edges, centroid
-                  Default is "vertices"
-
-        In case of location == 'centroid' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements in the mesh. Otherwise it must be of dimension Nx3
-
-        The values will be stored in elements following their
-        internal ordering.
-
-        If selected location is vertices, values for centroid and edges
-        will be assigned interpolated values.
-        In any other case, only values for the specified locations
-        will be assigned and the others will be left undefined.
-        """
-
-        from Numeric import array, Float
-
-        values = array(values).astype(Float)
-
-        N = self.centroid_values.shape[0]
-        
-        msg = 'Number of values must match number of elements'
-        assert values.shape[0] == N, msg
-        
-        if location == 'centroids':
-            assert len(values.shape) == 1, 'Values array must be 1d'
-            self.centroid_values = values                        
-        elif location == 'edges':
-            assert len(values.shape) == 2, 'Values array must be 2d'
-            msg = 'Array must be N x 3'            
-            assert values.shape[1] == 3, msg
-            
-            self.edge_values = values
-        else:
-            assert len(values.shape) == 2, 'Values array must be 2d'
-            msg = 'Array must be N x 3'            
-            assert values.shape[1] == 3, msg
-            
-            self.vertex_values = values
-

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -352 +352 @@
         Q.interpolate_from_vertices_to_edges()            
 
+
 def protect_against_infinitesimal_heights(domain):
     """Protect against infinitesimal heights and high speeds    
@@ -397 +398 @@
     #Update
     for k in range(domain.number_of_elements):
-
-
         if domain.newstyle:
             if hc[k] < domain.minimum_allowed_height:        

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -87 +87 @@
 
     def test_boundary_allocation(self):
-        quantity = Quantity(self.mesh4,
+        quantity = Conserved_quantity(self.mesh4,
                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
 
@@ -189 +189 @@
 
     def test_gradient(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)
 
         #Set up for a gradient of (3,0) at mid triangle 
@@ -218 +218 @@
 
     def test_second_order_extrapolation2(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)        
 
         #Set up for a gradient of (3,1), f(x) = 3x+y
@@ -327 +327 @@
 
     def test_first_order_extrapolator(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)                
 
         #Test centroids
@@ -342 +342 @@
 
     def test_second_order_extrapolator(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)                        
 
         #Set up for a gradient of (3,0) at mid triangle 
@@ -375 +375 @@
         
 
-    def test_limiter(self):        
-        quantity = Quantity(self.mesh4)
+    def test_limiter(self):
+        quantity = Conserved_quantity(self.mesh4)
 
         #Create a deliberate overshoot (e.g. from gradient computation)
@@ -406 +406 @@
 
     def test_distribute_first_order(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)        
 
         #Test centroids
@@ -427 +427 @@
 
     def test_update_explicit(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids
@@ -443 +443 @@
 
     def test_update_semi_implicit(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids
@@ -459 +459 @@
 
     def test_both_updates(self):
-        quantity = Quantity(self.mesh4)
+        quantity = Conserved_quantity(self.mesh4)
 
         #Test centroids