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;
190: PetscTryMethod(snes, "SNESNewtonALSetDiagonalScaling_C", (SNES, Vec), (snes, v));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode SNESNewtonALSetDiagonalScaling_NEWTONAL(SNES snes, Vec v)
195: {
196: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
198: PetscFunctionBegin;
199: if (v) PetscCall(PetscObjectReference((PetscObject)v));
200: PetscCall(MatDestroy(&al->mat_diag_scaling));
201: if (v) {
202: PetscCall(MatCreateDiagonal(v, &al->mat_diag_scaling));
203: PetscCall(PetscObjectDereference((PetscObject)v));
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes, PetscReal *lambda)
209: {
210: SNES_NEWTONAL *al;
212: PetscFunctionBeginHot;
213: al = (SNES_NEWTONAL *)snes->data;
214: *lambda = al->lambda;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@C
219: SNESNewtonALGetLoadParameter - Get the value of the load parameter `lambda` for the arc-length continuation method.
221: Logically Collective
223: Input Parameter:
224: . snes - the nonlinear solver object
226: Output Parameter:
227: . lambda - the arc-length parameter
229: Level: intermediate
231: Notes:
232: This function should be used in the functions provided to `SNESSetFunction()` and `SNESNewtonALSetFunction()`
233: to compute the residual and tangent load vectors for a given value of `lambda` (0 <= lambda <= 1).
235: Usually, `lambda` is used to scale the external force vector in the residual function, i.e. proportional loading,
236: in which case the tangent load vector is the full external force vector.
238: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
239: @*/
240: PetscErrorCode SNESNewtonALGetLoadParameter(SNES snes, PetscReal *lambda)
241: {
242: PetscFunctionBeginHot;
244: PetscAssertPointer(lambda, 2);
245: PetscUseMethod(snes, "SNESNewtonALGetLoadParameter_C", (SNES, PetscReal *), (snes, lambda));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode SNESNewtonALComputeFunction_NEWTONAL(SNES snes, Vec X, Vec Q)
250: {
251: void *ctx = NULL;
252: SNESFunctionFn *computealfunction = NULL;
253: SNES_NEWTONAL *al;
255: PetscFunctionBegin;
256: al = (SNES_NEWTONAL *)snes->data;
257: PetscCall(SNESNewtonALGetFunction(snes, &computealfunction, &ctx));
259: PetscCall(VecZeroEntries(Q));
260: 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");
261: if (computealfunction) {
262: PetscCall(VecLockReadPush(X));
263: PetscCallBack("SNES callback NewtonAL tangent load function", (*computealfunction)(snes, X, Q, ctx));
264: PetscCall(VecLockReadPop(X));
265: }
266: if (al->scale_rhs && snes->vec_rhs) {
267: /* Save original RHS vector values, then scale `snes->vec_rhs` by load parameter */
268: if (!al->vec_rhs_orig) PetscCall(VecDuplicate(snes->vec_rhs, &al->vec_rhs_orig));
269: if (!al->copied_rhs) {
270: PetscCall(VecCopy(snes->vec_rhs, al->vec_rhs_orig));
271: al->copied_rhs = PETSC_TRUE;
272: }
273: PetscCall(VecAXPBY(snes->vec_rhs, al->lambda, 0.0, al->vec_rhs_orig));
274: PetscCall(VecAXPY(Q, 1, al->vec_rhs_orig));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@C
280: SNESNewtonALComputeFunction - Calls the function that has been set with `SNESNewtonALSetFunction()`.
282: Collective
284: Input Parameters:
285: + snes - the `SNES` context
286: - X - input vector
288: Output Parameter:
289: . Q - tangent load vector, as set by `SNESNewtonALSetFunction()`
291: Level: developer
293: Note:
294: `SNESNewtonALComputeFunction()` is typically used within nonlinear solvers
295: implementations, so users would not generally call this routine themselves.
297: .seealso: [](ch_snes), `SNES`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`
298: @*/
299: PetscErrorCode SNESNewtonALComputeFunction(SNES snes, Vec X, Vec Q)
300: {
301: PetscFunctionBegin;
305: PetscCheckSameComm(snes, 1, X, 2);
306: PetscCheckSameComm(snes, 1, Q, 3);
307: PetscCall(VecValidValues_Internal(X, 2, PETSC_TRUE));
308: PetscCall(PetscLogEventBegin(SNES_NewtonALEval, snes, X, Q, 0));
309: PetscTryMethod(snes, "SNESNewtonALComputeFunction_C", (SNES, Vec, Vec), (snes, X, Q));
310: PetscCall(PetscLogEventEnd(SNES_NewtonALEval, snes, X, Q, 0));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*
315: SNESSolve_NEWTONAL - Solves a nonlinear system with Newton's method with arc length continuation.
317: Input Parameter:
318: . snes - the `SNES` context
320: Application Interface Routine: SNESSolve()
321: */
322: static PetscErrorCode SNESSolve_NEWTONAL(SNES snes)
323: {
324: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
325: PetscInt maxits, maxincs, lits;
326: PetscReal fnorm, xnorm, ynorm, stepSize;
327: Vec DeltaX, deltaX, X, R, Q, deltaX_Q, deltaX_R;
328: Mat W;
330: PetscFunctionBegin;
331: 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);
333: /* Register citations */
334: PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
335: if (data->correction_type == SNES_NEWTONAL_CORRECTION_EXACT) {
336: PetscCall(PetscCitationsRegister(NewtonALExactCitation, &NewtonALExactCitationSet));
337: } else if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
338: PetscCall(PetscCitationsRegister(NewtonALNormalCitation, &NewtonALNormalCitationSet));
339: }
341: data->copied_rhs = PETSC_FALSE;
342: data->lambda_update = 0.0;
343: data->lambda = 0.0;
344: snes->numFailures = 0;
345: snes->numLinearSolveFailures = 0;
346: snes->reason = SNES_CONVERGED_ITERATING;
347: snes->iter = 0;
349: maxits = snes->max_its; /* maximum number of iterations */
350: maxincs = data->max_continuation_steps; /* maximum number of increments */
351: X = snes->vec_sol; /* solution vector */
352: R = snes->vec_func; /* residual vector */
353: Q = snes->work[0]; /* tangent load vector */
354: W = data->mat_diag_scaling; /* diagonal scaling for DoFs */
355: deltaX_Q = snes->work[1]; /* variation of X with respect to lambda */
356: deltaX_R = snes->work[2]; /* linearized error correction */
357: DeltaX = snes->work[3]; /* step from equilibrium */
358: deltaX = snes->vec_sol_update; /* full newton step */
359: stepSize = data->step_size; /* initial step size */
361: PetscCall(VecZeroEntries(DeltaX));
363: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
364: /* set snes->max_its for convergence test */
365: snes->max_its = maxits * maxincs;
366: snes->iter = 0;
367: snes->norm = 0.0;
368: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
370: /* main incremental-iterative loop */
371: for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
372: PetscReal deltaLambda;
374: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
375: snes->norm = 0.0;
376: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
377: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
378: PetscCall(SNESComputeFunction(snes, X, R));
379: PetscCall(VecAXPY(R, 1, Q)); /* R <- R + Q */
380: PetscCall(VecNorm(R, NORM_2, &fnorm)); /* fnorm <- ||R|| */
381: SNESCheckFunctionDomainError(snes, fnorm);
383: /* Monitor convergence */
384: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
385: PetscCall(SNESMonitor(snes, snes->iter, fnorm));
386: if (snes->reason) break;
387: for (PetscInt j = 0; j < maxits; j++) {
388: PetscReal normsqX_Q, deltaS = 1;
390: /* Call general purpose update function */
391: PetscTryTypeMethod(snes, update, snes->iter);
393: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
394: SNESCheckJacobianDomainError(snes);
395: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
396: /* Solve J deltaX_Q = Q, where J is Jacobian matrix */
397: PetscCall(KSPSolve(snes->ksp, Q, deltaX_Q));
398: SNESCheckKSPSolve(snes);
399: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
400: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", tangent load linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
401: /* Compute load parameter variation */
402: if (W) PetscCall(MatANorm(W, deltaX_Q, &normsqX_Q));
403: else PetscCall(VecNorm(deltaX_Q, NORM_2, &normsqX_Q));
404: normsqX_Q *= normsqX_Q;
405: /* On first iter, use predictor. This is the same regardless of corrector scheme. */
406: if (j == 0) {
407: PetscReal sign = 1.0;
408: if (i > 0) {
409: PetscScalar dot;
410: if (W) PetscCall(MatADot(W, DeltaX, deltaX_Q, &dot));
411: else PetscCall(VecDot(DeltaX, deltaX_Q, &dot));
412: sign = PetscRealPart(dot) + data->psisq * data->lambda_update;
413: sign = sign >= 0 ? 1.0 : -1.0;
414: }
415: data->lambda_update = 0.0;
416: PetscCall(VecZeroEntries(DeltaX));
417: deltaLambda = sign * stepSize / PetscSqrtReal(normsqX_Q + data->psisq);
418: } else {
419: /* Solve J deltaX_R = -R */
420: PetscCall(KSPSolve(snes->ksp, R, deltaX_R));
421: SNESCheckKSPSolve(snes);
422: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
423: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", residual linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
424: PetscCall(VecScale(deltaX_R, -1));
426: if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
427: /*
428: Take a step orthogonal to the current incremental update DeltaX.
429: Note, this approach is cheaper than the exact correction, but may exhibit convergence
430: issues due to the iterative trial points not being on the quadratic constraint surface.
431: On the bright side, we always have a real and unique solution for deltaLambda.
432: */
433: PetscScalar coefs[2];
434: if (W) {
435: PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs[0]));
436: PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs[1]));
437: } else {
438: PetscCall(VecDot(DeltaX, deltaX_R, &coefs[0]));
439: PetscCall(VecDot(DeltaX, deltaX_Q, &coefs[1]));
440: }
441: deltaLambda = -PetscRealPart(coefs[0]) / (PetscRealPart(coefs[1]) + data->psisq * data->lambda_update);
442: } else {
443: /*
444: Solve
445: a*deltaLambda^2 + b*deltaLambda + c = 0 (*)
446: where
447: a = a0
448: b = b0 + b1*deltaS
449: c = c0 + c1*deltaS + c2*deltaS^2
450: and deltaS is either 1, or the largest value in (0, 1) that satisfies
451: b^2 - 4*a*c = as*deltaS^2 + bs*deltaS + cs >= 0
452: where
453: as = b1^2 - 4*a0*c2
454: bs = 2*b1*b0 - 4*a0*c1
455: cs = b0^2 - 4*a0*c0
456: These "partial corrections" prevent (*) from having complex roots.
457: */
458: PetscReal psisqLambdaUpdate, discriminant;
459: PetscReal as, bs, cs;
460: PetscReal a0, b0, b1, c0, c1, c2;
461: PetscScalar coefs1[3]; /* coefs[0] = deltaX_Q*DeltaX, coefs[1] = deltaX_R*DeltaX, coefs[2] = DeltaX*DeltaX */
462: PetscScalar coefs2[2]; /* coefs[0] = deltaX_Q*deltaX_R, coefs[1] = deltaX_R*deltaX_R */
464: psisqLambdaUpdate = data->psisq * data->lambda_update;
465: if (W) {
466: PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs1[0]));
467: PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs1[1]));
468: PetscCall(MatADot(W, DeltaX, DeltaX, &coefs1[2]));
469: PetscCall(MatADot(W, deltaX_R, deltaX_Q, &coefs2[0]));
470: PetscCall(MatADot(W, deltaX_R, deltaX_R, &coefs2[1]));
471: } else {
472: PetscCall(VecDot(DeltaX, deltaX_Q, &coefs1[0]));
473: PetscCall(VecDot(DeltaX, deltaX_R, &coefs1[1]));
474: PetscCall(VecDot(DeltaX, DeltaX, &coefs1[2]));
475: PetscCall(VecDot(deltaX_R, deltaX_Q, &coefs2[0]));
476: PetscCall(VecDot(deltaX_R, deltaX_R, &coefs2[1]));
477: }
479: a0 = normsqX_Q + data->psisq;
480: b0 = 2 * (PetscRealPart(coefs1[0]) + psisqLambdaUpdate);
481: b1 = 2 * PetscRealPart(coefs2[0]);
482: c0 = PetscRealPart(coefs1[2]) + psisqLambdaUpdate * data->lambda_update - stepSize * stepSize;
483: c1 = 2 * PetscRealPart(coefs1[1]);
484: c2 = PetscRealPart(coefs2[1]);
486: as = b1 * b1 - 4 * a0 * c2;
487: bs = 2 * (b1 * b0 - 2 * a0 * c1);
488: cs = b0 * b0 - 4 * a0 * c0;
490: discriminant = cs + bs * deltaS + as * deltaS * deltaS;
492: if (discriminant < 0) {
493: /* Take deltaS < 1 with the unique root -b/(2*a) */
494: PetscReal x1;
496: /* Compute deltaS to be the largest root of (as * x^2 + bs * x + cs = 0) */
497: PetscQuadraticRoots(as, bs, cs, &x1, &deltaS);
498: 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));
499: deltaLambda = -0.5 * (b0 + b1 * deltaS) / a0;
500: } else {
501: /* Use deltaS = 1, pick root that is closest to the last point to prevent doubling back */
502: PetscReal dlambda1, dlambda2;
504: PetscQuadraticRoots(a0, b0 + b1, c0 + c1 + c2, &dlambda1, &dlambda2);
505: deltaLambda = (b0 * dlambda1) > (b0 * dlambda2) ? dlambda1 : dlambda2;
506: }
507: }
508: }
509: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
510: data->lambda = data->lambda + deltaLambda;
511: if (data->lambda > data->lambda_max) {
512: /* Ensure that lambda = lambda_max exactly at the end of incremental process. This ensures the final solution matches the problem we want to solve. */
513: deltaLambda = deltaLambda - (data->lambda - data->lambda_max);
514: data->lambda = data->lambda_max;
515: }
516: if (data->lambda < data->lambda_min) {
517: // LCOV_EXCL_START
518: /* Ensure that lambda >= lambda_min. This prevents some potential oscillatory behavior. */
519: deltaLambda = deltaLambda - (data->lambda - data->lambda_min);
520: data->lambda = data->lambda_min;
521: // LCOV_EXCL_STOP
522: }
523: data->lambda_update = data->lambda_update + deltaLambda;
524: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
525: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", lambda=%18.16e, lambda_update=%18.16e\n", snes->iter, (double)data->lambda, (double)data->lambda_update));
526: if (j == 0) {
527: /* deltaX = deltaLambda*deltaX_Q */
528: PetscCall(VecCopy(deltaX_Q, deltaX));
529: PetscCall(VecScale(deltaX, deltaLambda));
530: } else {
531: /* deltaX = deltaS*deltaX_R + deltaLambda*deltaX_Q */
532: PetscCall(VecAXPBYPCZ(deltaX, deltaS, deltaLambda, 0, deltaX_R, deltaX_Q));
533: }
534: PetscCall(VecAXPY(DeltaX, 1, deltaX));
535: PetscCall(VecAXPY(X, 1, deltaX));
536: /* Q = -dF/dlambda(X, lambda)*/
537: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
538: /* R = F(X, lambda) */
539: PetscCall(SNESComputeFunction(snes, X, R));
540: PetscCall(VecNormBegin(R, NORM_2, &fnorm));
541: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
542: PetscCall(VecNormBegin(deltaX, NORM_2, &ynorm));
543: PetscCall(VecNormEnd(R, NORM_2, &fnorm));
544: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
545: PetscCall(VecNormEnd(deltaX, NORM_2, &ynorm));
546: SNESCheckFunctionDomainError(snes, fnorm);
547: if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));
549: /* Monitor convergence */
550: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
551: snes->iter++;
552: snes->norm = fnorm;
553: snes->ynorm = ynorm;
554: snes->xnorm = xnorm;
555: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
556: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
557: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
558: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
559: if (!snes->reason && j == maxits - 1) snes->reason = SNES_DIVERGED_MAX_IT;
560: if (snes->reason) break;
561: }
562: if (snes->reason < 0) break;
563: if (data->lambda >= data->lambda_max) {
564: break;
565: } else if (maxincs > 0 && i == maxincs - 1) {
566: snes->reason = SNES_DIVERGED_MAX_IT;
567: break;
568: } else {
569: snes->reason = SNES_CONVERGED_ITERATING;
570: }
571: }
572: /* Reset RHS vector, if changed */
573: if (data->copied_rhs) {
574: PetscCall(VecCopy(data->vec_rhs_orig, snes->vec_rhs));
575: data->copied_rhs = PETSC_FALSE;
576: }
577: snes->max_its = maxits; /* reset snes->max_its */
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*
582: SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
583: of the SNESNEWTONAL nonlinear solver.
585: Input Parameter:
586: . snes - the SNES context
587: . x - the solution vector
589: Application Interface Routine: SNESSetUp()
590: */
591: static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
592: {
593: PetscFunctionBegin;
594: PetscCall(SNESSetWorkVecs(snes, 4));
595: PetscCall(SNESSetUpMatrices(snes));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*
600: SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.
602: Input Parameter:
603: . snes - the SNES context
605: Application Interface Routine: SNESSetFromOptions()
606: */
607: static PetscErrorCode SNESSetFromOptions_NEWTONAL(SNES snes, PetscOptionItems PetscOptionsObject)
608: {
609: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
610: SNESNewtonALCorrectionType correction_type = data->correction_type;
612: PetscFunctionBegin;
613: PetscOptionsHeadBegin(PetscOptionsObject, "SNES Newton Arc Length options");
614: PetscCall(PetscOptionsReal("-snes_newtonal_step_size", "Initial arc length increment step size", "SNESNewtonAL", data->step_size, &data->step_size, NULL));
615: PetscCall(PetscOptionsInt("-snes_newtonal_max_continuation_steps", "Maximum number of increment steps", "SNESNewtonAL", data->max_continuation_steps, &data->max_continuation_steps, NULL));
616: PetscCall(PetscOptionsReal("-snes_newtonal_psisq", "Regularization parameter for arc length continuation, 0 for cylindrical", "SNESNewtonAL", data->psisq, &data->psisq, NULL));
617: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_min", "Minimum value of the load parameter lambda", "SNESNewtonAL", data->lambda_min, &data->lambda_min, NULL));
618: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_max", "Maximum value of the load parameter lambda", "SNESNewtonAL", data->lambda_max, &data->lambda_max, NULL));
619: 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));
620: 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));
621: PetscCall(SNESNewtonALSetCorrectionType(snes, correction_type));
622: PetscOptionsHeadEnd();
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
627: {
628: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
630: PetscFunctionBegin;
631: PetscCall(VecDestroy(&al->vec_rhs_orig));
632: PetscCall(MatDestroy(&al->mat_diag_scaling));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*
637: SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
638: with SNESCreate_NEWTONAL().
640: Input Parameter:
641: . snes - the SNES context
643: Application Interface Routine: SNESDestroy()
644: */
645: static PetscErrorCode SNESDestroy_NEWTONAL(SNES snes)
646: {
647: PetscFunctionBegin;
648: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", NULL));
649: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", NULL));
650: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", NULL));
651: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", NULL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", NULL));
653: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", NULL));
654: PetscCall(PetscFree(snes->data));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*MC
659: SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.
661: Options Database Keys:
662: + -snes_newtonal_step_size step - Initial arc length increment step size
663: . -snes_newtonal_max_continuation_steps max - Maximum number of continuation steps, or negative for no limit (not recommended)
664: . -snes_newtonal_psisq psisq - Regularization parameter for arc length continuation, 0 for cylindrical. Larger values generally lead to more steps.
665: . -snes_newtonal_lambda_min lambda_min - Minimum value of the load parameter lambda
666: . -snes_newtonal_lambda_max lambda_max - Maximum value of the load parameter lambda
667: . -snes_newtonal_scale_rhs (true|false) - Scale the constant vector passed to `SNESSolve` by the load parameter `lambda`
668: - -snes_newtonal_correction_type (exact|normal) - Type of correction to use in the arc-length continuation method
670: Level: intermediate
672: Note:
673: The exact correction scheme with partial updates is detailed in {cite}`Ritto-CorreaCamotim2008` and the implementation of the
674: normal correction scheme is based on {cite}`LeonPaulinoPereiraMenezesLages_2011`.
676: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
677: M*/
678: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
679: {
680: SNES_NEWTONAL *arclengthParameters;
682: PetscFunctionBegin;
683: snes->ops->setup = SNESSetUp_NEWTONAL;
684: snes->ops->solve = SNESSolve_NEWTONAL;
685: snes->ops->destroy = SNESDestroy_NEWTONAL;
686: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONAL;
687: snes->ops->reset = SNESReset_NEWTONAL;
689: snes->usesksp = PETSC_TRUE;
690: snes->usesnpc = PETSC_FALSE;
692: PetscCall(SNESParametersInitialize(snes));
693: PetscObjectParameterSetDefault(snes, max_funcs, 30000);
694: PetscObjectParameterSetDefault(snes, max_its, 50);
696: snes->alwayscomputesfinalresidual = PETSC_TRUE;
698: PetscCall(PetscNew(&arclengthParameters));
699: arclengthParameters->lambda = 0.0;
700: arclengthParameters->lambda_update = 0.0;
701: arclengthParameters->step_size = 1.0;
702: arclengthParameters->max_continuation_steps = 100;
703: arclengthParameters->psisq = 1.0;
704: arclengthParameters->lambda_min = 0.0;
705: arclengthParameters->lambda_max = 1.0;
706: arclengthParameters->scale_rhs = PETSC_TRUE;
707: arclengthParameters->correction_type = SNES_NEWTONAL_CORRECTION_EXACT;
708: arclengthParameters->mat_diag_scaling = NULL;
709: snes->data = (void *)arclengthParameters;
711: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", SNESNewtonALSetCorrectionType_NEWTONAL));
712: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", SNESNewtonALGetLoadParameter_NEWTONAL));
713: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", SNESNewtonALSetFunction_NEWTONAL));
714: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", SNESNewtonALGetFunction_NEWTONAL));
715: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", SNESNewtonALSetDiagonalScaling_NEWTONAL));
716: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", SNESNewtonALComputeFunction_NEWTONAL));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }