Changeset 198 for inundation/ga

Timestamp:

Aug 20, 2004, 3:28:34 PM (20 years ago)

Author:

ole

Message:

More work on limiters
Started working on gradients

Location:

inundation/ga/storm_surge/pyvolution

Files:

• config.py (modified) (1 diff)
• mesh.py (modified) (2 diffs)
• quantity.py (modified) (2 diffs)
• shallow_water.py (modified) (6 diffs)
• test_quantity.py (modified) (2 diffs)
• test_shallow_water.py (modified) (4 diffs)
• util.py (modified) (2 diffs)
• util_ext.c (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/config.py

@@ -57 +57 @@
 
 use_extensions = True   #Try to use C-extensions
-#use_extensions = False   #Do not use C-extensions
+use_extensions = False   #Do not use C-extensions
 
 use_psyco = True  #Use psyco optimisations

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -201 +201 @@
     def build_neighbour_structure(self):
         """Update all registered triangles to point to their neighbours.
+        
+        Also, keep a tally of the number of boundaries for each triangle
 
         Postconditions:
@@ -297 +299 @@
         if they exist. Otherwise point to the triangle itself.
 
-        Also, keep a tally of the number of boundaries for each triangle
-
         The surrogate neighbour structure is useful for computing gradients
         based on centroid values of neighbours.

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -106 +106 @@
             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.
+
+                k0 = k  #Self
+                
+                #Find index of the one neighbour
+                for k1 in self.domain.neighbours[k,:]:
+                    if k1 >= 0:
+                        break
+                assert k1 != k0
+                assert k1 >= 0                     
+                    
+                #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)
+                                      #array([q0]), array([q1]), array([q2]))
+        
+        return a, b
+        
 
     def second_order_limiter(self):
@@ -189 +267 @@
             self.first_order_limiter()
         elif order == 2:
+            a, b = self.compute_gradients()
+            
             self.second_order_limiter()
         else:

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -267 +267 @@
 ####################################
 # Functions for gradient limiting (distribute_to_vertices_and_edges)
-def first_order_limiter(domain):
-    """First order limiter function, specific to the shallow water wave
-    equation.
-    
-    It will ensure that h (w-z) is always non-negative even in the
-    presence of steep bed-slopes (see comment in code)
-    In addition, momemtums get distributed as constant values.
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
+
+def protect_against_infinitesimal_heights_centroid(domain):
+    """Adjust height and momentum at centroid if height is less than
+    minimum allowed height
     """
-
+    
     #Water levels at centroids
     wc = domain.quantities['level'].centroid_values
@@ -296 +286 @@
     ymomc = domain.quantities['ymomentum'].centroid_values        
 
-    #Water levels at vertices
-    wv = domain.quantities['level'].vertex_values
-
-    #Bed elevations at vertices
-    zv = domain.quantities['elevation'].vertex_values
-    
     for k in range(domain.number_of_elements):
         #Protect against infinitesimal heights and high velocities
@@ -312 +296 @@
             xmomc[k] = ymomc[k] = 0.0
                 
+
+    
+
+def first_order_limiter(domain):
+    """First order limiter function, specific to the shallow water wave
+    equation.
+    
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes (see comment in code)
+    In addition, momemtums get distributed as constant values.
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+    """
+
+    #FIXME: Might go elsewhere
+    protect_against_infinitesimal_heights_centroid(domain)
+    
+    #Update conserved quantities using straight first order
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        Q.first_order_limiter()
+
+        
+    #Water levels at centroids
+    wc = domain.quantities['level'].centroid_values
+
+    #Bed elevations at centroids
+    zc = domain.quantities['elevation'].centroid_values
+
+    #Water depths at centroids    
+    hc = wc - zc
+
+    #Water levels at vertices
+    wv = domain.quantities['level'].vertex_values
+
+    #Bed elevations at vertices
+    zv = domain.quantities['elevation'].vertex_values
                 
+    #Water depths at vertices
+    hv = wv-zv
+    
+    
+    #Computed weighted balance between constant levels and and
+    #levels parallel to the bed elevation.     
+    for k in range(domain.number_of_elements):                
                 
         #Compute maximal variation in bed elevation
@@ -327 +360 @@
         #  wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) = 
         #  (1-alpha)*zvi + alpha*zc + hc = 
-        #  zvi + hc + alpha*(zc-zvi);
+        #  zvi + hc + alpha*(zc-zvi)
         #
-        #where alpha in [0,1] and defined as the ration between hc and 
+        #where alpha in [0,1] and defined as the ratio between hc and 
         #the maximal difference from zc to zv0, zv1 and zv2
+        #
+        #Mathematically the following can be continued on using hc as
+        #  wvi = 
+        #  zvi + hc + alpha*(zc+hc-zvi-hc) =
+        #  zvi + hc + alpha*(hvi-hc)
+        #since hvi = zc+hc-zvi in the constant case
+
         
         if z_range > 0.0:
@@ -339 +379 @@
         #Update stage
         for i in range(3):
-            wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
-
-    #Update momentum using straight first order
-    for name in ['xmomentum', 'ymomentum']:
-        Q = domain.quantities[name]
-        for i in range(3):
-            Q.vertex_values[:, i] = Q.centroid_values
+            #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
+            wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+
 
         
@@ -370 +406 @@
     
     """
+
+    #FIXME: Might go elsewhere
+    protect_against_infinitesimal_heights_centroid(domain)
+    
+    #Update conserved quantities using straight second order
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        Q.second_order_limiter()
+
 
     #Water levels at centroids

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -145 +145 @@
 
 
+    def test_gradient(self):
+        quantity = Quantity(self.mesh4)
+
+        #Set up for a gradient of (1,0)
+        quantity.set_values([2.0/3, 4.0/3, 8.0/3, 2.0/3],
+                            location = 'centroids')
+        
+        #Gradients
+        quantity.compute_gradients()
+
+        #FIXME:
+
+        #HERTIL
+        
+        #Gradient of fitted pwl surface    
+        #a, b = compute_gradient(v2.id)
+
+        #assert a == 1.0
+        #assert b == 0.0
+
+
+        #And now for the second order stuff        
+        # - the full second order extrapolation
+        #domain.order = 2
+        #domain.limiter = dummy_limiter
+        #distribute_to_vertices_and_edges(domain)       
+    
+        #assert v2.conserved_quantities_vertex0 == 0.0
+        #assert v2.conserved_quantities_vertex1 == 2.0        
+        #assert v2.conserved_quantities_vertex2 == 2.0        
+
+
+
+#     def test_second_order_extrapolation2(self):
+
+#         initialise_consecutive_datastructure(points=6+4, elements=4)                
+        
+#         a = Point (0.0, 0.0)
+#         b = Point (0.0, 2.0)
+#         c = Point (2.0, 0.0)
+#         d = Point (0.0, 4.0)
+#         e = Point (2.0, 2.0)
+#         f = Point (4.0, 0.0)
+
+#         #Set up for a gradient of (3,1), f(x) = 3x+y
+#         v1 = Volume(b,a,c,array([2.0+2.0/3,0,0]))        
+#         v2 = Volume(b,c,e,array([4.0+4.0/3,0,0]))
+#         v3 = Volume(e,c,f,array([8.0+2.0/3,0,0]))
+#         v4 = Volume(d,b,e,array([2.0+8.0/3,0,0]))
+
+#         #Setup neighbour structure
+#         domain = Domain([v1,v2,v3,v4])
+#         domain.precompute()
+#       domain.check_integrity()
+        
+#         #Lets's check first order first, hey
+#       domain.order = 1
+#       domain.limiter = dummy_limiter
+#         distribute_to_vertices_and_edges(domain)              
+        
+#         assert allclose(v2.conserved_quantities_vertex0,
+#                         v2.conserved_quantities_centroid)
+#         assert allclose(v2.conserved_quantities_vertex1,
+#                         v2.conserved_quantities_centroid)
+#         assert allclose(v2.conserved_quantities_vertex2,
+#                         v2.conserved_quantities_centroid)                
+
+
+#         #Flux across right edge of volume 1
+#         #Outward pointing normal vector
+#         from shallow_water import flux_using_stage as flux_function
+#         normals = Volume.normals
+        
+#         normal = Vector.coordinates[normals[1][2]]
+#         ql = Volume.conserved_quantities_face2[1]
+#         qr = Volume.conserved_quantities_face1[0]
+#         fl = array([0.,0.])
+#         fr = array([0.,0.])
+#         flux0, max_speed = flux_function(normal, ql, qr, fl, fr)
+
+#         #print flux0, max_speed        
+
+#         #print
+#         #print v1.conserved_quantities_face0,\
+#         #      v2.conserved_quantities_face0,\
+#         #      v3.conserved_quantities_face0,\
+#         #      v4.conserved_quantities_face0              
+#         #print
+#         #edgelengths = Volume.geometric[:,2:]
+#         #print
+#         #print
+        
+#         from python_versions import compute_flux
+#         compute_flux(domain, 100)
+#         F1 = Volume.explicit_update
+
+#         from domain import compute_flux
+#         compute_flux(domain, 100)
+#         F2 = Volume.explicit_update        
+
+#         assert allclose(F1, F2)
+        
+#         #print F1
+#         #print F2
+        
+
+
+#         #Gradient of fitted pwl surface    
+#       a, b = compute_gradient(v2.id)  
+
+#         assert abs(a[0] - 3.0) < epsilon
+#         assert abs(b[0] - 1.0) < epsilon
+#         #assert qmin[0] == 2.0 + 2.0/3
+#         #assert qmax[0] == 8.0 + 2.0/3                
+
+#         #And now for the second order stuff        
+#         # - the full second order extrapolation
+#       domain.order = 2
+#       domain.limiter = dummy_limiter
+#         distribute_to_vertices_and_edges(domain)              
+
+#         assert allclose(v2.conserved_quantities_vertex0[0], 2.0)
+#         assert allclose(v2.conserved_quantities_vertex1[0], 6.0)        
+#         assert allclose(v2.conserved_quantities_vertex2[0], 8.0)
+
+
+
+
+#     def test_limiter(self):
+
+#         initialise_consecutive_datastructure(points=6+4, elements=4)          
+        
+#         a = Point (0.0, 0.0)
+#         b = Point (0.0, 2.0)
+#         c = Point (2.0, 0.0)
+#         d = Point (0.0, 4.0)
+#         e = Point (2.0, 2.0)
+#         f = Point (4.0, 0.0)
+
+#         #Set up for a gradient of (3,1), f(x) = 3x+y
+#         v1 = Volume(b,a,c,array([0.0,0,0]))        
+#         v2 = Volume(b,c,e,array([1.0,0,0]))
+#         v3 = Volume(e,c,f,array([10.0,0,0]))
+#         v4 = Volume(d,b,e,array([0.0,0,0]))
+
+#         #Setup neighbour structure
+#         domain = Domain([v1,v2,v3,v4])
+#         domain.precompute()
+        
+#         #Lets's check first order first, hey
+#       domain.order = 1
+#       domain.limiter = None
+#         distribute_to_vertices_and_edges(domain)
+#         assert allclose(v2.conserved_quantities_vertex0,
+#                         v2.conserved_quantities_centroid)
+#         assert allclose(v2.conserved_quantities_vertex1,
+#                         v2.conserved_quantities_centroid)
+#         assert allclose(v2.conserved_quantities_vertex2,
+#                         v2.conserved_quantities_centroid)
+
+
+#         #Gradient of fitted pwl surface
+#         a, b = compute_gradient(v2.id)
+        
+
+#         assert abs(a[0] - 5.0) < epsilon
+#         assert abs(b[0]) < epsilon        
+#         #assert qminr[0] == 0.0
+#         #assert qmaxr[0] == 10.0
+        
+#         #And now for the second order stuff        
+#         # - the full second order extrapolation
+#       domain.order = 2        
+#         distribute_to_vertices_and_edges(domain)
+
+
+#         qmin = qmax = v2.conserved_quantities_centroid
+
+#         qmin = minimum(qmin, v1.conserved_quantities_centroid)
+#         qmax = maximum(qmax, v1.conserved_quantities_centroid)
+
+#         qmin = minimum(qmin, v3.conserved_quantities_centroid)
+#         qmax = maximum(qmax, v3.conserved_quantities_centroid)
+
+#         qmin = minimum(qmin, v4.conserved_quantities_centroid)
+#         qmax = maximum(qmax, v4.conserved_quantities_centroid)                    
+#         #assert qminr == qmin
+#         #assert qmaxr == qmax
+
+#         assert v2.conserved_quantities_vertex0 <= qmax
+#         assert v2.conserved_quantities_vertex0 >= qmin
+#         assert v2.conserved_quantities_vertex1 <= qmax
+#         assert v2.conserved_quantities_vertex1 >= qmin
+#         assert v2.conserved_quantities_vertex2 <= qmax
+#         assert v2.conserved_quantities_vertex2 >= qmin                    
+
+
+#         #Check that volume has been preserved    
+
+#         q = v2.conserved_quantities_centroid[0]
+#         w = (v2.conserved_quantities_vertex0[0] +
+#              v2.conserved_quantities_vertex1[0] +
+#              v2.conserved_quantities_vertex2[0])/3
+
+#         assert allclose(q, w)
+
+
+        
+
+
     def test_first_order_limiter(self):
         quantity = Quantity(self.mesh4)
@@ -161 +371 @@
 
     def test_second_order_limiter(self):
-
-        #FIXME:!!!!!!!
-        
-        quantity = Quantity(self.mesh4)
-
-        quantity.set_values([2., 4., 5., 5.], location = 'centroids')
-        #Assume gradient
-
-        #Extrapolate
+        quantity = Quantity(self.mesh4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,3,3], [6,2,2], [3,8,5], [3,5,8]])
+
+        #Limit and extrapolate
         quantity.second_order_limiter()
-        #print quantity.vertex_values
-
-
+
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.vertex_values[1,0] <= quantity.vertex_values[2,2]
+        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
+        
+        assert quantity.vertex_values[1,1] <= quantity.vertex_values[3,0]
+        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
+        
+        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
+        assert quantity.vertex_values[1,2] >= quantity.vertex_values[0,1]
+
+        
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/3)
+        
         
         #Check vertices but not edge values

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -22 +22 @@
         ql = zeros( 3, Float )
         qr = zeros( 3, Float )
-        normal = zeros( 3, Float )
+        normal = zeros( 2, Float )
         zl = zr = 0.
         flux, max_speed = flux_function(normal, ql, qr, zl, zr)
@@ -110 +110 @@
 
         assert domain.get_conserved_quantities(0, edge=1) == 0.    
-
 
 
@@ -602 +601 @@
                                       [val3, val3, val3]])
 
+        E = domain.quantities['elevation'].vertex_values
         L = domain.quantities['level'].vertex_values        
-        E = domain.quantities['elevation'].vertex_values
 
         
@@ -612 +611 @@
         domain.first_order_limiter()
 
-        #Check that all levels are above elevation (within eps) 
+        #Check that all levels are above elevation (within eps)
+
+        
         assert alltrue(alltrue(greater_equal(L,E-epsilon)))                
         

inundation/ga/storm_surge/pyvolution/util.py

@@ -100 +100 @@
            'Vector of conserved quantities must have length 3'\
            'for 2D shallow water equation'
+
     try:
         l = len(normal)
@@ -105 +106 @@
         raise 'Normal vector must be an Numeric array'
 
+    #FIXME: Put this test into C-extension as well
     assert l == 2, 'Normal vector must have 2 components'
 

inundation/ga/storm_surge/pyvolution/util_ext.c

@@ -117 +117 @@
     return NULL;
   }  
-    
+  
   //Allocate space for return vectors a and b