Actual source code: al.c

      1: #include <../src/snes/impls/al/alimpl.h>

      3: /*
      4:      This file implements a truncated Newton method with arc length continuation,
      5:      for solving a system of nonlinear equations, using the KSP, Vec,
      6:      and Mat interfaces for linear solvers, vectors, and matrices,
      7:      respectively.
      8: */

     10: const char NewtonALExactCitation[]   = "@article{Ritto-CorreaCamotim2008,\n"
     11:                                        "  title={On the arc-length and other quadratic control methods: Established, less known and new implementation procedures},\n"
     12:                                        "  volume={86},\n"
     13:                                        "  ISSN={0045-7949},\n"
     14:                                        "  DOI={10.1016/j.compstruc.2007.08.003},\n"
     15:                                        "  number={11},\n"
     16:                                        "  journal={Computers & Structures},\n"
     17:                                        "  author={Ritto-Corr{\\^{e}}a, Manuel and Camotim, Dinar},\n"
     18:                                        "  year={2008},\n"
     19:                                        "  month=jun,\n"
     20:                                        "  pages={1353-1368},\n"
     21:                                        "}\n";
     22: PetscBool  NewtonALExactCitationSet  = PETSC_FALSE;
     23: const char NewtonALNormalCitation[]  = "@article{LeonPaulinoPereiraMenezesLages_2011,\n"
     24:                                        "  title={A Unified Library of Nonlinear Solution Schemes},\n"
     25:                                        "  volume={64},\n"
     26:                                        "  ISSN={0003-6900, 2379-0407},\n"
     27:                                        "  DOI={10.1115/1.4006992},\n"
     28:                                        "  number={4},\n"
     29:                                        "  journal={Applied Mechanics Reviews},\n"
     30:                                        "  author={Leon, Sofie E. and Paulino, Glaucio H. and Pereira, Anderson and Menezes, Ivan F. M. and Lages, Eduardo N.},\n"
     31:                                        "  year={2011},\n"
     32:                                        "  month=jul,\n"
     33:                                        "  pages={040803},\n"
     34:                                        "  language={en}\n"
     35:                                        "}\n";
     36: PetscBool  NewtonALNormalCitationSet = PETSC_FALSE;

     38: const char *const SNESNewtonALCorrectionTypes[] = {"EXACT", "NORMAL", "SNESNewtonALCorrectionType", "SNES_NEWTONAL_CORRECTION_", NULL};

     40: static PetscErrorCode SNESNewtonALCheckArcLength(SNES snes, Vec XStep, PetscReal lambdaStep, PetscReal stepSize)
     41: {
     42:   PetscReal      arcLength, arcLengthError;
     43:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

     45:   PetscFunctionBegin;
     46:   if (al->mat_diag_scaling) PetscCall(MatANorm(al->mat_diag_scaling, XStep, &arcLength));
     47:   else PetscCall(VecNorm(XStep, NORM_2, &arcLength));
     48:   arcLength      = PetscSqrtReal(arcLength * arcLength + al->psisq * lambdaStep * lambdaStep);
     49:   arcLengthError = PetscAbsReal(arcLength - stepSize);

     51:   if (arcLengthError > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscInfo(snes, "Arc length differs from specified step size: computed=%18.16e, expected=%18.16e, error=%18.16e \n", (double)arcLength, (double)stepSize, (double)arcLengthError));
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
     56: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
     57: {
     58:   PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
     59:   PetscReal x1   = temp / a;
     60:   PetscReal x2   = c / temp;
     61:   *xm            = PetscMin(x1, x2);
     62:   *xp            = PetscMax(x1, x2);
     63: }

     65: static PetscErrorCode SNESNewtonALSetCorrectionType_NEWTONAL(SNES snes, SNESNewtonALCorrectionType ctype)
     66: {
     67:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

     69:   PetscFunctionBegin;
     70:   al->correction_type = ctype;
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*@
     75:   SNESNewtonALSetCorrectionType - Set the type of correction to use in the arc-length continuation method.

     77:   Logically Collective

     79:   Input Parameters:
     80: + snes  - the nonlinear solver object
     81: - ctype - the type of correction to use

     83:   Options Database Key:
     84: . -snes_newtonal_correction_type type - Set the type of correction to use; use -help for a list of available types

     86:   Level: intermediate

     88: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALCorrectionType`
     89: @*/
     90: PetscErrorCode SNESNewtonALSetCorrectionType(SNES snes, SNESNewtonALCorrectionType ctype)
     91: {
     92:   PetscFunctionBegin;
     95:   PetscTryMethod(snes, "SNESNewtonALSetCorrectionType_C", (SNES, SNESNewtonALCorrectionType), (snes, ctype));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode SNESNewtonALSetFunction_NEWTONAL(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
    100: {
    101:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

    103:   PetscFunctionBegin;
    104:   al->computealfunction = func;
    105:   al->alctx             = ctx;
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: /*@C
    110:   SNESNewtonALSetFunction - Sets a user function that is called at each function evaluation to
    111:   compute the tangent load vector for the arc-length continuation method.

    113:   Logically Collective

    115:   Input Parameters:
    116: + snes - the nonlinear solver object
    117: . func - [optional] tangent load function evaluation routine, see `SNESFunctionFn` for the calling sequence. `U` is the current solution vector, `Q` is the output tangent load vector
    118: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

    120:   Level: intermediate

    122:   Notes:
    123:   If the current value of the load parameter is needed in `func`, it can be obtained with `SNESNewtonALGetLoadParameter()`.

    125:   The tangent load vector is the partial derivative of external load with respect to the load parameter.
    126:   In the case of proportional loading, the tangent load vector is the full external load vector at the end of the load step.

    128: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
    129: @*/
    130: PetscErrorCode SNESNewtonALSetFunction(SNES snes, SNESFunctionFn *func, PetscCtx ctx)
    131: {
    132:   PetscFunctionBegin;
    134:   PetscTryMethod(snes, "SNESNewtonALSetFunction_C", (SNES, SNESFunctionFn *, void *), (snes, func, ctx));
    135:   PetscFunctionReturn(PETSC_SUCCESS);
    136: }

    138: static PetscErrorCode SNESNewtonALGetFunction_NEWTONAL(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
    139: {
    140:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

    142:   PetscFunctionBegin;
    143:   if (func) *func = al->computealfunction;
    144:   if (ctx) *(void **)ctx = al->alctx;
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: /*@C
    149:   SNESNewtonALGetFunction - Get the user function and context set with `SNESNewtonALSetFunction`

    151:   Logically Collective

    153:   Input Parameters:
    154: + snes - the nonlinear solver object
    155: . func - [optional] tangent load function evaluation routine, see `SNESNewtonALSetFunction()` for the call sequence
    156: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

    158:   Level: intermediate

    160: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
    161: @*/
    162: PetscErrorCode SNESNewtonALGetFunction(SNES snes, SNESFunctionFn **func, PetscCtxRt ctx)
    163: {
    164:   PetscFunctionBegin;
    166:   PetscUseMethod(snes, "SNESNewtonALGetFunction_C", (SNES, SNESFunctionFn **, PetscCtxRt), (snes, func, ctx));
    167:   PetscFunctionReturn(PETSC_SUCCESS);
    168: }

    170: /*@C
    171:   SNESNewtonALSetDiagonalScaling - Set the global vector used to rescale DoFs for computation of arc length.

    173:   Logically Collective

    175:   Input Parameters:
    176: + snes - the nonlinear solver object
    177: - v    - the `Vec` containing diagonal scaling for each DoF, must be the same size as the solution vector (may be `NULL`)

    179:   Note:
    180:   This function stores a reference to `v`. Any changes to the vector will be reflected automatically in the arc length computation.

    182:   Level: intermediate

    184: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetLoadParameter()`
    185: @*/
    186: PetscErrorCode SNESNewtonALSetDiagonalScaling(SNES snes, Vec v)
    187: {
    188:   PetscFunctionBegin;
    191:   PetscTryMethod(snes, "SNESNewtonALSetDiagonalScaling_C", (SNES, Vec), (snes, v));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: static PetscErrorCode SNESNewtonALSetDiagonalScaling_NEWTONAL(SNES snes, Vec v)
    196: {
    197:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

    199:   PetscFunctionBegin;
    200:   PetscCall(PetscObjectReference((PetscObject)v));
    201:   PetscCall(MatDestroy(&al->mat_diag_scaling));
    202:   if (v) {
    203:     PetscCall(MatCreateDiagonal(v, &al->mat_diag_scaling));
    204:     PetscCall(PetscObjectDereference((PetscObject)v));
    205:   }
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: static PetscErrorCode SNESNewtonALGetLoadParameter_NEWTONAL(SNES snes, PetscReal *lambda)
    210: {
    211:   SNES_NEWTONAL *al;

    213:   PetscFunctionBeginHot;
    214:   al      = (SNES_NEWTONAL *)snes->data;
    215:   *lambda = al->lambda;
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }

    219: /*@C
    220:   SNESNewtonALGetLoadParameter - Get the value of the load parameter `lambda` for the arc-length continuation method.

    222:   Logically Collective

    224:   Input Parameter:
    225: . snes - the nonlinear solver object

    227:   Output Parameter:
    228: . lambda - the arc-length parameter

    230:   Level: intermediate

    232:   Notes:
    233:   This function should be used in the functions provided to `SNESSetFunction()` and `SNESNewtonALSetFunction()`
    234:   to compute the residual and tangent load vectors for a given value of `lambda` (0 <= lambda <= 1).

    236:   Usually, `lambda` is used to scale the external force vector in the residual function, i.e. proportional loading,
    237:   in which case the tangent load vector is the full external force vector.

    239: .seealso: [](ch_snes), `SNES`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`
    240: @*/
    241: PetscErrorCode SNESNewtonALGetLoadParameter(SNES snes, PetscReal *lambda)
    242: {
    243:   PetscFunctionBeginHot;
    245:   PetscAssertPointer(lambda, 2);
    246:   PetscUseMethod(snes, "SNESNewtonALGetLoadParameter_C", (SNES, PetscReal *), (snes, lambda));
    247:   PetscFunctionReturn(PETSC_SUCCESS);
    248: }

    250: static PetscErrorCode SNESNewtonALComputeFunction_NEWTONAL(SNES snes, Vec X, Vec Q)
    251: {
    252:   void           *ctx               = NULL;
    253:   SNESFunctionFn *computealfunction = NULL;
    254:   SNES_NEWTONAL  *al;

    256:   PetscFunctionBegin;
    257:   al = (SNES_NEWTONAL *)snes->data;
    258:   PetscCall(SNESNewtonALGetFunction(snes, &computealfunction, &ctx));

    260:   PetscCall(VecZeroEntries(Q));
    261:   PetscCheck(computealfunction || (snes->vec_rhs && al->scale_rhs), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No tangent load function or rhs vector has been set");
    262:   if (computealfunction) {
    263:     PetscCall(VecLockReadPush(X));
    264:     PetscCallBack("SNES callback NewtonAL tangent load function", (*computealfunction)(snes, X, Q, ctx));
    265:     PetscCall(VecLockReadPop(X));
    266:   }
    267:   if (al->scale_rhs && snes->vec_rhs) {
    268:     /* Save original RHS vector values, then scale `snes->vec_rhs` by load parameter */
    269:     if (!al->vec_rhs_orig) PetscCall(VecDuplicate(snes->vec_rhs, &al->vec_rhs_orig));
    270:     if (!al->copied_rhs) {
    271:       PetscCall(VecCopy(snes->vec_rhs, al->vec_rhs_orig));
    272:       al->copied_rhs = PETSC_TRUE;
    273:     }
    274:     PetscCall(VecAXPBY(snes->vec_rhs, al->lambda, 0.0, al->vec_rhs_orig));
    275:     PetscCall(VecAXPY(Q, 1, al->vec_rhs_orig));
    276:   }
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: /*@C
    281:   SNESNewtonALComputeFunction - Calls the function that has been set with `SNESNewtonALSetFunction()`.

    283:   Collective

    285:   Input Parameters:
    286: + snes - the `SNES` context
    287: - X    - input vector

    289:   Output Parameter:
    290: . Q - tangent load vector, as set by `SNESNewtonALSetFunction()`

    292:   Level: developer

    294:   Note:
    295:   `SNESNewtonALComputeFunction()` is typically used within nonlinear solvers
    296:   implementations, so users would not generally call this routine themselves.

    298: .seealso: [](ch_snes), `SNES`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`
    299: @*/
    300: PetscErrorCode SNESNewtonALComputeFunction(SNES snes, Vec X, Vec Q)
    301: {
    302:   PetscFunctionBegin;
    306:   PetscCheckSameComm(snes, 1, X, 2);
    307:   PetscCheckSameComm(snes, 1, Q, 3);
    308:   PetscCall(VecValidValues_Internal(X, 2, PETSC_TRUE));
    309:   PetscCall(PetscLogEventBegin(SNES_NewtonALEval, snes, X, Q, 0));
    310:   PetscTryMethod(snes, "SNESNewtonALComputeFunction_C", (SNES, Vec, Vec), (snes, X, Q));
    311:   PetscCall(PetscLogEventEnd(SNES_NewtonALEval, snes, X, Q, 0));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: /*
    316:   SNESSolve_NEWTONAL - Solves a nonlinear system with Newton's method with arc length continuation.

    318:   Input Parameter:
    319: . snes - the `SNES` context

    321:   Application Interface Routine: SNESSolve()
    322: */
    323: static PetscErrorCode SNESSolve_NEWTONAL(SNES snes)
    324: {
    325:   SNES_NEWTONAL *data = (SNES_NEWTONAL *)snes->data;
    326:   PetscInt       maxits, maxincs, lits;
    327:   PetscReal      fnorm, xnorm, ynorm, stepSize;
    328:   Vec            DeltaX, deltaX, X, R, Q, deltaX_Q, deltaX_R;
    329:   Mat            W;

    331:   PetscFunctionBegin;
    332:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

    334:   /* Register citations */
    335:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
    336:   if (data->correction_type == SNES_NEWTONAL_CORRECTION_EXACT) {
    337:     PetscCall(PetscCitationsRegister(NewtonALExactCitation, &NewtonALExactCitationSet));
    338:   } else if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
    339:     PetscCall(PetscCitationsRegister(NewtonALNormalCitation, &NewtonALNormalCitationSet));
    340:   }

    342:   data->copied_rhs             = PETSC_FALSE;
    343:   data->lambda_update          = 0.0;
    344:   data->lambda                 = 0.0;
    345:   snes->numFailures            = 0;
    346:   snes->numLinearSolveFailures = 0;
    347:   snes->reason                 = SNES_CONVERGED_ITERATING;
    348:   snes->iter                   = 0;

    350:   maxits   = snes->max_its;                /* maximum number of iterations */
    351:   maxincs  = data->max_continuation_steps; /* maximum number of increments */
    352:   X        = snes->vec_sol;                /* solution vector */
    353:   R        = snes->vec_func;               /* residual vector */
    354:   Q        = snes->work[0];                /* tangent load vector */
    355:   W        = data->mat_diag_scaling;       /* diagonal scaling for DoFs */
    356:   deltaX_Q = snes->work[1];                /* variation of X with respect to lambda */
    357:   deltaX_R = snes->work[2];                /* linearized error correction */
    358:   DeltaX   = snes->work[3];                /* step from equilibrium */
    359:   deltaX   = snes->vec_sol_update;         /* full newton step */
    360:   stepSize = data->step_size;              /* initial step size */

    362:   PetscCall(VecZeroEntries(DeltaX));

    364:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    365:   /* set snes->max_its for convergence test */
    366:   snes->max_its = maxits * maxincs;
    367:   snes->iter    = 0;
    368:   snes->norm    = 0.0;
    369:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

    371:   /* main incremental-iterative loop */
    372:   for (PetscInt i = 0; i < maxincs || maxincs < 0; i++) {
    373:     PetscReal deltaLambda;

    375:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    376:     snes->norm = 0.0;
    377:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    378:     PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
    379:     PetscCall(SNESComputeFunction(snes, X, R));
    380:     PetscCall(VecAXPY(R, 1, Q));           /* R <- R + Q */
    381:     PetscCall(VecNorm(R, NORM_2, &fnorm)); /* fnorm <- ||R|| */
    382:     SNESCheckFunctionDomainError(snes, fnorm);

    384:     /* Monitor convergence */
    385:     PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
    386:     PetscCall(SNESMonitor(snes, snes->iter, fnorm));
    387:     if (snes->reason) break;
    388:     for (PetscInt j = 0; j < maxits; j++) {
    389:       PetscReal normsqX_Q, deltaS = 1;

    391:       /* Call general purpose update function */
    392:       PetscTryTypeMethod(snes, update, snes->iter);

    394:       PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
    395:       SNESCheckJacobianDomainError(snes);
    396:       PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
    397:       /* Solve J deltaX_Q = Q, where J is Jacobian matrix */
    398:       PetscCall(KSPSolve(snes->ksp, Q, deltaX_Q));
    399:       SNESCheckKSPSolve(snes);
    400:       PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
    401:       PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", tangent load linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
    402:       /* Compute load parameter variation */
    403:       if (W) PetscCall(MatANorm(W, deltaX_Q, &normsqX_Q));
    404:       else PetscCall(VecNorm(deltaX_Q, NORM_2, &normsqX_Q));
    405:       normsqX_Q *= normsqX_Q;
    406:       /* On first iter, use predictor. This is the same regardless of corrector scheme. */
    407:       if (j == 0) {
    408:         PetscReal sign = 1.0;
    409:         if (i > 0) {
    410:           PetscScalar dot;
    411:           if (W) PetscCall(MatADot(W, DeltaX, deltaX_Q, &dot));
    412:           else PetscCall(VecDot(DeltaX, deltaX_Q, &dot));
    413:           sign = PetscRealPart(dot) + data->psisq * data->lambda_update;
    414:           sign = sign >= 0 ? 1.0 : -1.0;
    415:         }
    416:         data->lambda_update = 0.0;
    417:         PetscCall(VecZeroEntries(DeltaX));
    418:         deltaLambda = sign * stepSize / PetscSqrtReal(normsqX_Q + data->psisq);
    419:       } else {
    420:         /* Solve J deltaX_R = -R */
    421:         PetscCall(KSPSolve(snes->ksp, R, deltaX_R));
    422:         SNESCheckKSPSolve(snes);
    423:         PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
    424:         PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", residual linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
    425:         PetscCall(VecScale(deltaX_R, -1));

    427:         if (data->correction_type == SNES_NEWTONAL_CORRECTION_NORMAL) {
    428:           /*
    429:             Take a step orthogonal to the current incremental update DeltaX.
    430:             Note, this approach is cheaper than the exact correction, but may exhibit convergence
    431:             issues due to the iterative trial points not being on the quadratic constraint surface.
    432:             On the bright side, we always have a real and unique solution for deltaLambda.
    433:           */
    434:           PetscScalar coefs[2];
    435:           if (W) {
    436:             PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs[0]));
    437:             PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs[1]));
    438:           } else {
    439:             PetscCall(VecDot(DeltaX, deltaX_R, &coefs[0]));
    440:             PetscCall(VecDot(DeltaX, deltaX_Q, &coefs[1]));
    441:           }
    442:           deltaLambda = -PetscRealPart(coefs[0]) / (PetscRealPart(coefs[1]) + data->psisq * data->lambda_update);
    443:         } else {
    444:           /*
    445:             Solve
    446:               a*deltaLambda^2 + b*deltaLambda + c = 0  (*)
    447:             where
    448:               a = a0
    449:               b = b0 + b1*deltaS
    450:               c = c0 + c1*deltaS + c2*deltaS^2
    451:             and deltaS is either 1, or the largest value in (0, 1) that satisfies
    452:               b^2 - 4*a*c = as*deltaS^2 + bs*deltaS + cs >= 0
    453:             where
    454:               as = b1^2 - 4*a0*c2
    455:               bs = 2*b1*b0 - 4*a0*c1
    456:               cs = b0^2 - 4*a0*c0
    457:             These "partial corrections" prevent (*) from having complex roots.
    458:           */
    459:           PetscReal   psisqLambdaUpdate, discriminant;
    460:           PetscReal   as, bs, cs;
    461:           PetscReal   a0, b0, b1, c0, c1, c2;
    462:           PetscScalar coefs1[3]; /* coefs[0] = deltaX_Q*DeltaX, coefs[1] = deltaX_R*DeltaX, coefs[2] = DeltaX*DeltaX */
    463:           PetscScalar coefs2[2]; /* coefs[0] = deltaX_Q*deltaX_R, coefs[1] = deltaX_R*deltaX_R */

    465:           psisqLambdaUpdate = data->psisq * data->lambda_update;
    466:           if (W) {
    467:             PetscCall(MatADot(W, DeltaX, deltaX_Q, &coefs1[0]));
    468:             PetscCall(MatADot(W, DeltaX, deltaX_R, &coefs1[1]));
    469:             PetscCall(MatADot(W, DeltaX, DeltaX, &coefs1[2]));
    470:             PetscCall(MatADot(W, deltaX_R, deltaX_Q, &coefs2[0]));
    471:             PetscCall(MatADot(W, deltaX_R, deltaX_R, &coefs2[1]));
    472:           } else {
    473:             PetscCall(VecDot(DeltaX, deltaX_Q, &coefs1[0]));
    474:             PetscCall(VecDot(DeltaX, deltaX_R, &coefs1[1]));
    475:             PetscCall(VecDot(DeltaX, DeltaX, &coefs1[2]));
    476:             PetscCall(VecDot(deltaX_R, deltaX_Q, &coefs2[0]));
    477:             PetscCall(VecDot(deltaX_R, deltaX_R, &coefs2[1]));
    478:           }

    480:           a0 = normsqX_Q + data->psisq;
    481:           b0 = 2 * (PetscRealPart(coefs1[0]) + psisqLambdaUpdate);
    482:           b1 = 2 * PetscRealPart(coefs2[0]);
    483:           c0 = PetscRealPart(coefs1[2]) + psisqLambdaUpdate * data->lambda_update - stepSize * stepSize;
    484:           c1 = 2 * PetscRealPart(coefs1[1]);
    485:           c2 = PetscRealPart(coefs2[1]);

    487:           as = b1 * b1 - 4 * a0 * c2;
    488:           bs = 2 * (b1 * b0 - 2 * a0 * c1);
    489:           cs = b0 * b0 - 4 * a0 * c0;

    491:           discriminant = cs + bs * deltaS + as * deltaS * deltaS;

    493:           if (discriminant < 0) {
    494:             /* Take deltaS < 1 with the unique root -b/(2*a) */
    495:             PetscReal x1;

    497:             /* Compute deltaS to be the largest root of (as * x^2 + bs * x + cs = 0) */
    498:             PetscQuadraticRoots(as, bs, cs, &x1, &deltaS);
    499:             PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", discriminant=%18.16e < 0, shrinking residual update size to deltaS = %18.16e\n", snes->iter, (double)discriminant, (double)deltaS));
    500:             deltaLambda = -0.5 * (b0 + b1 * deltaS) / a0;
    501:           } else {
    502:             /* Use deltaS = 1, pick root that is closest to the last point to prevent doubling back */
    503:             PetscReal dlambda1, dlambda2;

    505:             PetscQuadraticRoots(a0, b0 + b1, c0 + c1 + c2, &dlambda1, &dlambda2);
    506:             deltaLambda = (b0 * dlambda1) > (b0 * dlambda2) ? dlambda1 : dlambda2;
    507:           }
    508:         }
    509:       }
    510:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    511:       data->lambda = data->lambda + deltaLambda;
    512:       if (data->lambda > data->lambda_max) {
    513:         /* Ensure that lambda = lambda_max exactly at the end of incremental process. This ensures the final solution matches the problem we want to solve. */
    514:         deltaLambda  = deltaLambda - (data->lambda - data->lambda_max);
    515:         data->lambda = data->lambda_max;
    516:       }
    517:       if (data->lambda < data->lambda_min) {
    518:         // LCOV_EXCL_START
    519:         /* Ensure that lambda >= lambda_min. This prevents some potential oscillatory behavior. */
    520:         deltaLambda  = deltaLambda - (data->lambda - data->lambda_min);
    521:         data->lambda = data->lambda_min;
    522:         // LCOV_EXCL_STOP
    523:       }
    524:       data->lambda_update = data->lambda_update + deltaLambda;
    525:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    526:       PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", lambda=%18.16e, lambda_update=%18.16e\n", snes->iter, (double)data->lambda, (double)data->lambda_update));
    527:       if (j == 0) {
    528:         /* deltaX = deltaLambda*deltaX_Q */
    529:         PetscCall(VecCopy(deltaX_Q, deltaX));
    530:         PetscCall(VecScale(deltaX, deltaLambda));
    531:       } else {
    532:         /* deltaX = deltaS*deltaX_R + deltaLambda*deltaX_Q */
    533:         PetscCall(VecAXPBYPCZ(deltaX, deltaS, deltaLambda, 0, deltaX_R, deltaX_Q));
    534:       }
    535:       PetscCall(VecAXPY(DeltaX, 1, deltaX));
    536:       PetscCall(VecAXPY(X, 1, deltaX));
    537:       /* Q = -dF/dlambda(X, lambda)*/
    538:       PetscCall(SNESNewtonALComputeFunction(snes, X, Q));
    539:       /* R = F(X, lambda) */
    540:       PetscCall(SNESComputeFunction(snes, X, R));
    541:       PetscCall(VecNormBegin(R, NORM_2, &fnorm));
    542:       PetscCall(VecNormBegin(X, NORM_2, &xnorm));
    543:       PetscCall(VecNormBegin(deltaX, NORM_2, &ynorm));
    544:       PetscCall(VecNormEnd(R, NORM_2, &fnorm));
    545:       PetscCall(VecNormEnd(X, NORM_2, &xnorm));
    546:       PetscCall(VecNormEnd(deltaX, NORM_2, &ynorm));
    547:       SNESCheckFunctionDomainError(snes, fnorm);
    548:       if (PetscLogPrintInfo) PetscCall(SNESNewtonALCheckArcLength(snes, DeltaX, data->lambda_update, stepSize));

    550:       /* Monitor convergence */
    551:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    552:       snes->iter++;
    553:       snes->norm  = fnorm;
    554:       snes->ynorm = ynorm;
    555:       snes->xnorm = xnorm;
    556:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    557:       PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
    558:       PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
    559:       PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    560:       if (!snes->reason && j == maxits - 1) snes->reason = SNES_DIVERGED_MAX_IT;
    561:       if (snes->reason) break;
    562:     }
    563:     if (snes->reason < 0) break;
    564:     if (data->lambda >= data->lambda_max) {
    565:       break;
    566:     } else if (maxincs > 0 && i == maxincs - 1) {
    567:       snes->reason = SNES_DIVERGED_MAX_IT;
    568:       break;
    569:     } else {
    570:       snes->reason = SNES_CONVERGED_ITERATING;
    571:     }
    572:   }
    573:   /* Reset RHS vector, if changed */
    574:   if (data->copied_rhs) {
    575:     PetscCall(VecCopy(data->vec_rhs_orig, snes->vec_rhs));
    576:     data->copied_rhs = PETSC_FALSE;
    577:   }
    578:   snes->max_its = maxits; /* reset snes->max_its */
    579:   PetscFunctionReturn(PETSC_SUCCESS);
    580: }

    582: /*
    583:    SNESSetUp_NEWTONAL - Sets up the internal data structures for the later use
    584:    of the SNESNEWTONAL nonlinear solver.

    586:    Input Parameter:
    587: .  snes - the SNES context
    588: .  x - the solution vector

    590:    Application Interface Routine: SNESSetUp()
    591:  */
    592: static PetscErrorCode SNESSetUp_NEWTONAL(SNES snes)
    593: {
    594:   PetscFunctionBegin;
    595:   PetscCall(SNESSetWorkVecs(snes, 4));
    596:   PetscCall(SNESSetUpMatrices(snes));
    597:   PetscFunctionReturn(PETSC_SUCCESS);
    598: }

    600: /*
    601:    SNESSetFromOptions_NEWTONAL - Sets various parameters for the SNESNEWTONAL method.

    603:    Input Parameter:
    604: .  snes - the SNES context

    606:    Application Interface Routine: SNESSetFromOptions()
    607: */
    608: static PetscErrorCode SNESSetFromOptions_NEWTONAL(SNES snes, PetscOptionItems PetscOptionsObject)
    609: {
    610:   SNES_NEWTONAL             *data            = (SNES_NEWTONAL *)snes->data;
    611:   SNESNewtonALCorrectionType correction_type = data->correction_type;

    613:   PetscFunctionBegin;
    614:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES Newton Arc Length options");
    615:   PetscCall(PetscOptionsReal("-snes_newtonal_step_size", "Initial arc length increment step size", "SNESNewtonAL", data->step_size, &data->step_size, NULL));
    616:   PetscCall(PetscOptionsInt("-snes_newtonal_max_continuation_steps", "Maximum number of increment steps", "SNESNewtonAL", data->max_continuation_steps, &data->max_continuation_steps, NULL));
    617:   PetscCall(PetscOptionsReal("-snes_newtonal_psisq", "Regularization parameter for arc length continuation, 0 for cylindrical", "SNESNewtonAL", data->psisq, &data->psisq, NULL));
    618:   PetscCall(PetscOptionsReal("-snes_newtonal_lambda_min", "Minimum value of the load parameter lambda", "SNESNewtonAL", data->lambda_min, &data->lambda_min, NULL));
    619:   PetscCall(PetscOptionsReal("-snes_newtonal_lambda_max", "Maximum value of the load parameter lambda", "SNESNewtonAL", data->lambda_max, &data->lambda_max, NULL));
    620:   PetscCall(PetscOptionsBool("-snes_newtonal_scale_rhs", "Scale the constant vector passed to `SNESSolve` by the load parameter lambda", "SNESNewtonAL", data->scale_rhs, &data->scale_rhs, NULL));
    621:   PetscCall(PetscOptionsEnum("-snes_newtonal_correction_type", "Type of correction to use in the arc-length continuation method", "SNESNewtonALCorrectionType", SNESNewtonALCorrectionTypes, (PetscEnum)correction_type, (PetscEnum *)&correction_type, NULL));
    622:   PetscCall(SNESNewtonALSetCorrectionType(snes, correction_type));
    623:   PetscOptionsHeadEnd();
    624:   PetscFunctionReturn(PETSC_SUCCESS);
    625: }

    627: static PetscErrorCode SNESReset_NEWTONAL(SNES snes)
    628: {
    629:   SNES_NEWTONAL *al = (SNES_NEWTONAL *)snes->data;

    631:   PetscFunctionBegin;
    632:   PetscCall(VecDestroy(&al->vec_rhs_orig));
    633:   PetscCall(MatDestroy(&al->mat_diag_scaling));
    634:   PetscFunctionReturn(PETSC_SUCCESS);
    635: }

    637: /*
    638:    SNESDestroy_NEWTONAL - Destroys the private SNES_NEWTONAL context that was created
    639:    with SNESCreate_NEWTONAL().

    641:    Input Parameter:
    642: .  snes - the SNES context

    644:    Application Interface Routine: SNESDestroy()
    645:  */
    646: static PetscErrorCode SNESDestroy_NEWTONAL(SNES snes)
    647: {
    648:   PetscFunctionBegin;
    649:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", NULL));
    650:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", NULL));
    651:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", NULL));
    652:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", NULL));
    653:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", NULL));
    654:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", NULL));
    655:   PetscCall(PetscFree(snes->data));
    656:   PetscFunctionReturn(PETSC_SUCCESS);
    657: }

    659: /*MC
    660:   SNESNEWTONAL - Newton based nonlinear solver that uses a arc-length continuation method to solve the nonlinear system.

    662:   Options Database Keys:
    663: +   -snes_newtonal_step_size step                 - Initial arc length increment step size
    664: .   -snes_newtonal_max_continuation_steps max     - Maximum number of continuation steps, or negative for no limit (not recommended)
    665: .   -snes_newtonal_psisq psisq                    - Regularization parameter for arc length continuation, 0 for cylindrical. Larger values generally lead to more steps.
    666: .   -snes_newtonal_lambda_min lambda_min          - Minimum value of the load parameter lambda
    667: .   -snes_newtonal_lambda_max lambda_max          - Maximum value of the load parameter lambda
    668: .   -snes_newtonal_scale_rhs (true|false)         - Scale the constant vector passed to `SNESSolve` by the load parameter `lambda`
    669: -   -snes_newtonal_correction_type (exact|normal) - Type of correction to use in the arc-length continuation method

    671:   Level: intermediate

    673:   Note:
    674:   The exact correction scheme with partial updates is detailed in {cite}`Ritto-CorreaCamotim2008` and the implementation of the
    675:   normal correction scheme is based on {cite}`LeonPaulinoPereiraMenezesLages_2011`.

    677: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONAL`, `SNESNewtonALSetFunction()`, `SNESNewtonALGetFunction()`, `SNESNewtonALGetLoadParameter()`
    678: M*/
    679: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONAL(SNES snes)
    680: {
    681:   SNES_NEWTONAL *arclengthParameters;

    683:   PetscFunctionBegin;
    684:   snes->ops->setup          = SNESSetUp_NEWTONAL;
    685:   snes->ops->solve          = SNESSolve_NEWTONAL;
    686:   snes->ops->destroy        = SNESDestroy_NEWTONAL;
    687:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONAL;
    688:   snes->ops->reset          = SNESReset_NEWTONAL;

    690:   snes->usesksp = PETSC_TRUE;
    691:   snes->usesnpc = PETSC_FALSE;

    693:   PetscCall(SNESParametersInitialize(snes));
    694:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
    695:   PetscObjectParameterSetDefault(snes, max_its, 50);

    697:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    699:   PetscCall(PetscNew(&arclengthParameters));
    700:   arclengthParameters->lambda                 = 0.0;
    701:   arclengthParameters->lambda_update          = 0.0;
    702:   arclengthParameters->step_size              = 1.0;
    703:   arclengthParameters->max_continuation_steps = 100;
    704:   arclengthParameters->psisq                  = 1.0;
    705:   arclengthParameters->lambda_min             = 0.0;
    706:   arclengthParameters->lambda_max             = 1.0;
    707:   arclengthParameters->scale_rhs              = PETSC_TRUE;
    708:   arclengthParameters->correction_type        = SNES_NEWTONAL_CORRECTION_EXACT;
    709:   arclengthParameters->mat_diag_scaling       = NULL;
    710:   snes->data                                  = (void *)arclengthParameters;

    712:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetCorrectionType_C", SNESNewtonALSetCorrectionType_NEWTONAL));
    713:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetLoadParameter_C", SNESNewtonALGetLoadParameter_NEWTONAL));
    714:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetFunction_C", SNESNewtonALSetFunction_NEWTONAL));
    715:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALGetFunction_C", SNESNewtonALGetFunction_NEWTONAL));
    716:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALSetDiagonalScaling_C", SNESNewtonALSetDiagonalScaling_NEWTONAL));
    717:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonALComputeFunction_C", SNESNewtonALComputeFunction_NEWTONAL));
    718:   PetscFunctionReturn(PETSC_SUCCESS);
    719: }