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: }