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;
308: maxits = snes->max_its; /* maximum number of iterations */
309: maxincs = data->max_continuation_steps; /* maximum number of increments */
310: X = snes->vec_sol; /* solution vector */
311: R = snes->vec_func; /* residual vector */
312: Q = snes->work[0]; /* tangent load vector */
313: deltaX_Q = snes->work[1]; /* variation of X with respect to lambda */
314: deltaX_R = snes->work[2]; /* linearized error correction */
315: DeltaX = snes->work[3]; /* step from equilibrium */
316: deltaX = snes->vec_sol_update; /* full newton step */
317: stepSize = data->step_size; /* initial step size */
319: PetscCall(VecZeroEntries(DeltaX));
321: /* set snes->max_its for convergence test */
322: snes->max_its = maxits * maxincs;
324: /* main incremental-iterative loop */
325: for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
326: PetscReal deltaLambda;
328: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
329: snes->norm = 0.0;
330: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
331: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
332: PetscCall(SNESComputeFunction(snes, X, R));
333: PetscCall(VecAXPY(R, 1, Q)); /* R <- R + Q */
334: PetscCall(VecNorm(R, NORM_2, &fnorm)); /* fnorm <- ||R|| */
335: SNESCheckFunctionNorm(snes, fnorm);
337: /* Monitor convergence */
338: PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
339: PetscCall(SNESMonitor(snes, snes->iter, fnorm));
340: if (i == 0 && snes->reason) break;
341: for (PetscInt j = 0; j < maxits; j++) {
342: PetscReal normsqX_Q, deltaS = 1;
344: /* Call general purpose update function */
345: PetscTryTypeMethod(snes, update, snes->iter);
347: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
348: SNESCheckJacobianDomainerror(snes);
349: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
350: /* Solve J deltaX_Q = Q, where J is Jacobian matrix */
351: PetscCall(KSPSolve(snes->ksp, Q, deltaX_Q));
352: SNESCheckKSPSolve(snes);
353: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
354: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", tangent load linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
355: /* Compute load parameter variation */
356: PetscCall(VecNorm(deltaX_Q, NORM_2, &normsqX_Q));
357: normsqX_Q *= normsqX_Q;
358: /* On first iter, use predictor. This is the same regardless of corrector scheme. */
359: if (j == 0) {
360: PetscReal sign = 1.0;
361: if (i > 0) {
362: PetscCall(VecDotRealPart(DeltaX, deltaX_Q, &sign));
363: sign += data->psisq * data->lambda_update;
364: sign = sign >= 0 ? 1.0 : -1.0;
365: }
366: data->lambda_update = 0.0;
367: PetscCall(VecZeroEntries(DeltaX));
368: deltaLambda = sign * stepSize / PetscSqrtReal(normsqX_Q + data->psisq);
369: } else {
370: /* Solve J deltaX_R = -R */
371: PetscCall(KSPSolve(snes->ksp, R, deltaX_R));
372: SNESCheckKSPSolve(snes);
373: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
374: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", residual linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
375: PetscCall(VecScale(deltaX_R, -1));
377: if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
378: /*
379: Take a step orthogonal to the current incremental update DeltaX.
380: Note, this approach is cheaper than the exact correction, but may exhibit convergence
381: issues due to the iterative trial points not being on the quadratic contraint surface.
382: On the bright side, we always have a real and unique solution for deltaLambda.
383: */
384: PetscScalar coefs[2];
385: Vec rhs[] = {deltaX_R, deltaX_Q};
387: PetscCall(VecMDot(DeltaX, 2, rhs, coefs));
388: deltaLambda = -PetscRealPart(coefs[0]) / (PetscRealPart(coefs[1]) + data->psisq * data->lambda_update);
389: } else {
390: /*
391: Solve
392: a*deltaLambda^2 + b*deltaLambda + c = 0 (*)
393: where
394: a = a0
395: b = b0 + b1*deltaS
396: c = c0 + c1*deltaS + c2*deltaS^2
397: and deltaS is either 1, or the largest value in (0, 1) that satisfies
398: b^2 - 4*a*c = as*deltaS^2 + bs*deltaS + cs >= 0
399: where
400: as = b1^2 - 4*a0*c2
401: bs = 2*b1*b0 - 4*a0*c1
402: cs = b0^2 - 4*a0*c0
403: These "partial corrections" prevent (*) from having complex roots.
404: */
405: PetscReal psisqLambdaUpdate, discriminant;
406: PetscReal as, bs, cs;
407: PetscReal a0, b0, b1, c0, c1, c2;
408: PetscScalar coefs1[3]; /* coefs[0] = deltaX_Q*DeltaX, coefs[1] = deltaX_R*DeltaX, coefs[2] = DeltaX*DeltaX */
409: PetscScalar coefs2[2]; /* coefs[0] = deltaX_Q*deltaX_R, coefs[1] = deltaX_R*deltaX_R */
410: const Vec rhs1[3] = {deltaX_Q, deltaX_R, DeltaX};
411: const Vec rhs2[2] = {deltaX_Q, deltaX_R};
413: psisqLambdaUpdate = data->psisq * data->lambda_update;
414: PetscCall(VecMDotBegin(DeltaX, 3, rhs1, coefs1));
415: PetscCall(VecMDotBegin(deltaX_R, 2, rhs2, coefs2));
416: PetscCall(VecMDotEnd(DeltaX, 3, rhs1, coefs1));
417: PetscCall(VecMDotEnd(deltaX_R, 2, rhs2, coefs2));
419: a0 = normsqX_Q + data->psisq;
420: b0 = 2 * (PetscRealPart(coefs1[0]) + psisqLambdaUpdate);
421: b1 = 2 * PetscRealPart(coefs2[0]);
422: c0 = PetscRealPart(coefs1[2]) + psisqLambdaUpdate * data->lambda_update - stepSize * stepSize;
423: c1 = 2 * PetscRealPart(coefs1[1]);
424: c2 = PetscRealPart(coefs2[1]);
426: as = b1 * b1 - 4 * a0 * c2;
427: bs = 2 * (b1 * b0 - 2 * a0 * c1);
428: cs = b0 * b0 - 4 * a0 * c0;
430: discriminant = cs + bs * deltaS + as * deltaS * deltaS;
432: if (discriminant < 0) {
433: /* Take deltaS < 1 with the unique root -b/(2*a) */
434: PetscReal x1;
436: /* Compute deltaS to be the largest root of (as * x^2 + bs * x + cs = 0) */
437: PetscQuadraticRoots(as, bs, cs, &x1, &deltaS);
438: 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));
439: deltaLambda = -0.5 * (b0 + b1 * deltaS) / a0;
440: } else {
441: /* Use deltaS = 1, pick root that is closest to the last point to prevent doubling back */
442: PetscReal dlambda1, dlambda2;
444: PetscQuadraticRoots(a0, b0 + b1, c0 + c1 + c2, &dlambda1, &dlambda2);
445: deltaLambda = (b0 * dlambda1) > (b0 * dlambda2) ? dlambda1 : dlambda2;
446: }
447: }
448: }
449: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
450: data->lambda = data->lambda + deltaLambda;
451: if (data->lambda > data->lambda_max) {
452: /* Ensure that lambda = lambda_max exactly at the end of incremental process. This ensures the final solution matches the problem we want to solve. */
453: deltaLambda = deltaLambda - (data->lambda - data->lambda_max);
454: data->lambda = data->lambda_max;
455: }
456: if (data->lambda < data->lambda_min) {
457: // LCOV_EXCL_START
458: /* Ensure that lambda >= lambda_min. This prevents some potential oscillatory behavior. */
459: deltaLambda = deltaLambda - (data->lambda - data->lambda_min);
460: data->lambda = data->lambda_min;
461: // LCOV_EXCL_STOP
462: }
463: data->lambda_update = data->lambda_update + deltaLambda;
464: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
465: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", lambda=%18.16e, lambda_update=%18.16e\n", snes->iter, (double)data->lambda, (double)data->lambda_update));
466: if (j == 0) {
467: /* deltaX = deltaLambda*deltaX_Q */
468: PetscCall(VecCopy(deltaX_Q, deltaX));
469: PetscCall(VecScale(deltaX, deltaLambda));
470: } else {
471: /* deltaX = deltaS*deltaX_R + deltaLambda*deltaX_Q */
472: PetscCall(VecAXPBYPCZ(deltaX, deltaS, deltaLambda, 0, deltaX_R, deltaX_Q));
473: }
474: PetscCall(VecAXPY(DeltaX, 1, deltaX));
475: PetscCall(VecAXPY(X, 1, deltaX));
476: /* Q = -dF/dlambda(X, lambda)*/
477: PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
478: /* R = F(X, lambda) */
479: PetscCall(SNESComputeFunction(snes, X, R));
480: PetscCall(VecNormBegin(R, NORM_2, &fnorm));
481: PetscCall(VecNormBegin(X, NORM_2, &xnorm));
482: PetscCall(VecNormBegin(deltaX, NORM_2, &ynorm));
483: PetscCall(VecNormEnd(R, NORM_2, &fnorm));
484: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
485: PetscCall(VecNormEnd(deltaX, NORM_2, &ynorm));
487: if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));
489: /* Monitor convergence */
490: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
491: snes->iter++;
492: snes->norm = fnorm;
493: snes->ynorm = ynorm;
494: snes->xnorm = xnorm;
495: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
496: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
497: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
498: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
499: if (snes->reason) break;
500: if (j == maxits - 1) {
501: snes->reason = SNES_DIVERGED_MAX_IT;
502: break;
503: }
504: }
505: if (snes->reason < 0) break;
506: if (data->lambda >= data->lambda_max) {
507: break;
508: } else if (maxincs > 0 && i == maxincs - 1) {
509: snes->reason = SNES_DIVERGED_MAX_IT;
510: break;
511: } else {
512: snes->reason = SNES_CONVERGED_ITERATING;
513: }
514: }
515: /* Reset RHS vector, if changed */
516: if (data->copied_rhs) {
517: PetscCall(VecCopy(data->vec_rhs_orig, snes->vec_rhs));
518: data->copied_rhs = PETSC_FALSE;
519: }
520: snes->max_its = maxits; /* reset snes->max_its */
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*
525: SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
526: of the SNESNEWTONAL nonlinear solver.
528: Input Parameter:
529: . snes - the SNES context
530: . x - the solution vector
532: Application Interface Routine: SNESSetUp()
533: */
534: static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
535: {
536: PetscFunctionBegin;
537: PetscCall(SNESSetWorkVecs(snes, 4));
538: PetscCall(SNESSetUpMatrices(snes));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*
543: SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.
545: Input Parameter:
546: . snes - the SNES context
548: Application Interface Routine: SNESSetFromOptions()
549: */
550: static PetscErrorCode SNESSetFromOptions_NEWTONAL(SNES snes, PetscOptionItems *PetscOptionsObject)
551: {
552: SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
553: SNESNewtonALCorrectionType correction_type = data->correction_type;
555: PetscFunctionBegin;
556: PetscOptionsHeadBegin(PetscOptionsObject, "SNES Newton Arc Length options");
557: PetscCall(PetscOptionsReal("-snes_newtonal_step_size", "Initial arc length increment step size", "SNESNewtonAL", data->step_size, &data->step_size, NULL));
558: PetscCall(PetscOptionsInt("-snes_newtonal_max_continuation_steps", "Maximum number of increment steps", "SNESNewtonAL", data->max_continuation_steps, &data->max_continuation_steps, NULL));
559: PetscCall(PetscOptionsReal("-snes_newtonal_psisq", "Regularization parameter for arc length continuation, 0 for cylindrical", "SNESNewtonAL", data->psisq, &data->psisq, NULL));
560: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_min", "Minimum value of the load parameter lambda", "SNESNewtonAL", data->lambda_min, &data->lambda_min, NULL));
561: PetscCall(PetscOptionsReal("-snes_newtonal_lambda_max", "Maximum value of the load parameter lambda", "SNESNewtonAL", data->lambda_max, &data->lambda_max, NULL));
562: 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));
563: 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));
564: PetscCall(SNESNewtonALSetCorrectionType(snes, correction_type));
565: PetscOptionsHeadEnd();
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
570: {
571: SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
573: PetscFunctionBegin;
574: PetscCall(VecDestroy(&al->vec_rhs_orig));
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: /*
579: SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
580: with SNESCreate_NEWTONAL().
582: Input Parameter:
583: . snes - the SNES context
585: Application Interface Routine: SNESDestroy()
586: */
587: static PetscErrorCode SNESDestroy_NEWTONAL(SNES snes)
588: {
589: PetscFunctionBegin;
590: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", NULL));
591: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", NULL));
592: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", NULL));
593: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", NULL));
594: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", NULL));
595: PetscCall(PetscFree(snes->data));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*MC
600: SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.
602: Options Database Keys:
603: + -snes_newtonal_step_size <1.0> - Initial arc length increment step size
604: . -snes_newtonal_max_continuation_steps <100> - Maximum number of continuation steps, or negative for no limit (not recommended)
605: . -snes_newtonal_psisq <1.0> - Regularization parameter for arc length continuation, 0 for cylindrical. Larger values generally lead to more steps.
606: . -snes_newtonal_lambda_min <0.0> - Minimum value of the load parameter lambda
607: . -snes_newtonal_lambda_max <1.0> - Maximum value of the load parameter lambda
608: . -snes_newtonal_scale_rhs <true> - Scale the constant vector passed to `SNESSolve` by the load parameter lambda
609: - -snes_newtonal_correction_type <exact> - Type of correction to use in the arc-length continuation method, `exact` or `normal`
611: Level: intermediate
613: Note:
614: The exact correction scheme with partial updates is detailed in {cite}`Ritto-CorreaCamotim2008` and the implementation of the
615: normal correction scheme is based on {cite}`LeonPaulinoPereiraMenezesLages_2011`.
617: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
618: M*/
619: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
620: {
621: SNES_NEWTONAL *arclengthParameters;
623: PetscFunctionBegin;
624: snes->ops->setup = SNESSetUp_NEWTONAL;
625: snes->ops->solve = SNESSolve_NEWTONAL;
626: snes->ops->destroy = SNESDestroy_NEWTONAL;
627: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONAL;
628: snes->ops->reset = SNESReset_NEWTONAL;
630: snes->usesksp = PETSC_TRUE;
631: snes->usesnpc = PETSC_FALSE;
633: PetscCall(SNESParametersInitialize(snes));
634: PetscObjectParameterSetDefault(snes, max_funcs, 30000);
635: PetscObjectParameterSetDefault(snes, max_its, 50);
637: snes->alwayscomputesfinalresidual = PETSC_TRUE;
639: PetscCall(PetscNew(&arclengthParameters));
640: arclengthParameters->lambda = 0.0;
641: arclengthParameters->lambda_update = 0.0;
642: arclengthParameters->step_size = 1.0;
643: arclengthParameters->max_continuation_steps = 100;
644: arclengthParameters->psisq = 1.0;
645: arclengthParameters->lambda_min = 0.0;
646: arclengthParameters->lambda_max = 1.0;
647: arclengthParameters->scale_rhs = PETSC_TRUE;
648: arclengthParameters->correction_type = SNES_NEWTONAL_CORRECTION_EXACT;
649: snes->data = (void *)arclengthParameters;
651: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", SNESNewtonALSetCorrectionType_NEWTONAL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", SNESNewtonALGetLoadParameter_NEWTONAL));
653: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", SNESNewtonALSetFunction_NEWTONAL));
654: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", SNESNewtonALGetFunction_NEWTONAL));
655: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", SNESNewtonALComputeFunction_NEWTONAL));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }