Changeset 8154

Timestamp:

Mar 17, 2011, 5:51:42 PM (14 years ago)

Author:

steve

Message:

Adding in kinematic viscosity. Added in some procedures to change boundary_values of
quantities

Location:

trunk/anuga_core/source/anuga

Files:

• __init__.py (modified) (2 diffs)
• abstract_2d_finite_volumes/generic_domain.py (modified) (3 diffs)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_quantity.py (modified) (5 diffs)
• operators/kinematic_viscosity.py (modified) (13 diffs)
• operators/test_kinematic_viscosity.py (modified) (12 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (3 diffs)
• utilities/cg_solve.py (modified) (4 diffs)
• utilities/sparse.py (modified) (1 diff)

trunk/anuga_core/source/anuga/__init__.py

@@ -29 +29 @@
 
 from anuga.shallow_water.shallow_water_domain import Domain
+
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
 
 from anuga.abstract_2d_finite_volumes.util import file_function, \
@@ -132 +134 @@
 
 #---------------------------
-# Structures
+# Structure Operators
 #---------------------------
 from anuga.structures.structure_operator import Structure_operator

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -100 +100 @@
         self.vertex_coordinates = self.mesh.vertex_coordinates
         self.boundary = self.mesh.boundary
+        self.boundary_enumeration = self.mesh.boundary_enumeration
         self.neighbours = self.mesh.neighbours
         self.surrogate_neighbours = self.mesh.surrogate_neighbours
@@ -708 +709 @@
                 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.neighbours[vol_id, edge_id] = \
+                                        #-len(self.boundary_objects)
                 else:
                     pass
@@ -721 +722 @@
                 raise Exception(msg)
 
+    ##
+    # @brief Set quantities based on a regional tag.
+    # @param args
+    # @param kwargs
     def set_region(self, *args, **kwargs):
         """Set quantities based on a regional tag.

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -178 +178 @@
         if verbose: log.critical('Mesh: Building boundary dictionary')
         self.build_boundary_dictionary(boundary)
+
+        #Update boundary_enumeration
+        self.build_boundary_neighbours()
 
         #Build tagged element  dictionary mapping (tag) to array of elements
@@ -441 +444 @@
 
             self.neighbours[id, edge] = index
+
+            self.boundary_enumeration[id,edge] = index
+
             index -= 1
+
+    def build_boundary_neighbours(self):
+        """Traverse boundary and
+        enumerate neighbour indices from -1 and
+        counting down.
+
+        Precondition:
+            self.boundary is defined.
+        Post condition:
+            neighbours array has unique negative indices for boundary
+            boundary_segments array imposes an ordering on segments
+            (not otherwise available from the dictionary)
+
+        """
+
+        if self.boundary is None:
+            msg = 'Boundary dictionary must be defined before '
+            msg += 'building boundary structure'
+            raise Exception(msg)
+
+        self.boundary_enumeration = {}
+
+        X = self.boundary.keys()
+        X.sort()
+
+        index = -1
+        for id, edge in X:
+            self.neighbours[id, edge] = index
+
+            self.boundary_enumeration[id,edge] = -index -1 
+
+            index -= 1
+
 
 

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -234 +234 @@
         self.beta = beta
 
+    ##
+    # @brief Get the current beta value.
+    # @return The current beta value.
     def get_beta(self):
         """Get default beta value for limiting"""
@@ -239 +242 @@
         return self.beta
 
+
+
+    ##
+    # @brief Set boundary values using a function or array or scalar
+    # @param numeric: function or array or scalar
+    def set_boundary_values(self, numeric = 0.0):
+        """Set boundary values """
+
+        if isinstance(numeric, (list, num.ndarray)):
+            self._set_boundary_values_from_array(numeric)
+        elif callable(numeric):
+            self._set_boundary_values_from_function(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_boundary_values_from_constant(numeric)
+
+
+    ##
+    # @brief Set boundary values using a function
+    # @param numeric: function
+    def _set_boundary_values_from_function(self, function):
+        """Set boundary values from function of x,y """
+
+        for (vol_id, edge_id) , j in self.domain.boundary_enumeration.items():
+            [x,y] = self.domain.get_edge_midpoint_coordinates(vol_id)[edge_id]
+            self.boundary_values[j] = function(x,y)
+
+    ##
+    # @brief Set boundary values using a scalar
+    # @param numeric: scalar
+    def _set_boundary_values_from_constant(self, scalar):
+        """Set boundary values from scalar """
+
+        for (vol_id, edge_id) , j in self.domain.boundary_enumeration.items():
+            [x,y] = self.domain.get_edge_midpoint_coordinates(vol_id)[edge_id]
+            self.boundary_values[j] = scalar
+
+    ##
+    # @brief Set boundary values using a scalar
+    # @param numeric: scalar
+    def _set_boundary_values_from_array(self, array):
+        """Set boundary values from array """
+
+        assert len(array)  == len(self.domain.boundary_enumeration)
+        
+        for (vol_id, edge_id) , j in self.domain.boundary_enumeration.items():
+            [x,y] = self.domain.get_edge_midpoint_coordinates(vol_id)[edge_id]
+            self.boundary_values[j] = array[j]
+
+    ##
+    # @brief Compute interpolated values at edges and centroid.
+    # @note vertex_values must be set before calling this.
     def interpolate(self):
         """Compute interpolated values at edges and centroid
@@ -583 +643 @@
                 assert values.shape[0] == N, msg
 
-                self.centroid_values = values
+                self.centroid_values[:] = values
             else:
                 msg = 'Number of values must match number of indices'
@@ -610 +670 @@
 
                 if indices is None:
-                    self.vertex_values = values
+                    self.vertex_values[:] = values
                 else:
                     for element_index, value in map(None, indices, values):

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1233 +1233 @@
 
         neighbours = mesh.get_triangle_neighbours(0)
-        assert num.allclose(neighbours, [-1,1,-1])
+        assert num.allclose(neighbours, [-1,1,-2])
         neighbours = mesh.get_triangle_neighbours(-10)
         assert neighbours == []

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -97 +97 @@
                                                      [0.,0.,0.], [0.,0.,0.]])
 
