La copie PETSc Matrix déclenche une exception: pourquoi?

Question

J'utilise petsc4py, et obtenir une exception que je ne comprends pas. Je définis la fonction suivante:

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() 

puis essayez ce qui suit. (pA a déjà été défini comme un PETSc.Mat 48x48 ailleurs, et assemblé.ksdg est une instance d'une classe sur laquelle je travaille, dont elle sera éventuellement une fonction membre, si je peux l'obtenir):

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

Ce lève l'exception suivante:

--------------------------------------------------------------------------- 
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 

Regarder petscerror.h.html, 63 est PETSC_ERR_ARG_OUTOFRANGE 63 /* input argument, out of range */.

Si quelqu'un comprend pourquoi PETSc ne me laisse pas copier les matrices A à B, j'apprécierais une explication. Merci.

Answer 1

OK, je l'ai compris. L'erreur 63 exception est levée lorsque l'on tente de créer un nouvel élément non nul dans une matrice creuse, si l'option

PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR 

est True. Ainsi, si j'exécute

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

immédiatement après la création de B, aucune exception est soulevée dans tsIJacobian. Ce que je ne comprends pas, c'est pourquoi je ne reçois pas la même exception pour J.