[cig-commits] r16735 - short/3D/PyLith/trunk/libsrc/problems

[brad at geodynamics.org]

Author: brad
Date: 2010-05-17 17:30:19 -0700 (Mon, 17 May 2010)
New Revision: 16735

Modified:
   short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
Log:
Added comments on fault preconditioner.

Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc	2010-05-17 09:07:25 UTC (rev 16734)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc	2010-05-18 00:30:19 UTC (rev 16735)
@@ -23,7 +23,7 @@
 
 #include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
 
-#define FIELD_SPLIT
+//#define FIELD_SPLIT
 #define NEW_FAULT_PRECONDITIONER
 
 #if defined(FIELD_SPLIT)
@@ -113,6 +113,45 @@
       // Create a mapping to indices for that space (might be in FS)
       err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
     } // if
+
+    /** We have J = [A C^T]
+     *              [C   0]
+     *
+     * We want to approximate C A^(-1) C^T.
+     *
+     * Consider Lagrange vertex L that constrains the relative
+     * displacement between vertex N on the negative side of the fault
+     * and vertex P on the positive side of the fault.
+     *
+     * If we approximate A(-1) by 1/diag(A), then we can write 
+     * C A^(-1) C^T for a 2-D case as
+     *
+     * [-R00 -R01  R00 R01][Ai_nx 0      0     0    ][-R00 -R10]
+     * [-R10 -R11  R10 R11][      Ai_ny  0     0    ][-R01 -R11]
+     *                     [      0      Ai_px 0    ][ R00  R10]
+     *                     [                   Ai_py][ R01  R11]
+     *
+     * where
+     *
+     * Ai_nx is the inverse of the diag(A) for DOF x of vertex N
+     * Ai_ny is the inverse of the diag(A) for DOF y of vertex N
+     * Ai_px is the inverse of the diag(A) for DOF x of vertex P
+     * Ai_py is the inverse of the diag(A) for DOF y of vertex P
+     *
+     * If Ai_nx == Ai_ny and Ai_px == Ai_py, then the result is
+     * diagonal. Otherwise, the off-diagonal terms will be nonzero,
+     * but we expect them to be small. Since we already approximate
+     * the inverse of A by the inverse of the diagonal, we drop the
+     * off-diagonal terms of C A^(-1) C^T:
+     *
+     * Term for DOF x of vertex L is: 
+     * R00^2 (Ai_nx + Ai_px) + R01^2 (Ai_ny + Ai_py)
+     *
+     * Term for DOF y of vertex L is: 
+     * R10^2 (Ai_nx + Ai_px) + R11^2 (Ai_ny + Ai_py)
+     *
+     */
+
 #endif
 #endif
   }