Changeset 1042

Timestamp:

Mar 8, 2005, 5:35:19 PM (19 years ago)

Author:

steve

Message:

 

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (103 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -21 +21 @@
 The quantities are
 
-symbol    variable name    explanation 
+symbol    variable name    explanation
 x         x                horizontal distance from origin [m]
-y         y                vertical distance from origin [m] 
+y         y                vertical distance from origin [m]
 z         elevation        elevation of bed on which flow is modelled [m]
 h         height           water height above z [m]
@@ -45 +45 @@
 
 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-Geoscience Australia, 2004   
+Geoscience Australia, 2004
 """
 
@@ -66 +66 @@
 
         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
-        other_quantities = ['elevation', 'friction'] 
-        
+        other_quantities = ['elevation', 'friction']
+
         Generic_domain.__init__(self, coordinates, vertices, boundary,
                                 conserved_quantities, other_quantities,
@@ -84 +84 @@
         #Stored output
         self.store = False
-        self.format = 'sww'    
+        self.format = 'sww'
         self.smooth = True
 
-        #Reduction operation for get_vertex_values              
+        #Reduction operation for get_vertex_values
         #from util import mean
-        #self.reduction = mean 
+        #self.reduction = mean
         self.reduction = min  #Looks better near steep slopes
-        
+
         self.quantities_to_be_stored = ['stage']
 
-        
+
         #Establish shortcuts to relevant quantities (for efficiency)
         #self.w = self.quantities['stage']
-        #self.uh = self.quantities['xmomentum']                 
-        #self.vh = self.quantities['ymomentum']                         
-        #self.z = self.quantities['elevation']          
-        #self.eta = self.quantities['friction']                 
+        #self.uh = self.quantities['xmomentum']
+        #self.vh = self.quantities['ymomentum']
+        #self.z = self.quantities['elevation']
+        #self.eta = self.quantities['friction']
 
     def check_integrity(self):
@@ -113 +113 @@
         msg = 'Third conserved quantity must be "ymomentum"'
         assert self.conserved_quantities[2] == 'ymomentum', msg
-        
+
 
     def compute_fluxes(self):
@@ -122 +122 @@
     def distribute_to_vertices_and_edges(self):
         #Call correct module function
-        #(either from this module or C-extension)        
+        #(either from this module or C-extension)
         distribute_to_vertices_and_edges(self)
 
 
-    #FIXME: Under construction    
+    #FIXME: Under construction
 #     def set_defaults(self):
 #         """Set default values for uninitialised quantities.
@@ -139 +139 @@
 #             if not self.quantities.has_key(name):
 #                 if name == 'stage':
-                    
+
 #                     if self.quantities.has_key('elevation'):
 #                         z = self.quantities['elevation'].vertex_values
-#                         self.set_quantity(name, z) 
+#                         self.set_quantity(name, z)
 #                     else:
-#                         self.set_quantity(name, 0.0)                        
-#                 else:    
+#                         self.set_quantity(name, 0.0)
+#                 else:
 #                     self.set_quantity(name, 0.0)
-                    
+
 
 
@@ -161 +161 @@
 #         #            w[k, i] = z[k, i]
 
-                    
+
 #         #self.quantities['stage'].interpolate()
-                
+
 
 
@@ -170 +170 @@
         """
 
-        #Call check integrity here rather than from user scripts 
+        #Call check integrity here rather than from user scripts
         #self.check_integrity()
-        
-        
-        #Initial update of vertex and edge values before any storage
-        #and or visualisation
+
+
+        #Initial update of vertex and edge values before any storage
+        #and or visualisation
         self.distribute_to_vertices_and_edges()
-        
-                
+
+
         #Initialise real time viz if requested
         if self.visualise is True and self.time == 0.0:
@@ -184 +184 @@
             visualise.create_surface(self)
 
-        #Store model data, e.g. for visualisation    
+        #Store model data, e.g. for visualisation
         if self.store is True and self.time == 0.0:
             self.initialise_storage()
@@ -201 +201 @@
                 visualise.update(self)
 
-            #Store model data, e.g. for subsequent visualisation    
+            #Store model data, e.g. for subsequent visualisation
             if self.store is True:
                 self.store_timestep(self.quantities_to_be_stored)
-                
-                #FIXME: Could maybe be taken from specified list 
-                #of 'store every step' quantities       
-
-            #Pass control on to outer loop for more specific actions    
+
+                #FIXME: Could maybe be taken from specified list
+                #of 'store every step' quantities
+
+            #Pass control on to outer loop for more specific actions
             yield(t)
-        
-
-    def initialise_storage(self):       
+
+
+    def initialise_storage(self):
         """Create and initialise self.writer object for storing data.
         Also, save x,y and bed elevation
-        """
+        """
 
         import data_manager
-        
+
         #Initialise writer
         self.writer = data_manager.get_dataobject(self, mode = 'w')
@@ -233 +233 @@
         """
         self.writer.store_timestep(name)
-        
+
 
 #Rotation of momentum vector
@@ -241 +241 @@
 
     If direction is negative the rotation is inverted.
-    
+
     Input vector is preserved
 
@@ -248 +248 @@
 
     from Numeric import zeros, Float
-    
+
     assert len(q) == 3,\
            'Vector of conserved quantities must have length 3'\
@@ -259 +259 @@
 
     assert l == 2, 'Normal vector must have 2 components'
-    
- 
+
+
     n1 = normal[0]
-    n2 = normal[1]    
-    
+    n2 = normal[1]
+
     r = zeros(len(q), Float) #Rotated quantities
     r[0] = q[0]              #First quantity, height, is not rotated
@@ -273 +273 @@
     r[1] =  n1*q[1] + n2*q[2]
     r[2] = -n2*q[1] + n1*q[2]
-    
+
     return r
 
@@ -279 +279 @@
 
 ####################################
-# Flux computation        
-def flux_function(normal, ql, qr, zl, zr): 
+# Flux computation
+def flux_function(normal, ql, qr, zl, zr):
     """Compute fluxes between volumes for the shallow water wave equation
     cast in terms of w = h+z using the 'central scheme' as described in
@@ -313 +313 @@
         u_left = 0.0  #Could have been negative
         h_left = 0.0
-    else:    
+    else:
         u_left  = uh_left/h_left
 
@@ -325 +325 @@
         u_right = 0.0  #Could have been negative
         h_right = 0.0
-    else:    
+    else:
         u_right  = uh_right/h_right
 
     vh_left  = q_left[2]
-    vh_right = q_right[2]        
-
-    soundspeed_left  = sqrt(g*h_left) 
+    vh_right = q_right[2]
+
+    soundspeed_left  = sqrt(g*h_left)
     soundspeed_right = sqrt(g*h_right)
 
     #Maximal wave speed
     s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    
-    #Minimal wave speed        
+
+    #Minimal wave speed
     s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
 
@@ -346 +346 @@
     flux_right = array([u_right*h_right,
                         u_right*uh_right + 0.5*g*h_right**2,
-                        u_right*vh_right])    
+                        u_right*vh_right])
 
     denom = s_max-s_min
