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

[Jed Brown]

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号
>
>
>
>
>