[cig-commits] r16741 - in short/3D/PyLith/trunk: libsrc/faults libsrc/problems modulesrc/problems pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Tue May 18 22:30:16 PDT 2010
Author: brad
Date: 2010-05-18 22:30:16 -0700 (Tue, 18 May 2010)
New Revision: 16741
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
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/modulesrc/problems/Formulation.i
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Work on fault preconditioner.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-19 05:30:16 UTC (rev 16741)
@@ -453,7 +453,7 @@
INSERT_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numVertices*(spaceDim*spaceDim*2));
+ PetscLogFlops(spaceDim*spaceDim*2);
_logger->eventBegin(updateEvent);
#endif
@@ -537,8 +537,12 @@
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields)
{ // calcPreconditioner
+ assert(0 != pc);
+ assert(0 != jacobian);
+ assert(0 != fields);
+ assert(0 != _fields);
+ assert(0 != _logger);
-
/** We have J = [A C^T]
* [C 0]
*
@@ -575,7 +579,153 @@
* Term for DOF y of vertex L is:
* 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.
*/
+
+ const int setupEvent = _logger->eventId("FaPr setup");
+ const int computeEvent = _logger->eventId("FaPr compute");
+ const int restrictEvent = _logger->eventId("FaPr restrict");
+ const int updateEvent = _logger->eventId("FaPr update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell information and setup storage for cell data
+ const int spaceDim = _quadrature->spaceDim();
+ const int orientationSize = spaceDim * spaceDim;
+
+ // Allocate vectors for vertex values
+ double_array orientationVertex(orientationSize);
+ double_array jacobianVertexP(spaceDim*spaceDim);
+ double_array jacobianVertexN(spaceDim*spaceDim);
+ double_array jacobianInvVertexP(spaceDim);
+ double_array jacobianInvVertexN(spaceDim);
+ double_array precondVertexL(spaceDim);
+ int_array indicesN(spaceDim);
+ int_array indicesP(spaceDim);
+ int_array indicesRel(spaceDim);
+ for (int i=0; i < spaceDim; ++i)
+ indicesRel[i] = i;
+
+ // Get sections
+ const ALE::Obj<RealSection>& solutionSection = fields->solution().section();
+ assert(!solutionSection.isNull());
+
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+ solutionSection);
+ assert(!globalOrder.isNull());
+
+ // Get Jacobian matrix
+ 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);
+#endif
+
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+ const int v_negative = _cohesiveVertices[iVertex].negative;
+ const int v_positive = _cohesiveVertices[iVertex].positive;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(restrictEvent);
+#endif
+
+ // Get orientations at fault cell's vertices.
+ orientationSection->restrictPoint(v_fault, &orientationVertex[0],
+ orientationVertex.size());
+
+ // Set global order indices
+ indicesN = indicesRel + globalOrder->getIndex(v_negative);
+ indicesP = indicesRel + globalOrder->getIndex(v_positive);
+
+ err = MatGetValues(jacobianMatrix,
+ indicesN.size(), &indicesN[0],
+ indicesN.size(), &indicesN[0],
+ &jacobianVertexN[0]); CHECK_PETSC_ERROR(err);
+ err = MatGetValues(jacobianMatrix,
+ indicesP.size(), &indicesP[0],
+ indicesP.size(), &indicesP[0],
+ &jacobianVertexP[0]); CHECK_PETSC_ERROR(err);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Compute inverse of Jacobian diagonals
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ jacobianInvVertexN[iDim] = 1.0/jacobianVertexN[iDim*spaceDim+iDim];
+ jacobianInvVertexP[iDim] = 1.0/jacobianVertexP[iDim*spaceDim+iDim];
+ } // for
+
+ // Compute [C] [Adiag]^(-1) [C]^T
+ precondVertexL = 0.0;
+ for (int kDim=0; kDim < spaceDim; ++kDim)
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ precondVertexL[kDim] +=
+ orientationVertex[kDim*spaceDim+iDim] *
+ orientationVertex[kDim*spaceDim+iDim] *
+ (jacobianInvVertexN[iDim] + jacobianInvVertexP[iDim]);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Set global index associated with Lagrange constraint vertex.
+ const int indexLglobal = globalOrder->getIndex(v_lagrange);
+
+ // Translate global index into index in preconditioner.
+ const int indexLprecond = 0; // MATT- FIX THIS.
+
+ // Set diagonal entries in preconditioned matrix.
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ MatSetValue(precondMatrix,
+ indexLprecond + iDim,
+ indexLprecond + iDim,
+ precondVertexL[iDim],
+ INSERT_VALUES);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(spaceDim*spaceDim*4);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ } // 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);
+#endif
} // calcPreconditioner
// ----------------------------------------------------------------------
@@ -962,6 +1112,12 @@
_logger->registerEvent("FaIJ compute");
_logger->registerEvent("FaIJ restrict");
_logger->registerEvent("FaIJ update");
+
+ _logger->registerEvent("FaPr setup");
+ _logger->registerEvent("FaPr geometry");
+ _logger->registerEvent("FaPr compute");
+ _logger->registerEvent("FaPr restrict");
+ _logger->registerEvent("FaPr update");
} // initializeLogger
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-19 05:30:16 UTC (rev 16741)
@@ -32,6 +32,7 @@
_jacobian(0),
_jacobianLumped(0),
_fields(0),
+ _pc(0),
_isJacobianSymmetric(false)
{ // constructor
} // constructor
@@ -51,6 +52,12 @@
_jacobian = 0; // :TODO: Use shared pointer.
_jacobianLumped = 0; // :TODO: Use shared pointer.
_fields = 0; // :TODO: Use shared pointer.
+
+ PetscErrorCode err = 0;
+ if (0 != _pc) {
+ err = PCDestroy(_pc); _pc = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
} // deallocate
// ----------------------------------------------------------------------
@@ -70,20 +77,20 @@
} // splitFields
// ----------------------------------------------------------------------
-// Set flag for using custom fault preconditioner.
+// Set flag for using custom preconditioner for Lagrange constraints.
void
-pylith::problems::Formulation::useCustomFaultPC(const bool flag)
-{ // useCustomFaultPC
- _useCustomFaultPC = flag;
-} // useCustomFaultPC
+pylith::problems::Formulation::useCustomConstraintPC(const bool flag)
+{ // useCustomConstraintPC
+ _useCustomConstraintPC = flag;
+} // useCustomConstraintPC
// ----------------------------------------------------------------------
-// Get flag indicating use of custom fault conditioner.
+// Get flag indicating use of custom conditioner for Lagrange constraints.
bool
-pylith::problems::Formulation::useCustomFaultPC(void) const
-{ // useCustomFaultPC
- return _useCustomFaultPC;
-} // useCustomFaultPC
+pylith::problems::Formulation::useCustomConstraintPC(void) const
+{ // useCustomConstraintPC
+ return _useCustomConstraintPC;
+} // useCustomConstraintPC
// ----------------------------------------------------------------------
// Return the fields
@@ -128,6 +135,17 @@
} // submeshIntegrators
// ----------------------------------------------------------------------
+// Set handle to preconditioner.
+void
+pylith::problems::Formulation::preconditioner(PetscPC& pc)
+{ // preconditioner
+ _pc = pc;
+
+ PetscErrorCode err = 0;
+ err = PetscObjectReference((PetscObject) _pc); CHECK_PETSC_ERROR(err);
+} // preconditioner
+
+// ----------------------------------------------------------------------
// Initialize formulation.
void
pylith::problems::Formulation::initialize(void)
@@ -342,18 +360,15 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
-
-#if 0
- // Recalculate preconditioner.
- if (_useCustomFaultPC) {
+ if (0 != _pc) {
+ // Recalculate preconditioner.
numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _meshIntegrators[i]->calcPreconditioner(pc, _jacobian, _fields);
+ _meshIntegrators[i]->calcPreconditioner(&_pc, _jacobian, _fields);
numIntegrators = _submeshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
- _submeshIntegrators[i]->calcPreconditioner(pc, _jacobian, _fields);
+ _submeshIntegrators[i]->calcPreconditioner(&_pc, _jacobian, _fields);
} // if
-#endif
} // reformJacobian
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-05-19 05:30:16 UTC (rev 16741)
@@ -67,17 +67,17 @@
*/
bool splitFields(void) const;
- /** Set flag for using custom fault preconditioner.
+ /** Set flag for using custom preconditioner for Lagrange constraints.
*
- * @param flag True if splitting fields, false otherwise.
+ * @param flag True if using custom fault preconditioner, false otherwise.
*/
- void useCustomFaultPC(const bool flag);
+ void useCustomConstraintPC(const bool flag);
- /** Get flag indicating use of custom fault conditioner.
+ /** Get flag indicating use of custom conditioner for Lagrange constraints.
*
- * @returns flag True if using custom fault preconditioner, false otherwise.
+ * @returns True if using custom fault preconditioner, false otherwise.
*/
- bool useCustomFaultPC(void) const;
+ bool useCustomConstraintPC(void) const;
/** Get solution fields.
*
@@ -107,6 +107,12 @@
void submeshIntegrators(IntegratorSubMesh** integrators,
const int numIntegrators);
+ /** Set handle to preconditioner.
+ *
+ * @param pc PETSc preconditioner.
+ */
+ void preconditioner(PetscPC& pc);
+
/// Initialize formulation.
virtual
void initialize(void);
@@ -191,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.
+ PetscPC _pc; ///< Handle to PETSc preconditioner.
topology::Field<topology::Mesh>* _jacobianLumped; ///< Handle to lumped Jacobian of system.
topology::SolutionFields* _fields; ///< Handle to solution fields for system.
@@ -203,8 +209,10 @@
bool _isJacobianSymmetric; ///< Is system Jacobian symmetric?
bool _splitFields; ///< True if splitting fields.
- bool _useCustomFaultPC; ///< True if using custom fault preconditioner.
+ /// True if using custom preconditioner for Lagrange constraints.
+ bool _useCustomConstraintPC;
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2010-05-19 05:30:16 UTC (rev 16741)
@@ -30,6 +30,10 @@
#endif
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Constructor
pylith::problems::Solver::Solver(void) :
_formulation(0),
@@ -72,30 +76,35 @@
Formulation* const formulation,
const topology::SolutionFields& fields)
{ // _setupFieldSplit
+ assert(0 != pc);
+ assert(0 != formulation);
PetscErrorCode err = 0;
+ const ALE::Obj<SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
const topology::Field<topology::Mesh>& solution = fields.solution();
- const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh =
- fields.mesh().sieveMesh();
+ const ALE::Obj<RealSection>& solutionSection = solution.section();
+ assert(!solutionSection.isNull());
+
err = PCSetType(*pc, PCFIELDSPLIT); CHECK_PETSC_ERROR(err);
err = PCSetOptionsPrefix(*pc, "fs_"); CHECK_PETSC_ERROR(err);
err = PCSetFromOptions(*pc); CHECK_PETSC_ERROR(err);
#if defined(FIELD_SPLIT)
- constructFieldSplit(solution.section(),
+ constructFieldSplit(solutionSection,
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- solution.section()),
+ solutionSection),
solution.vector(), *pc);
- if (formulation->useCustomFaultPC()) {
+ if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
+ formulation->useCustomConstraintPC()) {
PetscKSP *ksps = 0;
PetscMat A, P;
PetscInt num, m, n;
err = PCFieldSplitGetSubKSP(*pc, &num, &ksps); CHECK_PETSC_ERROR(err);
- // Put in PC matrix for fault
+ // Put in PC matrix for additional space (fault).
MatStructure flag;
err = KSPGetOperators(ksps[num-1], &A,
PETSC_NULL, &flag); CHECK_PETSC_ERROR(err);
@@ -106,9 +115,7 @@
err = MatSetSizes(P, m, n,
PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
- // Allocate just the diagonal. This should really go in the Fault
- // object that calculates the preconditioner since it is
- // implementation dependent.
+ // Allocate just the diagonal.
err = MatSeqAIJSetPreallocation(P, 1, PETSC_NULL); CHECK_PETSC_ERROR(err);
err = MatMPIAIJSetPreallocation(P, 1, PETSC_NULL,
0, PETSC_NULL); CHECK_PETSC_ERROR(err);
@@ -117,6 +124,9 @@
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);
} // if
#endif
Modified: short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-05-19 05:30:16 UTC (rev 16741)
@@ -46,18 +46,19 @@
*/
bool splitFields(void) const;
- /** Set flag for using custom fault preconditioner.
+ /** Set flag for using custom Lagrange constraint preconditioner.
*
- * @param flag True if splitting fields, false otherwise.
+ * @param flag True if using custom constraint precondition,
+ * false otherwise.
*/
- void useCustomFaultPC(const bool flag);
+ void useCustomConstraintPC(const bool flag);
- /** Get flag indicating use of custom fault conditioner.
+ /** Get flag indicating use of custom Lagrange constraint conditioner.
*
- * @returns flag True if using custom fault preconditioner,
+ * @returns flag True if using custom constraint preconditioner,
* false otherwise.
*/
- bool useCustomFaultPC(void) const;
+ bool useCustomConstraintPC(void) const;
/** Get solution fields.
*
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-05-18 23:15:13 UTC (rev 16740)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-05-19 05:30:16 UTC (rev 16741)
@@ -58,10 +58,10 @@
## \b Properties
## @li \b matrix_type Type of PETSc sparse matrix.
## @li \b split_fields Split solution fields into displacements
- ## @li \b use_custom_pc Use custom preconditioner for fault.
- ## and Lagrange multipliers for separate preconditioning.
+ ## @li \b use_custom_constraint_pc Use custom preconditioner for
+ ## Lagrange constraints.
## @li \b view_jacobian Flag to output Jacobian matrix when it is
- ## reformed.
+ ## reformed.
##
## \b Facilities
## @li \b time_step Time step size manager.
@@ -78,8 +78,10 @@
useSplitFields.meta['tip'] = "Split solution fields into displacements "\
"and Lagrange multipliers for separate preconditioning."
- useCustomFaultPC = pyre.inventory.bool("use_custom_fault_pc", default=True)
- useCustomFaultPC.meta['tip'] = "Use custom preconditioner for fault."
+ useCustomConstraintPC = pyre.inventory.bool("use_custom_constraint_pc",
+ default=True)
+ useCustomConstraintPC.meta['tip'] = "Use custom preconditioner for " \
+ "Lagrange constraints."
viewJacobian = pyre.inventory.bool("view_jacobian", default=False)
viewJacobian.meta['tip'] = "Write Jacobian matrix to binary file."
@@ -293,13 +295,16 @@
import journal
self._debug = journal.debug(self.name)
- if self.useCustomFaultPC and not self.inventory.useSplitFields:
- print "WARNING: Request to use custom fault preconditioner without " \
- "splitting fields. Setting split fields flag to 'True'."
+ if self.inventory.useCustomConstraintPC and \
+ not self.inventory.useSplitFields:
+ print "WARNING: Request to use custom preconditioner for Lagrange " \
+ "constraints without splitting fields. " \
+ "Setting split fields flag to 'True'."
self.inventory.useSplitFields = True
ModuleFormulation.splitFields(self, self.inventory.useSplitFields)
- ModuleFormulation.useCustomFaultPC(self, self.inventory.useCustomFaultPC)
+ ModuleFormulation.useCustomConstraintPC(self,
+ self.inventory.useCustomConstraintPC)
return
More information about the CIG-COMMITS
mailing list