Changeset 2716

Timestamp:

Apr 13, 2006, 4:34:22 PM (17 years ago)

Author:

jakeman

Message:

Adding new files

Location:

development/pyvolution-1d

Files:

• domain.py (modified) (21 diffs)
• generic_boundary_conditions.py (added)
• quantity.py (added)
• shallow_water_1d.py (modified) (18 diffs)
• shallow_water_ext.c (added)
• test_all.py (added)
• test_quantity.py (modified) (2 diffs)
• test_shallow_water_1d.py (added)
• util.py (added)
• util_ext.c (added)
• util_ext.h (added)

development/pyvolution-1d/domain.py

@@ -8 +8 @@
 """
 from generic_boundary_conditions import *
+from coordinate_transforms.geo_reference import Geo_reference
 
 class Domain:
 
     def __init__(self, coordinates, boundary = None,
-                 conserved_quantities = None,
-                 other_quantities = None,
-                 tagged_elements = None):
+                 conserved_quantities = None, other_quantities = None,
+                 tagged_elements = None, geo_reference = None):
         """
         Build 1D elements from x coordinates
@@ -25 +25 @@
         self.coordinates = array(coordinates)
 
+        if geo_reference is None:
+            self.geo_reference = Geo_reference() #Use defaults 
+        else:
+            self.geo_reference = geo_reference
+
         #Register number of Elements
         self.number_of_elements = N = len(self.coordinates)-1
@@ -30 +35 @@
         #Allocate space for neighbour and boundary structures
         self.neighbours = zeros((N, 2), Int)
-        self.neighbour_edges = zeros((N, 2), Int)
+        #self.neighbour_edges = zeros((N, 2), Int)
+        self.neighbour_vertices = zeros((N, 2), Int)
         self.number_of_boundaries = zeros(N, Int)
         self.surrogate_neighbours = zeros((N, 2), Int)
@@ -60 +66 @@
             self.neighbours[i, :] = [-1, -1]
             #Initialise edge ids of neighbours
+            #Initialise vertex ids of neighbours
             #In case of boundaries this slot is not used
-            self.neighbour_edges[i, :] = [-1, -1]
+            #self.neighbour_edges[i, :] = [-1, -1]
+            self.neighbour_vertices[i, :] = [-1, -1]
+
+        self.build_vertexlist()
 
         #Build neighbour structure
@@ -70 +80 @@
 
         #Build boundary dictionary mapping (id, edge) to symbolic tags
+        #Build boundary dictionary mapping (id, vertex) to symbolic tags
         self.build_boundary_dictionary(boundary)
 
         #Build tagged element  dictionary mapping (tag) to array of elements
         self.build_tagged_elements_dictionary(tagged_elements)
-
-
         
         from quantity import Quantity, Conserved_quantity
@@ -105 +114 @@
         self.forcing_terms = []
 
+        #Defaults
+        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
+        self.beta_w = beta_w
+        self.beta_h = beta_h
+        self.epsilon = epsilon
+
+        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
+        #Or maybe get rid of order altogether and use beta_w and beta_h
+        self.default_order = 1
+        self.order = self.default_order
+        self.smallsteps = 0
+        self.max_smallsteps = max_smallsteps
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+        self.CFL = CFL
+
+        #Model time
+        self.time = 0.0
+        self.finaltime = None
+        self.min_timestep = self.max_timestep = 0.0
+        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        #Checkpointing and storage
+        from config import default_datadir
+        self.datadir = default_datadir
+        self.filename = 'domain'
+        self.checkpoint = False
+
+    def __len__(self):
+        return self.number_of_elements
+
+    def build_vertexlist(self):
+        """Build vertexlist index by vertex ids and for each entry (point id)
+        build a list of (triangles, vertex_id) pairs that use the point
+        as vertex.
+
+        Preconditions:
+          self.coordinates and self.triangles are defined
+
+        Postcondition:
+          self.vertexlist is built
+        """
+        from Numeric import array
+
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(self.number_of_elements):
+
+            #a = self.triangles[i, 0]
+            #b = self.triangles[i, 1]
+            #c = self.triangles[i, 2]
+            a = i
+            b = i + 1
+
+            #Register the vertices v as lists of
+            #(triangle_id, vertex_id) tuples associated with them
+            #This is used for smoothing
+            #for vertex_id, v in enumerate([a,b,c]):
+            for vertex_id, v in enumerate([a,b]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )
+
+        self.vertexlist = vertexlist
+
+        
     def build_neighbour_structure(self):
         """Update all registered triangles to point to their neighbours.
@@ -112 +186 @@
         Postconditions:
           neighbours and neighbour_edges is populated
+          neighbours and neighbour_vertices is populated
           number_of_boundaries integer array is defined.
         """
