1c6db04a5SJed Brown #include <../src/snes/impls/tr/trimpl.h> /*I "petscsnes.h" I*/ 24800dd8cSBarry Smith 3971273eeSBarry Smith typedef struct { 4971273eeSBarry Smith SNES snes; 5df8705c3SBarry Smith /* Information on the regular SNES convergence test; which may have been user provided */ 6df8705c3SBarry Smith PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 7df8705c3SBarry Smith PetscErrorCode (*convdestroy)(void*); 8df8705c3SBarry Smith void *convctx; 9971273eeSBarry Smith } SNES_TR_KSPConverged_Ctx; 10971273eeSBarry Smith 11470880abSPatrick Sanan static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx) 12df60cc22SBarry Smith { 13971273eeSBarry Smith SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 14971273eeSBarry Smith SNES snes = ctx->snes; 1504d7464bSBarry Smith SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 16df60cc22SBarry Smith Vec x; 17064f8208SBarry Smith PetscReal nrm; 18df60cc22SBarry Smith 193a40ed3dSBarry Smith PetscFunctionBegin; 209566063dSJacob Faibussowitsch PetscCall((*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx)); 21329f5518SBarry Smith if (*reason) { 2263a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n",n,(double)rnorm)); 239b04593eSLois Curfman McInnes } 24a935fc98SLois Curfman McInnes /* Determine norm of solution */ 259566063dSJacob Faibussowitsch PetscCall(KSPBuildSolution(ksp,NULL,&x)); 269566063dSJacob Faibussowitsch PetscCall(VecNorm(x,NORM_2,&nrm)); 27064f8208SBarry Smith if (nrm >= neP->delta) { 289566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm)); 29329f5518SBarry Smith *reason = KSP_CONVERGED_STEP_LENGTH; 30df60cc22SBarry Smith } 313a40ed3dSBarry Smith PetscFunctionReturn(0); 32df60cc22SBarry Smith } 3382bf6240SBarry Smith 34470880abSPatrick Sanan static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx) 35971273eeSBarry Smith { 36971273eeSBarry Smith SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 37971273eeSBarry Smith 38971273eeSBarry Smith PetscFunctionBegin; 399566063dSJacob Faibussowitsch PetscCall((*ctx->convdestroy)(ctx->convctx)); 409566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 41971273eeSBarry Smith PetscFunctionReturn(0); 42971273eeSBarry Smith } 43971273eeSBarry Smith 4485385478SLisandro Dalcin /* ---------------------------------------------------------------- */ 4585385478SLisandro Dalcin /* 46470880abSPatrick Sanan SNESTR_Converged_Private -test convergence JUST for 4785385478SLisandro Dalcin the trust region tolerance. 4885385478SLisandro Dalcin 4985385478SLisandro Dalcin */ 50470880abSPatrick Sanan static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 5185385478SLisandro Dalcin { 5204d7464bSBarry Smith SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 5385385478SLisandro Dalcin 5485385478SLisandro Dalcin PetscFunctionBegin; 5585385478SLisandro Dalcin *reason = SNES_CONVERGED_ITERATING; 5685385478SLisandro Dalcin if (neP->delta < xnorm * snes->deltatol) { 579566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol)); 581c6b2ff8SBarry Smith *reason = SNES_DIVERGED_TR_DELTA; 59e71169deSBarry Smith } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 6063a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n",snes->max_funcs)); 6185385478SLisandro Dalcin *reason = SNES_DIVERGED_FUNCTION_COUNT; 6285385478SLisandro Dalcin } 6385385478SLisandro Dalcin PetscFunctionReturn(0); 6485385478SLisandro Dalcin } 6585385478SLisandro Dalcin 667cb011f5SBarry Smith /*@C 67c9368356SGlenn Hammond SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined. 68c9368356SGlenn Hammond Allows the user a chance to change or override the decision of the line search routine. 69c9368356SGlenn Hammond 70c9368356SGlenn Hammond Logically Collective on snes 71c9368356SGlenn Hammond 72c9368356SGlenn Hammond Input Parameters: 73c9368356SGlenn Hammond + snes - the nonlinear solver object 74c9368356SGlenn Hammond . func - [optional] function evaluation routine, see SNESNewtonTRPreCheck() for the calling sequence 75c9368356SGlenn Hammond - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 76c9368356SGlenn Hammond 77c9368356SGlenn Hammond Level: intermediate 78c9368356SGlenn Hammond 79c9368356SGlenn Hammond Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver. 80c9368356SGlenn Hammond 81*db781477SPatrick Sanan .seealso: `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()` 82c9368356SGlenn Hammond @*/ 83c9368356SGlenn Hammond PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx) 84c9368356SGlenn Hammond { 85c9368356SGlenn Hammond SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 86c9368356SGlenn Hammond 87c9368356SGlenn Hammond PetscFunctionBegin; 88c9368356SGlenn Hammond PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 89c9368356SGlenn Hammond if (func) tr->precheck = func; 90c9368356SGlenn Hammond if (ctx) tr->precheckctx = ctx; 91c9368356SGlenn Hammond PetscFunctionReturn(0); 92c9368356SGlenn Hammond } 93c9368356SGlenn Hammond 94c9368356SGlenn Hammond /*@C 95c9368356SGlenn Hammond SNESNewtonTRGetPreCheck - Gets the pre-check function 96c9368356SGlenn Hammond 97c9368356SGlenn Hammond Not collective 98c9368356SGlenn Hammond 99c9368356SGlenn Hammond Input Parameter: 100c9368356SGlenn Hammond . snes - the nonlinear solver context 101c9368356SGlenn Hammond 102c9368356SGlenn Hammond Output Parameters: 103c9368356SGlenn Hammond + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck() 104c9368356SGlenn Hammond - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 105c9368356SGlenn Hammond 106c9368356SGlenn Hammond Level: intermediate 107c9368356SGlenn Hammond 108*db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()` 109c9368356SGlenn Hammond @*/ 110c9368356SGlenn Hammond PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx) 111c9368356SGlenn Hammond { 112c9368356SGlenn Hammond SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 113c9368356SGlenn Hammond 114c9368356SGlenn Hammond PetscFunctionBegin; 115c9368356SGlenn Hammond PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 116c9368356SGlenn Hammond if (func) *func = tr->precheck; 117c9368356SGlenn Hammond if (ctx) *ctx = tr->precheckctx; 118c9368356SGlenn Hammond PetscFunctionReturn(0); 119c9368356SGlenn Hammond } 120c9368356SGlenn Hammond 121c9368356SGlenn Hammond /*@C 1227cb011f5SBarry Smith SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next 1237cb011f5SBarry Smith function evaluation. Allows the user a chance to change or override the decision of the line search routine 1247cb011f5SBarry Smith 1257cb011f5SBarry Smith Logically Collective on snes 1267cb011f5SBarry Smith 1277cb011f5SBarry Smith Input Parameters: 1287cb011f5SBarry Smith + snes - the nonlinear solver object 1297cb011f5SBarry Smith . func - [optional] function evaluation routine, see SNESNewtonTRPostCheck() for the calling sequence 1307cb011f5SBarry Smith - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 1317cb011f5SBarry Smith 1327cb011f5SBarry Smith Level: intermediate 1337cb011f5SBarry Smith 134c9368356SGlenn Hammond Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in 1357cb011f5SBarry Smith SNESLineSearchSetPostCheck() is called AFTER the function evaluation. 1367cb011f5SBarry Smith 137*db781477SPatrick Sanan .seealso: `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()` 1387cb011f5SBarry Smith @*/ 139c9368356SGlenn Hammond PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx) 1407cb011f5SBarry Smith { 1417cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 1427cb011f5SBarry Smith 1437cb011f5SBarry Smith PetscFunctionBegin; 1447cb011f5SBarry Smith PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1457cb011f5SBarry Smith if (func) tr->postcheck = func; 1467cb011f5SBarry Smith if (ctx) tr->postcheckctx = ctx; 1477cb011f5SBarry Smith PetscFunctionReturn(0); 1487cb011f5SBarry Smith } 1497cb011f5SBarry Smith 1507cb011f5SBarry Smith /*@C 1517cb011f5SBarry Smith SNESNewtonTRGetPostCheck - Gets the post-check function 1527cb011f5SBarry Smith 1537cb011f5SBarry Smith Not collective 1547cb011f5SBarry Smith 1557cb011f5SBarry Smith Input Parameter: 1567cb011f5SBarry Smith . snes - the nonlinear solver context 1577cb011f5SBarry Smith 1587cb011f5SBarry Smith Output Parameters: 1597cb011f5SBarry Smith + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck() 1607cb011f5SBarry Smith - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 1617cb011f5SBarry Smith 1627cb011f5SBarry Smith Level: intermediate 1637cb011f5SBarry Smith 164*db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()` 1657cb011f5SBarry Smith @*/ 166c9368356SGlenn Hammond PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx) 1677cb011f5SBarry Smith { 1687cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 1697cb011f5SBarry Smith 1707cb011f5SBarry Smith PetscFunctionBegin; 1717cb011f5SBarry Smith PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1727cb011f5SBarry Smith if (func) *func = tr->postcheck; 1737cb011f5SBarry Smith if (ctx) *ctx = tr->postcheckctx; 1747cb011f5SBarry Smith PetscFunctionReturn(0); 1757cb011f5SBarry Smith } 1767cb011f5SBarry Smith 1777cb011f5SBarry Smith /*@C 178c9368356SGlenn Hammond SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR 179c9368356SGlenn Hammond 180c9368356SGlenn Hammond Logically Collective on snes 181c9368356SGlenn Hammond 182c9368356SGlenn Hammond Input Parameters: 183c9368356SGlenn Hammond + snes - the solver 184c9368356SGlenn Hammond . X - The last solution 185c9368356SGlenn Hammond - Y - The step direction 186c9368356SGlenn Hammond 187c9368356SGlenn Hammond Output Parameters: 188c9368356SGlenn Hammond . changed_Y - Indicator that the step direction Y has been changed. 189c9368356SGlenn Hammond 190c9368356SGlenn Hammond Level: developer 191c9368356SGlenn Hammond 192*db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()` 193c9368356SGlenn Hammond @*/ 194c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y) 195c9368356SGlenn Hammond { 196c9368356SGlenn Hammond SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 197c9368356SGlenn Hammond 198c9368356SGlenn Hammond PetscFunctionBegin; 199c9368356SGlenn Hammond *changed_Y = PETSC_FALSE; 200c9368356SGlenn Hammond if (tr->precheck) { 2019566063dSJacob Faibussowitsch PetscCall((*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx)); 202c9368356SGlenn Hammond PetscValidLogicalCollectiveBool(snes,*changed_Y,4); 203c9368356SGlenn Hammond } 204c9368356SGlenn Hammond PetscFunctionReturn(0); 205c9368356SGlenn Hammond } 206c9368356SGlenn Hammond 207c9368356SGlenn Hammond /*@C 2087cb011f5SBarry Smith SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation 2097cb011f5SBarry Smith 2107cb011f5SBarry Smith Logically Collective on snes 2117cb011f5SBarry Smith 2127cb011f5SBarry Smith Input Parameters: 2136b867d5aSJose E. Roman + snes - the solver 2146b867d5aSJose E. Roman . X - The last solution 2157cb011f5SBarry Smith . Y - The full step direction 2163312a946SBarry Smith - W - The updated solution, W = X - Y 2177cb011f5SBarry Smith 2187cb011f5SBarry Smith Output Parameters: 2193312a946SBarry Smith + changed_Y - indicator if step has been changed 2203312a946SBarry Smith - changed_W - Indicator if the new candidate solution W has been changed. 2217cb011f5SBarry Smith 2227cb011f5SBarry Smith Notes: 2233312a946SBarry Smith If Y is changed then W is recomputed as X - Y 2247cb011f5SBarry Smith 2257cb011f5SBarry Smith Level: developer 2267cb011f5SBarry Smith 227*db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()` 2287cb011f5SBarry Smith @*/ 229c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W) 2307cb011f5SBarry Smith { 2317cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 2327cb011f5SBarry Smith 2337cb011f5SBarry Smith PetscFunctionBegin; 234c9368356SGlenn Hammond *changed_Y = PETSC_FALSE; 2357cb011f5SBarry Smith *changed_W = PETSC_FALSE; 2367cb011f5SBarry Smith if (tr->postcheck) { 2379566063dSJacob Faibussowitsch PetscCall((*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx)); 238c9368356SGlenn Hammond PetscValidLogicalCollectiveBool(snes,*changed_Y,5); 2397cb011f5SBarry Smith PetscValidLogicalCollectiveBool(snes,*changed_W,6); 2407cb011f5SBarry Smith } 2417cb011f5SBarry Smith PetscFunctionReturn(0); 2427cb011f5SBarry Smith } 24385385478SLisandro Dalcin 244df60cc22SBarry Smith /* 24504d7464bSBarry Smith SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust 246ddbbdb52SLois Curfman McInnes region approach for solving systems of nonlinear equations. 2474800dd8cSBarry Smith 2484800dd8cSBarry Smith */ 24904d7464bSBarry Smith static PetscErrorCode SNESSolve_NEWTONTR(SNES snes) 2504800dd8cSBarry Smith { 25104d7464bSBarry Smith SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 252c9368356SGlenn Hammond Vec X,F,Y,G,Ytmp,W; 253a7cc72afSBarry Smith PetscInt maxits,i,lits; 25485385478SLisandro Dalcin PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 25579f36460SBarry Smith PetscScalar cnorm; 256df60cc22SBarry Smith KSP ksp; 2573f149594SLisandro Dalcin SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 258df8705c3SBarry Smith PetscBool breakout = PETSC_FALSE; 259df8705c3SBarry Smith SNES_TR_KSPConverged_Ctx *ctx; 2605e28bcb6Sprj- PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*); 2615e28bcb6Sprj- void *convctx; 2624800dd8cSBarry Smith 2633a40ed3dSBarry Smith PetscFunctionBegin; 2640b121fc5SBarry Smith PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 265c579b300SPatrick Farrell 266fbe28522SBarry Smith maxits = snes->max_its; /* maximum number of iterations */ 267fbe28522SBarry Smith X = snes->vec_sol; /* solution vector */ 26839e2f89bSBarry Smith F = snes->vec_func; /* residual vector */ 269fbe28522SBarry Smith Y = snes->work[0]; /* work vectors */ 270fbe28522SBarry Smith G = snes->work[1]; 2716b5873e3SBarry Smith Ytmp = snes->work[2]; 272c9368356SGlenn Hammond W = snes->work[3]; 2734800dd8cSBarry Smith 2749566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 2754c49b128SBarry Smith snes->iter = 0; 2769566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 2774800dd8cSBarry Smith 278df8705c3SBarry Smith /* Set the linear stopping criteria to use the More' trick. */ 2799566063dSJacob Faibussowitsch PetscCall(SNESGetKSP(snes,&ksp)); 2809566063dSJacob Faibussowitsch PetscCall(KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy)); 281df8705c3SBarry Smith if (convtest != SNESTR_KSPConverged_Private) { 2829566063dSJacob Faibussowitsch PetscCall(PetscNew(&ctx)); 283df8705c3SBarry Smith ctx->snes = snes; 2849566063dSJacob Faibussowitsch PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy)); 2859566063dSJacob Faibussowitsch PetscCall(KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy)); 2869566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n")); 287df8705c3SBarry Smith } 288df8705c3SBarry Smith 289e4ed7901SPeter Brune if (!snes->vec_func_init_set) { 2909566063dSJacob Faibussowitsch PetscCall(SNESComputeFunction(snes,X,F)); /* F(X) */ 2911aa26658SKarl Rupp } else snes->vec_func_init_set = PETSC_FALSE; 292e4ed7901SPeter Brune 2939566063dSJacob Faibussowitsch PetscCall(VecNorm(F,NORM_2,&fnorm)); /* fnorm <- || F || */ 294422a814eSBarry Smith SNESCheckFunctionNorm(snes,fnorm); 2959566063dSJacob Faibussowitsch PetscCall(VecNorm(X,NORM_2,&xnorm)); /* xnorm <- || X || */ 2969566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 297fbe28522SBarry Smith snes->norm = fnorm; 2989566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 2998b84755bSBarry Smith delta = xnorm ? neP->delta0*xnorm : neP->delta0; 3004800dd8cSBarry Smith neP->delta = delta; 3019566063dSJacob Faibussowitsch PetscCall(SNESLogConvergenceHistory(snes,fnorm,0)); 3029566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes,0,fnorm)); 303b37302c6SLois Curfman McInnes 30485385478SLisandro Dalcin /* test convergence */ 3059566063dSJacob Faibussowitsch PetscCall((*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP)); 30685385478SLisandro Dalcin if (snes->reason) PetscFunctionReturn(0); 3073f149594SLisandro Dalcin 3084800dd8cSBarry Smith for (i=0; i<maxits; i++) { 30999a96b7cSMatthew Knepley 31099a96b7cSMatthew Knepley /* Call general purpose update function */ 311e7788613SBarry Smith if (snes->ops->update) { 3129566063dSJacob Faibussowitsch PetscCall((*snes->ops->update)(snes, snes->iter)); 31399a96b7cSMatthew Knepley } 31499a96b7cSMatthew Knepley 31585385478SLisandro Dalcin /* Solve J Y = F, where J is Jacobian matrix */ 3169566063dSJacob Faibussowitsch PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 31707b62357SFande Kong SNESCheckJacobianDomainerror(snes); 3189566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre)); 3199566063dSJacob Faibussowitsch PetscCall(KSPSolve(snes->ksp,F,Ytmp)); 3209566063dSJacob Faibussowitsch PetscCall(KSPGetIterationNumber(snes->ksp,&lits)); 321329f5518SBarry Smith snes->linear_its += lits; 3221aa26658SKarl Rupp 3239566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits)); 3249566063dSJacob Faibussowitsch PetscCall(VecNorm(Ytmp,NORM_2,&nrm)); 325064f8208SBarry Smith norm1 = nrm; 326284fb49fSHeeho Park 3276b5873e3SBarry Smith while (1) { 328c9368356SGlenn Hammond PetscBool changed_y; 3297cb011f5SBarry Smith PetscBool changed_w; 3309566063dSJacob Faibussowitsch PetscCall(VecCopy(Ytmp,Y)); 331064f8208SBarry Smith nrm = norm1; 3326b5873e3SBarry Smith 3332deea200SLois Curfman McInnes /* Scale Y if need be and predict new value of F norm */ 334064f8208SBarry Smith if (nrm >= delta) { 335064f8208SBarry Smith nrm = delta/nrm; 336064f8208SBarry Smith gpnorm = (1.0 - nrm)*fnorm; 337064f8208SBarry Smith cnorm = nrm; 3389566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Scaling direction by %g\n",(double)nrm)); 3399566063dSJacob Faibussowitsch PetscCall(VecScale(Y,cnorm)); 340064f8208SBarry Smith nrm = gpnorm; 341fbe28522SBarry Smith ynorm = delta; 342fbe28522SBarry Smith } else { 343fbe28522SBarry Smith gpnorm = 0.0; 3449566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Direction is in Trust Region\n")); 345064f8208SBarry Smith ynorm = nrm; 346fbe28522SBarry Smith } 347c9368356SGlenn Hammond /* PreCheck() allows for updates to Y prior to W <- X - Y */ 3489566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPreCheck(snes,X,Y,&changed_y)); 3499566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-1.0,Y,X)); /* W <- X - Y */ 3509566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w)); 3519566063dSJacob Faibussowitsch if (changed_y) PetscCall(VecWAXPY(W,-1.0,Y,X)); 3529566063dSJacob Faibussowitsch PetscCall(VecCopy(Y,snes->vec_sol_update)); 3539566063dSJacob Faibussowitsch PetscCall(SNESComputeFunction(snes,W,G)); /* F(X-Y) = G */ 3549566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); /* gnorm <- || g || */ 3557551ef16SBarry Smith SNESCheckFunctionNorm(snes,gnorm); 3564800dd8cSBarry Smith if (fnorm == gpnorm) rho = 0.0; 357fbe28522SBarry Smith else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 3584800dd8cSBarry Smith 3594800dd8cSBarry Smith /* Update size of trust region */ 3604800dd8cSBarry Smith if (rho < neP->mu) delta *= neP->delta1; 3614800dd8cSBarry Smith else if (rho < neP->eta) delta *= neP->delta2; 3624800dd8cSBarry Smith else delta *= neP->delta3; 3639566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm)); 3649566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta)); 3651aa26658SKarl Rupp 3664800dd8cSBarry Smith neP->delta = delta; 3674800dd8cSBarry Smith if (rho > neP->sigma) break; 3689566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Trying again in smaller region\n")); 369284fb49fSHeeho Park 3706b5873e3SBarry Smith /* check to see if progress is hopeless */ 371ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 3729566063dSJacob Faibussowitsch PetscCall(SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP)); 3739566063dSJacob Faibussowitsch if (!reason) PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP)); 3747cb011f5SBarry Smith if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER; 375184914b5SBarry Smith if (reason) { 37652392280SLois Curfman McInnes /* We're not progressing, so return with the current iterate */ 3779566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes,i+1,fnorm)); 378a7cc72afSBarry Smith breakout = PETSC_TRUE; 379454a90a3SBarry Smith break; 38052392280SLois Curfman McInnes } 38150ffb88aSMatthew Knepley snes->numFailures++; 3824800dd8cSBarry Smith } 3831acabf8cSLois Curfman McInnes if (!breakout) { 38485385478SLisandro Dalcin /* Update function and solution vectors */ 3854800dd8cSBarry Smith fnorm = gnorm; 3869566063dSJacob Faibussowitsch PetscCall(VecCopy(G,F)); 3879566063dSJacob Faibussowitsch PetscCall(VecCopy(W,X)); 38885385478SLisandro Dalcin /* Monitor convergence */ 3899566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 390c509a162SBarry Smith snes->iter = i+1; 391fbe28522SBarry Smith snes->norm = fnorm; 392c1e67a49SFande Kong snes->xnorm = xnorm; 393c1e67a49SFande Kong snes->ynorm = ynorm; 3949566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 3959566063dSJacob Faibussowitsch PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits)); 3969566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes,snes->iter,snes->norm)); 39785385478SLisandro Dalcin /* Test for convergence, xnorm = || X || */ 398ef87ac0dSKris Buschelman neP->itflag = PETSC_TRUE; 3999566063dSJacob Faibussowitsch if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm)); 4009566063dSJacob Faibussowitsch PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP)); 4013f149594SLisandro Dalcin if (reason) break; 4021aa26658SKarl Rupp } else break; 40338442cffSBarry Smith } 404284fb49fSHeeho Park 40552392280SLois Curfman McInnes if (i == maxits) { 4069566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits)); 40785385478SLisandro Dalcin if (!reason) reason = SNES_DIVERGED_MAX_IT; 40852392280SLois Curfman McInnes } 4099566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 410184914b5SBarry Smith snes->reason = reason; 4119566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 4125e28bcb6Sprj- if (convtest != SNESTR_KSPConverged_Private) { 4139566063dSJacob Faibussowitsch PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy)); 4149566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 4159566063dSJacob Faibussowitsch PetscCall(KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy)); 4165e28bcb6Sprj- } 4173a40ed3dSBarry Smith PetscFunctionReturn(0); 4184800dd8cSBarry Smith } 419284fb49fSHeeho Park 4204800dd8cSBarry Smith /*------------------------------------------------------------*/ 42104d7464bSBarry Smith static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) 4224800dd8cSBarry Smith { 4233a40ed3dSBarry Smith PetscFunctionBegin; 4249566063dSJacob Faibussowitsch PetscCall(SNESSetWorkVecs(snes,4)); 4259566063dSJacob Faibussowitsch PetscCall(SNESSetUpMatrices(snes)); 4263a40ed3dSBarry Smith PetscFunctionReturn(0); 4274800dd8cSBarry Smith } 4286b8b9a38SLisandro Dalcin 42904d7464bSBarry Smith PetscErrorCode SNESReset_NEWTONTR(SNES snes) 4306b8b9a38SLisandro Dalcin { 4316b8b9a38SLisandro Dalcin PetscFunctionBegin; 4326b8b9a38SLisandro Dalcin PetscFunctionReturn(0); 4336b8b9a38SLisandro Dalcin } 4346b8b9a38SLisandro Dalcin 43504d7464bSBarry Smith static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) 4364800dd8cSBarry Smith { 4373a40ed3dSBarry Smith PetscFunctionBegin; 4389566063dSJacob Faibussowitsch PetscCall(SNESReset_NEWTONTR(snes)); 4399566063dSJacob Faibussowitsch PetscCall(PetscFree(snes->data)); 4403a40ed3dSBarry Smith PetscFunctionReturn(0); 4414800dd8cSBarry Smith } 4424800dd8cSBarry Smith /*------------------------------------------------------------*/ 4434800dd8cSBarry Smith 4444416b707SBarry Smith static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes) 4454800dd8cSBarry Smith { 44604d7464bSBarry Smith SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data; 4474800dd8cSBarry Smith 4483a40ed3dSBarry Smith PetscFunctionBegin; 449d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject,"SNES trust region options for nonlinear equations"); 4509566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL)); 4519566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL)); 4529566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL)); 4539566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL)); 4549566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL)); 4559566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL)); 4569566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL)); 4579566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL)); 458d0609cedSBarry Smith PetscOptionsHeadEnd(); 4593a40ed3dSBarry Smith PetscFunctionReturn(0); 4604800dd8cSBarry Smith } 4614800dd8cSBarry Smith 46204d7464bSBarry Smith static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer) 463a935fc98SLois Curfman McInnes { 46404d7464bSBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 465ace3abfcSBarry Smith PetscBool iascii; 466a935fc98SLois Curfman McInnes 4673a40ed3dSBarry Smith PetscFunctionBegin; 4689566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 46932077d6dSBarry Smith if (iascii) { 47063a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," Trust region tolerance %g (-snes_trtol)\n",(double)snes->deltatol)); 4719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma)); 4729566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3)); 47319bcc07fSBarry Smith } 4743a40ed3dSBarry Smith PetscFunctionReturn(0); 475a935fc98SLois Curfman McInnes } 47652392280SLois Curfman McInnes /* ------------------------------------------------------------ */ 477ebe3b25bSBarry Smith /*MC 47804d7464bSBarry Smith SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 479ebe3b25bSBarry Smith 480ebe3b25bSBarry Smith Options Database: 4818e434772SBarry Smith + -snes_trtol <tol> - trust region tolerance 4828e434772SBarry Smith . -snes_tr_mu <mu> - trust region parameter 4838e434772SBarry Smith . -snes_tr_eta <eta> - trust region parameter 4848e434772SBarry Smith . -snes_tr_sigma <sigma> - trust region parameter 4858b84755bSBarry Smith . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 4868e434772SBarry Smith . -snes_tr_delta1 <delta1> - trust region parameter 4878e434772SBarry Smith . -snes_tr_delta2 <delta2> - trust region parameter 4888e434772SBarry Smith - -snes_tr_delta3 <delta3> - trust region parameter 489ebe3b25bSBarry Smith 490ebe3b25bSBarry Smith The basic algorithm is taken from "The Minpack Project", by More', 491ebe3b25bSBarry Smith Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 492ebe3b25bSBarry Smith of Mathematical Software", Wayne Cowell, editor. 493ebe3b25bSBarry Smith 494ee3001cbSBarry Smith Level: intermediate 495ee3001cbSBarry Smith 496*db781477SPatrick Sanan .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()` 497ebe3b25bSBarry Smith 498ebe3b25bSBarry Smith M*/ 4998cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) 5004800dd8cSBarry Smith { 50104d7464bSBarry Smith SNES_NEWTONTR *neP; 5024800dd8cSBarry Smith 5033a40ed3dSBarry Smith PetscFunctionBegin; 50404d7464bSBarry Smith snes->ops->setup = SNESSetUp_NEWTONTR; 50504d7464bSBarry Smith snes->ops->solve = SNESSolve_NEWTONTR; 50604d7464bSBarry Smith snes->ops->destroy = SNESDestroy_NEWTONTR; 50704d7464bSBarry Smith snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 50804d7464bSBarry Smith snes->ops->view = SNESView_NEWTONTR; 50904d7464bSBarry Smith snes->ops->reset = SNESReset_NEWTONTR; 510fbe28522SBarry Smith 51142f4f86dSBarry Smith snes->usesksp = PETSC_TRUE; 512efd4aadfSBarry Smith snes->usesnpc = PETSC_FALSE; 51342f4f86dSBarry Smith 5144fc747eaSLawrence Mitchell snes->alwayscomputesfinalresidual = PETSC_TRUE; 5154fc747eaSLawrence Mitchell 5169566063dSJacob Faibussowitsch PetscCall(PetscNewLog(snes,&neP)); 517fbe28522SBarry Smith snes->data = (void*)neP; 518fbe28522SBarry Smith neP->mu = 0.25; 519fbe28522SBarry Smith neP->eta = 0.75; 520fbe28522SBarry Smith neP->delta = 0.0; 521fbe28522SBarry Smith neP->delta0 = 0.2; 522fbe28522SBarry Smith neP->delta1 = 0.3; 523fbe28522SBarry Smith neP->delta2 = 0.75; 524fbe28522SBarry Smith neP->delta3 = 2.0; 525fbe28522SBarry Smith neP->sigma = 0.0001; 526ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 52775567043SBarry Smith neP->rnorm0 = 0.0; 52875567043SBarry Smith neP->ttol = 0.0; 5293a40ed3dSBarry Smith PetscFunctionReturn(0); 5304800dd8cSBarry Smith } 531