+    def test_set_boundary_values(self):
+
+        quantity = Quantity(self.mesh1)
+
+        quantity.set_boundary_values()
+
+        assert  num.allclose(quantity.boundary_values, [0.0, 0.0, 0.0])
+
+
+    def test_set_boundary_values_with_function(self):
+
+        quantity = Quantity(self.mesh1)
+        #assert num.allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+        def simple(x,y):
+            return x+3*y
+
+        quantity.set_boundary_values(simple)
+
+        assert  num.allclose(quantity.boundary_values, [1.0, 4.0, 3.0])
+
+    def test_set_boundary_values_with_constant(self):
+
+        quantity = Quantity(self.mesh1)
+        #assert num.allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+        quantity.set_boundary_values(10.0)
+
+        assert  num.allclose(quantity.boundary_values, [10.0, 10.0, 10.0])
+
+
+    def test_set_boundary_values_with_array(self):
+
+        quantity = Quantity(self.mesh1)
+        #assert num.allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+
+        quantity.set_boundary_values([10.0, 4.0, 5.0])
+
+        assert  num.allclose(quantity.boundary_values, [10.0, 4.0, 5.0])
+
+    def test_set_boundary_values_with_wrong_sized_array(self):
+
+        quantity = Quantity(self.mesh1)
+        #assert num.allclose(quantity.vertex_values, [[0.,0.,0.]])
+
+        try:
+            quantity.set_boundary_values([10.0, 4.0, 5.0, 8.0])
+        except:
+            pass
+        else:
+            msg = 'Should have caught this'
+            raise Exception(msg)
 
     def test_interpolation(self):
@@ -247 +300 @@
         quantity = Quantity(self.mesh4)
 
+        # get referece to data arrays
+        centroid_values = quantity.centroid_values
+        vertex_values = quantity.vertex_values
 
         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
