Actual source code: ntrdc.c


  2: #include <../src/snes/impls/ntrdc/ntrdcimpl.h>

  4: typedef struct {
  5:   SNES           snes;
  6:   /*  Information on the regular SNES convergence test; which may have been user provided
  7:       Copied from tr.c (maybe able to disposed, but this is a private function) - Heeho
  8:       Same with SNESTR_KSPConverged_Private, SNESTR_KSPConverged_Destroy, and SNESTR_Converged_Private
  9:  */

 11:   PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 12:   PetscErrorCode (*convdestroy)(void*);
 13:   void           *convctx;
 14: } SNES_TRDC_KSPConverged_Ctx;

 16: static PetscErrorCode SNESTRDC_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
 17: {
 18:   SNES_TRDC_KSPConverged_Ctx  *ctx = (SNES_TRDC_KSPConverged_Ctx*)cctx;
 19:   SNES                        snes = ctx->snes;
 20:   SNES_NEWTONTRDC             *neP = (SNES_NEWTONTRDC*)snes->data;
 21:   Vec                         x;
 22:   PetscReal                   nrm;

 24:   (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);
 25:   if (*reason) {
 26:     PetscInfo(snes,"Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n",n,(double)rnorm);
 27:   }
 28:   /* Determine norm of solution */
 29:   KSPBuildSolution(ksp,NULL,&x);
 30:   VecNorm(x,NORM_2,&nrm);
 31:   if (nrm >= neP->delta) {
 32:     PetscInfo(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
 33:     *reason = KSP_CONVERGED_STEP_LENGTH;
 34:   }
 35:   return 0;
 36: }

 38: static PetscErrorCode SNESTRDC_KSPConverged_Destroy(void *cctx)
 39: {
 40:   SNES_TRDC_KSPConverged_Ctx *ctx = (SNES_TRDC_KSPConverged_Ctx*)cctx;

 42:   (*ctx->convdestroy)(ctx->convctx);
 43:   PetscFree(ctx);

 45:   return 0;
 46: }

 48: /* ---------------------------------------------------------------- */
 49: /*
 50:    SNESTRDC_Converged_Private -test convergence JUST for
 51:    the trust region tolerance.

 53: */
 54: static PetscErrorCode SNESTRDC_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
 55: {
 56:   SNES_NEWTONTRDC  *neP = (SNES_NEWTONTRDC*)snes->data;

 58:   *reason = SNES_CONVERGED_ITERATING;
 59:   if (neP->delta < xnorm * snes->deltatol) {
 60:     PetscInfo(snes,"Diverged due to too small a trust region %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
 61:     *reason = SNES_DIVERGED_TR_DELTA;
 62:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
 63:     PetscInfo(snes,"Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n",snes->max_funcs);
 64:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
 65:   }
 66:   return 0;
 67: }

 69: /*@
 70:   SNESNewtonTRDCGetRhoFlag - Get whether the solution update is within the trust-region.

 72:   Input Parameters:
 73: . snes - the nonlinear solver object

 75:   Output Parameters:
 76: . rho_flag: PETSC_TRUE if the solution update is in the trust-region; otherwise, PETSC_FALSE

 78:   Level: developer

 80: @*/
 81: PetscErrorCode  SNESNewtonTRDCGetRhoFlag(SNES snes,PetscBool *rho_flag)
 82: {
 83:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

 87:   *rho_flag = tr->rho_satisfied;
 88:   return 0;
 89: }

 91: /*@C
 92:    SNESNewtonTRDCSetPreCheck - Sets a user function that is called before the search step has been determined.
 93:        Allows the user a chance to change or override the trust region decision.

 95:    Logically Collective on snes

 97:    Input Parameters:
 98: +  snes - the nonlinear solver object
 99: .  func - [optional] function evaluation routine, see SNESNewtonTRDCPreCheck()  for the calling sequence
100: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

102:    Level: intermediate

104:    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTRDC solver.

106: .seealso: SNESNewtonTRDCPreCheck(), SNESNewtonTRDCGetPreCheck(), SNESNewtonTRDCSetPostCheck(), SNESNewtonTRDCGetPostCheck()
107: @*/
108: PetscErrorCode  SNESNewtonTRDCSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
109: {
110:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

113:   if (func) tr->precheck    = func;
114:   if (ctx)  tr->precheckctx = ctx;
115:   return 0;
116: }

118: /*@C
119:    SNESNewtonTRDCGetPreCheck - Gets the pre-check function

121:    Not collective

123:    Input Parameter:
124: .  snes - the nonlinear solver context

126:    Output Parameters:
127: +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRDCPreCheck()
128: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

130:    Level: intermediate

132: .seealso: SNESNewtonTRDCSetPreCheck(), SNESNewtonTRDCPreCheck()
133: @*/
134: PetscErrorCode  SNESNewtonTRDCGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
135: {
136:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

139:   if (func) *func = tr->precheck;
140:   if (ctx)  *ctx  = tr->precheckctx;
141:   return 0;
142: }

144: /*@C
145:    SNESNewtonTRDCSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
146:        function evaluation. Allows the user a chance to change or override the decision of the line search routine

148:    Logically Collective on snes

150:    Input Parameters:
151: +  snes - the nonlinear solver object
152: .  func - [optional] function evaluation routine, see SNESNewtonTRDCPostCheck()  for the calling sequence
153: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

155:    Level: intermediate

157:    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTRDC solver while the function set in
158:    SNESLineSearchSetPostCheck() is called AFTER the function evaluation.

160: .seealso: SNESNewtonTRDCPostCheck(), SNESNewtonTRDCGetPostCheck()
161: @*/
162: PetscErrorCode  SNESNewtonTRDCSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
163: {
164:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

167:   if (func) tr->postcheck    = func;
168:   if (ctx)  tr->postcheckctx = ctx;
169:   return 0;
170: }

172: /*@C
173:    SNESNewtonTRDCGetPostCheck - Gets the post-check function

175:    Not collective

177:    Input Parameter:
178: .  snes - the nonlinear solver context

180:    Output Parameters:
181: +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRDCPostCheck()
182: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

184:    Level: intermediate

186: .seealso: SNESNewtonTRDCSetPostCheck(), SNESNewtonTRDCPostCheck()
187: @*/
188: PetscErrorCode  SNESNewtonTRDCGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
189: {
190:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

193:   if (func) *func = tr->postcheck;
194:   if (ctx)  *ctx  = tr->postcheckctx;
195:   return 0;
196: }

198: /*@C
199:    SNESNewtonTRDCPreCheck - Called before the step has been determined in SNESNEWTONTRDC

201:    Logically Collective on snes

203:    Input Parameters:
204: +  snes - the solver
205: .  X - The last solution
206: -  Y - The step direction

208:    Output Parameters:
209: .  changed_Y - Indicator that the step direction Y has been changed.

211:    Level: developer

213: .seealso: SNESNewtonTRDCSetPreCheck(), SNESNewtonTRDCGetPreCheck()
214: @*/
215: static PetscErrorCode SNESNewtonTRDCPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
216: {
217:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

219:   *changed_Y = PETSC_FALSE;
220:   if (tr->precheck) {
221:     (*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx);
223:   }
224:   return 0;
225: }

227: /*@C
228:    SNESNewtonTRDCPostCheck - Called after the step has been determined in SNESNEWTONTRDC but before the function evaluation

230:    Logically Collective on snes

232:    Input Parameters:
233: +  snes - the solver
234: .  X - The last solution
235: .  Y - The full step direction
236: -  W - The updated solution, W = X - Y

238:    Output Parameters:
239: +  changed_Y - indicator if step has been changed
240: -  changed_W - Indicator if the new candidate solution W has been changed.

242:    Notes:
243:      If Y is changed then W is recomputed as X - Y

245:    Level: developer

247: .seealso: SNESNewtonTRDCSetPostCheck(), SNESNewtonTRDCGetPostCheck()
248: @*/
249: static PetscErrorCode SNESNewtonTRDCPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
250: {
251:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;

253:   *changed_Y = PETSC_FALSE;
254:   *changed_W = PETSC_FALSE;
255:   if (tr->postcheck) {
256:     (*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx);
259:   }
260:   return 0;
261: }

263: /*
264:    SNESSolve_NEWTONTRDC - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
265:    (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
266:    nonlinear equations

268: */
269: static PetscErrorCode SNESSolve_NEWTONTRDC(SNES snes)
270: {
271:   SNES_NEWTONTRDC            *neP = (SNES_NEWTONTRDC*)snes->data;
272:   Vec                        X,F,Y,G,W,GradF,YNtmp;
273:   Vec                        YCtmp;
274:   Mat                        jac;
275:   PetscInt                   maxits,i,j,lits,inner_count,bs;
276:   PetscReal                  rho,fnorm,gnorm,xnorm=0,delta,ynorm,temp_xnorm,temp_ynorm;  /* TRDC inner iteration */
277:   PetscReal                  inorms[99]; /* need to make it dynamic eventually, fixed max block size of 99 for now */
278:   PetscReal                  deltaM,ynnorm,f0,mp,gTy,g,yTHy;  /* rho calculation */
279:   PetscReal                  auk,gfnorm,ycnorm,c0,c1,c2,tau,tau_pos,tau_neg,gTBg;  /* Cauchy Point */
280:   KSP                        ksp;
281:   SNESConvergedReason        reason = SNES_CONVERGED_ITERATING;
282:   PetscBool                  breakout = PETSC_FALSE;
283:   SNES_TRDC_KSPConverged_Ctx *ctx;
284:   PetscErrorCode             (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
285:   void                       *convctx;

287:   maxits = snes->max_its;               /* maximum number of iterations */
288:   X      = snes->vec_sol;               /* solution vector */
289:   F      = snes->vec_func;              /* residual vector */
290:   Y      = snes->work[0];               /* update vector */
291:   G      = snes->work[1];               /* updated residual */
292:   W      = snes->work[2];               /* temporary vector */
293:   GradF  = snes->work[3];               /* grad f = J^T F */
294:   YNtmp  = snes->work[4];               /* Newton solution */
295:   YCtmp  = snes->work[5];               /* Cauchy solution */


299:   VecGetBlockSize(YNtmp,&bs);

301:   PetscObjectSAWsTakeAccess((PetscObject)snes);
302:   snes->iter = 0;
303:   PetscObjectSAWsGrantAccess((PetscObject)snes);

305:   /* Set the linear stopping criteria to use the More' trick. From tr.c */
306:   SNESGetKSP(snes,&ksp);
307:   KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy);
308:   if (convtest != SNESTRDC_KSPConverged_Private) {
309:     PetscNew(&ctx);
310:     ctx->snes             = snes;
311:     KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
312:     KSPSetConvergenceTest(ksp,SNESTRDC_KSPConverged_Private,ctx,SNESTRDC_KSPConverged_Destroy);
313:     PetscInfo(snes,"Using Krylov convergence test SNESTRDC_KSPConverged_Private\n");
314:   }

316:   if (!snes->vec_func_init_set) {
317:     SNESComputeFunction(snes,X,F);          /* F(X) */
318:   } else snes->vec_func_init_set = PETSC_FALSE;

320:   VecNorm(F,NORM_2,&fnorm);             /* fnorm <- || F || */
321:   SNESCheckFunctionNorm(snes,fnorm);
322:   VecNorm(X,NORM_2,&xnorm);             /* xnorm <- || X || */

324:   PetscObjectSAWsTakeAccess((PetscObject)snes);
325:   snes->norm = fnorm;
326:   PetscObjectSAWsGrantAccess((PetscObject)snes);
327:   delta      = xnorm ? neP->delta0*xnorm : neP->delta0;  /* initial trust region size scaled by xnorm */
328:   deltaM     = xnorm ? neP->deltaM*xnorm : neP->deltaM;  /* maximum trust region size scaled by xnorm */
329:   neP->delta = delta;
330:   SNESLogConvergenceHistory(snes,fnorm,0);
331:   SNESMonitor(snes,0,fnorm);

333:   neP->rho_satisfied = PETSC_FALSE;

335:   /* test convergence */
336:   (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
337:   if (snes->reason) return 0;

339:   for (i=0; i<maxits; i++) {
340:     PetscBool changed_y;
341:     PetscBool changed_w;

343:     /* dogleg method */
344:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
345:     SNESCheckJacobianDomainerror(snes);
346:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian);
347:     KSPSolve(snes->ksp,F,YNtmp);   /* Quasi Newton Solution */
348:     SNESCheckKSPSolve(snes);  /* this is necessary but old tr.c did not have it*/
349:     KSPGetIterationNumber(snes->ksp,&lits);
350:     SNESGetJacobian(snes,&jac,NULL,NULL,NULL);

352:     /* rescale Jacobian, Newton solution update, and re-calculate delta for multiphase (multivariable)
353:        for inner iteration and Cauchy direction calculation
354:     */
355:     if (bs > 1 && neP->auto_scale_multiphase) {
356:       VecStrideNormAll(YNtmp,NORM_INFINITY,inorms);
357:       for (j=0; j<bs; j++) {
358:         if (neP->auto_scale_max > 1.0) {
359:           if (inorms[j] < 1.0/neP->auto_scale_max) {
360:             inorms[j] = 1.0/neP->auto_scale_max;
361:           }
362:         }
363:         VecStrideSet(W,j,inorms[j]);
364:         VecStrideScale(YNtmp,j,1.0/inorms[j]);
365:         VecStrideScale(X,j,1.0/inorms[j]);
366:       }
367:       VecNorm(X,NORM_2,&xnorm);
368:       if (i == 0) {
369:         delta = neP->delta0*xnorm;
370:       } else {
371:         delta = neP->delta*xnorm;
372:       }
373:       deltaM = neP->deltaM*xnorm;
374:       MatDiagonalScale(jac,PETSC_NULL,W);
375:     }

377:     /* calculating GradF of minimization function */
378:     MatMultTranspose(jac,F,GradF);  /* grad f = J^T F */
379:     VecNorm(YNtmp,NORM_2,&ynnorm);  /* ynnorm <- || Y_newton || */

381:     inner_count = 0;
382:     neP->rho_satisfied = PETSC_FALSE;
383:     while (1) {
384:       if (ynnorm <= delta) {  /* see if the Newton solution is within the trust region */
385:         VecCopy(YNtmp,Y);
386:       } else if (neP->use_cauchy) { /* use Cauchy direction if enabled */
387:         MatMult(jac,GradF,W);
388:         VecDotRealPart(W,W,&gTBg);  /* completes GradF^T J^T J GradF */
389:         VecNorm(GradF,NORM_2,&gfnorm);  /* grad f norm <- || grad f || */
390:         if (gTBg <= 0.0) {
391:           auk = PETSC_MAX_REAL;
392:         } else {
393:           auk = PetscSqr(gfnorm)/gTBg;
394:         }
395:         auk  = PetscMin(delta/gfnorm,auk);
396:         VecCopy(GradF,YCtmp);  /* this could be improved */
397:         VecScale(YCtmp,auk);  /* YCtmp, Cauchy solution*/
398:         VecNorm(YCtmp,NORM_2,&ycnorm);  /* ycnorm <- || Y_cauchy || */
399:         if (ycnorm >= delta) {  /* see if the Cauchy solution meets the criteria */
400:             VecCopy(YCtmp,Y);
401:             PetscInfo(snes,"DL evaluated. delta: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n",(double)delta,(double)ynnorm,(double)ycnorm);
402:         } else {  /* take ratio, tau, of Cauchy and Newton direction and step */
403:           VecAXPY(YNtmp,-1.0,YCtmp);  /* YCtmp = A, YNtmp = B */
404:           VecNorm(YNtmp,NORM_2,&c0);  /* this could be improved */
405:           c0      = PetscSqr(c0);
406:           VecDotRealPart(YCtmp,YNtmp,&c1);
407:           c1      = 2.0*c1;
408:           VecNorm(YCtmp,NORM_2,&c2);  /* this could be improved */
409:           c2      = PetscSqr(c2) - PetscSqr(delta);
410:           tau_pos = (c1 + PetscSqrtReal(PetscSqr(c1) - 4.*c0*c2))/(2.*c0); /* quadratic formula */
411:           tau_neg = (c1 - PetscSqrtReal(PetscSqr(c1) - 4.*c0*c2))/(2.*c0);
412:           tau     = PetscMax(tau_pos, tau_neg);  /* can tau_neg > tau_pos? I don't think so, but just in case. */
413:           PetscInfo(snes,"DL evaluated. tau: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n",(double)tau,(double)ynnorm,(double)ycnorm);
414:           VecWAXPY(W,tau,YNtmp,YCtmp);
415:           VecAXPY(W,-tau,YCtmp);
416:           VecCopy(W, Y); /* this could be improved */
417:         }
418:       } else {
419:         /* if Cauchy is disabled, only use Newton direction */
420:         auk = delta/ynnorm;
421:         VecScale(YNtmp,auk);
422:         VecCopy(YNtmp,Y); /* this could be improved (many VecCopy, VecNorm)*/
423:       }

425:       VecNorm(Y,NORM_2,&ynorm);  /* compute the final ynorm  */
426:       f0 = 0.5*PetscSqr(fnorm);  /* minimizing function f(X) */
427:       MatMult(jac,Y,W);
428:       VecDotRealPart(W,W,&yTHy);  /* completes GradY^T J^T J GradY */
429:       VecDotRealPart(GradF,Y,&gTy);
430:       mp = f0 - gTy + 0.5*yTHy;  /* quadratic model to satisfy, -gTy because our update is X-Y*/

432:       /* scale back solution update */
433:       if (bs > 1 && neP->auto_scale_multiphase) {
434:         for (j=0; j<bs; j++) {
435:           VecStrideScale(Y,j,inorms[j]);
436:           if (inner_count == 0) {
437:             /* TRDC inner algorithm does not need scaled X after calculating delta in the outer iteration */
438:             /* need to scale back X to match Y and provide proper update to the external code */
439:             VecStrideScale(X,j,inorms[j]);
440:           }
441:         }
442:         if (inner_count == 0) VecNorm(X,NORM_2,&temp_xnorm);  /* only in the first iteration */
443:         VecNorm(Y,NORM_2,&temp_ynorm);
444:       } else {
445:         temp_xnorm = xnorm;
446:         temp_ynorm = ynorm;
447:       }
448:       inner_count++;

450:       /* Evaluate the solution to meet the improvement ratio criteria */
451:       SNESNewtonTRDCPreCheck(snes,X,Y,&changed_y);
452:       VecWAXPY(W,-1.0,Y,X);
453:       SNESNewtonTRDCPostCheck(snes,X,Y,W,&changed_y,&changed_w);
454:       if (changed_y) VecWAXPY(W,-1.0,Y,X);
455:       VecCopy(Y,snes->vec_sol_update);
456:       SNESComputeFunction(snes,W,G); /*  F(X-Y) = G */
457:       VecNorm(G,NORM_2,&gnorm);      /* gnorm <- || g || */
458:       SNESCheckFunctionNorm(snes,gnorm);
459:       g = 0.5*PetscSqr(gnorm); /* minimizing function g(W) */
460:       if (f0 == mp) rho = 0.0;
461:       else rho = (f0 - g)/(f0 - mp);  /* actual improvement over predicted improvement */

463:       if (rho < neP->eta2) {
464:         delta *= neP->t1;  /* shrink the region */
465:       } else if (rho > neP->eta3) {
466:         delta = PetscMin(neP->t2*delta,deltaM); /* expand the region, but not greater than deltaM */
467:       }

469:       neP->delta = delta;
470:       if (rho >= neP->eta1) {
471:         /* unscale delta and xnorm before going to the next outer iteration */
472:         if (bs > 1 && neP->auto_scale_multiphase) {
473:           neP->delta = delta/xnorm;
474:           xnorm      = temp_xnorm;
475:           ynorm      = temp_ynorm;
476:         }
477:         neP->rho_satisfied = PETSC_TRUE;
478:         break;  /* the improvement ratio is satisfactory */
479:       }
480:       PetscInfo(snes,"Trying again in smaller region\n");

482:       /* check to see if progress is hopeless */
483:       neP->itflag = PETSC_FALSE;
484:       /* both delta, ynorm, and xnorm are either scaled or unscaled */
485:       SNESTRDC_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
486:       if (!reason) {
487:          /* temp_xnorm, temp_ynorm is always unscaled */
488:          /* also the inner iteration already calculated the Jacobian and solved the matrix */
489:          /* therefore, it should be passing iteration number of iter+1 instead of iter+0 in the first iteration and after */
490:          (*snes->ops->converged)(snes,snes->iter+1,temp_xnorm,temp_ynorm,fnorm,&reason,snes->cnvP);
491:       }
492:       /* if multiphase state changes, break out inner iteration */
493:       if (reason == SNES_BREAKOUT_INNER_ITER) {
494:         if (bs > 1 && neP->auto_scale_multiphase) {
495:           /* unscale delta and xnorm before going to the next outer iteration */
496:           neP->delta = delta/xnorm;
497:           xnorm      = temp_xnorm;
498:           ynorm      = temp_ynorm;
499:         }
500:         reason = SNES_CONVERGED_ITERATING;
501:         break;
502:       }
503:       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
504:       if (reason) {
505:         if (reason < 0) {
506:             /* We're not progressing, so return with the current iterate */
507:             SNESMonitor(snes,i+1,fnorm);
508:             breakout = PETSC_TRUE;
509:             break;
510:         } else if (reason > 0) {
511:             /* We're converged, so return with the current iterate and update solution */
512:             SNESMonitor(snes,i+1,fnorm);
513:             breakout = PETSC_FALSE;
514:             break;
515:         }
516:       }
517:       snes->numFailures++;
518:     }
519:     if (!breakout) {
520:       /* Update function and solution vectors */
521:       fnorm       = gnorm;
522:       VecCopy(G,F);
523:       VecCopy(W,X);
524:       /* Monitor convergence */
525:       PetscObjectSAWsTakeAccess((PetscObject)snes);
526:       snes->iter  = i+1;
527:       snes->norm  = fnorm;
528:       snes->xnorm = xnorm;
529:       snes->ynorm = ynorm;
530:       PetscObjectSAWsGrantAccess((PetscObject)snes);
531:       SNESLogConvergenceHistory(snes,snes->norm,lits);
532:       SNESMonitor(snes,snes->iter,snes->norm);
533:       /* Test for convergence, xnorm = || X || */
534:       neP->itflag = PETSC_TRUE;
535:       if (snes->ops->converged != SNESConvergedSkip) VecNorm(X,NORM_2,&xnorm);
536:       (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
537:       if (reason) break;
538:     } else break;
539:   }

541:   /* PetscFree(inorms); */
542:   if (i == maxits) {
543:     PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits);
544:     if (!reason) reason = SNES_DIVERGED_MAX_IT;
545:   }
546:   PetscObjectSAWsTakeAccess((PetscObject)snes);
547:   snes->reason = reason;
548:   PetscObjectSAWsGrantAccess((PetscObject)snes);
549:   if (convtest != SNESTRDC_KSPConverged_Private) {
550:     KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
551:     PetscFree(ctx);
552:     KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy);
553:   }
554:   return 0;
555: }

557: /*------------------------------------------------------------*/
558: static PetscErrorCode SNESSetUp_NEWTONTRDC(SNES snes)
559: {
560:   SNESSetWorkVecs(snes,6);
561:   SNESSetUpMatrices(snes);
562:   return 0;
563: }

565: PetscErrorCode SNESReset_NEWTONTRDC(SNES snes)
566: {
567:   return 0;
568: }

570: static PetscErrorCode SNESDestroy_NEWTONTRDC(SNES snes)
571: {
572:   SNESReset_NEWTONTRDC(snes);
573:   PetscFree(snes->data);
574:   return 0;
575: }
576: /*------------------------------------------------------------*/

578: static PetscErrorCode SNESSetFromOptions_NEWTONTRDC(PetscOptionItems *PetscOptionsObject,SNES snes)
579: {
580:   SNES_NEWTONTRDC  *ctx = (SNES_NEWTONTRDC*)snes->data;

582:   PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
583:   PetscOptionsReal("-snes_trdc_tol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
584:   PetscOptionsReal("-snes_trdc_eta1","eta1","None",ctx->eta1,&ctx->eta1,NULL);
585:   PetscOptionsReal("-snes_trdc_eta2","eta2","None",ctx->eta2,&ctx->eta2,NULL);
586:   PetscOptionsReal("-snes_trdc_eta3","eta3","None",ctx->eta3,&ctx->eta3,NULL);
587:   PetscOptionsReal("-snes_trdc_t1","t1","None",ctx->t1,&ctx->t1,NULL);
588:   PetscOptionsReal("-snes_trdc_t2","t2","None",ctx->t2,&ctx->t2,NULL);
589:   PetscOptionsReal("-snes_trdc_deltaM","deltaM","None",ctx->deltaM,&ctx->deltaM,NULL);
590:   PetscOptionsReal("-snes_trdc_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
591:   PetscOptionsReal("-snes_trdc_auto_scale_max","auto_scale_max","None",ctx->auto_scale_max,&ctx->auto_scale_max,NULL);
592:   PetscOptionsBool("-snes_trdc_use_cauchy","use_cauchy","use Cauchy step and direction",ctx->use_cauchy,&ctx->use_cauchy,NULL);
593:   PetscOptionsBool("-snes_trdc_auto_scale_multiphase","auto_scale_multiphase","Auto scaling for proper cauchy direction",ctx->auto_scale_multiphase,&ctx->auto_scale_multiphase,NULL);
594:   PetscOptionsTail();
595:   return 0;
596: }

598: static PetscErrorCode SNESView_NEWTONTRDC(SNES snes,PetscViewer viewer)
599: {
600:   SNES_NEWTONTRDC  *tr = (SNES_NEWTONTRDC*)snes->data;
601:   PetscBool        iascii;

603:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
604:   if (iascii) {
605:     PetscViewerASCIIPrintf(viewer,"  Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
606:     PetscViewerASCIIPrintf(viewer,"  eta1=%g, eta2=%g, eta3=%g\n",(double)tr->eta1,(double)tr->eta2,(double)tr->eta3);
607:     PetscViewerASCIIPrintf(viewer,"  delta0=%g, t1=%g, t2=%g, deltaM=%g\n",(double)tr->delta0,(double)tr->t1,(double)tr->t2,(double)tr->deltaM);
608:   }
609:   return 0;
610: }
611: /* ------------------------------------------------------------ */
612: /*MC
613:       SNESNEWTONTRDC - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction

615:    Options Database:
616: +   -snes_trdc_tol <tol> - trust region tolerance
617: .   -snes_trdc_eta1 <eta1> - trust region parameter 0.0 <= eta1 <= eta2, rho >= eta1 breaks out of the inner iteration (default: eta1=0.001)
618: .   -snes_trdc_eta2 <eta2> - trust region parameter 0.0 <= eta1 <= eta2, rho <= eta2 shrinks the trust region (default: eta2=0.25)
619: .   -snes_trdc_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
620: .   -snes_trdc_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
621: .   -snes_trdc_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
622: .   -snes_trdc_deltaM <deltaM> - trust region parameter, max size of trust region, deltaM*norm2(x) (default: 0.5)
623: .   -snes_trdc_delta0 <delta0> - trust region parameter, initial size of trust region, delta0*norm2(x) (default: 0.1)
624: .   -snes_trdc_auto_scale_max <auto_scale_max> - used with auto_scale_multiphase, caps the maximum auto-scaling factor
625: .   -snes_trdc_use_cauchy <use_cauchy> - True uses dogleg Cauchy (Steepest Descent direction) step & direction in the trust region algorithm
626: -   -snes_trdc_auto_scale_multiphase <auto_scale_multiphase> - True turns on auto-scaling for multivariable block matrix for Cauchy and trust region

628:     Notes:
629:     The algorithm is taken from "Linear and Nonlinear Solvers for Simulating Multiphase Flow
630:     within Large-Scale Engineered Subsurface Systems" by Heeho D. Park, Glenn E. Hammond,
631:     Albert J. Valocchi, Tara LaForce.

633:    Level: intermediate

635: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance(), SNESNEWTONTRDC

637: M*/
638: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTRDC(SNES snes)
639: {
640:   SNES_NEWTONTRDC  *neP;

642:   snes->ops->setup          = SNESSetUp_NEWTONTRDC;
643:   snes->ops->solve          = SNESSolve_NEWTONTRDC;
644:   snes->ops->destroy        = SNESDestroy_NEWTONTRDC;
645:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTRDC;
646:   snes->ops->view           = SNESView_NEWTONTRDC;
647:   snes->ops->reset          = SNESReset_NEWTONTRDC;

649:   snes->usesksp = PETSC_TRUE;
650:   snes->usesnpc = PETSC_FALSE;

652:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

654:   PetscNewLog(snes,&neP);
655:   snes->data  = (void*)neP;
656:   neP->delta  = 0.0;
657:   neP->delta0 = 0.1;
658:   neP->eta1   = 0.001;
659:   neP->eta2   = 0.25;
660:   neP->eta3   = 0.75;
661:   neP->t1     = 0.25;
662:   neP->t2     = 2.0;
663:   neP->deltaM = 0.5;
664:   neP->sigma  = 0.0001;
665:   neP->itflag = PETSC_FALSE;
666:   neP->rnorm0 = 0.0;
667:   neP->ttol   = 0.0;
668:   neP->use_cauchy            = PETSC_TRUE;
669:   neP->auto_scale_multiphase = PETSC_FALSE;
670:   neP->auto_scale_max        = -1.0;
671:   neP->rho_satisfied         = PETSC_FALSE;
672:   snes->deltatol             = 1.e-12;

674:   /* for multiphase (multivariable) scaling */
675:   /* may be used for dynamic allocation of inorms, but it fails snes_tutorials-ex3_13
676:      on test forced DIVERGED_JACOBIAN_DOMAIN test. I will use static array for now.
677:   VecGetBlockSize(snes->work[0],&neP->bs);
678:   PetscCalloc1(neP->bs,&neP->inorms);
679:   */

681:   return 0;
682: }