ODE schemes: a list of schemes to be implemented

[szaghi]

Just a list of schemes that we would like to be implemented into the library:

• [x] explicit schemes:
  • [x] Adams-Bashforth schemes see [7]:
  • [x] 1 step, namely the forward explicit Euler scheme, 1st order accurate;
  • [x] 2 steps, 2nd accurate;
  • [x] 3 steps, 3rd accurate;
  • [x] 4 steps, 4th accurate;
  • [x] Euler scheme, 1st order accurate;
  • [x] Leapfrog, 2nd order accurate:
  • [x] unfiltered leapfrog, 2nd order accurate, mostly unstable, see [4];
  • [x] Robert-Asselin filtered leapfrog, 1st order accurate, see [4, 5, 6];
  • [x] Robert-Asselin-Williams filtered leapfrog, 3rd order accurate, see [5, 6];
  • [ ] Runge-Kutta schemes:
  • [ ] embedded (adaptive) schemes:
    • [x] Heun-Euler, 2 stages, 2nd order accurate;
    • [ ] Runge-Kutta-Fehlberg, 5 stages, 4th order accurate, see [8];
    • [x] Cash-Karp, 6 stages, 4th order accurate, see[14];
    • [x] Prince-Dormand, 7 stages, 4th order accurate, see[10];
    • [x] Calvo et al., 9 stages, 6th order accurate, see[13];
    • [x] Feagin, 17 stages, 10th order accurate, see[15];
  • [ ] low-storage schemes, see [1, 2, 3]:
    • [x] 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
    • [ ] 2 stages;
    • [ ] 3 stages;
    • [ ] 4 stages;
    • [x] 5 stages, 4th order accurate, 2N registers, see [3];
    • [x] 6 stages, 4th order accurate, 2N registers, see [11];
    • [x] 7 stages, 4th order accurate, 2N registers, see [11];
    • [x] 12 stages, 4th order accurate, 2N registers, see [12];
    • [x] 13 stages, 4th order accurate, 2N registers, see [12];
    • [x] 14 stages, 4th order accurate, 2N registers, see [12];
  • [ ] TVD or SSP schemes, see [1]:
    • [x] 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
    • [x] 2 stages, 2nd order accurate;
    • [x] 3 stages, 3rd order accurate;
    • [ ] 4 stages;
    • [x] 5 stages, 4th order accurate;
• [ ] implicit schemes:
  • [x] Adams-Moulton schemes:
  • [x] 0 step, namely the backward implicit Euler scheme, 1st order accurate;
  • [x] 1 step, 2nd accurate;
  • [x] 2 steps, 3rd accurate;
  • [x] 3 steps, 4th accurate;
• [ ] predictor-corrector schemes:
  • [x] Adams-Bashforth-Moulton schemes:
  • [x] 1 step, AB(1)-AM(0), 1st order accurate;
  • [x] 2 steps, AB(2)-AM(1), 2nd accurate;
  • [x] 3 steps, AB(3)-AM(2), 3rd accurate;
  • [x] 4 steps, AB(4)-AM(3), 4th accurate;
  • [ ] Runge-Kutta schemes;

    Bibliography

[1] High Order Strong Stability Preserving Time Discretizations, Gottlieb, S., Ketcheson, D. I., Shu, C.W., Journal of Scientific Computing, vol. 38, N. 3, 2009, pp. 251-289.

[2] Low-Storage Runge-Kutta Schemes, J. H. Williamson, Journal of Computational Physics, vol. 35, 1980, pp. 48--56.

[3] Fourth-Order 2N-Storage Runge-Kutta Schemes, Mark H. Carpenter, Christopher A. Kennedy, NASA Technical Memorandum 109112, June 1994.

[4] Mesinger F. and A. Arakawa, 1976: Numerical methods used in atmospheric models. Global Atmospheric Research Programme (GARP) Technical Report. http://twister.ou.edu/CFD2003/Mesinger_ArakawaGARP.pdf

[5] Williams, P. D., 2009: A Proposed Modification to the Robert–Asselin Time Filter. Mon. Wea. Rev., 137, 2538–2546. doi: http://dx.doi.org/10.1175/2009MWR2724.1

[6] The RAW filter: An improvement to the Robert–Asselin filter in semi-implicit integrations, Williams, P.D., Monthly Weather Review, vol. 139(6), pages 1996--2007, June 2011.

