[petsc-users] finite difference jacobian errors when given non-constant initial condition

Mark Lohry mlohry at gmail.com
Sun Apr 21 12:34:54 CDT 2024


The coloring I'm fairly confident is correct -- I use the same process for
3D unstructured grids and everything seems to work. The residual function
is also validated.

As a test I did as you suggested -- assume the matrix is dense -- and I get
the same bad results, just now the zero blocks are filled.

Assuming dense, giving it a constant vector, all is good:

 4.23516e-16     -1.10266      0.31831   -0.0852909            0
 0     -0.31831      1.18795
     1.10266 -4.23516e-16     -1.18795      0.31831            0
 0    0.0852909     -0.31831
    -0.31831      1.18795  2.11758e-16     -1.10266      0.31831
-0.0852909            0            0
   0.0852909     -0.31831      1.10266 -4.23516e-16     -1.18795
 0.31831            0            0
           0            0     -0.31831      1.18795  2.11758e-16
-1.10266      0.31831   -0.0852909
           0            0    0.0852909     -0.31831      1.10266
-4.23516e-16     -1.18795      0.31831
     0.31831   -0.0852909            0            0     -0.31831
 1.18795  4.23516e-16     -1.10266


Assuming dense, giving it sin(x), all is bad:

-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08
 6.07287e+07  6.07287e+07  1.76177e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08
 4.52116e+07  4.52116e+07  1.31161e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08
-4.52116e+07 -4.52116e+07 -1.31161e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08
-6.07287e+07 -6.07287e+07 -1.76177e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08
-6.07287e+07 -6.07287e+07 -1.76177e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08
-4.52116e+07 -4.52116e+07 -1.31161e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08
 4.52116e+07  4.52116e+07  1.31161e+08
-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08
 6.07287e+07  6.07287e+07  1.76177e+08

Scratching my head over here... I've been using these routines successfully
for years in much more complex code.

On Sun, Apr 21, 2024 at 12:36 PM Zou, Ling <lzou at anl.gov> wrote:

