[cig-commits] r18245 - short/3D/PyLith/trunk/libsrc/problems
brad at geodynamics.org
brad at geodynamics.org
Tue Apr 19 09:00:47 PDT 2011
Author: brad
Date: 2011-04-19 09:00:47 -0700 (Tue, 19 Apr 2011)
New Revision: 18245
Modified:
short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
short/3D/PyLith/trunk/libsrc/problems/Solver.cc
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
Log:
Fixed order of SNES preconditioning setup calls. Renabled use of custom preconditioner for linear solver.
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2011-04-19 14:12:37 UTC (rev 18244)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2011-04-19 16:00:47 UTC (rev 18245)
@@ -38,7 +38,7 @@
_jacobian(0),
_jacobianLumped(0),
_fields(0),
- _precondMatrix(0),
+ _customConstraintPCMat(0),
_isJacobianSymmetric(false)
{ // constructor
} // constructor
@@ -61,12 +61,12 @@
#if 0 // :KLUDGE: Assume Solver deallocates matrix.
PetscErrorCode err = 0;
- if (0 != _precondMatrix) {
- err = PetscObjectDereference((PetscObject) _precondMatrix);
- _precondMatrix = 0; CHECK_PETSC_ERROR(err);
+ if (0 != _customConstraintPCMat) {
+ err = PetscObjectDereference((PetscObject) _customConstraintPCMat);
+ _customConstraintPCMat = 0; CHECK_PETSC_ERROR(err);
} // if
#else
- _precondMatrix = 0;
+ _customConstraintPCMat = 0;
#endif
} // deallocate
@@ -149,7 +149,7 @@
void
pylith::problems::Formulation::customPCMatrix(PetscMat& mat)
{ // preconditioner
- _precondMatrix = mat;
+ _customConstraintPCMat = mat;
#if 0 // :KLUDGE: Assume solver deallocates matrix
PetscErrorCode err = 0;
@@ -319,23 +319,23 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
- if (0 != _precondMatrix) {
+ if (0 != _customConstraintPCMat) {
// Recalculate preconditioner.
numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _meshIntegrators[i]->calcPreconditioner(&_precondMatrix,
+ _meshIntegrators[i]->calcPreconditioner(&_customConstraintPCMat,
_jacobian, _fields);
numIntegrators = _submeshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _submeshIntegrators[i]->calcPreconditioner(&_precondMatrix,
+ _submeshIntegrators[i]->calcPreconditioner(&_customConstraintPCMat,
_jacobian, _fields);
- MatAssemblyBegin(_precondMatrix, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_precondMatrix, MAT_FINAL_ASSEMBLY);
+ MatAssemblyBegin(_customConstraintPCMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(_customConstraintPCMat, MAT_FINAL_ASSEMBLY);
#if 0 // debugging
std::cout << "Preconditioner Matrix" << std::endl;
- MatView(_precondMatrix, PETSC_VIEWER_STDOUT_WORLD);
+ MatView(_customConstraintPCMat, PETSC_VIEWER_STDOUT_WORLD);
#endif
} // if
} // reformJacobian
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2011-04-19 14:12:37 UTC (rev 18244)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2011-04-19 16:00:47 UTC (rev 18245)
@@ -192,7 +192,7 @@
double _t; ///< Current time (nondimensional).
double _dt; ///< Current time step (nondimensional).
topology::Jacobian* _jacobian; ///< Handle to Jacobian of system.
- PetscMat _precondMatrix; ///< Custom PETSc preconditioning matrix.
+ PetscMat _customConstraintPCMat; ///< Custom PETSc preconditioning matrix for constraints.
topology::Field<topology::Mesh>* _jacobianLumped; ///< Handle to lumped Jacobian of system.
topology::SolutionFields* _fields; ///< Handle to solution fields for system.
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2011-04-19 14:12:37 UTC (rev 18244)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2011-04-19 16:00:47 UTC (rev 18245)
@@ -102,6 +102,35 @@
{ // initialize
assert(0 != formulation);
_formulation = formulation;
+
+ // Make global preconditioner matrix
+ PetscMat jacobianMat = jacobian.matrix();
+
+ const ALE::Obj<RealSection>& solutionSection = fields.solution().section();
+ assert(!solutionSection.isNull());
+ const ALE::Obj<SieveMesh>& sieveMesh = fields.solution().mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ if (formulation->splitFields() &&
+ formulation->useCustomConstraintPC() &&
+ solutionSection->getNumSpaces() > sieveMesh->getDimension()) {
+ // We have split fields with a custom constraint preconditioner
+ // and constraints exist.
+
+ PetscInt M, N, m, n;
+ PetscErrorCode err = 0;
+ err = MatGetSize(jacobianMat, &M, &N);CHECK_PETSC_ERROR(err);
+ err = MatGetLocalSize(jacobianMat, &m, &n);CHECK_PETSC_ERROR(err);
+ err = MatCreateShell(fields.mesh().comm(), m, n, M, N, &_ctx, &_jacobianPC);
+ CHECK_PETSC_ERROR(err);
+ err = MatShellSetOperation(_jacobianPC, MATOP_GET_SUBMATRIX,
+ (void (*)(void)) MyMatGetSubMatrix);
+ CHECK_PETSC_ERROR(err);
+ } else {
+ _jacobianPC = jacobianMat;
+ PetscErrorCode err = PetscObjectReference((PetscObject) jacobianMat);
+ CHECK_PETSC_ERROR(err);
+ } // if/else
} // initialize
// ----------------------------------------------------------------------
@@ -133,8 +162,12 @@
solution.vector(), *pc);
const int spaceDim = sieveMesh->getDimension();
- if (solutionSection->getNumSpaces() > spaceDim &&
- formulation->useCustomConstraintPC()) {
+ if (formulation->splitFields() &&
+ formulation->useCustomConstraintPC() &&
+ solutionSection->getNumSpaces() > sieveMesh->getDimension()) {
+ // 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",
@@ -171,27 +204,11 @@
// Set preconditioning matrix in formulation
formulation->customPCMatrix(_jacobianPCFault);
- // Make global preconditioner matrix
- PetscMat jacobianMat = jacobian.matrix();
-
- PetscInt M, N, m, n;
- err = MatGetSize(jacobianMat, &M, &N);CHECK_PETSC_ERROR(err);
- err = MatGetLocalSize(jacobianMat, &m, &n);CHECK_PETSC_ERROR(err);
- err = MatCreateShell(fields.mesh().comm(), m, n, M, N, &_ctx, &_jacobianPC);
- CHECK_PETSC_ERROR(err);
- err = MatShellSetOperation(_jacobianPC, MATOP_GET_SUBMATRIX,
- (void (*)(void)) MyMatGetSubMatrix);
- CHECK_PETSC_ERROR(err);
- _ctx.A = jacobianMat;
+ _ctx.pc = *pc;
+ _ctx.A = jacobian.matrix();
_ctx.faultFieldName = "3";
_ctx.faultA = _jacobianPCFault;
- } else {
- // Make global preconditioner matrix
- PetscMat jacobianMat = jacobian.matrix();
- _jacobianPC = jacobianMat;
- err = PetscObjectReference((PetscObject) jacobianMat);CHECK_PETSC_ERROR(err);
- } // if/else
-
+ } // if
} // _setupFieldSplit
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2011-04-19 14:12:37 UTC (rev 18244)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2011-04-19 16:00:47 UTC (rev 18245)
@@ -85,7 +85,6 @@
if (formulation->splitFields()) {
PetscPC pc = 0;
err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
- _ctx.pc = pc;
_setupFieldSplit(&pc, formulation, jacobian, fields);
} // if
} // initialize
@@ -129,9 +128,10 @@
assert(!solutionSection.isNull());
const ALE::Obj<SieveMesh>& sieveMesh = solution->mesh().sieveMesh();
assert(!sieveMesh.isNull());
-#if 0
+
+#if 1 // :TODO: Does this get replaced by the use of the MatShell like SNES?
if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
- 0 != _jacobianPCFault) {
+ _jacobianPCFault) {
PetscPC pc = 0;
PetscKSP *ksps = 0;
@@ -153,6 +153,7 @@
} // if
#endif
+
const PetscVec residualVec = residual.vector();
const PetscVec solutionVec = solution->vector();
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-04-19 14:12:37 UTC (rev 18244)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-04-19 16:00:47 UTC (rev 18245)
@@ -125,7 +125,6 @@
PetscPC pc = 0;
err = SNESGetKSP(_snes, &ksp); CHECK_PETSC_ERROR(err);
err = KSPGetPC(ksp, &pc); CHECK_PETSC_ERROR(err);
- _ctx.pc = pc;
_setupFieldSplit(&pc, formulation, jacobian, fields);
} // if
} // initialize
More information about the CIG-COMMITS
mailing list