[7] Linear multistep method, wikipedia article

[8] Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems, Erwin Fehlberg (1969), NASA Technical Report 315.

[9] A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, J. R. Cash, A. H. Karp, ACM Transactions on Mathematical Software 16: 201-222, 1990. doi:10.1145/79505.79507

[10] A family of embedded Runge-Kutta formulae, Dormand, J. R.; Prince, P. J. (1980), , Journal of Computational and Applied Mathematics 6 (1): 19--26, doi:10.1016/0771-050X(80)90013-3

[11] High-accuracy large-step explicit Runge–Kutta (HALE-RK) schemes for computational aeroacoustics, Vasanth Allampalli and Ray Hixon and M. Nallasamy and Scott D. Sawyer, Journal of Computational Physics, vol. 228, 2009, pp. 3837--3850.

[12] Efficient low-storage Runge–Kutta schemes with optimized stability regions, Jens Niegemann and Richard Diehl and Kurt Busch, Journal of Computational Physics, vol. 231, 2012, pp. 364--372.

[13] A New Embedded Pair of Runge-Kutta Formulas of orders 5 and 6, M. Calvo, J.I. Montijano, L. Randez, Computers & Mathematics with Applications, Volume 20, Issue 1, 1990, Pages 15--24, ISSN 0898-1221, http://dx.doi.org/10.1016/0898-1221(90)90064-Q.

[14] A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, J. R. Cash, A. H. Karp, ACM Transactions on Mathematical Software, vol. 16, pp. 201--222, 1990, doi:10.1145/79505.79507.

[15] A tenth-order Runge-Kutta method with error estimate, Feagin, T., Proceedings of the IAENG Conf. on Scientific Computing. 2007.

[zbeekman]

Can I add low-storage explicit RK schemes? http://www.sciencedirect.com/science/article/pii/0021999180900339 and maybe these: http://www.ece.uvic.ca/~bctill/papers/numacoust/Carpenter_Kennedy_1994.pdf

[szaghi]

Sure, you can add all you want!

I'll update references soon.

[milancurcic]

Some suggestions for 2nd and 3rd order schemes:

• [ ] Leapfrog, 2nd order accurate
  • [ ] Robert-Asselin filter, for stability and leapfrog -> 1st order accuracy
  • [ ] Robert-Asselin-Williams filter, for stability and leapfrog -> 3rd order accuracy
• [ ] Adams-Bashforth variants of 2nd and 3rd order accuracy

[szaghi]

Thank you very much @milancurcic

[szaghi]

Hi all,

I have done some steps over. The class of RK schemes that I usually use has been implemented (documentation coming very soon).

I am now extending the tests suite with the Burgers equation (and maybe some others). Then I would like to try to implement other schemes that you suggested: the low-storage class is quite in my radar. On the contrary for the schemes proposed by @milancurcic I need help: I have some implementation of those schemes from my old university references, but I prefer to have a more well known and opinionated references.

@milancurcic can you suggest some references?

[milancurcic]

@szaghi Sorry I didn't get back to you sooner on this - somehow I missed the notification.

This is the text that we used in our Numerical Methods course in college: http://twister.ou.edu/CFD2003/Mesinger_ArakawaGARP.pdf , and it describes the basics of numerical analysis, and the most basic integration and differentiation schemes used in atmospheric and oceanic models. In the text they are applied to the oscillation equation, and this may be a good candidate to include in FOODiE testing set as well.

I also like Dale Durran's book: http://www.atmos.washington.edu/numerical_methods/, but it is fairly expensive. Most academic libraries should have it though.

If you're not in a rush to implement these few basic schemes I suggested, let me work on it over the next week or two. What do you think?

[milancurcic]

The above references are more of a review/textbook format kind of texts. Specific schemes have more appropriate references, and I will include them as I go in the initial implementation.

[szaghi]

@milancurcic Thank you very much! Do not worry about delayed responses... I am not so responsive in general.

Mesinger et al. lectures seems nice (not so up-to-date...) I will read it under the sun... Instead, I should have a copy of the Durran's book, I will check. Anyhow, for your solvers I will refer to the schemes described into these two references.

