[petsc-dev] Inquiry Regarding Arc-Length Method Implementation in SNESNEWTONAL Solver

[Zach Atkins]

Hi David,

Thanks for reaching out! At a high-level, I would recommend reading through the docs and source code if you haven't already, as the tutorial does not go into much detail on the details of the arc length method implementation.

Docs: https://urldefense.us/v3/__https://petsc.org/main/manual/snes/*newton-with-arc-length-continuation__;Iw!!G_uCfscf7eWS!fkLQY6YXo3t6v6Di1MslMhnd9IaaagxBJ2eJStyMTIsHbXqM_XF9fAyqs0nggXxm9xDT2C63tge0xCpoFzV8wizyeLSBsA$ 
Source: https://urldefense.us/v3/__https://petsc.org/main/src/snes/impls/al/al.c.html__;!!G_uCfscf7eWS!fkLQY6YXo3t6v6Di1MslMhnd9IaaagxBJ2eJStyMTIsHbXqM_XF9fAyqs0nggXxm9xDT2C63tge0xCpoFzV8wiwoS3XDBQ$  (in particular, SNESSolve_NEWTONAL)

> ‎1. According to the tutorial, the arc-length constraint function `ComputeTangentLoad` is used to compute ( Q = dF^{ext}/d\lambda ), the derivative of the external load with respect to lambda. In the case of a linear load, i.e., ( F^{ext} = \lambda F^{ext} ), ( Q = F^{ext} ). In this scenario, is it correct to simply return ( F^{ext} ) inside `ComputeTangentLoad`?

Yes, this is correct. If you have a solution-independent load, you can also just pass it as the RHS vector to SNESSolve and it will automatically be used as a linear load in terms of lambda.

> 2. From my experiments, I observed that the `ComputeTangentLoad` function is called once before assembling the residual at every iteration. Interestingly, it seems that during the first iteration, ( Q ) contributes to the residual calculation, but from the second iteration onward, it no longer does. Could you help me understand why this happens?

We only include Q in the full residual in the predictor step, the corrector steps use a special newton iteration over lambda. If you are using a solution-dependent external force, i.e. the force is computed in the function passed to `SNESSetFunction`, then you will need to scale the force by the current value of lambda, which can be queried via `SNESNewtonALGetLoadParameter<https://urldefense.us/v3/__https://petsc.org/main/manualpages/SNES/SNESNewtonALGetLoadParameter/__;!!G_uCfscf7eWS!fkLQY6YXo3t6v6Di1MslMhnd9IaaagxBJ2eJStyMTIsHbXqM_XF9fAyqs0nggXxm9xDT2C63tge0xCpoFzV8wiwUW9aVcQ$ >()`.

> 3. I have attached my own buckling example of a Lee Frame structure where I implemented the arc-length method as the nonlinear solver. In my implementation, I divided the arc-length method into two steps:  
>    (1) A predictor step, which predicts the intersection of the constraint surface and the tangent;  
>    (2) A corrector step, which uses Newton-Raphson iterations to refine the intersection point between the constraint surface and the actual equilibrium path.  
>    This is the standard arc-length algorithm principle. I would like to know if `SNESNEWTONAL` follows a similar predictor-corrector procedure internally, or if I need to implement an outer load step loop and predictor step myself?

We also do a predictor/corrector loop, with two options for the corrector algorithm. The docs explain the selection of lambda in both cases in detail.

The arc length solver consists of an inner and outer solve, in an incremental-iterative style. Each increment is chosen via an inner iterative Newton solve which uses the predictor/corrector approach. The number of increments depends on the step size provided via -snes_newtonal_step_size, which is the 2-norm of the augmented increment vector [Delta x; Delta lambda].

> 4. I remember that the tutorial describes solving two linear systems during each corrector step to obtain ( \Delta x^F ) and ( \Delta x^Q ) separately. Where in the SNESNEWTONAL workflow are these two solves performed? Do they happen right after assembling the residual and Jacobian?
>
> 5. Additionally, is the Jacobian matrix used in SNESNEWTONAL the original system Jacobian, or is it augmented to include the constraint surface function? Similarly, is the input vector expanded to ( n+1 ) dimensions to account for ( \lambda ), or does it remain the original size ( n ) as in the standard Newton-Raphson implementation?

It is the original Jacobian — we avoid augmenting for efficiency and symmetry. The full process is described in the documentation as well as in the source code. Essentially, we factor the total SNES update into two parts, one which is determined by the the full residual at lambda, F(x, lambda), (this is \delta x^F ), and the other based on the variation in F with respect to lambda Q(x, lambda) (this is \delta x^Q). Both parts are computed once per inner Newton iteration right after the Jacobian is assembled. The residual and tangent load vectors are technically computed at the end of each inner Newton iteration, but that is mainly an implementation detail.

