[petsc-dev] TS confusion

[Barry Smith]

   I am ONLY providing RHSJacobian and RHSFunction and not IJacobian at all. I am using ARKIMEX but have called 
  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);
  ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr);

   You can run src/ts/examples/tutorials/advection-diffusion-reaction/ex10.c and put a breakpoint in it and you will see it is calling the MatAXPY() as I indicated. I suspect the logic of fully implicit TSARKIMEX is far from optimized when only RHSJacobian is provided.

    What TSType should I use?

   Barry

On Sep 15, 2013, at 10:01 AM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> Barry Smith <bsmith at mcs.anl.gov> writes:
> 
>> If no ijacobian is provided and not imex TS still seems to produce two matrices 
>> 
>>    Mat Arhs = NULL,Brhs = NULL;
>>    MatStructure flg2;
>>    if (rhsjacobian) {
>>      ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr);
>>      ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr);
>>    }
>> 
>> while TSComputeRHSJacobian
>> 
>>  ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr);
>>  if (Arhs) {
>>    if (!ts->Arhs) {
>>      ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr);
>>    }
>>    *Arhs = ts->Arhs;
>>  }
> 
> Where is this code?  I don't see it in TSComputeRHSJacobian.
> 
>>    } else if (Arhs) {          /* Both IJacobian and RHSJacobian */
>>      MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
>>      if (!ijacobian) {         /* No IJacobian provided, but we have a separate RHS matrix */
>>        ierr = MatZeroEntries(*A);CHKERRQ(ierr);
>>        ierr = MatShift(*A,shift);CHKERRQ(ierr);
>>        if (*A != *B) {
>>          ierr = MatZeroEntries(*B);CHKERRQ(ierr);
>>          ierr = MatShift(*B,shift);CHKERRQ(ierr);
>>        }
>>      }
>>      ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr);
> 
> Do we agree that if the user provides both IJacobian and RHSJacobian,
> and if we are using non-imex assembly, then we need both matrices and
> need the MatAXPY?
> 
>> The end result is two sets of matrices. A and Arhs and then MatAXPY costs me
>> 
>> MatAXPY               47 1.0 7.1141e+00 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 22  0  0  0  0  22  0  0  0  0     0
>> 
>> when it doesn't even need to be done.  The code needs some refactoring but I don't understand enough about its convoluted logic.
>> What is going on here. Why does it feel it needs to create the non RHS Jacobian matrix?
> 
> I suspect the logic just wasn't optimized.  Can you think of a way to
> simplify the logic?