assembly of form with keep_diagonal fails.

[Evan Cummings]

Hello,

I'm receiving an error when assembling a form defined over facets when using the keep_diagonal=True option of assemble. The simple script below assembles the facet normal vectors over a 3D mesh:

from fenics import *

k    = 2  # works
k    = 4  # fails
mesh = UnitCubeMesh(k,k,k)
V    = VectorFunctionSpace(mesh, "CG", 1)
n    = FacetNormal(mesh)

trial = TrialFunction(V)
test  = TestFunction(V)

a = inner(trial, test)*ds
L = inner(n,     test)*ds

A = assemble(a, keep_diagonal=True)
A.ident_zeros()
b = assemble(L)

n = Function(V)
solve(A, n.vector(), b, 'cg', 'amg')

File('n.pvd') << n

This fails at A = assemble(a, keep_diagonal=True) with the message

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [e-mail]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValues'.
*** Reason:  PETSc error code is: 63.
*** Where:   This error was encountered inside /build/dolfin-i8W3lr/dolfin-2016.1.0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2016.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

[Evan Cummings]

minor error in code.

[Jan Blechta]

Trace from debugger looks really funny:

(gdb) f 4
#4  0x00007fffeeef5c1b in dolfin::AssemblerBase::init_global_tensor (this=this@entry=0x1bd0ee0, A=..., a=...)
    at ../../dolfin/fem/AssemblerBase.cpp:123
123         _matA.set(&block, 1, &_i, 1, &_i);
(gdb) p _i
$4 = 138
(gdb) p &_i
$5 = (dolfin::la_index *) 0x7fffffffc010
(gdb) f 3
#3  0x00007fffeeffd004 in dolfin::PETScMatrix::set (this=0x1967130, block=0x7fffffffc090, m=1, rows=0x7fffedcf3d88, n=1, cols=0x0)
    at ../../dolfin/la/PETScMatrix.cpp:233

The value of _i is fine (space dimension is 375). But its address &_i and actual parameters on frame 3 rows and cols do not match.

Could be a bug in GCC.

[Jan Blechta]

Could somebody else try to reproduce it?

1. Save the MWE above to test0.py
2. Run gdb -ex "set breakpoint pending on" -ex "b PetscError" -ex r -ex "f 4" -ex "p &_i" -ex "f 3" -ex "p rows" -ex "p cols" -args python test0.py | grep ^\\$
3. Wait for output like

Function "PetscError" not defined.
$1 = (dolfin::la_index *) 0x7fffffffc010
$2 = (const dolfin::la_index *) 0x7fffedcf3d88
$3 = (const dolfin::la_index *) 0x0

and type q <return> y <return> to quit gdb.

If adresses $1-$3 are not all same there is something wrong.

[Tormod Landet]

I get the same "PETSc error code is: 63" when running in just python and then the following from gdb+python:

$ gdb -ex "set breakpoint pending on" -ex "b PetscError" -ex r -ex "f 4" -ex "p &_i" -ex "f 3" -ex "p rows" -ex "p cols" -args python test0.py | grep ^\\$
Function "PetscError" not defined.
Detaching after fork from child process 2277.
Detaching after fork from child process 2310.
$1 = (dolfin::la_index *) 0x7fffffffc4c0
$2 = (const dolfin::la_index *) 0x7fffffffc4c0
$3 = (const dolfin::la_index *) 0x7fffffffc4c0

[Tormod Landet]

That is with gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-4) and a full reinstall run of fenics-install-all.sh from yesterday

[Jan Blechta]

Thanks but I came to a belief that the screwed up pointers are probably just a code generation artifact. In fact, the issue will be a bit less difficult to debug once pull request #289 is merged.

[Jan Blechta]

Another failing example by @hernanmella

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 2)

u0 = interpolate(Expression(("x[0]","x[1]"), degree=2), V)
u = TestFunction(V)
v = TrialFunction(V)

a = inner(u, v)*ds
L = inner(u0, v)*ds

A = assemble(a, keep_diagonal=True)
L = assemble(L)
A.ident_zeros()

[Jan Blechta]

Suggested fix in https://bitbucket.org/fenics-project/dolfin/branch/jan/fix-issue-721.

[Jan Blechta]

@garth-wells, do you think that the quick fix 3365a1c is now ok?

I regard the fix rather as a workaround because SparsityPattern interface does not know anything about blocks and PETScMatrix::init assumes that passed sparsity pattern is divisible by blocks. If merge this quick fix I'd file an issue describing the principle issue.

[Prof Garth Wells]

I'm not really fussed since I really dislike this option. A proper fix would be:

1. An option to add a diagonal entry to a sparsity pattern; and
2. An option asking the assembler to not finalise the matrix assembly

[Jan Blechta]

Temporarily fixed by 3365a1ccf50067b9df5a0e05138c361049935893. The real issue filed at #802.