@@ -352 +352 @@
         edgeflux = array([0.0, 0.0, 0.0])
         max_speed = 0.0
-    else:    
+    else:
         edgeflux = (s_max*flux_left - s_min*flux_right)/denom
         edgeflux += s_max*s_min*(q_right-q_left)/denom
@@ -359 +359 @@
         max_speed = max(abs(s_max), abs(s_min))
 
-    return edgeflux, max_speed        
+    return edgeflux, max_speed
 
 
@@ -370 +370 @@
     Fluxes across each edge are scaled by edgelengths and summed up
     Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities 
-    stage, xmomentum and ymomentum    
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
 
     The maximal allowable speed computed by the flux_function for each volume
@@ -379 +379 @@
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes. 
+      domain.timestep is set to the largest step satisfying all volumes.
     """
 
@@ -386 +386 @@
 
     N = domain.number_of_elements
-    
+
     #Shortcuts
     Stage = domain.quantities['stage']
     Xmom = domain.quantities['xmomentum']
     Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']        
+    Bed = domain.quantities['elevation']
 
     #Arrays
@@ -397 +397 @@
     xmom =  Xmom.edge_values
     ymom =  Ymom.edge_values
-    bed =   Bed.edge_values    
+    bed =   Bed.edge_values
 
     stage_bdry = Stage.boundary_values
     xmom_bdry =  Xmom.boundary_values
     ymom_bdry =  Ymom.boundary_values
-    
+
     flux = zeros(3, Float) #Work array for summing up fluxes
 
     #Loop
-    timestep = float(sys.maxint)    
+    timestep = float(sys.maxint)
     for k in range(N):
 
@@ -416 +416 @@
 
             #Quantities at neighbour on nearest face
-            n = domain.neighbours[k,i] 
+            n = domain.neighbours[k,i]
             if n < 0:
                 m = -n-1 #Convert negative flag to index
                 qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
                 zr = zl #Extend bed elevation to boundary
-            else:    
+            else:
                 m = domain.neighbour_edges[k,i]
                 qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                zr = bed[n, m]                
-
-                
-            #Outward pointing normal vector   
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
             normal = domain.normals[k, 2*i:2*i+2]
 
@@ -447 +447 @@
         Ymom.explicit_update[k] = flux[2]
 
-    
-    domain.timestep = timestep    
+
+    domain.timestep = timestep
 
 
@@ -459 +459 @@
 
     N = domain.number_of_elements
-    
+
     #Shortcuts
     Stage = domain.quantities['stage']
     Xmom = domain.quantities['xmomentum']
     Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']        
-
-    timestep = float(sys.maxint)    
+    Bed = domain.quantities['elevation']
+
+    timestep = float(sys.maxint)
     from shallow_water_ext import compute_fluxes
     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
@@ -472 +472 @@
                                      domain.neighbour_edges,
                                      domain.normals,
-                                     domain.edgelengths,                       
+                                     domain.edgelengths,
                                      domain.radii,
                                      domain.areas,
-                                     Stage.edge_values, 
-                                     Xmom.edge_values, 
-                                     Ymom.edge_values,  
-                                     Bed.edge_values,   
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
                                      Stage.boundary_values,
                                      Xmom.boundary_values,
@@ -485 +485 @@
                                      Xmom.explicit_update,
                                      Ymom.explicit_update)
-                                     
+
 
 ####################################
@@ -496 +496 @@
 
     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 
+    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
 
@@ -510 +510 @@
     Postcondition
       Conserved quantities defined at vertices
-    
+
     """
 
