[petsc-users] ARKIMEX produces incorrect values

Mon Aug 31 17:27:44 CDT 2020

On 8/31/20 12:17 PM, Ed Bueler wrote:

Emil --

Thanks for looking at this.

> Hi Ed, can you please add the following
> TSSetEquationType(ts,TS_EQ_IMPLICIT);
> before calling TSSolve and try again? This is described in Table 12 in the pdf doc.

Yep, that fixes it.  After setting the TS_EQ_IMPLICIT flag programmatically I get:

It is only programmatic because it has to do with the form of RHS and IFunctions.
$ ./ex54 -ts_type arkimex -ts_arkimex_fully_implicit
error norm at tf = 1.000000 from 12 steps:  |u-u_exact| =  1.34500e-02

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?

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.

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:

IFunction := u_dot and RHSFunction := M^{-1}*H(u) [or solve Mx=H(u) in the RHS function].

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.

So -ts_arkimex_fully_implicit does not set this flag?

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:

1- mass is identity: IFunction:= u_dot-H(u); RHSFunction:= W(u), or

2- mass is full rank, but not identity: IFunction:= M u_dot-H(u); RHSFunction:= M^{-1} * W(u)

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

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.

Emil

> So that we improve our user experience, can you tell us what are your usual sources/starting points
> when implementing a new problem:
> 1- PDF doc

Yes.  Looked briefly at the PDF manual.  E.g. I saw the tables for IMEX methods but my eyes glazed over.

> 2- tutorials (if you find a good match)

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.

> 3- own PETSc implementations

Yes.  I have my own diffusion-reaction system (https://github.com/bueler/p4pdes/blob/master/c/ch5/pattern.c) 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.

> 4- online function doc

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.

> 5- other

Not sure.

Thanks,

Ed

On Mon, Aug 31, 2020 at 6:09 AM Constantinescu, Emil M. <emconsta at anl.gov<mailto:emconsta at anl.gov>> wrote:

On 8/30/20 6:04 PM, Ed Bueler wrote:
Actually, ARKIMEX is not off the hook.  It still gets the wrong answer if told the whole thing is implicit:

$ ./ex54 -ts_type arkimex -ts_arkimex_fully_implicit    # WRONG  (AND REALLY SLOW)
error norm at tf = 1.000000 from 224 steps:  |u-u_exact| =  2.76636e+00

Hi Ed, can you please add the following

 TSSetEquationType<https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSSetEquationType.html#TSSetEquationType>(ts,TS_EQ_IMPLICIT<https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSEquationType.html#TSEquationType>);

before calling TSSolve and try again? This is described in Table 12 in the pdf doc.

So that we improve our user experience, can you tell us what are your usual sources/starting points when implementing a new problem:

1- PDF doc

2- tutorials (if you find a good match)

3- own PETSc implementations

4- online function doc

5- other

Thanks,
Emil

versus

$ ./ex54 -ts_type arkimex      # WRONG BUT IFunction IS OF FLAGGED FORM
error norm at tf = 1.000000 from 16 steps:  |u-u_exact| =  1.93229e+01

$ ./ex54 -ts_type bdf   # RIGHT
error norm at tf = 1.000000 from 33 steps:  |u-u_exact| =  9.29170e-02

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.

Ed

On Sun, Aug 30, 2020 at 2:57 PM Ed Bueler <elbueler at alaska.edu<mailto:elbueler at alaska.edu>> wrote:
Darn, sorry.

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?

Ed

On Sun, Aug 30, 2020 at 2:44 PM Ed Bueler <elbueler at alaska.edu<mailto:elbueler at alaska.edu>> wrote:
Dear PETSc --

I tried twice to make this an issue at the gitlab.com<http://gitlab.com> host site, but both times got "something went wrong (500)".  So this is a bug report by old-fashioned means.

I created a TS example, https://github.com/bueler/p4pdes-next/blob/master/c/fix-arkimex/ex54.c at my github, also attached.  It solves a 2D linear ODE
```
   x' + y' = 6 y
        y' = x
```
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.

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:
$ ./ex54
error norm at tf = 1.000000 from 33 steps:  |u-u_exact| =  9.29170e-02
$ ./ex54 -snes_fd
error norm at tf = 1.000000 from 33 steps:  |u-u_exact| =  9.29170e-02
$ ./ex54 -ts_rtol 1.0e-14 -ts_atol 1.0e-14 -ts_bdf_order 6
error norm at tf = 1.000000 from 388 steps:  |u-u_exact| =  4.23624e-11
$ ./ex54 -ts_type beuler
error norm at tf = 1.000000 from 100 steps:  |u-u_exact| =  6.71676e-01
$ ./ex54 -ts_type cn
error norm at tf = 1.000000 from 100 steps:  |u-u_exact| =  2.22839e-03
$ ./ex54 -ts_type rosw
error norm at tf = 1.000000 from 21 steps:  |u-u_exact| =  5.64012e-03

But it produces wrong values with ARKIMEX:
$ ./ex54 -ts_type arkimex
error norm at tf = 1.000000 from 16 steps:  |u-u_exact| =  1.93229e+01

Neither tightening tolerance nor changing type (`-ts_arkimex_type`) helps ARKIMEX.

Thanks!

Ed

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.

--
Ed Bueler
Dept of Mathematics and Statistics
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
306C Chapman

--
Ed Bueler
Dept of Mathematics and Statistics
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
306C Chapman

--
Ed Bueler
Dept of Mathematics and Statistics
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
306C Chapman

--
Emil M. Constantinescu, Ph.D.
Computational Mathematician
Argonne National Laboratory
Mathematics and Computer Science Division

Ph: 630-252-0926
http://www.mcs.anl.gov/~emconsta

--
Ed Bueler
Dept of Mathematics and Statistics
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
306C Chapman

--
Emil M. Constantinescu, Ph.D.
Computational Mathematician
Argonne National Laboratory
Mathematics and Computer Science Division

Ph: 630-252-0926
http://www.mcs.anl.gov/~emconsta
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200831/a22d876a/attachment-0001.html>