[petsc-users] Different precision from MatAssembly/MatView

[Noam T.]

Hello,

Looking at the assembly of the Jacobian using MatView, I noticed values were somewhat different from the computed ones in the jacobian function (DMSNESSetJacobian). For the test, I used a single Q1 element, with 4 quadrature points, so the exact Jacobian matrix can be computed analytically.

In the jacobian function, after looping through all quadrature points, printing the matrix to stdout shows:

----
0.531360000000 0.066670000000 0.197333333327
0.066670000000 0.535360000000 -0.003333333327
0.197333333327 -0.003333333327 0.527413333333
----

(showing only the 3 free DoFs). The 2x2 submatrix in [0:1] has exactly four/five nonzero digits, as shown above. The rest of the elements have a periodic 3 ending (a couple of digits at the end are off, but that's fine).

This same matrix is the one added to the global jacobian during the quadrature loop, with "DMPlexMatSetClosure". After looping, calling "MatAssemblyBegin/End", then MatView:

----

Mat Object: 1 MPI process
type: seqaij
row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338)
row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385)
row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744) ----

Values are close, but definitely not the same as computed. Are these values from "MatView" the ones actually stored in the matrix? It seems there is some precision loss in the process.

Interestingly, computing the jacobian through finite differences, shows the same result for the "hand-coded" jacobian, whereas the finite differences ones is the "exact" one (same as first one above):

----
||J - Jfd||_F/||J||_F = 4.90736e-05, ||J - Jfd||_F = 4.74265e-05
Hand-coded Jacobian ----------
Mat Object: 1 MPI process
type: seqaij
row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338)
row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385)
row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744)
Finite difference Jacobian ----------
Mat Object: 1 MPI process
type: seqaij
row 0: (0, 0.53136) (1, 0.06667) (2, 0.197333)
row 1: (0, 0.06667) (1, 0.53536) (2, -0.00333333)
row 2: (0, 0.197333) (1, -0.00333333) (2, 0.527413)
Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ----------
Mat Object: 1 MPI process
type: seqaij
row 0: (0, 2.66554e-05)
row 1: (1, 2.66554e-05)
row 2: (2, 2.6558e-05)  ----

The main code is in Fortran, using double precision variables, in case it is relevant.

Thank you.

Noam
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260315/421ee049/attachment.html>