"buffer source array is read-only" error with firedrake-adjoint

[ctjacobs]

With the default branch of my fork of dolfin-adjoint/firedrake-adjoint and up-to-date Firedrake, PyOP2, PETSc, FFC, etc, I have been running the following code:

from firedrake import *
op2.init(lazy_evaluation=False)
parameters['form_compiler']['quadrature_degree'] = 4
parameters["coffee"]["O2"] = False
from firedrake_adjoint import *
parameters["adjoint"]["test_derivative"] = True

mesh = UnitSquareMesh(10, 10)
U = VectorFunctionSpace(mesh, "CG", 2)
H = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([U, H])

def main(annotate=True):
    solution = Function(W, name="Solution", annotate=annotate)
    u, h = split(solution)
    w, v = TestFunctions(W)

    h_mean = Function(H).interpolate(Expression("50.0"))
    n = FacetNormal(mesh)
    F = -(h_mean+h)*inner(u, grad(v))*dx + (h_mean+h)*inner(u, n)*v*ds

    bcs = [DirichletBC(W.sub(0), Expression(("2", "0")), (3))]

    solve(F==0, solution, bcs=bcs, annotate=annotate)
    return solution

if __name__ == "__main__":

    solution = main(annotate=True)

    adj_html("forward.html", "forward")
    print "Replaying forward model"
    assert replay_dolfin(tol=1e-13, stop=True)

The first run of the forward model is successful, but replaying it with replay_dolfin results in the following error:

christian@elevate ~ $ python readonly_bug.py 
Replaying forward model

Python traceback: 
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/libadjoint/libadjoint.py", line 1294, in newfunc
    func(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/libadjoint/libadjoint.py", line 1647, in __mat_solve_callback__
    x = A.solve(Variable(var=adj_var), b)
  File "/data/dolfin-adjoint-test/dolfin_adjoint/adjlinalg.py", line 361, in solve
    x = self.basic_solve(var, b)
  File "/data/dolfin-adjoint-test/dolfin_adjoint/adjlinalg.py", line 296, in basic_solve
    compatibility.solve(F == 0, x.data, b.nonlinear_bcs, J=J, solver_parameters=self.solver_parameters)
  File "/data/firedrake/firedrake/solving.py", line 119, in solve
    _solve_varproblem(*args, **kwargs)
  File "/data/firedrake/firedrake/solving.py", line 154, in _solve_varproblem
    nullspace=nullspace)
  File "/data/firedrake/firedrake/variational_solver.py", line 251, in __init__
    ctx.set_jacobian(self.snes)
  File "/data/firedrake/firedrake/variational_solver.py", line 104, in set_jacobian
    snes.setJacobian(self.form_jacobian, J=self._jac._M.handle,
  File "/data/PyOP2/pyop2/petsc_base.py", line 438, in handle
    self._init()
  File "/data/PyOP2/pyop2/petsc_base.py", line 223, in _init
    self._init_nest()
  File "/data/PyOP2/pyop2/petsc_base.py", line 238, in _init_nest
    mat.createNest([[m.handle for m in row_] for row_ in self._blocks])
  File "/data/PyOP2/pyop2/petsc_base.py", line 438, in handle
    self._init()
  File "/data/PyOP2/pyop2/petsc_base.py", line 225, in _init
    self._init_block()
  File "/data/PyOP2/pyop2/petsc_base.py", line 297, in _init_block
    sparsity.fill_with_zeros(mat, self.sparsity.dims, self.sparsity.maps)
  File "pyop2/sparsity.pyx", line 266, in pyop2.sparsity.fill_with_zeros (pyop2/sparsity.cpp:4918)
    rmap = pair[0].values_with_halo
  File "stringsource", line 614, in View.MemoryView.memoryview_cwrapper (pyop2/sparsity.cpp:13431)
  File "stringsource", line 321, in View.MemoryView.memoryview.__cinit__ (pyop2/sparsity.cpp:9977)
ValueError: buffer source array is read-only

C traceback: 
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_call_unix64+0x4c)[0x7fba9ee48adc]
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_call+0x1fc)[0x7fba9ee4840c]
/usr/lib/python2.7/lib-dynload/_ctypes.x86_64-linux-gnu.so(_ctypes_callproc+0x48e)[0x7fba9f05f5fe]
/usr/lib/python2.7/lib-dynload/_ctypes.x86_64-linux-gnu.so(+0x15f9e)[0x7fba9f060f9e]
python(PyEval_EvalFrameEx+0x970)[0x52cc20]
python[0x56d0aa]
python(PyObject_CallObject+0x59)[0x55ab29]
/usr/lib/python2.7/lib-dynload/_ctypes.x86_64-linux-gnu.so(+0x13e96)[0x7fba9f05ee96]
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_closure_unix64_inner+0x1fb)[0x7fba9ee488cb]
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_closure_unix64+0x46)[0x7fba9ee48c44]
//usr/lib/libadjoint.so(adj_get_forward_solution+0x10a)[0x7fba926a4170]
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_call_unix64+0x4c)[0x7fba9ee48adc]
/usr/lib/x86_64-linux-gnu/libffi.so.6(ffi_call+0x1fc)[0x7fba9ee4840c]
/usr/lib/python2.7/lib-dynload/_ctypes.x86_64-linux-gnu.so(_ctypes_callproc+0x48e)[0x7fba9f05f5fe]
/usr/lib/python2.7/lib-dynload/_ctypes.x86_64-linux-gnu.so(+0x15f9e)[0x7fba9f060f9e]
python(PyEval_EvalFrameEx+0x970)[0x52cc20]
python(PyEval_EvalFrameEx+0xc82)[0x52cf32]
python(PyEval_EvalCodeEx+0x2a4)[0x55c594]
python(PyEval_EvalFrameEx+0x7dd)[0x52ca8d]
python(PyEval_EvalCodeEx+0x2a4)[0x55c594]
python(PyEval_EvalCode+0x32)[0x5b7392]
python[0x469663]
python(PyRun_FileExFlags+0x92)[0x4699e3]
python(PyRun_SimpleFileExFlags+0x2ee)[0x469f1c]
python(Py_Main+0xb5e)[0x46ab81]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7fbaa5b97ec5]
python[0x57497e]

