Quality of algorithms

[sourceryinstitute]

I just came from a meeting where I met a mathematics professor who specializes in ODE's. When I showed him FOODiE, he commented that the algorithms are not state-of-the-art and that he would not use any of them. Although he teaches the algorithms in FOODiE in an introductory undergraduate course, even that undergraduate course progresses to more sophisticated algorithms by the end. He indicated that more advanced algorithms vary both the time step and the order of accuracy adaptively. He also indicated that it's important to distinguish stiff from non-stiff solvers and to distinguish solvers that target ODE's generally from solvers developed specifically for the ODE systems the arise from semi-discrete PDE formulations.

His comments raise several important considerations:

1. It's important that specialists in the field recognize the library as valuable. In other conversations, I've heard criticisms of the algorithms in even one of the most widely used numerical methods textbooks so it's quite possible for poor algorithms to reach a very wide audience and thereby do more harm than good.
2. My hunch is that the defined operators that have been written for the current ODE solvers will enable the rapid implementation of more sophisticated solvers. If I'm right, then the cost of replacing the current algorithms with better ones is relatively low and the opportunity should not be missed before submitting the announcement paper. If I'm wrong, then the operators themselves are not as useful as one would hope.
3. There should be an ODE specialist involved in this project. I wasn't able to interest him in joining the project, but he did at least provide a few citations:

(a) http://www.netlib.org/ode/rksuite/
(b) http://www.unige.ch/~hairer/software.html (he specifically recommended the DOPRI5 and DOP853 algorithms from this page)
(c) Dormand, John R., and Peter J. Prince. "A family of embedded Runge-Kutta formulae." Journal of computational and applied mathematics 6.1 (1980): 19-26.

1. He also indicated that some or all of the algorithms in FOODiE don't allow for the right-hand-side function to depend on the solution. I haven't read any of the code in FOODiE and am unaware of whether this is true. Hopefully he missed something here. He made all of these comments after just a few minutes of consideration.

I also mentioned one class of adaptive schemes I've used: Runga-Kutta-Fehlberg (RKF) schemes. He indicated that even those have been superseded by better methods. Fehlberg published the first RKF scheme in 1969, where as Dormand and Prince published their scheme in reference (c) in 1980. I can at least say that in the case of RKF methods, the hunch I stated in item 2 above is true. I haven't investigated whether it's true for the Dormant-Prince algorithm.

Although I have taught graduate-level numerical methods and have written semi-discrete PDE solvers, I'm clearly not an ODE expert and sadly I have not had the time to contribute to this project even though I admire the effort and appreciate that it builds on software design ideas that I published.

I urge serious reconsideration of the algorithms being presented in FOODiE. It seems that more sophisticated algorithms should be presented as the primary and preferred interface of FOODiE while the current algorithms could be background material and cited for instructional purposes only.

[szaghi]

Dear Damian,
Thank you very much for these comments and suggestions. I try to replay quickly now, but monday I try to give more detailed answer:

• I cannot say anything about the your expert, but from a superficial first impression on his comments he seems an awarded professor who devised cutting-edge schemes, but never solved anything more complex than a uniform 1D smooth, subsonic heat pipe...
• I have not the authority to be defined an expert, but I know experts (@andreadimascio is one of the best who developed more than 10 years ago a CFD code that is still now a state-of-the-art CFD code): when you deal with real complex problem you start using well known schemes (RK54 and similar embeded version and bacward formulae for stiff ptoblem)
• many of the schemes your expert has citeted are planned to be included in FOODiE, see ODE schemes: a list of schemes to be implemented #1 embedded RK is one of the my priority; however it is false that the already implemented schemes are not currently used: in weather forecasting leapfrog, Adams-Bashforth and in general linear multi-step method are the defacto standard; in gasdynamic RK54 is probably the most diffused;
• it is a very superficial comment to say that the current FOODiE solvers cannot allow for the right-hand-side function to depend on the solution: what this means? FOODiE is already proven to solve the most general form of non linear ODE like Du_t = R(t,u(t)) whe the rhs is a residual function non linearly depending on the solution and on time; moreover there is already an example solving the compressible multi-specie Euler PDE system... FOODiE does not only allow but help on solving non linear PDE.
• the ADT integrand API will be improved soon: for embedded RK and implicit solvers some other operators and methods should be necessary, but I feel that it is quite stable.

I agree with your expert that more cutting-edge schemes must be implemented, but it is just a matter of time: let me few lunch breaks and many solvers will come (I am working on FOODIE in my lunch break), but I cannot agree that the already implemented ones are not useful (just see what is used in almost all CFD RANSE simulations of real complex problems such as a fully appended ship manouvring with a fully discretized propeller...). The current solvers cannot be relegated to some historical citation of what we do last decades, simply because it is not true: the fact that new more sofisticated solvers have been already presented cannot automatically relegate all other to the history. In my (modest) opinion, there is not a best scheme for all problems (the great debate on Riemann solvers is an emblematic example): I will not surprise that someone still use leapfrog today, even if IMEXRK methods are here...

