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