I've noticed that if I remove the ds integral and/or the DirichletBC, replay_dolfin succeeds. Any thoughts as to what might be causing this? @funsim?

[wence-]

This error occurs if the array assigned to a memory view is read only. The line:

rmap = pair[0].values_with_halo

But the values_with_halo array is not read only in normal operation. Can you run in the debugger:

ipython foo.py --pdb

And when you get that error go up the stack to the fill_with_zeros call and inspect the read only status of those arrays?

for pair in self.sparsity.maps:
     print pair[0].values_with_halo.flags
     print pair[1].values_with_halo.flags

[ctjacobs]

I added

        for pair in self.sparsity.maps:
           print "Pair:", pair
           print pair[0].values_with_halo.flags
           print pair[1].values_with_halo.flags

just before sparsity.fill_with_zeros(mat, self.sparsity.dims, self.sparsity.maps) in petsc_base.py which gave:

Pair: (DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((121, 121, 121, 121), 'None_nodes'), 3, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Replaying forward model
Pair: (DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((200, 200, 200, 200), 'plexmesh_cells'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_cell_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
Pair: (DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])), DecoratedMap(Map(Set((40, 40, 40, 40), 'plexmesh_exterior_facets'), Set((441, 441, 441, 441), 'None_nodes'), 6, None, 'None_exterior_facet_node'), frozenset(['ALL']), frozenset([])))
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : False
  ALIGNED : True
  UPDATEIFCOPY : False
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : False
  ALIGNED : True
  UPDATEIFCOPY : False

[ctjacobs]

Is there a way in PyOP2 to force the exterior_facets pair to be writable? It's probably not a good thing to do, but it might be interesting to see what happens.

[wence-]

just before the line that causes the error in sparsity.pyx do:

from utils import maybe_setflags
maybe_setflags(bad_array, write=True)

[ctjacobs]

Thanks. That works, and the results look sensible.

[kynan]

Any clues yet as to what causes the issue?

[ctjacobs]

No idea yet I'm afraid.

[wence-]

Please try with OP2/PyOP2#444, I think I figured out what was going on. When we build the exterior_facet_boundary_node_map we use the exterior_facet_node_map as a Dat, and the data therefore gets set to read-only, the sparsity builder barfs because memory views require writeable buffers. My workaround is to temporarily set them writeable in the sparsity builder.

[ctjacobs]

Works for me - thanks.

[wence-]

Merged in PyOP2, so closing.