[petsc-dev] Picard with TS

[Mark Adams]

> > 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
>

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.  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.

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);
> > >
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170103/ee1369f6/attachment.html>