[cig-commits] r16750 - in short/3D/PyLith/trunk/libsrc: faults feassemble problems
brad at geodynamics.org
brad at geodynamics.org
Wed May 19 15:47:45 PDT 2010
Author: brad
Date: 2010-05-19 15:47:45 -0700 (Wed, 19 May 2010)
New Revision: 16750
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
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/Solver.hh
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh
Log:
Revised fault preconditioner stuff for proper initialization.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-19 22:47:45 UTC (rev 16750)
@@ -533,11 +533,11 @@
// operator.
void
pylith::faults::FaultCohesiveLagrange::calcPreconditioner(
- PetscPC* const pc,
+ PetscMat* const precondMatrix,
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields)
{ // calcPreconditioner
- assert(0 != pc);
+ assert(0 != precondMatrix);
assert(0 != jacobian);
assert(0 != fields);
assert(0 != _fields);
@@ -626,17 +626,6 @@
const PetscMat jacobianMatrix = jacobian->matrix();
assert(0 != jacobianMatrix);
- // Get preconditioner matrix
- PetscKSP *ksps = 0;
- PetscMat precondMatrix = 0;
- PetscInt numKSP = 0;
- PetscErrorCode err = 0;
- err = PCFieldSplitGetSubKSP(*pc, &numKSP, &ksps); CHECK_PETSC_ERROR(err);
-
- MatStructure flag;
- err = KSPGetOperators(ksps[numKSP-1], PETSC_NULL,
- &precondMatrix, &flag); CHECK_PETSC_ERROR(err);
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
@@ -661,6 +650,8 @@
indicesN = indicesRel + globalOrder->getIndex(v_negative);
indicesP = indicesRel + globalOrder->getIndex(v_positive);
+ PetscErrorCode err = 0;
+
err = MatGetValues(jacobianMatrix,
indicesN.size(), &indicesN[0],
indicesN.size(), &indicesN[0],
@@ -698,6 +689,7 @@
// Set global index associated with Lagrange constraint vertex.
const int indexLglobal = globalOrder->getIndex(v_lagrange);
+#if 0
// Translate global index into index in preconditioner.
const int indexLprecond = 0; // MATT- FIX THIS.
@@ -708,6 +700,7 @@
indexLprecond + iDim,
precondVertexL[iDim],
INSERT_VALUES);
+#endif
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(spaceDim*spaceDim*4);
@@ -716,12 +709,6 @@
} // for
- // Flush assembled portion.
- MatAssemblyBegin(precondMatrix,MAT_FLUSH_ASSEMBLY);
- MatAssemblyEnd(precondMatrix,MAT_FLUSH_ASSEMBLY);
-
- err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
-
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numVertices*(spaceDim*spaceDim*4));
_logger->eventEnd(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-05-19 22:47:45 UTC (rev 16750)
@@ -166,7 +166,7 @@
* @param fields Solution fields
*/
virtual
- void calcPreconditioner(PetscPC* const pc,
+ void calcPreconditioner(PetscMat* const precondMatrix,
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-05-19 22:47:45 UTC (rev 16750)
@@ -25,7 +25,7 @@
#include "pylith/topology/topologyfwd.hh" // USES Mesh, Field, SolutionFields
#include "pylith/utils/utilsfwd.hh" // HOLDSA EventLogger
-#include "pylith/utils/petscfwd.h" // USES PetscPC
+#include "pylith/utils/petscfwd.h" // USES PetscMat
#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES GravityField
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
@@ -237,12 +237,12 @@
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
*
- * @param pc PETSc preconditioner structure.
+ * @param precondMatrix Custom preconditioning matrix.
* @param jacobian Sparse matrix for Jacobian of system.
* @param fields Solution fields
*/
virtual
- void calcPreconditioner(PetscPC* const pc,
+ void calcPreconditioner(PetscMat* const precondMatrix,
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-05-19 22:47:45 UTC (rev 16750)
@@ -158,7 +158,7 @@
inline
void
pylith::feassemble::Integrator<quadrature_type>::calcPreconditioner(
- PetscPC* const pc,
+ PetscMat* const precondMatrix,
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields) {
} // calcPreconditioner
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-19 22:47:45 UTC (rev 16750)
@@ -32,7 +32,7 @@
_jacobian(0),
_jacobianLumped(0),
_fields(0),
- _pc(0),
+ _precondMatrix(0),
_isJacobianSymmetric(false)
{ // constructor
} // constructor
@@ -54,8 +54,8 @@
_fields = 0; // :TODO: Use shared pointer.
PetscErrorCode err = 0;
- if (0 != _pc) {
- err = PCDestroy(_pc); _pc = 0;
+ if (0 != _precondMatrix) {
+ err = MatDestroy(_precondMatrix); _precondMatrix = 0;
CHECK_PETSC_ERROR(err);
} // if
} // deallocate
@@ -137,12 +137,12 @@
// ----------------------------------------------------------------------
// Set handle to preconditioner.
void
-pylith::problems::Formulation::preconditioner(PetscPC& pc)
+pylith::problems::Formulation::customPCMatrix(PetscMat& mat)
{ // preconditioner
- _pc = pc;
+ _precondMatrix = mat;
PetscErrorCode err = 0;
- err = PetscObjectReference((PetscObject) _pc); CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) _precondMatrix); CHECK_PETSC_ERROR(err);
} // preconditioner
// ----------------------------------------------------------------------
@@ -360,14 +360,22 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
- if (0 != _pc) {
+ if (0 != _precondMatrix) {
// Recalculate preconditioner.
numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _meshIntegrators[i]->calcPreconditioner(&_pc, _jacobian, _fields);
+ _meshIntegrators[i]->calcPreconditioner(&_precondMatrix,
+ _jacobian, _fields);
numIntegrators = _submeshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _submeshIntegrators[i]->calcPreconditioner(&_pc, _jacobian, _fields);
+ _submeshIntegrators[i]->calcPreconditioner(&_precondMatrix,
+ _jacobian, _fields);
+
+ // Flush assembled portion.
+ MatAssemblyBegin(_precondMatrix,MAT_FLUSH_ASSEMBLY);
+ MatAssemblyEnd(_precondMatrix,MAT_FLUSH_ASSEMBLY);
+ MatAssemblyBegin(_precondMatrix,MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(_precondMatrix,MAT_FINAL_ASSEMBLY);
} // if
} // reformJacobian
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-05-19 22:47:45 UTC (rev 16750)
@@ -111,7 +111,7 @@
*
* @param pc PETSc preconditioner.
*/
- void preconditioner(PetscPC& pc);
+ void customPCMatrix(PetscMat& mat);
/// Initialize formulation.
virtual
@@ -197,7 +197,7 @@
double _t; ///< Current time (nondimensional).
double _dt; ///< Current time step (nondimensional).
topology::Jacobian* _jacobian; ///< Handle to Jacobian of system.
- PetscPC _pc; ///< Handle to PETSc preconditioner.
+ PetscMat _precondMatrix; ///< Custom PETSc preconditioning matrix.
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 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2010-05-19 22:47:45 UTC (rev 16750)
@@ -73,10 +73,12 @@
// Setup preconditioner for preconditioning with split fields.
void
pylith::problems::Solver::_setupFieldSplit(PetscPC* const pc,
+ PetscMat* const precondMatrix,
Formulation* const formulation,
const topology::SolutionFields& fields)
{ // _setupFieldSplit
assert(0 != pc);
+ assert(0 != precondMatrix);
assert(0 != formulation);
PetscErrorCode err = 0;
@@ -96,37 +98,38 @@
solutionSection),
solution.vector(), *pc);
- if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
+ const int spaceDim = sieveMesh->getDimension();
+ if (solutionSection->getNumSpaces() > spaceDim &&
formulation->useCustomConstraintPC()) {
- PetscKSP *ksps = 0;
- PetscMat A, P;
- PetscInt num, m, n;
+ // Get total number of DOF associated with constraints field split
+ const ALE::Obj<RealSection>& constraintSection =
+ solutionSection->getFibration(spaceDim);
+ assert(!constraintSection.isNull());
- err = PCFieldSplitGetSubKSP(*pc, &num, &ksps); CHECK_PETSC_ERROR(err);
+ // :TODO: Is this assembled across processors?
+ const int ndof = constraintSection->size();
+
+ if (0 != *precondMatrix) {
+ err = MatDestroy(*precondMatrix); *precondMatrix = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
+ PetscInt nrows = ndof;
+ PetscInt ncols = ndof;
- // Put in PC matrix for additional space (fault).
- MatStructure flag;
- err = KSPGetOperators(ksps[num-1], &A,
- PETSC_NULL, &flag); CHECK_PETSC_ERROR(err);
+ err = MatCreate(sieveMesh->comm(), precondMatrix); CHECK_PETSC_ERROR(err);
+ err = MatSetSizes(*precondMatrix, nrows, ncols,
+ PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
+ err = MatSetType(*precondMatrix, MATSBAIJ);
+ err = MatSetFromOptions(*precondMatrix); CHECK_PETSC_ERROR(err);
- err = PetscObjectReference((PetscObject) A); CHECK_PETSC_ERROR(err);
- err = MatGetLocalSize(A, &m, &n); CHECK_PETSC_ERROR(err);
- err = MatCreate(sieveMesh->comm(), &P); CHECK_PETSC_ERROR(err);
- err = MatSetSizes(P, m, n,
- PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
-
// Allocate just the diagonal.
- err = MatSeqAIJSetPreallocation(P, 1, PETSC_NULL); CHECK_PETSC_ERROR(err);
- err = MatMPIAIJSetPreallocation(P, 1, PETSC_NULL,
+ err = MatSeqAIJSetPreallocation(*precondMatrix, 1,
+ PETSC_NULL); CHECK_PETSC_ERROR(err);
+ err = MatMPIAIJSetPreallocation(*precondMatrix, 1, PETSC_NULL,
0, PETSC_NULL); CHECK_PETSC_ERROR(err);
- err = MatSetFromOptions(P); CHECK_PETSC_ERROR(err);
- err = KSPSetOperators(ksps[num-1], A, P, flag); CHECK_PETSC_ERROR(err);
-
- err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
-
- // Set preconditioner in formulation
- formulation->preconditioner(*pc);
+ // Set preconditioning matrix in formulation
+ formulation->customPCMatrix(*precondMatrix);
} // if
#endif
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2010-05-19 22:47:45 UTC (rev 16750)
@@ -25,7 +25,7 @@
#include "pylith/topology/topologyfwd.hh" // USES SolutionFields
#include "pylith/utils/utilsfwd.hh" // USES EventLogger
-#include "pylith/utils/petscfwd.h" // USES PetscPC
+#include "pylith/utils/petscfwd.h" // USES PetscMat
// Solver ---------------------------------------------------------
/** @brief Abstract C++ base class for using PETSc linear and
@@ -64,14 +64,16 @@
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
- /** Setup preconditioner for preconditioning with split fields.8
+ /** Setup preconditioner for preconditioning using split fields.
*
* @param pc PETSc preconditioner.
+ * @param precondMatrix Matrix for custom preconditioner.
* @param formulation Formulation of system of equations.
* @param fields Solution fields.
*/
void
_setupFieldSplit(PetscPC* const pc,
+ PetscMat* const precondMatrix,
Formulation* const formulation,
const topology::SolutionFields& fields);
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2010-05-19 22:47:45 UTC (rev 16750)
@@ -24,9 +24,14 @@
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverLinear::SolverLinear(void) :
- _ksp(0)
+ _ksp(0),
+ _precondMatrix(0)
{ // constructor
} // constructor
@@ -48,6 +53,11 @@
PetscErrorCode err = KSPDestroy(_ksp); _ksp = 0;
CHECK_PETSC_ERROR(err);
} // if
+
+ if (0 != _precondMatrix) {
+ PetscErrorCode err = MatDestroy(_precondMatrix); _precondMatrix = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
} // deallocate
// ----------------------------------------------------------------------
@@ -75,7 +85,7 @@
if (formulation->splitFields()) {
PetscPC pc = 0;
err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
- _setupFieldSplit(&pc, formulation, fields);
+ _setupFieldSplit(&pc, &_precondMatrix, formulation, fields);
} // if
} // initialize
@@ -113,6 +123,34 @@
} // else
jacobian->resetValuesChanged();
+ // Update KSP operators with custom preconditioner if necessary.
+ const ALE::Obj<RealSection>& solutionSection = solution->section();
+ assert(!solutionSection.isNull());
+ const ALE::Obj<SieveMesh>& sieveMesh = solution->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
+ 0 != _precondMatrix) {
+
+ PetscPC pc = 0;
+ PetscKSP *ksps = 0;
+ PetscMat A = 0;
+ PetscInt num = 0;
+
+ err = KSPSetUp(_ksp); CHECK_PETSC_ERROR(err);
+ err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
+ err = PCFieldSplitGetSubKSP(pc, &num, &ksps); CHECK_PETSC_ERROR(err);
+ assert(solutionSection->getNumSpaces() == num);
+
+ 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, _precondMatrix,
+ flag); CHECK_PETSC_ERROR(err);
+ err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
+ } // if
+
+
const PetscVec residualVec = residual.vector();
const PetscVec solutionVec = solution->vector();
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh 2010-05-19 20:19:01 UTC (rev 16749)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh 2010-05-19 22:47:45 UTC (rev 16750)
@@ -79,6 +79,7 @@
private :
PetscKSP _ksp; ///< PETSc KSP linear solver.
+ PetscMat _precondMatrix; ///< Preconditioning matrix for Lagrange constraints.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list