@@ -121 +196 @@
         N = self.number_of_elements
         neighbourdict = {}
-        l_edge = 0
-        r_edge = 1
+        #l_edge = 0
+        #r_edge = 1
+        l_vertex = 0
+        r_vertex = 1
         for i in range(N):
 
@@ -149 +226 @@
             #neighbourdict[a,b] = (i, 1) #(id, edge)
             #neighbourdict[b,a] = (i, 0) #(id, edge)
-            neighbourdict[a,l_edge] = (i, 0) #(id, edge)
-            neighbourdict[b,r_edge] = (i, 1) #(id, edge)
+            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
+            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
+            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
+            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
 
 
@@ -166 +245 @@
             #self.number_of_boundaries[i] = 3
             self.number_of_boundaries[i] = 2
-            if neighbourdict.has_key((b,l_edge)):
-                self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
-                self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
+            #if neighbourdict.has_key((b,l_edge)):
+            if neighbourdict.has_key((b,l_vertex)):
+                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
+                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
+                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
+                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
                 self.number_of_boundaries[i] -= 1
 
-            if neighbourdict.has_key((a,r_edge)):
-                self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
-                self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
+            #if neighbourdict.has_key((a,r_edge)):
+            if neighbourdict.has_key((a,r_vertex)):
+                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
+                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
+                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
+                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
                 self.number_of_boundaries[i] -= 1
+                
             #if neighbourdict.has_key((b,a)):
             #    self.neighbours[i, 1] = neighbourdict[b,a][0]
@@ -231 +317 @@
             for vol_id in range(self.number_of_elements):
                 #for edge_id in range(0, 3):
-                for edge_id in range(0, 2):
-                    if self.neighbours[vol_id, edge_id] < 0:
-                        boundary[(vol_id, edge_id)] = default_boundary_tag
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:   
+                        #boundary[(vol_id, edge_id)] = default_boundary_tag
+                        boundary[(vol_id, vertex_id)] = default_boundary_tag
         else:
             #Check that all keys in given boundary exist
-            for vol_id, edge_id in boundary.keys():
-                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+            #for vol_id, edge_id in boundary.keys():
+            for vol_id, vertex_id in boundary.keys():    
+                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
                 a, b = self.neighbours.shape
-                assert vol_id < a and edge_id < b, msg
+                #assert vol_id < a and edge_id < b, msg
+                assert vol_id < a and vertex_id < b, msg
 
                 #FIXME: This assert violates internal boundaries (delete it)
@@ -248 +340 @@
             for vol_id in range(self.number_of_elements):
                 #for edge_id in range(0, 3):
-                for edge_id in range(0, 2):    
-                    if self.neighbours[vol_id, edge_id] < 0:
-                        if not boundary.has_key( (vol_id, edge_id) ):
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:    
+                        #if not boundary.has_key( (vol_id, edge_id) ):
+                        if not boundary.has_key( (vol_id, vertex_id) ):    
                             msg = 'WARNING: Given boundary does not contain '