@@ -252 +308 @@
         assert num.allclose(quantity.vertex_values,
                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+
+        assert id(vertex_values) == id(quantity.vertex_values)
         assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
         assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
@@ -281 +339 @@
 
 
+
         try:
             quantity.set_values([[1,2,3], [0,0,9]])
-        except AssertionError:
+        except ValueError:
             pass
         except:
-            raise Exception('should have raised Assertionerror')
+            raise Exception('should have raised ValueeError')
 
 
@@ -2368 +2427 @@
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')    
+    suite = unittest.makeSuite(Test_Quantity, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py

@@ -9 +9 @@
 import anuga.utilities.log as log
 
+
+
 class Kinematic_Viscosity:
     """
@@ -35 +37 @@
     """
 
-    def __init__(self, domain, diffusivity = None, use_triangle_areas=True, verbose=False):
+    def __init__(self, domain, use_triangle_areas=True, verbose=False):
         if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
         #Expose the domain attributes
@@ -44 +46 @@
         self.boundary_enumeration = domain.boundary_enumeration
         
-        # Setup diffusivity quantity
-        if diffusivity is None:
-            diffusivity = Quantity(domain)
-            diffusivity.set_values(1.0)
-            diffusivity.set_boundary_values(1.0)
-        #check that diffusivity is a quantity associated with domain
-        assert diffusivity.domain == domain
-
-        self.diffusivity = diffusivity
-        #self.diffusivity_bdry_data = self.diffusivity.boundary_values
-        #self.diffusivity_cell_data = self.diffusivity.centroid_values
-
+        # Dummy diffusivity quantity
+        diffusivity = Quantity(domain)
+        diffusivity.set_values(1.0)
+        diffusivity.set_boundary_values(1.0)
 
         self.n = len(self.domain)
-        self.dt = 1e-6 #Need to set to domain.timestep
+        self.dt = 1.0 #Need to set to domain.timestep
         self.boundary_len = len(domain.boundary)
         self.tot_len = self.n + self.boundary_len
@@ -95 +89 @@
 
         # Build matrix self.elliptic_matrix [A B]
-        self.build_elliptic_matrix()
-
-        # Build self.boundary_term
-        #self.build_elliptic_boundary_term()
-
-        self.parabolic_solve = False #Are we doing a parabolic solve at the moment?
+        self.build_elliptic_matrix(diffusivity)
+
+        self.boundary_term = num.zeros((self.n, ), num.float)
+
+        self.parabolic = False #Are we doing a parabolic solve at the moment?
 
         if verbose: log.critical('Kinematic Viscosity: Initialisation Done')
@@ -109 +102 @@
         
 
-    def set_qty_considered(self, qty):
-        # FIXME SR: Probably should just be set by quantity to which operation is applied
-        if qty == 1 or qty == 'u':
-            self.qty_considered = 1
-        elif qty == 2 or qty == 'v':
-            self.qty_considered = 2
-        else: #Raise an exception
-            msg = "Incorrect input qty"
-            assert 0 == 1, msg
-
-    def build_elliptic_matrix(self):
+    def set_parabolic_solve(self,flag):
+
+        self.parabolic = flag
+
+
+    def build_elliptic_matrix(self, a):
         """
         Builds matrix representing
 
-        div ( diffusivity grad )
+        div ( a grad )
 
         which has the form [ A B ]
@@ -131 +119 @@
         # are setup via this call
         kinematic_viscosity_ext.build_elliptic_matrix(self, \
-                self.diffusivity.centroid_values, \
-                self.diffusivity.boundary_values)
+                a.centroid_values, \
+                a.boundary_values)
 
         self.elliptic_matrix = Sparse_CSR(None, \
@@ -138 +126 @@
                 self.n, self.tot_len)
 
-        #print 'elliptic_matrix'
-        #print self.elliptic_matrix
-
-#        #Set up the scaling matrix
-#        data = h
-#        num.putmask(data, data != 0, 1 / data) #take the reciprocal of each entry unless it is zero
-#        self.stage_heights_scaling = \
-#            Sparse_CSR(None, self.diffusivity_data, num.arange(self.n), num.arange(self.n +1), self.n, self.n)
-
-    def update_elliptic_matrix(self):
+
+    def update_elliptic_matrix(self, a=None):
         """
         Updates the data values of matrix representing
 
-        div ( diffusivity grad )
-
-        (as when diffusivitiy has changed)
+        div ( a grad )
+
+        If a == None then we set a = quantity which is set to 1
         """
 
@@ -159 +139 @@
         # through to the Sparse_CSR matrix.
 
+        if a == None:
+            a = Quantity(self.domain)
+            a.set_values(1.0)
+            a.set_boundary_values(1.0)
+            
         kinematic_viscosity_ext.update_elliptic_matrix(self, \
-                self.diffusivity.centroid_values, \
-                self.diffusivity.boundary_values)
-        
-
-    def build_elliptic_boundary_term(self,quantity):
-        """
-        Operator has form [A B] and U = [ u ; b]
-
-        This procedure calculates B b which can be calculated as
+                a.centroid_values, \
+                a.boundary_values)
+        
+
+
+
+
+    def update_elliptic_boundary_term(self, boundary):
+
+
+        if isinstance(boundary, Quantity):
+
+            self._update_elliptic_boundary_term(boundary.boundary_values)
+
+        elif isinstance(boundary, num.ndarray):
+
+            self._update_elliptic_boundary_term(boundary.boundary_values)
+
+        else:
+
+            raise  TypeError('expecting quantity or numpy array')
+
+
+    def _update_elliptic_boundary_term(self, b):
+        """
+        Operator has form [A B] and u = [ u_1 ; b]
+
+        u_1 associated with centroid values of u
+        u_2 associated with boundary_values of u
+
+        This procedure calculates B u_2 which can be calculated as
 
         [A B] [ 0 ; b]
+
+        Assumes that update_elliptic_matrix has just been run.
         """
 
@@ -176 +185 @@
         tot_len = self.tot_len
 
-        self.boundary_term = num.zeros((n, ), num.float)
-        
         X = num.zeros((tot_len,), num.float)
 
-        X[n:] = quantity.boundary_values
+        X[n:] = b
         self.boundary_term[:] = self.elliptic_matrix * X
 
@@ -188 +195 @@
 
 
-    def elliptic_multiply(self, quantity_in, quantity_out=None, include_boundary=True):
+
+    def elliptic_multiply(self, input, output=None):
+
+
+        if isinstance(input, Quantity):
+
+            assert isinstance(output, Quantity) or output is None
+
+            output = self._elliptic_multiply_quantity(input, output)
+
+        elif isinstance(input, num.ndarray):
+
+            assert isinstance(output, num.ndarray) or output is None
+
+            output = self._elliptic_multiply_array(input, output)
+
+        else:
+
+            raise TypeError('expecting quantity or numpy array')
+        
+        return output
+
+
+    def _elliptic_multiply_quantity(self, quantity_in, quantity_out):
+        """
+        Assumes that update_elliptic_matrix has been run
+        """
+
+        if quantity_out is None:
+            quantity_out = Quantity(self.domain)
+
+        array_in = quantity_in.centroid_values
+        array_out = quantity_out.centroid_values
+
+        X = self._elliptic_multiply_array(array_in, array_out)
+
+        quantity_out.set_values(X, location = 'centroids')
+        
+        return quantity_out
+
+    def _elliptic_multiply_array(self, array_in, array_out):
+        """
+        calculates [A B] [array_in ; 0]
+        """
 
         n = self.n
@@ -194 +244 @@
 
         V = num.zeros((tot_len,), num.float)
-        X = num.zeros((n,), num.float)
+
+        assert len(array_in) == n
+
+        if array_out is None:
+            array_out = num.zeros_like(array_in)
+
+        V[0:n] = array_in
+        V[n:] = 0.0
+
+
+        if self.apply_triangle_areas:
+            V[0:n] = self.triangle_areas * V[0:n]
+
+
+        array_out[:] = self.elliptic_matrix * V
+
+
+        return array_out
+
+
+
+
+
+    def parabolic_multiply(self, input, output=None):
+
+
+        if isinstance(input, Quantity):
+
+            assert isinstance(output, Quantity) or output is None
+
+            output = self._parabolic_multiply_quantity(input, output)
+
+        elif isinstance(input, num.ndarray):
+
+            assert isinstance(output, num.ndarray) or output is None
+
+            output = self._parabolic_multiply_array(input, output)
+
+        else:
+
+            raise TypeError('expecting quantity or numpy array')
+
+        return output
+
+
+    def _parabolic_multiply_quantity(self, quantity_in, quantity_out):
+        """
+        Assumes that update_elliptic_matrix has been run
+        """
 
         if quantity_out is None:
             quantity_out = Quantity(self.domain)
 
-        V[0:n] = quantity_in.centroid_values
+        array_in = quantity_in.centroid_values
+        array_out = quantity_out.centroid_values
+
+        X = self._parabolic_multiply_array(array_in, array_out)
+
+        quantity_out.set_values(X, location = 'centroids')
+
+        return quantity_out
+
+    def _parabolic_multiply_array(self, array_in, array_out):
+        """
+        calculates ( [ I 0 ; 0  0] + dt [A B] ) [array_in ; 0]
+        """
+
+        n = self.n
+        tot_len = self.tot_len
+
+        V = num.zeros((tot_len,), num.float)
+
+        assert len(array_in) == n
+
+        if array_out is None:
+            array_out = num.zeros_like(array_in)
+
+        V[0:n] = array_in
         V[n:] = 0.0
-
-
-        # FIXME SR: These sparse matrix vector multiplications
-        # should be done in such a way to reuse memory, ie
-        # should have a procedure
-        # matrix.multiply(vector_in, vector_out)
 
 
@@ -213 +329 @@
 
 
-        X[:] = self.elliptic_matrix * V
-
-        if include_boundary:
-            self.build_elliptic_boundary_term(quantity_in)
-            X[:] += self.boundary_term
-
-
-        quantity_out.set_values(X, location = 'centroids')
-        return quantity_out
-
-    #parabolic_multiply(V) = identity - dt*elliptic_multiply(V)
-    #We do not need include boundary, as we will only use this
-    #method for solving (include_boundary = False)
-    #Here V is already either (u) or (v), not (uh) or (vh)
-    def parabolic_multiply(self, V):
-        msg = "(KV_Operator.parabolic_multiply) V vector has incorrect dimensions"
-        assert V.shape == (self.n, ) or V.shape == (self.n, 1), msg
-
-        D = V #The return value
-        X = V #a temporary array
-
-        #Apply triangle areas
-        if self.apply_triangle_areas:
-            X = self.triangle_areas * X
-
-        #Multiply out
-        X = num.append(X, num.zeros((self.boundary_len, )), axis=0)
-        D = D - self.dt * (self.elliptic_matrix * X)
-
-        return num.array(D).reshape(self.n, )
-
-    def __mul__(self, other):
-        try:
-            B = num.array(other)
-        except:
-            msg = "Trying to multiply the Kinematic Viscosity Operator against a non-numeric type"
-            raise msg
-
-        if len(B.shape) == 0:
-            #Scalar
-            R = B * self
-        elif len(B.shape) == 1 or (len(B.shape) == 2 and B.shape[1] == 1):
-            #Vector
-            if self.parabolic_solve:
-                R = self.parabolic_multiply(other)
-            else:
-                #include_boundary=False is this is *only* used for cg_solve()
-                R = self.elliptic_multiply(other, include_boundary=False)
+        array_out[:] = array_in - self.dt * (self.elliptic_matrix * V)
+
+        return array_out
+
+
+
+
+    def __mul__(self, vector):
+        
+        #Vector
+        if self.parabolic:
+            R = self.parabolic_multiply(vector)
         else:
-            raise ValueError, 'Dimension too high: d=%d' % len(B.shape)
+            #include_boundary=False is this is *only* used for cg_solve()
+            R = self.elliptic_multiply(vector)
+
         return R
     
@@ -271 +353 @@
         except:
             msg = 'Sparse matrix can only "right-multiply" onto a scalar'
-            raise TypeError, msg
+            raise TypeError(msg)
         else:
             new = self.elliptic_matrix * new
         return new
 
-    def cg_solve(self, B, qty_to_consider=None):
-        if len(B.shape) == 1:
-            return self.cg_solve_vector(B, qty_to_consider)
-        elif len(B.shape) == 2:
-            return self.cg_solve_matrix(B)
-        else:
-            raise ValueError, 'Dimension too high: d=%d' % len(B.shape)
-
-    def cg_solve_matrix(self, B):
-        assert B.shape[1] < 3, "Input matrix has too many columns (max 2)"
-
-        X = num.zeros(B.shape, num.float)
-        for i in range(B.shape[1]):
-            X[:, i] = self.cg_solve_vector(B[:, i], i + 1) #assuming B columns are (uh) and (vh)
-        return X
     
-    def cg_solve_vector(self, b, qty_to_consider=None):
-        if not qty_to_consider == None:
-            self.set_qty_considered(qty_to_consider)
-            #Call the ANUGA conjugate gradient utility
-        x = conjugate_gradient(self, b - self.boundary_vector[:, self.qty_considered-1])
-        return x
-
-    #Solve the parabolic equation to find u^(n+1), where u = u^n
-    #Here B = [uh vh] where uh,vh = n*1 column vectors
-    def parabolic_solver(self, B):
-        assert B.shape == (self.n, 2), "Input matrix has incorrect dimensions"
-
-        next_B = num.zeros((self.n, 2), num.float)
-        #The equation is self.parabolic_multiply*u^(n+1) = u^n + (dt)*self.boundary_vector
-        self.parabolic_solve = True
-        #Change (uh) and (vh) to (u) and (v)
-        b = self.stage_heights_scaling * B
-
-        next_B = conjugate_gradient(self, b + (self.dt * self.boundary_vector), iprint=1)
-
-        #Return (u) and (v) back to (uh) and (vh)
-        next_B = self.stage_heights_scaling * next_B
-
-        self.parabolic_solve = False
-
-        return next_B
+    def elliptic_solve(self, u_in, b, a = None, u_out = None, \
+                       imax=10000, tol=1.0e-8, atol=1.0e-8, iprint=None):
+        """ Solving div ( a grad u ) = b
+        u | boundary = g
+
+        u_in, u_out, f anf g are Quantity objects
+
+        Dirichlet BC g encoded into u_in boundary_values
+
+        Initial guess for iterative scheme is given by
+        centroid values of u_in
+
+        Centroid values of a and b provide diffusivity and rhs
+
+        Solution u is retruned in u_out
+        """
+
+        if u_out == None:
+            u_out = Quantity(self.domain)
+
+        
+        self.update_elliptic_matrix(a)        
+        self.update_elliptic_boundary_term(u_in)
+
+        # Pull out arrays and a matrix operator
+        A = self
+        rhs = b.centroid_values - self.boundary_term
+        x0 = u_in.centroid_values
+
+        x = conjugate_gradient(A,rhs,x0,imax=imax, tol=tol, atol=atol, iprint=iprint)
+
+        u_out.set_values(x, location='centroids')
+        u_out.set_boundary_values(u_in.boundary_values)
+
+        return u_out
+    
+
+    def parabolic_solve(self, u_in, b, a = None, u_out = None, \
+                       imax=10000, tol=1.0e-8, atol=1.0e-8, iprint=None):
+        """
+        Solve for u in the equation
+
+        ( I + dt div a grad ) u = b
+
+        u | boundary = g
+
+
+        u_in, u_out, f anf g are Quantity objects
+
+        Dirichlet BC g encoded into u_in boundary_values
+
+        Initial guess for iterative scheme is given by
+        centroid values of u_in
+
+        Centroid values of a and b provide diffusivity and rhs
+
+        Solution u is retruned in u_out
+
+        """
+
+        if u_out == None:
+            u_out = Quantity(self.domain)
+
+
+        self.update_elliptic_matrix(a)
+        self.update_elliptic_boundary_term(u_in)
+
+        self.set_parabolic_solve(True)
+
+
+        # Pull out arrays and a matrix operator
+        IdtA = self
+        rhs = b.centroid_values + (self.dt * self.boundary_term)
+        x0 = u_in.centroid_values
+
+        x = conjugate_gradient(IdtA,rhs,x0,imax=imax, tol=tol, atol=atol, iprint=iprint)
+
+        self.set_parabolic_solve(False)
+
+        u_out.set_values(x, location='centroids')
+        u_out.set_boundary_values(u_in.boundary_values)
+
+        return u_out
+
+ 

trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py

@@ +1 @@
+import operator
 from anuga import Domain
 from anuga import Quantity
@@ -104 +105 @@
         operator = self.operator1()
         domain = operator.domain
-        
-        operator.update_elliptic_matrix()
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+        
+        operator.update_elliptic_matrix(a)
 
         A = operator.elliptic_matrix
@@ -111 +116 @@
         assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
 
-        diffusivity = operator.diffusivity
-        diffusivity.set_values(10.0)
-        diffusivity.set_boundary_values(10.0)
-        
-        operator.update_elliptic_matrix()
+        a.set_values(10.0)
+        a.set_boundary_values(10.0)
+        
+        operator.update_elliptic_matrix(a)
 
         assert num.allclose(A.todense(), 10*num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
@@ -126 +130 @@
 
         domain = operator.domain
-        diffusivity = operator.diffusivity
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+        operator.update_elliptic_matrix(a)
 
         A = operator.elliptic_matrix
 
-        diffusivity.set_values(1.0)
-        diffusivity.set_boundary_values(1.0)
-        operator.update_elliptic_matrix()
-
+    
         A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],
                         [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]])
