[cig-commits] r16813 - short/3D/PyLith/trunk/playpen/faultpc

[brad at geodynamics.org]

Author: brad
Date: 2010-05-27 14:16:10 -0700 (Thu, 27 May 2010)
New Revision: 16813

Modified:
   short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py
Log:
Real check of fault preconditioner.

Modified: short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py
===================================================================
--- short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py	2010-05-27 20:45:43 UTC (rev 16812)
+++ short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py	2010-05-27 21:16:10 UTC (rev 16813)
@@ -1,36 +1,57 @@
+# Jacobian matrix for twocells/twoquad4/dislocation.cfg.
+
 import numpy
 import numpy.linalg as linalg
 
-Afull = numpy.array([[1,1, 1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 1,1, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 0,0, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [1,1, 0,0, 1,1, 1,1, 0,0, 0,0, 0,0, 0,0],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 0,0, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 0,0, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1, 1,1],
-                     [0,0, 0,0, 0,0, 0,0, 1,1, 1,1, 1,1, 1,1]],
-                    dtype=numpy.float64)
+A = numpy.array([[4/3.0, -0.5, 1/6.0, 0,   0,0,0,0],
+                 [-0.5, 4/3.0, 0, -5/6.0,  0,0,0,0],
+                 [1/6.0, 0, 4/3.0, 0.5,    0,0,0,0],
+                 [0, -5/6.0, 0.5, 4/3.0,   0,0,0,0],
+                 [0,0,0,0,  4/3.0, 0.5, 1/6.0, 0],
+                 [0,0,0,0,  0.5, 4/3.0, 0, -5/6.0],
+                 [0,0,0,0, 1/6.0, 0, 4/3.0, -0.5],
+                 [0,0,0,0, 0, -5/6.0, -0.5, 4/3.0]],
+                dtype=numpy.float64)
+Ai = linalg.inv(A)
 
-#Ainv = linalg.inv(Afull)
-
-C = numpy.array([[0,0, 0,-1, 0,0, 0,0, 0,1, 0,0, 0,0, 0,0],
-                 [0,0, -1,0, 0,0, 0,0, 1,0, 0,0, 0,0, 0,0],
-                 [0,0, 0,0, 0,-1, 0,0, 0,0, 0,1, 0,0, 0,0],
-                 [0,0, 0,0, -1,0, 0,0, 0,0, 1,0, 0,0, 0,0],
-                 [0,0, 0,0, 0,0, 0,-1, 0,0, 0,0, 0,1, 0,0],
-                 [0,0, 0,0, 0,0, -1,0, 0,0, 0,0, 1,0, 0,0]],
+C = numpy.array([[0,-1,  0,0,  0,1, 0,0],
+                 [-1,0,  0,0,  1,0, 0,0],
+                 [0,0,  0,-1,  0,0, 0,1],
+                 [0,0, -1,0,   0,0, 1,0]],
                 dtype=numpy.float64)
-CAC = numpy.dot(numpy.dot(C, Afull), C.transpose())
+Z = numpy.zeros( (4,4), dtype=numpy.float64)
 
+J = numpy.vstack( (numpy.hstack( (A, C.transpose()) ),
+                   numpy.hstack( (C, Z) ) ) )
+Jinv = linalg.inv(J)
 
-print "Afull",Afull
-#print "Afull inverse",Ainv
-print "CACfull",CAC
+# Compute [C] [A]^(-1) [C]^T and its inverse.
+CAC = numpy.dot(numpy.dot(C, Ai), C.transpose())
+CACi = numpy.linalg.inv(CAC)
+
+# Compute diagonal approximation of CAC and its inverse
+Aid = 1.0 / A.diagonal() * numpy.identity(A.shape[0])
+CACd = numpy.dot(numpy.dot(C, Aid), C.transpose())
+CACdi = 1.0 / CACd.diagonal() * numpy.identity(CACd.shape[0])
+
+# Compute preconditioner using full matrices (no approximations)
+P = numpy.zeros(J.shape)
+P[0:8,0:8] = Ai + numpy.dot(numpy.dot(numpy.dot(numpy.dot(Ai, C.transpose()), -CACi), C), Ai)
+P[0:8,8:12] = numpy.dot(numpy.dot(-Ai, C.transpose()), CACi)
+P[8:12,0:8] = numpy.dot(CACi, numpy.dot(C, Ai))
+P[8:12,8:12] = -CACi
+
+# Compute preconditioner using diagonal approximations (but full Ai)
+Pd = numpy.zeros(J.shape)
+Pd[0:8,0:8] = Ai + numpy.dot(numpy.dot(numpy.dot(numpy.dot(Aid, C.transpose()), -CACdi), C), Aid)
+Pd[0:8,8:12] = numpy.dot(numpy.dot(-Aid, C.transpose()), CACdi)
+Pd[8:12,0:8] = numpy.dot(CACdi, numpy.dot(C, Aid))
+Pd[8:12,8:12] = -CACdi
+
+evals, evecs = numpy.linalg.eig(numpy.dot(P, J))
+print numpy.abs(evals)
+print numpy.max(numpy.abs(evals))/numpy.min(numpy.abs(evals))
+
+evals, evecs = numpy.linalg.eig(numpy.dot(Pd, J))
+print numpy.abs(evals)
+print numpy.max(numpy.abs(evals))/numpy.min(numpy.abs(evals))