- edited description
assembly of form with keep_diagonal fails.
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
***
*** fenics-support@googlegroups.com
***
*** 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
*** -------------------------------------------------------------------------
Comments (13)
-
reporter -
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 3rows
andcols
do not match.Could be a bug in GCC.
-
Could somebody else try to reproduce it?
- Save the MWE above to
test0.py
- 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 ^\\$
- 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. - Save the MWE above to
-
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
-
That is with
gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-4)
and a full reinstall run offenics-install-all.sh
from yesterday -
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.
-
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()
-
-
- changed milestone to 2017.1
- changed component to linear algebra
-
assigned issue to
- marked as minor
-
@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 andPETScMatrix::init
assumes that passed sparsity pattern is divisible by blocks. If merge this quick fix I'd file an issue describing the principle issue. -
I'm not really fussed since I really dislike this option. A proper fix would be:
- An option to add a diagonal entry to a sparsity pattern; and
- An option asking the assembler to not finalise the matrix assembly
-
Temporarily fixed by 3365a1ccf50067b9df5a0e05138c361049935893. The real issue filed at #802.
-
- changed status to resolved
- Log in to comment
minor error in code.