@@ -144 +149 @@
         assert num.allclose(A.todense(), A0+A1+A2)
 
-        diffusivity.set_values([2.0, 1.0], location = 'centroids')
-        diffusivity.set_boundary_values(1.0)
-        operator.update_elliptic_matrix()
+        a.set_values([2.0, 1.0], location = 'centroids')
+        a.set_boundary_values(1.0)
+        operator.update_elliptic_matrix(a)
 
         A = operator.elliptic_matrix
@@ -154 +159 @@
         assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:])
 
-        diffusivity.set_values([-2.0, -2.0], location = 'centroids')
-        diffusivity.set_boundary_values(1.0)
-        operator.update_elliptic_matrix()
+        a.set_values([-2.0, -2.0], location = 'centroids')
+        a.set_boundary_values(1.0)
+        operator.update_elliptic_matrix(a)
 
         assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:])
@@ -166 +171 @@
 
         print operator.apply_triangle_areas
-        
-        n = operator.n
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        operator.update_elliptic_matrix()
+
+        
+
+        q_in = Quantity(operator.domain)
+        q_in.set_values(1.0)
+        q_in.set_boundary_values(1.0)
+        
+        n = operator.n
+        
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+
+
+        q_1 = operator.elliptic_multiply(q_in)
+
+        q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in)
+
+        assert id(q_in) == id(q_2)
+
+        assert num.allclose(q_1.centroid_values,q_2.centroid_values)
+
+        assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values )
+
+        #Now have different boundary values
+
+        q_in.set_values(1.0)
+        q_in.set_boundary_values(0.0)
+        operator.update_elliptic_matrix(a)
+
+
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+
+
+        q_1 = operator.elliptic_multiply(q_in)
+
+        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )
+        
+    def test_elliptic_multiply_exclude_boundary_one_triangle(self):
+        operator = self.operator1()
+        operator.set_triangle_areas(False)
+
+        print operator.apply_triangle_areas
+        #n = operator.n
 
         q_in = Quantity(operator.domain)