@@ -523 +523 @@
         elif domain.order == 2:
             Q.extrapolate_second_order()
-            Q.limit()                                
+            Q.limit()
         else:
             raise 'Unknown order'
-        
+
     #Take bed elevation into account when water heights are small
     balance_deep_and_shallow(domain)
@@ -534 +534 @@
         Q = domain.quantities[name]
         Q.interpolate_from_vertices_to_edges()
-                    
+
 
 
@@ -543 +543 @@
 
     #FIXME: Experimental (from old version). Not in use at the moment
-    
+
     #Shortcuts
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
     xmomv = domain.quantities['xmomentum'].vertex_values
-    ymomv = domain.quantities['ymomentum'].vertex_values            
+    ymomv = domain.quantities['ymomentum'].vertex_values
     hv = wv - zv  #Water depths at vertices
 
@@ -554 +554 @@
     for k in range(domain.number_of_elements):
         hmax = max(hv[k, :])
-        
+
         if hmax < domain.minimum_allowed_height:
             #Control stage
@@ -561 +561 @@
             #Control momentum
             xmomv[k,:] = ymomv[k,:] = 0.0
-    
+
 
 def protect_against_infinitesimal_and_negative_heights(domain):
     """Protect against infinitesimal heights and associated high velocities
     """
-    
+
     #FIXME: Look here for error
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].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    
+    ymomc = domain.quantities['ymomentum'].centroid_values
+    hc = wc - zc  #Water depths at centroids
 
     #print zc
-    #print '1', wc      
+    #print '1', wc
     #Update
     for k in range(domain.number_of_elements):
 
-        if hc[k] < domain.minimum_allowed_height: 
+        if hc[k] < domain.minimum_allowed_height:
             #Control stage
             wc[k] = zc[k]
@@ -587 +587 @@
             #Control momentum
             xmomc[k] = ymomc[k] = 0.0
-        
-    #print '2', wc              
+
+    #print '2', wc
 
 
@@ -594 +594 @@
     """Protect against infinitesimal heights and associated high velocities
     """
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values            
+    ymomc = domain.quantities['ymomentum'].centroid_values
 
     from shallow_water_ext import protect
 
     protect(domain.minimum_allowed_height, wc, zc, xmomc, ymomc)
-    
-
-
-def h_limiter(domain):        
+
+
+
+def h_limiter(domain):
     """Limit slopes for each volume to eliminate artificial variance
     introduced by e.g. second order extrapolator
-    
+
     limit on h = w-z
 
@@ -621 +621 @@
     N = domain.number_of_elements
     beta_h = domain.beta_h
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     hc = wc - zc
-    
+
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
     hv = wv-zv
 
-    hvbar = zeros(hv.shape, Float) #h-limited values 
-    
+    hvbar = zeros(hv.shape, Float) #h-limited values
+
     #Find min and max of this and neighbour's centroid values
     hmax = zeros(hc.shape, Float)
-    hmin = zeros(hc.shape, Float)        
-    
+    hmin = zeros(hc.shape, Float)
+
     for k in range(N):
         hmax[k] = hmin[k] = hc[k]
@@ -643 +643 @@
             if n >= 0:
                 hn = hc[n] #Neighbour's centroid value
-                
+
                 hmin[k] = min(hmin[k], hn)
                 hmax[k] = max(hmax[k], hn)
-      
-    
+
+
     #Diffences between centroids and maxima/minima
-    dhmax = hmax - hc 
+    dhmax = hmax - hc
     dhmin = hmin - hc
-        
+
     #Deltas between vertex and centroid values
     dh = zeros(hv.shape, Float)
@@ -657 +657 @@
         dh[:,i] = hv[:,i] - hc
 
-    #Phi limiter    
+    #Phi limiter
     for k in range(N):
-        
+
         #Find the gradient limiter (phi) across vertices
         phi = 1.0
@@ -666 +666 @@
             if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
             if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
-            
-            phi = min( min(r*beta_h, 1), phi )    
+
+            phi = min( min(r*beta_h, 1), phi )
 
         #Then update using phi limiter
-        for i in range(3):            
+        for i in range(3):
             hvbar[k,i] = hc[k] + phi*dh[k,i]
 
@@ -677 +677 @@
 
 
-def h_limiter_c(domain):        
+def h_limiter_c(domain):
     """Limit slopes for each volume to eliminate artificial variance
     introduced by e.g. second order extrapolator
