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) { 229566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Default or user provided convergence test KSP iterations=%D, 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) { 609566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum number of function evaluations: %D\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 813312a946SBarry Smith .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 108c9368356SGlenn Hammond .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 1377cb011f5SBarry Smith .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 1647cb011f5SBarry Smith .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 192c9368356SGlenn Hammond .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 2277cb011f5SBarry Smith .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; 2642c71b3e2SJacob Faibussowitsch PetscCheckFalse(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; 449*d0609cedSBarry 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)); 458*d0609cedSBarry 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) { 4709566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," Trust region tolerance (-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 49604d7464bSBarry Smith .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