-                            msg += 'tags for edge (%d, %d). '\
-                                   %(vol_id, edge_id)
+                            #msg += 'tags for edge (%d, %d). '\
+                            #       %(vol_id, edge_id)
+                            msg += 'tags for vertex (%d, %d). '\
+                                   %(vol_id, vertex_id)
                             msg += 'Assigning default tag (%s).'\
                                    %default_boundary_tag
@@ -264 +361 @@
                             #enable default boundary-tags where
                             #tags a not specified
-                            boundary[ (vol_id, edge_id) ] =\
+                            #boundary[ (vol_id, edge_id) ] =\
+                            boundary[ (vol_id, vertex_id) ] =\
                                       default_boundary_tag
 
@@ -304 +402 @@
         return tags.keys()
         
-
-    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
+    def get_vertex_coordinates(self, obj = False):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+
+        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
+        FIXME, we might make that the default.
+        FIXME Maybe use keyword: continuous = False for this condition?
+
+        
+        """
+
+        if obj is True:
+            from Numeric import concatenate, reshape
+            #V = self.vertex_coordinates
+            V = self.vertices
+            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
+
+            N = V.shape[0]
+            #return reshape(V, (3*N, 2))
+            return reshape(V, (N, 2))
+        else:    
+            #return self.vertex_coordinates
+            return self.vertices
+    
+    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
         """Get conserved quantities at volume vol_id
 
@@ -319 +441 @@
         from Numeric import zeros, Float
 
-        if not (vertex is None or edge is None):
-            msg = 'Values for both vertex and edge was specified.'
-            msg += 'Only one (or none) is allowed.'
-            raise msg
+        #if not (vertex is None):# or edge is None):
+        #    msg = 'Values for both vertex and edge was specified.'
+        #    msg += 'Only one (or none) is allowed.'
+       #     raise msg
 
         q = zeros( len(self.conserved_quantities), Float)
@@ -330 +452 @@
             if vertex is not None:
                 q[i] = Q.vertex_values[vol_id, vertex]
-            elif edge is not None:
-                q[i] = Q.edge_values[vol_id, edge]
+            #elif edge is not None:
+            #    q[i] = Q.edge_values[vol_id, edge]
             else:
                 q[i] = Q.centroid_values[vol_id]
@@ -368 +490 @@
 
         return self.areas[elem_id]
+
+    def get_quantity(self, name, location='vertices', indices = None):
+        """Get values for named quantity
+
+        name: Name of quantity
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3.
+
+        Indices is the set of element ids that the operation applies to.
+
+        The values will be stored in elements following their
+        internal ordering.
+        """
+
+        return self.quantities[name].get_values( location, indices = indices)
 
     def set_quantity(self, name, *args, **kwargs):
@@ -462 +601 @@
         #Loop through edges that lie on the boundary and associate them with
         #callable boundary objects depending on their tags
-        for k, (vol_id, edge_id) in enumerate(x):
-            tag = self.boundary[ (vol_id, edge_id) ]
+        #for k, (vol_id, edge_id) in enumerate(x):
+        for k, (vol_id, vertex_id) in enumerate(x):
+            #tag = self.boundary[ (vol_id, edge_id) ]
+            tag = self.boundary[ (vol_id, vertex_id) ]
 
             if boundary_map.has_key(tag):
@@ -469 +610 @@
 
                 if B is not None:
-                    self.boundary_objects.append( ((vol_id, edge_id), B) )
-                    self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
+                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
+                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
+                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
+                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
                 else:
                     pass
@@ -494 +637 @@
         ##assert hasattr(self, 'boundary_objects')
 
