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

 36: static PetscErrorCode SNESTRDC_KSPConverged_Destroy(void *cctx)
 37: {
 38:   SNES_TRDC_KSPConverged_Ctx *ctx = (SNES_TRDC_KSPConverged_Ctx *)cctx;

 40:   (*ctx->convdestroy)(ctx->convctx);
 41:   PetscFree(ctx);

 43:   return 0;
 44: }

 46: /*
 47:    SNESTRDC_Converged_Private -test convergence JUST for
 48:    the trust region tolerance.

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

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

 66: /*@
 67:   SNESNewtonTRDCGetRhoFlag - Get whether the current solution update is within the trust-region.

 69:   Input Parameter:
 70: . snes - the nonlinear solver object

 72:   Output Parameter:
 73: . rho_flag: `PETSC_TRUE` if the solution update is in the trust-region; otherwise, `PETSC_FALSE`

 75:   Level: developer

 77: .seealso: `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, , `SNESNewtonTRDCSetPreCheck()`,
 78:           `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`
 79: @*/
 80: PetscErrorCode SNESNewtonTRDCGetRhoFlag(SNES snes, PetscBool *rho_flag)
 81: {
 82:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

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

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

 94:    Logically Collective

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

101:    Level: intermediate

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

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

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

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

122:    Not collective

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

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

131:    Level: intermediate

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

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

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

149:    Logically Collective

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

156:    Level: intermediate

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

162: .seealso: `SNESNEWTONTRDC`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
163: @*/
164: PetscErrorCode SNESNewtonTRDCSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
165: {
166:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

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

174: /*@C
175:    SNESNewtonTRDCGetPostCheck - Gets the post-check function

177:    Not collective

179:    Input Parameter:
180: .  snes - the nonlinear solver context

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

186:    Level: intermediate

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

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

200: /*@C
201:    SNESNewtonTRDCPreCheck - Called before the step has been determined in `SNESNEWTONTRDC`

203:    Logically Collective

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

210:    Output Parameters:
211: .  changed_Y - Indicator that the step direction Y has been changed.

213:    Level: developer

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

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

229: /*@C
230:    SNESNewtonTRDCPostCheck - Called after the step has been determined in `SNESNEWTONTRDC` but before the function evaluation at that step

232:    Logically Collective

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

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

244:    Note:
245:      If Y is changed then W is recomputed as X - Y

247:    Level: developer

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

255:   *changed_Y = PETSC_FALSE;
256:   *changed_W = PETSC_FALSE;
257:   if (tr->postcheck) {
258:     (*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx);
261:   }
262:   return 0;
263: }

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

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

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


301:   VecGetBlockSize(YNtmp, &bs);

303:   PetscObjectSAWsTakeAccess((PetscObject)snes);
304:   snes->iter = 0;
305:   PetscObjectSAWsGrantAccess((PetscObject)snes);

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

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

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

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

335:   neP->rho_satisfied = PETSC_FALSE;

337:   /* test convergence */
338:   PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
339:   if (snes->reason) return 0;

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

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

354:     /* rescale Jacobian, Newton solution update, and re-calculate delta for multiphase (multivariable)
355:        for inner iteration and Cauchy direction calculation
356:     */
357:     if (bs > 1 && neP->auto_scale_multiphase) {
358:       VecStrideNormAll(YNtmp, NORM_INFINITY, inorms);
359:       for (j = 0; j < bs; j++) {
360:         if (neP->auto_scale_max > 1.0) {
361:           if (inorms[j] < 1.0 / neP->auto_scale_max) inorms[j] = 1.0 / neP->auto_scale_max;
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:       PetscUseTypeMethod(snes, converged, 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: static PetscErrorCode SNESSetUp_NEWTONTRDC(SNES snes)
558: {
559:   SNESSetWorkVecs(snes, 6);
560:   SNESSetUpMatrices(snes);
561:   return 0;
562: }

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

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

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

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

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

601:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
602:   if (iascii) {
603:     PetscViewerASCIIPrintf(viewer, "  Trust region tolerance %g (-snes_trtol)\n", (double)snes->deltatol);
604:     PetscViewerASCIIPrintf(viewer, "  eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3);
605:     PetscViewerASCIIPrintf(viewer, "  delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM);
606:   }
607:   return 0;
608: }

610: /*MC
611:       SNESNEWTONTRDC - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction

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

626:     Reference:
627: .   - * "Linear and Nonlinear Solvers for Simulating Multiphase Flow
628:     within Large-Scale Engineered Subsurface Systems" by Heeho D. Park, Glenn E. Hammond,
629:     Albert J. Valocchi, Tara LaForce.

631:    Level: intermediate

633: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`,
634:           `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
635:           `SNESNewtonTRDCGetRhoFlag()`, `SNESNewtonTRDCSetPreCheck()`
636: M*/
637: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTRDC(SNES snes)
638: {
639:   SNES_NEWTONTRDC *neP;

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

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

651:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

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

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

680:   return 0;
681: }