xref: /petsc/src/snes/impls/tr/tr.c (revision dbbe0bcd3f3a8fbab5a45420dc06f8387e5764c6)
1c6db04a5SJed Brown #include <../src/snes/impls/tr/trimpl.h>                /*I   "petscsnes.h"   I*/
24800dd8cSBarry Smith 
3971273eeSBarry Smith typedef struct {
4971273eeSBarry Smith   SNES           snes;
5df8705c3SBarry Smith   /*  Information on the regular SNES convergence test; which may have been user provided */
6df8705c3SBarry Smith   PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
7df8705c3SBarry Smith   PetscErrorCode (*convdestroy)(void*);
8df8705c3SBarry Smith   void           *convctx;
9971273eeSBarry Smith } SNES_TR_KSPConverged_Ctx;
10971273eeSBarry Smith 
11470880abSPatrick Sanan static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
12df60cc22SBarry Smith {
13971273eeSBarry Smith   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
14971273eeSBarry Smith   SNES                     snes = ctx->snes;
1504d7464bSBarry Smith   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
16df60cc22SBarry Smith   Vec                      x;
17064f8208SBarry Smith   PetscReal                nrm;
18df60cc22SBarry Smith 
193a40ed3dSBarry Smith   PetscFunctionBegin;
209566063dSJacob Faibussowitsch   PetscCall((*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx));
21329f5518SBarry Smith   if (*reason) {
2263a3b9bcSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n",n,(double)rnorm));
239b04593eSLois Curfman McInnes   }
24a935fc98SLois Curfman McInnes   /* Determine norm of solution */
259566063dSJacob Faibussowitsch   PetscCall(KSPBuildSolution(ksp,NULL,&x));
269566063dSJacob Faibussowitsch   PetscCall(VecNorm(x,NORM_2,&nrm));
27064f8208SBarry Smith   if (nrm >= neP->delta) {
289566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm));
29329f5518SBarry Smith     *reason = KSP_CONVERGED_STEP_LENGTH;
30df60cc22SBarry Smith   }
313a40ed3dSBarry Smith   PetscFunctionReturn(0);
32df60cc22SBarry Smith }
3382bf6240SBarry Smith 
34470880abSPatrick Sanan static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
35971273eeSBarry Smith {
36971273eeSBarry Smith   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
37971273eeSBarry Smith 
38971273eeSBarry Smith   PetscFunctionBegin;
399566063dSJacob Faibussowitsch   PetscCall((*ctx->convdestroy)(ctx->convctx));
409566063dSJacob Faibussowitsch   PetscCall(PetscFree(ctx));
41971273eeSBarry Smith   PetscFunctionReturn(0);
42971273eeSBarry Smith }
43971273eeSBarry Smith 
4485385478SLisandro Dalcin /* ---------------------------------------------------------------- */
4585385478SLisandro Dalcin /*
46470880abSPatrick Sanan    SNESTR_Converged_Private -test convergence JUST for
4785385478SLisandro Dalcin    the trust region tolerance.
4885385478SLisandro Dalcin 
4985385478SLisandro Dalcin */
50470880abSPatrick Sanan static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
5185385478SLisandro Dalcin {
5204d7464bSBarry Smith   SNES_NEWTONTR  *neP = (SNES_NEWTONTR*)snes->data;
5385385478SLisandro Dalcin 
5485385478SLisandro Dalcin   PetscFunctionBegin;
5585385478SLisandro Dalcin   *reason = SNES_CONVERGED_ITERATING;
5685385478SLisandro Dalcin   if (neP->delta < xnorm * snes->deltatol) {
579566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol));
581c6b2ff8SBarry Smith     *reason = SNES_DIVERGED_TR_DELTA;
59e71169deSBarry Smith   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
6063a3b9bcSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n",snes->max_funcs));
6185385478SLisandro Dalcin     *reason = SNES_DIVERGED_FUNCTION_COUNT;
6285385478SLisandro Dalcin   }
6385385478SLisandro Dalcin   PetscFunctionReturn(0);
6485385478SLisandro Dalcin }
6585385478SLisandro Dalcin 
667cb011f5SBarry Smith /*@C
67c9368356SGlenn Hammond    SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
68c9368356SGlenn Hammond        Allows the user a chance to change or override the decision of the line search routine.
69c9368356SGlenn Hammond 
70c9368356SGlenn Hammond    Logically Collective on snes
71c9368356SGlenn Hammond 
72c9368356SGlenn Hammond    Input Parameters:
73c9368356SGlenn Hammond +  snes - the nonlinear solver object
74c9368356SGlenn Hammond .  func - [optional] function evaluation routine, see SNESNewtonTRPreCheck()  for the calling sequence
75c9368356SGlenn Hammond -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
76c9368356SGlenn Hammond 
77c9368356SGlenn Hammond    Level: intermediate
78c9368356SGlenn Hammond 
79c9368356SGlenn Hammond    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver.
80c9368356SGlenn Hammond 
81db781477SPatrick Sanan .seealso: `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
82c9368356SGlenn Hammond @*/
83c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
84c9368356SGlenn Hammond {
85c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
86c9368356SGlenn Hammond 
87c9368356SGlenn Hammond   PetscFunctionBegin;
88c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
89c9368356SGlenn Hammond   if (func) tr->precheck    = func;
90c9368356SGlenn Hammond   if (ctx)  tr->precheckctx = ctx;
91c9368356SGlenn Hammond   PetscFunctionReturn(0);
92c9368356SGlenn Hammond }
93c9368356SGlenn Hammond 
94c9368356SGlenn Hammond /*@C
95c9368356SGlenn Hammond    SNESNewtonTRGetPreCheck - Gets the pre-check function
96c9368356SGlenn Hammond 
97c9368356SGlenn Hammond    Not collective
98c9368356SGlenn Hammond 
99c9368356SGlenn Hammond    Input Parameter:
100c9368356SGlenn Hammond .  snes - the nonlinear solver context
101c9368356SGlenn Hammond 
102c9368356SGlenn Hammond    Output Parameters:
103c9368356SGlenn Hammond +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck()
104c9368356SGlenn Hammond -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
105c9368356SGlenn Hammond 
106c9368356SGlenn Hammond    Level: intermediate
107c9368356SGlenn Hammond 
108db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
109c9368356SGlenn Hammond @*/
110c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
111c9368356SGlenn Hammond {
112c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
113c9368356SGlenn Hammond 
114c9368356SGlenn Hammond   PetscFunctionBegin;
115c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
116c9368356SGlenn Hammond   if (func) *func = tr->precheck;
117c9368356SGlenn Hammond   if (ctx)  *ctx  = tr->precheckctx;
118c9368356SGlenn Hammond   PetscFunctionReturn(0);
119c9368356SGlenn Hammond }
120c9368356SGlenn Hammond 
121c9368356SGlenn Hammond /*@C
1227cb011f5SBarry Smith    SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
1237cb011f5SBarry Smith        function evaluation. Allows the user a chance to change or override the decision of the line search routine
1247cb011f5SBarry Smith 
1257cb011f5SBarry Smith    Logically Collective on snes
1267cb011f5SBarry Smith 
1277cb011f5SBarry Smith    Input Parameters:
1287cb011f5SBarry Smith +  snes - the nonlinear solver object
1297cb011f5SBarry Smith .  func - [optional] function evaluation routine, see SNESNewtonTRPostCheck()  for the calling sequence
1307cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1317cb011f5SBarry Smith 
1327cb011f5SBarry Smith    Level: intermediate
1337cb011f5SBarry Smith 
134c9368356SGlenn Hammond    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in
1357cb011f5SBarry Smith    SNESLineSearchSetPostCheck() is called AFTER the function evaluation.
1367cb011f5SBarry Smith 
137db781477SPatrick Sanan .seealso: `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`
1387cb011f5SBarry Smith @*/
139c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
1407cb011f5SBarry Smith {
1417cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
1427cb011f5SBarry Smith 
1437cb011f5SBarry Smith   PetscFunctionBegin;
1447cb011f5SBarry Smith   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
1457cb011f5SBarry Smith   if (func) tr->postcheck    = func;
1467cb011f5SBarry Smith   if (ctx)  tr->postcheckctx = ctx;
1477cb011f5SBarry Smith   PetscFunctionReturn(0);
1487cb011f5SBarry Smith }
1497cb011f5SBarry Smith 
1507cb011f5SBarry Smith /*@C
1517cb011f5SBarry Smith    SNESNewtonTRGetPostCheck - Gets the post-check function
1527cb011f5SBarry Smith 
1537cb011f5SBarry Smith    Not collective
1547cb011f5SBarry Smith 
1557cb011f5SBarry Smith    Input Parameter:
1567cb011f5SBarry Smith .  snes - the nonlinear solver context
1577cb011f5SBarry Smith 
1587cb011f5SBarry Smith    Output Parameters:
1597cb011f5SBarry Smith +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck()
1607cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1617cb011f5SBarry Smith 
1627cb011f5SBarry Smith    Level: intermediate
1637cb011f5SBarry Smith 
164db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
1657cb011f5SBarry Smith @*/
166c9368356SGlenn Hammond PetscErrorCode  SNESNewtonTRGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
1677cb011f5SBarry Smith {
1687cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
1697cb011f5SBarry Smith 
1707cb011f5SBarry Smith   PetscFunctionBegin;
1717cb011f5SBarry Smith   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
1727cb011f5SBarry Smith   if (func) *func = tr->postcheck;
1737cb011f5SBarry Smith   if (ctx)  *ctx  = tr->postcheckctx;
1747cb011f5SBarry Smith   PetscFunctionReturn(0);
1757cb011f5SBarry Smith }
1767cb011f5SBarry Smith 
1777cb011f5SBarry Smith /*@C
178c9368356SGlenn Hammond    SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR
179c9368356SGlenn Hammond 
180c9368356SGlenn Hammond    Logically Collective on snes
181c9368356SGlenn Hammond 
182c9368356SGlenn Hammond    Input Parameters:
183c9368356SGlenn Hammond +  snes - the solver
184c9368356SGlenn Hammond .  X - The last solution
185c9368356SGlenn Hammond -  Y - The step direction
186c9368356SGlenn Hammond 
187c9368356SGlenn Hammond    Output Parameters:
188c9368356SGlenn Hammond .  changed_Y - Indicator that the step direction Y has been changed.
189c9368356SGlenn Hammond 
190c9368356SGlenn Hammond    Level: developer
191c9368356SGlenn Hammond 
192db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
193c9368356SGlenn Hammond @*/
194c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
195c9368356SGlenn Hammond {
196c9368356SGlenn Hammond   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
197c9368356SGlenn Hammond 
198c9368356SGlenn Hammond   PetscFunctionBegin;
199c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
200c9368356SGlenn Hammond   if (tr->precheck) {
2019566063dSJacob Faibussowitsch     PetscCall((*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx));
202c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes,*changed_Y,4);
203c9368356SGlenn Hammond   }
204c9368356SGlenn Hammond   PetscFunctionReturn(0);
205c9368356SGlenn Hammond }
206c9368356SGlenn Hammond 
207c9368356SGlenn Hammond /*@C
2087cb011f5SBarry Smith    SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation
2097cb011f5SBarry Smith 
2107cb011f5SBarry Smith    Logically Collective on snes
2117cb011f5SBarry Smith 
2127cb011f5SBarry Smith    Input Parameters:
2136b867d5aSJose E. Roman +  snes - the solver
2146b867d5aSJose E. Roman .  X - The last solution
2157cb011f5SBarry Smith .  Y - The full step direction
2163312a946SBarry Smith -  W - The updated solution, W = X - Y
2177cb011f5SBarry Smith 
2187cb011f5SBarry Smith    Output Parameters:
2193312a946SBarry Smith +  changed_Y - indicator if step has been changed
2203312a946SBarry Smith -  changed_W - Indicator if the new candidate solution W has been changed.
2217cb011f5SBarry Smith 
2227cb011f5SBarry Smith    Notes:
2233312a946SBarry Smith      If Y is changed then W is recomputed as X - Y
2247cb011f5SBarry Smith 
2257cb011f5SBarry Smith    Level: developer
2267cb011f5SBarry Smith 
227db781477SPatrick Sanan .seealso: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
2287cb011f5SBarry Smith @*/
229c9368356SGlenn Hammond static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
2307cb011f5SBarry Smith {
2317cb011f5SBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
2327cb011f5SBarry Smith 
2337cb011f5SBarry Smith   PetscFunctionBegin;
234c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
2357cb011f5SBarry Smith   *changed_W = PETSC_FALSE;
2367cb011f5SBarry Smith   if (tr->postcheck) {
2379566063dSJacob Faibussowitsch     PetscCall((*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx));
238c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes,*changed_Y,5);
2397cb011f5SBarry Smith     PetscValidLogicalCollectiveBool(snes,*changed_W,6);
2407cb011f5SBarry Smith   }
2417cb011f5SBarry Smith   PetscFunctionReturn(0);
2427cb011f5SBarry Smith }
24385385478SLisandro Dalcin 
244df60cc22SBarry Smith /*
24504d7464bSBarry Smith    SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
246ddbbdb52SLois Curfman McInnes    region approach for solving systems of nonlinear equations.
2474800dd8cSBarry Smith 
2484800dd8cSBarry Smith */
24904d7464bSBarry Smith static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
2504800dd8cSBarry Smith {
25104d7464bSBarry Smith   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
252c9368356SGlenn Hammond   Vec                      X,F,Y,G,Ytmp,W;
253a7cc72afSBarry Smith   PetscInt                 maxits,i,lits;
25485385478SLisandro Dalcin   PetscReal                rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
25579f36460SBarry Smith   PetscScalar              cnorm;
256df60cc22SBarry Smith   KSP                      ksp;
2573f149594SLisandro Dalcin   SNESConvergedReason      reason = SNES_CONVERGED_ITERATING;
258df8705c3SBarry Smith   PetscBool                breakout = PETSC_FALSE;
259df8705c3SBarry Smith   SNES_TR_KSPConverged_Ctx *ctx;
2605e28bcb6Sprj-   PetscErrorCode           (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
2615e28bcb6Sprj-   void                     *convctx;
2624800dd8cSBarry Smith 
2633a40ed3dSBarry Smith   PetscFunctionBegin;
2640b121fc5SBarry Smith   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
265c579b300SPatrick Farrell 
266fbe28522SBarry Smith   maxits = snes->max_its;               /* maximum number of iterations */
267fbe28522SBarry Smith   X      = snes->vec_sol;               /* solution vector */
26839e2f89bSBarry Smith   F      = snes->vec_func;              /* residual vector */
269fbe28522SBarry Smith   Y      = snes->work[0];               /* work vectors */
270fbe28522SBarry Smith   G      = snes->work[1];
2716b5873e3SBarry Smith   Ytmp   = snes->work[2];
272c9368356SGlenn Hammond   W      = snes->work[3];
2734800dd8cSBarry Smith 
2749566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
2754c49b128SBarry Smith   snes->iter = 0;
2769566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2774800dd8cSBarry Smith 
278df8705c3SBarry Smith   /* Set the linear stopping criteria to use the More' trick. */
2799566063dSJacob Faibussowitsch   PetscCall(SNESGetKSP(snes,&ksp));
2809566063dSJacob Faibussowitsch   PetscCall(KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy));
281df8705c3SBarry Smith   if (convtest != SNESTR_KSPConverged_Private) {
2829566063dSJacob Faibussowitsch     PetscCall(PetscNew(&ctx));
283df8705c3SBarry Smith     ctx->snes             = snes;
2849566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy));
2859566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy));
2869566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
287df8705c3SBarry Smith   }
288df8705c3SBarry Smith 
289e4ed7901SPeter Brune   if (!snes->vec_func_init_set) {
2909566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(snes,X,F));          /* F(X) */
2911aa26658SKarl Rupp   } else snes->vec_func_init_set = PETSC_FALSE;
292e4ed7901SPeter Brune 
2939566063dSJacob Faibussowitsch   PetscCall(VecNorm(F,NORM_2,&fnorm));             /* fnorm <- || F || */
294422a814eSBarry Smith   SNESCheckFunctionNorm(snes,fnorm);
2959566063dSJacob Faibussowitsch   PetscCall(VecNorm(X,NORM_2,&xnorm));             /* xnorm <- || X || */
2969566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
297fbe28522SBarry Smith   snes->norm = fnorm;
2989566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2998b84755bSBarry Smith   delta      = xnorm ? neP->delta0*xnorm : neP->delta0;
3004800dd8cSBarry Smith   neP->delta = delta;
3019566063dSJacob Faibussowitsch   PetscCall(SNESLogConvergenceHistory(snes,fnorm,0));
3029566063dSJacob Faibussowitsch   PetscCall(SNESMonitor(snes,0,fnorm));
303b37302c6SLois Curfman McInnes 
30485385478SLisandro Dalcin   /* test convergence */
305*dbbe0bcdSBarry Smith   PetscUseTypeMethod(snes,converged ,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
30685385478SLisandro Dalcin   if (snes->reason) PetscFunctionReturn(0);
3073f149594SLisandro Dalcin 
3084800dd8cSBarry Smith   for (i=0; i<maxits; i++) {
30999a96b7cSMatthew Knepley 
31099a96b7cSMatthew Knepley     /* Call general purpose update function */
311*dbbe0bcdSBarry Smith     PetscTryTypeMethod(snes,update, snes->iter);
31299a96b7cSMatthew Knepley 
31385385478SLisandro Dalcin     /* Solve J Y = F, where J is Jacobian matrix */
3149566063dSJacob Faibussowitsch     PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre));
31507b62357SFande Kong     SNESCheckJacobianDomainerror(snes);
3169566063dSJacob Faibussowitsch     PetscCall(KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre));
3179566063dSJacob Faibussowitsch     PetscCall(KSPSolve(snes->ksp,F,Ytmp));
3189566063dSJacob Faibussowitsch     PetscCall(KSPGetIterationNumber(snes->ksp,&lits));
319329f5518SBarry Smith     snes->linear_its += lits;
3201aa26658SKarl Rupp 
3219566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits));
3229566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ytmp,NORM_2,&nrm));
323064f8208SBarry Smith     norm1 = nrm;
324284fb49fSHeeho Park 
3256b5873e3SBarry Smith     while (1) {
326c9368356SGlenn Hammond       PetscBool changed_y;
3277cb011f5SBarry Smith       PetscBool changed_w;
3289566063dSJacob Faibussowitsch       PetscCall(VecCopy(Ytmp,Y));
329064f8208SBarry Smith       nrm  = norm1;
3306b5873e3SBarry Smith 
3312deea200SLois Curfman McInnes       /* Scale Y if need be and predict new value of F norm */
332064f8208SBarry Smith       if (nrm >= delta) {
333064f8208SBarry Smith         nrm    = delta/nrm;
334064f8208SBarry Smith         gpnorm = (1.0 - nrm)*fnorm;
335064f8208SBarry Smith         cnorm  = nrm;
3369566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes,"Scaling direction by %g\n",(double)nrm));
3379566063dSJacob Faibussowitsch         PetscCall(VecScale(Y,cnorm));
338064f8208SBarry Smith         nrm    = gpnorm;
339fbe28522SBarry Smith         ynorm  = delta;
340fbe28522SBarry Smith       } else {
341fbe28522SBarry Smith         gpnorm = 0.0;
3429566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes,"Direction is in Trust Region\n"));
343064f8208SBarry Smith         ynorm  = nrm;
344fbe28522SBarry Smith       }
345c9368356SGlenn Hammond       /* PreCheck() allows for updates to Y prior to W <- X - Y */
3469566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPreCheck(snes,X,Y,&changed_y));
3479566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W,-1.0,Y,X));         /* W <- X - Y */
3489566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w));
3499566063dSJacob Faibussowitsch       if (changed_y) PetscCall(VecWAXPY(W,-1.0,Y,X));
3509566063dSJacob Faibussowitsch       PetscCall(VecCopy(Y,snes->vec_sol_update));
3519566063dSJacob Faibussowitsch       PetscCall(SNESComputeFunction(snes,W,G)); /*  F(X-Y) = G */
3529566063dSJacob Faibussowitsch       PetscCall(VecNorm(G,NORM_2,&gnorm));      /* gnorm <- || g || */
3537551ef16SBarry Smith       SNESCheckFunctionNorm(snes,gnorm);
3544800dd8cSBarry Smith       if (fnorm == gpnorm) rho = 0.0;
355fbe28522SBarry Smith       else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
3564800dd8cSBarry Smith 
3574800dd8cSBarry Smith       /* Update size of trust region */
3584800dd8cSBarry Smith       if      (rho < neP->mu)  delta *= neP->delta1;
3594800dd8cSBarry Smith       else if (rho < neP->eta) delta *= neP->delta2;
3604800dd8cSBarry Smith       else                     delta *= neP->delta3;
3619566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm));
3629566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta));
3631aa26658SKarl Rupp 
3644800dd8cSBarry Smith       neP->delta = delta;
3654800dd8cSBarry Smith       if (rho > neP->sigma) break;
3669566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Trying again in smaller region\n"));
367284fb49fSHeeho Park 
3686b5873e3SBarry Smith       /* check to see if progress is hopeless */
369ef87ac0dSKris Buschelman       neP->itflag = PETSC_FALSE;
3709566063dSJacob Faibussowitsch       PetscCall(SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP));
371*dbbe0bcdSBarry Smith       if (!reason) PetscUseTypeMethod(snes,converged ,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
3727cb011f5SBarry Smith       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
373184914b5SBarry Smith       if (reason) {
37452392280SLois Curfman McInnes         /* We're not progressing, so return with the current iterate */
3759566063dSJacob Faibussowitsch         PetscCall(SNESMonitor(snes,i+1,fnorm));
376a7cc72afSBarry Smith         breakout = PETSC_TRUE;
377454a90a3SBarry Smith         break;
37852392280SLois Curfman McInnes       }
37950ffb88aSMatthew Knepley       snes->numFailures++;
3804800dd8cSBarry Smith     }
3811acabf8cSLois Curfman McInnes     if (!breakout) {
38285385478SLisandro Dalcin       /* Update function and solution vectors */
3834800dd8cSBarry Smith       fnorm = gnorm;
3849566063dSJacob Faibussowitsch       PetscCall(VecCopy(G,F));
3859566063dSJacob Faibussowitsch       PetscCall(VecCopy(W,X));
38685385478SLisandro Dalcin       /* Monitor convergence */
3879566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
388c509a162SBarry Smith       snes->iter = i+1;
389fbe28522SBarry Smith       snes->norm = fnorm;
390c1e67a49SFande Kong       snes->xnorm = xnorm;
391c1e67a49SFande Kong       snes->ynorm = ynorm;
3929566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
3939566063dSJacob Faibussowitsch       PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits));
3949566063dSJacob Faibussowitsch       PetscCall(SNESMonitor(snes,snes->iter,snes->norm));
39585385478SLisandro Dalcin       /* Test for convergence, xnorm = || X || */
396ef87ac0dSKris Buschelman       neP->itflag = PETSC_TRUE;
3979566063dSJacob Faibussowitsch       if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm));
398*dbbe0bcdSBarry Smith       PetscUseTypeMethod(snes,converged ,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
3993f149594SLisandro Dalcin       if (reason) break;
4001aa26658SKarl Rupp     } else break;
40138442cffSBarry Smith   }
402284fb49fSHeeho Park 
40352392280SLois Curfman McInnes   if (i == maxits) {
4049566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits));
40585385478SLisandro Dalcin     if (!reason) reason = SNES_DIVERGED_MAX_IT;
40652392280SLois Curfman McInnes   }
4079566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
408184914b5SBarry Smith   snes->reason = reason;
4099566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
4105e28bcb6Sprj-   if (convtest != SNESTR_KSPConverged_Private) {
4119566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy));
4129566063dSJacob Faibussowitsch     PetscCall(PetscFree(ctx));
4139566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy));
4145e28bcb6Sprj-   }
4153a40ed3dSBarry Smith   PetscFunctionReturn(0);
4164800dd8cSBarry Smith }
417284fb49fSHeeho Park 
4184800dd8cSBarry Smith /*------------------------------------------------------------*/
41904d7464bSBarry Smith static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
4204800dd8cSBarry Smith {
4213a40ed3dSBarry Smith   PetscFunctionBegin;
4229566063dSJacob Faibussowitsch   PetscCall(SNESSetWorkVecs(snes,4));
4239566063dSJacob Faibussowitsch   PetscCall(SNESSetUpMatrices(snes));
4243a40ed3dSBarry Smith   PetscFunctionReturn(0);
4254800dd8cSBarry Smith }
4266b8b9a38SLisandro Dalcin 
42704d7464bSBarry Smith PetscErrorCode SNESReset_NEWTONTR(SNES snes)
4286b8b9a38SLisandro Dalcin {
4296b8b9a38SLisandro Dalcin   PetscFunctionBegin;
4306b8b9a38SLisandro Dalcin   PetscFunctionReturn(0);
4316b8b9a38SLisandro Dalcin }
4326b8b9a38SLisandro Dalcin 
43304d7464bSBarry Smith static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
4344800dd8cSBarry Smith {
4353a40ed3dSBarry Smith   PetscFunctionBegin;
4369566063dSJacob Faibussowitsch   PetscCall(SNESReset_NEWTONTR(snes));
4379566063dSJacob Faibussowitsch   PetscCall(PetscFree(snes->data));
4383a40ed3dSBarry Smith   PetscFunctionReturn(0);
4394800dd8cSBarry Smith }
4404800dd8cSBarry Smith /*------------------------------------------------------------*/
4414800dd8cSBarry Smith 
442*dbbe0bcdSBarry Smith static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes,PetscOptionItems *PetscOptionsObject)
4434800dd8cSBarry Smith {
44404d7464bSBarry Smith   SNES_NEWTONTR  *ctx = (SNES_NEWTONTR*)snes->data;
4454800dd8cSBarry Smith 
4463a40ed3dSBarry Smith   PetscFunctionBegin;
447d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject,"SNES trust region options for nonlinear equations");
4489566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL));
4499566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL));
4509566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL));
4519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL));
4529566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL));
4539566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL));
4549566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL));
4559566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL));
456d0609cedSBarry Smith   PetscOptionsHeadEnd();
4573a40ed3dSBarry Smith   PetscFunctionReturn(0);
4584800dd8cSBarry Smith }
4594800dd8cSBarry Smith 
46004d7464bSBarry Smith static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
461a935fc98SLois Curfman McInnes {
46204d7464bSBarry Smith   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
463ace3abfcSBarry Smith   PetscBool      iascii;
464a935fc98SLois Curfman McInnes 
4653a40ed3dSBarry Smith   PetscFunctionBegin;
4669566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
46732077d6dSBarry Smith   if (iascii) {
46863a3b9bcSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer,"  Trust region tolerance %g (-snes_trtol)\n",(double)snes->deltatol));
4699566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer,"  mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma));
4709566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer,"  delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3));
47119bcc07fSBarry Smith   }
4723a40ed3dSBarry Smith   PetscFunctionReturn(0);
473a935fc98SLois Curfman McInnes }
47452392280SLois Curfman McInnes /* ------------------------------------------------------------ */
475ebe3b25bSBarry Smith /*MC
47604d7464bSBarry Smith       SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
477ebe3b25bSBarry Smith 
478ebe3b25bSBarry Smith    Options Database:
4798e434772SBarry Smith +    -snes_trtol <tol> - trust region tolerance
4808e434772SBarry Smith .    -snes_tr_mu <mu> - trust region parameter
4818e434772SBarry Smith .    -snes_tr_eta <eta> - trust region parameter
4828e434772SBarry Smith .    -snes_tr_sigma <sigma> - trust region parameter
4838b84755bSBarry Smith .    -snes_tr_delta0 <delta0> -  initial size of the trust region is delta0*norm2(x)
4848e434772SBarry Smith .    -snes_tr_delta1 <delta1> - trust region parameter
4858e434772SBarry Smith .    -snes_tr_delta2 <delta2> - trust region parameter
4868e434772SBarry Smith -    -snes_tr_delta3 <delta3> - trust region parameter
487ebe3b25bSBarry Smith 
488ebe3b25bSBarry Smith    The basic algorithm is taken from "The Minpack Project", by More',
489ebe3b25bSBarry Smith    Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
490ebe3b25bSBarry Smith    of Mathematical Software", Wayne Cowell, editor.
491ebe3b25bSBarry Smith 
492ee3001cbSBarry Smith    Level: intermediate
493ee3001cbSBarry Smith 
494db781477SPatrick Sanan .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`
495ebe3b25bSBarry Smith 
496ebe3b25bSBarry Smith M*/
4978cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
4984800dd8cSBarry Smith {
49904d7464bSBarry Smith   SNES_NEWTONTR  *neP;
5004800dd8cSBarry Smith 
5013a40ed3dSBarry Smith   PetscFunctionBegin;
50204d7464bSBarry Smith   snes->ops->setup          = SNESSetUp_NEWTONTR;
50304d7464bSBarry Smith   snes->ops->solve          = SNESSolve_NEWTONTR;
50404d7464bSBarry Smith   snes->ops->destroy        = SNESDestroy_NEWTONTR;
50504d7464bSBarry Smith   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
50604d7464bSBarry Smith   snes->ops->view           = SNESView_NEWTONTR;
50704d7464bSBarry Smith   snes->ops->reset          = SNESReset_NEWTONTR;
508fbe28522SBarry Smith 
50942f4f86dSBarry Smith   snes->usesksp = PETSC_TRUE;
510efd4aadfSBarry Smith   snes->usesnpc = PETSC_FALSE;
51142f4f86dSBarry Smith 
5124fc747eaSLawrence Mitchell   snes->alwayscomputesfinalresidual = PETSC_TRUE;
5134fc747eaSLawrence Mitchell 
5149566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(snes,&neP));
515fbe28522SBarry Smith   snes->data  = (void*)neP;
516fbe28522SBarry Smith   neP->mu     = 0.25;
517fbe28522SBarry Smith   neP->eta    = 0.75;
518fbe28522SBarry Smith   neP->delta  = 0.0;
519fbe28522SBarry Smith   neP->delta0 = 0.2;
520fbe28522SBarry Smith   neP->delta1 = 0.3;
521fbe28522SBarry Smith   neP->delta2 = 0.75;
522fbe28522SBarry Smith   neP->delta3 = 2.0;
523fbe28522SBarry Smith   neP->sigma  = 0.0001;
524ef87ac0dSKris Buschelman   neP->itflag = PETSC_FALSE;
52575567043SBarry Smith   neP->rnorm0 = 0.0;
52675567043SBarry Smith   neP->ttol   = 0.0;
5273a40ed3dSBarry Smith   PetscFunctionReturn(0);
5284800dd8cSBarry Smith }
529