Actual source code: al.c
1: #include <../src/snes/impls/al/alimpl.h>
3: /*
4: This file implements a truncated Newton method with arc length continuation,
5: for solving a system of nonlinear equations, using the KSP, Vec,
6: and Mat interfaces for linear solvers, vectors, and matrices,
7: respectively.
8: */
10: const char NewtonALExactCitation[] = "@article{Ritto-CorreaCamotim2008,\n"
11: " title={On the arc-length and other quadratic control methods: Established, less known and new implementation procedures},\n"
12: " volume={86},\n"
13: " ISSN={0045-7949},\n"
14: " DOI={10.1016/j.compstruc.2007.08.003},\n"
15: " number={11},\n"
16: " journal={Computers & Structures},\n"
17: " author={Ritto-Corr{\\^{e}}a, Manuel and Camotim, Dinar},\n"
18: " year={2008},\n"
19: " month=jun,\n"
20: " pages={1353-1368},\n"
21: "}\n";
22: PetscBool NewtonALExactCitationSet = PETSC_FALSE;
23: const char NewtonALNormalCitation[] = "@article{LeonPaulinoPereiraMenezesLages_2011,\n"
24: " title={A Unified Library of Nonlinear Solution Schemes},\n"
25: " volume={64},\n"
26: " ISSN={0003-6900, 2379-0407},\n"
27: " DOI={10.1115/1.4006992},\n"
28: " number={4},\n"
29: " journal={Applied Mechanics Reviews},\n"
30: " author={Leon, Sofie E. and Paulino, Glaucio H. and Pereira, Anderson and Menezes, Ivan F. M. and Lages, Eduardo N.},\n"
31: " year={2011},\n"
32: " month=jul,\n"
33: " pages={040803},\n"
34: " language={en}\n"
35: "}\n";
36: PetscBool NewtonALNormalCitationSet = PETSC_FALSE;
38: const char *const SNESNewtonALCorrectionTypes[] = {"EXACT", "NORMAL", "SNESNewtonALCorrectionType", "SNES_NEWTONAL_CORRECTION_", NULL};
40: static PetscErrorCode SNESNewtonALCheckArcLength(SNES snes, Vec XStep, PetscReal lambdaStep, PetscReal stepSize)
41: {
42: PetscReal arcLength, arcLengthError;
43: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
45: PetscFunctionBegin;
46: if (al->mat_diag_scaling) PetscCall(MatANorm(al->mat_diag_scaling, XStep, &arcLength));
47: else PetscCall(VecNorm(XStep, NORM_2, &arcLength));
48: arcLength = PetscSqrtReal(arcLength * arcLength + al->psisq * lambdaStep * lambdaStep);
49: arcLengthError = PetscAbsReal(arcLength - stepSize);
51: if (arcLengthError > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscInfo(snes, "Arc length differs from specified step size: computed=%18.16e, expected=%18.16e, error=%18.16e \n", (double)arcLength, (double)stepSize, (double)arcLengthError));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
56: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
57: {
58: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
59: PetscReal x1 = temp / a;
60: PetscReal x2 = c / temp;
61: *xm = PetscMin(x1, x2);
62: *xp = PetscMax(x1, x2);
63: }
65: static PetscErrorCode SNESNewtonALSetCorrectionType_NEWTONAL(SNES snes, SNESNewtonALCorrectionType ctype)
66: {
67: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
69: PetscFunctionBegin;
70: al->correction_type = ctype;
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*@
75: SNESNewtonALSetCorrectionType - Set the type of correction to use in the arc-length continuation method.
77: Logically Collective
79: Input Parameters:
80: + snes - the nonlinear solver object
81: - ctype - the type of correction to use
83: Options Database Key:
84: . -snes_newtonal_correction_type type - Set the type of correction to use; use -help for a list of available types
86: Level: intermediate
88: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALCorrectionType`
89: @*/
90: PetscErrorCode SNESNewtonALSetCorrectionType(SNES snes, SNESNewtonALCorrectionType ctype)
91: {
92: PetscFunctionBegin;
95: PetscTryMethod(snes, "SNESNewtonALSetCorrectionType_C", (SNES, SNESNewtonALCorrectionType), (snes, ctype));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode SNESNewtonALSetFunction_NEWTONAL(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
100: {
101: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
103: PetscFunctionBegin;
104: al->computealfunction = func;
105: al->alctx = ctx;
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*@C
110: SNESNewtonALSetFunction - Sets a user function that is called at each function evaluation to
111: compute the tangent load vector for the arc-length continuation method.
113: Logically Collective
115: Input Parameters:
116: + snes - the nonlinear solver object
117: . func - [optional] tangent load function evaluation routine, see `SNESFunctionFn` for the calling sequence. `U` is the current solution vector, `Q` is the output tangent load vector
118: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
120: Level: intermediate
122: Notes:
123: If the current value of the load parameter is needed in `func`, it can be obtained with `SNESNewtonALGetLoadParameter()`.
125: The tangent load vector is the partial derivative of external load with respect to the load parameter.
126: In the case of proportional loading, the tangent load vector is the full external load vector at the end of the load step.
128: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
129: @*/
130: PetscErrorCode SNESNewtonALSetFunction(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
131: {
132: PetscFunctionBegin;
134: PetscTryMethod(snes, "SNESNewtonALSetFunction_C", (SNES, SNESFunctionFn *, void *), (snes, func, ctx));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode SNESNewtonALGetFunction_NEWTONAL(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
139: {
140: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
142: PetscFunctionBegin;
143: if (func) *func = al->computealfunction;
144: if (ctx) *(void **)ctx = al->alctx;
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*@C
149: SNESNewtonALGetFunction - Get the user function and context set with `SNESNewtonALSetFunction`
151: Logically Collective
153: Input Parameters:
154: + snes - the nonlinear solver object
155: . func - [optional] tangent load function evaluation routine, see `SNESNewtonALSetFunction()` for the call sequence
156: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
158: Level: intermediate
160: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
161: @*/
162: PetscErrorCode SNESNewtonALGetFunction(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
163: {
164: PetscFunctionBegin;
166: PetscUseMethod(snes, "SNESNewtonALGetFunction_C", (SNES, SNESFunctionFn **, PetscCtxRt), (snes, func, ctx));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*@C
171: SNESNewtonALSetDiagonalScaling - Set the global vector used to rescale DoFs for computation of arc length.
173: Logically Collective
175: Input Parameters:
176: + snes - the nonlinear solver object
177: - v - the `Vec` containing diagonal scaling for each DoF, must be the same size as the solution vector (may be `NULL`)
179: Note:
180: This function stores a reference to `v`. Any changes to the vector will be reflected automatically in the arc length computation.
182: Level: intermediate
184: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetLoadParameter()`
185: @*/
186: PetscErrorCode SNESNewtonALSetDiagonalScaling(SNES snes, Vec v)
187: {
188: PetscFunctionBegin;
191: PetscTryMethod(snes, "SNESNewtonALSetDiagonalScaling_C", (SNES, Vec), (snes, v));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: static PetscErrorCode SNESNewtonALSetDiagonalScaling_NEWTONAL(SNES snes, Vec v)
196: {
197: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
199: PetscFunctionBegin;
200: PetscCall(PetscObjectReference((PetscObject)v));
201: PetscCall(MatDestroy(&al->mat_diag_scaling));
202: if (v) {
203: PetscCall(MatCreateDiagonal(v, &al->mat_diag_scaling));
204: PetscCall(PetscObjectDereference((PetscObject)v));
205: }
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes, PetscReal *lambda)
210: {
211: SNES_NEWTONAL *al;
213: PetscFunctionBeginHot;
214: al = (SNES_NEWTONAL *)snes->data;
215: *lambda = al->lambda;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@C
220: SNESNewtonALGetLoadParameter - Get the value of the load parameter `lambda` for the arc-length continuation method.
222: Logically Collective
224: Input Parameter:
225: . snes - the nonlinear solver object
227: Output Parameter:
228: . lambda - the arc-length parameter
230: Level: intermediate
232: Notes:
233: This function should be used in the functions provided to `SNESSetFunction()` and `SNESNewtonALSetFunction()`
234: to compute the residual and tangent load vectors for a given value of `lambda` (0 <= lambda <= 1).
236: Usually, `lambda` is used to scale the external force vector in the residual function, i.e. proportional loading,
237: in which case the tangent load vector is the full external force vector.
239: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
240: @*/
241: PetscErrorCode SNESNewtonALGetLoadParameter(SNES snes, PetscReal *lambda)
242: {
243: PetscFunctionBeginHot;
245: PetscAssertPointer(lambda, 2);
246: PetscUseMethod(snes, "SNESNewtonALGetLoadParameter_C", (SNES, PetscReal *), (snes, lambda));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode SNESNewtonALComputeFunction_NEWTONAL(SNES snes, Vec X, Vec Q)
251: {
252: void *ctx = NULL;
253: SNESFunctionFn *computealfunction = NULL;
254: SNES_NEWTONAL *al;
256: PetscFunctionBegin;
257: al = (SNES_NEWTONAL *)snes->data;
258: PetscCall(SNESNewtonALGetFunction(snes, &computealfunction, &ctx));
260: PetscCall(VecZeroEntries(Q));
261: PetscCheck(computealfunction || (snes->vec_rhs && al->scale_rhs), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No tangent load function or rhs vector has been set");
262: if (computealfunction) {
263: PetscCall(VecLockReadPush(X));
264: PetscCallBack("SNES callback NewtonAL tangent load function", (*computealfunction)(snes, X, Q, ctx));
265: PetscCall(VecLockReadPop(X));
266: }
267: if (al->scale_rhs && snes->vec_rhs) {
268: /* Save original RHS vector values, then scale `snes->vec_rhs` by load parameter */
269: if (!al->vec_rhs_orig) PetscCall(VecDuplicate(snes->vec_rhs, &al->vec_rhs_orig));
270: if (!al->copied_rhs) {
271: PetscCall(VecCopy(snes->vec_rhs, al->vec_rhs_orig));
272: al->copied_rhs = PETSC_TRUE;
273: }
274: PetscCall(VecAXPBY(snes->vec_rhs, al->lambda, 0.0, al->vec_rhs_orig));
275: PetscCall(VecAXPY(Q, 1, al->vec_rhs_orig));
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@C
281: SNESNewtonALComputeFunction - Calls the function that has been set with `SNESNewtonALSetFunction()`.
283: Collective
285: Input Parameters:
286: + snes - the `SNES` context
287: - X - input vector
289: Output Parameter:
290: . Q - tangent load vector, as set by `SNESNewtonALSetFunction()`
292: Level: developer
294: Note:
295: `SNESNewtonALComputeFunction()` is typically used within nonlinear solvers
296: implementations, so users would not generally call this routine themselves.
298: .seealso: [](ch_snes), `SNES`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`
299: @*/
300: PetscErrorCode SNESNewtonALComputeFunction(SNES snes, Vec X, Vec Q)
301: {
302: PetscFunctionBegin;
306: PetscCheckSameComm(snes, 1, X, 2);
307: PetscCheckSameComm(snes, 1, Q, 3);
308: PetscCall(VecValidValues_Internal(X, 2, PETSC_TRUE));
309: PetscCall(PetscLogEventBegin(SNES_NewtonALEval, snes, X, Q, 0));
310: PetscTryMethod(snes, "SNESNewtonALComputeFunction_C", (SNES, Vec, Vec), (snes, X, Q));
311: PetscCall(PetscLogEventEnd(SNES_NewtonALEval, snes, X, Q, 0));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*
316: SNESSolve_NEWTONAL - Solves a nonlinear system with Newton's method with arc length continuation.
318: Input Parameter:
319: . snes - the `SNES` context
321: Application Interface Routine: SNESSolve()
322: */
323: static PetscErrorCode SNESSolve_NEWTONAL(SNES snes)
324: {
325: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
326: PetscInt maxits, maxincs, lits;
327: PetscReal fnorm, xnorm, ynorm, stepSize;
328: Vec DeltaX, deltaX, X, R, Q, deltaX_Q, deltaX_R;
329: Mat W;
331: PetscFunctionBegin;
332: PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
334: /* Register citations */
335: PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
336: if (data->correction_type == SNES_NEWTONAL_CORRECTION_EXACT) {
337: PetscCall(PetscCitationsRegister(NewtonALExactCitation, &NewtonALExactCitationSet));
338: } else if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
339: PetscCall(PetscCitationsRegister(NewtonALNormalCitation, &NewtonALNormalCitationSet));
340: }
342: data->copied_rhs = PETSC_FALSE;
343: data->lambda_update = 0.0;
344: data->lambda = 0.0;
345: snes->numFailures = 0;
346: snes->numLinearSolveFailures = 0;
347: snes->reason = SNES_CONVERGED_ITERATING;
348: snes->iter = 0;
350: maxits = snes->max_its; /* maximum number of iterations */
351: maxincs = data->max_continuation_steps; /* maximum number of increments */
352: X = snes->vec_sol; /* solution vector */
353: R = snes->vec_func; /* residual vector */
354: Q = snes->work[0]; /* tangent load vector */
355: W = data->mat_diag_scaling; /* diagonal scaling for DoFs */
356: deltaX_Q = snes->work[1]; /* variation of X with respect to lambda */
357: deltaX_R = snes->work[2]; /* linearized error correction */
358: DeltaX = snes->work[3]; /* step from equilibrium */
359: deltaX = snes->vec_sol_update; /* full newton step */
360: stepSize = data->step_size; /* initial step size */
362: PetscCall(VecZeroEntries(DeltaX));
364: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
365: /* set snes->max_its for convergence test */
366: snes->max_its = maxits * maxincs;
367: snes->iter = 0;
368: snes->norm = 0.0;
369: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
371: /* main incremental-iterative loop */
372: for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
373: PetscReal deltaLambda;
375: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
376: snes->norm = 0.0;
377: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
378: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
379: PetscCall(SNESComputeFunction(snes, X, R));
380: PetscCall(VecAXPY(R, 1, Q)); /* R <- R + Q */
381: PetscCall(VecNorm(R, NORM_2, &fnorm)); /* fnorm <- ||R|| */
382: SNESCheckFunctionDomainError(snes, fnorm);
384: /* Monitor convergence */
385: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
386: PetscCall(SNESMonitor(snes, snes->iter, fnorm));
387: if (snes->reason) break;
388: for (PetscInt j = 0; j < maxits; j++) {
389: PetscReal normsqX_Q, deltaS = 1;
391: /* Call general purpose update function */
392: PetscTryTypeMethod(snes, update, snes->iter);
394: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
395: SNESCheckJacobianDomainError(snes);
396: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
397: /* Solve J deltaX_Q = Q, where J is Jacobian matrix */
398: PetscCall(KSPSolve(snes->ksp, Q, deltaX_Q));
399: SNESCheckKSPSolve(snes);
400: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
401: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", tangent load linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
402: /* Compute load parameter variation */
403: if (W) PetscCall(MatANorm(W, deltaX_Q, &normsqX_Q));
404: else PetscCall(VecNorm(deltaX_Q, NORM_2, &normsqX_Q));
405: normsqX_Q *= normsqX_Q;
406: /* On first iter, use predictor. This is the same regardless of corrector scheme. */
407: if (j == 0) {
408: PetscReal sign = 1.0;
409: if (i > 0) {
410: PetscScalar dot;
411: if (W) PetscCall(MatADot(W, DeltaX, deltaX_Q, &dot));
412: else PetscCall(VecDot(DeltaX, deltaX_Q, &dot));
413: sign = PetscRealPart(dot) + data->psisq * data->lambda_update;
414: sign = sign >= 0 ? 1.0 : -1.0;
415: }
416: data->lambda_update = 0.0;
417: PetscCall(VecZeroEntries(DeltaX));
418: deltaLambda = sign * stepSize / PetscSqrtReal(normsqX_Q + data->psisq);
419: } else {
420: /* Solve J deltaX_R = -R */
421: PetscCall(KSPSolve(snes->ksp, R, deltaX_R));
422: SNESCheckKSPSolve(snes);
423: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
424: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", residual linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
425: PetscCall(VecScale(deltaX_R, -1));
427: if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
428: /*
429: Take a step orthogonal to the current incremental update DeltaX.
430: Note, this approach is cheaper than the exact correction, but may exhibit convergence
431: issues due to the iterative trial points not being on the quadratic constraint surface.
432: On the bright side, we always have a real and unique solution for deltaLambda.
433: */
434: PetscScalar coefs[2];
435: if (W) {
436: PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs[0]));
437: PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs[1]));
438: } else {
439: PetscCall(VecDot(DeltaX, deltaX_R, &coefs[0]));
440: PetscCall(VecDot(DeltaX, deltaX_Q, &coefs[1]));
441: }
442: deltaLambda = -PetscRealPart(coefs[0]) / (PetscRealPart(coefs[1]) + data->psisq * data->lambda_update);
443: } else {
444: /*
445: Solve
446: a*deltaLambda^2 + b*deltaLambda + c = 0 (*)
447: where
448: a = a0
449: b = b0 + b1*deltaS
450: c = c0 + c1*deltaS + c2*deltaS^2
451: and deltaS is either 1, or the largest value in (0, 1) that satisfies
452: b^2 - 4*a*c = as*deltaS^2 + bs*deltaS + cs >= 0
453: where
454: as = b1^2 - 4*a0*c2
455: bs = 2*b1*b0 - 4*a0*c1
456: cs = b0^2 - 4*a0*c0
457: These "partial corrections" prevent (*) from having complex roots.
458: */
459: PetscReal psisqLambdaUpdate, discriminant;
460: PetscReal as, bs, cs;
461: PetscReal a0, b0, b1, c0, c1, c2;
462: PetscScalar coefs1[3]; /* coefs[0] = deltaX_Q*DeltaX, coefs[1] = deltaX_R*DeltaX, coefs[2] = DeltaX*DeltaX */
463: PetscScalar coefs2[2]; /* coefs[0] = deltaX_Q*deltaX_R, coefs[1] = deltaX_R*deltaX_R */
465: psisqLambdaUpdate = data->psisq * data->lambda_update;
466: if (W) {
467: PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs1[0]));
468: PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs1[1]));
469: PetscCall(MatADot(W, DeltaX, DeltaX, &coefs1[2]));
470: PetscCall(MatADot(W, deltaX_R, deltaX_Q, &coefs2[0]));
471: PetscCall(MatADot(W, deltaX_R, deltaX_R, &coefs2[1]));
472: } else {
473: PetscCall(VecDot(DeltaX, deltaX_Q, &coefs1[0]));
474: PetscCall(VecDot(DeltaX, deltaX_R, &coefs1[1]));
475: PetscCall(VecDot(DeltaX, DeltaX, &coefs1[2]));
476: PetscCall(VecDot(deltaX_R, deltaX_Q, &coefs2[0]));
477: PetscCall(VecDot(deltaX_R, deltaX_R, &coefs2[1]));
478: }
480: a0 = normsqX_Q + data->psisq;
481: b0 = 2 * (PetscRealPart(coefs1[0]) + psisqLambdaUpdate);
482: b1 = 2 * PetscRealPart(coefs2[0]);
483: c0 = PetscRealPart(coefs1[2]) + psisqLambdaUpdate * data->lambda_update - stepSize * stepSize;
484: c1 = 2 * PetscRealPart(coefs1[1]);
485: c2 = PetscRealPart(coefs2[1]);
487: as = b1 * b1 - 4 * a0 * c2;
488: bs = 2 * (b1 * b0 - 2 * a0 * c1);
489: cs = b0 * b0 - 4 * a0 * c0;
491: discriminant = cs + bs * deltaS + as * deltaS * deltaS;
493: if (discriminant < 0) {
494: /* Take deltaS < 1 with the unique root -b/(2*a) */
495: PetscReal x1;
497: /* Compute deltaS to be the largest root of (as * x^2 + bs * x + cs = 0) */
498: PetscQuadraticRoots(as, bs, cs, &x1, &deltaS);
499: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", discriminant=%18.16e < 0, shrinking residual update size to deltaS = %18.16e\n", snes->iter, (double)discriminant, (double)deltaS));
500: deltaLambda = -0.5 * (b0 + b1 * deltaS) / a0;
501: } else {
502: /* Use deltaS = 1, pick root that is closest to the last point to prevent doubling back */
503: PetscReal dlambda1, dlambda2;
505: PetscQuadraticRoots(a0, b0 + b1, c0 + c1 + c2, &dlambda1, &dlambda2);
506: deltaLambda = (b0 * dlambda1) > (b0 * dlambda2) ? dlambda1 : dlambda2;
507: }
508: }
509: }
510: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
511: data->lambda = data->lambda + deltaLambda;
512: if (data->lambda > data->lambda_max) {
513: /* Ensure that lambda = lambda_max exactly at the end of incremental process. This ensures the final solution matches the problem we want to solve. */
514: deltaLambda = deltaLambda - (data->lambda - data->lambda_max);
515: data->lambda = data->lambda_max;
516: }
517: if (data->lambda < data->lambda_min) {
518: // LCOV_EXCL_START
519: /* Ensure that lambda >= lambda_min. This prevents some potential oscillatory behavior. */
520: deltaLambda = deltaLambda - (data->lambda - data->lambda_min);
521: data->lambda = data->lambda_min;
522: // LCOV_EXCL_STOP
523: }
524: data->lambda_update = data->lambda_update + deltaLambda;
525: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
526: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", lambda=%18.16e, lambda_update=%18.16e\n", snes->iter, (double)data->lambda, (double)data->lambda_update));
527: if (j == 0) {
528: /* deltaX = deltaLambda*deltaX_Q */
529: PetscCall(VecCopy(deltaX_Q, deltaX));
530: PetscCall(VecScale(deltaX, deltaLambda));
531: } else {
532: /* deltaX = deltaS*deltaX_R + deltaLambda*deltaX_Q */
533: PetscCall(VecAXPBYPCZ(deltaX, deltaS, deltaLambda, 0, deltaX_R, deltaX_Q));
534: }
535: PetscCall(VecAXPY(DeltaX, 1, deltaX));
536: PetscCall(VecAXPY(X, 1, deltaX));
537: /* Q = -dF/dlambda(X, lambda)*/
538: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
539: /* R = F(X, lambda) */
540: PetscCall(SNESComputeFunction(snes, X, R));
541: PetscCall(VecNormBegin(R, NORM_2, &fnorm));
542: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
543: PetscCall(VecNormBegin(deltaX, NORM_2, &ynorm));
544: PetscCall(VecNormEnd(R, NORM_2, &fnorm));
545: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
546: PetscCall(VecNormEnd(deltaX, NORM_2, &ynorm));
547: SNESCheckFunctionDomainError(snes, fnorm);
548: if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));
550: /* Monitor convergence */
551: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
552: snes->iter++;
553: snes->norm = fnorm;
554: snes->ynorm = ynorm;
555: snes->xnorm = xnorm;
556: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
557: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
558: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
559: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
560: if (!snes->reason && j == maxits - 1) snes->reason = SNES_DIVERGED_MAX_IT;
561: if (snes->reason) break;
562: }
563: if (snes->reason < 0) break;
564: if (data->lambda >= data->lambda_max) {
565: break;
566: } else if (maxincs > 0 && i == maxincs - 1) {
567: snes->reason = SNES_DIVERGED_MAX_IT;
568: break;
569: } else {
570: snes->reason = SNES_CONVERGED_ITERATING;
571: }
572: }
573: /* Reset RHS vector, if changed */
574: if (data->copied_rhs) {
575: PetscCall(VecCopy(data->vec_rhs_orig, snes->vec_rhs));
576: data->copied_rhs = PETSC_FALSE;
577: }
578: snes->max_its = maxits; /* reset snes->max_its */
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*
583: SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
584: of the SNESNEWTONAL nonlinear solver.
586: Input Parameter:
587: . snes - the SNES context
588: . x - the solution vector
590: Application Interface Routine: SNESSetUp()
591: */
592: static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
593: {
594: PetscFunctionBegin;
595: PetscCall(SNESSetWorkVecs(snes, 4));
596: PetscCall(SNESSetUpMatrices(snes));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*
601: SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.
603: Input Parameter:
604: . snes - the SNES context
606: Application Interface Routine: SNESSetFromOptions()
607: */
608: static PetscErrorCode SNESSetFromOptions_NEWTONAL(SNES snes, PetscOptionItems PetscOptionsObject)
609: {
610: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
611: SNESNewtonALCorrectionType correction_type = data->correction_type;
613: PetscFunctionBegin;
614: PetscOptionsHeadBegin(PetscOptionsObject, "SNES Newton Arc Length options");
615: PetscCall(PetscOptionsReal("-snes_newtonal_step_size", "Initial arc length increment step size", "SNESNewtonAL", data->step_size, &data->step_size, NULL));
616: PetscCall(PetscOptionsInt("-snes_newtonal_max_continuation_steps", "Maximum number of increment steps", "SNESNewtonAL", data->max_continuation_steps, &data->max_continuation_steps, NULL));
617: PetscCall(PetscOptionsReal("-snes_newtonal_psisq", "Regularization parameter for arc length continuation, 0 for cylindrical", "SNESNewtonAL", data->psisq, &data->psisq, NULL));
618: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_min", "Minimum value of the load parameter lambda", "SNESNewtonAL", data->lambda_min, &data->lambda_min, NULL));
619: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_max", "Maximum value of the load parameter lambda", "SNESNewtonAL", data->lambda_max, &data->lambda_max, NULL));
620: PetscCall(PetscOptionsBool("-snes_newtonal_scale_rhs", "Scale the constant vector passed to `SNESSolve` by the load parameter lambda", "SNESNewtonAL", data->scale_rhs, &data->scale_rhs, NULL));
621: PetscCall(PetscOptionsEnum("-snes_newtonal_correction_type", "Type of correction to use in the arc-length continuation method", "SNESNewtonALCorrectionType", SNESNewtonALCorrectionTypes, (PetscEnum)correction_type, (PetscEnum *)&correction_type, NULL));
622: PetscCall(SNESNewtonALSetCorrectionType(snes, correction_type));
623: PetscOptionsHeadEnd();
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
628: {
629: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
631: PetscFunctionBegin;
632: PetscCall(VecDestroy(&al->vec_rhs_orig));
633: PetscCall(MatDestroy(&al->mat_diag_scaling));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /*
638: SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
639: with SNESCreate_NEWTONAL().
641: Input Parameter:
642: . snes - the SNES context
644: Application Interface Routine: SNESDestroy()
645: */
646: static PetscErrorCode SNESDestroy_NEWTONAL(SNES snes)
647: {
648: PetscFunctionBegin;
649: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", NULL));
650: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", NULL));
651: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", NULL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", NULL));
653: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", NULL));
654: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", NULL));
655: PetscCall(PetscFree(snes->data));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*MC
660: SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.
662: Options Database Keys:
663: + -snes_newtonal_step_size step - Initial arc length increment step size
664: . -snes_newtonal_max_continuation_steps max - Maximum number of continuation steps, or negative for no limit (not recommended)
665: . -snes_newtonal_psisq psisq - Regularization parameter for arc length continuation, 0 for cylindrical. Larger values generally lead to more steps.
666: . -snes_newtonal_lambda_min lambda_min - Minimum value of the load parameter lambda
667: . -snes_newtonal_lambda_max lambda_max - Maximum value of the load parameter lambda
668: . -snes_newtonal_scale_rhs (true|false) - Scale the constant vector passed to `SNESSolve` by the load parameter `lambda`
669: - -snes_newtonal_correction_type (exact|normal) - Type of correction to use in the arc-length continuation method
671: Level: intermediate
673: Note:
674: The exact correction scheme with partial updates is detailed in {cite}`Ritto-CorreaCamotim2008` and the implementation of the
675: normal correction scheme is based on {cite}`LeonPaulinoPereiraMenezesLages_2011`.
677: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
678: M*/
679: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
680: {
681: SNES_NEWTONAL *arclengthParameters;
683: PetscFunctionBegin;
684: snes->ops->setup = SNESSetUp_NEWTONAL;
685: snes->ops->solve = SNESSolve_NEWTONAL;
686: snes->ops->destroy = SNESDestroy_NEWTONAL;
687: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONAL;
688: snes->ops->reset = SNESReset_NEWTONAL;
690: snes->usesksp = PETSC_TRUE;
691: snes->usesnpc = PETSC_FALSE;
693: PetscCall(SNESParametersInitialize(snes));
694: PetscObjectParameterSetDefault(snes, max_funcs, 30000);
695: PetscObjectParameterSetDefault(snes, max_its, 50);
697: snes->alwayscomputesfinalresidual = PETSC_TRUE;
699: PetscCall(PetscNew(&arclengthParameters));
700: arclengthParameters->lambda = 0.0;
701: arclengthParameters->lambda_update = 0.0;
702: arclengthParameters->step_size = 1.0;
703: arclengthParameters->max_continuation_steps = 100;
704: arclengthParameters->psisq = 1.0;
705: arclengthParameters->lambda_min = 0.0;
706: arclengthParameters->lambda_max = 1.0;
707: arclengthParameters->scale_rhs = PETSC_TRUE;
708: arclengthParameters->correction_type = SNES_NEWTONAL_CORRECTION_EXACT;
709: arclengthParameters->mat_diag_scaling = NULL;
710: snes->data = (void *)arclengthParameters;
712: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", SNESNewtonALSetCorrectionType_NEWTONAL));
713: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", SNESNewtonALGetLoadParameter_NEWTONAL));
714: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", SNESNewtonALSetFunction_NEWTONAL));
715: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", SNESNewtonALGetFunction_NEWTONAL));
716: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", SNESNewtonALSetDiagonalScaling_NEWTONAL));
717: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", SNESNewtonALComputeFunction_NEWTONAL));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }