[cig-commits] r16034 - in short/3D/PyLith/branches/pylith-friction/libsrc: feassemble problems
knepley at geodynamics.org
knepley at geodynamics.org
Tue Nov 24 17:53:45 PST 2009
Author: knepley
Date: 2009-11-24 17:53:44 -0800 (Tue, 24 Nov 2009)
New Revision: 16034
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
Log:
Added a custom line seach (just cubic now) and some support code for Laplacian preconditioning
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicit.cc 2009-11-24 22:44:31 UTC (rev 16033)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicit.cc 2009-11-25 01:53:44 UTC (rev 16034)
@@ -355,19 +355,26 @@
// Set variables dependent on dimension of cell
totalStrain_fn_type calcTotalStrainFn;
elasticityJacobian_fn_type elasticityJacobianFn;
+ elasticityJacobian_fn_type elasticityPreconFn;
if (1 == cellDim) {
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian1D;
+ elasticityPreconFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityPrecon1D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian2D;
+ elasticityPreconFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityPrecon2D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian3D;
+ elasticityPreconFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityPrecon3D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
} else
@@ -523,6 +530,21 @@
&_cellMatrix[0], ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
_logger->eventEnd(updateEvent);
+#if 0
+ _logger->eventBegin(computeEvent);
+ // Get laplacian matrix at quadrature points for this cell
+ CALL_MEMBER_FN(*this, elasticityPreconFn)(elasticConsts);
+ _logger->eventEnd(computeEvent);
+
+ // Assemble cell contribution into PETSc preconditioner matrix.
+ _logger->eventBegin(updateEvent);
+ jacobianVisitor.clear();
+ PetscErrorCode err = updateOperator(preconMat, *sieveMesh->getSieve(),
+ jacobianVisitor, *c_iter,
+ &_cellMatrix[0], ADD_VALUES);
+ CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+ _logger->eventEnd(updateEvent);
+#endif
} // for
_needNewJacobian = false;
_material->resetNeedNewJacobian();
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.cc 2009-11-24 22:44:31 UTC (rev 16033)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.cc 2009-11-25 01:53:44 UTC (rev 16034)
@@ -759,6 +759,38 @@
} // _elasticityJacobian1D
// ----------------------------------------------------------------------
+// Integrate laplacian term in Jacobian preconditioner for 1-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityPrecon1D(
+ const double_array& elasticConsts)
+{ // _elasticityPrecon1D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(1 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basisDeriv[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basisDeriv[iQ+jBasis];
+ const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
+ const int jBlock = jBasis*spaceDim;
+ _cellMatrix[iBlock+jBlock] += valIJ;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*2)));
+} // _elasticityPrecon1D
+
+// ----------------------------------------------------------------------
// Integrate elasticity term in Jacobian for 2-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityJacobian2D(
@@ -823,6 +855,49 @@
} // _elasticityJacobian2D
// ----------------------------------------------------------------------
+// Integrate laplacian term in Jacobian preconditioner for 2-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityPrecon2D(
+ const double_array& elasticConsts)
+{ // _elasticityPrecon2D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(2 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int numConsts = 6;
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ // Delta_ij = C_ijkl * e_kl
+ // = C_ijlk * 0.5 (u_k,l + u_l,k)
+ // = 0.5 * C_ijkl * (u_k,l + u_l,k)
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const double Nip = wt*basisDeriv[iQ+iBasis*spaceDim ];
+ const double Niq = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ const int iBlock = (iBasis*spaceDim ) * (numBasis*spaceDim);
+ const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double Njp = basisDeriv[iQ+jBasis*spaceDim ];
+ const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
+ const int jBlock = (jBasis*spaceDim );
+ const int jBlock1 = (jBasis*spaceDim+1);
+ _cellMatrix[iBlock +jBlock ] += Nip*Njp;
+ _cellMatrix[iBlock1+jBlock1] += Niq*Njq;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(4))));
+} // _elasticityPrecon2D
+
+// ----------------------------------------------------------------------
// Integrate elasticity term in Jacobian for 3-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityJacobian3D(
@@ -936,7 +1011,53 @@
} // _elasticityJacobian3D
// ----------------------------------------------------------------------
+// Integrate laplacian term in Jacobian preconditioner for 3-D cells.
void
+pylith::feassemble::IntegratorElasticity::_elasticityPrecon3D(
+ const double_array& elasticConsts)
+{ // _elasticityPrecon3D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(3 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int numConsts = 21;
+
+ // Compute Jacobian for consistent tangent matrix
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const double Nip = wt*basisDeriv[iQ+iBasis*spaceDim+0];
+ const double Niq = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ const double Nir = wt*basisDeriv[iQ+iBasis*spaceDim+2];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double Njp = basisDeriv[iQ+jBasis*spaceDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
+ const double Njr = basisDeriv[iQ+jBasis*spaceDim+2];
+ const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
+ const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
+ const int iBlock2 = (iBasis*spaceDim+2) * (numBasis*spaceDim);
+ const int jBlock = jBasis*spaceDim;
+ const int jBlock1 = jBasis*spaceDim+1;
+ const int jBlock2 = jBasis*spaceDim+2;
+ _cellMatrix[iBlock +jBlock ] += Nip*Njp;
+ _cellMatrix[iBlock1+jBlock1] += Niq*Njq;
+ _cellMatrix[iBlock2+jBlock2] += Nir*Njr;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(6))));
+} // _elasticityPrecon3D
+
+// ----------------------------------------------------------------------
+void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D(
double_array* strain,
const double_array& basisDeriv,
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.hh 2009-11-24 22:44:31 UTC (rev 16033)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticity.hh 2009-11-25 01:53:44 UTC (rev 16034)
@@ -186,6 +186,24 @@
*/
void _elasticityJacobian3D(const double_array& elasticConsts);
+ /** Integrate laplacian term in Jacobian preconditioner for 1-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityPrecon1D(const double_array& elasticConsts);
+
+ /** Integrate laplacian term in Jacobian preconditioner for 2-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityPrecon2D(const double_array& elasticConsts);
+
+ /** Integrate laplacian term in Jacobian preconditioner for 3-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityPrecon3D(const double_array& elasticConsts);
+
/** Compute total strain in at quadrature points of a cell.
*
* @param strain Strain tensor at quadrature points.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-24 22:44:31 UTC (rev 16033)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-25 01:53:44 UTC (rev 16034)
@@ -22,6 +22,218 @@
#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
+
+ 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
+
+ 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.
+
+ 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] )
+
+
+ Notes:
+ This line search is taken from "Numerical Methods for Unconstrained
+ Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
+
+ Level: advanced
+
+.keywords: SNES, nonlinear, line search, cubic
+
+.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.
+ */
+
+ PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
+ PetscReal minlambda,lambda,lambdatemp;
+#if defined(PETSC_USE_COMPLEX)
+ PetscScalar cinitslope;
+#endif
+ PetscErrorCode ierr;
+ PetscInt count;
+ SNES_LS *neP = (SNES_LS*)snes->data;
+ PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;
+
+ PetscFunctionBegin;
+ ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ *flag = PETSC_TRUE;
+
+ ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr);
+ if (!*ynorm) {
+ ierr = PetscInfo(snes,"Search direction and size is 0\n");CHKERRQ(ierr);
+ *gnorm = fnorm;
+ ierr = VecCopy(x,w);CHKERRQ(ierr);
+ ierr = VecCopy(f,g);CHKERRQ(ierr);
+ *flag = PETSC_FALSE;
+ goto theend1;
+ }
+ if (*ynorm > neP->maxstep) { /* Step too big, so scale back */
+ ierr = PetscInfo2(snes,"Scaling step by %G old ynorm %G\n",neP->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr);
+ ierr = VecScale(y,neP->maxstep/(*ynorm));CHKERRQ(ierr);
+ *ynorm = neP->maxstep;
+ }
+ ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr);
+ minlambda = neP->minlambda/rellength;
+ ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr);
+#if defined(PETSC_USE_COMPLEX)
+ ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr);
+ initslope = PetscRealPart(cinitslope);
+#else
+ ierr = VecDot(f,w,&initslope);CHKERRQ(ierr);
+#endif
+ if (initslope > 0.0) initslope = -initslope;
+ if (initslope == 0.0) initslope = -1.0;
+
+ ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr);
+ if (snes->nfuncs >= snes->max_funcs) {
+ ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
+ *flag = PETSC_FALSE;
+ snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
+ goto theend1;
+ }
+ ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
+ if (snes->domainerror) {
+ ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+ }
+ ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
+ if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr);
+ if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + neP->alpha*initslope) { /* Sufficient reduction */
+ ierr = PetscInfo2(snes,"Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr);
+ goto theend1;
+ }
+
+ /* Fit points with quadratic */
+ lambda = 1.0;
+ lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
+ lambdaprev = lambda;
+ gnormprev = *gnorm;
+ if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
+ if (lambdatemp <= .1*lambda) lambda = .1*lambda;
+ else lambda = lambdatemp;
+
+ ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr);
+ if (snes->nfuncs >= snes->max_funcs) {
+ ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
+ *flag = PETSC_FALSE;
+ snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
+ goto theend1;
+ }
+ ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
+ if (snes->domainerror) {
+ ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+ }
+ ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
+ if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ ierr = PetscInfo1(snes,"gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr);
+ if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*neP->alpha*initslope) { /* sufficient reduction */
+ ierr = PetscInfo1(snes,"Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr);
+ goto theend1;
+ }
+
+ /* Fit points with cubic */
+ count = 1;
+ while (PETSC_TRUE) {
+ if (lambda <= minlambda) {
+ ierr = PetscInfo1(snes,"Unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
+ ierr = PetscInfo6(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,minlambda,lambda,initslope);CHKERRQ(ierr);
+ *flag = PETSC_FALSE;
+ break;
+ }
+ t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
+ t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
+ a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
+ b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
+ d = b*b - 3*a*initslope;
+ if (d < 0.0) d = 0.0;
+ if (a == 0.0) {
+ lambdatemp = -initslope/(2.0*b);
+ } else {
+ lambdatemp = (-b + sqrt(d))/(3.0*a);
+ }
+ lambdaprev = lambda;
+ gnormprev = *gnorm;
+ if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
+ if (lambdatemp <= .1*lambda) lambda = .1*lambda;
+ else lambda = lambdatemp;
+ ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr);
+ if (snes->nfuncs >= snes->max_funcs) {
+ ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
+ ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr);
+ *flag = PETSC_FALSE;
+ snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
+ break;
+ }
+ ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
+ if (snes->domainerror) {
+ ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+ }
+ ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr);
+ if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*neP->alpha*initslope) { /* is reduction enough? */
+ ierr = PetscInfo2(snes,"Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr);
+ break;
+ } else {
+ ierr = PetscInfo2(snes,"Cubic step no good, shrinking lambda, current gnorem %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr);
+ }
+ count++;
+ }
+ theend1:
+ /* Optional user-defined check for line search step validity */
+ if (neP->postcheckstep && *flag) {
+ ierr = (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);CHKERRQ(ierr);
+ if (changed_y) {
+ ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr);
+ }
+ if (changed_y || changed_w) { /* recompute the function if the step has changed */
+ ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
+ if (snes->domainerror) {
+ ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+ }
+ ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr);
+ if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
+ ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr);
+ ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr);
+ ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr);
+ }
+ }
+ ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
// ----------------------------------------------------------------------
// Constructor
pylith::problems::SolverNonlinear::SolverNonlinear(void) :
@@ -80,6 +292,7 @@
CHECK_PETSC_ERROR(err);
err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
+ err = SNESLineSearchSet(_snes, FrictionLineSearchCubic, this); CHECK_PETSC_ERROR(err);
} // initialize
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list