Changeset 8154 for trunk/anuga_core/source/anuga/utilities/cg_solve.py

Timestamp:

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

Author:

steve

Message:

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

File:

• trunk/anuga_core/source/anuga/utilities/cg_solve.py (modified) (4 diffs)

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