> 6. In general, I would like to understand more about the internal iteration details of `SNESSetType(snes, SNESNEWTONAL);` compared to `SNESSetType(snes, SNESNEWTONLS);` so that I can more accurately replicate and improve upon the arc-length algorithm I have implemented in the attached code using SNES.

I would recommend again referring to the documentation and source for easier comparison. If you have additional specific questions about the implementation, please feel free to ask!

Best,
Zach Atkins

________________________________
From: Jed Brown <jed at jedbrown.org>
Sent: Tuesday, March 18, 2025 11:34 AM
To: David Jiawei LUO LIANG <12431140 at mail.sustech.edu.cn>; petsc-dev <petsc-dev at mcs.anl.gov>
Cc: Zach Atkins <Zach.Atkins at colorado.edu>
Subject: Re: [petsc-dev] Inquiry Regarding Arc-Length Method Implementation in SNESNEWTONAL Solver

[jed at jedbrown.org appears similar to someone who previously sent you email, but may not be that person. Learn why this could be a risk at https://urldefense.us/v3/__https://aka.ms/LearnAboutSenderIdentification__;!!G_uCfscf7eWS!fkLQY6YXo3t6v6Di1MslMhnd9IaaagxBJ2eJStyMTIsHbXqM_XF9fAyqs0nggXxm9xDT2C63tge0xCpoFzV8wiwz13ACMQ$  ]

[External email - use caution]


Hi, David. I'm copying Zach, who developed this implementation and will be able to more directly answer your questions.

"David Jiawei LUO LIANG"        <12431140 at mail.sustech.edu.cn> writes:

> Dear PETSc Developers,
>
>
> I am currently using the SNESNEWTONAL nonlinear solver and its arc-length method module. I noticed from the tutorial that, besides the standard residual and Jacobian definitions:  
>
>
> SNESSetFunction(snes, r, FormFunction2nd, &user);  
> SNESSetJacobian(snes, J, J, FormJacobian2nd, &user);  
>
>
> it also requires:  
>
>
> SNESNewtonALSetFunction(snes, ComputeTangentLoad, &user);  
>
>
> which is used to define the arc-length constraint function.
>
>
> I have several questions regarding the arc-length implementation, and I would greatly appreciate your clarification:
>
>
> 1. According to the tutorial, the arc-length constraint function `ComputeTangentLoad` is used to compute ( Q = dF^{ext}/d\lambda ), the derivative of the external load with respect to lambda. In the case of a linear load, i.e., ( F^{ext} = \lambda F^{ext} ), ( Q = F^{ext} ). In this scenario, is it correct to simply return ( F^{ext} ) inside `ComputeTangentLoad`?
>
>
> 2. From my experiments, I observed that the `ComputeTangentLoad` function is called once before assembling the residual at every iteration. Interestingly, it seems that during the first iteration, ( Q ) contributes to the residual calculation, but from the second iteration onward, it no longer does. Could you help me understand why this happens?
>
>
> 3. I have attached my own buckling example of a Lee Frame structure where I implemented the arc-length method as the nonlinear solver. In my implementation, I divided the arc-length method into two steps:  
>    (1) A predictor step, which predicts the intersection of the constraint surface and the tangent;  
>    (2) A corrector step, which uses Newton-Raphson iterations to refine the intersection point between the constraint surface and the actual equilibrium path.  
>    This is the standard arc-length algorithm principle. I would like to know if `SNESNEWTONAL` follows a similar predictor-corrector procedure internally, or if I need to implement an outer load step loop and predictor step myself?
>
>
> 4. I remember that the tutorial describes solving two linear systems during each corrector step to obtain ( \Delta x^F ) and ( \Delta x^Q ) separately. Where in the SNESNEWTONAL workflow are these two solves performed? Do they happen right after assembling the residual and Jacobian?
>
>
> 5. Additionally, is the Jacobian matrix used in SNESNEWTONAL the original system Jacobian, or is it augmented to include the constraint surface function? Similarly, is the input vector expanded to ( n+1 ) dimensions to account for ( \lambda ), or does it remain the original size ( n ) as in the standard Newton-Raphson implementation?
>
>
> 6. In general, I would like to understand more about the internal iteration details of `SNESSetType(snes, SNESNEWTONAL);` compared to `SNESSetType(snes, SNESNEWTONLS);` so that I can more accurately replicate and improve upon the arc-length algorithm I have implemented in the attached code using SNES.
>
>
> Thank you very much for your time and assistance!
>
>
> Best regards,  
>
>
>
> David
>
>
>
>
>
> David Jiawei LUO LIANG
>
>
>
> 南方科技大学/学生/研究生/2024
>
>
>
> 广东省深圳市南山区学苑大道1088号
>
>
>
>
>  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20250318/bc3eed3d/attachment-0001.html>