I am trying to conduct a simulation with the rate and state model but I have convergence problems. I found in the v2.2.2 manual (section 6.4.5.1) that “The solution of these two linear systems gives the increment in slip assuming all the degrees of freedom except those imme-
diately adjacent to the fault remain fixed. In real applications where the deformation associated with fault slip is localized
around the fault, this provides good enough approximations so that the nonlinear solver converges quickly. In problems where
deformation associated with slip on the fault is not localized (as in the case in some of the example problems), the increment in
slip computed by solving these two linear systems is not a good approximation and the nonlinear solve requires a large number
of iterations.”
I do not know whether this could be the issue (deformation is not localized on the fault surface in my case). I have checked the examples and YouTube tutorials and applied the suggested settings:
pc_type = asm
sub_pc_factor_shift_type = nonzero
ksp_rtol = 1.0e-30
ksp_atol = 1.0e-12
snes_rtol = 1.0e-30
snes_atol = 1.0e-10
zero_tolerance = 1.0e-11
zero_tolerance_normal = 1.0e-11
friction.linear_slip_rate = 1.0e-9
friction_pc_type = asm

friction_sub_pc_factor_shift_type = nonzero

friction_ksp_max_it = 100

friction_ksp_gmres_restart = 30

friction_ksp_error_if_not_converged = true

friction_ksp_monitor = true

#friction_ksp_view = true

friction_ksp_converged_reason = true

The linear solver takes thousands of iterations every time step, and eventually the non linear solver fails to converge. The model has about 30k elements.
Are there any specific settings for cases where deformation is not localized in the fault plane?

It looks like you are using the asm preconditioner. It is known not to work well for problems with faults. As discussed in PyLith v2.2.2 manual section 4.1 and the 2013 Aagard, Knepley, and Williams JGR paper, the custom fault preconditioner with algebriac multigrid preconditioning on the displacement block performs much better. See Table 4.4 with the corresponding settings provided in share/settings/solver_fault_fieldsplit.cfg.

Thanks for your reply. I’ve tried those settings before and I also get convergence problems (at an early time step). Could you please explain what does fs_* stand for?. Also, is it OK to use friction_* settings combined with fs_* settings ? When I set the *_error_if_not_converged & * _monitor to true, the friction solves converge within a few iterations, but SNES solver does not converge.
Is there a way to display which node has the highest residual?

PyLith v2.2.2 with friction uses two solvers. The “outer” solve is for the entire domain, while the “inner” solve is for the fault friction. The friction_ solve settings apply only to the inner solve.

When trying to improve the performance, the first thing to do is to get the outer linear solve to converge reasonably fast. This should happen with the field split solver and custom fault preconditioner that I mentioned. The inner (friction) solve should always converge quickly with the ASM preconditioner.

Diagnosing failure of the nonlinear solve to converge can be tricky. You should turn on the SNES monitor to snes_monitor = true to see what the residual is doing. Is it decreasing but doing so very slowly, does it approach a constant value, or does it increase to a large value?

Thanks for your answer. I’ve turned on the SNES monitor and the residual is very large the first time (0 SNES Function norm 4.877127797937e+05), then it gets stuck at about 5 e-7 until it eventually fails. I wonder why it fails at a number of iterations smaller than snes_max_it (message is Nonlinear solve did not converge due to DIVERGED_FUNCTION_COUNT iterations 258).
I’ve set the cohesion to 100 MPa to make sure there is no slip on the fault, and the problem persists, so it seems the problem is outside the fault, but I don’t know where.

If there is no slip on the fault and the bulk rheology is linear elastic, then the solution should converge in 1 SNES iteration. I would start by checking the mesh quality. Make sure there are no cells with condition numbers greater than about 2-5. I would then omit the fault from the simulation settings (you should change your preconditioner settings to simply use pc_type = ml and not use the field split and custom fault preconditioner settings). If the problem works without the fault, then put it back in and run a problem with a prescribed slip of zero (use the “fault” preconditioning settings). If prescribed slip works, then use a very large compressive normal traction like 1.0e+12 Pa and make sure the fault doesn’t slip and that the solution converges in one SNES iteration. If that works, slowly adjust the parameters to get to the desired ones. Hopefully, this will lead you to identify what is going wrong.

I was wondering how much tight should the absolute tolerances of ksp and snes be when using rate and state. I know there are constraints between zero_tolerance, snes_atol and ksp_atol (and also the linear slip rate), but can we just shift all the values?
Also, the value of the reference slip rate V0 can significantly affect the value of the computed friction coefficient. How do you estimate it?

Yes, you can just shift/scale the KSP absolute tolerance, SNES absolute tolerance, and zero_tolerance values. Note that these are all nondimensional values. That is, they apply directly to the internal values solver values that have been nondimensionalized.

The reference slip rate V0 in the rate-state friction model should be set according to the desired friction model values. The value should be chosen in accordance with the coefficient of friction and the characteristic slip distance.

Thanks for your answer. In choosing the reference slip rate, should one also consider the expected values of the slip rate V (and thus the time step size), given the log dependence?
As suggested, I checked the condition number of the mesh, it is between 1 and about 2.5. When I try the field split settings that you suggested, the simulation is extremely slow and it does not converge in a reasonable time (the KSP within SNES errates, the norms decrease and then increase and so on, globally decreasing but after thousands of iterations). The simulation runs with ASM. The conceptual model is basically a box with rollers on x=0, y=0 and zbottom, and Neuman boundary conditions on the other three boundaries. Gravity is on, and there is a quasi-vertical fault in the middle.

If the KSP (linear) solver does not converge in a reasonable number of iterations (usually less than 100) with the field split settings, then something is wrong in the model setup. Is this true for prescribed slip?

This suggests there is something wrong or unusual in the model setup. If you provide all of the input files for prescribed slip, we can take a look as see if we can find the problem. Along with the input files, please include a diagram / sketch showing the problem you are trying to solve.