Changeset 296 for inundation/ga/storm_surge/pyvolution-1d

Timestamp:

Sep 14, 2004, 5:41:22 PM (21 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge/pyvolution-1d

Files:

• quantity.py (modified) (2 diffs)
• test_quantity.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution-1d/quantity.py

@@ -180 +180 @@
 
 
-if __name__ == "__main__":
-    from domain import Domain
-
-    points1 = [0.0, 1.0, 2.0, 3.0]
-    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
-    
-    D1 = Domain(points1)
-
-    Q1 = Quantity(D1, vertex_values)
-
-    print Q1.vertex_values
-    print Q1.centroid_values
-
-    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
-
-    Q1.set_values(new_vertex_values)
-    
-   
-    print Q1.vertex_values
-    print Q1.centroid_values
-
-    new_centroid_values = [20,30,40]
-    Q1.set_values(new_centroid_values,'centroids')
-            
-    print Q1.vertex_values
-    print Q1.centroid_values
-
-    def fun(x):
-
-        return x**2
-
-    Q1.set_values(fun,'vertices')    
-            
-    print Q1.vertex_values
-    print Q1.centroid_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):
-##         """Update centroid values based on values stored in
-##         explicit_update and semi_implicit_update as well as given timestep
-##         """
-
-##         from Numeric import sum, equal, ones, Float
-        
-##         N = self.centroid_values.shape[0]
-        
-##         #Explicit updates
-##         self.centroid_values += timestep*self.explicit_update
-            
-##         #Semi implicit updates
-##         denominator = ones(N, Float)-timestep*self.semi_implicit_update
-
-##         if sum(equal(denominator, 0.0)) > 0.0:
-##             msg = 'Zero division in semi implicit update. Call Stephen :-)'
-##             raise msg
-##         else:
-##             #Update conserved_quantities from semi implicit updates
-##             self.centroid_values /= denominator
-
-
-##     def compute_gradients(self):
-##         """Compute gradients of triangle surfaces defined by centroids of
-##         neighbouring volumes.
-##         If one face is on the boundary, use own centroid as neighbour centroid.
-##         If two or more are on the boundary, fall back to first order scheme.
-        
-##         Also return minimum and maximum of conserved quantities 
-##         """
-
-        
-##         from Numeric import array, zeros, Float
-##         from util import gradient
-
-##         N = self.centroid_values.shape[0]
-
-##         a = zeros(N, Float)
-##         b = zeros(N, Float)
-        
-##         for k in range(N):
-        
-##             number_of_boundaries = self.domain.number_of_boundaries[k]
-
-##             if number_of_boundaries == 3:                 
-##                 #We have zero neighbouring volumes -
-##                 #Fall back to first order scheme
-##                 pass
-            
-##             elif number_of_boundaries == 2:
-##                 #Special case where we have only one neighbouring volume.
-
-##                 #Find index of the one neighbour
-##                 for k0 in self.domain.neighbours[k,:]:
-##                     if k0 >= 0:
-##                         break
-
-##                 assert k0 != k
-##                 assert k0 >= 0
-
-##                 k1 = k  #Self
-
-##                 #Get data
-##                 q0 = self.centroid_values[k0]
-##                 q1 = self.centroid_values[k1]
-
-##                 x0, y0 = self.domain.centroids[k0] #V0 centroid
-##                 x1, y1 = self.domain.centroids[k1] #V1 centroid        
-
-##                 #Gradient
-##                 det = x0*y1 - x1*y0
-##                 if det != 0.0:
-##                     a[k] = (y1*q0 - y0*q1)/det
-##                     b[k] = (x0*q1 - x1*q0)/det
-
-##             else:
-##                 #One or zero missing neighbours
-##                 #In case of one boundary - own centroid
-##                 #has been inserted as a surrogate for the one
-##                 #missing neighbour in neighbour_surrogates
-
-##                 #Get data        
-##                 k0 = self.domain.surrogate_neighbours[k,0]
-##                 k1 = self.domain.surrogate_neighbours[k,1]
-##                 k2 = self.domain.surrogate_neighbours[k,2]
-
-##                 q0 = self.centroid_values[k0]
-##                 q1 = self.centroid_values[k1]
-##                 q2 = self.centroid_values[k2]                    
-
-##                 x0, y0 = self.domain.centroids[k0] #V0 centroid
-##                 x1, y1 = self.domain.centroids[k1] #V1 centroid
-##                 x2, y2 = self.domain.centroids[k2] #V2 centroid
-
-##                 #Gradient
-##                 a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        
-##         return a, b
+
+class Conserved_quantity(Quantity):
+    """Class conserved quantity adds to Quantity:
+
+    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 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):
+        """Update centroid values based on values stored in
+        explicit_update and semi_implicit_update as well as given timestep
+        """
+
+        from Numeric import sum, equal, ones, Float
+        
+        N = self.centroid_values.shape[0]
+        
+        #Explicit updates
+        self.centroid_values += timestep*self.explicit_update
+            
+        #Semi implicit updates
+        denominator = ones(N, Float)-timestep*self.semi_implicit_update
+
+        if sum(equal(denominator, 0.0)) > 0.0:
+            msg = 'Zero division in semi implicit update. Call Stephen :-)'
+            raise msg
+        else:
+            #Update conserved_quantities from semi implicit updates
+            self.centroid_values /= denominator
+
+
+    def compute_gradients(self):
+        """Compute gradients of triangle surfaces defined by centroids of
+        neighbouring volumes.
+        If one face is on the boundary, use own centroid as neighbour centroid.
+        If two or more are on the boundary, fall back to first order scheme.
+        
+        Also return minimum and maximum of conserved quantities 
+        """
+
+        
+        from Numeric import array, zeros, Float
+        from util import gradient
+
+        N = self.centroid_values.shape[0]
+
+        a = zeros(N, Float)
+        
+        for k in range(N):
+
+            # first and last elements have boundaries
+
+            if k == 1:
+
+                #Get data
+                k0 = k
+                k1 = k+1
+                
+                q0 = self.centroid_values[k0]
+                q1 = self.centroid_values[k1]
+
+                x0 = self.domain.centroids[k0] #V0 centroid
+                x1 = self.domain.centroids[k1] #V1 centroid        
+
+                #Gradient
+                a[k] = (q1 - q0)/(x1 - x0)
+
+            elif k == N-1:
+                
+                #Get data
+                k0 = k
+                k1 = k-1
+                
+                q0 = self.centroid_values[k0]
+                q1 = self.centroid_values[k1]
+
+                x0 = self.domain.centroids[k0] #V0 centroid
+                x1 = self.domain.centroids[k1] #V1 centroid        
+
+                #Gradient
+                a[k] = (q1 - q0)/(x1 - x0)
+                
+            else
+                #Interior Volume (2 neighbours)
+
+                #Get data
+                k0 = self.domain.neighbours[k,0]
+                k2 = self.domain.neighbours[k,1]
+
+                q0 = self.centroid_values[k0]
+                q1 = self.centroid_values[k]
+                q2 = self.centroid_values[k2]                    
+
+                x0 = self.domain.centroids[k0] #V0 centroid
+                x1 = self.domain.centroids[k1] #V1 centroid (Self)
+                x2 = self.domain.centroids[k2] #V2 centroid
+
+                #Gradient
+                a[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+        
+        return a
         
 
@@ -452 +406 @@
 
 
+
+if __name__ == "__main__":
+    from domain import Domain
+
+    points1 = [0.0, 1.0, 2.0, 3.0]
+    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
+    
+    D1 = Domain(points1)
+
+    Q1 = Quantity(D1, vertex_values)
+
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
+
+    Q1.set_values(new_vertex_values)
+    
+   
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    new_centroid_values = [20,30,40]
+    Q1.set_values(new_centroid_values,'centroids')
+            
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+    def fun(x):
+
+        return x**2
+
+    Q1.set_values(fun,'vertices')    
+            
+    print Q1.vertex_values
+    print Q1.centroid_values
+
+
+
+

inundation/ga/storm_surge/pyvolution-1d/test_quantity.py

@@ -76 +76 @@
                             [[1,2], [5,5], [0,9], [-6, 3], [3,4]])
         assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid
-
-
-##     def test_boundary_allocation(self):
-##         quantity = Conserved_quantity(self.mesh4,
-##                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
-
-##         assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
-
 
     def test_set_values(self):
@@ -275 +267 @@
 ## #         assert v2.conserved_quantities_vertex2 <= qmax
 ## #         assert v2.conserved_quantities_vertex2 >= qmin                    
-163.968025
+
 
 ## #         #Check that volume has been preserved    
@@ -390 +382 @@
 
 
-##     def test_update_explicit(self):
-##         quantity = Conserved_quantity(self.mesh4)
-
-##         #Test centroids
-##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-##         #Set explicit_update
-##         quantity.explicit_update = array( [1.,1.,1.,1.] )
-
-##         #Update with given timestep
-##         quantity.update(0.1)
-
-##         x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
-##         assert allclose( quantity.centroid_values, x)
-
-##     def test_update_semi_implicit(self):
-##         quantity = Conserved_quantity(self.mesh4)
-
-##         #Test centroids
-##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-##         #Set semi implicit update
-##         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
-
-##         #Update with given timestep
-##         quantity.update(0.1)
-
-##         x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
-##         assert allclose( quantity.centroid_values, x)
-
-##     def test_both_updates(self):
-##         quantity = Conserved_quantity(self.mesh4)
-
-##         #Test centroids
-##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-##         #Set explicit_update
-##         quantity.explicit_update = array( [4.,3.,2.,1.] )
-        
-##         #Set semi implicit update
-##         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
-
-##         #Update with given timestep
-##         quantity.update(0.1)
-
-##         x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
-##         x /= array( [.9,.9,.9,.9] )
-##         assert allclose( quantity.centroid_values, x)        
+    def test_update_explicit(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_update_semi_implicit(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)
+
+    def test_both_updates(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [4.,3.,2.,1.,0.0] )
+        
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
+
+        #Update with given timestep
+        quantity.update(0.1)
+
+        x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] )
+        x /= array( [.9,.9,.9,.9,.9] )
+        assert allclose( quantity.centroid_values, x)