@@ -174 +225 @@
         operator.update_elliptic_matrix()
 
-        
+
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+
+
+        q_1 = operator.elliptic_multiply(q_in, include_boundary=False)
+
+
+        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )
+
+    def test_elliptic_multiply_include_boundary_one_triangle(self):
+        operator = self.operator1()
+        operator.set_triangle_areas(True)
+
+        n = operator.n
+
+        q_in = Quantity(operator.domain)
+        q_in.set_values(1.0)
+        q_in.set_boundary_values(1.0)
+        
+        operator.update_elliptic_matrix()
+
+
         A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
 
@@ -180 +252 @@
         q_1 = operator.elliptic_multiply(q_in)
 
-        q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in)
+        q_2 = operator.elliptic_multiply(q_in, output = q_in)
 
         assert id(q_in) == id(q_2)
@@ -186 +258 @@
         assert num.allclose(q_1.centroid_values,q_2.centroid_values)
 
-        assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values )
+        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )
 
         #Now have different boundary values
@@ -200 +272 @@
         q_1 = operator.elliptic_multiply(q_in)
 
-        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )
-        
+        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )
+
     def test_elliptic_multiply_exclude_boundary_one_triangle(self):
         operator = self.operator1()
-        operator.set_triangle_areas(False)
-
-        print operator.apply_triangle_areas
-        #n = operator.n
+        operator.set_triangle_areas(True)
 
         q_in = Quantity(operator.domain)
