[petsc-users] finite difference jacobian errors when given non-constant initial condition
Mark Lohry
mlohry at gmail.com
Fri Jul 12 15:09:59 CDT 2024
The root cause of this turned out to be that I took code that i had
historically used for JFNK solves via TS and SNESSolve interfaces. Here
instead of calling SNESSolve I only wanted to compute the jacobian using
what i had set up as to be computed via finite differences via coloring.
when calling SNESSolve i confirmed it does generate the expected matrix,
but if instead i directly call
SNESComputeJacobian(ctx.snes_, petscsoln, ctx.JPre_, ctx.JPre_); // this is
computing without coloring
This appears to produce a matrix that is erroneous in some spots due to not
using coloring. Is there some other way to get the jacobian only? i tried
to follow similar steps as what is in the actual SNESSolve but i seem to
have botched it somehow.
On Sun, Apr 21, 2024 at 3:36 PM Zou, Ling <lzou at anl.gov> wrote:
> The other symptom is the same:
>
> - Using coloring, finite differencing respects the specified non-zero
> pattern, but gives wrong (very large) Jacobian entries (J_ij)
> - Using dense matrix assumption, finite differencing does not respect
> the non-zero pattern determined by your numeric, which is a clear sign of
> residual function code bug (your residual function does not respect your
> numeric).
>
> -Ling
>
>
>
> *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 2:28 PM
> *To: *Mark Lohry <mlohry at gmail.com>
> *Cc: *PETSc <petsc-users at mcs.anl.gov>
> *Subject: *Re: [petsc-users] finite difference jacobian errors when given
> non-constant initial condition
>
> Very interesting.
>
> I happened to encounter something very similar a couple of days ago,
> which, of course, was due to a code bug I introduced. The code bug was in
> the residual function. I used a local vector to track ‘heat flux’, which
> should be zero-ed out at the beginning of each residual function
> evaluation. I did not zero it, and I observed very similar results, the
> Jacobian is completely wrong, with large values (J_ij keeps increasing
> after each iteration), and non-zero values are observed in locations which
> should be perfect zero. The symptom is very much like what you are seeing
> here. I suspect a similar bug. (Maybe you forgot zero the coefficients of
> P1 re-construction? Using constant value 1, reconstructed dphi/dx = 0, so
> however many iterations, still zero).
>
>
>
> -Ling
>
>
>
> *From: *Mark Lohry <mlohry at gmail.com>
> *Date: *Sunday, April 21, 2024 at 12:35 PM
> *To: *Zou, Ling <lzou at anl.gov>
> *Cc: *PETSc <petsc-users at mcs.anl.gov>
> *Subject: *Re: [petsc-users] finite difference jacobian errors when given
> non-constant initial condition
>
> 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
>
> ZjQcmQRYFpfptBannerStart
>
> *This Message Is From an External Sender *
>
> This message came from outside your organization.
>
>
>
> ZjQcmQRYFpfptBannerEnd
>
> 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/20240712/58d4e25d/attachment-0001.html>
More information about the petsc-users
mailing list