[petsc-users] Fwd: Some guidance on setting an arc-length solver up in PETSc with TS/SNES

[Francesc Levrero-Florencio]

Dear PETSc team and users,

I am trying to implement a “non-consistent arc-length method” (i.e.
non-consistent as in the Jacobian from a traditional load-controlled method
is used instead of the “augmented one”, the latter would need an extra/row
column for the constraint terms; the non-consistent method requires two
linear solves, with different RHS, at every nonlinear iteration).

I am trying to do this in the context of TS/SNES. Computation of the
Jacobian remains the same, but I am thinking of doing most of the work
within the function to compute the residual. The rough non-consistent
arc-length algorithm is (delta here refers to the correction per nonlinear
iteration):

- Solve K * delta u^k_1 = - R

- Solve K * delta u^k_2 = - F_ext

- Find the correct root of delta lambda

- lambda^k = lambda^{k-1} + delta lambda^k

- delta u^k = delta u^k_1 + delta lambda^k * delta u^k_2

- u^{k} = u^{k-1} + delta u^k, u here is the solution to the mechanical
problem (displacement).


The rest is the same as a load-controlled scheme (for a more complete
arc-length algorithm, please see “A simple extrapolated predictor for
overcoming the starting and tracking issues in the arc-length method for
nonlinear structural mechanics”, Kadapa, C. (2021), *Eng Struct*). I see
two potential issues with an arc-length implementation by using TS/SNES:

- Needing to modify the solution vector within each nonlinear iteration. I
believe this could be sorted out by allowing TS/SNES to obtain the solution
to the traditional load-controlled problem (u_1 above), and only keep track
of the correction per nonlinear iteration, which would correspond to delta
u^k_1 = u^k_1 (given by the TS at the current iteration) – u^{k-1}_1 (given
by the TS at the previous iteration). I imagine the nonlinear iterative
correction would be correct because the residual (and jacobian) would be
calculated by using an internally stored u (u above, stored outside the TS,
which is not the same as the TS internally stored u_1). I believe this
would be okay in quasistatic simulations since solution derivatives are not
required. How could this be amended for transient simulations? Is there any
other simpler way to achieve this? Can the solution vector be set to a
predictor value at the beginning of a time step (maybe prestep or poststep
related functions)?

- Solving two linear systems per iteration. I imagine the way around this
is to do TSGetKSP(ts, ksp), and then KSPSolve(ksp, - F_ext, delta u^k_2)
within the TS, since the TS would be using the correctly calculated
Jacobian K above. Would this an effective way to achieve this?

Thanks for your help in advance and please keep up the good work!


Regards,

Francesc.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210716/a6a257e4/attachment.html>