[cig-commits] r18125 - in short/3D/PyLith/trunk: . libsrc/meshio libsrc/problems libsrc/topology libsrc/utils
knepley at geodynamics.org
knepley at geodynamics.org
Tue Mar 22 15:39:08 PDT 2011
Author: knepley
Date: 2011-03-22 15:39:08 -0700 (Tue, 22 Mar 2011)
New Revision: 18125
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIOSieve.cc
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
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc
short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh
short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh
short/3D/PyLith/trunk/libsrc/utils/sievetypes.hh
Log:
Another try on nonlinear preconditioning
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/TODO 2011-03-22 22:39:08 UTC (rev 18125)
@@ -17,6 +17,10 @@
each processor having lots of gaps in the range of cells and
vertices it has?
+ - Can Fusion be leaking memory?
+
+ - How do we separate parallel stratification?
+
IntSections - Matt improved distribution to use allow IntSections
over vertices or cells rather than vertices + cells. This cuts
differences down to 75%.
@@ -30,12 +34,19 @@
set the KSP operator. The _precondMatrix in SolverNonlinear::solve
matches the one used in the linear solve.
+ - Brad verifies that _precondMatrix at line SolverNonlinear.cc:180 is correct. Check that
+ the same Mat pointer is used in ML.
+
+ - FOUND: This is caused by the KSPSetOperators() in SNESSolve_LS()
+
* Unform global refinement
(2) The overlap is not updated properly for hex8 refinement. Running
examples/twocells/twohex8/axialdisp.cfg with 2 processors and
uniform refinement is a good test case.
+ - Checking addition of MPI type
+
* Output to HDF5 files. [BRAD and MATT]
(1) The Sieve section to PETSc Vec is not done correctly for
@@ -47,7 +58,10 @@
(3) Add creation of .xmf metadata file for ParaView/Visit. [BRAD]
+ - Need a flag to use s->GetIndicesRaw() in scatter creation
+ - examples/bar_shearwave/tri3
+
* Fault preconditioner [BRAD]
FaultCohesiveLagrange::calcPreconditioner() [need unit test]
Modified: short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -18,7 +18,7 @@
#include <portinfo>
-#include <petscmesh_viewers.hh> // USES VTKViewer
+#include <petscdmmesh_viewers.hh> // USES VTKViewer
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOSieve.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOSieve.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOSieve.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -20,7 +20,7 @@
#include "MeshIOSieve.hh" // implementation of class methods
-#include "petscmesh.hh"
+#include <petscdmmesh.hh>
#include "MeshBuilder.hh" // USES MeshBuilder
#include "pylith/topology/Mesh.hh" // USES Mesh
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -23,6 +23,7 @@
#include "Formulation.hh" // USES Formulation
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/utils/EventLogger.hh" // USES EventLogger
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
@@ -32,18 +33,42 @@
#define FIELD_SPLIT
#if defined(FIELD_SPLIT)
-#include <petscmesh_solvers.hh> // USES constructFieldSplit()
+#include <petscdmmesh_solvers.hh> // USES constructFieldSplit()
#endif
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+EXTERN_C_BEGIN
+PetscErrorCode MyMatGetSubMatrix(Mat mat, IS isrow, IS iscol, MatReuse reuse, Mat *newmat) {
+ FaultPreconCtx *ctx;
+ IS faultIS;
+ PetscBool isFaultRow, isFaultCol;
+ PetscErrorCode ierr;
+
+ ierr = MatShellGetContext(mat, (void **) &ctx);CHKERRQ(ierr);
+ ierr = PCFieldSplitGetIS(ctx->pc, ctx->faultFieldName, &faultIS);CHKERRQ(ierr);
+ ierr = ISEqual(isrow, faultIS, &isFaultRow);CHKERRQ(ierr);
+ ierr = ISEqual(iscol, faultIS, &isFaultCol);CHKERRQ(ierr);
+ if (isFaultRow && isFaultCol) {
+ if (reuse == MAT_INITIAL_MATRIX) {
+ ierr = PetscObjectReference((PetscObject) ctx->faultA);CHKERRQ(ierr);
+ *newmat = ctx->faultA;
+ }
+ } else {
+ ierr = MatGetSubMatrix(ctx->A, isrow, iscol, reuse, newmat);CHKERRQ(ierr);
+ }
+}
+EXTERN_C_END
+
// ----------------------------------------------------------------------
// Constructor
pylith::problems::Solver::Solver(void) :
_formulation(0),
- _logger(0)
+ _logger(0),
+ _precondMatrix(0),
+ _jacobianPre(0)
{ // constructor
} // constructor
@@ -61,6 +86,14 @@
{ // deallocate
_formulation = 0; // Handle only, do not manage memory.
delete _logger; _logger = 0;
+ if (0 != _jacobianPre) {
+ PetscErrorCode err = MatDestroy(_jacobianPre); _jacobianPre = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
+ if (0 != _precondMatrix) {
+ PetscErrorCode err = MatDestroy(_precondMatrix); _precondMatrix = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
} // deallocate
// ----------------------------------------------------------------------
@@ -72,6 +105,26 @@
{ // initialize
assert(0 != formulation);
+ PetscMat jacobianMat = jacobian.matrix();
+ _jacobianPre = jacobianMat;
+ // Make global preconditioner matrix
+ const ALE::Obj<RealSection>& solutionSection = fields.solution().section();
+ assert(!solutionSection.isNull());
+ const ALE::Obj<SieveMesh>& sieveMesh = fields.solution().mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ if (solutionSection->getNumSpaces() > sieveMesh->getDimension() && 0 != _precondMatrix) {
+ PetscInt M, N, m, n;
+ PetscErrorCode err;
+
+ err = MatGetSize(jacobianMat, &M, &N);CHECK_PETSC_ERROR(err);
+ err = MatGetLocalSize(jacobianMat, &m, &n);CHECK_PETSC_ERROR(err);
+ err = MatCreateShell(fields.mesh().comm(), m, n, M, N, &_ctx, &_jacobianPre);CHECK_PETSC_ERROR(err);
+ err = MatShellSetOperation(_jacobianPre, MATOP_GET_SUBMATRIX, (void (*)(void)) MyMatGetSubMatrix);
+ _ctx.A = jacobianMat;
+ _ctx.faultFieldName = "3";
+ _ctx.faultA = _precondMatrix;
+ }
+
_formulation = formulation;
} // initialize
Modified: short/3D/PyLith/trunk/libsrc/problems/Solver.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/Solver.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -33,6 +33,14 @@
#include "pylith/utils/utilsfwd.hh" // USES EventLogger
#include "pylith/utils/petscfwd.h" // USES PetscMat
+typedef struct {
+ PetscPC pc;
+ PetscMat A;
+ const char *faultFieldName;
+ PetscMat faultA;
+} FaultPreconCtx;
+
+
// Solver ---------------------------------------------------------
/** @brief Abstract C++ base class for using PETSc linear and
* nonlinear solvers.
@@ -88,7 +96,9 @@
Formulation* _formulation; ///< Handle to formulation for system of eqns.
utils::EventLogger* _logger; ///< Event logger.
-
+ PetscMat _precondMatrix; ///< Preconditioning matrix for Lagrange constraints.
+ PetscMat _jacobianPre; ///< global fault preconditioning matrix
+ FaultPreconCtx _ctx; ///< Context for global fault preconditioning matrix
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -36,8 +36,7 @@
// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverLinear::SolverLinear(void) :
- _ksp(0),
- _precondMatrix(0)
+ _ksp(0)
{ // constructor
} // constructor
@@ -59,11 +58,6 @@
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
// ----------------------------------------------------------------------
@@ -91,6 +85,7 @@
if (formulation->splitFields()) {
PetscPC pc = 0;
err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
+ _ctx.pc = pc;
_setupFieldSplit(&pc, &_precondMatrix, formulation, fields);
} // if
} // initialize
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLinear.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -85,7 +85,6 @@
private :
PetscKSP _ksp; ///< PETSc KSP linear solver.
- PetscMat _precondMatrix; ///< Preconditioning matrix for Lagrange constraints.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -64,8 +64,7 @@
// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverNonlinear::SolverNonlinear(void) :
- _snes(0),
- _precondMatrix(0)
+ _snes(0)
{ // constructor
} // constructor
@@ -87,10 +86,6 @@
PetscErrorCode err = SNESDestroy(_snes); _snes = 0;
CHECK_PETSC_ERROR(err);
} // if
- if (0 != _precondMatrix) {
- PetscErrorCode err = MatDestroy(_precondMatrix); _precondMatrix = 0;
- CHECK_PETSC_ERROR(err);
- } // if
} // deallocate
// ----------------------------------------------------------------------
@@ -119,20 +114,18 @@
(void*) formulation);
CHECK_PETSC_ERROR(err);
- PetscMat jacobianMat = jacobian.matrix();
- err = SNESSetJacobian(_snes, jacobianMat, jacobianMat, reformJacobian,
- (void*) formulation);
+ err = SNESSetJacobian(_snes, jacobian.matrix(), _jacobianPre, reformJacobian, (void*) formulation);
CHECK_PETSC_ERROR(err);
err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
- err = SNESLineSearchSet(_snes, lineSearch,
- (void*) formulation); CHECK_PETSC_ERROR(err);
+ err = SNESLineSearchSet(_snes, lineSearch, (void*) formulation); CHECK_PETSC_ERROR(err);
if (formulation->splitFields()) {
PetscKSP ksp = 0;
PetscPC pc = 0;
err = SNESGetKSP(_snes, &ksp); CHECK_PETSC_ERROR(err);
err = KSPGetPC(ksp, &pc); CHECK_PETSC_ERROR(err);
+ _ctx.pc = pc;
_setupFieldSplit(&pc, &_precondMatrix, formulation, fields);
} // if
} // initialize
@@ -147,11 +140,12 @@
{ // solve
assert(0 != solution);
+#if 0
// 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());
+ assert(!sieveMesh.isNull());
if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
0 != _precondMatrix) {
PetscKSP ksp = 0;
@@ -181,6 +175,7 @@
flag); CHECK_PETSC_ERROR(err);
err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
} // if
+#endif
const int solveEvent = _logger->eventId("SoNl solve");
const int scatterEvent = _logger->eventId("SoNl scatter");
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -151,7 +151,6 @@
private :
PetscSNES _snes; ///< PETSc SNES nonlinear solver.
- PetscMat _precondMatrix; ///< Preconditioning matrix for Lagrange constraints.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -763,8 +763,7 @@
assert(!order.isNull());
// Create scatter
- err = MeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order,
- &sinfo.scatter);
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, &sinfo.scatter);
CHECK_PETSC_ERROR(err);
// Create scatterVec
@@ -837,8 +836,7 @@
assert(!order.isNull());
// Create scatter
- err = MeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order,
- &sinfo.scatter);
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, &sinfo.scatter);
CHECK_PETSC_ERROR(err);
// Create scatterVec
@@ -913,8 +911,7 @@
assert(!order.isNull());
// Create scatter
- err = MeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order,
- &sinfo.scatter);
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, &sinfo.scatter);
CHECK_PETSC_ERROR(err);
// Create scatterVec
@@ -988,8 +985,7 @@
order->view("GLOBAL ORDER");
// Create scatter
- err = MeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order,
- &sinfo.scatter);
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, &sinfo.scatter);
CHECK_PETSC_ERROR(err);
// Create scatterVec
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -32,7 +32,7 @@
#include "pylith/utils/arrayfwd.hh" // HASA int_array
#include "pylith/utils/petscfwd.h" // HASA PetscVec
-#include <petscmesh.hh>
+#include <petscdmmesh.hh>
#include <map> // USES std::map
#include <string> // USES std::string
Modified: short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/topology/Jacobian.cc 2011-03-22 22:39:08 UTC (rev 18125)
@@ -47,8 +47,8 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, fieldSection->getName(), fieldSection)->view("Global Order");
sieveMesh->setDebug(3);
#endif
- PetscErrorCode err = MeshCreateMatrix(sieveMesh, fieldSection,
- matrixType, &_matrix, blockFlag);
+ PetscErrorCode err = DMMeshCreateMatrix(sieveMesh, fieldSection,
+ matrixType, &_matrix, blockFlag);
CHECK_PETSC_ERROR_MSG(err, "Could not create PETSc sparse matrix "
"associated with system Jacobian.");
logger.stagePop();
@@ -73,8 +73,8 @@
// dimension, otherwise use a block size of 1.
const int blockFlag = (blockOkay) ? -1 : 1;
- PetscErrorCode err = MeshCreateMatrix(sieveMesh, fieldSection,
- matrixType, &_matrix, blockFlag);
+ PetscErrorCode err = DMMeshCreateMatrix(sieveMesh, fieldSection,
+ matrixType, &_matrix, blockFlag);
CHECK_PETSC_ERROR_MSG(err, "Could not create PETSc sparse matrix "
"associated with subsystem Jacobian.");
logger.stagePop();
Modified: short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/topology/Jacobian.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -31,7 +31,7 @@
#include "pylith/utils/petscfwd.h" // HOLDSA PetscMat
-#include <petscmesh.hh> // USES MPI_Comm
+#include <petscdmmesh.hh> // USES MPI_Comm
#include <string> // USES std::string
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -30,7 +30,7 @@
// Include directives ---------------------------------------------------
#include "topologyfwd.hh" // forward declarations
-#include <petscmesh.hh> // HASA ALE::IMesh
+#include <petscdmmesh.hh> // HASA ALE::IMesh
// MeshOrder ------------------------------------------------------------
/// Object for managing order of mesh entities.
Modified: short/3D/PyLith/trunk/libsrc/utils/sievetypes.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/sievetypes.hh 2011-03-22 21:52:38 UTC (rev 18124)
+++ short/3D/PyLith/trunk/libsrc/utils/sievetypes.hh 2011-03-22 22:39:08 UTC (rev 18125)
@@ -25,7 +25,7 @@
#if !defined(pylith_utils_sievetypes_hh)
#define pylith_utils_sievetypes_hh
-#include <petscmesh.hh> // PETSc Mesh
+#include <petscdmmesh.hh> // PETSc Mesh
namespace pylith {
More information about the CIG-COMMITS
mailing list