[petsc-users] Error when using TSSetRHSJacobian+DMSetMatrixPreallocateOnly

[Vijay S. Mahadevan]

> DMPlex does it.

I was just thinking about how you managed this in DMPlex. I will take
a look at this implementation.

Vijay

On Wed, Sep 18, 2013 at 12:10 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Wed, Sep 18, 2013 at 10:06 AM, Vijay S. Mahadevan <vijay.m at gmail.com>
> wrote:
>>
>> >    Until you can provide the preallocation  you need to simply call
>> > MatSetUp() and then MatAssemblyBegin/End() on your matrix inside your
>> > DMCreateMatrix_Moab(). The expectation is that the DM matrix provided is
>> > "ready to go", we didn't expect to handle matrices that were "incomplete" in
>> > other parts of the code.
>>
>> Interesting. Calling MatAssemblyBegin/End() without any MatSetValues*
>> wipes out the preallocation information. With -info turned on, I get
>> the following for some sample run.
>>
>> [0] MatAssemblyEnd_SeqAIJ(): Matrix size: 5192 X 5192; storage space:
>> 6906304 unneeded,0 used
>> [0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
>> [0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 0
>>
>> I do understand the expectation that matrix from DMCreateMatrix needs
>> to be complete but just that traversing through the mesh and setting
>> the right NNZ is difficult (I dont think this should exist in the
>> library code since this is physics/problem dependent). Somewhat
>> manageable for DMDA since it relies on a predefined stencil but I
>> couldn't wrap my head around doing this for a generic enough
>> implementation. I will think about this some more.
>
>
> DMPlex does it.
>
>    Matt
>
>>
>> > The "user matrix" must have been assembled otherwise it would have
>> > stopped above.
>>
>> It did.
>>
>> Vijay
>>
>> On Wed, Sep 18, 2013 at 11:17 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> >
>> > On Sep 18, 2013, at 10:52 AM, "Vijay S. Mahadevan" <vijay.m at gmail.com>
>> > wrote:
>> >
>> >>>    Why are you setting this flag?  If you don't set it then the code
>> >>> will go through.
>> >>
>> >> Yes it passes through with DMDA. But with my custom DM implementation
>> >> (MOAB based), this fails unless I explicitly set the NNZ pattern in
>> >> DMCreateMatrix. I currently don't do that and the only way to
>> >> replicate this behavior with DMDA was setting the flag explicitly.
>> >
>> >    Until you can provide the preallocation  you need to simply call
>> > MatSetUp() and then MatAssemblyBegin/End() on your matrix inside your
>> > DMCreateMatrix_Moab(). The expectation is that the DM matrix provided is
>> > "ready to go", we didn't expect to handle matrices that were "incomplete" in
>> > other parts of the code.
>> >
>> >>
>> >>>   There is no reason to use a user created matrix since the DM can
>> >>> provide the correct one. One of the "main jobs" of DM is to create that
>> >>> matrix so we really don't intend people to use a DM but then create the
>> >>> matrix themselves.
>> >>
>> >> Agreed, and that's what I expected when I passed in NULL args instead
>> >> of user matrix. But the DM created matrix errors out during Duplicate
>> >> while user provided matrix passes through, even though both of them
>> >> are unassembled.
>> >
>> > PetscErrorCode  MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
>> > {
>> >   PetscErrorCode ierr;
>> >   Mat            B;
>> >   PetscInt       i;
>> >
>> >   PetscFunctionBegin;
>> >   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
>> >   PetscValidType(mat,1);
>> >   PetscValidPointer(M,3);
>> >   if (!mat->assembled)
>> > SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for
>> > unassembled matrix");
>> >
>> > The "user matrix" must have been assembled otherwise it would have
>> > stopped above.
>> >
>> >
>> >> And this is what is confusing.
>> >>
>> >> Vijay
>> >>
>> >> On Wed, Sep 18, 2013 at 10:46 AM, Barry Smith <bsmith at mcs.anl.gov>
>> >> wrote:
>> >>>
>> >>>   The code did not "crash", it ended with an error message
>> >>>
>> >>>
>> >>> On Sep 18, 2013, at 10:32 AM, "Vijay S. Mahadevan" <vijay.m at gmail.com>
>> >>> wrote:
>> >>>
>> >>>> All,
>> >>>>
>> >>>> When I call DMSetMatrixPreallocateOnly
>> >>>
>> >>>    Why are you setting this flag?  If you don't set it then the code
>> >>> will go through.
>> >>>
>> >>>> and then TSSetRHSJacobian
>> >>>> without a matrix input argument (NULL instead), then the code crashes
>> >>>> during TSSolve and more specifically during MatDuplicate. This is
>> >>>> replicable with src/ts/examples/tutorials/ex2.c. Here's the error
>> >>>> with
>> >>>
>> >>>   Jed is reworking the logic here so the code will change.
>> >>>
>> >>>
>> >>>> stack trace.
>> >>>>
>> >>>> [0]PETSC ERROR: --------------------- Error Message
>> >>>>
>> >>>>
>> >>>> I hope this is not supposed to happen. The fix seems to be to call
>> >>>> TSSetRHSJacobian with a user created matrix instead of letting DM
>> >>>> create one,
>> >>>
>> >>>   There is no reason to use a user created matrix since the DM can
>> >>> provide the correct one. One of the "main jobs" of DM is to create that
>> >>> matrix so we really don't intend people to use a DM but then create the
>> >>> matrix themselves.
>> >>>
>> >>>   Barry
>> >>>
>> >>>
>> >>>> which rectifies this issue. This behavior feels
>> >>>> inconsistent and is there a way to fix this.
>> >>>>
>> >>>> Vijay
>> >>>>
>> >>>> PS: Just replace the following lines in ex2.c to replicate the error.
>> >>>>
>> >>>> /*  ierr =
>> >>>> TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);CHKERRQ(ierr); */
>> >>>> ierr =
>> >>>> TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
>> >>>> ierr =
>> >>>> DMSetMatrixPreallocateOnly(appctx.da,PETSC_TRUE);CHKERRQ(ierr);
>> >>>
>> >
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener