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.