@@ -218 +287 @@
 
 
-        q_1 = operator.elliptic_multiply(q_in, include_boundary=False)
-
-
-        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )
-
-    def test_elliptic_multiply_include_boundary_one_triangle(self):
-        operator = self.operator1()
-        operator.set_triangle_areas(True)
-
-        n = operator.n
-
-        q_in = Quantity(operator.domain)
-        q_in.set_values(1.0)
-        q_in.set_boundary_values(1.0)
+        q_1 = operator.elliptic_multiply(q_in)
+
+
+        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )
+
+    def test_mul_arg(self):
+        operator = self.operator1()
+
+        u = Quantity(operator.domain)
+        u.set_values(2.0)
+        #q boundary_values should equal 0.0
+
+
+
+        operator.update_elliptic_boundary_term(u)
+
+        r = 2.0
+
+        try:
+            q_out = operator * 2.0
+        except TypeError:
+            pass
+        else:
+            raise Exception('Should have caught an TypeError')
+
+
+    def test_mul(self):
+        operator = self.operator1()
+
+        u = Quantity(operator.domain)
+        u.set_values(2.0)
+        #q boundary_values should equal 0.0
+
         operator.update_elliptic_matrix()
 
-
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-
-
-        q_1 = operator.elliptic_multiply(q_in)
-
-        q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in)
-
-        assert id(q_in) == id(q_2)
-
-        assert num.allclose(q_1.centroid_values,q_2.centroid_values)
-
-        assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values )
-
-        #Now have different boundary values
-
-        q_in.set_values(1.0)
-        q_in.set_boundary_values(0.0)
-        operator.update_elliptic_matrix()
-
-
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-
-
-        q_1 = operator.elliptic_multiply(q_in)
-
-        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )
-
-    def test_elliptic_multiply_exclude_boundary_one_triangle(self):
-        operator = self.operator1()
-        operator.set_triangle_areas(True)
-
-        q_in = Quantity(operator.domain)
-        q_in.set_values(1.0)
-        q_in.set_boundary_values(1.0)
-        operator.update_elliptic_matrix()
-
-
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-
-
-        q_1 = operator.elliptic_multiply(q_in, include_boundary=False)
-
-
-        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )
-
-
-
-    def test_mul(self):
-        operator = self.operator1()
-
-        q = Quantity(operator.domain)
-        q.set_values(2.0)
-        #q boundary_values should equal 0.0
-
-        operator.build_elliptic_boundary_term()
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-        V1 = num.array([2.0]) #(uh)=2
+        operator.update_elliptic_boundary_term(u)
+
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+        V1 = num.array([2.0]) #u=2
         U1 = num.array([[2.0],[0.0],[0.0],[0.0]])