> Edit:
>
>    - how do you do the coloring when using PETSc finite differencing? An
>    incorrect coloring may give you wrong Jacobian. For debugging purpose,
>    the simplest way to avoid an incorrect coloring is to assume the matrix is
>    dense (slow but error proofing). If the numeric converges as expected,
>    then fine tune your coloring to make it right and fast.
>
>
>
>
>
> *From: *petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf of Zou,
> Ling via petsc-users <petsc-users at mcs.anl.gov>
> *Date: *Sunday, April 21, 2024 at 11:29 AM
> *To: *Mark Lohry <mlohry at gmail.com>, PETSc <petsc-users at mcs.anl.gov>
> *Subject: *Re: [petsc-users] finite difference jacobian errors when given
> non-constant initial condition
>
> Hi Mark, I am working on a project having similar numeric you have,
> one-dimensional finite volume method with second-order slope limiter TVD,
> and PETSc finite differencing gives perfect Jacobian even for complex
> problems.
>
> So, I tend to believe that your implementation may have some problem. Some
> lessons I learned during my code development:
>
>
>
>    - how do you do the coloring when using PETSc finite differencing? An
>    incorrect coloring may give you wrong Jacobian. The simplest way to avoid
>    an incorrect coloring is to assume the matrix is dense (slow but error
>    proofing).
>    - Residual function evaluation not correctly implemented can also lead
>    to incorrect Jacobian. In your case, you may want to take a careful look at
>    the order of execution, when to update your unknown vector, when to perform
>    P1 reconstruction, and when to evaluate the residual.
>
>
>
> -Ling
>
>
>
> *From: *petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf of Mark
> Lohry <mlohry at gmail.com>
> *Date: *Saturday, April 20, 2024 at 1:35 PM
> *To: *PETSc <petsc-users at mcs.anl.gov>
> *Subject: *[petsc-users] finite difference jacobian errors when given
> non-constant initial condition
>
> I have a 1-dimensional P1 discontinuous Galerkin discretization of the
> linear advection equation with 4 cells and periodic boundaries on
> [-pi,+pi]. I'm comparing the results from SNESComputeJacobian with a
> hand-written Jacobian. Being linear,
>
> ZjQcmQRYFpfptBannerStart
>
> *This Message Is From an External Sender *
>
> This message came from outside your organization.
>
>
>
> ZjQcmQRYFpfptBannerEnd
>
> I have a 1-dimensional P1 discontinuous Galerkin discretization of the
> linear advection equation with 4 cells and periodic boundaries on
> [-pi,+pi]. I'm comparing the results from SNESComputeJacobian with a
> hand-written Jacobian. Being linear, the Jacobian should be
> constant/independent of the solution.
>
>
>
> When I set the initial condition passed to SNESComputeJacobian as some
> constant, say f(x)=1 or 0, the petsc finite difference jacobian agrees with
> my hand coded-version. But when I pass it some non-constant value, e.g.
> f(x)=sin(x), something goes horribly wrong in the petsc jacobian.
> Implementing my own rudimentary finite difference approximation (similar to
> how I thought petsc computes it) it returns the correct jacobian to
> expected error. Any idea what could be going on?
>
>
>
> Analytically computed Jacobian:
>
>  4.44089e-16     -1.10266      0.31831   -0.0852909            0
>  0     -0.31831      1.18795
>      1.10266 -4.44089e-16     -1.18795      0.31831            0
>  0    0.0852909     -0.31831
>     -0.31831      1.18795  4.44089e-16     -1.10266      0.31831
> -0.0852909            0            0
>    0.0852909     -0.31831      1.10266 -4.44089e-16     -1.18795
>  0.31831            0            0
>            0            0     -0.31831      1.18795  4.44089e-16
> -1.10266      0.31831   -0.0852909
>            0            0    0.0852909     -0.31831      1.10266
> -4.44089e-16     -1.18795      0.31831
>      0.31831   -0.0852909            0            0     -0.31831
>  1.18795  4.44089e-16     -1.10266
>     -1.18795      0.31831            0            0    0.0852909
> -0.31831      1.10266 -4.44089e-16
>
>
>
>
>
> petsc finite difference jacobian when given f(x)=1:
>
>  4.44089e-16     -1.10266      0.31831   -0.0852909            0
>  0     -0.31831      1.18795
>      1.10266 -4.44089e-16     -1.18795      0.31831            0
>  0    0.0852909     -0.31831
>     -0.31831      1.18795  4.44089e-16     -1.10266      0.31831
> -0.0852909            0            0
>    0.0852909     -0.31831      1.10266 -4.44089e-16     -1.18795
>  0.31831            0            0
>            0            0     -0.31831      1.18795  4.44089e-16
> -1.10266      0.31831   -0.0852909
>            0            0    0.0852909     -0.31831      1.10266
> -4.44089e-16     -1.18795      0.31831
>      0.31831   -0.0852909            0            0     -0.31831
>  1.18795  4.44089e-16     -1.10266
>     -1.18795      0.31831            0            0    0.0852909
> -0.31831      1.10266 -4.44089e-16
>
>
>
> petsc finite difference jacobian when given f(x) = sin(x):
>
> -1.65547e+08 -3.31856e+08 -1.25427e+09   4.4844e+08            0
>  0  1.03206e+08  7.86375e+07
>  9.13788e+07  1.83178e+08  6.92336e+08  -2.4753e+08            0
>  0 -5.69678e+07 -4.34064e+07
>   3.7084e+07  7.43387e+07  2.80969e+08 -1.00455e+08  -5.0384e+07
> -2.99747e+07            0            0
>   3.7084e+07  7.43387e+07  2.80969e+08 -1.00455e+08  -5.0384e+07
> -2.99747e+07            0            0
>            0            0  2.80969e+08 -1.00455e+08  -5.0384e+07
> -2.99747e+07 -2.31191e+07 -1.76155e+07
>            0            0  2.80969e+08 -1.00455e+08  -5.0384e+07
> -2.99747e+07 -2.31191e+07 -1.76155e+07
>  9.13788e+07  1.83178e+08            0            0 -1.24151e+08
> -7.38608e+07 -5.69678e+07 -4.34064e+07
> -1.65547e+08 -3.31856e+08            0            0  2.24919e+08
> 1.3381e+08  1.03206e+08  7.86375e+07
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240421/0dc502b3/attachment-0001.html>


More information about the petsc-users mailing list