[cig-commits] r16036 - short/3D/PyLith/branches/pylith-friction/libsrc/problems
brad at geodynamics.org
brad at geodynamics.org
Wed Nov 25 08:46:26 PST 2009
Author: brad
Date: 2009-11-25 08:46:26 -0800 (Wed, 25 Nov 2009)
New Revision: 16036
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.hh
Log:
Added use of customized line search for nonlinear solver.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-11-25 14:25:45 UTC (rev 16035)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-11-25 16:46:26 UTC (rev 16036)
@@ -138,36 +138,12 @@
solution.scatterVectorToSection(*tmpSolutionVec);
} // if
- // ===================================================================
- // :KLUDGE: WAITING FOR MATT TO IMPLEMENT DIFFERENT LINE SEARCH
- // Need a customized replacement routine for line search that will be
- // set via SNES_SetLineSearch().
- int numIntegrators = _meshIntegrators.size();
- assert(numIntegrators > 0); // must have at least 1 bulk integrator
- for (int i=0; i < numIntegrators; ++i) {
- _meshIntegrators[i]->timeStep(_dt);
- _meshIntegrators[i]->constrainSolnSpace(_fields, _t, *_jacobian);
- } // for
- numIntegrators = _submeshIntegrators.size();
- for (int i=0; i < numIntegrators; ++i) {
- _submeshIntegrators[i]->timeStep(_dt);
- _submeshIntegrators[i]->constrainSolnSpace(_fields, _t, *_jacobian);
- } // for
-
- // Update PETScVec of solution for changes to Lagrange multipliers.
- if (0 != tmpSolutionVec) {
- topology::Field<topology::Mesh>& solution = _fields->solution();
- solution.scatterSectionToVector(*tmpSolutionVec);
- } // if
- // :KLUDGE: END
- // ===================================================================
-
// Set residual to zero.
topology::Field<topology::Mesh>& residual = _fields->get("residual");
residual.zero();
// Add in contributions that require assembly.
- numIntegrators = _meshIntegrators.size();
+ int numIntegrators = _meshIntegrators.size();
assert(numIntegrators > 0); // must have at least 1 bulk integrator
for (int i=0; i < numIntegrators; ++i) {
_meshIntegrators[i]->timeStep(_dt);
@@ -293,6 +269,39 @@
} // reformJacobian
// ----------------------------------------------------------------------
+// Constrain solution space.
+void
+pylith::problems::Formulation::constrainSolnSpace(const PetscVec* tmpSolutionVec)
+{ // constrainSolnSpace
+ assert(0 != tmpSolutionVec);
+ assert(0 != _fields);
+
+ // Update section view of field.
+ if (0 != tmpSolutionVec) {
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*tmpSolutionVec);
+ } // if
+
+ int numIntegrators = _meshIntegrators.size();
+ assert(numIntegrators > 0); // must have at least 1 bulk integrator
+ for (int i=0; i < numIntegrators; ++i) {
+ _meshIntegrators[i]->timeStep(_dt);
+ _meshIntegrators[i]->constrainSolnSpace(_fields, _t, *_jacobian);
+ } // for
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i) {
+ _submeshIntegrators[i]->timeStep(_dt);
+ _submeshIntegrators[i]->constrainSolnSpace(_fields, _t, *_jacobian);
+ } // for
+
+ // Update PETScVec of solution for changes to Lagrange multipliers.
+ if (0 != tmpSolutionVec) {
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterSectionToVector(*tmpSolutionVec);
+ } // if
+} // constrainSolnSpace
+
+// ----------------------------------------------------------------------
// Adjust solution from solver with lumped Jacobian to match Lagrange
// multiplier constraints.
void
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh 2009-11-25 14:25:45 UTC (rev 16035)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh 2009-11-25 16:46:26 UTC (rev 16036)
@@ -106,10 +106,10 @@
/** Reform system residual.
*
* @param tmpResidualVec Temporary PETSc vector for residual.
- * @param tmpSolveSolnVec Temporary PETSc vector for solution.
+ * @param tmpSolutionVec Temporary PETSc vector for solution.
*/
void reformResidual(const PetscVec* tmpResidualVec =0,
- const PetscVec* tmpSolveSolnVec =0);
+ const PetscVec* tmpSolutionVec =0);
/* Reform system Jacobian.
*
@@ -121,6 +121,12 @@
*/
void reformJacobianLumped(void);
+ /** Constrain solution space.
+ *
+ * @param tmpSolutionVec Temporary PETSc vector for solution.
+ */
+ void constrainSolnSpace(const PetscVec* tmpSolutionVec);
+
/** Adjust solution from solver with lumped Jacobian to match Lagrange
* multiplier constraints.
*/
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-25 14:25:45 UTC (rev 16035)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-25 16:46:26 UTC (rev 16036)
@@ -23,55 +23,147 @@
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
#include "../src/snes/impls/ls/ls.h"
-#undef __FUNCT__
-#define __FUNCT__ "FrictionLineSearchCubic"
-/*@C
- FrictionLineSearchCubic - Performs a cubic line search (default line search method).
- Collective on SNES
+// ----------------------------------------------------------------------
+// Constructor
+pylith::problems::SolverNonlinear::SolverNonlinear(void) :
+ _snes(0)
+{ // constructor
+} // constructor
- Input Parameters:
-+ snes - nonlinear context
-. lsctx - optional context for line search (not used here)
-. x - current iterate
-. f - residual evaluated at x
-. y - search direction
-. w - work vector
-. fnorm - 2-norm of f
-- xnorm - norm of x if known, otherwise 0
+// ----------------------------------------------------------------------
+// Destructor
+pylith::problems::SolverNonlinear::~SolverNonlinear(void)
+{ // destructor
+ deallocate();
+} // destructor
- Output Parameters:
-+ g - residual evaluated at new iterate y
-. w - new iterate
-. gnorm - 2-norm of g
-. ynorm - 2-norm of search length
-- flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.
+// ----------------------------------------------------------------------
+// Deallocate data structures.
+void
+pylith::problems::SolverNonlinear::deallocate(void)
+{ // deallocate
+ Solver::deallocate();
- Options Database Key:
-+ -snes_ls cubic - Activates SNESLineSearchCubic()
-. -snes_ls_alpha <alpha> - Sets alpha
-. -snes_ls_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y)
-- -snes_ls_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda/ max_i ( y[i]/x[i] )
+ if (0 != _snes) {
+ PetscErrorCode err = SNESDestroy(_snes); _snes = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Initialize solver.
+void
+pylith::problems::SolverNonlinear::initialize(
+ const topology::SolutionFields& fields,
+ const topology::Jacobian& jacobian,
+ Formulation* formulation)
+{ // initialize
+ assert(0 != formulation);
-
- Notes:
- This line search is taken from "Numerical Methods for Unconstrained
- Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
+ Solver::initialize(fields, jacobian, formulation);
- Level: advanced
+ PetscErrorCode err = 0;
+ if (0 != _snes) {
+ err = SNESDestroy(_snes); _snes = 0;
+ CHECK_PETSC_ERROR(err);
+ } // if
+ err = SNESCreate(fields.mesh().comm(), &_snes); CHECK_PETSC_ERROR(err);
-.keywords: SNES, nonlinear, line search, cubic
+ const topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const PetscVec residualVec = residual.vector();
+ err = SNESSetFunction(_snes, residualVec, reformResidual,
+ (void*) formulation);
+ CHECK_PETSC_ERROR(err);
-.seealso: SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms()
-@*/
-PetscErrorCode PETSCSNES_DLLEXPORT FrictionLineSearchCubic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
-{
- /*
- Note that for line search purposes we work with with the related
- minimization problem:
- min z(x): R^n -> R,
- where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
- */
+ PetscMat jacobianMat = jacobian.matrix();
+ err = SNESSetJacobian(_snes, jacobianMat, jacobianMat, reformJacobian,
+ (void*) formulation);
+ CHECK_PETSC_ERROR(err);
+
+ err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
+ err = SNESLineSearchSet(_snes, lineSearch,
+ (void*) formulation); CHECK_PETSC_ERROR(err);
+} // initialize
+
+// ----------------------------------------------------------------------
+// Solve the system.
+void
+pylith::problems::SolverNonlinear::solve(
+ topology::Field<topology::Mesh>* solution,
+ const topology::Jacobian& jacobian,
+ const topology::Field<topology::Mesh>& residual)
+{ // solve
+ assert(0 != solution);
+
+ PetscErrorCode err = 0;
+
+ const PetscVec solutionVec = solution->vector();
+ err = SNESSolve(_snes, PETSC_NULL, solutionVec); CHECK_PETSC_ERROR(err);
+
+ // Update section view of field.
+ solution->scatterVectorToSection();
+} // solve
+
+// ----------------------------------------------------------------------
+// Generic C interface for reformResidual for integration with
+// PETSc SNES solvers.
+PetscErrorCode
+pylith::problems::SolverNonlinear::reformResidual(PetscSNES snes,
+ PetscVec tmpSolutionVec,
+ PetscVec tmpResidualVec,
+ void* context)
+{ // reformResidual
+ assert(0 != context);
+ Formulation* formulation = (Formulation*) context;
+ assert(0 != formulation);
+
+ // Reform residual
+ formulation->reformResidual(&tmpResidualVec, &tmpSolutionVec);
+
+ return 0;
+} // reformResidual
+
+// ----------------------------------------------------------------------
+// Generic C interface for reformJacobian for integration with
+// PETSc SNES solvers.
+PetscErrorCode
+pylith::problems::SolverNonlinear::reformJacobian(PetscSNES snes,
+ PetscVec tmpSolutionVec,
+ PetscMat* jacobianMat,
+ PetscMat* preconditionerMat,
+ MatStructure* preconditionerLayout,
+ void* context)
+{ // reformJacobian
+ assert(0 != context);
+ Formulation* formulation = (Formulation*) context;
+ assert(0 != formulation);
+
+ formulation->reformJacobian(&tmpSolutionVec);
+
+ return 0;
+} // reformJacobian
+
+// ----------------------------------------------------------------------
+// Generic C interface for customized PETSc line search.
+PetscErrorCode
+pylith::problems::SolverNonlinear::lineSearch(PetscSNES snes,
+ void *lsctx,
+ PetscVec x,
+ PetscVec f,
+ PetscVec g,
+ PetscVec y,
+ PetscVec w,
+ PetscReal fnorm,
+ PetscReal xnorm,
+ PetscReal *ynorm,
+ PetscReal *gnorm,
+ PetscTruth *flag)
+{ // lineSearch
+ // Note that for line search purposes we work with with the related
+ // minimization problem:
+ // min z(x): R^n -> R,
+ // where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
PetscReal minlambda,lambda,lambdatemp;
@@ -87,6 +179,20 @@
ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
*flag = PETSC_TRUE;
+ // ======================================================================
+ // Code to constraint solution space.
+ //
+ // Matt- It seems like I should adjust both x (current iterate) and
+ // the search direction, with the Lagrange multipliers in x modified
+ // to conform to the constraints imposed by friction, and y adjusted
+ // so that the DOF associated with the Lagrange multipliers are zero
+ // (search direction doesn't change the Lagrange multipliers).
+ assert(0 != lsctx);
+ Formulation* formulation = (Formulation*) lsctx;
+ assert(0 != formulation);
+ formulation->constrainSolnSpace(&y);
+ // ======================================================================
+
ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr);
if (!*ynorm) {
ierr = PetscInfo(snes,"Search direction and size is 0\n");CHKERRQ(ierr);
@@ -231,127 +337,9 @@
}
}
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+
PetscFunctionReturn(0);
}
-// ----------------------------------------------------------------------
-// Constructor
-pylith::problems::SolverNonlinear::SolverNonlinear(void) :
- _snes(0)
-{ // constructor
-} // constructor
-// ----------------------------------------------------------------------
-// Destructor
-pylith::problems::SolverNonlinear::~SolverNonlinear(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate data structures.
-void
-pylith::problems::SolverNonlinear::deallocate(void)
-{ // deallocate
- Solver::deallocate();
-
- if (0 != _snes) {
- PetscErrorCode err = SNESDestroy(_snes); _snes = 0;
- CHECK_PETSC_ERROR(err);
- } // if
-} // deallocate
-
-// ----------------------------------------------------------------------
-// Initialize solver.
-void
-pylith::problems::SolverNonlinear::initialize(
- const topology::SolutionFields& fields,
- const topology::Jacobian& jacobian,
- Formulation* formulation)
-{ // initialize
- assert(0 != formulation);
-
- Solver::initialize(fields, jacobian, formulation);
-
- PetscErrorCode err = 0;
- if (0 != _snes) {
- err = SNESDestroy(_snes); _snes = 0;
- CHECK_PETSC_ERROR(err);
- } // if
- err = SNESCreate(fields.mesh().comm(), &_snes); CHECK_PETSC_ERROR(err);
-
- const topology::Field<topology::Mesh>& residual = fields.get("residual");
- const PetscVec residualVec = residual.vector();
- err = SNESSetFunction(_snes, residualVec, reformResidual,
- (void*) formulation);
- CHECK_PETSC_ERROR(err);
-
- PetscMat jacobianMat = jacobian.matrix();
- err = SNESSetJacobian(_snes, jacobianMat, jacobianMat, reformJacobian,
- (void*) formulation);
- CHECK_PETSC_ERROR(err);
-
- err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
- err = SNESLineSearchSet(_snes, FrictionLineSearchCubic, this); CHECK_PETSC_ERROR(err);
-} // initialize
-
-// ----------------------------------------------------------------------
-// Solve the system.
-void
-pylith::problems::SolverNonlinear::solve(
- topology::Field<topology::Mesh>* solution,
- const topology::Jacobian& jacobian,
- const topology::Field<topology::Mesh>& residual)
-{ // solve
- assert(0 != solution);
-
- PetscErrorCode err = 0;
-
- const PetscVec solutionVec = solution->vector();
- err = SNESSolve(_snes, PETSC_NULL, solutionVec); CHECK_PETSC_ERROR(err);
-
- // Update section view of field.
- solution->scatterVectorToSection();
-} // solve
-
-// ----------------------------------------------------------------------
-// Generic C interface for reformResidual for integration with
-// PETSc SNES solvers.
-PetscErrorCode
-pylith::problems::SolverNonlinear::reformResidual(PetscSNES snes,
- PetscVec tmpSolutionVec,
- PetscVec tmpResidualVec,
- void* context)
-{ // reformResidual
- assert(0 != context);
- Formulation* formulation = (Formulation*) context;
- assert(0 != formulation);
-
- // Reform residual
- formulation->reformResidual(&tmpResidualVec, &tmpSolutionVec);
-
- return 0;
-} // reformResidual
-
-// ----------------------------------------------------------------------
-// Generic C interface for reformJacobian for integration with
-// PETSc SNES solvers.
-PetscErrorCode
-pylith::problems::SolverNonlinear::reformJacobian(PetscSNES snes,
- PetscVec tmpSolutionVec,
- PetscMat* jacobianMat,
- PetscMat* preconditionerMat,
- MatStructure* preconditionerLayout,
- void* context)
-{ // reformJacobian
- assert(0 != context);
- Formulation* formulation = (Formulation*) context;
- assert(0 != formulation);
-
- formulation->reformJacobian(&tmpSolutionVec);
-
- return 0;
-} // reformJacobian
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.hh 2009-11-25 14:25:45 UTC (rev 16035)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.hh 2009-11-25 16:46:26 UTC (rev 16036)
@@ -104,6 +104,37 @@
MatStructure* preconditionerLayout,
void* context);
+ /** Generic C interface for customized PETSc line search.
+ *
+ * @param snes PETSc SNES solver.
+ * @param lsctx Optional context for line search (not used here).
+ * @param x Current iterate.
+ * @param f Residual evaluated at x.
+ * @param y Search direction.
+ * @param w Work vector
+ * @param f 2-norm of f.
+ * @param xnorm Norm of x if known, otherwise 0.
+ * @param g Residual evaluated at new iterate y.
+ * @param w New iterate.
+ * @param gnorm 2-norm of g.
+ * @param ynorm 2-norm of search length.
+ * @param PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.
+ * @returns PETSc error code.
+ */
+ static
+ PetscErrorCode lineSearch(PetscSNES snes,
+ void *lsctx,
+ PetscVec x,
+ PetscVec f,
+ PetscVec g,
+ PetscVec y,
+ PetscVec w,
+ PetscReal fnorm,
+ PetscReal xnorm,
+ PetscReal *ynorm,
+ PetscReal *gnorm,
+ PetscTruth *flag);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list