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: PetscCall(VecNorm(XStep, NORM_2, &arcLength));
47: arcLength = PetscSqrtReal(PetscSqr(arcLength) + al->psisq * lambdaStep * lambdaStep);
48: arcLengthError = PetscAbsReal(arcLength - stepSize);
50: 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));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
55: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
56: {
57: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
58: PetscReal x1 = temp / a;
59: PetscReal x2 = c / temp;
60: *xm = PetscMin(x1, x2);
61: *xp = PetscMax(x1, x2);
62: }
64: static PetscErrorCode SNESNewtonALSetCorrectionType_NEWTONAL(SNES snes, SNESNewtonALCorrectionType ctype)
65: {
66: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
68: PetscFunctionBegin;
69: al->correction_type = ctype;
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@
74: SNESNewtonALSetCorrectionType - Set the type of correction to use in the arc-length continuation method.
76: Logically Collective
78: Input Parameters:
79: + snes - the nonlinear solver object
80: - ctype - the type of correction to use
82: Options Database Key:
83: . -snes_newtonal_correction_type <type> - Set the type of correction to use; use -help for a list of available types
85: Level: intermediate
87: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALCorrectionType`
88: @*/
89: PetscErrorCode SNESNewtonALSetCorrectionType(SNES snes, SNESNewtonALCorrectionType ctype)
90: {
91: PetscFunctionBegin;
94: PetscTryMethod(snes, "SNESNewtonALSetCorrectionType_C", (SNES, SNESNewtonALCorrectionType), (snes, ctype));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode SNESNewtonALSetFunction_NEWTONAL(SNES snes, SNESFunctionFn *func, void *ctx)
99: {
100: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
102: PetscFunctionBegin;
103: al->computealfunction = func;
104: al->alctx = ctx;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: /*@C
109: SNESNewtonALSetFunction - Sets a user function that is called at each function evaluation to
110: compute the tangent load vector for the arc-length continuation method.
112: Logically Collective
114: Input Parameters:
115: + snes - the nonlinear solver object
116: . 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
117: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
119: Level: intermediate
121: Notes:
122: If the current value of the load parameter is needed in `func`, it can be obtained with `SNESNewtonALGetLoadParameter()`.
124: The tangent load vector is the partial derivative of external load with respect to the load parameter.
125: In the case of proportional loading, the tangent load vector is the full external load vector at the end of the load step.
127: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
128: @*/
129: PetscErrorCode SNESNewtonALSetFunction(SNES snes, SNESFunctionFn *func, void *ctx)
130: {
131: PetscFunctionBegin;
133: PetscTryMethod(snes, "SNESNewtonALSetFunction_C", (SNES, SNESFunctionFn *, void *), (snes, func, ctx));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode SNESNewtonALGetFunction_NEWTONAL(SNES snes, SNESFunctionFn **func, void **ctx)
138: {
139: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
141: PetscFunctionBegin;
142: if (func) *func = al->computealfunction;
143: if (ctx) *ctx = al->alctx;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@C
148: SNESNewtonALGetFunction - Get the user function and context set with `SNESNewtonALSetFunction`
150: Logically Collective
152: Input Parameters:
153: + snes - the nonlinear solver object
154: . func - [optional] tangent load function evaluation routine, see `SNESNewtonALSetFunction()` for the call sequence
155: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
157: Level: intermediate
159: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
160: @*/
161: PetscErrorCode SNESNewtonALGetFunction(SNES snes, SNESFunctionFn **func, void **ctx)
162: {
163: PetscFunctionBegin;
165: PetscUseMethod(snes, "SNESNewtonALGetFunction_C", (SNES, SNESFunctionFn **, void **), (snes, func, ctx));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes, PetscReal *lambda)
170: {
171: SNES_NEWTONAL *al;
173: PetscFunctionBeginHot;
174: al = (SNES_NEWTONAL *)snes->data;
175: *lambda = al->lambda;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@C
180: SNESNewtonALGetLoadParameter - Get the value of the load parameter `lambda` for the arc-length continuation method.
182: Logically Collective
184: Input Parameter:
185: . snes - the nonlinear solver object
187: Output Parameter:
188: . lambda - the arc-length parameter
190: Level: intermediate
192: Notes:
193: This function should be used in the functions provided to `SNESSetFunction()` and `SNESNewtonALSetFunction()`
194: to compute the residual and tangent load vectors for a given value of `lambda` (0 <= lambda <= 1).
196: Usually, `lambda` is used to scale the external force vector in the residual function, i.e. proportional loading,
197: in which case the tangent load vector is the full external force vector.
199: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
200: @*/
201: PetscErrorCode SNESNewtonALGetLoadParameter(SNES snes, PetscReal *lambda)
202: {
203: PetscFunctionBeginHot;
205: PetscAssertPointer(lambda, 2);
206: PetscUseMethod(snes, "SNESNewtonALGetLoadParameter_C", (SNES, PetscReal *), (snes, lambda));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static PetscErrorCode SNESNewtonALComputeFunction_NEWTONAL(SNES snes, Vec X, Vec Q)
211: {
212: void *ctx = NULL;
213: SNESFunctionFn *computealfunction = NULL;
214: SNES_NEWTONAL *al;
216: PetscFunctionBegin;
217: al = (SNES_NEWTONAL *)snes->data;
218: PetscCall(SNESNewtonALGetFunction(snes, &computealfunction, &ctx));
220: PetscCall(VecZeroEntries(Q));
221: 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");
222: if (computealfunction) {
223: PetscCall(VecLockReadPush(X));
224: PetscCallBack("SNES callback NewtonAL tangent load function", (*computealfunction)(snes, X, Q, ctx));
225: PetscCall(VecLockReadPop(X));
226: }
227: if (al->scale_rhs && snes->vec_rhs) {
228: /* Save original RHS vector values, then scale `snes->vec_rhs` by load parameter */
229: if (!al->vec_rhs_orig) PetscCall(VecDuplicate(snes->vec_rhs, &al->vec_rhs_orig));
230: if (!al->copied_rhs) {
231: PetscCall(VecCopy(snes->vec_rhs, al->vec_rhs_orig));
232: al->copied_rhs = PETSC_TRUE;
233: }
234: PetscCall(VecAXPBY(snes->vec_rhs, al->lambda, 0.0, al->vec_rhs_orig));
235: PetscCall(VecAXPY(Q, 1, al->vec_rhs_orig));
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@C
241: SNESNewtonALComputeFunction - Calls the function that has been set with `SNESNewtonALSetFunction()`.
243: Collective
245: Input Parameters:
246: + snes - the `SNES` context
247: - X - input vector
249: Output Parameter:
250: . Q - tangent load vector, as set by `SNESNewtonALSetFunction()`
252: Level: developer
254: Note:
255: `SNESNewtonALComputeFunction()` is typically used within nonlinear solvers
256: implementations, so users would not generally call this routine themselves.
258: .seealso: [](ch_snes), `SNES`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`
259: @*/
260: PetscErrorCode SNESNewtonALComputeFunction(SNES snes, Vec X, Vec Q)
261: {
262: PetscFunctionBegin;
266: PetscCheckSameComm(snes, 1, X, 2);
267: PetscCheckSameComm(snes, 1, Q, 3);
268: PetscCall(VecValidValues_Internal(X, 2, PETSC_TRUE));
269: PetscCall(PetscLogEventBegin(SNES_NewtonALEval, snes, X, Q, 0));
270: PetscTryMethod(snes, "SNESNewtonALComputeFunction_C", (SNES, Vec, Vec), (snes, X, Q));
271: PetscCall(PetscLogEventEnd(SNES_NewtonALEval, snes, X, Q, 0));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*
276: SNESSolve_NEWTONAL - Solves a nonlinear system with Newton's method with arc length continuation.
278: Input Parameter:
279: . snes - the `SNES` context
281: Application Interface Routine: SNESSolve()
282: */
283: static PetscErrorCode SNESSolve_NEWTONAL(SNES snes)
284: {
285: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
286: PetscInt maxits, maxincs, lits;
287: PetscReal fnorm, xnorm, ynorm, stepSize;
288: Vec DeltaX, deltaX, X, R, Q, deltaX_Q, deltaX_R;
290: PetscFunctionBegin;
291: 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);
293: /* Register citations */
294: PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
295: if (data->correction_type == SNES_NEWTONAL_CORRECTION_EXACT) {
296: PetscCall(PetscCitationsRegister(NewtonALExactCitation, &NewtonALExactCitationSet));
297: } else if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
298: PetscCall(PetscCitationsRegister(NewtonALNormalCitation, &NewtonALNormalCitationSet));
299: }
301: data->copied_rhs = PETSC_FALSE;
302: data->lambda_update = 0.0;
303: data->lambda = 0.0;
304: snes->numFailures = 0;
305: snes->numLinearSolveFailures = 0;
306: snes->reason = SNES_CONVERGED_ITERATING;
307: snes->iter = 0;
309: maxits = snes->max_its; /* maximum number of iterations */
310: maxincs = data->max_continuation_steps; /* maximum number of increments */
311: X = snes->vec_sol; /* solution vector */
312: R = snes->vec_func; /* residual vector */
313: Q = snes->work[0]; /* tangent load vector */
314: deltaX_Q = snes->work[1]; /* variation of X with respect to lambda */
315: deltaX_R = snes->work[2]; /* linearized error correction */
316: DeltaX = snes->work[3]; /* step from equilibrium */
317: deltaX = snes->vec_sol_update; /* full newton step */
318: stepSize = data->step_size; /* initial step size */
320: PetscCall(VecZeroEntries(DeltaX));
322: /* set snes->max_its for convergence test */
323: snes->max_its = maxits * maxincs;
325: /* main incremental-iterative loop */
326: for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
327: PetscReal deltaLambda;
329: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
330: snes->norm = 0.0;
331: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
332: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
333: PetscCall(SNESComputeFunction(snes, X, R));
334: PetscCall(VecAXPY(R, 1, Q)); /* R <- R + Q */
335: PetscCall(VecNorm(R, NORM_2, &fnorm)); /* fnorm <- ||R|| */
336: SNESCheckFunctionNorm(snes, fnorm);
338: /* Monitor convergence */
339: PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
340: PetscCall(SNESMonitor(snes, snes->iter, fnorm));
341: if (i == 0 && snes->reason) break;
342: for (PetscInt j = 0; j < maxits; j++) {
343: PetscReal normsqX_Q, deltaS = 1;
345: /* Call general purpose update function */
346: PetscTryTypeMethod(snes, update, snes->iter);
348: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
349: SNESCheckJacobianDomainerror(snes);
350: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
351: /* Solve J deltaX_Q = Q, where J is Jacobian matrix */
352: PetscCall(KSPSolve(snes->ksp, Q, deltaX_Q));
353: SNESCheckKSPSolve(snes);
354: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
355: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", tangent load linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
356: /* Compute load parameter variation */
357: PetscCall(VecNorm(deltaX_Q, NORM_2, &normsqX_Q));
358: normsqX_Q *= normsqX_Q;
359: /* On first iter, use predictor. This is the same regardless of corrector scheme. */
360: if (j == 0) {
361: PetscReal sign = 1.0;
362: if (i > 0) {
363: PetscCall(VecDotRealPart(DeltaX, deltaX_Q, &sign));
364: sign += data->psisq * data->lambda_update;
365: sign = sign >= 0 ? 1.0 : -1.0;
366: }
367: data->lambda_update = 0.0;
368: PetscCall(VecZeroEntries(DeltaX));
369: deltaLambda = sign * stepSize / PetscSqrtReal(normsqX_Q + data->psisq);
370: } else {
371: /* Solve J deltaX_R = -R */
372: PetscCall(KSPSolve(snes->ksp, R, deltaX_R));
373: SNESCheckKSPSolve(snes);
374: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
375: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", residual linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
376: PetscCall(VecScale(deltaX_R, -1));
378: if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
379: /*
380: Take a step orthogonal to the current incremental update DeltaX.
381: Note, this approach is cheaper than the exact correction, but may exhibit convergence
382: issues due to the iterative trial points not being on the quadratic constraint surface.
383: On the bright side, we always have a real and unique solution for deltaLambda.
384: */
385: PetscScalar coefs[2];
386: Vec rhs[] = {deltaX_R, deltaX_Q};
388: PetscCall(VecMDot(DeltaX, 2, rhs, coefs));
389: deltaLambda = -PetscRealPart(coefs[0]) / (PetscRealPart(coefs[1]) + data->psisq * data->lambda_update);
390: } else {
391: /*
392: Solve
393: a*deltaLambda^2 + b*deltaLambda + c = 0 (*)
394: where
395: a = a0
396: b = b0 + b1*deltaS
397: c = c0 + c1*deltaS + c2*deltaS^2
398: and deltaS is either 1, or the largest value in (0, 1) that satisfies
399: b^2 - 4*a*c = as*deltaS^2 + bs*deltaS + cs >= 0
400: where
401: as = b1^2 - 4*a0*c2
402: bs = 2*b1*b0 - 4*a0*c1
403: cs = b0^2 - 4*a0*c0
404: These "partial corrections" prevent (*) from having complex roots.
405: */
406: PetscReal psisqLambdaUpdate, discriminant;
407: PetscReal as, bs, cs;
408: PetscReal a0, b0, b1, c0, c1, c2;
409: PetscScalar coefs1[3]; /* coefs[0] = deltaX_Q*DeltaX, coefs[1] = deltaX_R*DeltaX, coefs[2] = DeltaX*DeltaX */
410: PetscScalar coefs2[2]; /* coefs[0] = deltaX_Q*deltaX_R, coefs[1] = deltaX_R*deltaX_R */
411: const Vec rhs1[3] = {deltaX_Q, deltaX_R, DeltaX};
412: const Vec rhs2[2] = {deltaX_Q, deltaX_R};
414: psisqLambdaUpdate = data->psisq * data->lambda_update;
415: PetscCall(VecMDotBegin(DeltaX, 3, rhs1, coefs1));
416: PetscCall(VecMDotBegin(deltaX_R, 2, rhs2, coefs2));
417: PetscCall(VecMDotEnd(DeltaX, 3, rhs1, coefs1));
418: PetscCall(VecMDotEnd(deltaX_R, 2, rhs2, coefs2));
420: a0 = normsqX_Q + data->psisq;
421: b0 = 2 * (PetscRealPart(coefs1[0]) + psisqLambdaUpdate);
422: b1 = 2 * PetscRealPart(coefs2[0]);
423: c0 = PetscRealPart(coefs1[2]) + psisqLambdaUpdate * data->lambda_update - stepSize * stepSize;
424: c1 = 2 * PetscRealPart(coefs1[1]);
425: c2 = PetscRealPart(coefs2[1]);
427: as = b1 * b1 - 4 * a0 * c2;
428: bs = 2 * (b1 * b0 - 2 * a0 * c1);
429: cs = b0 * b0 - 4 * a0 * c0;
431: discriminant = cs + bs * deltaS + as * deltaS * deltaS;
433: if (discriminant < 0) {
434: /* Take deltaS < 1 with the unique root -b/(2*a) */
435: PetscReal x1;
437: /* Compute deltaS to be the largest root of (as * x^2 + bs * x + cs = 0) */
438: PetscQuadraticRoots(as, bs, cs, &x1, &deltaS);
439: 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));
440: deltaLambda = -0.5 * (b0 + b1 * deltaS) / a0;
441: } else {
442: /* Use deltaS = 1, pick root that is closest to the last point to prevent doubling back */
443: PetscReal dlambda1, dlambda2;
445: PetscQuadraticRoots(a0, b0 + b1, c0 + c1 + c2, &dlambda1, &dlambda2);
446: deltaLambda = (b0 * dlambda1) > (b0 * dlambda2) ? dlambda1 : dlambda2;
447: }
448: }
449: }
450: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
451: data->lambda = data->lambda + deltaLambda;
452: if (data->lambda > data->lambda_max) {
453: /* Ensure that lambda = lambda_max exactly at the end of incremental process. This ensures the final solution matches the problem we want to solve. */
454: deltaLambda = deltaLambda - (data->lambda - data->lambda_max);
455: data->lambda = data->lambda_max;
456: }
457: if (data->lambda < data->lambda_min) {
458: // LCOV_EXCL_START
459: /* Ensure that lambda >= lambda_min. This prevents some potential oscillatory behavior. */
460: deltaLambda = deltaLambda - (data->lambda - data->lambda_min);
461: data->lambda = data->lambda_min;
462: // LCOV_EXCL_STOP
463: }
464: data->lambda_update = data->lambda_update + deltaLambda;
465: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
466: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", lambda=%18.16e, lambda_update=%18.16e\n", snes->iter, (double)data->lambda, (double)data->lambda_update));
467: if (j == 0) {
468: /* deltaX = deltaLambda*deltaX_Q */
469: PetscCall(VecCopy(deltaX_Q, deltaX));
470: PetscCall(VecScale(deltaX, deltaLambda));
471: } else {
472: /* deltaX = deltaS*deltaX_R + deltaLambda*deltaX_Q */
473: PetscCall(VecAXPBYPCZ(deltaX, deltaS, deltaLambda, 0, deltaX_R, deltaX_Q));
474: }
475: PetscCall(VecAXPY(DeltaX, 1, deltaX));
476: PetscCall(VecAXPY(X, 1, deltaX));
477: /* Q = -dF/dlambda(X, lambda)*/
478: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
479: /* R = F(X, lambda) */
480: PetscCall(SNESComputeFunction(snes, X, R));
481: PetscCall(VecNormBegin(R, NORM_2, &fnorm));
482: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
483: PetscCall(VecNormBegin(deltaX, NORM_2, &ynorm));
484: PetscCall(VecNormEnd(R, NORM_2, &fnorm));
485: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
486: PetscCall(VecNormEnd(deltaX, NORM_2, &ynorm));
488: if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));
490: /* Monitor convergence */
491: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
492: snes->iter++;
493: snes->norm = fnorm;
494: snes->ynorm = ynorm;
495: snes->xnorm = xnorm;
496: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
497: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
498: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
499: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
500: if (!snes->reason && j == maxits - 1) snes->reason = SNES_DIVERGED_MAX_IT;
501: if (snes->reason) break;
502: }
503: if (snes->reason < 0) break;
504: if (data->lambda >= data->lambda_max) {
505: break;
506: } else if (maxincs > 0 && i == maxincs - 1) {
507: snes->reason = SNES_DIVERGED_MAX_IT;
508: break;
509: } else {
510: snes->reason = SNES_CONVERGED_ITERATING;
511: }
512: }
513: /* Reset RHS vector, if changed */
514: if (data->copied_rhs) {
515: PetscCall(VecCopy(data->vec_rhs_orig, snes->vec_rhs));
516: data->copied_rhs = PETSC_FALSE;
517: }
518: snes->max_its = maxits; /* reset snes->max_its */
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*
523: SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
524: of the SNESNEWTONAL nonlinear solver.
526: Input Parameter:
527: . snes - the SNES context
528: . x - the solution vector
530: Application Interface Routine: SNESSetUp()
531: */
532: static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
533: {
534: PetscFunctionBegin;
535: PetscCall(SNESSetWorkVecs(snes, 4));
536: PetscCall(SNESSetUpMatrices(snes));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*
541: SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.
543: Input Parameter:
544: . snes - the SNES context
546: Application Interface Routine: SNESSetFromOptions()
547: */
548: static PetscErrorCode SNESSetFromOptions_NEWTONAL(SNES snes, PetscOptionItems PetscOptionsObject)
549: {
550: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
551: SNESNewtonALCorrectionType correction_type = data->correction_type;
553: PetscFunctionBegin;
554: PetscOptionsHeadBegin(PetscOptionsObject, "SNES Newton Arc Length options");
555: PetscCall(PetscOptionsReal("-snes_newtonal_step_size", "Initial arc length increment step size", "SNESNewtonAL", data->step_size, &data->step_size, NULL));
556: PetscCall(PetscOptionsInt("-snes_newtonal_max_continuation_steps", "Maximum number of increment steps", "SNESNewtonAL", data->max_continuation_steps, &data->max_continuation_steps, NULL));
557: PetscCall(PetscOptionsReal("-snes_newtonal_psisq", "Regularization parameter for arc length continuation, 0 for cylindrical", "SNESNewtonAL", data->psisq, &data->psisq, NULL));
558: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_min", "Minimum value of the load parameter lambda", "SNESNewtonAL", data->lambda_min, &data->lambda_min, NULL));
559: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_max", "Maximum value of the load parameter lambda", "SNESNewtonAL", data->lambda_max, &data->lambda_max, NULL));
560: 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));
561: 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));
562: PetscCall(SNESNewtonALSetCorrectionType(snes, correction_type));
563: PetscOptionsHeadEnd();
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
568: {
569: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
571: PetscFunctionBegin;
572: PetscCall(VecDestroy(&al->vec_rhs_orig));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*
577: SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
578: with SNESCreate_NEWTONAL().
580: Input Parameter:
581: . snes - the SNES context
583: Application Interface Routine: SNESDestroy()
584: */
585: static PetscErrorCode SNESDestroy_NEWTONAL(SNES snes)
586: {
587: PetscFunctionBegin;
588: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", NULL));
589: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", NULL));
590: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", NULL));
591: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", NULL));
592: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", NULL));
593: PetscCall(PetscFree(snes->data));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*MC
598: SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.
600: Options Database Keys:
601: + -snes_newtonal_step_size <1.0> - Initial arc length increment step size
602: . -snes_newtonal_max_continuation_steps <100> - Maximum number of continuation steps, or negative for no limit (not recommended)
603: . -snes_newtonal_psisq <1.0> - Regularization parameter for arc length continuation, 0 for cylindrical. Larger values generally lead to more steps.
604: . -snes_newtonal_lambda_min <0.0> - Minimum value of the load parameter lambda
605: . -snes_newtonal_lambda_max <1.0> - Maximum value of the load parameter lambda
606: . -snes_newtonal_scale_rhs <true> - Scale the constant vector passed to `SNESSolve` by the load parameter lambda
607: - -snes_newtonal_correction_type <exact> - Type of correction to use in the arc-length continuation method, `exact` or `normal`
609: Level: intermediate
611: Note:
612: The exact correction scheme with partial updates is detailed in {cite}`Ritto-CorreaCamotim2008` and the implementation of the
613: normal correction scheme is based on {cite}`LeonPaulinoPereiraMenezesLages_2011`.
615: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
616: M*/
617: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
618: {
619: SNES_NEWTONAL *arclengthParameters;
621: PetscFunctionBegin;
622: snes->ops->setup = SNESSetUp_NEWTONAL;
623: snes->ops->solve = SNESSolve_NEWTONAL;
624: snes->ops->destroy = SNESDestroy_NEWTONAL;
625: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONAL;
626: snes->ops->reset = SNESReset_NEWTONAL;
628: snes->usesksp = PETSC_TRUE;
629: snes->usesnpc = PETSC_FALSE;
631: PetscCall(SNESParametersInitialize(snes));
632: PetscObjectParameterSetDefault(snes, max_funcs, 30000);
633: PetscObjectParameterSetDefault(snes, max_its, 50);
635: snes->alwayscomputesfinalresidual = PETSC_TRUE;
637: PetscCall(PetscNew(&arclengthParameters));
638: arclengthParameters->lambda = 0.0;
639: arclengthParameters->lambda_update = 0.0;
640: arclengthParameters->step_size = 1.0;
641: arclengthParameters->max_continuation_steps = 100;
642: arclengthParameters->psisq = 1.0;
643: arclengthParameters->lambda_min = 0.0;
644: arclengthParameters->lambda_max = 1.0;
645: arclengthParameters->scale_rhs = PETSC_TRUE;
646: arclengthParameters->correction_type = SNES_NEWTONAL_CORRECTION_EXACT;
647: snes->data = (void *)arclengthParameters;
649: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", SNESNewtonALSetCorrectionType_NEWTONAL));
650: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", SNESNewtonALGetLoadParameter_NEWTONAL));
651: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", SNESNewtonALSetFunction_NEWTONAL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", SNESNewtonALGetFunction_NEWTONAL));
653: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", SNESNewtonALComputeFunction_NEWTONAL));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }