Changeset 274

Timestamp:

Sep 6, 2004, 1:57:41 PM (21 years ago)

Author:

ole

Message:

Cleaned up distribute

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (5 diffs)
• mesh.py (modified) (3 diffs)
• quantity.py (modified) (1 diff)
• shallow_water.py (modified) (2 diffs)
• test_quantity.py (modified) (1 diff)
• util.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -219 +219 @@
 
 
+    def get_name(self):
+        return self.filename    
+    
 
 
@@ -316 +319 @@
             
     def evolve_to_end(self, finaltime = 1.0):
-        """Iterate evolve generator all the way to the end
+        """Iterate evolve all the way to the end
         """
 
@@ -393 +396 @@
     def compute_forcing_terms(self):          
         """If there are any forcing functions driving the system
-        they should be defined in Domain subclass
+        they should be defined in Domain subclass and appended to
+        the list self.forcing_terms
         """
 
@@ -421 +425 @@
 
             #Clean up
+            #Note that Q.explicit_update is reset by compute_fluxes            
             Q.semi_implicit_update[:] = 0.0
-            Q.explicit_update[:] = 0.0 #Unnecessary as fluxes will set it
+
             
 
@@ -449 +454 @@
 ##############################################
 #Initialise module
-
-#C-extensions
-#import compile
-#if compile.can_use_C_extension('domain_ext.c'):
-#    compute_fluxes = compute_fluxes_c
-    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
-    #update_conserved_quantities = update_conserved_quantities_c
-#else:
-#    from shallow_water import compute_fluxes
-    #from python_versions import distribute_to_vertices_and_edges
-    #from python_versions import update_conserved_quantities
-
 
 #Optimisation with psyco

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -210 +210 @@
         #Build dictionary mapping from segments (2-tuple of points)
         #to left hand side edge (facing neighbouring triangle)
-        #Also build list of vertices containing indices of triangles  
 
         N = self.number_of_elements        
         neighbourdict = {}
-        vertexlist = [None]*len(self.coordinates)
         for i in range(N):
 
@@ -226 +224 @@
             neighbourdict[c,a] = (i, 1) #(id, edge)
 
-            #Register the vertices as keys mapping to 
-            #(triangle, edge) tuples associated with them
-            #This is used for smoothing
-            for vertex_id, v in enumerate([a,b,c]):
-                if vertexlist[v] is None:
-                    vertexlist[v] = []
-
-                vertexlist[v].append( (i, vertex_id) )    
 
         #Step 2:
@@ -260 +250 @@
 
 
-        self.vertexlist = vertexlist
     
     def build_vertexlist(self):

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -200 +200 @@
 
 
+    
+    def smooth_vertex_values(self, value_array='field_values',
+                             precision = None):
+        """ Smooths field_values or conserved_quantities data.
+        TODO: be able to smooth individual fields
+        NOTE:  This function does not have a test.
+        """
+
+        from Numeric import concatenate, zeros, Float, Int, array, reshape
+
+        
+        A,V = self.get_vertex_values(xy=False,
+                                     value_array=value_array,
+                                     smooth = True,
+                                     precision = precision)
+
+        #Set some field values
+        for volume in self:
+            for i,v in enumerate(volume.vertices):
+                if value_array == 'field_values':
+                    volume.set_field_values('vertex', i, A[v,:])
+                elif value_array == 'conserved_quantities':
+                    volume.set_conserved_quantities('vertex', i, A[v,:])
+
+        if value_array == 'field_values':
+            self.precompute()
+        elif value_array == 'conserved_quantities':
+            Volume.interpolate_conserved_quantities()
+        
+    
+    #Method for outputting model results
+    def get_vertex_values(self,
+                          xy=True, 
+                          smooth = None,
+                          precision = None,
+                          reduction = None):
+        """Return vertex values like an OBJ format
+
+        The vertex values are returned as one sequence in the 1D float array A.
+        If requested the coordinates will be returned in 1D arrays X and Y.
+        
+        The connectivity is represented as an integer array, V, of dimension
+        M x 3, where M is the number of volumes. Each row has three indices
+        into the X, Y, A arrays defining the triangle.
+
+        if smooth is True, vertex values corresponding to one common
+        coordinate set will be smoothed according to the given
+        reduction operator. In this case vertex coordinates will be
+        de-duplicated.
+        
+        If no smoothings is required, vertex coordinates and values will
+        be aggregated as a concatenation of of values at
+        vertices 0, vertices 1 and vertices 2
+
+        
+        Calling convention
+        if xy is True:
+           X,Y,A,V = get_vertex_values
+        else:
+           A,V = get_vertex_values         
+           
+        """
+
+        from Numeric import concatenate, zeros, Float, Int, array, reshape
+
+
+        if smooth is None:
+            smooth = self.domain.smooth
+
+        if precision is None:
+            precision = Float
+
+        if reduction is None:
+            reduction = self.domain.reduction
+
+        if smooth == True:
+            M = len(self.domain.vertexlist) #Number of unique vertices
+            
+            A = zeros((M,N), precision)
+            V = zeros((len(self),3), Int)            
+
+            #Rebuild triangle structure
+            for k, v in enumerate(self.vertexdict.keys()):            
+                for volume_id, vertex_id in self.vertexdict[v]:                
+                    V[volume_id, vertex_id] = k
+
+
+            #Smoothing loop
+            for j, i in enumerate(indices):
+                #For each quantity
+                #
+                #j is an index into returned quantities (0 <= j < N)
+                #i the given index into internal quantities
+                #(e.g. 0: bed elevation or 1: friction)
+
+                for k, v in enumerate(self.vertexdict.keys()):
+                    #Go through all contributions to vertex v
+                    #k is an index into returned vertices (0 <= k < M)
+                    #v is an index into internal vertices
+                    
+                    L = len(self.vertexdict[v])
+
+                    ulist = []
+                    for volume_id, vertex_id in self.vertexdict[v]:
+                        #For all contributing triange, vertex pairs
+                        
+                        cmd = 'u = Volume.%s_vertex%d[%d,%d]'\
+                              %(value_array, vertex_id, volume_id, i)
+                        exec(cmd)
+                        ulist.append(u)
+
+                    A[k,j] = reduction(ulist)    
+
+
+            if xy is True:
+                X = zeros(M, precision)
+                Y = zeros(M, precision)
+
+                for k, v in enumerate(self.vertexdict.keys()):            
+                    X[k], Y[k] = Point.coordinates[v].astype(precision)
+
+                return X, Y, A, V
+            else:
+                return A, V
+        else:
+            #Don't smooth
+
+            m = len(self)  #Number of volumes
+            M = 3*m        #Total number of unique vertices
+
+            A = self.vertex_values[:]
+
+            #Create connectivity
+            tmp = reshape(array(range(m)).astype(Int), (m,1))
+            V = concatenate( (tmp, tmp+m, tmp+2*m), axis = 1 )
+
+            #Do vertex coordinates    
+            if xy is True:                
+                V = self.domain.get_vertex_coordinates()
+                X = concatenate((V[:,0], V[:,2], V[:,4]), axis=0).\
+                    astype(precision)
+                Y = concatenate((V[:,1], V[:,3], V[:,5]), axis=0).\
+                    astype(precision)
+                
+                return X, Y, A , V           
+            else:
+                return A, V
+
+            
+            
+
+
 class Conserved_quantity(Quantity):
     """Class conserved quantity adds to Quantity:

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -386 +386 @@
 
 def distribute_to_vertices_and_edges(domain):
-
-    protect_against_infintesimal_and_negative_heights(domain)    
-    if domain.order == 1:
-        extrapolate_first_order(domain)
-    elif domain.order == 2:
-        extrapolate_second_order(domain)
-    else:
-        raise 'Unknown order'
-    
-    #Compute edge values
+    """Distribution from centroids to vertices 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 by taking a weighted average between shallow 
+    and deep cases.
+    
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+    
+    FIXME: more explanation about removal of artificial variability etc
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+    
+    """
+
+    #Remove very thin layers of water
+    protect_against_infintesimal_and_negative_heights(domain)
+
+    #Extrapolate all conserved quantities
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+            Q.limit()                                
+        else:
+            raise 'Unknown order'
+
+    #Take bed slevation into account when heights are small    
+    balance_deep_and_shallow(domain)
+    
+    #Compute edge values by interpolation
     for name in domain.conserved_quantities:
         Q = domain.quantities[name]
         Q.interpolate_from_vertices_to_edges()
                     
-    #domain.w.interpolate_from_vertices_to_edges()            
-    #domain.uh.interpolate_from_vertices_to_edges()            
-    #domain.vh.interpolate_from_vertices_to_edges()                    
     
 
@@ -443 +470 @@
     protect(domain.minimum_allowed_height, wc, zc, xmomc, ymomc)
     
-
-
-def extrapolate_first_order(domain):
-    """First order extrapolator 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
-    """
-
-    #Shortcuts
-    wc = domain.quantities['level'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values            
-    hc = wc - zc   #Water depths at centroids
-
-
-    #Update conserved quantities using straight first order
-    for name in domain.conserved_quantities:
-        Q = domain.quantities[name]
-        Q.extrapolate_first_order()
-
-
-
-
-    balance_deep_and_shallow(domain)
-
-
-    
-            
-        
-def extrapolate_second_order(domain):
-    """Second 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)
-
-    A weighted average between shallow 
-    and deep cases is as in the first order case.
-    
-    In addition, all conserved quantities get distributed as per a
-    piecewise linear function.
-    FIXME: more explanation about removal of artificial variability etc
-
-    Precondition:
-      All quantities defined at centroids and bed elevation defined at
-      vertices.
-
-    Postcondition
-      Conserved quantities defined at vertices
-    
-    """
-    
-
-    #Update conserved quantities using straight second order with limiter
-    for name in domain.conserved_quantities:
-        Q = domain.quantities[name]
-        Q.extrapolate_second_order()
-        Q.limit()        
-
-    balance_deep_and_shallow(domain)
-
 
 def balance_deep_and_shallow(domain):

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -254 +254 @@
 
 
-#     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_extrapolator(self):
         quantity = Conserved_quantity(self.mesh4)                

inundation/ga/storm_surge/pyvolution/util.py

@@ -42 +42 @@
 
 
+def mean(x):
+    from Numeric import sum
+    return sum(x)/len(x)
 
 ####################################################################