+    def get_name(self):
+        return self.filename
+
+    def set_name(self, name):
+        self.filename = name
+
+    def get_datadir(self):
+        return self.datadir
+
+    def set_datadir(self, name):
+        self.datadir = name
+
+#Main components of evolve
+
+    def evolve(self, yieldstep = None, finaltime = None,
+               skip_initial_step = False):
+        """Evolve model from time=0.0 to finaltime yielding results
+        every yieldstep.
+
+        Internally, smaller timesteps may be taken.
+
+        Evolve is implemented as a generator and is to be called as such, e.g.
+
+        for t in domain.evolve(timestep, yieldstep, finaltime):
+            <Do something with domain and t>
+
+        """
+
+        from config import min_timestep, max_timestep, epsilon
+
+        #FIXME: Maybe lump into a larger check prior to evolving
+        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
+        msg += 'e.g. using the method set_boundary.\n'
+        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
+        assert hasattr(self, 'boundary_objects'), msg
+
+        ##self.set_defaults()
+
+        if yieldstep is None:
+            yieldstep = max_timestep
+        else:
+            yieldstep = float(yieldstep)
+
+        self.order = self.default_order
+
+
+        self.yieldtime = 0.0 #Time between 'yields'
+
+        #Initialise interval of timestep sizes (for reporting only)
+        self.min_timestep = max_timestep
+        self.max_timestep = min_timestep
+        self.finaltime = finaltime
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+
+        #update ghosts
+        #self.update_ghosts()
+
+        #Initial update of vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Initial update boundary values
+        self.update_boundary()
+
+        #Or maybe restore from latest checkpoint
+        if self.checkpoint is True:
+            self.goto_latest_checkpoint()
+
+        if skip_initial_step is False:
+            yield(self.time)  #Yield initial values
+
+        while True:
+
+            #Compute fluxes across each element edge
+            self.compute_fluxes()
+
+            #Update timestep to fit yieldstep and finaltime
+            self.update_timestep(yieldstep, finaltime)
+
+            #Update conserved quantities
+            self.update_conserved_quantities()
+
+            #update ghosts
+            #self.update_ghosts()
+
+            #Update vertex and edge values
+            self.distribute_to_vertices_and_edges()
+
+            #Update boundary values
+            self.update_boundary()
+
+            #Update time
+            self.time += self.timestep
+            self.yieldtime += self.timestep
+            self.number_of_steps += 1
+            if self.order == 1:
+                self.number_of_first_order_steps += 1
+
+            #Yield results
+            if finaltime is not None and abs(self.time - finaltime) < epsilon:
+
+                #FIXME: There is a rare situation where the
+                #final time step is stored twice. Can we make a test?
+
+                # Yield final time and stop
+                yield(self.time)
+                break
+
+
+            if abs(self.yieldtime - yieldstep) < epsilon:
+                # Yield (intermediate) time and allow inspection of domain
+
+                if self.checkpoint is True:
+                    self.store_checkpoint()
+                    self.delete_old_checkpoints()
+
+                #Pass control on to outer loop for more specific actions
+                yield(self.time)
+
+                # Reinitialise
+                self.yieldtime = 0.0
+                self.min_timestep = max_timestep
+                self.max_timestep = min_timestep
+                self.number_of_steps = 0
+                self.number_of_first_order_steps = 0
+
+    def distribute_to_vertices_and_edges(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume
+
+        Default implementation is straight first order,
+        i.e. constant values throughout each element and
+        no reference to non-conserved quantities.
+        """
+
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            if self.order == 1:
+                Q.extrapolate_first_order()
+            elif self.order == 2:
+                Q.extrapolate_second_order()
+                Q.limit()
+            else:
+                raise 'Unknown order'
+            Q.interpolate_from_vertices_to_edges()
+
+
     def update_boundary(self):
         """Go through list of boundary objects and update boundary values
@@ -501 +791 @@
         #FIXME: Update only those that change (if that can be worked out)
         #FIXME: Boundary objects should not include ghost nodes.
-        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
-            q = B.evaluate(vol_id, edge_id)
+        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+        #    q = B.evaluate(vol_id, edge_id)
+        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
+            q = B.evaluate(vol_id, vertex_id)
 
             for j, name in enumerate(self.conserved_quantities):
                 Q = self.quantities[name]
                 Q.boundary_values[i] = q[j]
+
+    def compute_forcing_terms(self):
+        """If there are any forcing functions driving the system
+        they should be defined in Domain subclass and appended to
+        the list self.forcing_terms
+        """
+
+        for f in self.forcing_terms:
+            f(self)
 
 

development/pyvolution-1d/shallow_water_1d.py

@@ -49 +49 @@
 class Domain(Generic_Domain):
 
-    def __init__(self, coordinates, boundary = None, tagged_elements = None):
+    def __init__(self, coordinates, boundary = None, tagged_elements = None,
+                 geo_reference = None):
 
         conserved_quantities = ['stage', 'xmomentum']
@@ -55 +56 @@
         Generic_Domain.__init__(self, coordinates, boundary,
                                 conserved_quantities, other_quantities,
-                                tagged_elements)
+                                tagged_elements, geo_reference)
         
         from config import minimum_allowed_height, g