Ok this was the short answer :-) monday I will try to give more details :-)

As always, thank you very much for your help.

P.S. one expert that I would like to joing us is David Ketcheson (https://github.com/ketch), but I have not the prestige to contact him: maybe you can try to write him, your emails will not tagged to be spam as mine :-)

[milancurcic]

@sourceryinstitute Hi Damian, thank you so much for your and your colleague's comments. I think they are on point and well taken. It is true that the solvers implemented in FOODiE are textbook classics, and are not state-of-the-art methods, especially for mathematicians that specialize in ODEs. However, I do think that your hunch (2.) is correct - the software design of FOODiE will allow for straightforward implementation of more sophisticated methods.

In my opinion, FOODiE should not be advertised as a swiss-army knife collection of solvers to solve and do research on ODEs, at least not just yet. I am sure there exist sophisticated integration methods that neither you nor Stefano nor me have heard of before. Instead, I think FOODiE should be advertised as a library that facilitates implementation, testing, and analysis of sophisticated solvers. @szaghi is correct that different applications call for different methods. Solvers used by mathematicians are different than those used by aerodynamicists, which are different from those used by meteorologists and oceanographers (I'm one of those). Most models that we use and develop use 1st and 2nd order schemes like implicit Euler, leapfrog, AB2. RK4 is used but considered sophisticated. This is not because we don't know better but because errors from turbulence closure schemes, or cloud convection and microphysics parameterizations, far outweigh those from differentiation and integration algorithms, and computational efficiency is of prime importance.

Again, I think your hunch (2.) is correct, though it's possible that API may have to be revised a few more times to facilitate implementation of more sophisticated algorithms.

@szaghi Your comment may be a bit too defensive and assuming but I know you are well-meaning. Comments and criticisms like these from Damian's colleague are gold at this stage of FOODiE development and will help us a lot with going forward. I also recognize David Ketcheson as an expert whose insights would be of tremendous help.

@ketch David, when you get a chance, can you please take a quick glance at FOODiE and share some thoughts on the above? Thank you in advance!

[szaghi]

@milancurcic I am sorry, if my statements were too irrispectous, this was not my intention. I think that critics are more important than compliment, thus I will never stop to thank Damian, his expert, you and all others firstly for your critics, that not only in this stage, but ever are gold. I am very sorry if someone is hurted by my words.

The defense of obsolete schemes will come monday, but I like to point out for Damian and others one key point: almost all the current solvers have been implemented in few minutes, this the great merit of Damian idea. In my modest opinion, as Milan stated, FOODiE has great potential not as swiss-army, but for helping to develop new schemes and solve new problems in a rapid, concise and clear approach. I have not yet found the time to write a decent introduction on the why I need FOODiE. I depict here some thouths. Before reading Damian's book, I was a scissed man: I was attracted by python beacuse it is ridicolously simple to gain the abstraction for Rapid Application Development, but python (if not hackered mixing c and fortran) has very bad performance if compared to compiled language as fortran that was the other planet attracting me. After reading Damian's book I realized that also Fortran provide a sufficient level of abstraction with the plus to be the fastest language into the west.

Commonly there are two main categories of numerical specialists:

• mathematicians that use matlab, scilab, python for devising new great schemes; they need to test their idea quickly, without much effort, on simple handy cases;
• engineers, or the oompa loompa workers like me that use fortran or other inflexible but performant language for solving real complex problems; they typically hard-code matematicians great idea on each specific case;

I would like that FOODiE demonstrates that devising new idea is also possible with Fortran and in the meanwhile I would like to prove that solving a real ODE problem is very simple with the built-in arsenal of FOODiE.

[szaghi]

@ketch Milan speaks only on his benhalf... I would never disturb you for a toy like FFODiE :-) I am joking, but I seriously take your work with great ammiration.

[sourceryinstitute]

@szaghi and @milancurcic, thank you for taking my comments in the helpful manner in which they were intended. I don't interpret any of the responses as overly defensive. Your responses are valid ones. Given the physical uncertainties in complex meteorological models, for example, I can certainly understand that the time discretization errors are not first order.

For what it's worth, I think the person who made the comment is tackling very challenging ODE problems in an area of direct relevance to science, i.e., not toy problems such as 1D heat conduction. Nonetheless, he made it clear that he is not an expert in solving the ODEs that result from the semi-discretization of PDEs.

Based on your responses and on my own experience CFD and multiphysics modeling, I think the greatest value of the numerical methods you have implemented so far is their utility in solving very important PDE problems as evidenced by their use in the applications you've cited. It might be the case that the only change you need to make is a documentation change (and possibly a project name change, but that's radical and I'm just suggesting it -- not pushing for it). Most importantly, you might consider stating somewhere very prominently in the documentation what the range of applications are that inspire the methods your project. Let people know right away at the top of the README file that you are implementing schemes of interest in CFD, meteorology, and other multiphysics modeling areas involving the semi-discretization of PDEs. You can then make it clear that the methods being implemented first could also be used for ODEs that do not result from discretizing PDEs, but that the more challenging ODE problem domains warrant more sophisticated numerical methods. You could then emphasize that the operators you're defining on the simpler methods will make the future implementation of more sophisticated ODE integrators reasonably straightforward.

Keep up the great work. I think you're on the right path and I hope I can free up enough time to contribute at some point. For now, I do have one minor terminology correction to make in your README file's reference to my book. I'll post that in a separate Issue. Please also let me know how you would prefer to receive corrections if I were to make them myself directly. I'm a fan of the Forking Workflow described at https://www.atlassian.com/git/tutorials/comparing-workflows/forking-workflow. I recognize it's complicated and therefore not for ever project, but it's especially useful for "outside" developers, which is what I would currently consider myself. In the Forking Workflow, outside developers fork the repository, make changes and then submit pull requests.

Damian

[szaghi]

@sourceryinstitute agree on all. Tomorrow more detailed replay will come. Feel free to correct, update, change, rewrite, delete... all as you prefer (forking and pull is welcome): Milan directly correct my errors pushing his correction upstream, you should have also the same permissions of Milan and mine (the project belongs to all). Thank very very very much!

[rouson]

It just dawned on me that there's a convenient renaming that you might like. It's "Fortran Object-Oriented Differential-equation Integration Environment (FOODIE)", which has the following benefits:

1. Your background and interests (as well as mine) relate primarily to the ODEs that result from semi-discretization of PDEs. That puts the project in a middle ground where neither ODE nor PDE quite fits. Just saying "differentiate equations" without being more specific might be the best compromise.
2. You get to keep almost the same acronym.
3. As a bonus, the ordering of the letters in the acronym now matches the ordering of the words.

Thoughts?

:D

[milancurcic]

@rouson I love it. Of course, it will be Stefano's call.

[szaghi]

@sourceryinstitute @milancurcic you are great! I agree, just let Zaak to say his thougths.

[szaghi]

Hi all, yesterday I forget to announce that I have changed the project name as Damian suggested without waiting @zbeekman opinion: I hope Zaak does not matter about this.

[rouson]

Bravo!

:D

[szaghi]

Dear all,

the explicit embedded Runge-Kutta class of solvers is within us! As the Damian's expert anticipated, it works very well. For more details on what I have done read this post.

Other embedded RK solvers will come soon together the CAF enabled Euler 1D test.

See you soon.

[rouson]

This is fabulous progress.

[szaghi]

@rouson and all,

new embedded RK solvers up to 10th order are here!

See you soon.

[sourceryinstitute]

Great work!

D

Sent from my iPhone

On Oct 29, 2015, at 9:05 AM, Stefano Zaghi [email protected] wrote:

@rouson and all,

new embedded RK solvers up to 10th order are here!

See you soon.

—
Reply to this email directly or view it on GitHub.

[milancurcic]

@szaghi Amazing! My head is spinning from all these orders of accuracy :)