-    
+
     limit on h = w-z
 
@@ -693 +693 @@
     N = domain.number_of_elements
     beta_h = domain.beta_h
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     hc = wc - zc
-    
+
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
@@ -709 +709 @@
     return hvbar
 
-                                            
+
 def balance_deep_and_shallow(domain):
     """Compute linear combination between stage as computed by
-    gradient-limiters limiting using w, and stage computed as 
+    gradient-limiters limiting using w, and stage computed as
     constant height above bed and limited using h.
     The former takes precedence when heights are large compared to the
@@ -721 +721 @@
     modes.
 
-    The h-limiter is always applied irrespective of the order. 
+    The h-limiter is always applied irrespective of the order.
     """
 
     #New idea.
     #
-    # In the presence and near of bedslope it is necessary to 
-    # limit slopes based on differences in h rather than w 
+    # In the presence and near of bedslope it is necessary to
+    # limit slopes based on differences in h rather than w
     # (which is what is needed away from the bed).
     #
@@ -733 +733 @@
     # it will need to be balanced with a h-limited gradient.
     #
-    # For this we will use the quantity alpha as before 
+    # For this we will use the quantity alpha as before
     #
-        
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     hc = wc - zc
-    
+
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
@@ -760 +760 @@
                  abs(zv[k,2]-zc[k]))
 
-        
+
         hmin = min( hv[k,:] )
 
-        #Create alpha in [0,1], where alpha==0 means using the h-limited 
+        #Create alpha in [0,1], where alpha==0 means using the h-limited
         #stage and alpha==1 means using the w-limited stage as
         #computed by the gradient limiter (both 1st or 2nd order)
-    
+
         #If hmin > dz/2 then alpha = 1 and the bed will have no effect
         #If hmin < 0 then alpha = 0 reverting to constant height above bed.
@@ -773 +773 @@
             alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
         else:
-            #Flat bed 
+            #Flat bed
             alpha = 1.0
 
-        #Let 
+        #Let
         #
-        #  wvi be the w-limited stage (wvi = zvi + hvi)        
+        #  wvi be the w-limited stage (wvi = zvi + hvi)
         #  wvi- be the h-limited state (wvi- = zvi + hvi-)
-        #  
-        #    
+        #
+        #
         #where i=0,1,2 denotes the vertex ids
         #
-        #Weighted balance between w-limited and h-limited stage is  
+        #Weighted balance between w-limited and h-limited stage is
         #
-        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
+        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
         #
         #It follows that the updated wvi is
@@ -793 +793 @@
         # Momentum is balanced between constant and limited
 
-        
+
         #for i in range(3):
         #    wv[k,i] = zv[k,i] + hvbar[k,i]
-        
+
         #return
-        
-        if alpha < 1:         
-        
+
+        if alpha < 1:
+
             for i in range(3):
                 wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
-            
+
             #Momentums at centroids
             xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values        
+            ymomc = domain.quantities['ymomentum'].centroid_values
 
             #Momentums at vertices
             xmomv = domain.quantities['xmomentum'].vertex_values
-            ymomv = domain.quantities['ymomentum'].vertex_values         
-
-            # Update momentum as a linear combination of 
+            ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
             # xmomc and ymomc (shallow) and momentum
             # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-                                    
-            
+
+
 def balance_deep_and_shallow_c(domain):
     """Wrapper for C implementation
     """
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     hc = wc - zc
-    
+
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
@@ -834 +834 @@
     #Momentums at centroids
     xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values        
+    ymomc = domain.quantities['ymomentum'].centroid_values
 
     #Momentums at vertices
     xmomv = domain.quantities['xmomentum'].vertex_values
-    ymomv = domain.quantities['ymomentum'].vertex_values         
+    ymomv = domain.quantities['ymomentum'].vertex_values
 
     #Limit h
-    hvbar = h_limiter(domain)    
-    
+    hvbar = h_limiter(domain)
+
     #This is how one would make a first order h_limited value
     #as in the old balancer (pre 17 Feb 2005):
@@ -852 +852 @@
     from shallow_water_ext import balance_deep_and_shallow
     balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