@@ -281 +282 @@
 # Flux computation
 #def flux_function(normal, ql, qr, zl, zr):
-def flux_function(ql, qr, zl, zr):
+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
@@ -305 +306 @@
     #q_right = rotate(qr, normal, direction = 1)
     q_left = ql
+    q_left[1] = q_left[1]*normal
     q_right = qr
+    q_right[1] = q_right[1]*normal
 
     z = (zl+zr)/2 #Take average of field values
@@ -338 +341 @@
 
     #Maximal wave speed
-    #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right,0)
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
     
     #Minimal wave speed
-    #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right,0)
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
     
     #Flux computation
@@ -363 +364 @@
     denom = s_max-s_min
     if denom == 0.0:
-        edgeflux = array([0.0, 0.0, 0.0])
+        #edgeflux = array([0.0, 0.0, 0.0])
         edgeflux = array([0.0, 0.0])
         max_speed = 0.0
@@ -369 +370 @@
         edgeflux = (s_max*flux_left - s_min*flux_right)/denom
         edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = -edgeflux[1]*normal
 
         #edgeflux = rotate(edgeflux, normal, direction=-1)
@@ -407 +410 @@
 
     #Arrays
-    stage = Stage.edge_values
-    xmom =  Xmom.edge_values
+    #stage = Stage.edge_values
+    #xmom =  Xmom.edge_values
  #   ymom =  Ymom.edge_values
-    bed =   Bed.edge_values
+    #bed =   Bed.edge_values
+    
+    stage = Stage.vertex_values
+    xmom =  Xmom.vertex_values
+    bed =   Bed.vertex_values
+    
+    print 'stage edge values'
+    print stage
+    print 'xmom edge values'
+    print xmom
+    print 'bed values'
+    print bed
 
     stage_bdry = Stage.boundary_values
     xmom_bdry =  Xmom.boundary_values
+    print 'stage_bdry'
+    print stage_bdry
+    print 'xmom_bdry'
+    print xmom_bdry
 #    ymom_bdry =  Ymom.boundary_values
 
@@ -439 +457 @@
                 zr = zl #Extend bed elevation to boundary
             else:
-                m = domain.neighbour_edges[k,i]
+                #m = domain.neighbour_edges[k,i]
+                m = domain.neighbour_vertices[k,i]
                 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                 qr = [stage[n, m], xmom[n, m]]
@@ -446 +465 @@
 
             #Outward pointing normal vector
-            normal = domain.normals[k, 2*i:2*i+2]
+            #normal = domain.normals[k, 2*i:2*i+2]
             #CHECK THIS LINE [k, 2*i:2*i+1]
 
             #Flux computation using provided function
             #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            edgeflux, max_speed = flux_function(ql, qr, zl, zr)