[zbeekman]

Hi everyone, I'm trying to get caught up on this while I have a moment.

I like the renaming scheme @rouson proposed. I also think that, as we all discussed on the call, that our usage of FOODIE is primarily directed at solving PDEs/multi-physics, CFD, numerical meteorology/oceanography etc. and that we need to emphasize that in the documentation, presentations and papers.

Also, we should emphasize the flexibility, modularity and ease of adding new algorithms, potentially including those that vary the time-step as well as order of accuracy.

@szaghi thank you for your diligence and hard work.

I would like to try my hand at adding a second order accurate low-storage RK from the Williamson JCP article that we use for our turbulence simulations. After that, I would like to investigate whether or not it is possible to generalize and implement my PhD advisor's DP2 (data parallel, second order accurate implicit) scheme, as well as some of the data-parallel schemes preceding it like Yoon & Jameson's LU-SGS and Candler & Weir's DP-LUR. These are formulated for CFD, and I haven't gone about looking at whether or not they could be generalized. Also the implicit Data Parallel line relaxation time integration used in NASA's DPLR code. At any rate, I've just started thinking of this and I'm not sure whether it's feasible to implement them in an ADT Calculus pattern manner.

[szaghi]

@zbeekman great news! Do not worry about the, take your time: your solver is very interesting and welcome, but think firstly to you Ph.D. defense, we will not escape :-)

Very soon I will add the corrections we talk about.

See you soon.

[szaghi]

I this can be closed, feel free to re-open if necessary.