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 119371c9d4SSatish Balay static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx) { 12971273eeSBarry Smith SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx; 13971273eeSBarry Smith SNES snes = ctx->snes; 1404d7464bSBarry Smith SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data; 15df60cc22SBarry Smith Vec x; 16064f8208SBarry Smith PetscReal nrm; 17df60cc22SBarry Smith 183a40ed3dSBarry Smith PetscFunctionBegin; 199566063dSJacob Faibussowitsch PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx)); 2048a46eb9SPierre Jolivet if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm)); 21a935fc98SLois Curfman McInnes /* Determine norm of solution */ 229566063dSJacob Faibussowitsch PetscCall(KSPBuildSolution(ksp, NULL, &x)); 239566063dSJacob Faibussowitsch PetscCall(VecNorm(x, NORM_2, &nrm)); 24064f8208SBarry Smith if (nrm >= neP->delta) { 259566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm)); 26329f5518SBarry Smith *reason = KSP_CONVERGED_STEP_LENGTH; 27df60cc22SBarry Smith } 283a40ed3dSBarry Smith PetscFunctionReturn(0); 29df60cc22SBarry Smith } 3082bf6240SBarry Smith 319371c9d4SSatish Balay static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx) { 32971273eeSBarry Smith SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx; 33971273eeSBarry Smith 34971273eeSBarry Smith PetscFunctionBegin; 359566063dSJacob Faibussowitsch PetscCall((*ctx->convdestroy)(ctx->convctx)); 369566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 37971273eeSBarry Smith PetscFunctionReturn(0); 38971273eeSBarry Smith } 39971273eeSBarry Smith 4085385478SLisandro Dalcin /* 41470880abSPatrick Sanan SNESTR_Converged_Private -test convergence JUST for 4285385478SLisandro Dalcin the trust region tolerance. 4385385478SLisandro Dalcin */ 449371c9d4SSatish Balay static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy) { 4504d7464bSBarry Smith SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data; 4685385478SLisandro Dalcin 4785385478SLisandro Dalcin PetscFunctionBegin; 4885385478SLisandro Dalcin *reason = SNES_CONVERGED_ITERATING; 4985385478SLisandro Dalcin if (neP->delta < xnorm * snes->deltatol) { 509566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Converged due to trust region param %g<%g*%g\n", (double)neP->delta, (double)xnorm, (double)snes->deltatol)); 511c6b2ff8SBarry Smith *reason = SNES_DIVERGED_TR_DELTA; 52e71169deSBarry Smith } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 5363a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs)); 5485385478SLisandro Dalcin *reason = SNES_DIVERGED_FUNCTION_COUNT; 5585385478SLisandro Dalcin } 5685385478SLisandro Dalcin PetscFunctionReturn(0); 5785385478SLisandro Dalcin } 5885385478SLisandro Dalcin 597cb011f5SBarry Smith /*@C 60c9368356SGlenn Hammond SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined. 61c9368356SGlenn Hammond Allows the user a chance to change or override the decision of the line search routine. 62c9368356SGlenn Hammond 63f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 64f6dfbefdSBarry Smith 65c9368356SGlenn Hammond Logically Collective on snes 66c9368356SGlenn Hammond 67c9368356SGlenn Hammond Input Parameters: 68c9368356SGlenn Hammond + snes - the nonlinear solver object 69f6dfbefdSBarry Smith . func - [optional] function evaluation routine, see `SNESNewtonTRPreCheck()` for the calling sequence 70c9368356SGlenn Hammond - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 71c9368356SGlenn Hammond 72c9368356SGlenn Hammond Level: intermediate 73c9368356SGlenn Hammond 74f6dfbefdSBarry Smith Note: 75f6dfbefdSBarry Smith This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver. 76c9368356SGlenn Hammond 77f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()` 78c9368356SGlenn Hammond @*/ 799371c9d4SSatish Balay PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx) { 80c9368356SGlenn Hammond SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 81c9368356SGlenn Hammond 82c9368356SGlenn Hammond PetscFunctionBegin; 83c9368356SGlenn Hammond PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 84c9368356SGlenn Hammond if (func) tr->precheck = func; 85c9368356SGlenn Hammond if (ctx) tr->precheckctx = ctx; 86c9368356SGlenn Hammond PetscFunctionReturn(0); 87c9368356SGlenn Hammond } 88c9368356SGlenn Hammond 89c9368356SGlenn Hammond /*@C 90c9368356SGlenn Hammond SNESNewtonTRGetPreCheck - Gets the pre-check function 91c9368356SGlenn Hammond 92f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 93f6dfbefdSBarry Smith 94c9368356SGlenn Hammond Not collective 95c9368356SGlenn Hammond 96c9368356SGlenn Hammond Input Parameter: 97c9368356SGlenn Hammond . snes - the nonlinear solver context 98c9368356SGlenn Hammond 99c9368356SGlenn Hammond Output Parameters: 100f6dfbefdSBarry Smith + func - [optional] function evaluation routine, see for the calling sequence `SNESNewtonTRPreCheck()` 101c9368356SGlenn Hammond - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 102c9368356SGlenn Hammond 103c9368356SGlenn Hammond Level: intermediate 104c9368356SGlenn Hammond 105f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()` 106c9368356SGlenn Hammond @*/ 1079371c9d4SSatish Balay PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx) { 108c9368356SGlenn Hammond SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 109c9368356SGlenn Hammond 110c9368356SGlenn Hammond PetscFunctionBegin; 111c9368356SGlenn Hammond PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 112c9368356SGlenn Hammond if (func) *func = tr->precheck; 113c9368356SGlenn Hammond if (ctx) *ctx = tr->precheckctx; 114c9368356SGlenn Hammond PetscFunctionReturn(0); 115c9368356SGlenn Hammond } 116c9368356SGlenn Hammond 117c9368356SGlenn Hammond /*@C 1187cb011f5SBarry Smith SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next 1197cb011f5SBarry Smith function evaluation. Allows the user a chance to change or override the decision of the line search routine 1207cb011f5SBarry Smith 121f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 122f6dfbefdSBarry Smith 1237cb011f5SBarry Smith Logically Collective on snes 1247cb011f5SBarry Smith 1257cb011f5SBarry Smith Input Parameters: 1267cb011f5SBarry Smith + snes - the nonlinear solver object 127f6dfbefdSBarry Smith . func - [optional] function evaluation routine, see `SNESNewtonTRPostCheck()` for the calling sequence 1287cb011f5SBarry Smith - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 1297cb011f5SBarry Smith 1307cb011f5SBarry Smith Level: intermediate 1317cb011f5SBarry Smith 132f6dfbefdSBarry Smith Note: 133f6dfbefdSBarry Smith This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver while the function set in 134f6dfbefdSBarry Smith `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation. 1357cb011f5SBarry Smith 136f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()` 1377cb011f5SBarry Smith @*/ 1389371c9d4SSatish Balay PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx) { 1397cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 1407cb011f5SBarry Smith 1417cb011f5SBarry Smith PetscFunctionBegin; 1427cb011f5SBarry Smith PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1437cb011f5SBarry Smith if (func) tr->postcheck = func; 1447cb011f5SBarry Smith if (ctx) tr->postcheckctx = ctx; 1457cb011f5SBarry Smith PetscFunctionReturn(0); 1467cb011f5SBarry Smith } 1477cb011f5SBarry Smith 1487cb011f5SBarry Smith /*@C 1497cb011f5SBarry Smith SNESNewtonTRGetPostCheck - Gets the post-check function 1507cb011f5SBarry Smith 151f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 152f6dfbefdSBarry Smith 1537cb011f5SBarry Smith Not collective 1547cb011f5SBarry Smith 1557cb011f5SBarry Smith Input Parameter: 1567cb011f5SBarry Smith . snes - the nonlinear solver context 1577cb011f5SBarry Smith 1587cb011f5SBarry Smith Output Parameters: 159f6dfbefdSBarry 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 164f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()` 1657cb011f5SBarry Smith @*/ 1669371c9d4SSatish Balay PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx) { 1677cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 1687cb011f5SBarry Smith 1697cb011f5SBarry Smith PetscFunctionBegin; 1707cb011f5SBarry Smith PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 1717cb011f5SBarry Smith if (func) *func = tr->postcheck; 1727cb011f5SBarry Smith if (ctx) *ctx = tr->postcheckctx; 1737cb011f5SBarry Smith PetscFunctionReturn(0); 1747cb011f5SBarry Smith } 1757cb011f5SBarry Smith 1767cb011f5SBarry Smith /*@C 177f6dfbefdSBarry Smith SNESNewtonTRPreCheck - Called before the step has been determined in `SNESNEWTONTR` 178f6dfbefdSBarry Smith 179f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 180c9368356SGlenn Hammond 181c9368356SGlenn Hammond Logically Collective on snes 182c9368356SGlenn Hammond 183c9368356SGlenn Hammond Input Parameters: 184c9368356SGlenn Hammond + snes - the solver 185c9368356SGlenn Hammond . X - The last solution 186c9368356SGlenn Hammond - Y - The step direction 187c9368356SGlenn Hammond 188c9368356SGlenn Hammond Output Parameters: 189c9368356SGlenn Hammond . changed_Y - Indicator that the step direction Y has been changed. 190c9368356SGlenn Hammond 191c9368356SGlenn Hammond Level: developer 192c9368356SGlenn Hammond 193f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()` 194c9368356SGlenn Hammond @*/ 1959371c9d4SSatish Balay static PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y) { 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 208f6dfbefdSBarry Smith SNESNewtonTRPostCheck - Called after the step has been determined in `SNESNEWTONTR` but before the function evaluation 209f6dfbefdSBarry Smith 210f6dfbefdSBarry Smith Deprecated use `SNESNEWTONDCTRDC` 2117cb011f5SBarry Smith 2127cb011f5SBarry Smith Logically Collective on snes 2137cb011f5SBarry Smith 2147cb011f5SBarry Smith Input Parameters: 2156b867d5aSJose E. Roman + snes - the solver 2166b867d5aSJose E. Roman . X - The last solution 2177cb011f5SBarry Smith . Y - The full step direction 2183312a946SBarry Smith - W - The updated solution, W = X - Y 2197cb011f5SBarry Smith 2207cb011f5SBarry Smith Output Parameters: 2213312a946SBarry Smith + changed_Y - indicator if step has been changed 2223312a946SBarry Smith - changed_W - Indicator if the new candidate solution W has been changed. 2237cb011f5SBarry Smith 224f6dfbefdSBarry Smith Note: 2253312a946SBarry Smith If Y is changed then W is recomputed as X - Y 2267cb011f5SBarry Smith 2277cb011f5SBarry Smith Level: developer 2287cb011f5SBarry Smith 229f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()` 2307cb011f5SBarry Smith @*/ 2319371c9d4SSatish Balay static PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W) { 2327cb011f5SBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 2337cb011f5SBarry Smith 2347cb011f5SBarry Smith PetscFunctionBegin; 235c9368356SGlenn Hammond *changed_Y = PETSC_FALSE; 2367cb011f5SBarry Smith *changed_W = PETSC_FALSE; 2377cb011f5SBarry Smith if (tr->postcheck) { 2389566063dSJacob Faibussowitsch PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx)); 239c9368356SGlenn Hammond PetscValidLogicalCollectiveBool(snes, *changed_Y, 5); 2407cb011f5SBarry Smith PetscValidLogicalCollectiveBool(snes, *changed_W, 6); 2417cb011f5SBarry Smith } 2427cb011f5SBarry Smith PetscFunctionReturn(0); 2437cb011f5SBarry Smith } 24485385478SLisandro Dalcin 245df60cc22SBarry Smith /* 24604d7464bSBarry Smith SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust 247ddbbdb52SLois Curfman McInnes region approach for solving systems of nonlinear equations. 2484800dd8cSBarry Smith 2494800dd8cSBarry Smith */ 2509371c9d4SSatish Balay static PetscErrorCode SNESSolve_NEWTONTR(SNES snes) { 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 */ 305dbbe0bcdSBarry 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 /* Call general purpose update function */ 310dbbe0bcdSBarry Smith PetscTryTypeMethod(snes, update, snes->iter); 31199a96b7cSMatthew Knepley 31285385478SLisandro Dalcin /* Solve J Y = F, where J is Jacobian matrix */ 3139566063dSJacob Faibussowitsch PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre)); 31407b62357SFande Kong SNESCheckJacobianDomainerror(snes); 3159566063dSJacob Faibussowitsch PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre)); 3169566063dSJacob Faibussowitsch PetscCall(KSPSolve(snes->ksp, F, Ytmp)); 3179566063dSJacob Faibussowitsch PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); 318329f5518SBarry Smith snes->linear_its += lits; 3191aa26658SKarl Rupp 3209566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits)); 3219566063dSJacob Faibussowitsch PetscCall(VecNorm(Ytmp, NORM_2, &nrm)); 322064f8208SBarry Smith norm1 = nrm; 323284fb49fSHeeho Park 3246b5873e3SBarry Smith while (1) { 325c9368356SGlenn Hammond PetscBool changed_y; 3267cb011f5SBarry Smith PetscBool changed_w; 3279566063dSJacob Faibussowitsch PetscCall(VecCopy(Ytmp, Y)); 328064f8208SBarry Smith nrm = norm1; 3296b5873e3SBarry Smith 3302deea200SLois Curfman McInnes /* Scale Y if need be and predict new value of F norm */ 331064f8208SBarry Smith if (nrm >= delta) { 332064f8208SBarry Smith nrm = delta / nrm; 333064f8208SBarry Smith gpnorm = (1.0 - nrm) * fnorm; 334064f8208SBarry Smith cnorm = nrm; 3359566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Scaling direction by %g\n", (double)nrm)); 3369566063dSJacob Faibussowitsch PetscCall(VecScale(Y, cnorm)); 337064f8208SBarry Smith nrm = gpnorm; 338fbe28522SBarry Smith ynorm = delta; 339fbe28522SBarry Smith } else { 340fbe28522SBarry Smith gpnorm = 0.0; 3419566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Direction is in Trust Region\n")); 342064f8208SBarry Smith ynorm = nrm; 343fbe28522SBarry Smith } 344c9368356SGlenn Hammond /* PreCheck() allows for updates to Y prior to W <- X - Y */ 3459566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y)); 3469566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -1.0, Y, X)); /* W <- X - Y */ 3479566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w)); 3489566063dSJacob Faibussowitsch if (changed_y) PetscCall(VecWAXPY(W, -1.0, Y, X)); 3499566063dSJacob Faibussowitsch PetscCall(VecCopy(Y, snes->vec_sol_update)); 3509566063dSJacob Faibussowitsch PetscCall(SNESComputeFunction(snes, W, G)); /* F(X-Y) = G */ 3519566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); /* gnorm <- || g || */ 3527551ef16SBarry Smith SNESCheckFunctionNorm(snes, gnorm); 3534800dd8cSBarry Smith if (fnorm == gpnorm) rho = 0.0; 354fbe28522SBarry Smith else rho = (fnorm * fnorm - gnorm * gnorm) / (fnorm * fnorm - gpnorm * gpnorm); 3554800dd8cSBarry Smith 3564800dd8cSBarry Smith /* Update size of trust region */ 3574800dd8cSBarry Smith if (rho < neP->mu) delta *= neP->delta1; 3584800dd8cSBarry Smith else if (rho < neP->eta) delta *= neP->delta2; 3594800dd8cSBarry Smith else delta *= neP->delta3; 3609566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "fnorm=%g, gnorm=%g, ynorm=%g\n", (double)fnorm, (double)gnorm, (double)ynorm)); 3619566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "gpred=%g, rho=%g, delta=%g\n", (double)gpnorm, (double)rho, (double)delta)); 3621aa26658SKarl Rupp 3634800dd8cSBarry Smith neP->delta = delta; 3644800dd8cSBarry Smith if (rho > neP->sigma) break; 3659566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Trying again in smaller region\n")); 366284fb49fSHeeho Park 3676b5873e3SBarry Smith /* check to see if progress is hopeless */ 368ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 3699566063dSJacob Faibussowitsch PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP)); 370dbbe0bcdSBarry Smith if (!reason) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP); 3717cb011f5SBarry Smith if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER; 372184914b5SBarry Smith if (reason) { 37352392280SLois Curfman McInnes /* We're not progressing, so return with the current iterate */ 3749566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes, i + 1, fnorm)); 375a7cc72afSBarry Smith breakout = PETSC_TRUE; 376454a90a3SBarry Smith break; 37752392280SLois Curfman McInnes } 37850ffb88aSMatthew Knepley snes->numFailures++; 3794800dd8cSBarry Smith } 3801acabf8cSLois Curfman McInnes if (!breakout) { 38185385478SLisandro Dalcin /* Update function and solution vectors */ 3824800dd8cSBarry Smith fnorm = gnorm; 3839566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3849566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 38585385478SLisandro Dalcin /* Monitor convergence */ 3869566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 387c509a162SBarry Smith snes->iter = i + 1; 388fbe28522SBarry Smith snes->norm = fnorm; 389c1e67a49SFande Kong snes->xnorm = xnorm; 390c1e67a49SFande Kong snes->ynorm = ynorm; 3919566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 3929566063dSJacob Faibussowitsch PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 3939566063dSJacob Faibussowitsch PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 39485385478SLisandro Dalcin /* Test for convergence, xnorm = || X || */ 395ef87ac0dSKris Buschelman neP->itflag = PETSC_TRUE; 3969566063dSJacob Faibussowitsch if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm)); 397dbbe0bcdSBarry Smith PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP); 3983f149594SLisandro Dalcin if (reason) break; 3991aa26658SKarl Rupp } else break; 40038442cffSBarry Smith } 401284fb49fSHeeho Park 40252392280SLois Curfman McInnes if (i == maxits) { 4039566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 40485385478SLisandro Dalcin if (!reason) reason = SNES_DIVERGED_MAX_IT; 40552392280SLois Curfman McInnes } 4069566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 407184914b5SBarry Smith snes->reason = reason; 4089566063dSJacob Faibussowitsch PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 4095e28bcb6Sprj- if (convtest != SNESTR_KSPConverged_Private) { 4109566063dSJacob Faibussowitsch PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy)); 4119566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 4129566063dSJacob Faibussowitsch PetscCall(KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy)); 4135e28bcb6Sprj- } 4143a40ed3dSBarry Smith PetscFunctionReturn(0); 4154800dd8cSBarry Smith } 416284fb49fSHeeho Park 4179371c9d4SSatish Balay static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) { 4183a40ed3dSBarry Smith PetscFunctionBegin; 4199566063dSJacob Faibussowitsch PetscCall(SNESSetWorkVecs(snes, 4)); 4209566063dSJacob Faibussowitsch PetscCall(SNESSetUpMatrices(snes)); 4213a40ed3dSBarry Smith PetscFunctionReturn(0); 4224800dd8cSBarry Smith } 4236b8b9a38SLisandro Dalcin 4249371c9d4SSatish Balay PetscErrorCode SNESReset_NEWTONTR(SNES snes) { 4256b8b9a38SLisandro Dalcin PetscFunctionBegin; 4266b8b9a38SLisandro Dalcin PetscFunctionReturn(0); 4276b8b9a38SLisandro Dalcin } 4286b8b9a38SLisandro Dalcin 4299371c9d4SSatish Balay static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) { 4303a40ed3dSBarry Smith PetscFunctionBegin; 4319566063dSJacob Faibussowitsch PetscCall(SNESReset_NEWTONTR(snes)); 4329566063dSJacob Faibussowitsch PetscCall(PetscFree(snes->data)); 4333a40ed3dSBarry Smith PetscFunctionReturn(0); 4344800dd8cSBarry Smith } 4354800dd8cSBarry Smith 4369371c9d4SSatish Balay static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject) { 43704d7464bSBarry Smith SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data; 4384800dd8cSBarry Smith 4393a40ed3dSBarry Smith PetscFunctionBegin; 440d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations"); 4419566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_trtol", "Trust region tolerance", "SNESSetTrustRegionTolerance", snes->deltatol, &snes->deltatol, NULL)); 4429566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_mu", "mu", "None", ctx->mu, &ctx->mu, NULL)); 4439566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_eta", "eta", "None", ctx->eta, &ctx->eta, NULL)); 4449566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_sigma", "sigma", "None", ctx->sigma, &ctx->sigma, NULL)); 4459566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL)); 4469566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta1", "delta1", "None", ctx->delta1, &ctx->delta1, NULL)); 4479566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta2", "delta2", "None", ctx->delta2, &ctx->delta2, NULL)); 4489566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_tr_delta3", "delta3", "None", ctx->delta3, &ctx->delta3, NULL)); 449d0609cedSBarry Smith PetscOptionsHeadEnd(); 4503a40ed3dSBarry Smith PetscFunctionReturn(0); 4514800dd8cSBarry Smith } 4524800dd8cSBarry Smith 4539371c9d4SSatish Balay static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer) { 45404d7464bSBarry Smith SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data; 455ace3abfcSBarry Smith PetscBool iascii; 456a935fc98SLois Curfman McInnes 4573a40ed3dSBarry Smith PetscFunctionBegin; 4589566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 45932077d6dSBarry Smith if (iascii) { 46063a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g (-snes_trtol)\n", (double)snes->deltatol)); 4619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " mu=%g, eta=%g, sigma=%g\n", (double)tr->mu, (double)tr->eta, (double)tr->sigma)); 4629566063dSJacob 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)); 46319bcc07fSBarry Smith } 4643a40ed3dSBarry Smith PetscFunctionReturn(0); 465a935fc98SLois Curfman McInnes } 466f6dfbefdSBarry Smith 467ebe3b25bSBarry Smith /*MC 46804d7464bSBarry Smith SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 469ebe3b25bSBarry Smith 470f6dfbefdSBarry Smith Deprecated use `SNESNEWTONTRDC` 471f6dfbefdSBarry Smith 472f6dfbefdSBarry Smith Options Database Keys: 4738e434772SBarry Smith + -snes_trtol <tol> - trust region tolerance 4748e434772SBarry Smith . -snes_tr_mu <mu> - trust region parameter 4758e434772SBarry Smith . -snes_tr_eta <eta> - trust region parameter 4768e434772SBarry Smith . -snes_tr_sigma <sigma> - trust region parameter 4778b84755bSBarry Smith . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 4788e434772SBarry Smith . -snes_tr_delta1 <delta1> - trust region parameter 4798e434772SBarry Smith . -snes_tr_delta2 <delta2> - trust region parameter 4808e434772SBarry Smith - -snes_tr_delta3 <delta3> - trust region parameter 481ebe3b25bSBarry Smith 482f6dfbefdSBarry Smith Reference: 483f6dfbefdSBarry Smith . - * "The Minpack Project", by More', Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 484ebe3b25bSBarry Smith of Mathematical Software", Wayne Cowell, editor. 485ebe3b25bSBarry Smith 486ee3001cbSBarry Smith Level: intermediate 487ee3001cbSBarry Smith 488f6dfbefdSBarry Smith .seealso: `SNESNEWTONTRDC`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()` 489ebe3b25bSBarry Smith M*/ 4909371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) { 49104d7464bSBarry Smith SNES_NEWTONTR *neP; 4924800dd8cSBarry Smith 4933a40ed3dSBarry Smith PetscFunctionBegin; 49404d7464bSBarry Smith snes->ops->setup = SNESSetUp_NEWTONTR; 49504d7464bSBarry Smith snes->ops->solve = SNESSolve_NEWTONTR; 49604d7464bSBarry Smith snes->ops->destroy = SNESDestroy_NEWTONTR; 49704d7464bSBarry Smith snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 49804d7464bSBarry Smith snes->ops->view = SNESView_NEWTONTR; 49904d7464bSBarry Smith snes->ops->reset = SNESReset_NEWTONTR; 500fbe28522SBarry Smith 50142f4f86dSBarry Smith snes->usesksp = PETSC_TRUE; 502efd4aadfSBarry Smith snes->usesnpc = PETSC_FALSE; 50342f4f86dSBarry Smith 5044fc747eaSLawrence Mitchell snes->alwayscomputesfinalresidual = PETSC_TRUE; 5054fc747eaSLawrence Mitchell 506*4dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&neP)); 507fbe28522SBarry Smith snes->data = (void *)neP; 508fbe28522SBarry Smith neP->mu = 0.25; 509fbe28522SBarry Smith neP->eta = 0.75; 510fbe28522SBarry Smith neP->delta = 0.0; 511fbe28522SBarry Smith neP->delta0 = 0.2; 512fbe28522SBarry Smith neP->delta1 = 0.3; 513fbe28522SBarry Smith neP->delta2 = 0.75; 514fbe28522SBarry Smith neP->delta3 = 2.0; 515fbe28522SBarry Smith neP->sigma = 0.0001; 516ef87ac0dSKris Buschelman neP->itflag = PETSC_FALSE; 51775567043SBarry Smith neP->rnorm0 = 0.0; 51875567043SBarry Smith neP->ttol = 0.0; 5193a40ed3dSBarry Smith PetscFunctionReturn(0); 5204800dd8cSBarry Smith } 521