[petsc-dev] Picard with TS

[Barry Smith]

> On Jan 2, 2017, at 4:05 PM, Mark Adams <mfadams at lbl.gov> wrote:
> 
> 
> 
> On Sun, Jan 1, 2017 at 1:00 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>    What are LandSNESFunction() and LandSNESJacobian()?
> 
> These are the same methods that I give to TS (IFunction and IJacobian),  except there is no mass matrix in IJacobian. What do I do here? I do not have a shift.
> 
> But I see that IFunction is not correct. I am solving du/dt - C(u) = 0.

    Do you mean du/dt - c(u)* u = 0?

>  My IFunction applies C(u)*u, but LandSNESFunction should just return zero, maybe?

    If you are solving du/dt - c(u)*u = with backward Euler then it can be written as 

u^{n+1}  - u^n 
 ------------------       - c(u^{n+1})u^{n+1}  = 0
      dt 

  or     (I - dt*c(u^{n+1})u^{n+1})u^{n+1}  = u^n

   Which is the standard form for PETSc's Picard iteration; so you would provide I - dt*c(u^{n+1})u^{n+1} as the compute matrix J() and u^n as the b() Note that the b() is a constant independent of the unknown u^{n+1}.    

    But I do not think it is possible to use TS to solve it this way, you need to use SNES and manage the time-step yourself. Trying to call SNESSetPicard() inside the TS won't work I think. Perhaps Jed or Emil have better thoughts but I think it would require some additional plumbing in TS.

  You can also discretize du/dt - c(u)*u with lazy person's backward Euler as 


u^{n+1}  - u^n 
 ------------------       - c(u^n)u^{n+1}  = 0
      dt 

thus giving the linear system (I - dt c(u^n)) u^{n+1} = u^n but then you need to be sure it is a stable scheme with some accuracy, which it may not be. Again this you can do with KPS but not directly with TS I think.


  If you really are solving du/dt - c(u) = 0 then backward Euler gives you 

u^{n+1}  - u^n 
 ------------------       - c(u^{n+1}) = 0 
      dt 


(n^{n+1} - dt c(u^{n+1}) - u^n  and the logical solver is Newton because the Jacobian of this is 

I - dt J(u^{n+1}) .   This you can give directly to TS, but using Picard doesn't make much sense since there is no natural A(u^{n+1}) u^{n+1} term.


But now I remember something Emil showed me.  If you can write c(x) = A(x) x  + f(x)  then you can do another lazy person's backward Euler with 

u^{n+1}  - u^n 
 ------------------       - A(u^{n+1})u^{n+1}  - f(u^n) = 0
      dt 

Resulting in

(I - dt A(u^{n+1}))u^{n+1} = u^n + f(u^n)   and use SNESSetPicard to solve (again without using TS).

Or if you only want to solve a linear system using


u^{n+1}  - u^n 
 ------------------       - A(u^n)u^{n+1}  - f(u^n) = 0
      dt 

and solve 

(I - dt A(u^n))u^{n+1} = u^n + f(u^n).

But for these to be reasonable approaches you need to know something about A() and f(). 

I am making all these suggestions not because I think doing any of these things is reasonable, I'm just trying to guess at what you are trying to accomplish.

  Barry




 




> 
> 
>  
> 
>    Note that Picard
> 
>        Solves the equation A(x) x = b(x) via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = b(x^{n}) - A(x^{n})x^{n}
> 
>    The functions you need to provide are (compute the vector) b(x)  and (compute the matrix) A(x).   A() is most definitely not the Jacobian of b(). To use Picard you have to rewrite the nonlinear equations that come out of the definition of the time-step in the form A(x) x = b(x).  This may require additions to TS, I don't know.
> 
>    Barry
> 
> 
> 
> > On Jan 1, 2017, at 8:56 AM, Mark Adams <mfadams at lbl.gov> wrote:
> >
> > I tried using Picard inside of TS, and SNES diverged. Is this supposed to be possible?  I tried to do do this by wrapping my operator in a SNES version and simply adding it like this.
> >
> > Mark
> >
> >   ierr = TSSetIFunction(ts,NULL,LandIFunction,ctx);CHKERRQ(ierr);
> >   ierr = TSSetIJacobian(ts,J,J,LandIJacobian,ctx);CHKERRQ(ierr);
> >     SNES snes;
> >     Vec vec;
> >     ierr = VecDuplicate(uu,&vec);CHKERRQ(ierr);
> >     ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
> >     ierr = SNESSetPicard(snes, vec, LandSNESFunction, J, J, LandSNESJacobian, ctx);CHKERRQ(ierr);
> >
> 
>