[petsc-dev] Proposed changes to TS API

[Smith, Barry F.]

   Here is MY summary of the discussion so far.

1) the IFunction/IJacobian interface has its supporters. There is valid argument that for certain cases it can be more efficient than the proposed alternative; but this seems largely a theoretical believe at this time since there are no comparisons between nontrivial (or even trivial) codes that use the assembly directly the mass shift plus Jacobian vs the assembly separately and MatAXPY the two parts together.  As far as I can see this potential performance is the ONLY benefit to keeping the IFunction/IJacobian() interface? Please list additional benefits if there are any?

2) The IFunction/IJacobian approach makes it impossible for TS to access the mass matrix alone. 

3) But one can access the IJacobian matrix alone by passing a shift of zero

4) TSComputeIJacobian() is an ugly piece of shit code that has a confusing name since it also can incorporates the RHS Jacobian. 

4a) the TSComputeIJacobian() has issues for linear problems because it overwrites the user provided Jacobian and has hacks to deal with it

5) There is no standard way to solve M u = F() explicitly from the IFunction() form, (and cannot even with expressed with the explicit RHS approach) the user must manage the M solve themselves inside their RHSFunction.

6) There is some support for adding two new function callbacks that a) compute the mass matrix and b) compute the Jacobian part of the IJacobian. This appears particularly useful for implementing adjoints.

7) If we added the two new interfaces the internals of an already overly complicated TS become even more convoluted and unmaintainable.  For kicks I list the current TSComputeIJacobian() which I consider unmaintainable already.

  if (ijacobian) {
    PetscBool missing;
    PetscStackPush("TS user implicit Jacobian");
    ierr = (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);CHKERRQ(ierr);
    PetscStackPop;
    ierr = MatMissingDiagonal(A,&missing,NULL);CHKERRQ(ierr);
    if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Amat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
    if (B != A) {
      ierr = MatMissingDiagonal(B,&missing,NULL);CHKERRQ(ierr);
      if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Bmat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
    }
  }
  if (imex) {
    if (!ijacobian) {  /* system was written as Udot = G(t,U) */
      PetscBool assembled;
      ierr = MatZeroEntries(A);CHKERRQ(ierr);
      ierr = MatAssembled(A,&assembled);CHKERRQ(ierr);
      if (!assembled) {
        ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
        ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
      }
      ierr = MatShift(A,shift);CHKERRQ(ierr);
      if (A != B) {
        ierr = MatZeroEntries(B);CHKERRQ(ierr);
        ierr = MatAssembled(B,&assembled);CHKERRQ(ierr);
        if (!assembled) {
          ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
          ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
        }
        ierr = MatShift(B,shift);CHKERRQ(ierr);
      }
    }
  } else {
    Mat Arhs = NULL,Brhs = NULL;
    if (rhsjacobian) {
      ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr);
      ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr);
    }
    if (Arhs == A) {           /* No IJacobian, so we only have the RHS matrix */
      PetscBool flg;
      ts->rhsjacobian.scale = -1;
      ts->rhsjacobian.shift = shift;
      ierr = SNESGetUseMatrixFree(ts->snes,NULL,&flg);CHKERRQ(ierr);
      /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
      if (!flg) {
        ierr = MatScale(A,-1);CHKERRQ(ierr);
        ierr = MatShift(A,shift);CHKERRQ(ierr);
      }
      if (A != B) {
        ierr = MatScale(B,-1);CHKERRQ(ierr);
        ierr = MatShift(B,shift);CHKERRQ(ierr);
      }
    } 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);
      if (A != B) {
        ierr = MatAXPY(B,-1,Brhs,axpy);CHKERRQ(ierr);
      }
    }
  }

  Please comment and continue discussion.


> On May 11, 2018, at 5:09 PM, Jed Brown <jed at jedbrown.org> wrote:
> 
> "Smith, Barry F." <bsmith at mcs.anl.gov> writes:
> 
>>>> The current IJacobian is essentially SNESJacobian. And the single-matrix SNESJacobian interface is always there. Power users could set up the SNESJacobian directly if we pass a read-only shift parameter to them. So we are by no means prohibiting power users from doing what they want.
>>> 
>>> You mean the user would call TSGetSNES and SNESSetJacobian, then
>>> internally call TSGetIJacobianShift and some new function to create
>>> Udot?  That sounds way messier and error-prone.
>>> 
>>> And completely impossible to compose. We should be explicitly talking about subsystems that use TS. For example,
>>> I have a nonlinear system for plasticity. I want to use a preconditioner that introduces some elasticity and timesteps to
>>> steady state to provide a good Newton direction. I need to to be able to create the solver without all sorts of digging
>>> in and setting things. This idea that you "just set SNESJacobian" is a non-starter.
>> 
>>   But this is exactly what TSComputeIJacobian currently does, so is the current interface a non-starter?
> 
> It isn't at all the current interface.  If the user is calling
> SNESSetJacobian, then we would need to open up the bowels of every
> SNESTSFormJacobian_* and make stable public APIs out of those internals
> (including the wiring for nonlinear multigrid).  This sounds like a
> terrible thing to make public.