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

[Zou, Ling]

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/f9339966/attachment-0001.html>