-            flux -= edgeflux * domain.edgelengths[k,i]
+            print 'ql'
+            print ql
+            print 'qr'
+            print qr
+            
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+            print 'edgeflux'
+            print edgeflux
+            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
+            # flux = edgefluxleft - edgefluxright
+            flux -= edgeflux #* domain.edgelengths[k,i]
 
             #Update optimal_timestep
             try:
-                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
+                timestep = min(timestep, 0.5/max_speed)
+                print 'check time step in compute fluxes is ok'
             except ZeroDivisionError:
                 pass
@@ -514 +544 @@
     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
                                      domain.neighbours,
-                                     domain.neighbour_edges,
+                                     #domain.neighbour_edges,
+                                     domain.neighbour_vertices,
                                      domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
+                                     #domain.edgelengths,
+                                     #domain.radii,
                                      domain.areas,
-                                     Stage.edge_values,
-                                     Xmom.edge_values,
+                                     #Stage.edge_values,
+                                     #Xmom.edge_values,
                                      #Ymom.edge_values,
-                                     Bed.edge_values,
+                                     #Bed.edge_values,
+                                     Stage.vertex_values,
+                                     Xmom.vertex_values,
+                                     Bed.vertex_values,
                                      Stage.boundary_values,
                                      Xmom.boundary_values,
@@ -568 +602 @@
         #MH090605 if second order,
         #perform the extrapolation and limiting on
-        #all of the conserved quantitie
+        #all of the conserved quantities
 
         if (domain.order == 1):
@@ -595 +629 @@
 
     #Compute edge values by interpolation
-    for name in domain.conserved_quantities:
-        Q = domain.quantities[name]
-        Q.interpolate_from_vertices_to_edges()
+    #for name in domain.conserved_quantities:
+    #    Q = domain.quantities[name]
+    #    Q.interpolate_from_vertices_to_edges()
 
 
@@ -898 +932 @@
 
         #Handy shorthands
-        self.stage   = domain.quantities['stage'].edge_values
-        self.xmom    = domain.quantities['xmomentum'].edge_values
+        #self.stage   = domain.quantities['stage'].edge_values
+        #self.xmom    = domain.quantities['xmomentum'].edge_values
         #self.ymom    = domain.quantities['ymomentum'].edge_values
         #self.normals = domain.normals
+        self.stage   = domain.quantities['stage'].vertex_values
+        self.xmom    = domain.quantities['xmomentum'].vertex_values
 
         from Numeric import zeros, Float
@@ -931 +967 @@
         r[1] = -q[1]
         q = r
+        #For start interval there is no outward momentum so do not need to
+        #reverse direction in this case
 
         return q
@@ -950 +988 @@
     Stage = domain.quantities['stage']
     Elevation = domain.quantities['elevation']
-    h = Stage.edge_values - Elevation.edge_values
+    #h = Stage.edge_values - Elevation.edge_values
+    h = Stage.vertex_values - Elevation.vertex_values
     v = Elevation.vertex_values
 
@@ -1005 +1044 @@
 
     uh = domain.quantities['xmomentum'].centroid_values
-    vh = domain.quantities['ymomentum'].centroid_values
+    #vh = domain.quantities['ymomentum'].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
@@ -1018 +1057 @@
         if eta[k] >= eps:
             if h[k] >= eps:
-                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+                S = -g * eta[k]**2 * uh[k]
                 S /= h[k]**(7.0/3)
 
                 #Update momentum
                 xmom_update[k] += S*uh[k]
-                ymom_update[k] += S*vh[k]
+                #ymom_update[k] += S*vh[k]
 
 

development/pyvolution-1d/test_quantity.py

@@ -112 +112 @@
             raise 'should have raised AssertionError'
 
+
     def test_set_values_const(self):
         quantity = Quantity(self.domain5)
@@ -120 +121 @@
         assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
 
-
         quantity.set_values(2.0, location = 'centroids')
         assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
-
 
     def test_set_values_func(self):