I will try to implement the solvers you suggested not so quickly (two weeks I think, because I am going to take my holidays), but whatever I will do, you are free to do whatever you want with FOODiE: if at same point my implementations are bad, feel free to castigate me! Or start your implementations whenever you want, push on this repository and I will pull your work. I think that you should have write permissions.

Thank you again for your help!

[szaghi]

@milancurcic Adams-Bashforth class is here! It is a very basic implementation, I do not have the time to check literature for the best, state-of-the-art implementation, I have just followed the wikipedia reference... I am sorry for my laziness.

I have updated all the current 3 tests (burgers, lorenz and oscillation). Note that I followed your suggestion and I have embedded the previous time steps into the %state member of the 3 concrete extensions. The upper bound index of the second subscript is the up-to-date state, i.e. state(n+s) where s is the number of time steps(levels) you are using. So, considering your oscillation test, the type_oscillation is now:

type, extends(integrand) :: oscillation
  !< Oscillation equations field.
  !<
  !< It is a FOODiE integrand class.
  private
  real(R_P), dimension(:,:), allocatable :: state      !< Solution vector, [1:state_dims,1:time_steps_stored].
  real(R_P)                              :: f = 0._R_P !< Oscillation frequency (Hz).
  contains
    procedure, pass(self), public :: output                                                           !< Extract Oscillation field.
    procedure, pass(self), public :: update_previous_steps                                            !< Update previous time steps.
    procedure, pass(self), public :: t => dOscillation_dt                                             !< Time derivative, residuals.
    procedure, pass(lhs),  public :: integrand_multiply_integrand => oscillation_multiply_oscillation !< Oscillation * oscillation.
    procedure, pass(lhs),  public :: integrand_multiply_real => oscillation_multiply_real             !< Oscillation * real.
    procedure, pass(rhs),  public :: real_multiply_integrand => real_multiply_oscillation             !< Real * Oscillation.
    procedure, pass(lhs),  public :: add => add_oscillation                                           !< Oscillation + Oscillation.
    procedure, pass(lhs),  public :: assign_integrand => oscillation_assign_oscillation               !< Oscillation = Oscillation.
    procedure, pass(lhs),  public :: assign_real => oscillation_assign_real                           !< Oscillation = real.
endtype oscillation

Note that the type bound procedure update_previous_steps has been added. It is necessary inside the time-like loop when using an Adams-Bashforth integrator, e.g.

  do step = 1, num_steps
    if (s>step) then
      ! the time steps from 1 to s - 1 must be computed with other scheme...
      call euler_integrator%integrate(field=attractor, dt=dt)
    else
      call ab_integrator%integrate(field=attractor, steps=s, dt=dt)
    endif
    solution(0, step) = step * dt
    solution(1:space_dimension, step) = attractor%output()
    call attractor%update_previous_steps()
  enddo

Currently I have implemented only up to 3 steps, formal 3rd order accurate.

Please, review my implementation.

See you soon.

[szaghi]

Hi all,

implicit Adams-Moulton and predictor-corrector Adams-Bashforth-Moulton solver are within us!

[szaghi]

Hi all,

explicit embedded Runge-Kutta solvers class is here. Presently, there is only the Dormand and Prince 7 stages, 4th order scheme, but others will come soon.

The implementation of the embedded RK class has required to modify the ADT integrand API: a new operator is necessary, namely the method for the estimation of the local truncation error, see this. The only test I have updated is the oscillation one (note that I have tried to sanitize the tree of directories thus some files are in new places for you): the oscillation type shows a typical operator for estimating the local truncation error (a normalize norm L2 difference error between 2 oscillation fields) that you can use for other tests. For the oscillation test this new solver works incredibly well varying the time step: it obtain essentially the same accuracy of other with the most wide time size admissible (thus with less time step integration). Here is a plot of a result fixing the local truncation error to e-6.

I have a question for you: do how we can measure the accuracy of an adaptive time step schemes such that? The time resolution varies during the integration, thus an asymptotic time resolution reduction is not possible. I can analyze the solution varying the tolerance on the local truncation error that drives the time step adaption. but I am not sure that it is the right way. Any suggestions are very appreciated.

See you soon with CAF parallel test!!!

[szaghi]

Hi all,

some new embedded RK solvers are here. I go up to the Feagin 17 stages solvers, it being a RK of order 10th with error estimator of order 8th.

See you.