-                             xmomc, ymomc, xmomv, ymomv)    
-
-    
+                             xmomc, ymomc, xmomv, ymomv)
+
+
 
 
@@ -867 +867 @@
     and third conserved quantities.
     """
-    
-    def __init__(self, domain = None):        
+
+    def __init__(self, domain = None):
         Boundary.__init__(self)
 
@@ -874 +874 @@
             msg = 'Domain must be specified for reflective boundary'
             raise msg
-        
+
         #Handy shorthands
         self.stage = domain.quantities['stage'].edge_values
         self.xmom = domain.quantities['xmomentum'].edge_values
-        self.ymom = domain.quantities['ymomentum'].edge_values          
+        self.ymom = domain.quantities['ymomentum'].edge_values
         self.normals = domain.normals
-        
-        from Numeric import zeros, Float        
-        self.conserved_quantities = zeros(3, Float) 
-        
+
+        from Numeric import zeros, Float
+        self.conserved_quantities = zeros(3, Float)
+
     def __repr__(self):
         return 'Reflective_boundary'
@@ -895 +895 @@
         q = self.conserved_quantities
         q[0] = self.stage[vol_id, edge_id]
-        q[1] = self.xmom[vol_id, edge_id] 
-        q[2] = self.ymom[vol_id, edge_id] 
-        
-        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]     
-
-        
+        q[1] = self.xmom[vol_id, edge_id]
+        q[2] = self.ymom[vol_id, edge_id]
+
+        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+
+
         r = rotate(q, normal, direction = 1)
         r[1] = -r[1]
@@ -907 +907 @@
         return q
 
-        
+
 #class Spatio_temporal_boundary(Boundary):
-#    """The spatio-temporal boundary, reads values for the conserved 
-#    quantities from an sww NetCDF file, and returns interpolated values 
+#    """The spatio-temporal boundary, reads values for the conserved
+#    quantities from an sww NetCDF file, and returns interpolated values
 #    at the midpoints of each associated boundaty segment.
 #    Time dependency is interpolated linearly as in util.File_function.#
 #
 #    Example:
-#    Bf = Spatio_temporal_boundary('source_file.sww', domain) 
-#          
+#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
+#
 #    """
 Spatio_temporal_boundary = File_boundary
-        
-
-        
+
+
+
 
 #########################
@@ -937 +937 @@
 
     Stage = domain.quantities['stage']
-    Elevation = domain.quantities['elevation']    
+    Elevation = domain.quantities['elevation']
     h = Stage.edge_values - Elevation.edge_values
     v = Elevation.vertex_values
-    
+
     x = domain.get_vertex_coordinates()
     g = domain.g
-    
-    for k in range(domain.number_of_elements):    
+
+    for k in range(domain.number_of_elements):
         avg_h = sum( h[k,:] )/3
-        
+
         #Compute bed slope
-        x0, y0, x1, y1, x2, y2 = x[k,:]    
+        x0, y0, x1, y1, x2, y2 = x[k,:]
         z0, z1, z2 = v[k,:]
 
@@ -955 +955 @@
         #Update momentum
         xmom[k] += -g*zx*avg_h
-        ymom[k] += -g*zy*avg_h        
+        ymom[k] += -g*zy*avg_h
 
 
 def gravity_c(domain):
     """Wrapper calling C version
-    """    
+    """
 
     xmom = domain.quantities['xmomentum'].explicit_update
     ymom = domain.quantities['ymomentum'].explicit_update
 
-    Stage = domain.quantities['stage']    
-    Elevation = domain.quantities['elevation']    
+    Stage = domain.quantities['stage']
+    Elevation = domain.quantities['elevation']
     h = Stage.edge_values - Elevation.edge_values
     v = Elevation.vertex_values
-    
+
     x = domain.get_vertex_coordinates()
     g = domain.g
-    
-    
+
+
     from shallow_water_ext import gravity
     gravity(g, h, v, x, xmom, ymom)
-    
-        
+
+
 def manning_friction(domain):
     """Apply (Manning) friction to water momentum
@@ -987 +987 @@
     z = domain.quantities['elevation'].centroid_values
     h = w-z
-    
+
     uh = domain.quantities['xmomentum'].centroid_values
     vh = domain.quantities['ymomentum'].centroid_values
-    eta = domain.quantities['friction'].centroid_values    
-    
+    eta = domain.quantities['friction'].centroid_values
+
     xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    ymom_update = domain.quantities['ymomentum'].semi_implicit_update    
-    
+    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
     N = domain.number_of_elements
     eps = domain.minimum_allowed_height
     g = domain.g
-    
+
     for k in range(N):
         if eta[k] >= eps:
@@ -1008 +1008 @@
                 xmom_update[k] += S*uh[k]
                 ymom_update[k] += S*vh[k]
