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