[petsc-users] Strange behavior of TS after setting hand-coded Jacobian

[Zhang, Hong]

You can avoid the zero-diagonal elements in your case. Just remove the boundary points and the \hat{\Lambda}_r from the solution vector Y (eq(17) in your PDE file). You can calculate them separately in your RHS function (when computing the residual vector for 2 <=i<=N-1). It is not necessary to form a DAE. You can actually form an ODE instead.

Hong (Mr.)

> On Oct 31, 2021, at 5:38 AM, 仓宇 <yhcy1993 at gmail.com> wrote:
> 
> Thanks for your suggestions and experience sharing.
> 
> But it seems that you have missed my point. This kind of zero-diagonal
> element that I described earlier is somehow inevitable.
> If interested in this problem, you may refer to the PDF file that was
> attached in the previous post.
> 
> Regards
> 
> Yu
> 
> Zhang, Hong <hongzhang at anl.gov> 于2021年10月29日周五 下午10:05写道：
>> 
>> One way to avoid the zero element in Jacobian is to exclude the boundary point from the solution vector. I often do this for Dirichlet boundary conditions since the value at the boundary is given directly and does not need to be taken as a degree of freedom.
>> 
>> Hong (Mr.)
>> 
>> On Oct 28, 2021, at 9:49 PM, 仓宇 <yhcy1993 at gmail.com> wrote:
>> 
>> Thanks for your careful inspection and thoughtful suggestions.
>> 
>>>   finite differencing may put a small non-zero value in that location due to numerical round-off
>> 
>> I think your explanation is reasonable. This numerical round-off may somehow help to avoid this pivot issue.
>> 
>> The structure of my jacobian matrix looks like this (generated by '-mat_view draw'):
>> <jac_view.png>
>> 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.
>> 
>> Actually, the jacobian matrix is not singular, but I do believe this numerical difficulty is caused by the zero-element in diagonal.
>> In this regard, I've performed some trial and test. It seems that several methods have been worked out for this pivot issue:
>> 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.
>> b) By setting '-pc_type none', converged solution is also obtained, but it takes too many KSP iterations to converge per SNES step (usually hundreds), making the overall solution procedure very slow.
>> 
>> Do you think these methods really solved this kind of pivot issue? Not by chance like the numerical round-off in finite difference previously.
>> 
>> Regards
>> 
>> Yu Cang
>> 
>> Barry Smith <bsmith at petsc.dev> 于2021年10月27日周三 下午9:43写道：
>>> 
>>> 
>>>   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.
>>> 
>>>> The problem under consideration contains an eigen-value to be solved,
>>>> making the first diagonal element of the jacobian matrix being zero.
>>>> From these outputs, it seems that the PC failed to factorize, which is
>>>> caused by this 0 diagonal element.  But I'm wondering why it works
>>>> with jacobian matrix generated by finite-difference?
>>> 
>>>   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.
>>> 
>>>   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).
>>> 
>>>  Barry
>>> 
>>> 
>>>> On Oct 27, 2021, at 1:47 AM, 仓宇 <yhcy1993 at gmail.com> wrote:
>>>> 
>>>> Thanks for your kind reply.
>>>> 
>>>> Several comparison tests have been performed. Attached are execution
>>>> output files. Below are corresponding descriptions.
>>>> 
>>>> good.txt -- Run without hand-coded jacobian, solution converged, with
>>>> option '-ts_monitor -snes_monitor -snes_converged_reason
>>>> -ksp_monitor_true_residual -ksp_converged_reason';
>>>> jac1.txt -- Run with hand-coded jacobian, does not converge, with
>>>> option '-ts_monitor -snes_monitor -snes_converged_reason
>>>> -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian';
>>>> jac2.txt -- Run with hand-coded jacobian, does not converge, with
>>>> option '-ts_monitor -snes_monitor -snes_converged_reason
>>>> -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian
>>>> -ksp_view';
>>>> jac3.txt -- Run with hand-coded jacobian, does not converge, with
>>>> option '-ts_monitor -snes_monitor -snes_converged_reason
>>>> -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian
>>>> -ksp_view -ts_max_snes_failures -1 ';
>>>> 
>>>> The problem under consideration contains an eigen-value to be solved,
>>>> making the first diagonal element of the jacobian matrix being zero.
>>>> From these outputs, it seems that the PC failed to factorize, which is
>>>> caused by this 0 diagonal element.  But I'm wondering why it works
>>>> with jacobian matrix generated by finite-difference? Would employing
>>>> DMDA for discretization be helpful?
>>>> 
>>>> Regards
>>>> 
>>>> Yu Cang
>>>> 
>>>> Barry Smith <bsmith at petsc.dev> 于2021年10月25日周一 下午10:50写道：
>>>>> 
>>>>> 
>>>>> It is definitely unexpected that -snes_test_jacobian verifies the Jacobian as matching but the solve process is completely different.
>>>>> 
>>>>>  Please run with -snes_monitor -snes_converged_reason -ksp_monitor_true_residual -ksp_converged_reason -snes_test_jacobian and send all the output
>>>>> 
>>>>> Barry
>>>>> 
>>>>> 
>>>>>> On Oct 25, 2021, at 9:53 AM, 仓宇 <yhcy1993 at gmail.com> wrote:
>>>>>> 
>>>>>> I'm using TS to solve a set of DAE, which originates from a
>>>>>> one-dimensional problem. The grid points are uniformly distributed.
>>>>>> For simplicity, the DMDA is not employed for discretization.
>>>>>> 
>>>>>> At first, only the residual function is prescribed through
>>>>>> 'TSSetIFunction', and PETSC produces converged results. However, after
>>>>>> providing hand-coded Jacobian through 'TSSetIJacobian', the internal
>>>>>> SNES object fails (residual norm does not change), and TS reports
>>>>>> 'DIVERGED_STEP_REJECTED'.
>>>>>> 
>>>>>> I have tried to add the option '-snes_test_jacobian' to see if the
>>>>>> hand-coded jacobian is somewhere wrong, but it shows '||J -
>>>>>> Jfd||_F/||J||_F = 1.07488e-10, ||J - Jfd||_F = 2.14458e-07',
>>>>>> indicating that the hand-coded jacobian is correct.
>>>>>> 
>>>>>> Then, I added a monitor for the internal SNES object through
>>>>>> 'SNESMonitorSet', in which the solution vector will be displayed at
>>>>>> each iteration. It is interesting to find that, if the jacobian is not
>>>>>> provided, meaning finite-difference is utilized for jacobian
>>>>>> evaluation internally, the solution vector converges to steady
>>>>>> solution and the SNES residual norm is reduced continuously. However,
>>>>>> it turns out that, as long as the jacobian is provided, the solution
>>>>>> vector will NEVER get changed! So the solution procedure stucked!
>>>>>> 
>>>>>> This is quite strange!  Hope to get some advice.
>>>>>> PETSC version=3.14.6, program run in serial mode.
>>>>>> 
>>>>>> Regards
>>>>>> 
>>>>>> Yu Cang
>>>>> 
>>>> <jac3.txt><jac2.txt><jac1.txt><good.txt>
>>> 
>> <TFM.pdf>
>> 
>>