-            
-        
+
+
 def manning_friction_c(domain):
     """Wrapper for c version
@@ -1016 +1016 @@
 
     xmom = domain.quantities['xmomentum']
-    ymom = domain.quantities['ymomentum']    
-        
+    ymom = domain.quantities['ymomentum']
+
     w = domain.quantities['stage'].centroid_values
     z = domain.quantities['elevation'].centroid_values
-    
+
     uh = xmom.centroid_values
     vh = ymom.centroid_values
-    eta = domain.quantities['friction'].centroid_values    
-    
+    eta = domain.quantities['friction'].centroid_values
+
     xmom_update = xmom.semi_implicit_update
-    ymom_update = ymom.semi_implicit_update    
-    
+    ymom_update = ymom.semi_implicit_update
+
     N = domain.number_of_elements
     eps = domain.minimum_allowed_height
@@ -1050 +1050 @@
     uh = domain.quantities['xmomentum'].centroid_values
     vh = domain.quantities['ymomentum'].centroid_values
-    tau = domain.quantities['linear_friction'].centroid_values    
-    
+    tau = domain.quantities['linear_friction'].centroid_values
+
     xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    ymom_update = domain.quantities['ymomentum'].semi_implicit_update    
-    
+    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
+
     N = domain.number_of_elements
     eps = domain.minimum_allowed_height
     g = domain.g #Not necessary? Why was this added?
-    
+
     for k in range(N):
         if tau[k] >= eps:
@@ -1067 +1067 @@
                 xmom_update[k] += S*uh[k]
                 ymom_update[k] += S*vh[k]
-            
-        
-        
+
+
+
 def check_forcefield(f):
     """Check that f is either
@@ -1075 +1075 @@
        and that it returns an array or a list of same length
        as x and y
-    2: a scalar   
+    2: a scalar
     """
 
     from Numeric import ones, Float, array
 
-    
+
     if callable(f):
         N = 3
         x = ones(3, Float)
-        y = ones(3, Float)            
+        y = ones(3, Float)
         try:
             q = f(1.0, x, y)
@@ -1091 +1091 @@
             #FIXME: Reconsider this semantics
             raise msg
-                    
+
         try:
             q = array(q).astype(Float)
@@ -1097 +1097 @@
             msg = 'Return value from vector function %s could ' %f
             msg += 'not be converted into a Numeric array of floats.\n'
-            msg += 'Specified function should return either list or array.' 
+            msg += 'Specified function should return either list or array.'
             raise msg
 
@@ -1103 +1103 @@
         msg += 'must have same lenght as input vectors'
         assert len(q) == N, msg
-        
-    else:        
+
+    else:
         try:
             f = float(f)
@@ -1111 +1111 @@
             msg += ' or a vector function'
             raise msg
-    return f    
+    return f
 
 
 class Wind_stress:
     """Apply wind stress to water momentum in terms of
-    wind speed [m/s] and wind direction [degrees]    
+    wind speed [m/s] and wind direction [degrees]
     """
 
@@ -1122 +1122 @@
         """Initialise windfield from wind speed s [m/s]
         and wind direction phi [degrees]
-        
+
         Inputs v and phi can be either scalars or Python functions, e.g.
 
@@ -1128 +1128 @@
 
         #FIXME - 'normal' degrees are assumed for now, i.e. the
-        vector (1,0) has zero degrees. 
+        vector (1,0) has zero degrees.
         We may need to convert from 'compass' degrees later on and also
         map from True north to grid north.
@@ -1137 +1137 @@
             ...
             return s
-        
+
         def angle(t,x,y):
             ...
@@ -1154 +1154 @@
         can be applied, providing both quantities simultaneously.
         As in
-        W = Wind_stress(F), where returns (speed, angle) for each t.         
+        W = Wind_stress(F), where returns (speed, angle) for each t.
 
         domain.forcing_terms.append(W)
@@ -1160 +1160 @@
 
         from config import rho_a, rho_w, eta_w
-        from Numeric import array, Float                
+        from Numeric import array, Float
 
         if len(args) == 2:
@@ -1172 +1172 @@
         else:
            #Assume info is in 2 keyword arguments
-        
+
            if len(kwargs) == 2:
                s = kwargs['s']
                phi = kwargs['phi']
            else:
-               raise 'Assumes two keyword arguments: s=..., phi=....'   
-            
+               raise 'Assumes two keyword arguments: s=..., phi=....'
+
         self.speed = check_forcefield(s)
         self.phi = check_forcefield(phi)
-            
+
         self.const = eta_w*rho_a/rho_w
-        
+
 
     def __call__(self, domain):
@@ -1190 +1190 @@
 
         from math import pi, cos, sin, sqrt
-        from Numeric import Float, ones, ArrayType 
-        
+        from Numeric import Float, ones, ArrayType
+
         xmom_update = domain.quantities['xmomentum'].explicit_update