-        assert num.allclose(operator * q, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))
-
-    def test_cg_solve(self):
-        #cf self.test_mul()
-        operator1 = self.operator1()
-        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
-        operator1.set_qty_considered('u')
-        V = num.array([2.0]) #h=1, (uh)=2
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-        U = num.array([[2.0,2.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
-        test = 2*num.mat(A)*num.mat(U[:,0].reshape(4,1))
-        X = operator1.cg_solve(num.array(test).reshape(1,))
-        assert num.allclose(V, X)
-
-    def test_parabolic_solve(self):
-        operator1 = self.operator1()
-        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
-        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
-        U = num.array([[2.0,1.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
-        u = num.array([[2.0,1.0]])
-        U_new = operator1.parabolic_solver(u)
-        U_mod = num.array([[0.0,0.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
-        U_mod[0,:] = U_new
-        assert num.allclose(U_new - operator1.dt * 2 * num.mat(A)*num.mat(U_mod), U[0,:])
+
+        q_out = operator * u
+        
+        assert num.allclose(q_out.centroid_values, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))
+
+    def test_elliptic_solve_one_triangle(self):
+
+        operator = self.operator1()
+        n = operator.n
+        
+        U = num.array([2.0,2.0,1.0,1.0])
+
+        u_in = Quantity(operator.domain)
+        u_in.set_values(U[:1], location='centroids')
+        u_in.set_boundary_values(U[1:])
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        # Do this to get access to the matrix
+        # This is also called inside elliptic_solve
+        operator.update_elliptic_matrix(a)
+
+        V = num.array([2.0]) #h=1, u=2
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+        #U = num.array([[2.0,2.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
+
+        #Setup up rhs as b = A u
+        X = num.array(2*num.mat(A)*num.mat(U.reshape(4,1))).reshape(1,)
+        b = Quantity(operator.domain)
+        b.set_values(X, location='centroids')
+
+        u_in.set_values(0.0)
+
+        u_out = operator.elliptic_solve(u_in, b, a, iprint=1)
+
+        assert num.allclose(u_out.centroid_values, U[:n])
+
+
+    def test_elliptic_solve_two_triangle(self):
+
+        operator = self.operator2()
+        n = operator.n
+        
+        U = num.array([2.0,3.0,1.0,1.0,4.0,3.0])
+
+        u_in = Quantity(operator.domain)
+        u_in.set_values(U[:2], location='centroids')
+        u_in.set_boundary_values(U[2:])
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        # Do this to get access to the matrix
+        # This is also called inside elliptic_solve
+        operator.update_elliptic_matrix(a)
+
+        V1 = U[:n]
+        V2 = U[n:]
+
+        A = num.mat(operator.elliptic_matrix.todense())
+        U = num.mat(U.reshape(6,1))
+
+        #Setup up rhs as b = A u
+        X = num.array(2*A*U).reshape(2,)
+
+        b = Quantity(operator.domain)
+        b.set_values(X, location='centroids')
+
+        u_in.set_values(0.0)
+
+        u_out = operator.elliptic_solve(u_in, b, a, iprint=1)
+
+        assert num.allclose(u_out.centroid_values, V1)
+        assert num.allclose(u_out.boundary_values, V2)
+
+    def test_elliptic_solve_rectangular_cross(self):
+
+        from anuga import rectangular_cross_domain
+
+        m1 = 10
+        n1 = 10
+        domain = rectangular_cross_domain(m1,n1)
+
+        # Diffusivity
+        a = Quantity(domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        # Quantity to solve
+        u = Quantity(domain)
+        u.set_values(0.0)
+        u.set_boundary_values(1.0)
+
+        # Quantity for rhs
+        b = Quantity(domain)
+        b.set_values(0.0)
+        b.set_boundary_values(0.0)
+
+        operator = Kinematic_Viscosity(domain)
+
+        n = operator.n
+        tot_len = operator.tot_len
+
+        u_out = operator.elliptic_solve(u, b, a, iprint=1)
+    
+        assert num.allclose(u_out.centroid_values, num.ones_like(u_out.centroid_values))
+        assert num.allclose(u_out.boundary_values, num.ones_like(u_out.boundary_values))
+
+
+
+    def test_parabolic_solve_one_triangle(self):
+        operator = self.operator1()
+        n = operator.n
+        dt = operator.dt
+
+        U = num.array([2.0,2.0,1.0,1.0])
+        U_mod = num.array([10.0, 2.0, 1.0, 1.0])
+
+        u_in = Quantity(operator.domain)
+        u_in.set_values(U[:n], location='centroids')
+        u_in.set_boundary_values(U_mod[n:])
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+
+        V = num.array([2.0])
+        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+
+        #Setup up rhs
+        X = U_mod[:n] - dt*2*num.array(num.mat(A)*num.mat(U_mod.reshape(4,1))).reshape(n,)
+        b = Quantity(operator.domain)
+        b.set_values(X, location='centroids')
+
+
+        u_out = operator.parabolic_solve(u_in, b, a, iprint=1)
+
+
+        assert num.allclose(u_out.centroid_values, U_mod[:n])
+
+
+    def test_parabolic_solve_two_triangles(self):
+        operator = self.operator2()
+        n = operator.n
+        nt = operator.tot_len
+
+        dt = operator.dt
+
+        U = num.array([2.0,3.0,1.0,1.0,4.0,3.0])
+        U_mod = num.array([4.0,2.0,1.0,1.0,4.0,3.0])
+
+
+        u_in = Quantity(operator.domain)
+        u_in.set_values(U[:n], location='centroids')
+        u_in.set_boundary_values(U_mod[n:])
+
+        a = Quantity(operator.domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        operator.update_elliptic_matrix(a)
+
+
+        A = num.array([[-8.36656315,  3., 2.68328157,  2.68328157,  0.,  0. ],
+                       [ 3., -8.36656315 , 0. , 0. ,  2.68328157,  2.68328157]])
+
+        assert num.allclose(A,operator.elliptic_matrix.todense())
+
+
+
+        #Setup up rhs
+        X = U_mod[:n] - dt*2*num.array(num.mat(A)*num.mat(U_mod.reshape(nt,1))).reshape(n,)
+        b = Quantity(operator.domain)
+        b.set_values(X, location='centroids')
+
+
+        u_out = operator.parabolic_solve(u_in, b, a, iprint=1)
+
+
+        assert num.allclose(u_out.centroid_values, U_mod[:n])
+
+    def test_parabolic_solve_rectangular_cross(self):
+
+        from anuga import rectangular_cross_domain
+
+        m1 = 10
+        n1 = 10
+        domain = rectangular_cross_domain(m1,n1)
+
+        # Diffusivity
+        a = Quantity(domain)
+        a.set_values(1.0)
+        a.set_boundary_values(1.0)
+
+        # Quantity initial condition
+        u_in = Quantity(domain)
+        #u_in.set_values( 0.0 )
+        u_in.set_values(lambda x,y : 16.0*x*(1-x)*y*(1-y))
+        u_in.set_boundary_values(0.0)
+
+        # Quantity to solve
+        u_mod = Quantity(domain)
+        u_mod.set_values(lambda x,y : 15.9*x*(1-x)*y*(1-y) )
+        u_mod.set_boundary_values(0.0)
+
+        # Quantity for rhs
+        b = Quantity(domain)
+        b.set_values(0.0)
+        b.set_boundary_values(0.0)
+
+        operator = Kinematic_Viscosity(domain)
+
+        dt = 0.01
+        operator.dt = dt
+        n = operator.n
+        nt = operator.tot_len
+
+        operator.update_elliptic_matrix(a)
+
+        A = num.mat(operator.elliptic_matrix.todense())
+        D = num.mat(operator.triangle_areas.todense())
+        U_mod = num.concatenate( (u_mod.centroid_values, u_mod.boundary_values) )
+
+
+        #Setup up rhs
+        X = U_mod[:n] - dt*num.array(D*A*num.mat(U_mod.reshape(nt,1))).reshape(n,)
+        b = Quantity(operator.domain)
+        b.set_values(X, location='centroids')
+
+
+        u_out = operator.parabolic_solve(u_in, b, a, iprint=1)
+
+        assert num.allclose(u_out.centroid_values, U_mod[:n])
+
 
 ################################################################################

trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -512 +512 @@
                     (3, 1): 'Second',
                     (3, 2): 'Third',
-                    (0, 1): 'Internal'}
+                    (0, 1): 'Internal',
+                    (1, 2): 'Internal'}
 
         domain = Domain(points, vertices, boundary)
@@ -532 +533 @@
         domain.update_boundary()
         domain.check_integrity()
+
+    def test_boundary_conditions_inconsistent_internal(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        boundary = {(0, 0): 'Third',
+                    (0, 2): 'First',
+                    (2, 0): 'Second',
+                    (2, 1): 'Second',
+                    (3, 1): 'Second',
+                    (3, 2): 'Third',
+                    (0, 1): 'Internal'}
+
+        domain = Domain(points, vertices, boundary)
+
+        try:
+            domain.check_integrity()
+        except AssertionError:
+            pass
+        else:
+            msg = 'should have thrown assertion exception'
+            raise Exception(msg)
+
+
+
+
 
     def test_boundary_conditionsIII(self):
@@ -678 +712 @@
 
         assert num.allclose(domain.neighbours,
-                            [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
+                            [[-1,1,-2], [2,3,0], [-3,-4,1],[1,-5,-6]])
         assert num.allclose(domain.neighbour_edges,
                             [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])

trunk/anuga_core/source/anuga/utilities/cg_solve.py

@@ -8 +8 @@
 
 
-def conjugate_gradient(A,b,x0=None,imax=10000,tol=1.0e-8,atol=1.0e-14,iprint=None):
+def conjugate_gradient(A, b, x0=None, imax=10000, tol=1.0e-8, atol=1.0e-14, iprint=None):
     """
     Try to solve linear equation Ax = b using
@@ -23 +23 @@
 
     b  = num.array(b, dtype=num.float)
-    if len(b.shape) != 1 :
+    if len(b.shape) != 1:
        
         for i in range(b.shape[1]):
-            x0[:,i] = _conjugate_gradient(A, b[:,i], x0[:,i],
-                                          imax, tol, atol, iprint)
+            x0[:, i] = _conjugate_gradient(A, b[:, i], x0[:, i],
+                                           imax, tol, atol, iprint)
     else:
         x0 = _conjugate_gradient(A, b, x0, imax, tol, atol, iprint)
@@ -34 +34 @@
     
 def _conjugate_gradient(A, b, x0, 
-                        imax=10000, tol=1.0e-8, atol=1.0e-14, iprint=None):
-   """
+                        imax=10000, tol=1.0e-8, atol=1.0e-10, iprint=None):
+    """
    Try to solve linear equation Ax = b using
    conjugate gradient method
@@ -51 +51 @@
    """
 
+    b  = num.array(b, dtype=num.float)
+    if len(b.shape) != 1:
+        raise VectorShapeError, 'input vector should consist of only one column'
 
-   b  = num.array(b, dtype=num.float)
-   if len(b.shape) != 1 :
-      raise VectorShapeError, 'input vector should consist of only one column'
-
-   if x0 is None:
-      x0 = num.zeros(b.shape, dtype=num.float)
-   else:
-      x0 = num.array(x0, dtype=num.float)
+    if x0 is None:
+        x0 = num.zeros(b.shape, dtype=num.float)
+    else:
+        x0 = num.array(x0, dtype=num.float)
 
 
-   #FIXME: Should test using None
-   if iprint == None  or iprint == 0:
-      iprint = imax
+    if iprint == None  or iprint == 0:
+        iprint = imax
 
-   i=1
-   x = x0
-   r = b - A*x
-   d = r
-   rTr = num.dot(r,r)
-   rTr0 = rTr
+    i = 1
+    x = x0
+    r = b - A * x
+    d = r
+    rTr = num.dot(r, r)
+    rTr0 = rTr
 
-   #FIXME Let the iterations stop if starting with a small residual
-   while (i<imax and rTr>tol**2*rTr0 and rTr>atol**2 ):
-       q = A*d
-       alpha = rTr/num.dot(d,q)
-       x = x + alpha*d
-       if i%50 :
-           r = b - A*x
-       else:
-           r = r - alpha*q
-       rTrOld = rTr
-       rTr = num.dot(r,r)
-       bt = rTr/rTrOld
 
-       d = r + bt*d
-       i = i+1
-       if i%iprint == 0 :
-          log.info('i = %g rTr = %20.15e' %(i,rTr))
+    
+    #FIXME Let the iterations stop if starting with a small residual
+    while (i < imax and rTr > tol ** 2 * rTr0 and rTr > atol ** 2):
+        q = A * d
+        alpha = rTr / num.dot(d, q)
+        xold = x
+        x = x + alpha * d
 
-       if i==imax:
+        dx = num.linalg.norm(x-xold)
+        if dx < atol :
+            break
+            
+        if i % 50:
+            r = b - A * x
+        else:
+            r = r - alpha * q
+        rTrOld = rTr
+        rTr = num.dot(r, r)
+        bt = rTr / rTrOld
+
+        d = r + bt * d
+        i = i + 1
+        if i % iprint == 0:
+            log.info('i = %g rTr = %15.8e dx = %15.8e' % (i, rTr, dx))
+
+        if i == imax:
             log.warning('max number of iterations attained')
-            msg = 'Conjugate gradient solver did not converge: rTr==%20.15e' %rTr
+            msg = 'Conjugate gradient solver did not converge: rTr==%20.15e' % rTr
             raise ConvergenceError, msg
 
-   return x
+    #print x
+    return x
 

trunk/anuga_core/source/anuga/utilities/sparse.py

@@ -291 +291 @@
 
     def __repr__(self):
-        return '%d X %d sparse matrix:\n' %(self.M, self.N) + `self.data`
+        return '%d X %d sparse matrix:\n' %(self.M, self.N) + 'data '+ `self.data` + '\ncolind ' + \
+            `self.colind` + '\nrow_ptr ' + `self.row_ptr`
 
     def __len__(self):