xref: /petsc/src/snes/impls/tr/tr.c (revision 63a3b9bc7a1f24f247904ccba9383635fe6abade)
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) {
22*63a3b9bcSJacob 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) {
60*63a3b9bcSJacob 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 
813312a946SBarry Smith .seealso: SNESNewtonTRPreCheck(), SNESNewtonTRGetPreCheck(), SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
82c9368356SGlenn Hammond @*/
83c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
84c9368356SGlenn Hammond {
85c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
86c9368356SGlenn Hammond 
87c9368356SGlenn Hammond   PetscFunctionBegin;
88c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
89c9368356SGlenn Hammond   if (func) tr->precheck    = func;
90c9368356SGlenn Hammond   if (ctx)  tr->precheckctx = ctx;
91c9368356SGlenn Hammond   PetscFunctionReturn(0);
92c9368356SGlenn Hammond }
93c9368356SGlenn Hammond 
94c9368356SGlenn Hammond /*@C
95c9368356SGlenn Hammond    SNESNewtonTRGetPreCheck - Gets the pre-check function
96c9368356SGlenn Hammond 
97c9368356SGlenn Hammond    Not collective
98c9368356SGlenn Hammond 
99c9368356SGlenn Hammond    Input Parameter:
100c9368356SGlenn Hammond .  snes - the nonlinear solver context
101c9368356SGlenn Hammond 
102c9368356SGlenn Hammond    Output Parameters:
103c9368356SGlenn Hammond +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck()
104c9368356SGlenn Hammond -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
105c9368356SGlenn Hammond 
106c9368356SGlenn Hammond    Level: intermediate
107c9368356SGlenn Hammond 
108c9368356SGlenn Hammond .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRPreCheck()
109c9368356SGlenn Hammond @*/
110c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
111c9368356SGlenn Hammond {
112c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
113c9368356SGlenn Hammond 
114c9368356SGlenn Hammond   PetscFunctionBegin;
115c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
116c9368356SGlenn Hammond   if (func) *func = tr->precheck;
117c9368356SGlenn Hammond   if (ctx)  *ctx  = tr->precheckctx;
118c9368356SGlenn Hammond   PetscFunctionReturn(0);
119c9368356SGlenn Hammond }
120c9368356SGlenn Hammond 
121c9368356SGlenn Hammond /*@C
1227cb011f5SBarry Smith    SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
1237cb011f5SBarry Smith        function evaluation. Allows the user a chance to change or override the decision of the line search routine
1247cb011f5SBarry Smith 
1257cb011f5SBarry Smith    Logically Collective on snes
1267cb011f5SBarry Smith 
1277cb011f5SBarry Smith    Input Parameters:
1287cb011f5SBarry Smith +  snes - the nonlinear solver object
1297cb011f5SBarry Smith .  func - [optional] function evaluation routine, see SNESNewtonTRPostCheck()  for the calling sequence
1307cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1317cb011f5SBarry Smith 
1327cb011f5SBarry Smith    Level: intermediate
1337cb011f5SBarry Smith 
134c9368356SGlenn Hammond    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in
1357cb011f5SBarry Smith    SNESLineSearchSetPostCheck() is called AFTER the function evaluation.
1367cb011f5SBarry Smith 
1377cb011f5SBarry Smith .seealso: SNESNewtonTRPostCheck(), SNESNewtonTRGetPostCheck()
1387cb011f5SBarry Smith @*/
139c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
1407cb011f5SBarry Smith {
1417cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
1427cb011f5SBarry Smith 
1437cb011f5SBarry Smith   PetscFunctionBegin;
1447cb011f5SBarry Smith   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
1457cb011f5SBarry Smith   if (func) tr->postcheck    = func;
1467cb011f5SBarry Smith   if (ctx)  tr->postcheckctx = ctx;
1477cb011f5SBarry Smith   PetscFunctionReturn(0);
1487cb011f5SBarry Smith }
1497cb011f5SBarry Smith 
1507cb011f5SBarry Smith /*@C
1517cb011f5SBarry Smith    SNESNewtonTRGetPostCheck - Gets the post-check function
1527cb011f5SBarry Smith 
1537cb011f5SBarry Smith    Not collective
1547cb011f5SBarry Smith 
1557cb011f5SBarry Smith    Input Parameter:
1567cb011f5SBarry Smith .  snes - the nonlinear solver context
1577cb011f5SBarry Smith 
1587cb011f5SBarry Smith    Output Parameters:
1597cb011f5SBarry Smith +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck()
1607cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1617cb011f5SBarry Smith 
1627cb011f5SBarry Smith    Level: intermediate
1637cb011f5SBarry Smith 
1647cb011f5SBarry Smith .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRPostCheck()
1657cb011f5SBarry Smith @*/
166c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
1677cb011f5SBarry Smith {
1687cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
1697cb011f5SBarry Smith 
1707cb011f5SBarry Smith   PetscFunctionBegin;
1717cb011f5SBarry Smith   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
1727cb011f5SBarry Smith   if (func) *func = tr->postcheck;
1737cb011f5SBarry Smith   if (ctx)  *ctx  = tr->postcheckctx;
1747cb011f5SBarry Smith   PetscFunctionReturn(0);
1757cb011f5SBarry Smith }
1767cb011f5SBarry Smith 
1777cb011f5SBarry Smith /*@C
178c9368356SGlenn Hammond    SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR
179c9368356SGlenn Hammond 
180c9368356SGlenn Hammond    Logically Collective on snes
181c9368356SGlenn Hammond 
182c9368356SGlenn Hammond    Input Parameters:
183c9368356SGlenn Hammond +  snes - the solver
184c9368356SGlenn Hammond .  X - The last solution
185c9368356SGlenn Hammond -  Y - The step direction
186c9368356SGlenn Hammond 
187c9368356SGlenn Hammond    Output Parameters:
188c9368356SGlenn Hammond .  changed_Y - Indicator that the step direction Y has been changed.
189c9368356SGlenn Hammond 
190c9368356SGlenn Hammond    Level: developer
191c9368356SGlenn Hammond 
192c9368356SGlenn Hammond .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRGetPreCheck()
193c9368356SGlenn Hammond @*/
194c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
195c9368356SGlenn Hammond {
196c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
197c9368356SGlenn Hammond 
198c9368356SGlenn Hammond   PetscFunctionBegin;
199c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
200c9368356SGlenn Hammond   if (tr->precheck) {
2019566063dSJacob Faibussowitsch     PetscCall((*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx));
202c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes,*changed_Y,4);
203c9368356SGlenn Hammond   }
204c9368356SGlenn Hammond   PetscFunctionReturn(0);
205c9368356SGlenn Hammond }
206c9368356SGlenn Hammond 
207c9368356SGlenn Hammond /*@C
2087cb011f5SBarry Smith    SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation
2097cb011f5SBarry Smith 
2107cb011f5SBarry Smith    Logically Collective on snes
2117cb011f5SBarry Smith 
2127cb011f5SBarry Smith    Input Parameters:
2136b867d5aSJose E. Roman +  snes - the solver
2146b867d5aSJose E. Roman .  X - The last solution
2157cb011f5SBarry Smith .  Y - The full step direction
2163312a946SBarry Smith -  W - The updated solution, W = X - Y
2177cb011f5SBarry Smith 
2187cb011f5SBarry Smith    Output Parameters:
2193312a946SBarry Smith +  changed_Y - indicator if step has been changed
2203312a946SBarry Smith -  changed_W - Indicator if the new candidate solution W has been changed.
2217cb011f5SBarry Smith 
2227cb011f5SBarry Smith    Notes:
2233312a946SBarry Smith      If Y is changed then W is recomputed as X - Y
2247cb011f5SBarry Smith 
2257cb011f5SBarry Smith    Level: developer
2267cb011f5SBarry Smith 
2277cb011f5SBarry Smith .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
2287cb011f5SBarry Smith @*/
229c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
2307cb011f5SBarry Smith {
2317cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
2327cb011f5SBarry Smith 
2337cb011f5SBarry Smith   PetscFunctionBegin;
234c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
2357cb011f5SBarry Smith   *changed_W = PETSC_FALSE;
2367cb011f5SBarry Smith   if (tr->postcheck) {
2379566063dSJacob Faibussowitsch     PetscCall((*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx));
238c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes,*changed_Y,5);
2397cb011f5SBarry Smith     PetscValidLogicalCollectiveBool(snes,*changed_W,6);
2407cb011f5SBarry Smith   }
2417cb011f5SBarry Smith   PetscFunctionReturn(0);
2427cb011f5SBarry Smith }
24385385478SLisandro Dalcin 
244df60cc22SBarry Smith /*
24504d7464bSBarry Smith    SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
246ddbbdb52SLois Curfman McInnes    region approach for solving systems of nonlinear equations.
2474800dd8cSBarry Smith 
2484800dd8cSBarry Smith */
24904d7464bSBarry Smith static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
2504800dd8cSBarry Smith {
25104d7464bSBarry Smith   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
252c9368356SGlenn Hammond   Vec                      X,F,Y,G,Ytmp,W;
253a7cc72afSBarry Smith   PetscInt                 maxits,i,lits;
25485385478SLisandro Dalcin   PetscReal                rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
25579f36460SBarry Smith   PetscScalar              cnorm;
256df60cc22SBarry Smith   KSP                      ksp;
2573f149594SLisandro Dalcin   SNESConvergedReason      reason = SNES_CONVERGED_ITERATING;
258df8705c3SBarry Smith   PetscBool                breakout = PETSC_FALSE;
259df8705c3SBarry Smith   SNES_TR_KSPConverged_Ctx *ctx;
2605e28bcb6Sprj-   PetscErrorCode           (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
2615e28bcb6Sprj-   void                     *convctx;
2624800dd8cSBarry Smith 
2633a40ed3dSBarry Smith   PetscFunctionBegin;
2642c71b3e2SJacob Faibussowitsch   PetscCheckFalse(snes->xl || snes->xu || snes->ops->computevariablebounds,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
265c579b300SPatrick Farrell 
266fbe28522SBarry Smith   maxits = snes->max_its;               /* maximum number of iterations */
267fbe28522SBarry Smith   X      = snes->vec_sol;               /* solution vector */
26839e2f89bSBarry Smith   F      = snes->vec_func;              /* residual vector */
269fbe28522SBarry Smith   Y      = snes->work[0];               /* work vectors */
270fbe28522SBarry Smith   G      = snes->work[1];
2716b5873e3SBarry Smith   Ytmp   = snes->work[2];
272c9368356SGlenn Hammond   W      = snes->work[3];
2734800dd8cSBarry Smith 
2749566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
2754c49b128SBarry Smith   snes->iter = 0;
2769566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2774800dd8cSBarry Smith 
278df8705c3SBarry Smith   /* Set the linear stopping criteria to use the More' trick. */
2799566063dSJacob Faibussowitsch   PetscCall(SNESGetKSP(snes,&ksp));
2809566063dSJacob Faibussowitsch   PetscCall(KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy));
281df8705c3SBarry Smith   if (convtest != SNESTR_KSPConverged_Private) {
2829566063dSJacob Faibussowitsch     PetscCall(PetscNew(&ctx));
283df8705c3SBarry Smith     ctx->snes             = snes;
2849566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy));
2859566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy));
2869566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
287df8705c3SBarry Smith   }
288df8705c3SBarry Smith 
289e4ed7901SPeter Brune   if (!snes->vec_func_init_set) {
2909566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(snes,X,F));          /* F(X) */
2911aa26658SKarl Rupp   } else snes->vec_func_init_set = PETSC_FALSE;
292e4ed7901SPeter Brune 
2939566063dSJacob Faibussowitsch   PetscCall(VecNorm(F,NORM_2,&fnorm));             /* fnorm <- || F || */
294422a814eSBarry Smith   SNESCheckFunctionNorm(snes,fnorm);
2959566063dSJacob Faibussowitsch   PetscCall(VecNorm(X,NORM_2,&xnorm));             /* xnorm <- || X || */
2969566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
297fbe28522SBarry Smith   snes->norm = fnorm;
2989566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2998b84755bSBarry Smith   delta      = xnorm ? neP->delta0*xnorm : neP->delta0;
3004800dd8cSBarry Smith   neP->delta = delta;
3019566063dSJacob Faibussowitsch   PetscCall(SNESLogConvergenceHistory(snes,fnorm,0));
3029566063dSJacob Faibussowitsch   PetscCall(SNESMonitor(snes,0,fnorm));
303b37302c6SLois Curfman McInnes 
30485385478SLisandro Dalcin   /* test convergence */
3059566063dSJacob Faibussowitsch   PetscCall((*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP));
30685385478SLisandro Dalcin   if (snes->reason) PetscFunctionReturn(0);
3073f149594SLisandro Dalcin 
3084800dd8cSBarry Smith   for (i=0; i<maxits; i++) {
30999a96b7cSMatthew Knepley 
31099a96b7cSMatthew Knepley     /* Call general purpose update function */
311e7788613SBarry Smith     if (snes->ops->update) {
3129566063dSJacob Faibussowitsch       PetscCall((*snes->ops->update)(snes, snes->iter));
31399a96b7cSMatthew Knepley     }
31499a96b7cSMatthew Knepley 
31585385478SLisandro Dalcin     /* Solve J Y = F, where J is Jacobian matrix */
3169566063dSJacob Faibussowitsch     PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre));
31707b62357SFande Kong     SNESCheckJacobianDomainerror(snes);
3189566063dSJacob Faibussowitsch     PetscCall(KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre));
3199566063dSJacob Faibussowitsch     PetscCall(KSPSolve(snes->ksp,F,Ytmp));
3209566063dSJacob Faibussowitsch     PetscCall(KSPGetIterationNumber(snes->ksp,&lits));
321329f5518SBarry Smith     snes->linear_its += lits;
3221aa26658SKarl Rupp 
3239566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits));
3249566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ytmp,NORM_2,&nrm));
325064f8208SBarry Smith     norm1 = nrm;
326284fb49fSHeeho Park 
3276b5873e3SBarry Smith     while (1) {
328c9368356SGlenn Hammond       PetscBool changed_y;
3297cb011f5SBarry Smith       PetscBool changed_w;
3309566063dSJacob Faibussowitsch       PetscCall(VecCopy(Ytmp,Y));
331064f8208SBarry Smith       nrm  = norm1;
3326b5873e3SBarry Smith 
3332deea200SLois Curfman McInnes       /* Scale Y if need be and predict new value of F norm */
334064f8208SBarry Smith       if (nrm >= delta) {
335064f8208SBarry Smith         nrm    = delta/nrm;
336064f8208SBarry Smith         gpnorm = (1.0 - nrm)*fnorm;
337064f8208SBarry Smith         cnorm  = nrm;
3389566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes,"Scaling direction by %g\n",(double)nrm));
3399566063dSJacob Faibussowitsch         PetscCall(VecScale(Y,cnorm));
340064f8208SBarry Smith         nrm    = gpnorm;
341fbe28522SBarry Smith         ynorm  = delta;
342fbe28522SBarry Smith       } else {
343fbe28522SBarry Smith         gpnorm = 0.0;
3449566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes,"Direction is in Trust Region\n"));
345064f8208SBarry Smith         ynorm  = nrm;
346fbe28522SBarry Smith       }
347c9368356SGlenn Hammond       /* PreCheck() allows for updates to Y prior to W <- X - Y */
3489566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPreCheck(snes,X,Y,&changed_y));
3499566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W,-1.0,Y,X));         /* W <- X - Y */
3509566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w));
3519566063dSJacob Faibussowitsch       if (changed_y) PetscCall(VecWAXPY(W,-1.0,Y,X));
3529566063dSJacob Faibussowitsch       PetscCall(VecCopy(Y,snes->vec_sol_update));
3539566063dSJacob Faibussowitsch       PetscCall(SNESComputeFunction(snes,W,G)); /*  F(X-Y) = G */
3549566063dSJacob Faibussowitsch       PetscCall(VecNorm(G,NORM_2,&gnorm));      /* gnorm <- || g || */
3557551ef16SBarry Smith       SNESCheckFunctionNorm(snes,gnorm);
3564800dd8cSBarry Smith       if (fnorm == gpnorm) rho = 0.0;
357fbe28522SBarry Smith       else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
3584800dd8cSBarry Smith 
3594800dd8cSBarry Smith       /* Update size of trust region */
3604800dd8cSBarry Smith       if      (rho < neP->mu)  delta *= neP->delta1;
3614800dd8cSBarry Smith       else if (rho < neP->eta) delta *= neP->delta2;
3624800dd8cSBarry Smith       else                     delta *= neP->delta3;
3639566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm));
3649566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta));
3651aa26658SKarl Rupp 
3664800dd8cSBarry Smith       neP->delta = delta;
3674800dd8cSBarry Smith       if (rho > neP->sigma) break;
3689566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Trying again in smaller region\n"));
369284fb49fSHeeho Park 
3706b5873e3SBarry Smith       /* check to see if progress is hopeless */
371ef87ac0dSKris Buschelman       neP->itflag = PETSC_FALSE;
3729566063dSJacob Faibussowitsch       PetscCall(SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP));
3739566063dSJacob Faibussowitsch       if (!reason) PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP));
3747cb011f5SBarry Smith       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
375184914b5SBarry Smith       if (reason) {
37652392280SLois Curfman McInnes         /* We're not progressing, so return with the current iterate */
3779566063dSJacob Faibussowitsch         PetscCall(SNESMonitor(snes,i+1,fnorm));
378a7cc72afSBarry Smith         breakout = PETSC_TRUE;
379454a90a3SBarry Smith         break;
38052392280SLois Curfman McInnes       }
38150ffb88aSMatthew Knepley       snes->numFailures++;
3824800dd8cSBarry Smith     }
3831acabf8cSLois Curfman McInnes     if (!breakout) {
38485385478SLisandro Dalcin       /* Update function and solution vectors */
3854800dd8cSBarry Smith       fnorm = gnorm;
3869566063dSJacob Faibussowitsch       PetscCall(VecCopy(G,F));
3879566063dSJacob Faibussowitsch       PetscCall(VecCopy(W,X));
38885385478SLisandro Dalcin       /* Monitor convergence */
3899566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
390c509a162SBarry Smith       snes->iter = i+1;
391fbe28522SBarry Smith       snes->norm = fnorm;
392c1e67a49SFande Kong       snes->xnorm = xnorm;
393c1e67a49SFande Kong       snes->ynorm = ynorm;
3949566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
3959566063dSJacob Faibussowitsch       PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits));
3969566063dSJacob Faibussowitsch       PetscCall(SNESMonitor(snes,snes->iter,snes->norm));
39785385478SLisandro Dalcin       /* Test for convergence, xnorm = || X || */
398ef87ac0dSKris Buschelman       neP->itflag = PETSC_TRUE;
3999566063dSJacob Faibussowitsch       if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm));
4009566063dSJacob Faibussowitsch       PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP));
4013f149594SLisandro Dalcin       if (reason) break;
4021aa26658SKarl Rupp     } else break;
40338442cffSBarry Smith   }
404284fb49fSHeeho Park 
40552392280SLois Curfman McInnes   if (i == maxits) {
4069566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits));
40785385478SLisandro Dalcin     if (!reason) reason = SNES_DIVERGED_MAX_IT;
40852392280SLois Curfman McInnes   }
4099566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
410184914b5SBarry Smith   snes->reason = reason;
4119566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
4125e28bcb6Sprj-   if (convtest != SNESTR_KSPConverged_Private) {
4139566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy));
4149566063dSJacob Faibussowitsch     PetscCall(PetscFree(ctx));
4159566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy));
4165e28bcb6Sprj-   }
4173a40ed3dSBarry Smith   PetscFunctionReturn(0);
4184800dd8cSBarry Smith }
419284fb49fSHeeho Park 
4204800dd8cSBarry Smith /*------------------------------------------------------------*/
42104d7464bSBarry Smith static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
4224800dd8cSBarry Smith {
4233a40ed3dSBarry Smith   PetscFunctionBegin;
4249566063dSJacob Faibussowitsch   PetscCall(SNESSetWorkVecs(snes,4));
4259566063dSJacob Faibussowitsch   PetscCall(SNESSetUpMatrices(snes));
4263a40ed3dSBarry Smith   PetscFunctionReturn(0);
4274800dd8cSBarry Smith }
4286b8b9a38SLisandro Dalcin 
42904d7464bSBarry Smith PetscErrorCode SNESReset_NEWTONTR(SNES snes)
4306b8b9a38SLisandro Dalcin {
4316b8b9a38SLisandro Dalcin   PetscFunctionBegin;
4326b8b9a38SLisandro Dalcin   PetscFunctionReturn(0);
4336b8b9a38SLisandro Dalcin }
4346b8b9a38SLisandro Dalcin 
43504d7464bSBarry Smith static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
4364800dd8cSBarry Smith {
4373a40ed3dSBarry Smith   PetscFunctionBegin;
4389566063dSJacob Faibussowitsch   PetscCall(SNESReset_NEWTONTR(snes));
4399566063dSJacob Faibussowitsch   PetscCall(PetscFree(snes->data));
4403a40ed3dSBarry Smith   PetscFunctionReturn(0);
4414800dd8cSBarry Smith }
4424800dd8cSBarry Smith /*------------------------------------------------------------*/
4434800dd8cSBarry Smith 
4444416b707SBarry Smith static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
4454800dd8cSBarry Smith {
44604d7464bSBarry Smith   SNES_NEWTONTR  *ctx = (SNES_NEWTONTR*)snes->data;
4474800dd8cSBarry Smith 
4483a40ed3dSBarry Smith   PetscFunctionBegin;
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) {
470*63a3b9bcSJacob 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 
49604d7464bSBarry Smith .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
497ebe3b25bSBarry Smith 
498ebe3b25bSBarry Smith M*/
4998cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
5004800dd8cSBarry Smith {
50104d7464bSBarry Smith   SNES_NEWTONTR  *neP;
5024800dd8cSBarry Smith 
5033a40ed3dSBarry Smith   PetscFunctionBegin;
50404d7464bSBarry Smith   snes->ops->setup          = SNESSetUp_NEWTONTR;
50504d7464bSBarry Smith   snes->ops->solve          = SNESSolve_NEWTONTR;
50604d7464bSBarry Smith   snes->ops->destroy        = SNESDestroy_NEWTONTR;
50704d7464bSBarry Smith   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
50804d7464bSBarry Smith   snes->ops->view           = SNESView_NEWTONTR;
50904d7464bSBarry Smith   snes->ops->reset          = SNESReset_NEWTONTR;
510fbe28522SBarry Smith 
51142f4f86dSBarry Smith   snes->usesksp = PETSC_TRUE;
512efd4aadfSBarry Smith   snes->usesnpc = PETSC_FALSE;
51342f4f86dSBarry Smith 
5144fc747eaSLawrence Mitchell   snes->alwayscomputesfinalresidual = PETSC_TRUE;
5154fc747eaSLawrence Mitchell 
5169566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(snes,&neP));
517fbe28522SBarry Smith   snes->data  = (void*)neP;
518fbe28522SBarry Smith   neP->mu     = 0.25;
519fbe28522SBarry Smith   neP->eta    = 0.75;
520fbe28522SBarry Smith   neP->delta  = 0.0;
521fbe28522SBarry Smith   neP->delta0 = 0.2;
522fbe28522SBarry Smith   neP->delta1 = 0.3;
523fbe28522SBarry Smith   neP->delta2 = 0.75;
524fbe28522SBarry Smith   neP->delta3 = 2.0;
525fbe28522SBarry Smith   neP->sigma  = 0.0001;
526ef87ac0dSKris Buschelman   neP->itflag = PETSC_FALSE;
52775567043SBarry Smith   neP->rnorm0 = 0.0;
52875567043SBarry Smith   neP->ttol   = 0.0;
5293a40ed3dSBarry Smith   PetscFunctionReturn(0);
5304800dd8cSBarry Smith }
531