1 2 #include <../src/snes/impls/tr/trimpl.h> /*I "petscsnes.h" I*/ 3 4 typedef struct { 5 SNES snes; 6 /* Information on the regular SNES convergence test; which may have been user provided */ 7 PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 8 PetscErrorCode (*convdestroy)(void*); 9 void *convctx; 10 } SNES_TR_KSPConverged_Ctx; 11 12 static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx) 13 { 14 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 15 SNES snes = ctx->snes; 16 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 17 Vec x; 18 PetscReal nrm; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);CHKERRQ(ierr); 23 if (*reason) { 24 ierr = PetscInfo2(snes,"Default or user provided convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);CHKERRQ(ierr); 25 } 26 /* Determine norm of solution */ 27 ierr = KSPBuildSolution(ksp,NULL,&x);CHKERRQ(ierr); 28 ierr = VecNorm(x,NORM_2,&nrm);CHKERRQ(ierr); 29 if (nrm >= neP->delta) { 30 ierr = PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);CHKERRQ(ierr); 31 *reason = KSP_CONVERGED_STEP_LENGTH; 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx) 37 { 38 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = (*ctx->convdestroy)(ctx->convctx);CHKERRQ(ierr); 43 ierr = PetscFree(ctx);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 /* ---------------------------------------------------------------- */ 48 /* 49 SNESTR_Converged_Private -test convergence JUST for 50 the trust region tolerance. 51 52 */ 53 static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 54 { 55 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 *reason = SNES_CONVERGED_ITERATING; 60 if (neP->delta < xnorm * snes->deltatol) { 61 ierr = PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);CHKERRQ(ierr); 62 *reason = SNES_DIVERGED_TR_DELTA; 63 } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 64 ierr = PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);CHKERRQ(ierr); 65 *reason = SNES_DIVERGED_FUNCTION_COUNT; 66 } 67 PetscFunctionReturn(0); 68 } 69 70 /*@C 71 SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined. 72 Allows the user a chance to change or override the decision of the line search routine. 73 74 Logically Collective on snes 75 76 Input Parameters: 77 + snes - the nonlinear solver object 78 . func - [optional] function evaluation routine, see SNESNewtonTRPreCheck() for the calling sequence 79 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 80 81 Level: intermediate 82 83 Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver. 84 85 .seealso: SNESNewtonTRPreCheck(), SNESNewtonTRGetPreCheck(), SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck() 86 @*/ 87 PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx) 88 { 89 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 90 91 PetscFunctionBegin; 92 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 93 if (func) tr->precheck = func; 94 if (ctx) tr->precheckctx = ctx; 95 PetscFunctionReturn(0); 96 } 97 98 /*@C 99 SNESNewtonTRGetPreCheck - Gets the pre-check function 100 101 Not collective 102 103 Input Parameter: 104 . snes - the nonlinear solver context 105 106 Output Parameters: 107 + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck() 108 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 109 110 Level: intermediate 111 112 .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRPreCheck() 113 @*/ 114 PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx) 115 { 116 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 117 118 PetscFunctionBegin; 119 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 120 if (func) *func = tr->precheck; 121 if (ctx) *ctx = tr->precheckctx; 122 PetscFunctionReturn(0); 123 } 124 125 /*@C 126 SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next 127 function evaluation. Allows the user a chance to change or override the decision of the line search routine 128 129 Logically Collective on snes 130 131 Input Parameters: 132 + snes - the nonlinear solver object 133 . func - [optional] function evaluation routine, see SNESNewtonTRPostCheck() for the calling sequence 134 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 135 136 Level: intermediate 137 138 Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in 139 SNESLineSearchSetPostCheck() is called AFTER the function evaluation. 140 141 .seealso: SNESNewtonTRPostCheck(), SNESNewtonTRGetPostCheck() 142 @*/ 143 PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx) 144 { 145 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 146 147 PetscFunctionBegin; 148 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 149 if (func) tr->postcheck = func; 150 if (ctx) tr->postcheckctx = ctx; 151 PetscFunctionReturn(0); 152 } 153 154 /*@C 155 SNESNewtonTRGetPostCheck - Gets the post-check function 156 157 Not collective 158 159 Input Parameter: 160 . snes - the nonlinear solver context 161 162 Output Parameters: 163 + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck() 164 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 165 166 Level: intermediate 167 168 .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRPostCheck() 169 @*/ 170 PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx) 171 { 172 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 173 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 176 if (func) *func = tr->postcheck; 177 if (ctx) *ctx = tr->postcheckctx; 178 PetscFunctionReturn(0); 179 } 180 181 /*@C 182 SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR 183 184 Logically Collective on snes 185 186 Input Parameters: 187 + snes - the solver 188 . X - The last solution 189 - Y - The step direction 190 191 Output Parameters: 192 . changed_Y - Indicator that the step direction Y has been changed. 193 194 Level: developer 195 196 .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRGetPreCheck() 197 @*/ 198 static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y) 199 { 200 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 *changed_Y = PETSC_FALSE; 205 if (tr->precheck) { 206 ierr = (*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx);CHKERRQ(ierr); 207 PetscValidLogicalCollectiveBool(snes,*changed_Y,4); 208 } 209 PetscFunctionReturn(0); 210 } 211 212 /*@C 213 SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation 214 215 Logically Collective on snes 216 217 Input Parameters: 218 + snes - the solver 219 . X - The last solution 220 . Y - The full step direction 221 - W - The updated solution, W = X - Y 222 223 Output Parameters: 224 + changed_Y - indicator if step has been changed 225 - changed_W - Indicator if the new candidate solution W has been changed. 226 227 Notes: 228 If Y is changed then W is recomputed as X - Y 229 230 Level: developer 231 232 .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck() 233 @*/ 234 static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W) 235 { 236 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 *changed_Y = PETSC_FALSE; 241 *changed_W = PETSC_FALSE; 242 if (tr->postcheck) { 243 ierr = (*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx);CHKERRQ(ierr); 244 PetscValidLogicalCollectiveBool(snes,*changed_Y,5); 245 PetscValidLogicalCollectiveBool(snes,*changed_W,6); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 /* 251 SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust 252 region approach for solving systems of nonlinear equations. 253 254 */ 255 static PetscErrorCode SNESSolve_NEWTONTR(SNES snes) 256 { 257 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 258 Vec X,F,Y,G,Ytmp,W; 259 PetscErrorCode ierr; 260 PetscInt maxits,i,lits; 261 PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 262 PetscScalar cnorm; 263 KSP ksp; 264 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 265 PetscBool breakout = PETSC_FALSE; 266 SNES_TR_KSPConverged_Ctx *ctx; 267 PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*); 268 void *convctx; 269 270 PetscFunctionBegin; 271 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 272 273 maxits = snes->max_its; /* maximum number of iterations */ 274 X = snes->vec_sol; /* solution vector */ 275 F = snes->vec_func; /* residual vector */ 276 Y = snes->work[0]; /* work vectors */ 277 G = snes->work[1]; 278 Ytmp = snes->work[2]; 279 W = snes->work[3]; 280 281 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 282 snes->iter = 0; 283 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 284 285 /* Set the linear stopping criteria to use the More' trick. */ 286 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 287 ierr = KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy);CHKERRQ(ierr); 288 if (convtest != SNESTR_KSPConverged_Private) { 289 ierr = PetscNew(&ctx);CHKERRQ(ierr); 290 ctx->snes = snes; 291 ierr = KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);CHKERRQ(ierr); 292 ierr = KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);CHKERRQ(ierr); 293 ierr = PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");CHKERRQ(ierr); 294 } 295 296 if (!snes->vec_func_init_set) { 297 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 298 } else snes->vec_func_init_set = PETSC_FALSE; 299 300 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 301 SNESCheckFunctionNorm(snes,fnorm); 302 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 303 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 304 snes->norm = fnorm; 305 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 306 delta = xnorm ? neP->delta0*xnorm : neP->delta0; 307 neP->delta = delta; 308 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 309 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 310 311 /* test convergence */ 312 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 313 if (snes->reason) PetscFunctionReturn(0); 314 315 for (i=0; i<maxits; i++) { 316 317 /* Call general purpose update function */ 318 if (snes->ops->update) { 319 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 320 } 321 322 /* Solve J Y = F, where J is Jacobian matrix */ 323 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 324 SNESCheckJacobianDomainerror(snes); 325 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 326 ierr = KSPSolve(snes->ksp,F,Ytmp);CHKERRQ(ierr); 327 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 328 329 snes->linear_its += lits; 330 331 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 332 ierr = VecNorm(Ytmp,NORM_2,&nrm);CHKERRQ(ierr); 333 norm1 = nrm; 334 while (1) { 335 PetscBool changed_y; 336 PetscBool changed_w; 337 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 338 nrm = norm1; 339 340 /* Scale Y if need be and predict new value of F norm */ 341 if (nrm >= delta) { 342 nrm = delta/nrm; 343 gpnorm = (1.0 - nrm)*fnorm; 344 cnorm = nrm; 345 ierr = PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);CHKERRQ(ierr); 346 ierr = VecScale(Y,cnorm);CHKERRQ(ierr); 347 nrm = gpnorm; 348 ynorm = delta; 349 } else { 350 gpnorm = 0.0; 351 ierr = PetscInfo(snes,"Direction is in Trust Region\n");CHKERRQ(ierr); 352 ynorm = nrm; 353 } 354 /* PreCheck() allows for updates to Y prior to W <- X - Y */ 355 ierr = SNESNewtonTRPreCheck(snes,X,Y,&changed_y);CHKERRQ(ierr); 356 ierr = VecWAXPY(W,-1.0,Y,X);CHKERRQ(ierr); /* W <- X - Y */ 357 ierr = SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr); 358 if (changed_y) {ierr = VecWAXPY(W,-1.0,Y,X);CHKERRQ(ierr);} 359 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 360 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); /* F(X) */ 361 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 362 SNESCheckFunctionNorm(snes,gnorm); 363 if (fnorm == gpnorm) rho = 0.0; 364 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 365 366 /* Update size of trust region */ 367 if (rho < neP->mu) delta *= neP->delta1; 368 else if (rho < neP->eta) delta *= neP->delta2; 369 else delta *= neP->delta3; 370 ierr = PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);CHKERRQ(ierr); 371 ierr = PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);CHKERRQ(ierr); 372 373 neP->delta = delta; 374 if (rho > neP->sigma) break; 375 ierr = PetscInfo(snes,"Trying again in smaller region\n");CHKERRQ(ierr); 376 /* check to see if progress is hopeless */ 377 neP->itflag = PETSC_FALSE; 378 ierr = SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 379 if (!reason) {ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr);} 380 if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER; 381 if (reason) { 382 /* We're not progressing, so return with the current iterate */ 383 ierr = SNESMonitor(snes,i+1,fnorm);CHKERRQ(ierr); 384 breakout = PETSC_TRUE; 385 break; 386 } 387 snes->numFailures++; 388 } 389 if (!breakout) { 390 /* Update function and solution vectors */ 391 fnorm = gnorm; 392 ierr = VecCopy(G,F);CHKERRQ(ierr); 393 ierr = VecCopy(W,X);CHKERRQ(ierr); 394 /* Monitor convergence */ 395 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 396 snes->iter = i+1; 397 snes->norm = fnorm; 398 snes->xnorm = xnorm; 399 snes->ynorm = ynorm; 400 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 401 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 402 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 403 /* Test for convergence, xnorm = || X || */ 404 neP->itflag = PETSC_TRUE; 405 if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 406 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 407 if (reason) break; 408 } else break; 409 } 410 if (i == maxits) { 411 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 412 if (!reason) reason = SNES_DIVERGED_MAX_IT; 413 } 414 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 415 snes->reason = reason; 416 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 417 if (convtest != SNESTR_KSPConverged_Private) { 418 ierr = KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);CHKERRQ(ierr); 419 ierr = PetscFree(ctx);CHKERRQ(ierr); 420 ierr = KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy);CHKERRQ(ierr); 421 } 422 PetscFunctionReturn(0); 423 } 424 /*------------------------------------------------------------*/ 425 static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) 426 { 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 ierr = SNESSetWorkVecs(snes,4);CHKERRQ(ierr); 431 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 435 PetscErrorCode SNESReset_NEWTONTR(SNES snes) 436 { 437 438 PetscFunctionBegin; 439 PetscFunctionReturn(0); 440 } 441 442 static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) 443 { 444 PetscErrorCode ierr; 445 446 PetscFunctionBegin; 447 ierr = SNESReset_NEWTONTR(snes);CHKERRQ(ierr); 448 ierr = PetscFree(snes->data);CHKERRQ(ierr); 449 PetscFunctionReturn(0); 450 } 451 /*------------------------------------------------------------*/ 452 453 static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes) 454 { 455 SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data; 456 PetscErrorCode ierr; 457 458 PetscFunctionBegin; 459 ierr = PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");CHKERRQ(ierr); 460 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);CHKERRQ(ierr); 461 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);CHKERRQ(ierr); 462 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);CHKERRQ(ierr); 463 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);CHKERRQ(ierr); 464 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);CHKERRQ(ierr); 465 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);CHKERRQ(ierr); 466 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);CHKERRQ(ierr); 467 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);CHKERRQ(ierr); 468 ierr = PetscOptionsTail();CHKERRQ(ierr); 469 PetscFunctionReturn(0); 470 } 471 472 static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer) 473 { 474 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 475 PetscErrorCode ierr; 476 PetscBool iascii; 477 478 PetscFunctionBegin; 479 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 480 if (iascii) { 481 ierr = PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);CHKERRQ(ierr); 482 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);CHKERRQ(ierr); 483 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);CHKERRQ(ierr); 484 } 485 PetscFunctionReturn(0); 486 } 487 /* ------------------------------------------------------------ */ 488 /*MC 489 SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 490 491 Options Database: 492 + -snes_trtol <tol> - trust region tolerance 493 . -snes_tr_mu <mu> - trust region parameter 494 . -snes_tr_eta <eta> - trust region parameter 495 . -snes_tr_sigma <sigma> - trust region parameter 496 . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 497 . -snes_tr_delta1 <delta1> - trust region parameter 498 . -snes_tr_delta2 <delta2> - trust region parameter 499 - -snes_tr_delta3 <delta3> - trust region parameter 500 501 The basic algorithm is taken from "The Minpack Project", by More', 502 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 503 of Mathematical Software", Wayne Cowell, editor. 504 505 Level: intermediate 506 507 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance() 508 509 M*/ 510 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) 511 { 512 SNES_NEWTONTR *neP; 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 snes->ops->setup = SNESSetUp_NEWTONTR; 517 snes->ops->solve = SNESSolve_NEWTONTR; 518 snes->ops->destroy = SNESDestroy_NEWTONTR; 519 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 520 snes->ops->view = SNESView_NEWTONTR; 521 snes->ops->reset = SNESReset_NEWTONTR; 522 523 snes->usesksp = PETSC_TRUE; 524 snes->usesnpc = PETSC_FALSE; 525 526 snes->alwayscomputesfinalresidual = PETSC_TRUE; 527 528 ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr); 529 snes->data = (void*)neP; 530 neP->mu = 0.25; 531 neP->eta = 0.75; 532 neP->delta = 0.0; 533 neP->delta0 = 0.2; 534 neP->delta1 = 0.3; 535 neP->delta2 = 0.75; 536 neP->delta3 = 2.0; 537 neP->sigma = 0.0001; 538 neP->itflag = PETSC_FALSE; 539 neP->rnorm0 = 0.0; 540 neP->ttol = 0.0; 541 PetscFunctionReturn(0); 542 } 543