[petsc-dev] Picard with TS

[Barry Smith]

> On Jan 3, 2017, at 7:29 AM, Mark Adams <mfadams at lbl.gov> wrote:
> 
> 
> 
> 
> > 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?
> 
> yes,
>  
> 
> >  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

    I have an extra u^{n+1} in the above it should be 

      (I - dt*c(u^{n+1}))u^{n+1}  = u^n
> 
> OK, this makes sense, but I think I can use backward Euler as you say below. I am just trying to match what others are doing to make sure that I am not missing something algorithmically.  

     This is a good idea. Sometimes paper's don't give enough clear information to determine exactly what they are doing.

> They are using BE, as I recall now. I am now using C-N now. I think I can assure myself that BE with Newton is pretty much the same as Picard and BE.

   Not sure what you mean by "pretty much the same", but they are definitely IMHO pretty different. The matrix created in each case is different and the convergence of both the linear and nonlinear iteration will be different. Hopefully the Newton convergence will be quadratic while the Picard will not be.

   Barry

> 
> Thanks,
>  
> 
>    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);
> > >
> >
> >