> On Jan 3, 2017, at 7:29 AM, Mark Adams <[email protected]> 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 <[email protected]> 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);
> > >
> >
> >
