[cig-commits] r22112 - in short/3D/PyLith/trunk: examples/3d/tet4 libsrc/pylith/faults libsrc/pylith/problems modulesrc/problems pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Mon May 20 12:14:13 PDT 2013
Author: brad
Date: 2013-05-20 12:14:13 -0700 (Mon, 20 May 2013)
New Revision: 22112
Modified:
short/3D/PyLith/trunk/examples/3d/tet4/README
short/3D/PyLith/trunk/examples/3d/tet4/step02.cfg
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh
short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Simplify PC settings in case with not fault.
Modified: short/3D/PyLith/trunk/examples/3d/tet4/README
===================================================================
--- short/3D/PyLith/trunk/examples/3d/tet4/README 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/examples/3d/tet4/README 2013-05-20 19:14:13 UTC (rev 22112)
@@ -41,9 +41,9 @@
particular example. The example problems are briefly described below:
step01: Dirichlet BC (static) causing shear.
-step02: step01 with mesh refinement and field split preconditioning.
+step02: step01 with mesh refinement and algebraic multigrid preconditioner.
step03: DirichletBC with fault slip (static) Slip on a vertical fault
-step04: step03 with mesh refinement and field split preconditioning.
+step04: step03 with mesh refinement and algebraic multigrid preconditioner.
----------------------------------------
figures directory
Modified: short/3D/PyLith/trunk/examples/3d/tet4/step02.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/tet4/step02.cfg 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/examples/3d/tet4/step02.cfg 2013-05-20 19:14:13 UTC (rev 22112)
@@ -94,16 +94,17 @@
[pylithapp.petsc]
# Field split preconditioning. We use an algebraic multigrid
-# preconditioner from Trilinos/ML built via PETSc on each component of
-# the displacement field separately. In this way, we split the
-# fields. This method requires the 'aij' matrix type.
+# preconditioner from Trilinos/ML. This method requires the 'aij'
+# matrix type.
[pylithapp.timedependent.formulation]
-split_fields = True
+#split_fields = True
matrix_type = aij
[pylithapp.petsc]
-fs_pc_type = fieldsplit
-fs_pc_fieldsplit_real_diagonal = True
-fs_pc_fieldsplit_type = multiplicative
-fs_fieldsplit_0_pc_type = ml
-fs_fieldsplit_0_ksp_type = preonly
+pc_type = ml
+# Old settings
+#fs_pc_type = fieldsplit
+#fs_pc_fieldsplit_real_diagonal = True
+#fs_pc_fieldsplit_type = multiplicative
+#fs_fieldsplit_0_pc_type = ml
+#fs_fieldsplit_0_ksp_type = preonly
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-20 19:14:13 UTC (rev 22112)
@@ -137,9 +137,10 @@
PetscErrorCode err;
err = PetscSectionGetNumFields(fieldSection, &numFields);PYLITH_CHECK_ERROR(err);
// TODO: Does this make sense?
- if (!numFields)
+ if (!numFields) {
PYLITH_METHOD_END;
- assert(numFields == 2);
+ } // if
+ assert(2 == numFields);
err = PetscSectionGetFieldComponents(fieldSection, 0, &numComp);PYLITH_CHECK_ERROR(err);assert(numComp == spaceDim);
err = PetscSectionGetFieldComponents(fieldSection, 1, &numComp);PYLITH_CHECK_ERROR(err);assert(numComp == spaceDim);
@@ -148,8 +149,7 @@
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
PetscInt dof;
- err = PetscSectionGetDof(fieldSection, v_lagrange, &dof);PYLITH_CHECK_ERROR(err);
- assert(spaceDim == dof);
+ err = PetscSectionGetDof(fieldSection, v_lagrange, &dof);PYLITH_CHECK_ERROR(err);assert(spaceDim == dof);
err = PetscSectionSetFieldDof(fieldSection, v_lagrange, 1, dof);PYLITH_CHECK_ERROR(err);
} // for
err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);PYLITH_CHECK_ERROR(err);
@@ -159,7 +159,7 @@
if (!dof) {
err = PetscSectionSetFieldDof(fieldSection, p, 0, spaceDim);PYLITH_CHECK_ERROR(err);
} // if
- }
+ } // for
PYLITH_METHOD_END;
} // splitField
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2013-05-20 19:14:13 UTC (rev 22112)
@@ -40,8 +40,7 @@
_fields(0),
_customConstraintPCMat(0),
_isJacobianSymmetric(false),
- _splitFields(false),
- _splitFieldComponents(false)
+ _splitFields(false)
{ // constructor
} // constructor
@@ -93,22 +92,6 @@
} // splitFields
// ----------------------------------------------------------------------
-// Set flag for splitting field components.
-void
-pylith::problems::Formulation::splitFieldComponents(const bool flag)
-{ // splitFieldComponents
- _splitFieldComponents = flag;
-} // splitFieldComponents
-
-// ----------------------------------------------------------------------
-// Get flag for splitting field components.
-bool
-pylith::problems::Formulation::splitFieldComponents(void) const
-{ // splitFieldComponents
- return _splitFieldComponents;
-} // splitFieldComponents
-
-// ----------------------------------------------------------------------
// Set flag for using custom preconditioner for Lagrange constraints.
void
pylith::problems::Formulation::useCustomConstraintPC(const bool flag)
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh 2013-05-20 19:14:13 UTC (rev 22112)
@@ -67,18 +67,6 @@
*/
bool splitFields(void) const;
- /** Set flag for splitting field components.
- *
- * @param flag True if splitting fields comonents, false otherwise.
- */
- void splitFieldComponents(const bool flag);
-
- /** Get flag for splitting field components.
- *
- * @returns flag True if splitting field components, false otherwise.
- */
- bool splitFieldComponents(void) const;
-
/** Set flag for using custom preconditioner for Lagrange constraints.
*
* @param flag True if using custom fault preconditioner, false otherwise.
@@ -196,10 +184,8 @@
bool _isJacobianSymmetric; ///< Is system Jacobian symmetric?
bool _splitFields; ///< True if splitting fields.
- bool _splitFieldComponents; ///< True if splitting field components
- /// True if using custom preconditioner for Lagrange constraints.
- bool _useCustomConstraintPC;
+ bool _useCustomConstraintPC; ///< True if using custom preconditioner for Lagrange constraints.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2013-05-20 19:14:13 UTC (rev 22112)
@@ -166,19 +166,20 @@
PetscInt spaceDim, numFields;
PetscErrorCode err;
- const bool separateComponents = formulation->splitFieldComponents();
-
err = PetscObjectGetComm((PetscObject) dmMesh, &comm);PYLITH_CHECK_ERROR(err);
err = DMPlexGetDimension(dmMesh, &spaceDim);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetNumFields(solutionSection, &numFields);PYLITH_CHECK_ERROR(err);
err = PCSetDM(*pc, dmMesh);PYLITH_CHECK_ERROR(err);
err = PCSetType(*pc, PCFIELDSPLIT);PYLITH_CHECK_ERROR(err);
err = PCSetOptionsPrefix(*pc, "fs_");PYLITH_CHECK_ERROR(err);
err = PCSetFromOptions(*pc);PYLITH_CHECK_ERROR(err);
- if (formulation->splitFields() && formulation->useCustomConstraintPC() &&
- numFields > 1) {
+ if (formulation->splitFields() && formulation->useCustomConstraintPC()) {
+ err = PetscSectionGetNumFields(solutionSection, &numFields);PYLITH_CHECK_ERROR(err);
+ if (numFields < 2) {
+ throw std::logic_error("Cannot setup solution field for split fields. Solution field has only one field.");
+ } // if
+
// We have split fields with a custom constraint preconditioner
// and constraints exist.
@@ -206,98 +207,84 @@
_ctx.pc = *pc;
_ctx.A = jacobian.matrix();
- switch (spaceDim) {
- case 1 :
- _ctx.faultFieldName = "1";
- break;
- case 2 :
- _ctx.faultFieldName = (separateComponents) ? "2" : "1";
- break;
- case 3 :
- _ctx.faultFieldName = (separateComponents) ? "3" : "1";
- break;
- default:
- assert(0);
- throw std::logic_error("Unknown space dimension in "
- "Problems::_setupFieldSplit().");
- } // switch
+ _ctx.faultFieldName = "1";
_ctx.faultA = _jacobianPCFault;
} // if
+
+ // Create rigid body null space.
+ MatNullSpace nullsp = NULL;
+ PetscObject field = NULL;
+ PetscSection coordinateSection = NULL;
+ PetscVec coordinateVec = NULL;
+ PetscInt dim = spaceDim, vStart, vEnd;
+ PetscVec mode[6];
+
+ err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);PYLITH_CHECK_ERROR(err);
+ err = DMPlexGetCoordinateSection(dmMesh, &coordinateSection);PYLITH_CHECK_ERROR(err);assert(coordinateSection);
+ err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);PYLITH_CHECK_ERROR(err);assert(coordinateVec);
+ if (dim > 1) {
+ const int m = (dim * (dim + 1)) / 2;
+ for(int i = 0; i < m; ++i) {
+ err = VecDuplicate(solutionGlobalVec, &mode[i]);PYLITH_CHECK_ERROR(err);
+ } // for
+ // :KLUDGE: Assume P1
+ for(int d = 0; d < dim; ++d) {
+ PetscScalar values[3] = {0.0, 0.0, 0.0};
+ values[d] = 1.0;
- if (separateComponents) {
- throw std::logic_error("Separate components is no longer supported");
- } else {
- // Create rigid body null space.
- MatNullSpace nullsp = NULL;
- PetscObject field = NULL;
- PetscSection coordinateSection;
- PetscVec coordinateVec;
- PetscInt dim = spaceDim, vStart, vEnd;
- PetscVec mode[6];
-
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);PYLITH_CHECK_ERROR(err);
- err = DMPlexGetCoordinateSection(dmMesh, &coordinateSection);PYLITH_CHECK_ERROR(err);assert(coordinateSection);
- err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);PYLITH_CHECK_ERROR(err);assert(coordinateVec);
- if (dim > 1) {
- const int m = (dim * (dim + 1)) / 2;
- for(int i = 0; i < m; ++i) {
- err = VecDuplicate(solutionGlobalVec, &mode[i]);PYLITH_CHECK_ERROR(err);
+ err = VecSet(solutionVec, 0.0);PYLITH_CHECK_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = DMPlexVecSetClosure(dmMesh, solutionSection, solutionVec, v, values, INSERT_VALUES);PYLITH_CHECK_ERROR(err);
} // for
- // :KLUDGE: Assume P1
- for(int d = 0; d < dim; ++d) {
- PetscScalar values[3] = {0.0, 0.0, 0.0};
+ err = DMLocalToGlobalBegin(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ err = DMLocalToGlobalEnd(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ } // for
+ for(int d = dim; d < m; ++d) {
+ PetscInt k = (dim > 2) ? d - dim : d;
+
+ err = VecSet(solutionVec, 0.0);PYLITH_CHECK_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscScalar values[3] = {0.0, 0.0, 0.0};
+ PetscScalar *coords = NULL;
- values[d] = 1.0;
- err = VecSet(solutionVec, 0.0);PYLITH_CHECK_ERROR(err);
- for(PetscInt v = vStart; v < vEnd; ++v) {
- err = DMPlexVecSetClosure(dmMesh, solutionSection, solutionVec, v, values, INSERT_VALUES);PYLITH_CHECK_ERROR(err);
- } // for
- err = DMLocalToGlobalBegin(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
- err = DMLocalToGlobalEnd(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ err = DMPlexVecGetClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);PYLITH_CHECK_ERROR(err);
+ for(int i = 0; i < dim; ++i) {
+ for(int j = 0; j < dim; ++j) {
+ values[j] += _epsilon(i, j, k)*coords[i];
+ } // for
+ } // for
+ err = DMPlexVecRestoreClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);PYLITH_CHECK_ERROR(err);
+ err = DMPlexVecSetClosure(dmMesh, solutionSection, solutionVec, v, values, INSERT_VALUES);PYLITH_CHECK_ERROR(err);
} // for
- for(int d = dim; d < m; ++d) {
- PetscInt k = (dim > 2) ? d - dim : d;
-
- err = VecSet(solutionVec, 0.0);PYLITH_CHECK_ERROR(err);
- for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscScalar values[3] = {0.0, 0.0, 0.0};
- PetscScalar *coords = NULL;
-
- err = DMPlexVecGetClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);PYLITH_CHECK_ERROR(err);
- for(int i = 0; i < dim; ++i) {
- for(int j = 0; j < dim; ++j) {
- values[j] += _epsilon(i, j, k)*coords[i];
- } // for
- } // for
- err = DMPlexVecRestoreClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);PYLITH_CHECK_ERROR(err);
- err = DMPlexVecSetClosure(dmMesh, solutionSection, solutionVec, v, values, INSERT_VALUES);PYLITH_CHECK_ERROR(err);
- }
- err = DMLocalToGlobalBegin(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
- err = DMLocalToGlobalEnd(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ err = DMLocalToGlobalBegin(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ err = DMLocalToGlobalEnd(dmMesh, solutionVec, INSERT_VALUES, mode[d]);PYLITH_CHECK_ERROR(err);
+ } // for
+ for(int i = 0; i < dim; ++i) {
+ err = VecNormalize(mode[i], NULL);PYLITH_CHECK_ERROR(err);
+ } // for
+ // Orthonormalize system
+ for(int i = dim; i < m; ++i) {
+ PetscScalar dots[6];
+
+ err = VecMDot(mode[i], i, mode, dots);PYLITH_CHECK_ERROR(err);
+ for(int j = 0; j < i; ++j) {
+ dots[j] *= -1.0;
} // for
- for(int i = 0; i < dim; ++i) {
- err = VecNormalize(mode[i], NULL);PYLITH_CHECK_ERROR(err);
- } // for
- /* Orthonormalize system */
- for(int i = dim; i < m; ++i) {
- PetscScalar dots[6];
+ err = VecMAXPY(mode[i], i, dots, mode);PYLITH_CHECK_ERROR(err);
+ err = VecNormalize(mode[i], NULL);PYLITH_CHECK_ERROR(err);
+ } // for
+ err = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, &nullsp);PYLITH_CHECK_ERROR(err);
+ for(int i = 0; i< m; ++i) {
+ err = VecDestroy(&mode[i]);PYLITH_CHECK_ERROR(err);
+ } // for
+ } // if
+
+ err = DMSetNumFields(dmMesh, 2);PYLITH_CHECK_ERROR(err);
+ err = DMGetField(dmMesh, 0, &field);PYLITH_CHECK_ERROR(err);
+ err = PetscObjectCompose(field, "nearnullspace", (PetscObject) nullsp);PYLITH_CHECK_ERROR(err);
+ err = MatNullSpaceDestroy(&nullsp);PYLITH_CHECK_ERROR(err);
+ err = PetscObjectCompose(field, "pmat", (PetscObject) _jacobianPCFault);PYLITH_CHECK_ERROR(err);
- err = VecMDot(mode[i], i, mode, dots);PYLITH_CHECK_ERROR(err);
- for(int j = 0; j < i; ++j) dots[j] *= -1.0;
- err = VecMAXPY(mode[i], i, dots, mode);PYLITH_CHECK_ERROR(err);
- err = VecNormalize(mode[i], NULL);PYLITH_CHECK_ERROR(err);
- } // for
- err = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, &nullsp);PYLITH_CHECK_ERROR(err);
- for(int i = 0; i< m; ++i) {err = VecDestroy(&mode[i]);PYLITH_CHECK_ERROR(err);}
- } // if
-
- err = DMSetNumFields(dmMesh, 2);PYLITH_CHECK_ERROR(err);
- err = DMGetField(dmMesh, 0, &field);PYLITH_CHECK_ERROR(err);
- err = PetscObjectCompose(field, "nearnullspace", (PetscObject) nullsp);PYLITH_CHECK_ERROR(err);
- err = MatNullSpaceDestroy(&nullsp);PYLITH_CHECK_ERROR(err);
- err = PetscObjectCompose(field, "pmat", (PetscObject) _jacobianPCFault);PYLITH_CHECK_ERROR(err);
- } // if/else
-
PYLITH_METHOD_END;
} // _setupFieldSplit
Modified: short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2013-05-20 19:14:13 UTC (rev 22112)
@@ -52,18 +52,6 @@
*/
bool splitFields(void) const;
- /** Set flag for splitting field components.
- *
- * @param flag True if splitting fields comonents, false otherwise.
- */
- void splitFieldComponents(const bool flag);
-
- /** Get flag for splitting field components.
- *
- * @returns flag True if splitting field components, false otherwise.
- */
- bool splitFieldComponents(void) const;
-
/** Set flag for using custom Lagrange constraint preconditioner.
*
* @param flag True if using custom constraint precondition,
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2013-05-20 17:47:47 UTC (rev 22111)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2013-05-20 19:14:13 UTC (rev 22112)
@@ -74,13 +74,9 @@
## \b Properties
## @li \b use_cuda Enable use of CUDA for finite-element integrations.
## @li \b matrix_type Type of PETSc sparse matrix.
- ## @li \b split_fields Split solution fields into displacements
- ## and Lagrange constraints.
- ## @li \b split_field_components Split fields into components.
- ## @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.
+ ## @li \b split_fields Split solution fields into displacements and Lagrange constraints.
+ ## @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.
##
## \b Facilities
## @li \b time_step Time step size manager.
@@ -103,9 +99,6 @@
useSplitFields.meta['tip'] = "Split solution fields into displacements "\
"and Lagrange multipliers for separate preconditioning."
- useSplitFieldComponents = pyre.inventory.bool("split_field_components", default=False)
- useSplitFieldComponents.meta['tip'] = "Split solution fields into components for separate preconditioning."
-
useCustomConstraintPC = pyre.inventory.bool("use_custom_constraint_pc",
default=False)
useCustomConstraintPC.meta['tip'] = "Use custom preconditioner for " \
@@ -338,9 +331,7 @@
self.inventory.useSplitFields = True
ModuleFormulation.splitFields(self, self.inventory.useSplitFields)
- ModuleFormulation.splitFieldComponents(self, self.inventory.useSplitFieldComponents)
- ModuleFormulation.useCustomConstraintPC(self,
- self.inventory.useCustomConstraintPC)
+ ModuleFormulation.useCustomConstraintPC(self, self.inventory.useCustomConstraintPC)
return
@@ -535,7 +526,6 @@
solution.vectorFieldType(solution.VECTOR)
solution.scale(lengthScale.value)
if self.splitFields():
- solution.splitDefault()
for integrator in self.integrators:
integrator.splitField(solution)
for constraint in self.constraints:
More information about the CIG-COMMITS
mailing list