[cig-commits] r16039 - in short/3D/PyLith/branches/pylith-friction: libsrc/problems playpen/friction
brad at geodynamics.org
brad at geodynamics.org
Wed Nov 25 11:07:27 PST 2009
Author: brad
Date: 2009-11-25 11:07:26 -0800 (Wed, 25 Nov 2009)
New Revision: 16039
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4-axial.cfg
Log:
Modification to customized line search. Need to integrate constraining solution into line search.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-25 16:54:23 UTC (rev 16038)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverNonlinear.cc 2009-11-25 19:07:26 UTC (rev 16039)
@@ -179,20 +179,6 @@
ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
*flag = PETSC_TRUE;
- // ======================================================================
- // Code to constrain solution space.
- //
- // Matt- It seems like I should adjust both x (current iterate) and
- // y (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(&x); // Adjust x? y? both?
- // ======================================================================
-
ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr);
if (!*ynorm) {
ierr = PetscInfo(snes,"Search direction and size is 0\n");CHKERRQ(ierr);
@@ -203,7 +189,8 @@
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 = 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;
}
@@ -221,32 +208,42 @@
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);
+ 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;
- }
+ } // if
+ // TEMPORARY: update w?
+
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 (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);
+ 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;
- }
+ } // if
- /* Fit points with quadratic */
+ // 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;
+ 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) {
@@ -254,92 +251,142 @@
*flag = PETSC_FALSE;
snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
goto theend1;
- }
+ } // if
+ // TEMPORARY: update w?
+
ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
if (snes->domainerror) {
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
PetscFunctionReturn(0);
- }
+ } // if
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;
- }
+ } // if
- /* Fit points with cubic */
+ // 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);
+ 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);
+ 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);
- }
+ } // if/else
lambdaprev = lambda;
gnormprev = *gnorm;
- if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
- if (lambdatemp <= .1*lambda) lambda = .1*lambda;
- else lambda = lambdatemp;
+ 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);
+ 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;
- }
+ } // if
+ // TEMPORARY: update w?
+
ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
if (snes->domainerror) {
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
PetscFunctionReturn(0);
- }
+ } // if
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);
+ 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);
- }
+ ierr = PetscInfo2(snes,"Cubic step no good, shrinking lambda, "
+ "current gnorem %G lambda=%18.16e\n",
+ *gnorm,lambda);CHKERRQ(ierr);
+ } // if/else
count++;
- }
+ } // while
+
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 = (*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 */
+ // TEMPORARY: update w?
+
+ 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");
+ 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);
- }
- }
+ } // if
+ } // if
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
+ // ======================================================================
+ // Code to constrain solution space.
+ assert(0 != lsctx);
+ Formulation* formulation = (Formulation*) lsctx;
+ assert(0 != formulation);
+ formulation->constrainSolnSpace(&w);
+
+ 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);
+ // ======================================================================
+
PetscFunctionReturn(0);
-}
+} // lineSearch
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4-axial.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4-axial.cfg 2009-11-25 16:54:23 UTC (rev 16038)
+++ short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4-axial.cfg 2009-11-25 19:07:26 UTC (rev 16039)
@@ -67,7 +67,7 @@
# The properties for this material are given in the spatial database file
# 'matprops.spatialdb'.
-db_properties.iohandler.filename = matprops.spatialdb
+db_properties.iohandler.filename = matelastic2D.spatialdb
# Set cell type to quadrilateral (2-d Lagrange).
quadrature.cell = pylith.feassemble.FIATLagrange
More information about the CIG-COMMITS
mailing list