Копия PETSc Matrix вызывает исключение: почему?

Question

Я использую petsc4py и получаю исключение, которое не понимаю. Я определяю следующую функцию:

def tsIJacobian(self, ts, t, u, udot, shift, A, B):
    self.setup_problem()
    psol = fe.as_backend_type(self.sol.vector()).vec()
    pA = fe.as_backend_type(self.A).mat()
    u.copy(psol)
    JU = fe.assemble(fe.derivative(ksdg.U_terms, ksdg.U))
    Jrho = fe.assemble(fe.derivative(ksdg.rho_terms, ksdg.rho))
    pJU = fe.as_backend_type(JU).mat()
    pJrho = fe.as_backend_type(Jrho).mat()
    pA.copy(A)
    A.scale(shift)
    A.axpy(1.0, pJU)
    A.axpy(1.0, pJrho)
    A.assemble()
    if not (A is B):
        A.copy(B)
        B.assemble()

тогда попробуйте следующее. (pA уже был определен как 48x48 PETSc.Mat в другом месте и собран. ksdg - это экземпляр класса, над которым я работаю, и со временем это будет функция-член, если я смогу заставить его работать):

J = pA.duplicate()
B = pA.duplicate()
tsIJacobian(ksdg, ts, 0, psol, pdsol, 0.1, J, B)

Это вызывает следующее исключение:

---------------------------------------------------------------------------
Error                                     Traceback (most recent call last)
<ipython-input-24-2a71f7ccf0af> in <module>()
----> 1 tsIJacobian(ksdg, ts, 0, psol, pdsol, 0.1, J, B)

<ipython-input-22-14579c08d6ae> in tsIJacobian(self, ts, t, u, udot, shift, A, B)
     15     A.assemble()
     16     if not (A is B):
---> 17         A.copy(B)
     18         B.assemble()

PETSc/Mat.pyx in petsc4py.PETSc.Mat.copy (src/petsc4py.PETSc.c:118071)()

Error: error code 63

Смотря на petscerror.h.html63 - это PETSC_ERR_ARG_OUTOFRANGE 63 /* input argument, out of range */,

Если кто-то понимает, почему PETSc не позволяет мне копировать матрицу A в B, я был бы признателен за объяснение. Благодарю.

Answer 1

Хорошо, я вроде понял это. Исключение ошибки 63 возникает при попытке создать новый ненулевой элемент в разреженной матрице, если опция

PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR

является True, Таким образом, если я выполню

B.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)

сразу после создания BВ якобском языке не исключение. Я не понимаю, почему я не получаю то же исключение для J,