<div dir="ltr">Thanks for your careful inspection and thoughtful suggestions.<br><br>>    finite differencing may put a small non-zero value in that location due to numerical round-off<br><br>I think your explanation is reasonable. This numerical round-off may somehow help to avoid this pivot issue.<br><br>The structure of my jacobian matrix looks like this (generated by '-mat_view draw'):<div><img src="cid:ii_kvbn0pq30" alt="jac_view.png" width="280" height="282"></div><div>Analytically, the first diagonal element of the jacobian is indeed 0, as its corresponding residual function is solely determined from boundary condition of another variable. This seems a little bit wired but is mathematically well-posed. For more description about the background physics, please refer to attached PDF file for more detailed explanation on the discretization and boundary conditions.</div><div><br></div><div>Actually, the jacobian matrix is not singular, but I do believe this numerical difficulty is caused by the zero-element in diagonal. </div><div>In this regard, I've performed some trial and test. It seems that several methods have been worked out for this pivot issue:</div><div>a) By setting '-pc_type svd', PETSC does not panic any more with my hand-coded jacobian, and converged solution is obtained. Efficiency is also preserved.</div><div>b) By setting '-pc_type none', converged solution is also obtained, but it takes too many KSP iterations to converge per SNES step<span class="gmail-Apple-converted-space"> </span>(usually hundreds), making the overall solution procedure very slow.</div><div><br></div><div>Do you think these methods really solved this kind of pivot issue? Not by chance like the numerical round-off in finite difference previously.</div><div><br></div><div>Regards</div><div><br></div><div>Yu Cang</div><div><br>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>> 于2021年10月27日周三 下午9:43写道：<br>><br>><br>>    You can run with -ksp_error_if_not_converged to get it to stop as soon as a linear solve fails to help track down the exact breaking point.<br>><br>> > The problem under consideration contains an eigen-value to be solved,<br>> > making the first diagonal element of the jacobian matrix being zero.<br>> > From these outputs, it seems that the PC failed to factorize, which is<br>> > caused by this 0 diagonal element.  But I'm wondering why it works<br>> > with jacobian matrix generated by finite-difference?<br>><br>>    Presumably your "exact" Jacobian puts a zero on the diagonal while the finite differencing may put a small non-zero value in that location due to numerical round-off. In that case even if the factorization succeeds it may produce an inaccurate solution if the value on the diagonal is very small.<br>><br>>    If your matrix is singular or cannot be factored with LU then you need to use a different solver for the linear system that will be robust to the zero on the diagonal. What is the structure of your Jacobian? (The analytic form).<br>><br>>   Barry<br>><br>><br>> > On Oct 27, 2021, at 1:47 AM, 仓宇 <<a href="mailto:yhcy1993@gmail.com">yhcy1993@gmail.com</a>> wrote:<br>> ><br>> > Thanks for your kind reply.<br>> ><br>> > Several comparison tests have been performed. Attached are execution<br>> > output files. Below are corresponding descriptions.<br>> ><br>> > good.txt -- Run without hand-coded jacobian, solution converged, with<br>> > option '-ts_monitor -snes_monitor -snes_converged_reason<br>> > -ksp_monitor_true_residual -ksp_converged_reason';<br>> > jac1.txt -- Run with hand-coded jacobian, does not converge, with<br>> > option '-ts_monitor -snes_monitor -snes_converged_reason<br>> > -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian';<br>> > jac2.txt -- Run with hand-coded jacobian, does not converge, with<br>> > option '-ts_monitor -snes_monitor -snes_converged_reason<br>> > -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian<br>> > -ksp_view';<br>> > jac3.txt -- Run with hand-coded jacobian, does not converge, with<br>> > option '-ts_monitor -snes_monitor -snes_converged_reason<br>> > -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian<br>> > -ksp_view -ts_max_snes_failures -1 ';<br>> ><br>> > The problem under consideration contains an eigen-value to be solved,<br>> > making the first diagonal element of the jacobian matrix being zero.<br>> > From these outputs, it seems that the PC failed to factorize, which is<br>> > caused by this 0 diagonal element.  But I'm wondering why it works<br>> > with jacobian matrix generated by finite-difference? Would employing<br>> > DMDA for discretization be helpful?<br>> ><br>> > Regards<br>> ><br>> > Yu Cang<br>> ><br>> > Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>> 于2021年10月25日周一 下午10:50写道：<br>> >><br>> >><br>> >>  It is definitely unexpected that -snes_test_jacobian verifies the Jacobian as matching but the solve process is completely different.<br>> >><br>> >>   Please run with -snes_monitor -snes_converged_reason -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian and send all the output<br>> >><br>> >>  Barry<br>> >><br>> >><br>> >>> On Oct 25, 2021, at 9:53 AM, 仓宇 <<a href="mailto:yhcy1993@gmail.com">yhcy1993@gmail.com</a>> wrote:<br>> >>><br>> >>> I'm using TS to solve a set of DAE, which originates from a<br>> >>> one-dimensional problem. The grid points are uniformly distributed.<br>> >>> For simplicity, the DMDA is not employed for discretization.<br>> >>><br>> >>> At first, only the residual function is prescribed through<br>> >>> 'TSSetIFunction', and PETSC produces converged results. However, after<br>> >>> providing hand-coded Jacobian through 'TSSetIJacobian', the internal<br>> >>> SNES object fails (residual norm does not change), and TS reports<br>> >>> 'DIVERGED_STEP_REJECTED'.<br>> >>><br>> >>> I have tried to add the option '-snes_test_jacobian' to see if the<br>> >>> hand-coded jacobian is somewhere wrong, but it shows '||J -<br>> >>> Jfd||_F/||J||_F = 1.07488e-10, ||J - Jfd||_F = 2.14458e-07',<br>> >>> indicating that the hand-coded jacobian is correct.<br>> >>><br>> >>> Then, I added a monitor for the internal SNES object through<br>> >>> 'SNESMonitorSet', in which the solution vector will be displayed at<br>> >>> each iteration. It is interesting to find that, if the jacobian is not<br>> >>> provided, meaning finite-difference is utilized for jacobian<br>> >>> evaluation internally, the solution vector converges to steady<br>> >>> solution and the SNES residual norm is reduced continuously. However,<br>> >>> it turns out that, as long as the jacobian is provided, the solution<br>> >>> vector will NEVER get changed! So the solution procedure stucked!<br>> >>><br>> >>> This is quite strange!  Hope to get some advice.<br>> >>> PETSC version=3.14.6, program run in serial mode.<br>> >>><br>> >>> Regards<br>> >>><br>> >>> Yu Cang<br>> >><br>> > <jac3.txt><jac2.txt><jac1.txt><good.txt><br>></div></div>