[cig-commits] r15597 - in short/3D/PyLith/trunk: libsrc/faults libsrc/problems libsrc/topology modulesrc/problems pylith/topology
brad at geodynamics.org
brad at geodynamics.org
Wed Aug 26 09:38:56 PDT 2009
Author: brad
Date: 2009-08-26 09:38:55 -0700 (Wed, 26 Aug 2009)
New Revision: 15597
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh
short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc
short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh
short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc
short/3D/PyLith/trunk/modulesrc/problems/SolverLinear.i
short/3D/PyLith/trunk/modulesrc/problems/SolverNonlinear.i
short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
Log:
Fixed setting KSP operators (reuse preconditioner if sparse matrix values haven't changed).
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2009-08-26 16:38:55 UTC (rev 15597)
@@ -136,6 +136,8 @@
*firstFaultVertex, *firstFaultCell, useLagrangeConstraints());
} else {
+ //std::cout << "BEFORE ADJUSTING TOPOLOGY FOR FAULT '" << label() << "' firstFaultVertex: " << *firstFaultVertex << ", firstFaultCell: " << *firstFaultCell << std::endl;
+
const int faultDim = 2;
assert(3 == mesh->dimension());
@@ -154,6 +156,8 @@
assert(!groupField.isNull());
CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
*firstFaultVertex, *firstFaultCell, useLagrangeConstraints());
+
+ //std::cout << "AFTER ADJUSTING TOPOLOGY FOR FAULT '" << label() << "' firstFaultVertex: " << *firstFaultVertex << ", firstFaultCell: " << *firstFaultCell << std::endl;
} // if/else
} // adjustTopology
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2009-08-26 16:38:55 UTC (rev 15597)
@@ -98,19 +98,26 @@
void
pylith::problems::SolverLinear::solve(
topology::Field<topology::Mesh>* solution,
- const topology::Jacobian& jacobian,
+ topology::Jacobian* jacobian,
const topology::Field<topology::Mesh>& residual)
{ // solve
assert(0 != solution);
+ assert(0 != jacobian);
PetscErrorCode err = 0;
// Update PetscVector view of field.
residual.scatterSectionToVector();
- const PetscMat jacobianMat = jacobian.matrix();
- err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
- SAME_NONZERO_PATTERN); CHECK_PETSC_ERROR(err);
+ const PetscMat jacobianMat = jacobian->matrix();
+ if (!jacobian->valuesChanged()) {
+ err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
+ SAME_PRECONDITIONER); CHECK_PETSC_ERROR(err);
+ } else {
+ err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
+ DIFFERENT_NONZERO_PATTERN); CHECK_PETSC_ERROR(err);
+ } // else
+ jacobian->resetValuesChanged();
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 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh 2009-08-26 16:38:55 UTC (rev 15597)
@@ -62,7 +62,7 @@
* @param residual Residual field.
*/
void solve(topology::Field<topology::Mesh>* solution,
- const topology::Jacobian& jacobian,
+ topology::Jacobian* jacobian,
const topology::Field<topology::Mesh>& residual);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2009-08-26 16:38:55 UTC (rev 15597)
@@ -87,7 +87,7 @@
void
pylith::problems::SolverNonlinear::solve(
topology::Field<topology::Mesh>* solution,
- const topology::Jacobian& jacobian,
+ topology::Jacobian* jacobian,
const topology::Field<topology::Mesh>& residual)
{ // solve
assert(0 != solution);
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh 2009-08-26 16:38:55 UTC (rev 15597)
@@ -64,7 +64,7 @@
* @param residual Residual field.
*/
void solve(topology::Field<topology::Mesh>* solveSoln,
- const topology::Jacobian& jacobian,
+ topology::Jacobian* jacobian,
const topology::Field<topology::Mesh>& residual);
/** Generic C interface for reformResidual for integration with
Modified: short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc 2009-08-26 16:38:55 UTC (rev 15597)
@@ -26,7 +26,8 @@
const char* matrixType,
const bool blockOkay) :
_fields(fields),
- _matrix(0)
+ _matrix(0),
+ _valuesChanged(true)
{ // constructor
const ALE::Obj<Mesh::SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
const ALE::Obj<Mesh::RealSection>& solnSection = fields.solution().section();
@@ -97,6 +98,8 @@
} else
throw std::runtime_error("Unknown mode for assembly of sparse matrix "
"associated with system Jacobian.");
+
+ _valuesChanged = true;
} // assemble
// ----------------------------------------------------------------------
@@ -106,6 +109,7 @@
{ // zero
PetscErrorCode err = MatZeroEntries(_matrix);
CHECK_PETSC_ERROR(err);
+ _valuesChanged = true;
} // zero
// ----------------------------------------------------------------------
@@ -189,5 +193,22 @@
throw std::runtime_error("Jacobian matrix is not symmetric.");
} // verifySymmetry
+// ----------------------------------------------------------------------
+// Get flag indicating if sparse matrix values have been
+// updated.
+bool
+pylith::topology::Jacobian::valuesChanged(void) const
+{ // valuesChanged
+ return _valuesChanged;
+} // valuesChanged
+// ----------------------------------------------------------------------
+// Reset flag indicating if sparse matrix values have been updated.
+void
+pylith::topology::Jacobian::resetValuesChanged(void)
+{ // resetValuesChanged
+ _valuesChanged = false;
+} // resteValuesChanged
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh 2009-08-26 16:38:55 UTC (rev 15597)
@@ -82,12 +82,25 @@
/// Verify symmetry of matrix. For debugger purposes only.
void verifySymmetry(void) const;
+ /** Get flag indicating if sparse matrix values have been
+ * updated. This pertains to the values being changed, NOT the
+ * pattern of nonzero entries.
+ *
+ * @returns True if values have been updated/altered.
+ */
+ bool valuesChanged(void) const;
+
+ /// Reset flag indicating if sparse matrix values have been updated.
+ void resetValuesChanged(void);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
const SolutionFields& _fields; ///< Solution fields associated with problem.
PetscMat _matrix; ///< Sparse matrix for Jacobian of problem.
+ bool _valuesChanged; ///< Sparse matrix values have been updated.
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc 2009-08-26 16:38:55 UTC (rev 15597)
@@ -58,6 +58,9 @@
const ALE::Obj<SieveMesh::label_type>& materialsLabel =
sieveMesh->getLabel("material-id");
+ //mesh.view("===== MESH BEFORE CHECKING MATERIALS =====");
+ //materialsLabel->view("material-id");
+
int* matBegin = materialIds;
int* matEnd = materialIds + numMaterials;
std::sort(matBegin, matEnd);
Modified: short/3D/PyLith/trunk/modulesrc/problems/SolverLinear.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/SolverLinear.i 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/modulesrc/problems/SolverLinear.i 2009-08-26 16:38:55 UTC (rev 15597)
@@ -52,7 +52,7 @@
* @param residual Residual field.
*/
void solve(pylith::topology::Field<pylith::topology::Mesh>* solution,
- const pylith::topology::Jacobian& jacobian,
+ pylith::topology::Jacobian* jacobian,
const pylith::topology::Field<pylith::topology::Mesh>& residual);
}; // SolverLinear
Modified: short/3D/PyLith/trunk/modulesrc/problems/SolverNonlinear.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/SolverNonlinear.i 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/modulesrc/problems/SolverNonlinear.i 2009-08-26 16:38:55 UTC (rev 15597)
@@ -52,7 +52,7 @@
* @param residual Residual field.
*/
void solve(pylith::topology::Field<pylith::topology::Mesh>* solution,
- const pylith::topology::Jacobian& jacobian,
+ pylith::topology::Jacobian* jacobian,
const pylith::topology::Field<pylith::topology::Mesh>& residual);
}; // SolverNonlinear
Modified: short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2009-08-26 03:47:39 UTC (rev 15596)
+++ short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2009-08-26 16:38:55 UTC (rev 15597)
@@ -91,10 +91,14 @@
logEvent = "%sadjTopo" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ #self._info.activate()
+ #mesh.view("===== MESH BEFORE ADJUSTING TOPOLOGY =====")
+
if not interfaces is None:
firstFaultVertex = 0
firstFaultCell = 0
for interface in interfaces:
+ self._info.log("Counting vertices for fault '%s'." % interface.label())
nvertices = interface.numVertices(mesh)
firstFaultCell += nvertices
if interface.useLagrangeConstraints():
@@ -105,6 +109,9 @@
(interface.label(), nvertices))
firstFaultVertex, firstFaultCell = \
interface.adjustTopology(mesh, firstFaultVertex, firstFaultCell)
+
+ #mesh.view("===== MESH AFTER ADJUSTING TOPOLOGY =====")
+ #self._info.deactivate()
self._eventLogger.eventEnd(logEvent)
return
More information about the CIG-COMMITS
mailing list