<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body>
<p>On 8/31/20 12:17 PM, Ed Bueler wrote:<br>
</p>
<blockquote type="cite" cite="mid:CAOHboJ_cZ2xs2trjRnrLdWHUZ2BXn6mN0TKsbbrYGe94iohnUA@mail.gmail.com">
<div dir="ltr">Emil --
<div><br>
</div>
<div>Thanks for looking at this.</div>
<div><br>
</div>
> Hi Ed, can you please add the following<br>
> TSSetEquationType(ts,TS_EQ_IMPLICIT);<br>
> before calling TSSolve and try again? This is described in Table 12 in the pdf doc.
<div><br>
</div>
<div>Yep, that fixes it. After setting the TS_EQ_IMPLICIT flag programmatically I get:</div>
<div><br>
</div>
</div>
</blockquote>
It is only programmatic because it has to do with the form of RHS and IFunctions.<br>
<blockquote type="cite" cite="mid:CAOHboJ_cZ2xs2trjRnrLdWHUZ2BXn6mN0TKsbbrYGe94iohnUA@mail.gmail.com">
<div dir="ltr">
<div>
<div>$ ./ex54 -ts_type arkimex -ts_arkimex_fully_implicit<br>
error norm at tf = 1.000000 from 12 steps: |u-u_exact| = 1.34500e-02<br>
</div>
<div><br>
</div>
</div>
<div>Without -ts_arkimex_fully_implicit we still get the wrong answer, but, as I understand it, we expect the wrong answer because dF/d(dudt) != I, correct?</div>
<div><br>
</div>
</div>
</blockquote>
<p>Yes. I keep mixing F and G, but if you want to solve Mu'=H(u), then define the IFunction := M u_dot - H(u) then it should work with all time steppers.
<br>
</p>
<p>If you want to set the RHS of your ODE in the RHS function (so that you can use explicit integrators, too) you have to provide:
<br>
</p>
<p>IFunction := u_dot and RHSFunction := M^{-1}*H(u) [or solve Mx=H(u) in the RHS function].
<br>
</p>
<p>Note that M u_dot - H(u) can only be solved by implicit solvers directly so IFunction := M u_dot and RHSFunction := H(u). Table 12 in the PDF doc explains these cases, but that can be improved as well.
</p>
<blockquote type="cite" cite="mid:CAOHboJ_cZ2xs2trjRnrLdWHUZ2BXn6mN0TKsbbrYGe94iohnUA@mail.gmail.com">
<div dir="ltr">
<div>So -ts_arkimex_fully_implicit does not set this flag?</div>
<div><br>
</div>
</div>
</blockquote>
No, its use is for when you have both IFunction (for stiff) and RHSfunction (for nonstiff) defined to solve Mu'=H(u) + W(u) and:<br>
<p>1- mass is identity: IFunction:= u_dot-H(u); RHSFunction:= W(u), or<br>
</p>
<p>2- mass is full rank, but not identity: IFunction:= M u_dot-H(u); RHSFunction:= M^{-1} * W(u)</p>
<p>and you have a choice of using either an IMEX scheme [-ts_arkimex_fully_implicit false] or just the implicit part [-ts_arkimex_fully_implicit true].<br>
</p>
<p>Thank you for your feedback on our short survey - it is very valuable in helping us crafting a less painful path to using all these options.</p>
<p>Emil<br>
</p>
<blockquote type="cite" cite="mid:CAOHboJ_cZ2xs2trjRnrLdWHUZ2BXn6mN0TKsbbrYGe94iohnUA@mail.gmail.com">
<div dir="ltr">> So that we improve our user experience, can you tell us what are your usual sources/starting points
<div>> when implementing a new problem:<br>
> 1- PDF doc</div>
<div><br>
Yes. Looked briefly at the PDF manual. E.g. I saw the tables for IMEX methods but my eyes glazed over.</div>
<div><br>
> 2- tutorials (if you find a good match)<br>
<br>
</div>
<div>Yes. Looked at various html pages including the one for TSARKIMEX. But I missed the sentence "Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X)." I did not expect that ARKIMEX
had this restriction, and did not pick it up.<br>
<br>
> 3- own PETSc implementations<br>
<br>
Yes. I have my own diffusion-reaction system (<a href="https://github.com/bueler/p4pdes/blob/master/c/ch5/pattern.c" moz-do-not-send="true">https://github.com/bueler/p4pdes/blob/master/c/ch5/pattern.c</a>) in which ARKIMEX works well. (Or at least as far
as I can tell. I don't have a manufactured solution for it, for example.) I am in the midst of tracking down a different kind of error, probably from DMDA callbacks, when I got distracted by the current issue.<br>
<br>
> 4- online function doc<br>
<br>
Yes. See above comment on TSARKIMEX page. By my memory I also looked at the TSSet{I,RHS}Jacobian() pages, for example, and probably others.</div>
<div><br>
> 5- other</div>
<div><br>
</div>
<div>Not sure.</div>
<div><br>
</div>
<div>Thanks,</div>
<div><br>
</div>
<div>Ed</div>
<div><br>
</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Mon, Aug 31, 2020 at 6:09 AM Constantinescu, Emil M. <<a href="mailto:emconsta@anl.gov" moz-do-not-send="true">emconsta@anl.gov</a>> wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<p><br>
</p>
<div>On 8/30/20 6:04 PM, Ed Bueler wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">Actually, ARKIMEX is not off the hook. It still gets the wrong answer if told the whole thing is implicit:
<div><br>
</div>
<div>$ ./ex54 -ts_type arkimex -ts_arkimex_fully_implicit # WRONG (AND REALLY SLOW)<br>
error norm at tf = 1.000000 from 224 steps: |u-u_exact| = 2.76636e+00<br>
<br>
</div>
</div>
</blockquote>
Hi Ed, can you please add the following
<pre width="80" style="color:rgb(0,0,0);font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;word-spacing:0px;text-decoration-style:initial;text-decoration-color:initial"> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSSetEquationType.html#TSSetEquationType" target="_blank" moz-do-not-send="true">TSSetEquationType</a>(ts,<a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSEquationType.html#TSEquationType" target="_blank" moz-do-not-send="true">TS_EQ_IMPLICIT</a>);</pre>
<p>before calling TSSolve and try again? This is described in Table 12 in the pdf doc.<br>
</p>
<p><br>
</p>
<p>So that we improve our user experience, can you tell us what are your usual sources/starting points when implementing a new problem:</p>
<p>1- PDF doc</p>
<p>2- tutorials (if you find a good match)<br>
</p>
<p>3- own PETSc implementations</p>
<p>4- online function doc</p>
<p>5- other<br>
</p>
<p>Thanks,<br>
Emil<br>
</p>
<blockquote type="cite">
<div dir="ltr">
<div>versus</div>
<div><br>
</div>
<div>$ ./ex54 -ts_type arkimex # WRONG BUT IFunction IS OF FLAGGED FORM<br>
error norm at tf = 1.000000 from 16 steps: |u-u_exact| = 1.93229e+01<br>
</div>
<div><br>
</div>
<div>$ ./ex54 -ts_type bdf # RIGHT<br>
error norm at tf = 1.000000 from 33 steps: |u-u_exact| = 9.29170e-02<br>
</div>
<div><br>
</div>
<div>So I am not sure what "Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X)." means.</div>
<div><br>
</div>
<div>Ed</div>
<div><br>
</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Sun, Aug 30, 2020 at 2:57 PM Ed Bueler <<a href="mailto:elbueler@alaska.edu" target="_blank" moz-do-not-send="true">elbueler@alaska.edu</a>> wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px
0px 0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div dir="ltr">Darn, sorry.
<div><br>
</div>
<div>I realize the ARKIMEX page does say "Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X)." So my example does not do that. Is there a way for ARKIMEX to detect that dG/d(Xdot) = I?
<div><br>
</div>
<div>Ed</div>
</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Sun, Aug 30, 2020 at 2:44 PM Ed Bueler <<a href="mailto:elbueler@alaska.edu" target="_blank" moz-do-not-send="true">elbueler@alaska.edu</a>> wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px
0px 0px 0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div dir="ltr">
<div>Dear PETSc --</div>
<div><br>
</div>
<div>I tried twice to make this an issue at the <a href="http://gitlab.com" target="_blank" moz-do-not-send="true">gitlab.com</a> host site, but both times got "something went wrong (500)". So this is a bug report by old-fashioned means.</div>
<div><br>
</div>
I created a TS example, <a href="https://github.com/bueler/p4pdes-next/blob/master/c/fix-arkimex/ex54.c" target="_blank" moz-do-not-send="true">
https://github.com/bueler/p4pdes-next/blob/master/c/fix-arkimex/ex54.c</a> at my github, also attached. It solves a 2D linear ODE<br>
```<br>
x' + y' = 6 y<br>
y' = x<br>
```<br>
Pretty basic; the known exact solution is just exponentials. The code writes it as F(t,u,u')=G(t,u) and supplies all the pieces, namely IFunction,IJacobian,RHSFunction,RHSJacobian. Note both F and G must be seen by TS to get the correct solution. In summary,
a boring (and valgrind-clean ;-)) example.
<div><br>
For current master branch it runs fine for the fully-implicit methods (e.g. BDF, CN, ROSW) which can use the IFunction F, including with finite-differenced Jacobians. With BDF2, BDF2+-snes_fd, BDF6+tight tol., CN, BEULER, ROSW:<br>
$ ./ex54<br>
error norm at tf = 1.000000 from 33 steps: |u-u_exact| = 9.29170e-02<br>
$ ./ex54 -snes_fd<br>
error norm at tf = 1.000000 from 33 steps: |u-u_exact| = 9.29170e-02<br>
$ ./ex54 -ts_rtol 1.0e-14 -ts_atol 1.0e-14 -ts_bdf_order 6<br>
error norm at tf = 1.000000 from 388 steps: |u-u_exact| = 4.23624e-11<br>
$ ./ex54 -ts_type beuler<br>
error norm at tf = 1.000000 from 100 steps: |u-u_exact| = 6.71676e-01<br>
$ ./ex54 -ts_type cn<br>
error norm at tf = 1.000000 from 100 steps: |u-u_exact| = 2.22839e-03<br>
$ ./ex54 -ts_type rosw<br>
error norm at tf = 1.000000 from 21 steps: |u-u_exact| = 5.64012e-03<br>
<br>
But it produces wrong values with ARKIMEX:<br>
$ ./ex54 -ts_type arkimex<br>
error norm at tf = 1.000000 from 16 steps: |u-u_exact| = 1.93229e+01<br>
<div><br>
Neither tightening tolerance nor changing type (`-ts_arkimex_type`) helps ARKIMEX.<br>
<br>
Thanks!</div>
<div><br>
Ed</div>
<div><br>
</div>
<div>PS My book is at a late proofs stage, and out of my hands. It should appear SIAM Press in a couple of months. In all the examples in my book, only my diffusion-reaction system example using F(t,u,u') = G(t,u) is broken. Thus the motivation for a trivial
ODE example as above.</div>
<div><br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>Ed Bueler<br>
Dept of Mathematics and Statistics<br>
University of Alaska Fairbanks<br>
Fairbanks, AK 99775-6660<br>
306C Chapman<br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>Ed Bueler<br>
Dept of Mathematics and Statistics<br>
University of Alaska Fairbanks<br>
Fairbanks, AK 99775-6660<br>
306C Chapman<br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>Ed Bueler<br>
Dept of Mathematics and Statistics<br>
University of Alaska Fairbanks<br>
Fairbanks, AK 99775-6660<br>
306C Chapman<br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
<pre cols="72">--
Emil M. Constantinescu, Ph.D.
Computational Mathematician
Argonne National Laboratory
Mathematics and Computer Science Division
Ph: 630-252-0926
<a href="http://www.mcs.anl.gov/~emconsta" target="_blank" moz-do-not-send="true">http://www.mcs.anl.gov/~emconsta</a></pre>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr" class="gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>Ed Bueler<br>
Dept of Mathematics and Statistics<br>
University of Alaska Fairbanks<br>
Fairbanks, AK 99775-6660<br>
306C Chapman<br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
<pre class="moz-signature" cols="72">--
Emil M. Constantinescu, Ph.D.
Computational Mathematician
Argonne National Laboratory
Mathematics and Computer Science Division
Ph: 630-252-0926
<a class="moz-txt-link-freetext" href="http://www.mcs.anl.gov/~emconsta">http://www.mcs.anl.gov/~emconsta</a></pre>
</body>
</html>