> On Jan 2, 2017, at 4:05 PM, Mark Adams <[email protected]> wrote:
>
>
>
> On Sun, Jan 1, 2017 at 1:00 PM, Barry Smith <[email protected]> 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 <[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);
> >
>
>
