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

[仓宇]

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'):
[image: 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>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211029/5e79cd18/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: jac_view.png
Type: image/png
Size: 1998 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211029/5e79cd18/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: TFM.pdf
Type: application/pdf
Size: 44074 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211029/5e79cd18/attachment-0001.pdf>