Changeset 7823

Timestamp:

Jun 11, 2010, 6:14:08 PM (14 years ago)

Author:

steve

Message:

Adding old files as well as workinng files from anuga_1d

Location:

anuga_work/development/pipeflow

Files:

• generic_domain.py (modified) (2 diffs)
• old_quantity.py (added)
• old_sww_comp_flux_ext.c (added)
• old_test_quantity.py (copied) (copied from anuga_work/development/pipeflow/test_quantity.py)
• old_test_sww_domain_vel.py (copied) (copied from anuga_work/development/pipeflow/test_sww_domain_vel.py)
• quantity.py (modified) (35 diffs)
• sww_comp_flux_ext.c (modified) (16 diffs)
• sww_domain.py (modified) (9 diffs)
• sww_python.py (modified) (1 diff)
• test_quantity.py (added)
• test_sww_domain.py (added)
• test_sww_domain_vel.py (added)

anuga_work/development/pipeflow/generic_domain.py

@@ -12 +12 @@
 
 
-class Generic_Domain:
+class Generic_domain:
 
     def __init__(self,
@@ -587 +587 @@
             return self.vertices[elem_id,vertex]
 
-    def get_area(self, elem_id):
+    def get_area(self, elem_id=None):
         """Return area of element id
         """
 
-        return self.areas[elem_id]
+        if elem_id is None:
+            return sum(self.areas)
+        else:
+            return self.areas[elem_id]
+
 
     def get_quantity(self, name, location='vertices', indices = None):

anuga_work/development/pipeflow/quantity.py

@@ -1 @@
-"""Class Quantity - Implements values at each 1d element
+"""
+Class Quantity - Implements values at each 1d element
 
 To create:
@@ -17 +18 @@
 import numpy
 
+
+
 class Quantity:
 
@@ -23 +26 @@
         #Initialise Quantity using optional vertex values.
         
-        from generic_domain import Generic_Domain
+        from generic_domain import Generic_domain
 
         msg = 'First argument in Quantity.__init__ '
         msg += 'must be of class Domain (or a subclass thereof)'
-        assert isinstance(domain, Generic_Domain), msg
+        assert isinstance(domain, Generic_domain), msg
 
         if vertex_values is None:
@@ -52 +55 @@
         self.centroid_values = numpy.zeros(N, numpy.float)
         self.centroid_backup_values = numpy.zeros(N, numpy.float)
-        #self.edge_values = zeros((N, 2), numpy.float)
+        #self.edge_values = numpy.zeros((N, 2), numpy.float)
         #edge values are values of the ends of each interval
       
@@ -58 +61 @@
         self.interpolate()
 
+
         
         #Allocate space for boundary values
@@ -82 +86 @@
         """
         return self.centroid_values.shape[0]
-    
+
+    def __neg__(self):
+        """Negate all values in this quantity giving meaning to the
+        expression -Q where Q is an instance of class Quantity
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(-self.vertex_values)
+        return Q
+
+
+
+    def __add__(self, other):
+        """Add to self anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
+        return result
+
+    def __radd__(self, other):
+        """Handle cases like 7+Q, where Q is an instance of class Quantity
+        """
+
+        return self + other
+
+    def __sub__(self, other):
+        return self + -other            # Invoke self.__neg__()
+
+    def __mul__(self, other):
+        """Multiply self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The product of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be products from q1 and q2
+        result.vertex_values = self.vertex_values * Q.vertex_values
+        result.centroid_values = self.centroid_values * Q.centroid_values
+
+        return result
+
+    def __rmul__(self, other):
+        """Handle cases like 3*Q, where Q is an instance of class Quantity
+        """
+
+        return self * other
+
+    def __div__(self, other):
+        """Divide self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+
+        Zero division is dealt with by adding an epsilon to the divisore
+        FIXME (Ole): Replace this with native INF once we migrate to NumPy
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The quotient of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be quotient of q1 and q2
+        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
+        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
+
+        return result
+
+    def __rdiv__(self, other):
+        """Handle cases like 3/Q, where Q is an instance of class Quantity
+        """
+
+        return self / other
+
+    def __pow__(self, other):
+        """Raise quantity to (numerical) power
+
+        As with __mul__ vertex values are processed entry by entry
+        while centroid and edge values are re-interpolated.
+
+        Example using __pow__:
+          Q = (Q1**2 + Q2**2)**0.5
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The power of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be correct
+        result.vertex_values = self.vertex_values ** other
+        result.centroid_values = self.centroid_values ** other
+
+        return result
+
+    def set_values_from_numeric(self, numeric):
+        
+        x = numpy.array([1.0,2.0])
+        y = [1.0,2.0]
+
+        if type(numeric) == type(y):
+            self.set_values_from_array(numeric)
+        elif type(numeric) == type(x):
+            self.set_values_from_array(numeric)
+        elif callable(numeric):
+            self.set_values_from_function(numeric)
+        elif isinstance(numeric, Quantity):
+            self.set_values_from_quantity(numeric)
+        else:   # see if it's coercible to a float (float, int or long, etc)
+            try:
+                numeric = float(numeric)
+            except ValueError:
+                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
+                raise Exception(msg)
+            self.set_values_from_constant(numeric)
+
+    def set_values_from_constant(self,numeric):
+
+        self.vertex_values[:,:]   = numeric
+        self.centroid_values[:,] = numeric
+
+
+    def set_values_from_array(self,numeric):
+
+        self.vertex_values[:,:]   = numeric
+        self.interpolate()
+
+
+    def set_values_from_quantity(self,quantity):
+
+        self.vertex_values[:,:]   = quantity.vertex_values
+        self.centroid_values[:,] = quantity.centroid_values
+
+    def set_values_from_function(self,function):
+
+        self.vertex_values[:,:]   = map(function, self.domain.vertices)
+        self.centroid_values[:,] = map(function, self.domain.centroids)
+
+
     def interpolate(self):
         """
@@ -96 +271 @@
             self.centroid_values[i] = (v0 + v1)/2.0
 
+
+
+
     def set_values(self, X, location='vertices'):
         """Set values for quantity
 
-        X: Compatible list, numpy array (see below), constant or function
+        X: Compatible list, Numeric array (see below), constant or function
         location: Where values are to be stored.
                   Permissible options are: vertices, centroid
@@ -105 +283 @@
 
         In case of location == 'centroid' the dimension values must
-        be a list of a numpyal array of length N, N being the number
+        be a list of a Numerical array of length N, N being the number
         of elements in the mesh. Otherwise it must be of dimension Nx2
 
@@ -119 +297 @@
         """
 
-        if location not in ['vertices', 'centroids']:
-            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
-            raise msg
+        if location not in ['vertices', 'centroids', 'unique_vertices']:
+            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
+            raise Exception(msg)
 
         if X is None:
             msg = 'Given values are None'
-            raise msg
+            raise Exception(msg)
 
         import types
@@ -134 +312 @@
           
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
-            if location == 'centroids':
-                self.centroid_values[:] = X
-            else:
-                self.vertex_values[:] = X
+            
+            self.centroid_values[:,] = float(X)
+            self.vertex_values[:,:] = float(X)
+
+        elif isinstance(X, Quantity):
+            self.set_array_values(X.vertex_values, location)
 
         else:
@@ -143 +323 @@
             self.set_array_values(X, location)
 
-        if location == 'vertices':
-            #Intialise centroid and edge values
+        if location == 'vertices' or location == 'unique_vertices':
+            #Intialise centroid
             self.interpolate()
+
+        if location == 'centroid':
+            self.extrapolate_first_order()
 
 
@@ -176 +359 @@
         """Set values for quantity
 
-        values: numpy array
+        values: Numeric array
         location: Where values are to be stored.
-                  Permissible options are: vertices, centroid, edges
+                  Permissible options are: vertices, centroid, unique_vertices
                   Default is "vertices"
 
         In case of location == 'centroid' the dimension values must
-        be a list of a numpyal array of length N, N being the number
-        of elements in the mesh. Otherwise it must be of dimension Nx2
+        be a list of a Numerical array of length N, N being the number
+        of elements in the mesh. If location == 'unique_vertices' the
+        dimension values must be a list or a Numeric array of length N+1.
+        Otherwise it must be of dimension Nx2
 
         The values will be stored in elements following their
@@ -195 +380 @@
 
 
-        values = numpy.array(values, numpy.float)
+        values = numpy.array(values).astype(numpy.float)
 
         N = self.centroid_values.shape[0]
 
-        msg = 'Number of values must match number of elements'
-        assert values.shape[0] == N, msg
 
         if location == 'centroids':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             assert len(values.shape) == 1, 'Values array must be 1d'
-            self.centroid_values = values
-        #elif location == 'edges':
-        #    assert len(values.shape) == 2, 'Values array must be 2d'
-        #    msg = 'Array must be N x 2'
-        #    self.edge_values = values
-        else:
+
+            for i in range(N):
+                self.centroid_values[i] = values[i]
+
+            self.vertex_values[:,0] = values
+            self.vertex_values[:,1] = values
+ 
+        if location == 'vertices':
+            msg = 'Number of values must match number of elements'
+            assert values.shape[0] == N, msg
             assert len(values.shape) == 2, 'Values array must be 2d'
             msg = 'Array must be N x 2'
             assert values.shape[1] == 2, msg
 
-            self.vertex_values = values
+            self.vertex_values[:,:] = values[:,:]
+
+        if location == 'unique_vertices':
+            msg = 'Number of values must match number of elements +1'
+            assert values.shape[0] == N+1, msg
+            assert len(values.shape) == 1, 'Values array must be 1d'
+
+            self.vertex_values[:,0] = values[0:-1]
+            self.vertex_values[:,1] = values[1:N+1]  
+
+
+            
 
 
@@ -220 +420 @@
         """get values for quantity
 
-        return X, Compatible list, numpy array (see below)
+        return X, Compatible list, Numeric array (see below)
         location: Where values are to be stored.
                   Permissible options are: vertices, centroid
@@ -226 +426 @@
 
         In case of location == 'centroids' the dimension values must
-        be a list of a numpyal array of length N, N being the number
+        be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3
 
@@ -236 +436 @@
 
         """
-        from numpy import take
+
 
         if location not in ['vertices', 'centroids', 'unique vertices']:
             msg = 'Invalid location: %s' %location
-            raise msg
-
-        import types
+            raise Exception(msg)
+
+        import types, numpy
         assert type(indices) in [types.ListType, types.NoneType,
                                  numpy.ArrayType],\
@@ -262 +462 @@
                 if cells is None:
                     msg = 'Unique vertex not associated with cells'
-                    raise msg
+                    raise Exception(msg)
 
                 # Go through all cells, vertex pairs
@@ -274 +474 @@
             if (indices ==  None):
                 indices = range(len(self))
-            return numpy.take(self.vertex_values,indices)
+            return take(self.vertex_values,indices)
 
 
@@ -284 +484 @@
         """Return vertex values like an OBJ format
 
-        The vertex values are returned as one sequence in the 1D numpy.float array A.
+        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.
 
@@ -308 +508 @@
 
         """
+
+
 
 
@@ -399 +601 @@
         """
 
+
+
         N = self.centroid_values.shape[0]
 
@@ -407 +611 @@
         denominator = numpy.ones(N, numpy.float)-timestep*self.semi_implicit_update
 
-        if numpy.sum(numpy.equal(denominator, 0.0)) > 0.0:
+        if sum(numpy.equal(denominator, 0.0)) > 0.0:
             msg = 'Zero division in semi implicit update. Call Stephen :-)'
-            raise msg
+            raise Exception(msg)
         else:
             #Update conserved_quantities from semi implicit updates
@@ -420 +624 @@
         """
 
-        #print 'compute_gradient'
+
+
 
         N = self.centroid_values.shape[0]
@@ -510 +715 @@
 
         #print 'compute_minmod_gradients'
+        from numpy import sign
+
         
         def xmin(a,b):
@@ -516 +723 @@
         def xmic(t,a,b):
             return xmin(t*xmin(a,b), 0.50*(a+b) )
+
+
 
         N = self.centroid_values.shape[0]
@@ -702 +911 @@
         """
 
+
         N = self.domain.number_of_elements
         beta = self.domain.beta
@@ -729 +939 @@
 
         #Deltas between vertex and centroid values
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] = qv[:,i] - qc
@@ -753 +963 @@
 
         from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
+
         limiter = self.domain.limiter
         #print limiter
         
         #print 'limit_range'
-        N  = self.domain.number_of_elements
+        N = self.domain.number_of_elements
         qc = self.centroid_values
         qv = self.vertex_values
-        C  = self.domain.centroids
-        X  = self.domain.vertices
+        C = self.domain.centroids
+        X = self.domain.vertices
         beta_p = numpy.zeros(N,numpy.float)
         beta_m = numpy.zeros(N,numpy.float)
@@ -777 +988 @@
                 beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
                 
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
@@ -817 +1028 @@
 
         import sys
+
         from util import minmod, minmod_kurganov, maxmod, vanleer
 
@@ -841 +1053 @@
                 # Check denominator not zero
                 if (qc[k+1]-qc[k]) == 0.0:
-                    beta_p[k] = numpy.float(sys.maxint)
-                    beta_m[k] = numpy.float(sys.maxint)
+                    beta_p[k] = float(sys.maxint)
+                    beta_m[k] = float(sys.maxint)
                 else:
                     #STEVE LIMIT
@@ -849 +1061 @@
 
         #Deltas between vertex and centroid values
-        dq = zeros(qv.shape, numpy.float)
+        dq = numpy.zeros(qv.shape, numpy.float)
         for i in range(2):
             dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
@@ -885 +1097 @@
         #backup_centroid_values(self)
 
-        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
+        self.centroid_backup_values[:,] = (self.centroid_values).astype('f')
 
     def saxpy_centroid_values(self,a,b):
         # Call correct module function
         # (either from this module or C-extension)
-        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
+        self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
         
 
@@ -914 +1126 @@
 if __name__ == "__main__":
     #from domain import Domain
-    from shallow_water_domain import Domain     
-    from numpy import arange
+    from sww_domain import Domain
     
     points1 = [0.0, 1.0, 2.0, 3.0]
@@ -968 +1179 @@
     raw_input('press return')
 
-    points2 = arange(10)
+    points2 = numpy.arange(10)
     D2 = Domain(points2)
 

anuga_work/development/pipeflow/sww_comp_flux_ext.c

@@ -10 +10 @@
 
 
+/* double max(double a, double b) { */
+/*      double z; */
+/*      z=(a>b)?a:b; */
+/*      return z;} */
+        
+/* double min(double a, double b) { */
+/*      double z; */
+/*      z=(a<b)?a:b; */
+/*      return z;} */
+
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+                      double *h, 
+                      double epsilon, 
+                      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+     u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);    
+  }
+  
+
+  // Adjust momentum to be consistent with speed
+  *uh = u * *h;
+  
+  return u;
+}
+
+
+
+
 //Innermost flux function (using w=z+h)
 int _flux_function(double *q_left, double *q_right,
+        double z_left, double z_right,
                    double normals, double g, double epsilon, double h0,
-                   double *edgeflux, double *max_speed) {
+                double *edgeflux, double *max_speed) {
                 
                 int i;
-                double flux_left[2], flux_right[2];
-                double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
-                double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
-                double z, s_max, s_min, denom;
+                double ql[2], qr[2], flux_left[2], flux_right[2];
+                double z, w_left, h_left, uh_left, soundspeed_left, u_left;
+                double w_right, h_right, uh_right, soundspeed_right, u_right;
+                double s_max, s_min, denom;
                 
                 //printf("h0 = %f \n",h0);
-                w_left  = q_left[0];
-                uh_left = q_left[1]*normals;
-                h_left  = q_left[2];
-                z_left  = q_left[3];
-                u_left  = q_left[4]*normals;
-        
-                w_right  = q_right[0];
-                uh_right = q_right[1]*normals;
-                h_right  = q_right[2];
-                z_right  = q_right[3];
-                u_right  = q_right[4]*normals;          
+                ql[0] = q_left[0];
+                ql[1] = q_left[1];
+                ql[1] = ql[1]*normals;
+        
+                qr[0] = q_right[0];
+                qr[1] = q_right[1];
+                qr[1] = qr[1]*normals;
           
                 z = (z_left+z_right)/2.0;                  
                                    
+                //w_left = ql[0];
+                //h_left = w_left-z;
+                //uh_left = ql[1];
+                
+
+
+                // Compute speeds in x-direction
+                w_left = ql[0];          
+                h_left = w_left-z;
+                uh_left = ql[1];
+
+                u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
+
+                w_right = qr[0];
+                h_right = w_right-z;
+                uh_right = qr[1];
+
+                u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
+  
                 soundspeed_left = sqrt(g*h_left);
                 soundspeed_right = sqrt(g*h_right);
@@ -48 +102 @@
                 // Flux formulas
                 flux_left[0] = u_left*h_left;
-                flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
+                flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
 
                 flux_right[0] = u_right*h_right;
-                flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
+                flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
 
                 // Flux computation
@@ -60 +114 @@
                 } else {
                         edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
-                        edgeflux[0] += s_max*s_min*(q_right[0]-q_left[0]);
+                        edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
                         edgeflux[0] /= denom;
                         edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
-                        edgeflux[1] += s_max*s_min*(q_right[1]-q_left[1]);
+                        edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
                         edgeflux[1] /= denom;
                         edgeflux[1] *= normals;
     
-                        // Maximal wavespeed
-                        *max_speed = max(fabs(s_max), fabs(s_min));
+                // Maximal wavespeed
+        *max_speed = max(fabs(s_max), fabs(s_min));
                 }  
-                return 0;               
-}
+    return 0;           
+        }
                 
                 
@@ -90 +144 @@
                            double* xmom_edge_values,
                            double* bed_edge_values,
-                           double* height_edge_values,
-                           double* velocity_edge_values,
                            double* stage_boundary_values,
                            double* xmom_boundary_values,
-                           double* bed_boundary_values,
-                           double* height_boundary_values,
-                           double* velocity_boundary_values,
                            double* stage_explicit_update,
                            double* xmom_explicit_update,
@@ -102 +151 @@
                            double* max_speed_array) {
                 
-                double flux[2], ql[5], qr[5], edgeflux[2];
-                double max_speed, normal;
+                double flux[2], ql[2], qr[2], edgeflux[2];
+                double zl, zr, max_speed, normal;
                 int k, i, ki, n, m, nm=0;
                 
@@ -116 +165 @@
                                 ql[0] = stage_edge_values[ki];
                                 ql[1] = xmom_edge_values[ki];
-                                ql[2] = bed_edge_values[ki];
-                                ql[3] = height_edge_values[ki];
-                                ql[4] = velocity_edge_values[ki];
+                                zl = bed_edge_values[ki];
                                 
                                 n = neighbours[ki];
@@ -125 +172 @@
                                         qr[0] = stage_boundary_values[m];
                                         qr[1] = xmom_boundary_values[m];
-                                        qr[2] = bed_edge_values[m];
-                                        qr[3] = height_edge_values[m];
-                                        qr[4] = velocity_edge_values[m];
+                                        zr = zl;
                                 } else {
                                         m = neighbour_vertices[ki];
@@ -133 +178 @@
                                         qr[0] = stage_edge_values[nm];
                                         qr[1] = xmom_edge_values[nm];
-                                        qr[2] = bed_edge_values[nm];
-                                        qr[3] = height_edge_values[nm];
-                                        qr[4] = velocity_edge_values[nm];
+                                        zr = bed_edge_values[nm];                               
                                 }
                                 
                                 normal = normals[ki];
-                                _flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
+                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
                                 flux[0] -= edgeflux[0];
                                 flux[1] -= edgeflux[1];
@@ -152 +195 @@
                                     }
                                 }
-                        } // End edge i (and neighbour n)
+            } // End edge i (and neighbour n)
                         flux[0] /= areas[k];
                         stage_explicit_update[k] = flux[0];
@@ -163 +206 @@
                 return timestep;        
         }
+        
+        
+        
+        
+        
+
+
 
 //=========================================================================
@@ -168 +218 @@
 //=========================================================================
 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-
-    PyObject
+  
+    PyArrayObject 
+        *neighbours, 
+        *neighbour_vertices,
+        *normals, 
+        *areas,
+        *stage_edge_values,
+        *xmom_edge_values,
+        *bed_edge_values,
+        *stage_boundary_values,
+        *xmom_boundary_values,
+        *stage_explicit_update,
+        *xmom_explicit_update,
+        *max_speed_array;
+    
+  double timestep, epsilon, g, h0, cfl;
+  int number_of_elements;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddddOOOOOOOOOOOiO",
+                        &cfl,
+                        &timestep,
+                        &epsilon,
+                        &g,
+                        &h0,
+                        &neighbours,
+                        &neighbour_vertices,
+                        &normals,
+                        &areas,
+                        &stage_edge_values,
+                        &xmom_edge_values,
+                        &bed_edge_values,
+                        &stage_boundary_values,
+                        &xmom_boundary_values,
+                        &stage_explicit_update,
+                        &xmom_explicit_update,
+                        &number_of_elements,
+                        &max_speed_array)) {
+    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
+    return NULL;
+  }
+
+ 
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_ext(cfl,
+                                 timestep,
+                                 epsilon,
+                                 g,
+                                 h0,
+                                 (long*) neighbours -> data,
+                                 (long*) neighbour_vertices -> data,
+                                 (double*) normals -> data,
+                                 (double*) areas -> data,
+                                 (double*) stage_edge_values -> data,
+                                 (double*) xmom_edge_values -> data,
+                                 (double*) bed_edge_values -> data,
+                                 (double*) stage_boundary_values -> data,
+                                 (double*) xmom_boundary_values -> data,
+                                 (double*) stage_explicit_update -> data,
+                                 (double*) xmom_explicit_update -> data,
+                                 number_of_elements,
+                                 (double*) max_speed_array -> data);
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+//-----------------------
+PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
+  
+   PyObject 
         *domain,
-        *stage,
-        *xmom,
-        *bed,
-        *height,
-        *velocity;
-
-    PyArrayObject
-        *neighbours,
+        *stage, 
+        *xmom, 
+        *bed;
+
+    PyArrayObject 
+        *neighbours, 
         *neighbour_vertices,
-        *normals,
+        *normals, 
         *areas,
         *stage_vertex_values,
         *xmom_vertex_values,
         *bed_vertex_values,
-        *height_vertex_values,
-        *velocity_vertex_values,
         *stage_boundary_values,
         *xmom_boundary_values,
-        *bed_boundary_values,
-        *height_boundary_values,
-        *velocity_boundary_values,
         *stage_explicit_update,
         *xmom_explicit_update,
         *max_speed_array;
-
+    
   double timestep, epsilon, g, h0, cfl;
   int number_of_elements;
 
-
+    
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dOOOOOO",
+  if (!PyArg_ParseTuple(args, "dOOOO",
                         &timestep,
                         &domain,
                         &stage,
                         &xmom,
-                        &bed,
-                        &height,
-                        &velocity)) {
+                        &bed)) {
       PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
       return NULL;
@@ -217 +329 @@
     g                 = get_python_double(domain,"g");
     h0                = get_python_double(domain,"h0");
-    cfl               = get_python_double(domain,"cfl");
-
-
+    cfl               = get_python_double(domain,"CFL");
+ 
+    
     neighbours        = get_consecutive_array(domain, "neighbours");
-    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
     normals           = get_consecutive_array(domain, "normals");
-    areas             = get_consecutive_array(domain, "areas");
+    areas             = get_consecutive_array(domain, "areas");    
     max_speed_array   = get_consecutive_array(domain, "max_speed_array");
-
-    stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");
-    xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");
-    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
-    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
-    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
-
-    stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");
-    xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");
-    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
-    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
-    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
-
-
-    stage_explicit_update = get_consecutive_array(stage, "explicit_update");
-    xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");
+    
+    stage_vertex_values = get_consecutive_array(stage, "vertex_values");    
+    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");    
+    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+
+
 
     number_of_elements = stage_vertex_values -> dimensions[0];
 
-    // Call underlying flux computation routine and update
-    // the explicit update arrays
+
+ 
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
     timestep = _compute_fluxes_ext(
         cfl,
@@ -259 +370 @@
         (double*) xmom_vertex_values -> data,
         (double*) bed_vertex_values -> data,
-        (double*) height_vertex_values -> data,
-        (double*) velocity_vertex_values -> data,
         (double*) stage_boundary_values -> data,
         (double*) xmom_boundary_values -> data,
-        (double*) bed_boundary_values -> data,
-        (double*) height_boundary_values -> data,
-        (double*) velocity_boundary_values -> data,
         (double*) stage_explicit_update -> data,
         (double*) xmom_explicit_update -> data,
@@ -279 +385 @@
   Py_DECREF(xmom_vertex_values);
   Py_DECREF(bed_vertex_values);
-  Py_DECREF(height_vertex_values);
-  Py_DECREF(velocity_vertex_values);
   Py_DECREF(stage_boundary_values);
   Py_DECREF(xmom_boundary_values);
-  Py_DECREF(bed_boundary_values);
-  Py_DECREF(height_boundary_values);
-  Py_DECREF(velocity_boundary_values);
   Py_DECREF(stage_explicit_update);
   Py_DECREF(xmom_explicit_update);
@@ -291 +392 @@
 
 
+
+
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
 
+ 
 
 
@@ -303 +407 @@
 static struct PyMethodDef MethodTable[] = {
   {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
   {NULL, NULL}
 };

anuga_work/development/pipeflow/sww_domain.py

@@ -50 +50 @@
 
 #Shallow water domain
-class Domain(Generic_Domain):
+class Domain(Generic_domain):
 
     def __init__(self, coordinates, boundary = None, tagged_elements = None):
@@ -57 +57 @@
         evolved_quantities   = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
         other_quantities     = ['friction']
-        Generic_Domain.__init__(self,
+        Generic_domain.__init__(self,
                                 coordinates          = coordinates,
                                 boundary             = boundary,
@@ -143 +143 @@
 
 
-        timestep = numpy.float(sys.maxint)
+        timestep = float(sys.maxint)
 
         stage    = self.quantities['stage']
@@ -152 +152 @@
 
 
-        from sww_comp_flux_ext import compute_fluxes_ext
-
-        self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
+        from sww_comp_flux_ext import compute_fluxes_ext_short
+
+        #self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
+
+    
+        self.flux_timestep = compute_fluxes_ext_short(timestep,self,stage,xmom,bed)
+
 
 
@@ -168 +172 @@
 
         #Call basic machinery from parent class
-        for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
+        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
                                        skip_initial_step):
 
@@ -266 +270 @@
 ##      else:
 
-    u_C[:]  = uh_C/(h_C + h0/h_C)
+    u_C[:,]  = uh_C/(h_C + h0/h_C)
         
     for name in [ 'velocity', 'stage' ]:
@@ -286 +290 @@
     xmom_V  = domain.quantities['xmomentum'].vertex_values
 
-    h_V[:]    = stage_V - bed_V
+    h_V[:,:]    = stage_V - bed_V
     for i in range(len(h_C)):
         for j in range(2):
@@ -297 +301 @@
                 stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
                 
-    xmom_V[:] = u_V * h_V
+    xmom_V[:,:] = u_V * h_V
     
     return
@@ -314 +318 @@
 
     #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
+    wc    = domain.quantities['stage'].centroid_values
+    zc    = domain.quantities['elevation'].centroid_values
     xmomc = domain.quantities['xmomentum'].centroid_values
     hc = wc - zc  #Water depths at centroids

anuga_work/development/pipeflow/sww_python.py

@@ -58 +58 @@
     for k in range(N):
 
-        flux[:] = 0.  #Reset work numpy.array
+        flux[:,] = 0.  #Reset work numpy.array
         #for i in range(3):
         for i in range(2):