<div dir="ltr">I'm solving a "linear elasticity with contact" problem using SNES. The matrix is symmetric, and it corresponds to the standard linear elasticity matrix, plus some extra nonlinear "penalty spring" terms due to the contact conditions.<div><br></div><div>The solve works fine if I use a direct solver (e.g. MUMPS). I've also found that CG usually works fine too, but in one specific case it fails to converge. (I thought that penalty springs + iterative solvers may not be a good combination, but it does seem to work fine in most cases.)<div><br></div><div>Below I've pasted the (last part of the) output of "-ksp_type cg -ksp_norm_type natural -snes_view -ksp_monitor" for the case that fails.</div><div><br></div><div>If anyone has some suggestions about other solver/preconditioner options that would be worth trying, that would be most appreciated.</div><div><br></div><div>Thanks,</div><div>David</div><div><br></div><div>-------------------------------------------------------------------------------</div><div><br></div><div><div>  102 KSP Residual norm 8.910585826247e-05 </div><div>  103 KSP Residual norm 8.395178827803e-05 </div><div>  104 KSP Residual norm 7.758478472168e-05 </div><div>  105 KSP Residual norm 7.099534889984e-05 </div><div>  106 KSP Residual norm 6.615351528474e-05 </div><div>  107 KSP Residual norm 6.303646443688e-05 </div><div>  NL step  5, |residual|_2 = 1.996487e+01</div><div>    0 KSP Residual norm 5.844456247091e+00 </div><div>    1 KSP Residual norm 1.308756943973e+00 </div><div>    2 KSP Residual norm 7.545133187463e-01 </div><div>    3 KSP Residual norm 1.346068942607e+00 </div><div>SNES Object: 1 MPI processes</div><div>  type: newtonls</div><div>  maximum iterations=50, maximum function evaluations=10000</div><div>  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08</div><div>  total number of linear solver iterations=535</div><div>  total number of function evaluations=11</div><div>  norm schedule ALWAYS</div><div>  SNESLineSearch Object:   1 MPI processes</div><div>    type: bt</div><div>      interpolation: cubic</div><div>      alpha=1.000000e-04</div><div>    maxstep=1.000000e+08, minlambda=1.000000e-12</div><div>    tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08</div><div>    maximum iterations=40</div><div>  KSP Object:   1 MPI processes</div><div>    type: cg</div><div>    maximum iterations=10000, initial guess is zero</div><div>    tolerances:  relative=1e-05, absolute=1e-50, divergence=10000</div><div>    left preconditioning</div><div>    using NATURAL norm type for convergence test</div><div>  PC Object:   1 MPI processes</div><div>    type: ilu</div><div>      ILU: out-of-place factorization</div><div>      0 levels of fill</div><div>      tolerance for zero pivot 2.22045e-14</div><div>      matrix ordering: natural</div><div>      factor fill ratio given 1, needed 1</div><div>        Factored matrix follows:</div><div>          Mat Object:           1 MPI processes</div><div>            type: seqaij</div><div>            rows=967608, cols=967608</div><div>            package used to perform factorization: petsc</div><div>            total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07</div><div>            total number of mallocs used during MatSetValues calls =0</div><div>              using I-node routines: found 322536 nodes, limit used is 5</div><div>    linear system matrix = precond matrix:</div><div>    Mat Object:    ()     1 MPI processes</div><div>      type: seqaij</div><div>      rows=967608, cols=967608</div><div>      total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07</div><div>      total number of mallocs used during MatSetValues calls =0</div><div>        using I-node routines: found 322536 nodes, limit used is 5</div><div>Number of nonlinear iterations: 5</div><div>Nonlinear solver convergence/divergence reason: DIVERGED_LINEAR_SOLVE</div></div></div></div>