<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
> But I see that IFunction is not correct. I am solving du/dt - C(u) = 0.<br>
<br>
</span> Do you mean du/dt - c(u)* u = 0?<br></blockquote><div><br></div><div>yes,</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<span class=""><br>
> My IFunction applies C(u)*u, but LandSNESFunction should just return zero, maybe?<br>
<br>
</span> If you are solving du/dt - c(u)*u = with backward Euler then it can be written as<br>
<br>
u^{n+1} - u^n<br>
------------------ - c(u^{n+1})u^{n+1} = 0<br>
dt<br>
<br>
or (I - dt*c(u^{n+1})u^{n+1})u^{n+1} = u^n<br></blockquote><div><br></div><div>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.</div><div><br></div><div>Thanks,</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
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}.<br>
<br>
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.<br>
<br>
You can also discretize du/dt - c(u)*u with lazy person's backward Euler as<br>
<br>
<br>
u^{n+1} - u^n<br>
------------------ - c(u^n)u^{n+1} = 0<br>
dt<br>
<br>
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.<br>
<br>
<br>
If you really are solving du/dt - c(u) = 0 then backward Euler gives you<br>
<br>
u^{n+1} - u^n<br>
------------------ - c(u^{n+1}) = 0<br>
dt<br>
<br>
<br>
(n^{n+1} - dt c(u^{n+1}) - u^n and the logical solver is Newton because the Jacobian of this is<br>
<br>
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.<br>
<br>
<br>
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<br>
<br>
u^{n+1} - u^n<br>
------------------ - A(u^{n+1})u^{n+1} - f(u^n) = 0<br>
dt<br>
<br>
Resulting in<br>
<br>
(I - dt A(u^{n+1}))u^{n+1} = u^n + f(u^n) and use SNESSetPicard to solve (again without using TS).<br>
<br>
Or if you only want to solve a linear system using<br>
<br>
<br>
u^{n+1} - u^n<br>
------------------ - A(u^n)u^{n+1} - f(u^n) = 0<br>
dt<br>
<br>
and solve<br>
<br>
(I - dt A(u^n))u^{n+1} = u^n + f(u^n).<br>
<br>
But for these to be reasonable approaches you need to know something about A() and f().<br>
<br>
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.<br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
><br>
><br>
><br>
><br>
> Note that Picard<br>
><br>
> 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}<br>
><br>
> 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.<br>
><br>
> Barry<br>
><br>
><br>
><br>
> > On Jan 1, 2017, at 8:56 AM, Mark Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>> wrote:<br>
> ><br>
> > 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.<br>
> ><br>
> > Mark<br>
> ><br>
> > ierr = TSSetIFunction(ts,NULL,<wbr>LandIFunction,ctx);CHKERRQ(<wbr>ierr);<br>
> > ierr = TSSetIJacobian(ts,J,J,<wbr>LandIJacobian,ctx);CHKERRQ(<wbr>ierr);<br>
> > SNES snes;<br>
> > Vec vec;<br>
> > ierr = VecDuplicate(uu,&vec);CHKERRQ(<wbr>ierr);<br>
> > ierr = TSGetSNES(ts,&snes);CHKERRQ(<wbr>ierr);<br>
> > ierr = SNESSetPicard(snes, vec, LandSNESFunction, J, J, LandSNESJacobian, ctx);CHKERRQ(ierr);<br>
> ><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div></div>