1 #include <../src/snes/impls/al/alimpl.h> /*I "petscsnes.h" I*/
2
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 */
9
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;
37
38 const char *const SNESNewtonALCorrectionTypes[] = {"EXACT", "NORMAL", "SNESNewtonALCorrectionType", "SNES_NEWTONAL_CORRECTION_", NULL};
39
SNESNewtonALCheckArcLength(SNES snes,Vec XStep,PetscReal lambdaStep,PetscReal stepSize)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;
44
45 PetscFunctionBegin;
46 PetscCall(VecNorm(XStep, NORM_2, &arcLength));
47 arcLength = PetscSqrtReal(PetscSqr(arcLength) + al->psisq * lambdaStep * lambdaStep);
48 arcLengthError = PetscAbsReal(arcLength - stepSize);
49
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 }
53
54 /* stable implementation of roots of a*x^2 + b*x + c = 0 */
PetscQuadraticRoots(PetscReal a,PetscReal b,PetscReal c,PetscReal * xm,PetscReal * xp)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 }
63
SNESNewtonALSetCorrectionType_NEWTONAL(SNES snes,SNESNewtonALCorrectionType ctype)64 static PetscErrorCode SNESNewtonALSetCorrectionType_NEWTONAL(SNES snes, SNESNewtonALCorrectionType ctype)
65 {
66 SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
67
68 PetscFunctionBegin;
69 al->correction_type = ctype;
70 PetscFunctionReturn(PETSC_SUCCESS);
71 }
72
73 /*@
74 SNESNewtonALSetCorrectionType - Set the type of correction to use in the arc-length continuation method.
75
76 Logically Collective
77
78 Input Parameters:
79 + snes - the nonlinear solver object
80 - ctype - the type of correction to use
81
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
84
85 Level: intermediate
86
87 .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALCorrectionType`
88 @*/
SNESNewtonALSetCorrectionType(SNES snes,SNESNewtonALCorrectionType ctype)89 PetscErrorCode SNESNewtonALSetCorrectionType(SNES snes, SNESNewtonALCorrectionType ctype)
90 {
91 PetscFunctionBegin;
92 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
93 PetscValidLogicalCollectiveEnum(snes, ctype, 2);
94 PetscTryMethod(snes, "SNESNewtonALSetCorrectionType_C", (SNES, SNESNewtonALCorrectionType), (snes, ctype));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
SNESNewtonALSetFunction_NEWTONAL(SNES snes,SNESFunctionFn * func,PetscCtx ctx)98 static PetscErrorCode SNESNewtonALSetFunction_NEWTONAL(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
99 {
100 SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
101
102 PetscFunctionBegin;
103 al->computealfunction = func;
104 al->alctx = ctx;
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
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.
111
112 Logically Collective
113
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`)
118
119 Level: intermediate
120
121 Notes:
122 If the current value of the load parameter is needed in `func`, it can be obtained with `SNESNewtonALGetLoadParameter()`.
123
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.
126
127 .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
128 @*/
SNESNewtonALSetFunction(SNES snes,SNESFunctionFn * func,PetscCtx ctx)129 PetscErrorCode SNESNewtonALSetFunction(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
130 {
131 PetscFunctionBegin;
132 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
133 PetscTryMethod(snes, "SNESNewtonALSetFunction_C", (SNES, SNESFunctionFn *, void *), (snes, func, ctx));
134 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136
SNESNewtonALGetFunction_NEWTONAL(SNES snes,SNESFunctionFn ** func,PetscCtxRt ctx)137 static PetscErrorCode SNESNewtonALGetFunction_NEWTONAL(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
138 {
139 SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
140
141 PetscFunctionBegin;
142 if (func) *func = al->computealfunction;
143 if (ctx) *(void **)ctx = al->alctx;
144 PetscFunctionReturn(PETSC_SUCCESS);
145 }
146
147 /*@C
148 SNESNewtonALGetFunction - Get the user function and context set with `SNESNewtonALSetFunction`
149
150 Logically Collective
151
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`)
156
157 Level: intermediate
158
159 .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
160 @*/
SNESNewtonALGetFunction(SNES snes,SNESFunctionFn ** func,PetscCtxRt ctx)161 PetscErrorCode SNESNewtonALGetFunction(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
162 {
163 PetscFunctionBegin;
164 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
165 PetscUseMethod(snes, "SNESNewtonALGetFunction_C", (SNES, SNESFunctionFn **, PetscCtxRt), (snes, func, ctx));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes,PetscReal * lambda)169 static PetscErrorCode SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes, PetscReal *lambda)
170 {
171 SNES_NEWTONAL *al;
172
173 PetscFunctionBeginHot;
174 al = (SNES_NEWTONAL *)snes->data;
175 *lambda = al->lambda;
176 PetscFunctionReturn(PETSC_SUCCESS);
177 }
178
179 /*@C
180 SNESNewtonALGetLoadParameter - Get the value of the load parameter `lambda` for the arc-length continuation method.
181
182 Logically Collective
183
184 Input Parameter:
185 . snes - the nonlinear solver object
186
187 Output Parameter:
188 . lambda - the arc-length parameter
189
190 Level: intermediate
191
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).
195
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.
198
199 .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
200 @*/
SNESNewtonALGetLoadParameter(SNES snes,PetscReal * lambda)201 PetscErrorCode SNESNewtonALGetLoadParameter(SNES snes, PetscReal *lambda)
202 {
203 PetscFunctionBeginHot;
204 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
205 PetscAssertPointer(lambda, 2);
206 PetscUseMethod(snes, "SNESNewtonALGetLoadParameter_C", (SNES, PetscReal *), (snes, lambda));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
SNESNewtonALComputeFunction_NEWTONAL(SNES snes,Vec X,Vec Q)210 static PetscErrorCode SNESNewtonALComputeFunction_NEWTONAL(SNES snes, Vec X, Vec Q)
211 {
212 void *ctx = NULL;
213 SNESFunctionFn *computealfunction = NULL;
214 SNES_NEWTONAL *al;
215
216 PetscFunctionBegin;
217 al = (SNES_NEWTONAL *)snes->data;
218 PetscCall(SNESNewtonALGetFunction(snes, &computealfunction, &ctx));
219
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 }
239
240 /*@C
241 SNESNewtonALComputeFunction - Calls the function that has been set with `SNESNewtonALSetFunction()`.
242
243 Collective
244
245 Input Parameters:
246 + snes - the `SNES` context
247 - X - input vector
248
249 Output Parameter:
250 . Q - tangent load vector, as set by `SNESNewtonALSetFunction()`
251
252 Level: developer
253
254 Note:
255 `SNESNewtonALComputeFunction()` is typically used within nonlinear solvers
256 implementations, so users would not generally call this routine themselves.
257
258 .seealso: [](ch_snes), `SNES`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`
259 @*/
SNESNewtonALComputeFunction(SNES snes,Vec X,Vec Q)260 PetscErrorCode SNESNewtonALComputeFunction(SNES snes, Vec X, Vec Q)
261 {
262 PetscFunctionBegin;
263 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
264 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
265 PetscValidHeaderSpecific(Q, VEC_CLASSID, 3);
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 }
274
275 /*
276 SNESSolve_NEWTONAL - Solves a nonlinear system with Newton's method with arc length continuation.
277
278 Input Parameter:
279 . snes - the `SNES` context
280
281 Application Interface Routine: SNESSolve()
282 */
SNESSolve_NEWTONAL(SNES snes)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;
289
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);
292
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 }
300
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;
308
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 */
319
320 PetscCall(VecZeroEntries(DeltaX));
321
322 /* set snes->max_its for convergence test */
323 snes->max_its = maxits * maxincs;
324
325 /* main incremental-iterative loop */
326 for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
327 PetscReal deltaLambda;
328
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 SNESCheckFunctionDomainError(snes, fnorm);
337
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;
344
345 /* Call general purpose update function */
346 PetscTryTypeMethod(snes, update, snes->iter);
347
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));
377
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};
387
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};
413
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));
419
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]);
426
427 as = b1 * b1 - 4 * a0 * c2;
428 bs = 2 * (b1 * b0 - 2 * a0 * c1);
429 cs = b0 * b0 - 4 * a0 * c0;
430
431 discriminant = cs + bs * deltaS + as * deltaS * deltaS;
432
433 if (discriminant < 0) {
434 /* Take deltaS < 1 with the unique root -b/(2*a) */
435 PetscReal x1;
436
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;
444
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));
487 SNESCheckFunctionDomainError(snes, fnorm);
488 if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));
489
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 }
521
522 /*
523 SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
524 of the SNESNEWTONAL nonlinear solver.
525
526 Input Parameter:
527 . snes - the SNES context
528 . x - the solution vector
529
530 Application Interface Routine: SNESSetUp()
531 */
SNESSetUp_NEWTONAL(SNES snes)532 static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
533 {
534 PetscFunctionBegin;
535 PetscCall(SNESSetWorkVecs(snes, 4));
536 PetscCall(SNESSetUpMatrices(snes));
537 PetscFunctionReturn(PETSC_SUCCESS);
538 }
539
540 /*
541 SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.
542
543 Input Parameter:
544 . snes - the SNES context
545
546 Application Interface Routine: SNESSetFromOptions()
547 */
SNESSetFromOptions_NEWTONAL(SNES snes,PetscOptionItems PetscOptionsObject)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;
552
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 }
566
SNESReset_NEWTONAL(SNES snes)567 static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
568 {
569 SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;
570
571 PetscFunctionBegin;
572 PetscCall(VecDestroy(&al->vec_rhs_orig));
573 PetscFunctionReturn(PETSC_SUCCESS);
574 }
575
576 /*
577 SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
578 with SNESCreate_NEWTONAL().
579
580 Input Parameter:
581 . snes - the SNES context
582
583 Application Interface Routine: SNESDestroy()
584 */
SNESDestroy_NEWTONAL(SNES snes)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 }
596
597 /*MC
598 SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.
599
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`
608
609 Level: intermediate
610
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`.
614
615 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
616 M*/
SNESCreate_NEWTONAL(SNES snes)617 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
618 {
619 SNES_NEWTONAL *arclengthParameters;
620
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;
627
628 snes->usesksp = PETSC_TRUE;
629 snes->usesnpc = PETSC_FALSE;
630
631 PetscCall(SNESParametersInitialize(snes));
632 PetscObjectParameterSetDefault(snes, max_funcs, 30000);
633 PetscObjectParameterSetDefault(snes, max_its, 50);
634
635 snes->alwayscomputesfinalresidual = PETSC_TRUE;
636
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;
648
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 }
656