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

[brad at geodynamics.org]

Author: brad
Date: 2010-05-29 15:02:06 -0700 (Sat, 29 May 2010)
New Revision: 16831

Modified:
   short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py
   short/3D/PyLith/trunk/playpen/faultpc/notes.tex
Log:
Updated notes.

Modified: short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py
===================================================================
--- short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py	2010-05-29 20:21:54 UTC (rev 16830)
+++ short/3D/PyLith/trunk/playpen/faultpc/checkfaultpc.py	2010-05-29 22:02:06 UTC (rev 16831)
@@ -35,38 +35,32 @@
 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
+P = J
+Pi = numpy.linalg.inv(P)
 
 # Compute condition number
-evals, evecs = numpy.linalg.eig(numpy.dot(P, J))
+evals, evecs = numpy.linalg.eig(numpy.dot(J, Pi))
 print numpy.abs(evals)
 print numpy.max(numpy.abs(evals))/numpy.min(numpy.abs(evals))
 
-# Compute preconditioner using diagonal approximations (but full Ai)
+# Compute preconditioner using diagonal approximations (but full A)
 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
+Pd[0:8,0:8] = A
+Pd[0:8,8:12] = C.transpose()
+Pd[8:12,0:8] = C
+Pd[8:12,8:12] = -1.0/(CACdi*CACdi)
 
+# Compute expected inverse of preconditioner
+Pdi = numpy.zeros(Pd.shape)
+Pdi[0:8,0:8] = Ai
+Pdi[0:8,8:12] = numpy.dot(-Ai, C.transpose())
+Pdi[8:12,0:8] = 0.0
+#Pdi[8:12,8:12] = -numpy.dot(numpy.dot(numpy.dot(C, Ai), C.transpose()), CACdi*CACdi)
+Pdi[8:12,8:12] = numpy.eye(4)
+
 # Compute condition number for diagonal approximations
-evals, evecs = numpy.linalg.eig(numpy.dot(Pd, J))
+evals, evecs = numpy.linalg.eig(numpy.dot(Pdi, J))
 print numpy.abs(evals)
 print numpy.max(numpy.abs(evals))/numpy.min(numpy.abs(evals))
 
-# Print preconditioner formed with diagonal approximations
-print "Pd00:", numpy.dot(numpy.dot(numpy.dot(numpy.dot(Aid, C.transpose()), -CACdi), C), Aid)
-print "Pd01:", numpy.dot(numpy.dot(-Aid, C.transpose()), -CACdi)
-print "Pd10:", numpy.dot(CACdi, numpy.dot(C, Aid))
-print "Pd11:", -CACdi
-
-# Print simplified terms for preconditioner formed with diagonal approximations
-CCti = numpy.linalg.inv(numpy.dot(C, C.transpose()))
-print "Pd00:", -numpy.dot(Aid, numpy.dot(numpy.dot(C.transpose(), CCti), C))
-print "Pd01:", 0.5*C.transpose()
-print "Pd10:", 0.5*C
-print "Pd11:", -CACdi
+print numpy.dot(Pdi, J)

Modified: short/3D/PyLith/trunk/playpen/faultpc/notes.tex
===================================================================
--- short/3D/PyLith/trunk/playpen/faultpc/notes.tex	2010-05-29 20:21:54 UTC (rev 16830)
+++ short/3D/PyLith/trunk/playpen/faultpc/notes.tex	2010-05-29 22:02:06 UTC (rev 16831)
@@ -40,6 +40,50 @@
     C & 0
   \end{array} \right).
 \end{equation}
+We provide PETSc with $P$ so that it can create $P^{-1}$. Using the
+field split preconditioner, we form
+\begin{equation}
+  P = \left( \begin{array}{cc}
+    A_\mathit{ml} & 0 \\
+    0 & P_f
+  \end{array} \right),
+\end{equation}
+where we use the ML package to form $A_\mathit{ml}$ and we create a
+custom matrix for the portion of the preconditioner associated with
+the Lagrange constraints, $P_f$. Let $n$ be the number of conventional
+degrees of freedom and $l$ be the number of Lagrange constraints. This
+means $A$ and $P$ are $(n+l) \times (n+l)$, $A$ and $A_\mathit{ml}$
+are $n \times n$, $C$ is $l \times n$, and $P_f$ is $l \times l$.
+
+Using the {\tt multiplicative} field split type, PETSc will form
+$P^{-1}$ as
+\begin{equation}
+  P^{-1} = \left( \begin{array}{cc}
+    A_\mathit{ml}^{-1} & -A_\mathit{ml}^{-1} C^T \\
+    0 & C A_\mathit{ml}^{-1} C^T P_f^{-1}
+  \end{array} \right).
+\end{equation}
+
+\subsection*{QUESTIONS FOR MATT}
+We are solving $J u = b$.
+\begin{itemize}
+\item Does PETSc solve the right-preconditioned system using
+  \begin{gather}
+    (J P^{-1}) v = b \\
+    u = P^{-1} v
+  \end{gather}
+  or the left-preconditioned system using
+  \begin{gather}
+    P^{-1}(J u - b) = 0.
+  \end{gather}
+\item
+  Do we want to reduce the conditioner number (ratio of largest to
+  smallest eigenvalue) of $J P^{-1}$ or $P^{-1} J$?
+\end{itemize}
+      
+
+\subsection*{Old, wrong stuff}
+
 Using the Schur complement, we want to form the preconditioner
 \begin{equation}
   P = \left( \begin{array}{cc}