[cig-commits] r20085 - short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems

[brad at geodynamics.org]

Author: brad
Date: 2012-05-11 14:17:14 -0700 (Fri, 11 May 2012)
New Revision: 20085

Modified:
   short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
   short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc
Log:
Merge from trunk (v1.7-trunk).

Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc	2012-05-11 21:15:58 UTC (rev 20084)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc	2012-05-11 21:17:14 UTC (rev 20085)
@@ -160,12 +160,64 @@
   assert(!solutionSection.isNull());
   const int spaceDim = sieveMesh->getDimension();
   const int numSpaces = solutionSection->getNumSpaces();
+  const bool separateComponents = formulation->splitFieldComponents();
 
   err = PCSetType(*pc, PCFIELDSPLIT); CHECK_PETSC_ERROR(err);
   err = PCSetOptionsPrefix(*pc, "fs_"); CHECK_PETSC_ERROR(err);
   err = PCSetFromOptions(*pc); CHECK_PETSC_ERROR(err);
 
-  bool separateComponents = formulation->splitFieldComponents();
+  if (formulation->splitFields() && 
+      formulation->useCustomConstraintPC() &&
+      numSpaces > spaceDim) {
+    // We have split fields with a custom constraint preconditioner
+    // and constraints exist.
+
+    // Get total number of DOF associated with constraints field split
+    const ALE::Obj<SieveMesh::order_type>& lagrangeGlobalOrder =
+      sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "faultDefault",
+                                              solutionSection, spaceDim);
+    assert(!lagrangeGlobalOrder.isNull());
+
+    err = MatDestroy(&_jacobianPCFault); CHECK_PETSC_ERROR(err);
+    PylithInt nrows = lagrangeGlobalOrder->getLocalSize();
+    PylithInt ncols = nrows;
+
+    err = MatCreate(sieveMesh->comm(), &_jacobianPCFault); CHECK_PETSC_ERROR(err);
+    err = MatSetSizes(_jacobianPCFault, nrows, ncols, 
+		      PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
+    err = MatSetType(_jacobianPCFault, MATAIJ);
+    err = MatSetFromOptions(_jacobianPCFault); CHECK_PETSC_ERROR(err);
+    
+    // Allocate just the diagonal.
+    err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1, 
+				    PETSC_NULL); CHECK_PETSC_ERROR(err);
+    err = MatMPIAIJSetPreallocation(_jacobianPCFault, 1, PETSC_NULL, 
+				    0, PETSC_NULL); CHECK_PETSC_ERROR(err);
+    // Set preconditioning matrix in formulation
+    formulation->customPCMatrix(_jacobianPCFault);
+
+    assert(_jacobianPCFault);
+
+    _ctx.pc = *pc;
+    _ctx.A = jacobian.matrix();
+    switch (spaceDim) {
+    case 1 :
+      _ctx.faultFieldName = "1";
+      break;
+    case 2 :
+      _ctx.faultFieldName = (separateComponents) ? "2" : "1";
+      break;
+    case 3 :
+      _ctx.faultFieldName = (separateComponents) ? "3" : "1";
+      break;
+    default:
+      assert(0);
+      throw std::logic_error("Unknown space dimension in "
+			     "Problems::_setupFieldSplit().");
+    } // switch
+    _ctx.faultA = _jacobianPCFault;
+  } // if
+
   if (separateComponents) {
     PetscMat* precon = new PetscMat[numSpaces];
     for (int i=0; i < numSpaces; ++i) {
@@ -256,58 +308,6 @@
     err = MatNullSpaceDestroy(&nullsp[0]);CHECK_PETSC_ERROR(err);
     delete[] fields;
   } // if/else
-
-  if (formulation->splitFields() && 
-      formulation->useCustomConstraintPC() &&
-      numSpaces > spaceDim) {
-    // We have split fields with a custom constraint preconditioner
-    // and constraints exist.
-
-    // Get total number of DOF associated with constraints field split
-    const ALE::Obj<SieveMesh::order_type>& lagrangeGlobalOrder =
-      sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "faultDefault",
-                                              solutionSection, spaceDim);
-    assert(!lagrangeGlobalOrder.isNull());
-
-    err = MatDestroy(&_jacobianPCFault); CHECK_PETSC_ERROR(err);
-    PylithInt nrows = lagrangeGlobalOrder->getLocalSize();
-    PylithInt ncols = nrows;
-
-    err = MatCreate(sieveMesh->comm(), &_jacobianPCFault); CHECK_PETSC_ERROR(err);
-    err = MatSetSizes(_jacobianPCFault, nrows, ncols, 
-		      PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
-    err = MatSetType(_jacobianPCFault, MATAIJ);
-    err = MatSetFromOptions(_jacobianPCFault); CHECK_PETSC_ERROR(err);
-    
-    // Allocate just the diagonal.
-    err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1, 
-				    PETSC_NULL); CHECK_PETSC_ERROR(err);
-    err = MatMPIAIJSetPreallocation(_jacobianPCFault, 1, PETSC_NULL, 
-				    0, PETSC_NULL); CHECK_PETSC_ERROR(err);
-    // Set preconditioning matrix in formulation
-    formulation->customPCMatrix(_jacobianPCFault);
-
-    assert(_jacobianPCFault);
-
-    _ctx.pc = *pc;
-    _ctx.A = jacobian.matrix();
-    switch (spaceDim) {
-    case 1 :
-      _ctx.faultFieldName = "1";
-      break;
-    case 2 :
-      _ctx.faultFieldName = (separateComponents) ? "2" : "1";
-      break;
-    case 3 :
-      _ctx.faultFieldName = (separateComponents) ? "3" : "1";
-      break;
-    default:
-      assert(0);
-      throw std::logic_error("Unknown space dimension in "
-			     "Problems::_setupFieldSplit().");
-    } // switch
-    _ctx.faultA = _jacobianPCFault;
-  } // if
 } // _setupFieldSplit
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc	2012-05-11 21:15:58 UTC (rev 20084)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc	2012-05-11 21:17:14 UTC (rev 20085)
@@ -130,20 +130,23 @@
       _jacobianPCFault) {
     
     PetscPC pc = 0;
-    PetscKSP *ksps = 0;
-    PetscMat A = 0;
-    PylithInt num = 0;
+    PCCompositeType ctype;
     
     err = KSPSetUp(_ksp); CHECK_PETSC_ERROR(err);
     err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
-    err = PCFieldSplitGetSubKSP(pc, &num, &ksps); CHECK_PETSC_ERROR(err);
-    // Now only true if splitting components assert(solutionSection->getNumSpaces() == num);
+    err = PCFieldSplitGetType(pc, &ctype); CHECK_PETSC_ERROR(err);
+    if (ctype != PC_COMPOSITE_SCHUR) {
+      PetscKSP    *ksps = 0;
+      PylithInt    num  = 0;
+      PetscMat     A    = 0;
+      MatStructure flag;
 
-    MatStructure flag;
-    err = KSPGetOperators(ksps[num-1], &A, PETSC_NULL, &flag);CHECK_PETSC_ERROR(err);
-    err = PetscObjectReference((PetscObject) A);CHECK_PETSC_ERROR(err);
-    err = KSPSetOperators(ksps[num-1], A, _jacobianPCFault, flag);CHECK_PETSC_ERROR(err);
-    err = PetscFree(ksps);CHECK_PETSC_ERROR(err);
+      err = PCFieldSplitGetSubKSP(pc, &num, &ksps); CHECK_PETSC_ERROR(err);
+      err = KSPGetOperators(ksps[num-1], &A, PETSC_NULL, &flag);CHECK_PETSC_ERROR(err);
+      err = PetscObjectReference((PetscObject) A);CHECK_PETSC_ERROR(err);
+      err = KSPSetOperators(ksps[num-1], A, _jacobianPCFault, flag);CHECK_PETSC_ERROR(err);
+      err = PetscFree(ksps);CHECK_PETSC_ERROR(err);
+    }
   } // if
 #endif