DirichletBC has not a proper way checking that tensors are compatible
Description: More careful compatibility checks in DirichletBC will be introduced by 670c1f385c27d5ce64e9123114baa33f945761a4 when solving #717. But the problem is that in combination with #719 user codes are prone to subtle bugs:
# Create tensors on one place
def create_data(A=None):
"This function creates data used in the tests below"
mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
if A is None:
A = assemble(u*v*dx)
else:
assemble(u*v*dx, tensor=A)
b = assemble(v*dx)
x = Vector()
return A, x, b
# Create bc on other place from different dofmap!!!
def test_p31_box_2():
mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, "CG", 1)
A, x, b = create_data()
u_0 = Expression("sin(x[0])", degree=1)
bc = DirichletBC(V, u_0, "x[0] > 0.5 && on_boundary")
u = Function(V)
bc.apply(A, b)
bc.apply(u.vector())
Although incorrect, this will worl fine in sequential and will produce random strange errors in parallel due to #719. In fact, this is wrong code from test/unit/python/book/test_chapter10.py
which was causing random failures for a long time and was not understood until #717 was being resolved.
Possible resolution: One possibility would be that tensors would carry std::size_t GenericTensor::_index_map_id
which would be compared to function_space->dofmap()->index_map()->id()
in DirichletBC
.
Comments (7)
-
-
reporter I agree. But the fact that we haven't been able to debug this for years demonstrates how DirichletBC accepting a wrong input is dangerous.
-
reporter I have of course fixed the wrong code in the test.
-
Since the latest SCOTCH behaves as expected and the old version of SCOTCH isn't causing a serious problem, I wouldn't bother with any workaround.
-
reporter Do you mean a workaround for
#719? This is a feature (rather than a bug) which still might be worth resolving. It helps users not to write an incorrect code. -
I mean do not bother implementing any work-arounds for SCOTCH 6.0.3. It's more code, and we'd need to remember to remove it later.
We've tended to not add work-arounds for bugs in old version of libraries because it clutters the code and isn't great for maintenance.
-
reporter Ok, I agree.
- Log in to comment
This is a 'bug' in the test. If we need to be sure that two meshes are partitioned in the same way, the mesh object should be re-used.