[cig-commits] r16810 - in short/3D/PyLith/trunk: . examples/twocells/twoquad4 libsrc/faults libsrc/problems

[brad at geodynamics.org]

Author: brad
Date: 2010-05-27 13:21:23 -0700 (Thu, 27 May 2010)
New Revision: 16810

Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/trunk/libsrc/problems/Solver.cc
Log:
Switch default behavior back to diagonal fault preconditioner. Fix for diagonal preconditioner. Want -(CAiC^T)^(-1) not CAiC^T.

Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2010-05-27 18:59:55 UTC (rev 16809)
+++ short/3D/PyLith/trunk/TODO	2010-05-27 20:21:23 UTC (rev 16810)
@@ -107,11 +107,10 @@
   Write up description of Savage and Prescott (1978) benchmark and
   distribute to Greg Lyzenga and Jay Parker
 
+    Add time step information
+
 CODE
 
-  * Fix parallel bugs
-    + Need to limit MatGetValues to local.
-
   * Better preconditioning [Matt/Brad]
 
     Diagonal preconditioner does not work very well. It looks like the

Modified: short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg	2010-05-27 18:59:55 UTC (rev 16809)
+++ short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg	2010-05-27 20:21:23 UTC (rev 16810)
@@ -96,3 +96,22 @@
 
 # start_in_debugger = true
 # debugger_timeout = 100
+
+# Field split
+[pylithapp.timedependent.formulation]
+split_fields = True
+use_custom_constraint_pc = True
+matrix_type = aij
+
+[pylithapp.petsc]
+fs_pc_type = fieldsplit
+fs_pc_fieldsplit_real_diagonal =
+fs_pc_fieldsplit_type = multiplicative
+fs_fieldsplit_0_pc_type = ml
+fs_fieldsplit_1_pc_type = ml
+fs_fieldsplit_2_pc_type = ilu
+#fs_fieldsplit_2_pc_type = jacobi
+fs_fieldsplit_0_ksp_type = gmres
+fs_fieldsplit_1_ksp_type = gmres
+fs_fieldsplit_2_ksp_type = richardson
+

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc	2010-05-27 18:59:55 UTC (rev 16809)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc	2010-05-27 20:21:23 UTC (rev 16810)
@@ -575,7 +575,7 @@
   /** We have J = [A C^T]
    *              [C   0]
    *
-   * We want to approximate C A^(-1) C^T.
+   * We want to approximate -( C A^(-1) C^T)^(-1)
    *
    * Consider Lagrange vertex L that constrains the relative
    * displacement between vertex N on the negative side of the fault
@@ -603,16 +603,16 @@
    * 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)
+   * -1.0 / (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)
+   * -1.0 / (R10^2 (Ai_nx + Ai_px) + R11^2 (Ai_ny + Ai_py))
    *
    * Translate DOF for global vertex L into DOF in split field for
    * preconditioner.
    */
 
-#if 0 // DIAGONAL PRECONDITIONER
+#if 1 // DIAGONAL PRECONDITIONER
   const int setupEvent = _logger->eventId("FaPr setup");
   const int computeEvent = _logger->eventId("FaPr compute");
   const int restrictEvent = _logger->eventId("FaPr restrict");
@@ -712,12 +712,17 @@
     //   Adiag^{-1}_{ii} = jacobianInvVertexN[i] + jacobianInvVertexP[i] (BRAD: Are you sure its not a minus sign here?)
     //  \sum_{j} C_{ij} Adiag^{-1}_{jj} C^T_{ji}
     precondVertexL = 0.0;
-    for (int kDim=0; kDim < spaceDim; ++kDim)
+    for (int kDim=0; kDim < spaceDim; ++kDim) {
+      double value = 0.0;
       for (int iDim=0; iDim < spaceDim; ++iDim)
-        precondVertexL[kDim] += 
+	value += 
           orientationVertex[kDim*spaceDim+iDim] * 
           orientationVertex[kDim*spaceDim+iDim] * 
           (jacobianInvVertexN[iDim] + jacobianInvVertexP[iDim]);
+      if (fabs(value) > 1.0e-8)
+	precondVertexL[kDim] = -1.0 / value;
+    } // for
+    
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(computeEvent);
@@ -1078,12 +1083,25 @@
     _logger->eventBegin(updateEvent);
 #endif
 
+#if 0
     err = MatSetValues(*precondMatrix,
                  indicesLagrange.size(), &indicesLagrange[0],
                  indicesLagrange.size(), &indicesLagrange[0],
                  &preconditionerCell[0],
-                 INSERT_VALUES);
+                 ADD_VALUES);
     CHECK_PETSC_ERROR_MSG(err, "Setting values in fault preconditioner failed.");
+#else
+    for (int iLagrange=0; iLagrange < numBasis; ++iLagrange)
+      for (int iDim=0; iDim < spaceDim; ++iDim) {
+	const int iL = iLagrange*spaceDim + iDim;
+	const int indexL = indicesLagrange[iL];
+	MatSetValue(*precondMatrix,
+		    indexL, indexL,
+		    preconditionerCell[iL*nrowsF+iL],
+		    ADD_VALUES);
+      } // for
+    
+#endif
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(updateEvent);

Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.cc	2010-05-27 18:59:55 UTC (rev 16809)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc	2010-05-27 20:21:23 UTC (rev 16810)
@@ -122,7 +122,7 @@
     err = MatSetType(*precondMatrix, MATAIJ);
     err = MatSetFromOptions(*precondMatrix); CHECK_PETSC_ERROR(err);
     
-#if 0
+#if 1
     // Allocate just the diagonal.
     err = MatSeqAIJSetPreallocation(*precondMatrix, 1, 
 				    PETSC_NULL); CHECK_PETSC_ERROR(err);