-        ymom_update = domain.quantities['ymomentum'].explicit_update    
-        
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
         N = domain.number_of_elements
-        t = domain.time        
-        
+        t = domain.time
+
         if callable(self.speed):
-            xc = domain.get_centroid_coordinates()            
+            xc = domain.get_centroid_coordinates()
             s_vec = self.speed(t, xc[:,0], xc[:,1])
         else:
@@ -1212 +1212 @@
 
         if callable(self.phi):
-            xc = domain.get_centroid_coordinates()            
+            xc = domain.get_centroid_coordinates()
             phi_vec = self.phi(t, xc[:,0], xc[:,1])
         else:
@@ -1218 +1218 @@
 
             try:
-                phi_vec = self.phi * ones(N, Float)            
+                phi_vec = self.phi * ones(N, Float)
             except:
                 msg = 'Angle must be either callable or a scalar: %s' %self.phi
@@ -1233 +1233 @@
     """
     from math import pi, cos, sin, sqrt
-    
+
     N = len(s_vec)
     for k in range(N):
@@ -1245 +1245 @@
         u = s*cos(phi)
         v = s*sin(phi)
-        
+
         #Compute wind stress
         S = const * sqrt(u**2 + v**2)
         xmom_update[k] += S*u
-        ymom_update[k] += S*v        
-            
+        ymom_update[k] += S*v
+
 
 
 ##############################
 #OBSOLETE STUFF
-    
+
 def balance_deep_and_shallow_old(domain):
     """Compute linear combination between stage as computed by
@@ -1268 +1268 @@
 
     #OBSOLETE
-    
+
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
     zc = domain.quantities['elevation'].centroid_values
     hc = wc - zc
-    
+
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
@@ -1280 +1280 @@
 
     #Computed linear combination between constant stages and and
-    #stages parallel to the bed elevation.     
+    #stages parallel to the bed elevation.
     for k in range(domain.number_of_elements):
         #Compute maximal variation in bed elevation
@@ -1293 +1293 @@
                  abs(zv[k,2]-zc[k]))
 
-        
+
         hmin = min( hv[k,:] )
 
-        #Create alpha in [0,1], where alpha==0 means using shallow 
+        #Create alpha in [0,1], where alpha==0 means using shallow
         #first order scheme and alpha==1 means using the stage w as
         #computed by the gradient limiter (1st or 2nd order)
@@ -1302 +1302 @@
         #If hmin > dz/2 then alpha = 1 and the bed will have no effect
         #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-      
+
         if dz > 0.0:
             alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
         else:
-            #Flat bed 
+            #Flat bed
             alpha = 1.0
 
-        #Weighted balance between stage parallel to bed elevation 
-        #(wvi = zvi + hc) and stage as computed by 1st or 2nd 
+        #Weighted balance between stage parallel to bed elevation
+        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
         #order gradient limiter
         #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
         #
         #It follows that the updated wvi is
-        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
+        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
         #  zvi + hc + alpha*(hvi - hc)
         #
         #Note that hvi = zc+hc-zvi in the first order case (constant).
 
-        if alpha < 1:         
+        if alpha < 1:
             for i in range(3):
                 wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-            
+
 
             #Momentums at centroids
             xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values        
+            ymomc = domain.quantities['ymomentum'].centroid_values
 
             #Momentums at vertices
             xmomv = domain.quantities['xmomentum'].vertex_values
-            ymomv = domain.quantities['ymomentum'].vertex_values         
-
-            # Update momentum as a linear combination of 
+            ymomv = domain.quantities['ymomentum'].vertex_values
+
+            # Update momentum as a linear combination of
             # xmomc and ymomc (shallow) and momentum
             # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-            
