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 81db781477SPatrick 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 108db781477SPatrick 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 137db781477SPatrick 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 164db781477SPatrick 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 192db781477SPatrick 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 227db781477SPatrick 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 */ 305*dbbe0bcdSBarry Smith PetscUseTypeMethod(snes,converged ,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 */ 311*dbbe0bcdSBarry Smith PetscTryTypeMethod(snes,update, snes->iter); 31299a96b7cSMatthew Knepley 31385385478SLisandro Dalcin /* Solve J Y = F, where J is Jacobian matrix */ 3149566063dSJacob Faibussowitsch PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 31507b62357SFande Kong SNESCheckJacobianDomainerror(snes); 3169566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre)); 3179566063dSJacob Faibussowitsch PetscCall(KSPSolve(snes->ksp,F,Ytmp)); 3189566063dSJacob Faibussowitsch PetscCall(KSPGetIterationNumber(snes->ksp,&lits)); 319329f5518SBarry Smith snes->linear_its += lits; 3201aa26658SKarl Rupp 3219566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits)); 3229566063dSJacob Faibussowitsch PetscCall(VecNorm(Ytmp,NORM_2,&nrm)); 323064f8208SBarry Smith norm1 = nrm; 324284fb49fSHeeho Park 3256b5873e3SBarry Smith while (1) { 326c9368356SGlenn Hammond PetscBool changed_y; 3277cb011f5SBarry Smith PetscBool changed_w; 3289566063dSJacob Faibussowitsch PetscCall(VecCopy(Ytmp,Y)); 329064f8208SBarry Smith nrm = norm1; 3306b5873e3SBarry Smith 3312deea200SLois Curfman McInnes /* Scale Y if need be and predict new value of F norm */ 332064f8208SBarry Smith if (nrm >= delta) { 333064f8208SBarry Smith nrm = delta/nrm; 334064f8208SBarry Smith gpnorm = (1.0 - nrm)*fnorm; 335064f8208SBarry Smith cnorm = nrm; 3369566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Scaling direction by %g\n",(double)nrm)); 3379566063dSJacob Faibussowitsch PetscCall(VecScale(Y,cnorm)); 338064f8208SBarry Smith nrm = gpnorm; 339fbe28522SBarry Smith ynorm = delta; 340fbe28522SBarry Smith } else { 341fbe28522SBarry Smith gpnorm = 0.0; 3429566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Direction is in Trust Region\n")); 343064f8208SBarry Smith ynorm = nrm; 344fbe28522SBarry Smith } 345c9368356SGlenn Hammond /* PreCheck() allows for updates to Y prior to W <- X - Y */ 3469566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPreCheck(snes,X,Y,&changed_y)); 3479566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-1.0,Y,X)); /* W <- X - Y */ 3489566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w)); 3499566063dSJacob Faibussowitsch if (changed_y) PetscCall(VecWAXPY(W,-1.0,Y,X)); 3509566063dSJacob Faibussowitsch PetscCall(VecCopy(Y,snes->vec_sol_update)); 3519566063dSJacob Faibussowitsch PetscCall(SNESComputeFunction(snes,W,G)); /* F(X-Y) = G */ 3529566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); /* gnorm <- || g || */ 3537551ef16SBarry Smith SNESCheckFunctionNorm(snes,gnorm); 3544800dd8cSBarry Smith if (fnorm == gpnorm) rho = 0.0; 355fbe28522SBarry Smith else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 3564800dd8cSBarry Smith 3574800dd8cSBarry Smith /* Update size of trust region */ 3584800dd8cSBarry Smith if (rho < neP->mu) delta *= neP->delta1; 3594800dd8cSBarry Smith else if (rho < neP->eta) delta *= neP->delta2; 3604800dd8cSBarry Smith else delta *= neP->delta3; 3619566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm)); 3629566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta)); 3631aa26658SKarl Rupp 3644800dd8cSBarry Smith neP->delta = delta; 3654800dd8cSBarry Smith if (rho > neP->sigma) break; 3669566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Trying again in smaller region\n")); 367284fb49fSHeeho Park 3686b5873e3SBarry Smith /* check to see if progress is hopeless */ 369ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 3709566063dSJacob Faibussowitsch PetscCall(SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP)); 371*dbbe0bcdSBarry Smith if (!reason) PetscUseTypeMethod(snes,converged ,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); 3727cb011f5SBarry Smith if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER; 373184914b5SBarry Smith if (reason) { 37452392280SLois Curfman McInnes /* We're not progressing, so return with the current iterate */ 3759566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes,i+1,fnorm)); 376a7cc72afSBarry Smith breakout = PETSC_TRUE; 377454a90a3SBarry Smith break; 37852392280SLois Curfman McInnes } 37950ffb88aSMatthew Knepley snes->numFailures++; 3804800dd8cSBarry Smith } 3811acabf8cSLois Curfman McInnes if (!breakout) { 38285385478SLisandro Dalcin /* Update function and solution vectors */ 3834800dd8cSBarry Smith fnorm = gnorm; 3849566063dSJacob Faibussowitsch PetscCall(VecCopy(G,F)); 3859566063dSJacob Faibussowitsch PetscCall(VecCopy(W,X)); 38685385478SLisandro Dalcin /* Monitor convergence */ 3879566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 388c509a162SBarry Smith snes->iter = i+1; 389fbe28522SBarry Smith snes->norm = fnorm; 390c1e67a49SFande Kong snes->xnorm = xnorm; 391c1e67a49SFande Kong snes->ynorm = ynorm; 3929566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 3939566063dSJacob Faibussowitsch PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits)); 3949566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes,snes->iter,snes->norm)); 39585385478SLisandro Dalcin /* Test for convergence, xnorm = || X || */ 396ef87ac0dSKris Buschelman neP->itflag = PETSC_TRUE; 3979566063dSJacob Faibussowitsch if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm)); 398*dbbe0bcdSBarry Smith PetscUseTypeMethod(snes,converged ,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); 3993f149594SLisandro Dalcin if (reason) break; 4001aa26658SKarl Rupp } else break; 40138442cffSBarry Smith } 402284fb49fSHeeho Park 40352392280SLois Curfman McInnes if (i == maxits) { 4049566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits)); 40585385478SLisandro Dalcin if (!reason) reason = SNES_DIVERGED_MAX_IT; 40652392280SLois Curfman McInnes } 4079566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 408184914b5SBarry Smith snes->reason = reason; 4099566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 4105e28bcb6Sprj- if (convtest != SNESTR_KSPConverged_Private) { 4119566063dSJacob Faibussowitsch PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy)); 4129566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 4139566063dSJacob Faibussowitsch PetscCall(KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy)); 4145e28bcb6Sprj- } 4153a40ed3dSBarry Smith PetscFunctionReturn(0); 4164800dd8cSBarry Smith } 417284fb49fSHeeho Park 4184800dd8cSBarry Smith /*------------------------------------------------------------*/ 41904d7464bSBarry Smith static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) 4204800dd8cSBarry Smith { 4213a40ed3dSBarry Smith PetscFunctionBegin; 4229566063dSJacob Faibussowitsch PetscCall(SNESSetWorkVecs(snes,4)); 4239566063dSJacob Faibussowitsch PetscCall(SNESSetUpMatrices(snes)); 4243a40ed3dSBarry Smith PetscFunctionReturn(0); 4254800dd8cSBarry Smith } 4266b8b9a38SLisandro Dalcin 42704d7464bSBarry Smith PetscErrorCode SNESReset_NEWTONTR(SNES snes) 4286b8b9a38SLisandro Dalcin { 4296b8b9a38SLisandro Dalcin PetscFunctionBegin; 4306b8b9a38SLisandro Dalcin PetscFunctionReturn(0); 4316b8b9a38SLisandro Dalcin } 4326b8b9a38SLisandro Dalcin 43304d7464bSBarry Smith static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) 4344800dd8cSBarry Smith { 4353a40ed3dSBarry Smith PetscFunctionBegin; 4369566063dSJacob Faibussowitsch PetscCall(SNESReset_NEWTONTR(snes)); 4379566063dSJacob Faibussowitsch PetscCall(PetscFree(snes->data)); 4383a40ed3dSBarry Smith PetscFunctionReturn(0); 4394800dd8cSBarry Smith } 4404800dd8cSBarry Smith /*------------------------------------------------------------*/ 4414800dd8cSBarry Smith 442*dbbe0bcdSBarry Smith static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes,PetscOptionItems *PetscOptionsObject) 4434800dd8cSBarry Smith { 44404d7464bSBarry Smith SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data; 4454800dd8cSBarry Smith 4463a40ed3dSBarry Smith PetscFunctionBegin; 447d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject,"SNES trust region options for nonlinear equations"); 4489566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL)); 4499566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL)); 4509566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL)); 4519566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL)); 4529566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL)); 4539566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL)); 4549566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL)); 4559566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL)); 456d0609cedSBarry Smith PetscOptionsHeadEnd(); 4573a40ed3dSBarry Smith PetscFunctionReturn(0); 4584800dd8cSBarry Smith } 4594800dd8cSBarry Smith 46004d7464bSBarry Smith static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer) 461a935fc98SLois Curfman McInnes { 46204d7464bSBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 463ace3abfcSBarry Smith PetscBool iascii; 464a935fc98SLois Curfman McInnes 4653a40ed3dSBarry Smith PetscFunctionBegin; 4669566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 46732077d6dSBarry Smith if (iascii) { 46863a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," Trust region tolerance %g (-snes_trtol)\n",(double)snes->deltatol)); 4699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma)); 4709566063dSJacob 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)); 47119bcc07fSBarry Smith } 4723a40ed3dSBarry Smith PetscFunctionReturn(0); 473a935fc98SLois Curfman McInnes } 47452392280SLois Curfman McInnes /* ------------------------------------------------------------ */ 475ebe3b25bSBarry Smith /*MC 47604d7464bSBarry Smith SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 477ebe3b25bSBarry Smith 478ebe3b25bSBarry Smith Options Database: 4798e434772SBarry Smith + -snes_trtol <tol> - trust region tolerance 4808e434772SBarry Smith . -snes_tr_mu <mu> - trust region parameter 4818e434772SBarry Smith . -snes_tr_eta <eta> - trust region parameter 4828e434772SBarry Smith . -snes_tr_sigma <sigma> - trust region parameter 4838b84755bSBarry Smith . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 4848e434772SBarry Smith . -snes_tr_delta1 <delta1> - trust region parameter 4858e434772SBarry Smith . -snes_tr_delta2 <delta2> - trust region parameter 4868e434772SBarry Smith - -snes_tr_delta3 <delta3> - trust region parameter 487ebe3b25bSBarry Smith 488ebe3b25bSBarry Smith The basic algorithm is taken from "The Minpack Project", by More', 489ebe3b25bSBarry Smith Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 490ebe3b25bSBarry Smith of Mathematical Software", Wayne Cowell, editor. 491ebe3b25bSBarry Smith 492ee3001cbSBarry Smith Level: intermediate 493ee3001cbSBarry Smith 494db781477SPatrick Sanan .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()` 495ebe3b25bSBarry Smith 496ebe3b25bSBarry Smith M*/ 4978cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) 4984800dd8cSBarry Smith { 49904d7464bSBarry Smith SNES_NEWTONTR *neP; 5004800dd8cSBarry Smith 5013a40ed3dSBarry Smith PetscFunctionBegin; 50204d7464bSBarry Smith snes->ops->setup = SNESSetUp_NEWTONTR; 50304d7464bSBarry Smith snes->ops->solve = SNESSolve_NEWTONTR; 50404d7464bSBarry Smith snes->ops->destroy = SNESDestroy_NEWTONTR; 50504d7464bSBarry Smith snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 50604d7464bSBarry Smith snes->ops->view = SNESView_NEWTONTR; 50704d7464bSBarry Smith snes->ops->reset = SNESReset_NEWTONTR; 508fbe28522SBarry Smith 50942f4f86dSBarry Smith snes->usesksp = PETSC_TRUE; 510efd4aadfSBarry Smith snes->usesnpc = PETSC_FALSE; 51142f4f86dSBarry Smith 5124fc747eaSLawrence Mitchell snes->alwayscomputesfinalresidual = PETSC_TRUE; 5134fc747eaSLawrence Mitchell 5149566063dSJacob Faibussowitsch PetscCall(PetscNewLog(snes,&neP)); 515fbe28522SBarry Smith snes->data = (void*)neP; 516fbe28522SBarry Smith neP->mu = 0.25; 517fbe28522SBarry Smith neP->eta = 0.75; 518fbe28522SBarry Smith neP->delta = 0.0; 519fbe28522SBarry Smith neP->delta0 = 0.2; 520fbe28522SBarry Smith neP->delta1 = 0.3; 521fbe28522SBarry Smith neP->delta2 = 0.75; 522fbe28522SBarry Smith neP->delta3 = 2.0; 523fbe28522SBarry Smith neP->sigma = 0.0001; 524ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 52575567043SBarry Smith neP->rnorm0 = 0.0; 52675567043SBarry Smith neP->ttol = 0.0; 5273a40ed3dSBarry Smith PetscFunctionReturn(0); 5284800dd8cSBarry Smith } 529