[cig-commits] r21588 - short/3D/PyLith/trunk/libsrc/pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Tue Mar 19 19:13:07 PDT 2013
Author: brad
Date: 2013-03-19 19:13:06 -0700 (Tue, 19 Mar 2013)
New Revision: 21588
Modified:
short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.hh
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLinear.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
Log:
Code cleanup in libsrc/problems.
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2013-03-20 01:37:35 UTC (rev 21587)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2013-03-20 02:13:06 UTC (rev 21588)
@@ -34,22 +34,22 @@
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "MyMatGetSubMatrix"
-PetscErrorCode MyMatGetSubMatrix(Mat mat, IS isrow, IS iscol, MatReuse reuse, Mat *newmat) {
- FaultPreconCtx *ctx;
- IS faultIS;
- PetscBool isFaultRow, isFaultCol;
+PetscErrorCode MyMatGetSubMatrix(PetscMat mat,
+ PetscIS isrow,
+ PetscIS iscol,
+ MatReuse reuse,
+ PetscMat *newmat) {
+ FaultPreconCtx *ctx = NULL;
+ PetscIS faultIS = NULL;
+ PetscBool isFaultRow, isFaultCol;
PetscErrorCode ierr = 0;
PetscFunctionBegin;
ierr = MatShellGetContext(mat, (void **) &ctx);CHKERRQ(ierr);
- ierr = PCFieldSplitGetIS(ctx->pc, ctx->faultFieldName, &faultIS);CHKERRQ(ierr);
- assert(faultIS);
+ ierr = PCFieldSplitGetIS(ctx->pc, ctx->faultFieldName, &faultIS);CHKERRQ(ierr);assert(faultIS);
ierr = ISEqual(isrow, faultIS, &isFaultRow);CHKERRQ(ierr);
ierr = ISEqual(iscol, faultIS, &isFaultCol);CHKERRQ(ierr);
if (isFaultRow && isFaultCol) {
@@ -112,12 +112,11 @@
// Make global preconditioner matrix
PetscMat jacobianMat = jacobian.matrix();
- PetscSection solutionSection = fields.solution().petscSection();
- DM dmMesh = fields.solution().mesh().dmMesh();
- PetscInt numFields;
+ PetscSection solutionSection = fields.solution().petscSection();assert(solutionSection);
+ PetscDM dmMesh = fields.solution().mesh().dmMesh();assert(dmMesh);
+ PetscInt numFields;
PetscErrorCode err;
- assert(solutionSection);assert(dmMesh);
err = DMGetNumFields(dmMesh, &numFields);CHECK_PETSC_ERROR(err);
if (formulation->splitFields() && formulation->useCustomConstraintPC() &&
numFields > 1) {
@@ -125,78 +124,19 @@
// and constraints exist.
PylithInt M, N, m, n;
- PetscErrorCode err = 0;
err = MatDestroy(&_jacobianPC);CHECK_PETSC_ERROR(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, &_jacobianPC);
- CHECK_PETSC_ERROR(err);
+ err = MatCreateShell(fields.mesh().comm(), m, n, M, N, &_ctx, &_jacobianPC);CHECK_PETSC_ERROR(err);
err = MatShellSetOperation(_jacobianPC, MATOP_GET_SUBMATRIX,
- (void (*)(void)) MyMatGetSubMatrix);
- CHECK_PETSC_ERROR(err);
+ (void (*)(void)) MyMatGetSubMatrix);CHECK_PETSC_ERROR(err);
} else {
_jacobianPC = jacobianMat;
- PetscErrorCode err = PetscObjectReference((PetscObject) jacobianMat);
- CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) jacobianMat);CHECK_PETSC_ERROR(err);
} // if/else
} // initialize
-static inline PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
-{
- switch(i) {
- case 0:
- switch(j) {
- case 0: return 0;
- case 1:
- switch(k) {
- case 0: return 0;
- case 1: return 0;
- case 2: return 1;
- }
- case 2:
- switch(k) {
- case 0: return 0;
- case 1: return -1;
- case 2: return 0;
- }
- }
- case 1:
- switch(j) {
- case 0:
- switch(k) {
- case 0: return 0;
- case 1: return 0;
- case 2: return -1;
- }
- case 1: return 0;
- case 2:
- switch(k) {
- case 0: return 1;
- case 1: return 0;
- case 2: return 0;
- }
- }
- case 2:
- switch(j) {
- case 0:
- switch(k) {
- case 0: return 0;
- case 1: return 1;
- case 2: return 0;
- }
- case 1:
- switch(k) {
- case 0: return -1;
- case 1: return 0;
- case 2: return 0;
- }
- case 2: return 0;
- }
- }
- return 0;
-}
-
// ----------------------------------------------------------------------
// Setup preconditioner for preconditioning with split fields.
void
@@ -208,19 +148,16 @@
assert(pc);
assert(formulation);
- DM dmMesh = fields.mesh().dmMesh();
+ PetscDM dmMesh = fields.mesh().dmMesh();assert(dmMesh);
MPI_Comm comm;
- assert(dmMesh);
- PetscSection solutionSection = fields.solution().petscSection();
- Vec solutionVec = fields.solution().localVector();
- Vec solutionGlobalVec = fields.solution().globalVector();
- PetscInt spaceDim, numFields;
+ PetscSection solutionSection = fields.solution().petscSection();assert(solutionSection);
+ PetscVec solutionVec = fields.solution().localVector();assert(solutionVec);
+ PetscVec solutionGlobalVec = fields.solution().globalVector();assert(solutionGlobalVec);
+ PetscInt spaceDim, numFields;
PetscErrorCode err;
-
const bool separateComponents = formulation->splitFieldComponents();
- assert(solutionSection);
err = PetscObjectGetComm((PetscObject) dmMesh, &comm);CHECK_PETSC_ERROR(err);
err = DMPlexGetDimension(dmMesh, &spaceDim);CHECK_PETSC_ERROR(err);
err = PetscSectionGetNumFields(solutionSection, &numFields);CHECK_PETSC_ERROR(err);
@@ -236,32 +173,27 @@
// and constraints exist.
// Get total number of DOF associated with constraints field split
- DM lagrangeDM;
- PetscSection lagrangeSection;
- PetscInt lagrangeFields[1] = {1}, nrows, ncols;
+ PetscDM lagrangeDM = NULL;
+ PetscSection lagrangeSection = NULL;
+ PetscInt lagrangeFields[1] = {1}, nrows, ncols;
err = MatDestroy(&_jacobianPCFault);CHECK_PETSC_ERROR(err);
- err = DMCreateSubDM(dmMesh, 1, lagrangeFields, PETSC_NULL, &lagrangeDM);CHECK_PETSC_ERROR(err);
+ err = DMCreateSubDM(dmMesh, 1, lagrangeFields, NULL, &lagrangeDM);CHECK_PETSC_ERROR(err);
err = DMGetDefaultGlobalSection(lagrangeDM, &lagrangeSection);CHECK_PETSC_ERROR(err);
err = PetscSectionGetStorageSize(lagrangeSection, &nrows);CHECK_PETSC_ERROR(err);
ncols = nrows;
err = MatCreate(comm, &_jacobianPCFault);CHECK_PETSC_ERROR(err);
- err = MatSetSizes(_jacobianPCFault, nrows, ncols,
- PETSC_DECIDE, PETSC_DECIDE);CHECK_PETSC_ERROR(err);
+ err = MatSetSizes(_jacobianPCFault, nrows, ncols, PETSC_DECIDE, PETSC_DECIDE);CHECK_PETSC_ERROR(err);
err = MatSetType(_jacobianPCFault, MATAIJ);CHECK_PETSC_ERROR(err);
err = MatSetFromOptions(_jacobianPCFault);CHECK_PETSC_ERROR(err);
// Allocate just the diagonal.
- err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1,
- PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = MatMPIAIJSetPreallocation(_jacobianPCFault, 1, PETSC_NULL,
- 0, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1, NULL);CHECK_PETSC_ERROR(err);
+ err = MatMPIAIJSetPreallocation(_jacobianPCFault, 1, NULL, 0, NULL);CHECK_PETSC_ERROR(err);
// Set preconditioning matrix in formulation
- formulation->customPCMatrix(_jacobianPCFault);
+ formulation->customPCMatrix(_jacobianPCFault);assert(_jacobianPCFault);
- assert(_jacobianPCFault);
-
_ctx.pc = *pc;
_ctx.A = jacobian.matrix();
switch (spaceDim) {
@@ -283,22 +215,19 @@
} // if
if (separateComponents) {
- throw std::runtime_error("Separate components is no longer supported");
- // constructFieldSplit(solutionSection, PETSC_DETERMINE, PETSC_NULL, PETSC_NULL,
- // sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionSection), precon, PETSC_NULL, solution.vector(), *pc);
+ throw std::logic_error("Separate components is no longer supported");
} else {
// Create rigid body null space.
- MatNullSpace nullsp;
- PetscObject field;
+ MatNullSpace nullsp = NULL;
+ PetscObject field = NULL;
PetscSection coordinateSection;
- Vec coordinateVec;
- PetscInt dim = spaceDim, vStart, vEnd;
- PetscVec mode[6];
+ PetscVec coordinateVec;
+ PetscInt dim = spaceDim, vStart, vEnd;
+ PetscVec mode[6];
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexGetCoordinateSection(dmMesh, &coordinateSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);CHECK_PETSC_ERROR(err);
- assert(coordinateSection);assert(coordinateVec);
+ err = DMPlexGetCoordinateSection(dmMesh, &coordinateSection);CHECK_PETSC_ERROR(err);assert(coordinateSection);
+ err = DMGetCoordinatesLocal(dmMesh, &coordinateVec);CHECK_PETSC_ERROR(err);assert(coordinateVec);
if (dim > 1) {
const int m = (dim * (dim + 1)) / 2;
for(int i = 0; i < m; ++i) {
@@ -322,22 +251,22 @@
err = VecSet(solutionVec, 0.0);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
PetscScalar values[3] = {0.0, 0.0, 0.0};
- PetscScalar *coords;
+ PetscScalar *coords = NULL;
- err = DMPlexVecGetClosure(dmMesh, coordinateSection, coordinateVec, v, PETSC_NULL, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecGetClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);CHECK_PETSC_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, PETSC_NULL, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecRestoreClosure(dmMesh, coordinateSection, coordinateVec, v, NULL, &coords);CHECK_PETSC_ERROR(err);
err = DMPlexVecSetClosure(dmMesh, solutionSection, solutionVec, v, values, INSERT_VALUES);CHECK_PETSC_ERROR(err);
}
err = DMLocalToGlobalBegin(dmMesh, solutionVec, INSERT_VALUES, mode[d]);CHECK_PETSC_ERROR(err);
err = DMLocalToGlobalEnd(dmMesh, solutionVec, INSERT_VALUES, mode[d]);CHECK_PETSC_ERROR(err);
} // for
for(int i = 0; i < dim; ++i) {
- err = VecNormalize(mode[i], PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = VecNormalize(mode[i], NULL);CHECK_PETSC_ERROR(err);
} // for
/* Orthonormalize system */
for(int i = dim; i < m; ++i) {
@@ -346,7 +275,7 @@
err = VecMDot(mode[i], i, mode, dots);CHECK_PETSC_ERROR(err);
for(int j = 0; j < i; ++j) dots[j] *= -1.0;
err = VecMAXPY(mode[i], i, dots, mode);CHECK_PETSC_ERROR(err);
- err = VecNormalize(mode[i], PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = VecNormalize(mode[i], NULL);CHECK_PETSC_ERROR(err);
} // for
err = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, &nullsp);CHECK_PETSC_ERROR(err);
for(int i = 0; i< m; ++i) {err = VecDestroy(&mode[i]);CHECK_PETSC_ERROR(err);}
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.hh 2013-03-20 01:37:35 UTC (rev 21587)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.hh 2013-03-20 02:13:06 UTC (rev 21588)
@@ -34,10 +34,10 @@
#include "pylith/utils/petscfwd.h" // USES PetscMat
typedef struct {
- PetscPC pc;
- PetscMat A;
+ PetscPC pc;
+ PetscMat A;
const char *faultFieldName;
- PetscMat faultA;
+ PetscMat faultA;
} FaultPreconCtx;
@@ -91,11 +91,10 @@
const topology::Jacobian& jacobian,
const topology::SolutionFields& fields);
- /**
+ /** :MATT: :TODO: DOCUMENT THIS.
*/
static
- int
- _epsilon(int i,
+ int _epsilon(int i,
int j,
int k);
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLinear.cc 2013-03-20 01:37:35 UTC (rev 21587)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLinear.cc 2013-03-20 02:13:06 UTC (rev 21588)
@@ -30,10 +30,6 @@
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverLinear::SolverLinear(void) :
_ksp(0)
@@ -60,10 +56,9 @@
// ----------------------------------------------------------------------
// Initialize solver.
void
-pylith::problems::SolverLinear::initialize(
- const topology::SolutionFields& fields,
- const topology::Jacobian& jacobian,
- Formulation* formulation)
+pylith::problems::SolverLinear::initialize(const topology::SolutionFields& fields,
+ const topology::Jacobian& jacobian,
+ Formulation* formulation)
{ // initialize
assert(formulation);
@@ -86,10 +81,9 @@
// ----------------------------------------------------------------------
// Solve the system.
void
-pylith::problems::SolverLinear::solve(
- topology::Field<topology::Mesh>* solution,
- topology::Jacobian* jacobian,
- const topology::Field<topology::Mesh>& residual)
+pylith::problems::SolverLinear::solve(topology::Field<topology::Mesh>* solution,
+ topology::Jacobian* jacobian,
+ const topology::Field<topology::Mesh>& residual)
{ // solve
assert(solution);
assert(jacobian);
@@ -108,15 +102,12 @@
_logger->eventEnd(scatterEvent);
_logger->eventBegin(setupEvent);
-
PetscErrorCode err = 0;
const PetscMat jacobianMat = jacobian->matrix();
if (!jacobian->valuesChanged()) {
- err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
- SAME_PRECONDITIONER);CHECK_PETSC_ERROR(err);
+ err = KSPSetOperators(_ksp, jacobianMat, jacobianMat, SAME_PRECONDITIONER);CHECK_PETSC_ERROR(err);
} else {
- err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
- DIFFERENT_NONZERO_PATTERN);CHECK_PETSC_ERROR(err);
+ err = KSPSetOperators(_ksp, jacobianMat, jacobianMat, DIFFERENT_NONZERO_PATTERN);CHECK_PETSC_ERROR(err);
} // else
jacobian->resetValuesChanged();
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc 2013-03-20 01:37:35 UTC (rev 21587)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc 2013-03-20 02:13:06 UTC (rev 21588)
@@ -22,15 +22,13 @@
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
#include "pylith/problems/Formulation.hh" // USES Formulation
#include "pylith/utils/EventLogger.hh" // USES EventLogger
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverLumped::SolverLumped(void)
{ // constructor
@@ -54,12 +52,11 @@
// ----------------------------------------------------------------------
// Initialize solver.
void
-pylith::problems::SolverLumped::initialize(
- const topology::SolutionFields& fields,
- const topology::Field<topology::Mesh>& jacobian,
- Formulation* formulation)
+pylith::problems::SolverLumped::initialize(const topology::SolutionFields& fields,
+ const topology::Field<topology::Mesh>& jacobian,
+ Formulation* formulation)
{ // initialize
- assert(0 != formulation);
+ assert(formulation);
_initializeLogger();
@@ -69,13 +66,12 @@
// ----------------------------------------------------------------------
// Solve the system.
void
-pylith::problems::SolverLumped::solve(
- topology::Field<topology::Mesh>* solution,
- const topology::Field<topology::Mesh>& jacobian,
- const topology::Field<topology::Mesh>& residual)
+pylith::problems::SolverLumped::solve(topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residual)
{ // solve
- assert(0 != solution);
- assert(0 != _formulation);
+ assert(solution);
+ assert(_formulation);
// solution = residual / jacobian
@@ -84,66 +80,44 @@
const int adjustEvent = _logger->eventId("SoLu adjust");
_logger->eventBegin(setupEvent);
- const spatialdata::geocoords::CoordSys* cs = solution->mesh().coordsys();
- assert(0 != cs);
+ const spatialdata::geocoords::CoordSys* cs = solution->mesh().coordsys();assert(cs);
const int spaceDim = cs->spaceDim();
// Get mesh vertices.
- DM dmMesh = solution->mesh().dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
-
- assert(dmMesh);
- err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ PetscDM dmMesh = solution->mesh().dmMesh(); assert(dmMesh);
+ topology::Stratum depthStratum(dmMesh, topology::Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
// Get sections.
- PetscSection solutionSection = solution->petscSection();
- Vec solutionVec = solution->localVector();
- PetscScalar *solutionArray;
- assert(solutionSection);assert(solutionVec);
-
- PetscSection jacobianSection = jacobian.petscSection();
- Vec jacobianVec = jacobian.localVector();
- PetscScalar *jacobianArray;
- assert(jacobianSection);assert(jacobianVec);
+ topology::VecVisitorMesh solutionVisitor(*solution);
+ PetscScalar* solutionArray = solutionVisitor.localArray();
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
- PetscScalar *residualArray;
- assert(residualSection);assert(residualVec);
-
+ topology::VecVisitorMesh jacobianVisitor(jacobian);
+ PetscScalar* jacobianArray = jacobianVisitor.localArray();
+
+ topology::VecVisitorMesh residualVisitor(residual);
+ PetscScalar* residualArray = residualVisitor.localArray();
+
_logger->eventEnd(setupEvent);
_logger->eventBegin(solveEvent);
- err = VecGetArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
- PetscInt jdof, joff;
+ const PetscInt joff = jacobianVisitor.sectionOffset(v);
+ assert(spaceDim == jacobianVisitor.sectionDof(v));
- err = PetscSectionGetDof(jacobianSection, v, &jdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(jacobianSection, v, &joff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == jdof);
- PetscInt rdof, roff;
+ const PetscInt roff = residualVisitor.sectionOffset(v);
+ assert(spaceDim == residualVisitor.sectionDof(v));
- err = PetscSectionGetDof(residualSection, v, &rdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(residualSection, v, &roff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == rdof);
- PetscInt sdof, soff;
+ const PetscInt soff = solutionVisitor.sectionOffset(v);
+ assert(spaceDim == solutionVisitor.sectionDof(v));
- err = PetscSectionGetDof(solutionSection, v, &sdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(solutionSection, v, &soff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == sdof);
-
for (int i=0; i < spaceDim; ++i) {
assert(jacobianArray[joff+i] != 0.0);
solutionArray[soff+i] = residualArray[roff+i] / jacobianArray[joff+i];
} // for
} // for
PetscLogFlops((vEnd - vStart) * spaceDim);
- err = VecRestoreArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
_logger->eventEnd(solveEvent);
_logger->eventBegin(adjustEvent);
@@ -165,7 +139,7 @@
pylith::problems::SolverLumped::_initializeLogger(void)
{ // initializeLogger
delete _logger; _logger = new utils::EventLogger;
- assert(0 != _logger);
+ assert(_logger);
_logger->className("SolverLumped");
_logger->initialize();
_logger->registerEvent("SoLu setup");
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2013-03-20 01:37:35 UTC (rev 21587)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2013-03-20 02:13:06 UTC (rev 21588)
@@ -47,10 +47,6 @@
#define PYLITH_CUSTOM_LINESEARCH 1
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverNonlinear::SolverNonlinear(void) :
_snes(0)
@@ -77,18 +73,17 @@
// ----------------------------------------------------------------------
// Initialize solver.
void
-pylith::problems::SolverNonlinear::initialize(
- const topology::SolutionFields& fields,
- const topology::Jacobian& jacobian,
- Formulation* formulation)
+pylith::problems::SolverNonlinear::initialize(const topology::SolutionFields& fields,
+ const topology::Jacobian& jacobian,
+ Formulation* formulation)
{ // initialize
- assert(0 != formulation);
+ assert(formulation);
_initializeLogger();
Solver::initialize(fields, jacobian, formulation);
PetscErrorCode err = 0;
- if (0 != _snes) {
+ if (_snes) {
err = SNESDestroy(&_snes); _snes = 0;
CHECK_PETSC_ERROR(err);
} // if
@@ -96,12 +91,10 @@
const topology::Field<topology::Mesh>& residual = fields.get("residual");
const PetscVec residualVec = residual.globalVector();
- err = SNESSetFunction(_snes, residualVec, reformResidual,
- (void*) formulation);
+ err = SNESSetFunction(_snes, residualVec, reformResidual, (void*) formulation);
CHECK_PETSC_ERROR(err);
- err = SNESSetJacobian(_snes, jacobian.matrix(), _jacobianPC, reformJacobian, (void*) formulation);
- CHECK_PETSC_ERROR(err);
+ err = SNESSetJacobian(_snes, jacobian.matrix(), _jacobianPC, reformJacobian, (void*) formulation);CHECK_PETSC_ERROR(err);
// Set default line search type to SNESSHELL and use our custom line search
PetscSNESLineSearch ls;
@@ -126,12 +119,11 @@
// ----------------------------------------------------------------------
// Solve the system.
void
-pylith::problems::SolverNonlinear::solve(
- topology::Field<topology::Mesh>* solution,
- topology::Jacobian* jacobian,
- const topology::Field<topology::Mesh>& residual)
+pylith::problems::SolverNonlinear::solve(topology::Field<topology::Mesh>* solution,
+ topology::Jacobian* jacobian,
+ const topology::Field<topology::Mesh>& residual)
{ // solve
- assert(0 != solution);
+ assert(solution);
const int solveEvent = _logger->eventId("SoNl solve");
const int scatterEvent = _logger->eventId("SoNl scatter");
@@ -163,9 +155,9 @@
PetscVec tmpResidualVec,
void* context)
{ // reformResidual
- assert(0 != context);
+ assert(context);
Formulation* formulation = (Formulation*) context;
- assert(0 != formulation);
+ assert(formulation);
// Make sure we have an admissible Lagrange multiplier (\lambda)
formulation->constrainSolnSpace(&tmpSolutionVec);
@@ -187,9 +179,9 @@
MatStructure* preconditionerLayout,
void* context)
{ // reformJacobian
- assert(0 != context);
+ assert(context);
Formulation* formulation = (Formulation*) context;
- assert(0 != formulation);
+ assert(formulation);
formulation->reformJacobian(&tmpSolutionVec);
@@ -213,8 +205,8 @@
PetscBool changed_y =PETSC_FALSE, changed_w =PETSC_FALSE;
PetscErrorCode ierr;
- Vec X, F, Y, W, G;
- SNES snes;
+ PetscVec X, F, Y, W, G;
+ PetscSNES snes;
PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev;
PetscReal lambda, lambdatemp, lambdaprev, minlambda, maxstep, initslope, alpha, stol;
PetscReal t1, t2, a, b, d;
@@ -225,7 +217,7 @@
PetscViewer monitor;
PetscInt max_its, count;
PetscSNESLineSearch_BT *bt;
- Mat jac;
+ PetscMat jac;
PetscFunctionBegin;
@@ -512,7 +504,7 @@
PetscVec initialGuessVec,
void *lsctx)
{ // initialGuess
- VecSet(initialGuessVec, 0.0);
+ PetscErrorCode err = VecSet(initialGuessVec, 0.0);CHECK_PETSC_ERROR(err);
} // initialGuess
// ----------------------------------------------------------------------
@@ -521,7 +513,7 @@
pylith::problems::SolverNonlinear::_initializeLogger(void)
{ // initializeLogger
delete _logger; _logger = new utils::EventLogger;
- assert(0 != _logger);
+ assert(_logger);
_logger->className("SolverNonlinear");
_logger->initialize();
_logger->registerEvent("SoNl setup");
More information about the CIG-COMMITS
mailing list