Scalar::apply() is wrong

[Jan Blechta]

Following MWE shows the bug

from dolfin import *
mesh = UnitSquareMesh(3, 3)
S = Scalar()
assemble(Constant(1.0)*dx(mesh), tensor=S, add_values=True)
assemble(Constant(1.0)*dx(mesh), tensor=S, add_values=True)
print S.getval() # 3 when run by two processes

which happens because of

    /// Finalize assembly of tensor
    void apply(std::string mode)
    { _value = MPI::sum(_mpi_comm, _value); }

in dolfin/la/Scalar.cpp.

[Martin Sandve Alnæs]

I have a fix in martinal/fix-scalar-apply introducing a dual local/global value of Scalar.

However I wonder if the behaviour of all the functions are correct. Scalar is probably undertested.

In particular, should operator=(double) in C++ and s.assign(value) in python assign the global or local value? I've implemented it now to set the local value for the tiny existing tests to pass. However this means that we get the behaviour:

double d = 3.0;
Scalar s;
s = d;
// at this point (double)s != d
s.apply("insert");
// at this point (double)s == d

which I'd strongly argue is not expected behaviour of operator=.

Rather more confusing would be code like this:

double d = 3.0;
Scalar s;
s = d;
d = s;
// at this point d == 0.0!

Can we just remove operator= and operator double()? That means .assign in python as well.

[Jan Blechta]

I think that Scalar::insert*, Scalar::operator= and Scalar::apply('insert') are not well-defined operations and could be implemented by { dolfin_not_implemented(); }. It's not obvious what they should do when you assign different values across ranks. This should be sufficient as assemblers need just add_local and apply('add').

I'd also suggest that Scalar::_local_value is renamed to something like _local_increment or something similar.

[Prof Garth Wells]

I agree with @blechta that the operations are not well defined in parallel, and that an informative error message could be added.

Also, maybe we can remove the string argument to GenericTensor::apply. It was added to handle the peculiar behaviour of Epetra. We no longer support Epetra.

[Martin Sandve Alnæs]

Good!

[Jan Blechta]

@garth-wells There's still apply('flush') which we officially support - not that anybody is using it much. BTW, you've forgotten implementing it for Scalar,@martinal. And now I see that PETScVector does not implement it - only PETScMatrix.

[Jan Blechta]

VecAssemblyBegin does not support flush in a comparison to MatAssemblyBegin. So it would be probably better to remove it from GenericTensor.

[Martin Sandve Alnæs]

Fix in next.

[Jan Blechta]

Fixed by 39177f2.

[Prof Garth Wells]

Removing milestone: 1.5 (automated comment)