+
 
 
@@ -1353 +1353 @@
     x,y are assumed to be in the unit square
     """
-    
+
     def __init__(self, stage):
         self.inflow_stage = stage
 
-    def __call__(self, x, y):    
+    def __call__(self, x, y):
         from Numeric import zeros, Float
         from math import sqrt
-    
+
         N = len(x)
-        assert N == len(y)    
+        assert N == len(y)
 
         z = zeros(N, Float)
@@ -1370 +1370 @@
             #Flattish bit to the left
             if x[i] < 0.3:
-                z[i] = -x[i]/10 
-            
+                z[i] = -x[i]/10
+
             #Weir
             if x[i] >= 0.3 and x[i] < 0.4:
@@ -1381 +1381 @@
             depth = -1.0
             #plateaux = -0.9
-            plateaux = -0.6            
+            plateaux = -0.6
             if y[i] < 0.7:
                 if x[i] > x0 and x[i] < 0.9:
@@ -1389 +1389 @@
                 if x[i] >= 0.9:
                     z[i] = plateaux
-                    
-                    
+
+
             elif y[i] >= 0.7 and y[i] < 1.5:
                 #Restrict and deepen
@@ -1397 +1397 @@
                     #z[i] = depth-y[i]/5
                     #z[i] = depth
-                elif x[i] >= 0.8:    
+                elif x[i] >= 0.8:
                     #RHS plateaux
                     z[i] = plateaux
-                    
+
             elif y[i] >= 1.5:
                 if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
                     #Widen up and stay at constant depth
                     z[i] = depth-1.5/5
-                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:                        
+                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
                     #RHS plateaux
-                    z[i] = plateaux                                    
-                    
+                    z[i] = plateaux
+
 
             #Hole in weir (slightly higher than inflow condition)
@@ -1418 +1418 @@
             if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
                 z[i] = -x[i]+self.inflow_stage + 0.02
-                
+
             if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
                 #Flatten it out towards the end
@@ -1432 +1432 @@
             if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
                 z[i] = -0.9 #North south
-                
+
             if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
                 z[i] = -1.0 + (x[i]-0.9)/3 #East west
-            
-            
+
+
 
             #Stuff not in use
-            
+
             #Upward slope at inlet to the north west
             #if x[i] < 0.0: # and y[i] > 0.5:
@@ -1460 +1460 @@
     x,y are assumed to be in the unit square
     """
-    
+
     def __init__(self, stage):
         self.inflow_stage = stage
 
-    def __call__(self, x, y):    
-        from Numeric import zeros, Float 
-    
+    def __call__(self, x, y):
+        from Numeric import zeros, Float
+
         N = len(x)
-        assert N == len(y)    
+        assert N == len(y)
 
         z = zeros(N, Float)
@@ -1477 +1477 @@
             if x[i] < 0.3:
                 z[i] = -x[i]/10  #General slope
-            
+
             #Weir
             if x[i] > 0.3 and x[i] < 0.4:
@@ -1488 +1488 @@
             #Hole in weir (slightly higher than inflow condition)
             if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i]+self.inflow_stage + 0.05       
-            
+                z[i] = -x[i]+self.inflow_stage + 0.05
+
 
         return z/2
-        
-
-            
+
+
+
 class Constant_height:
     """Set an initial condition with constant water height, e.g
@@ -1521 +1521 @@
     Fluxes across each edge are scaled by edgelengths and summed up
     Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities 
-    stage, xmomentum and ymomentum    
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
 
     The maximal allowable speed computed by the flux_function for each volume
@@ -1530 +1530 @@
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes. 
+      domain.timestep is set to the largest step satisfying all volumes.
     """
 
@@ -1537 +1537 @@
 
     N = domain.number_of_elements
-    
+
     #Shortcuts
     Stage = domain.quantities['stage']
     Xmom = domain.quantities['xmomentum']
     Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']        
+    Bed = domain.quantities['elevation']
 
     #Arrays
@@ -1548 +1548 @@
     xmom =  Xmom.edge_values
     ymom =  Ymom.edge_values
-    bed =   Bed.edge_values    
+    bed =   Bed.edge_values
 
     stage_bdry = Stage.boundary_values
     xmom_bdry =  Xmom.boundary_values
     ymom_bdry =  Ymom.boundary_values
-    
+
     flux = zeros((N,3), Float) #Work array for summing up fluxes
 
     #Loop
-    timestep = float(sys.maxint)    
+    timestep = float(sys.maxint)
     for k in range(N):
 
@@ -1566 +1566 @@
 
             #Quantities at neighbour on nearest face
-            n = domain.neighbours[k,i] 
+            n = domain.neighbours[k,i]
             if n < 0:
                 m = -n-1 #Convert negative flag to index
                 qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
                 zr = zl #Extend bed elevation to boundary
-            else:    
+            else:
                 m = domain.neighbour_edges[k,i]
                 qr = [stage[n, m], xmom[n, m], ymom[n, m]]
-                zr = bed[n, m]                
-
-                
-            #Outward pointing normal vector   
+                zr = bed[n, m]
+
+
+            #Outward pointing normal vector
             normal = domain.normals[k, 2*i:2*i+2]
 
             #Flux computation using provided function
             edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-        
+
             flux[k,:] = edgeflux
 
     return flux
-    
-  
+
+
 
 
@@ -1600 +1600 @@
 if compile.can_use_C_extension('shallow_water_ext.c'):
     #Replace python version with c implementations
-    
+
     from shallow_water_ext import rotate, assign_windfield_values
     compute_fluxes = compute_fluxes_c
@@ -1609 +1609 @@
     protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
 
-    
+
     #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
 
@@ -1631 +1631 @@
         psyco.bind(Domain.distribute_to_vertices_and_edges)
         psyco.bind(Domain.compute_fluxes)
-        
+
 if __name__ == "__main__":
     pass