[cig-commits] r16740 - in short/3D/PyLith/trunk: . libsrc/faults libsrc/feassemble libsrc/meshio libsrc/problems libsrc/utils modulesrc/problems pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Tue May 18 16:15:13 PDT 2010
Author: brad
Date: 2010-05-18 16:15:13 -0700 (Tue, 18 May 2010)
New Revision: 16740
Modified:
short/3D/PyLith/trunk/TODO
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/meshio/DataWriterVTK.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/libsrc/problems/Solver.hh
short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
short/3D/PyLith/trunk/libsrc/utils/petscfwd.h
short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Reorganization of field split stuff. Started working on hooks for computing custom fault preconditioner.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/TODO 2010-05-18 23:15:13 UTC (rev 16740)
@@ -122,6 +122,14 @@
Need field split working for both SolverLinear and SolverNonlinear.
+ Preconditioner Object
+
+ Owned by formulation?
+ Solver::initialize() sets pc object in Preconditioner
+
+ use_field_split
+ use_custom_fault_pc
+
* Optimization [Brad]
+ Specialized elasticity integrator objects for Tri3, Quad4, Hex8
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -529,6 +529,56 @@
} // integrateJacobianAssembled
// ----------------------------------------------------------------------
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator.
+void
+pylith::faults::FaultCohesiveLagrange::calcPreconditioner(
+ PetscPC* const pc,
+ topology::Jacobian* const jacobian,
+ topology::SolutionFields* const fields)
+{ // calcPreconditioner
+
+
+ /** We have J = [A C^T]
+ * [C 0]
+ *
+ * We want to approximate C A^(-1) C^T.
+ *
+ * Consider Lagrange vertex L that constrains the relative
+ * displacement between vertex N on the negative side of the fault
+ * and vertex P on the positive side of the fault.
+ *
+ * If we approximate A(-1) by 1/diag(A), then we can write
+ * C A^(-1) C^T for a 2-D case as
+ *
+ * [-R00 -R01 R00 R01][Ai_nx 0 0 0 ][-R00 -R10]
+ * [-R10 -R11 R10 R11][ Ai_ny 0 0 ][-R01 -R11]
+ * [ 0 Ai_px 0 ][ R00 R10]
+ * [ Ai_py][ R01 R11]
+ *
+ * where
+ *
+ * Ai_nx is the inverse of the diag(A) for DOF x of vertex N
+ * Ai_ny is the inverse of the diag(A) for DOF y of vertex N
+ * Ai_px is the inverse of the diag(A) for DOF x of vertex P
+ * Ai_py is the inverse of the diag(A) for DOF y of vertex P
+ *
+ * If Ai_nx == Ai_ny and Ai_px == Ai_py, then the result is
+ * diagonal. Otherwise, the off-diagonal terms will be nonzero,
+ * but we expect them to be small. Since we already approximate
+ * the inverse of A by the inverse of the diagonal, we drop the
+ * off-diagonal terms of C A^(-1) C^T:
+ *
+ * Term for DOF x of vertex L is:
+ * R00^2 (Ai_nx + Ai_px) + R01^2 (Ai_ny + Ai_py)
+ *
+ * Term for DOF y of vertex L is:
+ * R10^2 (Ai_nx + Ai_px) + R11^2 (Ai_ny + Ai_py)
+ *
+ */
+} // calcPreconditioner
+
+// ----------------------------------------------------------------------
// Adjust solution from solver with lumped Jacobian to match Lagrange
// multiplier constraints.
void
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.hh 2010-05-18 23:15:13 UTC (rev 16740)
@@ -154,6 +154,22 @@
const double t,
topology::SolutionFields* const fields);
+ /** Compute custom fault precoditioner using Schur complement.
+ *
+ * We have J = [A C^T]
+ * [C 0]
+ *
+ * We approximate C A^(-1) C^T.
+ *
+ * @param pc PETSc preconditioner structure.
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param fields Solution fields
+ */
+ virtual
+ void calcPreconditioner(PetscPC* const pc,
+ topology::Jacobian* const jacobian,
+ topology::SolutionFields* const fields);
+
/** Adjust solution from solver with lumped Jacobian to match Lagrange
* multiplier constraints.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-05-18 23:15:13 UTC (rev 16740)
@@ -25,6 +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 "spatialdata/spatialdb/spatialdbfwd.hh" // USES GravityField
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
@@ -233,6 +234,18 @@
const double t,
topology::SolutionFields* const fields);
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param pc PETSc preconditioner structure.
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param fields Solution fields
+ */
+ virtual
+ void calcPreconditioner(PetscPC* const pc,
+ topology::Jacobian* const jacobian,
+ topology::SolutionFields* const fields);
+
/** Update state variables as needed.
*
* @param t Current time
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -152,6 +152,17 @@
topology::SolutionFields* const fields) {
} // integrateJacobianAssembled
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::calcPreconditioner(
+ PetscPC* const pc,
+ topology::Jacobian* const jacobian,
+ topology::SolutionFields* const fields) {
+} // calcPreconditioner
+
// Update state variables as needed.
template<typename quadrature_type>
inline
Modified: short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -121,7 +121,7 @@
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- err = VTKViewer::writeHeader(_viewer);
+ err = VTKViewer::writeHeader(sieveMesh, _viewer);
CHECK_PETSC_ERROR(err);
//std::cout << "Wrote header for " << filename << std::endl;
err = VTKViewer::writeVertices(sieveMesh, _viewer);
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -54,6 +54,38 @@
} // deallocate
// ----------------------------------------------------------------------
+// Set flag for splitting fields.
+void
+pylith::problems::Formulation::splitFields(const bool flag)
+{ // splitFields
+ _splitFields = flag;
+} // splitFields
+
+// ----------------------------------------------------------------------
+// Get flag for splitting fields.
+bool
+pylith::problems::Formulation::splitFields(void) const
+{ // splitFields
+ return _splitFields;
+} // splitFields
+
+// ----------------------------------------------------------------------
+// Set flag for using custom fault preconditioner.
+void
+pylith::problems::Formulation::useCustomFaultPC(const bool flag)
+{ // useCustomFaultPC
+ _useCustomFaultPC = flag;
+} // useCustomFaultPC
+
+// ----------------------------------------------------------------------
+// Get flag indicating use of custom fault conditioner.
+bool
+pylith::problems::Formulation::useCustomFaultPC(void) const
+{ // useCustomFaultPC
+ return _useCustomFaultPC;
+} // useCustomFaultPC
+
+// ----------------------------------------------------------------------
// Return the fields
const pylith::topology::SolutionFields&
pylith::problems::Formulation::fields(void) const
@@ -310,6 +342,18 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
+
+#if 0
+ // Recalculate preconditioner.
+ if (_useCustomFaultPC) {
+ numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->calcPreconditioner(pc, _jacobian, _fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _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:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-05-18 23:15:13 UTC (rev 16740)
@@ -26,7 +26,7 @@
#include "pylith/feassemble/feassemblefwd.hh" // USES Integrator
#include "pylith/topology/topologyfwd.hh" // USES Mesh, Field, SolutionFields
-#include "pylith/utils/petscfwd.h" // USES PetscVec, PetscMat, PetscSNES
+#include "pylith/utils/petscfwd.h" // USES PetscVec, PetscMat
#include "pylith/utils/array.hh" // HASA std::vector
@@ -55,6 +55,30 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
+ /** Set flag for splitting fields.
+ *
+ * @param flag True if splitting fields, false otherwise.
+ */
+ void splitFields(const bool flag);
+
+ /** Get flag for splitting fields.
+ *
+ * @returns flag True if splitting fields, false otherwise.
+ */
+ bool splitFields(void) const;
+
+ /** Set flag for using custom fault preconditioner.
+ *
+ * @param flag True if splitting fields, false otherwise.
+ */
+ void useCustomFaultPC(const bool flag);
+
+ /** Get flag indicating use of custom fault conditioner.
+ *
+ * @returns flag True if using custom fault preconditioner, false otherwise.
+ */
+ bool useCustomFaultPC(void) const;
+
/** Get solution fields.
*
* @returns solution fields.
@@ -167,6 +191,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.
topology::Field<topology::Mesh>* _jacobianLumped; ///< Handle to lumped Jacobian of system.
topology::SolutionFields* _fields; ///< Handle to solution fields for system.
@@ -177,6 +202,8 @@
std::vector<IntegratorSubMesh*> _submeshIntegrators;
bool _isJacobianSymmetric; ///< Is system Jacobian symmetric?
+ bool _splitFields; ///< True if splitting fields.
+ bool _useCustomFaultPC; ///< True if using custom fault preconditioner.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -14,10 +14,21 @@
#include "Solver.hh" // implementation of class methods
+#include "Formulation.hh" // USES Formulation
+
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
#include "pylith/utils/EventLogger.hh" // USES EventLogger
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
#include <cassert> // USES assert()
+#define FIELD_SPLIT
+
+#if defined(FIELD_SPLIT)
+#include <petscmesh_solvers.hh> // USES constructFieldSplit()
+#endif
+
// ----------------------------------------------------------------------
// Constructor
pylith::problems::Solver::Solver(void) :
@@ -54,5 +65,62 @@
_formulation = formulation;
} // initialize
+// ----------------------------------------------------------------------
+// Setup preconditioner for preconditioning with split fields.
+void
+pylith::problems::Solver::_setupFieldSplit(PetscPC* const pc,
+ Formulation* const formulation,
+ const topology::SolutionFields& fields)
+{ // _setupFieldSplit
+ PetscErrorCode err = 0;
+
+ const topology::Field<topology::Mesh>& solution = fields.solution();
+ const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh =
+ fields.mesh().sieveMesh();
+ 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(),
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+ solution.section()),
+ solution.vector(), *pc);
+
+ if (formulation->useCustomFaultPC()) {
+ 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
+ MatStructure flag;
+ err = KSPGetOperators(ksps[num-1], &A,
+ PETSC_NULL, &flag); 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. This should really go in the Fault
+ // object that calculates the preconditioner since it is
+ // implementation dependent.
+ err = MatSeqAIJSetPreallocation(P, 1, PETSC_NULL); CHECK_PETSC_ERROR(err);
+ err = MatMPIAIJSetPreallocation(P, 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);
+ } // if
+
+#endif
+} // _setupFieldSplit
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2010-05-18 23:15:13 UTC (rev 16740)
@@ -25,6 +25,7 @@
#include "pylith/topology/topologyfwd.hh" // USES SolutionFields
#include "pylith/utils/utilsfwd.hh" // USES EventLogger
+#include "pylith/utils/petscfwd.h" // USES PetscPC
// Solver ---------------------------------------------------------
/** @brief Abstract C++ base class for using PETSc linear and
@@ -60,6 +61,20 @@
const topology::Jacobian& jacobian,
Formulation* const formulation);
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+ /** Setup preconditioner for preconditioning with split fields.8
+ *
+ * @param pc PETSc preconditioner.
+ * @param formulation Formulation of system of equations.
+ * @param fields Solution fields.
+ */
+ void
+ _setupFieldSplit(PetscPC* const pc,
+ Formulation* const formulation,
+ const topology::SolutionFields& fields);
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2010-05-18 23:15:13 UTC (rev 16740)
@@ -23,13 +23,6 @@
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
-//#define FIELD_SPLIT
-#define NEW_FAULT_PRECONDITIONER
-
-#if defined(FIELD_SPLIT)
-#include <petscmesh_solvers.hh> // USES constructFieldSplit()
-#endif
-
// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverLinear::SolverLinear(void) :
@@ -79,82 +72,11 @@
err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE); CHECK_PETSC_ERROR(err);
err = KSPSetFromOptions(_ksp); CHECK_PETSC_ERROR(err);
- const topology::Field<topology::Mesh>& residual = fields.get("residual");
-
- // Check for fibration
- if (residual.section()->getNumSpaces() > 0) {
- const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
- PC pc;
-
+ if (formulation->splitFields()) {
+ PetscPC pc = 0;
err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
- 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(residual.section(), sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", residual.section()), residual.vector(), pc);
-#if defined(NEW_FAULT_PRECONDITIONER)
- if (residual.section()->getNumSpaces() > sieveMesh->getDimension()) {
- KSP *ksps;
- Mat A, M;
- PetscInt num, m, n;
-
- err = PCFieldSplitGetSubKSP(pc, &num, &ksps); CHECK_PETSC_ERROR(err);
- // Put in PC matrix for fault
- MatStructure flag;
- err = KSPGetOperators(ksps[num-1], &A, PETSC_NULL, &flag); 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(), &M); CHECK_PETSC_ERROR(err);
- err = MatSetSizes(M, m, n, PETSC_DECIDE, PETSC_DECIDE); CHECK_PETSC_ERROR(err);
- err = MatSeqAIJSetPreallocation(M, 1, PETSC_NULL); CHECK_PETSC_ERROR(err);
- err = MatMPIAIJSetPreallocation(M, 1, PETSC_NULL, 0, PETSC_NULL); CHECK_PETSC_ERROR(err);
- err = MatSetFromOptions(M); CHECK_PETSC_ERROR(err);
- err = KSPSetOperators(ksps[num-1], A, M, flag); CHECK_PETSC_ERROR(err);
- // Create a mapping to indices for that space (might be in FS)
- err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
- } // if
-
- /** We have J = [A C^T]
- * [C 0]
- *
- * We want to approximate C A^(-1) C^T.
- *
- * Consider Lagrange vertex L that constrains the relative
- * displacement between vertex N on the negative side of the fault
- * and vertex P on the positive side of the fault.
- *
- * If we approximate A(-1) by 1/diag(A), then we can write
- * C A^(-1) C^T for a 2-D case as
- *
- * [-R00 -R01 R00 R01][Ai_nx 0 0 0 ][-R00 -R10]
- * [-R10 -R11 R10 R11][ Ai_ny 0 0 ][-R01 -R11]
- * [ 0 Ai_px 0 ][ R00 R10]
- * [ Ai_py][ R01 R11]
- *
- * where
- *
- * Ai_nx is the inverse of the diag(A) for DOF x of vertex N
- * Ai_ny is the inverse of the diag(A) for DOF y of vertex N
- * Ai_px is the inverse of the diag(A) for DOF x of vertex P
- * Ai_py is the inverse of the diag(A) for DOF y of vertex P
- *
- * If Ai_nx == Ai_ny and Ai_px == Ai_py, then the result is
- * diagonal. Otherwise, the off-diagonal terms will be nonzero,
- * but we expect them to be small. Since we already approximate
- * the inverse of A by the inverse of the diagonal, we drop the
- * off-diagonal terms of C A^(-1) C^T:
- *
- * Term for DOF x of vertex L is:
- * R00^2 (Ai_nx + Ai_px) + R01^2 (Ai_ny + Ai_py)
- *
- * Term for DOF y of vertex L is:
- * R10^2 (Ai_nx + Ai_px) + R11^2 (Ai_ny + Ai_py)
- *
- */
-
-#endif
-#endif
- }
+ _setupFieldSplit(&pc, formulation, fields);
+ } // if
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/utils/petscfwd.h
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/petscfwd.h 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/libsrc/utils/petscfwd.h 2010-05-18 23:15:13 UTC (rev 16740)
@@ -38,6 +38,9 @@
/// forward declaration for PETSc SNES
typedef struct _p_SNES* PetscSNES;
+/// forward declaration for PETSc PC
+typedef struct _p_PC* PetscPC;
+
/// forward declaration for PETSc PetscErrorCode
typedef int PetscErrorCode;
Modified: short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-05-18 23:15:13 UTC (rev 16740)
@@ -34,6 +34,31 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
+ /** Set flag for splitting fields.
+ *
+ * @param flag True if splitting fields, false otherwise.
+ */
+ void splitFields(const bool flag);
+
+ /** Get flag for splitting fields.
+ *
+ * @returns flag True if splitting fields, false otherwise.
+ */
+ bool splitFields(void) const;
+
+ /** Set flag for using custom fault preconditioner.
+ *
+ * @param flag True if splitting fields, false otherwise.
+ */
+ void useCustomFaultPC(const bool flag);
+
+ /** Get flag indicating use of custom fault conditioner.
+ *
+ * @returns flag True if using custom fault preconditioner,
+ * false otherwise.
+ */
+ bool useCustomFaultPC(void) const;
+
/** Get solution fields.
*
* @returns solution fields.
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-05-18 23:03:03 UTC (rev 16739)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2010-05-18 23:15:13 UTC (rev 16740)
@@ -58,6 +58,7 @@
## \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 view_jacobian Flag to output Jacobian matrix when it is
## reformed.
@@ -73,10 +74,13 @@
matrixType = pyre.inventory.str("matrix_type", default="unknown")
matrixType.meta['tip'] = "Type of PETSc sparse matrix."
- splitFields = pyre.inventory.bool("split_fields", default=False)
- splitFields.meta['tip'] = "Split solution fields into displacements "\
+ useSplitFields = pyre.inventory.bool("split_fields", default=False)
+ 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."
+
viewJacobian = pyre.inventory.bool("view_jacobian", default=False)
viewJacobian.meta['tip'] = "Write Jacobian matrix to binary file."
@@ -279,7 +283,6 @@
"""
PetscComponent._configure(self)
self.matrixType = self.inventory.matrixType
- self.splitFields = self.inventory.splitFields
self.timeStep = self.inventory.timeStep
self.solver = self.inventory.solver
self.output = self.inventory.output
@@ -289,6 +292,14 @@
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'."
+ self.inventory.useSplitFields = True
+
+ ModuleFormulation.splitFields(self, self.inventory.useSplitFields)
+ ModuleFormulation.useCustomFaultPC(self, self.inventory.useCustomFaultPC)
return
More information about the CIG-COMMITS
mailing list