We have a long-standing design flaw with GenericTensor::init (and
GenericVector::resize) that has raised its head with the fix to
https://bitbucket.org/fenics-project/dolfin/issue/132. The flaw is that we
allow init() to be called more than once. GenericFoo objects store a reference
to the underlying backend implementation, e.g. a PETSc Mat. Other objects, e.g.
a PETSc KSP, will often store a reference to the matrix or vector object too.
However, if init() is called a second time, the reference to the underlying
object is released and a new object is created. This leads to a mis-match
between objects. In the case of the fix to the SNES interface bug (Issue 132),
the SNES object has a reference to the operator A and the preconditioner matrix
P. In most case, these are the same. However, because init() is called by
default in the assembler, the two SNES matrix pointers no longer point to the
same object. The result is that the solver fails. To use the fixed SNES
interface, a user must now tell the assembler to not reset the Jacobian.
Allowing GenericFoo::init() to be called more than once also misleads the user
and does not help the user write good code. Since init() destroys the
underlying object and creates a new one, we should just have the user create a
new object if this is required (which is very rare). This way it’s clear what’s
happening. GenericFoo::init was conceived to do nothing if resizing was not
required, but it is not reasonably possible to perform the necessary checks as
the size, parallel layout and sparsity pattern all need to be checked. The
first is easy to check efficiently, the latter two are not. Hence, the
assemblers call init() by default in DOLFIN.
I propose that we permit a matrix/vector to be sized/initialised only once.
This is:
- consistent with the backends (PETSc and Epetra);
- does not obscure what happens behind-the-scene;
- avoids errors and inefficiencies when other objects hold references to the
underlying data; and
- will help users write faster code (considerably faster in many cases) by
re-using matrix data structures rather than frequent behind-the-scenes
destruction and rebuilding.
If init() is called twice, we can throw an error. During a transition, this can
be a warning. The assemblers will not need the ‘reset_tensor’ option because
the assembler will only initialise an empty matrix.
Garth
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
