fortran-foss-programmers/foodie

Fortran Object-Oriented Differential-equations Integration Environment, FOODIE

FOODIE

FOODIE, Fortran Object-Oriented Differential-equations Integration Environment

• FOODIE is a pure Fortran (KISS) library providing an awesome environment for the numerical integration of Differential-equations (ODE, PDE);
• FOODIE is Fortran 2008+ standard compliant;
• FOODIE is OOP designed;
• FOODIE is a Free, Open Source Project.

Table of Contents

• What is FOODIE?
• Main features
• Status
• Copyrights
• Documentation

What is FOODIE?

Modern Fortran standards (2003+) have introduced support for Object-Oriented Programming (OOP). Exploiting new features like Abstract Data Type (ADT) is now possible to develop a KISS library providing an awesome environment for the numerical integration of Differential-equations such as Ordinary and Partial Differential Equations (ODE, PDE). FOODIE is tailored to the systems arising from the semi-discretization of PDEs, but it is not limited to them.

The FOODIE environment allows the (numerical) solution of general, non linear differential equations system of the form:

where:

• U_t = dU/dt;
• U is the vector of state variables being a function of the time-like independent variable t;
• R is the (vectorial) residual function, it could be a non linear function of the solution U itself;
• F is the (vectorial) initial conditions function.

The FOODIE has two main purposes:

• for developers devising new schemes for the numerical integrations of differential equations (DE):
  • provide a concise, clear, robust and comprehensive abstract environment by means of which:
    • express the solvers formulae with a very high-level language, it being close as much as possible to their natual mathematical formulations; this ensures:
      • clearness, conciseness and robustness;
      • fast-developing;
• for clients that must solve a differential equations system:
  • provide a simple, standard API for many built-in DE solvers out-of-the-box, thus allowing:
    • fast-solution of new problems;
    • robustness: the same DE solver is applied to different problems, i.e. cross-validation;

Go to Top

Main features

FOODIE is aimed to be a KISS-pure-Fortran library for integrating Ordinary Differential Equations (ODE), it being:

• Pure Fortran implementation;
• KISS and user-friendly:
  • simple API, presently based on the Rouson's Abstract Data Type Pattern [8];
  • easy building and porting on heterogeneous architectures;
• comprehensive solvers set out-of-the-box:
  • explicit schemes:
    • Adams-Bashforth schemes see [7, 12]:
      • 1 step, namely the forward explicit Euler scheme, 1st order accurate;
      • 2 to 16 steps, 2nd to 16th accurate, respectively;
    • Euler (forward explicit) scheme, 1st order accurate;
    • Leapfrog, 2nd order accurate:
      • unfiltered leapfrog, 2nd order accurate, mostly unstable, see [4];
      • Robert-Asselin filtered leapfrog, 1st order accurate, see [4, 5, 6];
      • Robert-Asselin-Williams filtered leapfrog, 3rd order accurate, see [5, 6];
    • Linear Multistep Methods, SSP schemes see [16]:
      • 3 steps, 2nd order accurate;
      • 4 steps, 3rd order accurate;
      • 5 steps, 3rd order accurate;
    • Linear Multistep Methods, SSP with Variable Step Size (VSS) schemes see [17]:
      • 2 steps, 2nd order accurate;
      • 3 steps, 2nd order accurate;
      • 3 steps, 3rd order accurate;
      • 4 steps, 3rd order accurate;
      • 5 steps, 3rd order accurate;
    • Linear Multistep Runge-Kutta Methods SSP schemes see [18]:
      • 2 steps, 2 stages, 3rd order accurate;
      • 3 steps, 2 stages, 3rd order accurate;
      • 4 steps, 3 stages, 5th order accurate;
      • 3 steps, 6 stages, 5th order accurate;
      • 3 steps, 5 stages, 6th order accurate;
      • 3 steps, 7 stages, 7th order accurate;
      • 4 steps, 7 stages, 7th order accurate;
      • 4 steps, 5 stages, 8th order accurate;
      • 5 steps, 9 stages, 8th order accurate;
      • 4 steps, 9 stages, 9th order accurate;
      • 3 steps, 20 stages, 10th order accurate;
    • Runge-Kutta schemes:
      • [+] Linear SSP (of any order) schemes, see [16]:
        • generic s-stages of order (s-1)-th;
        • generic s-stages of order s-th;
      • low-storage schemes, see [1, 2, 3]:
        • 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
        • 2 stages;
        • 3 stages;
        • 4 stages;
        • 5 stages, 4th order accurate, 2N registers, see [3];
        • 6 stages, 4th order accurate, 2N registers, see [9];
        • 7 stages, 4th order accurate, 2N registers, see [9];
        • 12 stages, 4th order accurate, 2N registers, see [10];
        • 13 stages, 4th order accurate, 2N registers, see [10];
        • 14 stages, 4th order accurate, 2N registers, see [10];
      • TVD/SSP schemes, see [1]:
        • 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
        • 2 stages, 2nd order accurate;
        • 3 stages, 3rd order accurate;
        • 4 stages;
        • 5 stages, 4th order accurate;
      • embedded (adaptive) schemes:
        • Heun-Euler, 2 stages, 2nd order accurate;
        • Runge-Kutta-Fehlberg, 5 stages, 4th order accurate;
        • Runge-Kutta-Cash-Karp, 6 stages, 5th order accurate, see [13];
        • Prince-Dormand, 7 stages, 4th order accurate, see [11];
        • Calvo, 9 stages, 6th order accurate, see [14];
        • Feagin, 17 stages, 10th order accurate, see [15];
  • implicit schemes:
    • Runge-Kutta schemes;
    • Adams-Moulton schemes:
      • 0 step, 1st order accurate;
      • 1 step, 2nd accurate;
      • 2 steps, 3rd accurate;
      • 3 steps, 4th accurate;
    • Backward Differentiation Formula schemes:
      • 1 step, namely the backward implicit Euler scheme, 1st order accurate;
      • 2 to 6 steps, 2nd to 6th accurate, respectively;
  • predictor-corrector schemes:
    • Adams-Bashforth-Moulton schemes:
      • 1 step, AB(1)-AM(0), 1st order accurate;
      • 2 steps, AB(2)-AM(1), 2nd accurate;
      • 3 steps, AB(3)-AM(2), 3rd accurate;
      • 4 steps, AB(4)-AM(3), 4th accurate;
• efficient and non intrusive:
  • FOODIE environment is unaware of any eventual parallel paradigms the clients used, but it is proved to preserve high scalability on parallel architectures such as:
    • OpenMP directive-based codes on shared memory multi/many cores architectures;
    • CoArray Fortran (CAF) based codes for Partitioned Global Address Space (PGAS) programming model;
    • MPI based code on distributed memory clusters;
    • GPGPU/accelerators device enabled codes;
• Tests-Driven Developed (TDD):
• well documented:
  • clear documentation of schemes implementations, e.g. see Adams-Bashforth API documentation;
  • complete API reference;
  • comprehensive wiki:
• collaborative developed on GitHub;
• FOSS licensed;

Any feature request is welcome.

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] Numerical methods used in atmospheric models, Mesinger F. and A. Arakawa, Global Atmospheric Research Programme (GARP), Technical Report, 1976.

[5] A Proposed Modification to the Robert-Asselin Time Filter, Williams, P. D., Mon. Wea. Rev., vol. 137, pp. 2538--2546, 2009, 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] Scientific Software Design: The Object-Oriented Way, Rouson, Damian and Xia, Jim and Xu, Xiaofeng, 2011, ISBN 9780521888134, Cambridge University Press, New York, NY, USA.

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

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

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

[12] Cowell Type Numerical Integration As Applied to Satellite Orbit Computation, J. L. Maury Jr., G. P. Segal, X-553-69-46, April 1969, NASA-TM-X-63542.

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

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

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

[16] Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations, S. Gottlieb, D. Ketcheson, C.W. Shu, 2011, 978-981-4289-26-9, doi:10.1142/7498, World Scientific Publishing Co. Pte. Ltd.

[17] Strong Stability Preserving Explicit Linear Multistep Methods with Variable Step Size, Y. Hadjimichael, D. Ketcheson, L. Lóczi, A. Németh, 2016, SIAM J. Numer. Anal., 54(5), 2799–2832.DOI:10.1137/15M101717X.

[18] Explicit Strong Stability Preserving Multistep Runge-Kutta Methods, C. Bresten, S. Gottlieb, Z. Grant, D. Higgs, D. Ketcheson, A. Nemeth, Mathematics of Computation, 2016, http://dx.doi.org/10.1090/mcom/3115.

Go to Top

Status

FOODIE project is young, but developed with love. Many integrators have been implemented using the Rouson's Abstract Data Type Pattern and tested with complex problems, but the library API is still in beta testing status. Nevertheless, FOODIE is already proven to be able to integrate a wide range of different ODE problems, from pure ODEs (Lorenz and inertial oscillations equations) to complex PDEs (Burgers and Euler equations), see the documentation.

We are searching for Fortraners enthusiast joining our team!

Copyrights

FOODIE is an open source project, it is distributed under a multi-licensing system:

• for FOSS projects:
  • GPL v3;
• for closed source/commercial projects:
  • BSD 2-Clause;
  • BSD 3-Clause;
  • MIT.

Anyone is interest to use, to develop or to contribute to FOODIE is welcome, feel free to select the license that best matches your soul!

More details can be found on wiki.

Go to Top

Documentation

Besides this README file the FOODIE documentation is contained into its own wiki. Detailed documentation of the API is contained into the GitHub Pages that can also be created locally by means of ford tool.

Go to Top

foodie's Issues

Error in Adams-Bashforth integrator implementation

The integrate method of the adams_bashforth_integrator class is implemented as:

  subroutine integrate(self, field, dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with Adams-Bashforth class scheme.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(adams_bashforth_integrator), intent(IN)    :: self  !< Actual AB integrator.
  class(integrand),                  intent(INOUT) :: field !< Field to be integrated.
  real(R_P),                         intent(in)    :: dt    !< Time step.
  integer(I_P)                                     :: s     !< Steps counter.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  do s=1, self%steps
    field = field + field%t(n=s) * (dt * self%b(s))
  enddo
  call field%update_previous_steps
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate

For simplicity, let's consider the 2nd order Adams-Bashforth (AB2). field%state then has the second dimension of size 2, which corresponds to the number of stored time levels. The solution to the ODE du/dt = F(u,t) using AB2 is:

u(n+1) = u(n) + dt*[-0.5*F(n-1)+1.5*F(n)]

In the code, this is calculated with the loop:

  do s=1, self%steps
    field = field + field%t(n=s) * (dt * self%b(s))
  enddo

which is equivalent to:

field = field + field%t(n=1) * (dt * self%b(1))
field = field + field%t(n=2) * (dt * self%b(2))

The problem that arises here is that field%t(n=1) * (dt * self%b(1)) is added to field using the overloaded operator + which performs element-wise sum for each respective array element. So, tendencies from each time level (n-2, n-1) are always added to the same time level, instead of time level n. The implementation is correct only for number of steps = 1, where the scheme degenerates to forward Euler, and there is only time level n in both the state and the tendency.

Please discuss what would be the best way to solve this. Somehow we need to add the tendency from previous time levels to the field at the present time level, something like:

  do s=1, self%steps
    field(self%steps) = field(self%steps) + field%t(n=s) * (dt * self%b(s))
  enddo

but this is not allowed because integrand is an abstract derived type, and we cannot access its state at this level in the code. Perhaps a solution would be to implement another method, perhaps field%update(), that would act to add tendencies from different time levels to state at the desired time level?

How to implement multiple time level storage in FOODiE integrands?

With hope of including some multilevel integrators into FOODiE, such as leapfrog or Adams-Bashforth, the integrand type needs to have a capability of storing multiple levels. Currently, the integrand type is defined as:

type, abstract :: integrand

  private
  real(R_P), dimension(:), allocatable :: state
  ...
endtype integrand

where the dimension of state represents exactly that - the number of dimensions, that is number of equations in the system. For example, in case of Lorenz equations, the size of state is 3 (see https://github.com/Fortran-FOSS-Programmers/FOODiE/blob/master/src/tests/lorenz/lorenz.f90), and in case of oscillation equation, it is 2 (see https://github.com/Fortran-FOSS-Programmers/FOODiE/blob/master/src/tests/oscillation/oscillation.f90).

In the current implementation of integrand type, its state is stored only at one time level, the present one (say, n). However, in order to implement multi-level schemes, it will need to be able to store states at multiple time levels - n, n-1, n-2... etc., and perhaps even n+1, n+2...etc. for implicit integrators.

I suggest that the time storage be implemented as a second rank in the array of states:

type, abstract :: integrand

  private
  real(R_P), dimension(:,:), allocatable :: state
  ...
endtype integrand

The first index would still reference the equation of the system. The second index would now reference the time level. The time derivative function, for example in the case of oscillation equations:

  pure function dOscillation_dt(self) result(dState_dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Time derivative of Oscillation field.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(oscillation), intent(IN) :: self      !< Oscillation field.
  class(integrand),  allocatable :: dState_dt !< Oscillation field time derivative.
  type(oscillation), allocatable :: delta     !< Delta state used as temporary variables.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! preparing temporary delta
  allocate(delta)
  allocate(delta%state(size(self%state)))
  ! Oscillation equations
  delta%state(1) = -self%f * self%state(2)
  delta%state(2) =  self%f * self%state(1)
  ! hold Oscillation parameters constant over time
  delta%f = 0.
  call move_alloc (delta, dState_dt)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endfunction dOscillation_dt

would need to be re-implemented to accept the time level as an argument, something like this:

  pure function dOscillation_dt(self,n) result(dState_dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Time derivative of Oscillation field.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(oscillation), intent(IN) :: self      !< Oscillation field.
  integer,intent(IN) :: n !< Time level
  class(integrand),  allocatable :: dState_dt !< Oscillation field time derivative.
  type(oscillation), allocatable :: delta     !< Delta state used as temporary variables.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! preparing temporary delta
  allocate(delta)
  allocate(delta%state(size(self%state)))
  ! Oscillation equations
  delta%state(1,n) = -self%f * self%state(2,n)
  delta%state(2,n) =  self%f * self%state(1,n)
  ! hold Oscillation parameters constant over time
  delta%f = 0.
  call move_alloc (delta, dState_dt)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endfunction dOscillation_dt

This issue aims to open discussion about the most sensible way to implement this. Let me know if you have any ideas. I will try to implement some multi-level logic and integrators next week.

Benchmark of ADT-based solvers against procedural ones

Exactly what the title says. I think it will be important to test and benchmark the performance of FOODiE's solvers against their equivalents written in a procedural style, essentially testing the efficiency of abstract data types and overloaded operators. The results will have to be satisfactory in order for FOODiE to be considered in HPC applications. I am not sure yet what satisfactory means in this case but I will think about it.

Comparing the performance of FOODiE between major Fortran compilers will be elucidating as well.

Please discuss. I will be happy to make some tests in the next few weeks.

Tracking download statistics

One thing I think is very important when going after funding is having data on the usage and impact of the software. This could take several forms:

1. Download statistics.
2. A list of installations at major sites, e.g., ask CINECA if they would install it for users.
3. Support letters for grant proposals (if you know users, take the time get to know them and how FOODIE is helping them so that you cultivate the relationships required to solicit such letters and also of course for lots of other reasons).
4. Citations via Google Scholar or a similar tool. This is where it will be very useful to get an initial conference paper in the pipeline early if that gets the work into the literature sooner than a journal article. Let's keep an eye out for workshops Jeff hosts next year. He submits the workshops as proposals to larger conferences so the schedule is only set when the conferences announce acceptance of the proposals. I think his calls for papers usually go out at least 6 months before the conference with submission deadlines at least 3 months before the conference (these are very rough estimates from memory but I don't recall for sure).

Item 1 above is something the most challenging and is the subject of the remainder of this post.

One is that we post our own tar ball with releases, e.g., see "opencoarrays-1.1.1.tar.gz" at https://github.com/sourceryinstitute/opencoarrays/releases/tag/1.1.1. Then we use the following external tool to track downloads of the aforementioned tar ball: http://www.somsubhra.com/github-release-stats/. Enter "sourceryinstitute" and "opencoarrays" into the two fields on the latter page. Sadly even this only provides incomplete data because if someone downloads either of the tar balls that GitHub automatically posts to the aforementioned URL, we get no information on those tar balls. I've even thought about renaming the tar ball "download-this-one.tar.gz" or something similar. Very frustrating.

To help make installation easier for users, I also asked Alessandro to develop a Portfile that enables installation via MacPorts package management software on OS X. That Portfile now exists, but it turns out that MacPorts will only report downloads if the user also installs the mpstats port. Arrggghh. We ask people to do so, but we have no control over it and the corresponding web page shows only one download: http://stats.macports.neverpanic.de/categories/38/ports/25626#installs_over_time, which I'm almost certain is wrong. Our tar ball has been downloaded from GitHub over 275 times. It doesn't seem believable to me that the MacPorts installations would be less than more than 100 times smaller than the tar ball downloads now that we mention MacPorts prominently in our installation instructions.

I hope people can push on the maintainers of GitHub to do a better job with tracking such data, but at least it's good the we have something. That's better than nothing, which is what BitBucket offers.

Also, I have a web developer investigating whether we can use Google Analytics to track clicks on the link to the OpenCoarrays release tar ball from the main page on www.opencoarrays.org, but I don't have an answer on that yet either. I'm really amazed at how difficult it can be to get data that I would think would be easily accessible. Nonetheless, I believe strongly it's worth the effort.

No rule to make target 'pyplot_module.f90'

Hello,

I wanted to give the package a spin and ran into what appears to be a
Makefile problem:

ig25@linux-fd1f:~/Downloads/FOODIE-master> make
make FC=gfortran FCFLAGS="-cpp -g -O0 -C -fbacktrace" --directory=src/lib
make[1]: Entering directory '/home/ig25/Downloads/FOODIE-master/src/lib'
ln -sfv ../../external/pyplot-fortran/src/pyplot_module.f90 .
„./pyplot_module.f90“ -> „../../external/pyplot-fortran/src/pyplot_module.f90“
make foodie
make[2]: Entering directory '/home/ig25/Downloads/FOODIE-master/src/lib'
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_kinds.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_adt_integrand.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_adams_bashforth.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_adams_moulton.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_adams_bashforth_moulton.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_euler_explicit.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_leapfrog.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_low_storage_runge_kutta.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_embedded_runge_kutta.f90
gfortran -c -cpp -g -O0 -C -fbacktrace foodie_integrator_tvd_runge_kutta.f90
make[2]: *** No rule to make target 'pyplot_module.f90', needed by 'pyplot_module.o'. Schluss.
make[2]: Leaving directory '/home/ig25/Downloads/FOODIE-master/src/lib'
Makefile:25: recipe for target 'all' failed
make[1]: *** [all] Error 2
make[1]: Leaving directory '/home/ig25/Downloads/FOODIE-master/src/lib'
Makefile:35: recipe for target 'foodie' failed
make: *** [foodie] Error 2

Did I miss anything? Is there something wrong with my environment?

This is on x86_64-pc-linux-gnu using a patched gfortran.

Move toward OpenFoam...

From a great Idea of Fra...

How we can implement an intermediate abstract type that solves some tasks that are COMMON to many U=R(U)...

Benchmarks

Here we can discuss the benchmarks analysis.

It could be useful to compare with the performance of other non abstract implementations.

Architecture	Comparison with other non ACP	Conclusion
Serial
Parallel shared memory
Parallel distributed memory

OO Design of the time_derivative method

Hic sunt leones

The OOP abstract calculus as the one of Rouson already implemented and tested with the Lorenz equations example could be limiting or problematic...

As a matter of facts, in the Lorenz equations the time derivative (hereafter the residuals) function is very simple to compute. In a more complex scenario of non-linear PDE system the residuals is a very complex function (involving computations on large stencils). Generally, solving a PDE means applay a spatial operator that builds the residuals that are the advanced in time by an ODE solver like FOODiE.

Consequently, the design of residuals function is crucial. Moreover, for multi-steps/multi-stages like RK ones the residuals function is used many times for each single time step. This involves also a smart boundary conditions handling.

In my previous experience, the residuals function is a procedure operating on a whole set of integrand field (a block of finite volumes being self-consistent for each time step, i.e. boundary conditions, space operator...). Now (in the present integrand definition), the residual function seems to be designed to operate on a single-valued field. Obviously, one can assume that the Lorenz field example is just the time integration of a single space location and applay it to tge whole domain. For the Lorenz example this is simple, but for non linear PDE each space lication residuals depend on other space locations...

In such a case, one solution could be to reduce the integrand time derivative class function to a simple lookup array where the residuals are computed previously in the whole domain taking into account the non linearity. This could work for a simple scheme like the explicit Euler one, but I think that it not feasible for more comple schemes like the RK ones.

Any suggestions?

I will try study better the Burgers solution of Rouson it being a more complex example than the Lorenz one.

Enabling FOODIE compilation on Fermi

Currently, IBM XL Fortran complains with a lot of issues...

I am solving it right now.

Auxiliary methods for integrators

Dear all,

some of you (@jacobwilliams and @zbeekman) have expressed the idea that the FOODiE's integrators should provide some auxiliary methods:

It would be wonderful to have an object oriented library of time integrators that could be used, for integration of abstract data type (ADT) calculus (See Rouson et al.) with procedures to query their properties, like their stability locus plots sigma roots, etc. Users could then be guided to picking a good time integrator for their PDE/ODE of interest and it would take a lot of the headache out of picking a suitable time integration scheme to match your spatial scheme when making a PDE solver.

I am now quite confident with the ADT of Rouson (Euler and explicit TVD RK schemes have been implemented) thus, in the meanwhile other schemes are under development, I would like to implement some of the auxiliary methods that you asked for. Here is list tracking issue:

• stability locus plot;
• event finding... @jacobwilliams what do you mean?
• supported stages/steps inquiring methods, e.g. tvd_rk%is_supported(s=7);
  • lower and upper bound limits of supported stages/steps, e.g. tvd_rk%max_stages();
• provide a descriptive summary (varying string?) of the solver or of the class of solvers, e.g. tvd_rk%description();
• ...

update compile options

The currently compile steps do not reflect the current state of the project:

FoBiS.py build -lmodes does not result in:

The fobos file defines the following modes:
  - "shared-gnu"
  - "shared-intel"
  - "static-gnu"
  - "static-intel"
  - "tests-gnu"
  - "tests-intel"
  - "tests-gnu-debug"
  - "tests-intel-debug"

as stated in the docs. On openSUSE 42.1 with gfortran 5.3.1, I get:

The fobos file defines the following modes:
  - "foodie-static-opencoarrays-gnu" Build FOODIE library in static release mode with OpenCoarrays-GNU gfortran
  - "foodie-static-opencoarrays-gnu-debug" Build FOODIE library in static release mode with OpenCoarrays-GNU gfortran
  - "foodie-static-gnu" Build FOODIE library in static release mode with GNU gfortran
  - "foodie-static-gnu-debug" Build FOODIE library in static debug mode with GNU gfortran
  - "foodie-shared-gnu" Build FOODIE library in shared release mode with GNU gfortran
  - "foodie-shared-gnu-debug" Build FOODIE library in shared debug mode with GNU gfortran
  - "foodie-static-intel" Build FOODIE library in static release mode with Intel Fortran
  - "foodie-static-intel-debug" Build FOODIE library in static debug mode with Intel Fortran
  - "foodie-shared-intel" Build FOODIE library in shared release mode with Intel Fortran
  - "foodie-shared-intel-debug" Build FOODIE library in shared debug mode with Intel Fortran
  - "foodie-static-gnu-coverage" Build FOODIE library in static release coverage mode with GNU gfortran
  - "foodie-static-gnu-pure" Build FOODIE library in static release mode with GNU gfortran
  - "foodie-static-gnu-debug-pure" Build FOODIE library in static debug mode with GNU gfortran
  - "foodie-static-gnu-coverage-pure" Build FOODIE library in static release coverage mode with GNU gfortran
  - "foodie-static-intel-pure" Build FOODIE library in static release mode with Intel Fortran
  - "foodie-static-intel-debug-pure" Build FOODIE library in static debug mode with Intel Fortran

Should the docs be updated to reflect this or is it compiler specific? I updated the wiki myself, but further explanation from a maintainer is probably helpful

Collaborative workflow: which model(s) you like?

Dear all,
a nice discussion with @zbeekman and @rouson is started here about which workflow model adopt to develop FOODiE. I created this room for not polluting the other, this topic being interesting too.

Workflow, what?

The topic simply concerns the approach we will use to develop FOODiE in a collaborative way, obviously using Git and GitHub. Presently, FOODiE has an official GitHub repository owned by our group where I pushed my commits, but all of you (I hope) will push yours, thus it could be useful to agree on a common workflow.

My approach

Currently, I use the Driessen's branching model (by means of gitflow help) for my local developments and pushing only the master branch on the FOODiE GitHub repository. Essentially, my local workflow is based on a hopefully sane branching approach: I work on the develop branch where I create new branch for each new feature; when a feature is correctly implemented its branch is merged into the develop one; when a set of features have been implemented into the develop branch it is merged into the master and finally pushed to GitHub remote. Besides, I use branch also for hot-fixes and releases.

Zaak's approach

Zaak, if I understand right, prefers the forking that should be more light than the previous one. There is no formal branching naming model for driving the implementation of new feature or fixing bugs, etc... In few words, just fork, branch as you like, push your branches and merge pull requests (maybe I am wrong, I have not deeply read the forking model documentation).

My first thought

The gitflow and forking models are not in contrast: I think we can use both at the same time. I suppose I can (locally) continue to chain up my chaotic development approach on a sane gitflow branching model, while I can (remotely) still collaborate with you by means of a more light and flexible forking model.

Well known workflow models

Besides the Zaak's consideration and mine, it can be useful to review the models that others smart developers use. From this atlassian tutorial, well known workflow models are:

• centralized;
• feature-branch;
• gitflow;
• forking;
• suggest others.

The vote is started

Post your detailed opinions and your current workflow, please :-)

Coarray enabling

Hi all,

I have found 1 hour to play with CAF (note that FoBiS now supports a new compiler OpenCoarrays-GNU gfortran). I faced with some problems and reading the work of @rouson I think that the current integrand ADT API must be slightly modified, but I appreciate your opinions.

I think that the best way (for reasons that I will explain later) is to let concrete extensions of integrand to embedded CAF members. However, it seems that a valid extension of an abstract type is allowed to contain a coarray only if the abstract parent does too. For example you can see the @rouson co_object and its concrete extensiin global_field implemented for his Burgers CAF test. As a consequence, I propose to change our integrand with this:

type, abstract :: integrand
#ifdef CAF
class(*), allocatable :: dummy_to_allow_extension[:]
#endif
  !< Abstract type for building FOODIE ODE integrators.
  contains
    ! public deferred procedures that concrete integrand-field must implement
    procedure(time_derivative),      pass(self), deferred, public :: t !< Time derivative, residuals.
    ! operators
    ....
endtype integrand

If I understand right, this should allow a concrete extension of the type

type, extends(integrand) :: integrand_concrete
#ifdef CAF
  real, allocatable :: U(:)[:]
#else
  real, allocatable :: U(:)
#endif
  contains
    ...
endtype integrand_concrete

Why embedd Coarrays?

Because solving PDE on multi-block split domain involves global comunication of boundaries values for example. Consequently, when FOODIE solvers invoke %t() a global comunication is necessary.

What do you think?

Implementation of low-storage Runge-Kutta class of solvers

This is just the roomwhere we discuss the implementation of low-storage RK class. This is in my radar for the next week.

is a "foodie solver constructor" class desirable?

Following a thread on google group CLF, I realized that we have never talked about the possibility to add a solver constructor class in order to provide a unique foodie solver class that is then specialized into the concrete solvers provided by means of the factory pattern (I hope to have used the tags correctly in the OOP slang), e.g.

use foodie, only : foodie_emd_rk_constructor, foodie_factory, foodie_integrator

type(foodie_factory) :: factory
class(foodie_integrator), allocatable :: integrator

call factory%create(constructor=foodie_emd_rk_constructor(S=10), integrator=integrator)

do 
  call integrator%integrate(...)
enddo

What do you think about this?

P.S. I have hopefully destroyed the Rock fallen on my shoulders in the two last weeks, I hope tomorrow I will come back 😄

Announcement Paper

As anticipated on gitter, I am planning to write an announcement paper to be published on CPC or similar.

I have just uploaded a skeleton tex with the list of the authors that (for the moment) I would like to include into the paper. No matter the order of names, we can modify it in any order.

I hope others will join us.

See you soon.

Implementation of Leapfrog class of solvers

This is just the roomwhere we discuss the implementation of Leapfrog class. This is in my radar for the next 2 weeks.

foodie_integrator_leapfrog.f90 does not compile with ifort 15.0.0

foodie_integrator_leapfrog.f90 fails to compile with Intel Fortran compiler 15.0.0 with an internal compiler error:

$ make foodie
make FC=ifort FCFLAGS="-cpp -g -O0 -C -traceback -assume realloc_lhs" --directory=src/lib
make[1]: Entering directory `/home/milan/FOODiE/src/lib'
ln -sfv ../../external/pyplot-fortran/src/pyplot_module.f90 .
‘./pyplot_module.f90’ -> ‘../../external/pyplot-fortran/src/pyplot_module.f90’
make foodie
make[2]: Entering directory `/home/milan/FOODiE/src/lib'
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_kinds.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs type_integrand.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_integrator_adams_bashforth.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_integrator_adams_moulton.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_integrator_adams_bashforth_moulton.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_integrator_euler_explicit.f90
ifort -c -cpp -g -O0 -C -traceback -assume realloc_lhs foodie_integrator_leapfrog.f90
foodie_integrator_leapfrog.f90(100): catastrophic error: **Internal compiler error: internal abort** Please report this error along with the circumstances in which it occurred in a Software Problem Report.  Note: File and line given may not be explicit cause of this error.
compilation aborted for foodie_integrator_leapfrog.f90 (code 1)

Can anybody try to reproduce this with a later version? I am on a x86-64 linux.

How to parallelize (OpenMP) procedural version of Euler 1D test

Hi all,

answering to a question of @rouson I have noted that for the procedural version of the OpenMP Euler 1D test I have missed to parallalize an important loop. See the following

Procedural integrate

This is the hard-coded-Euler-1D RK integrator

  subroutine integrate_rk(self, U, stage, Dt, t)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with explicit TVD (or SSP) Runge-Kutta scheme.
  !<
  !< @note This method can be used **after** the integrator is created (i.e. the RK coeficients are initialized).
  !---------------------------------------------------------------------------------------------------------------------------------
  class(tvd_runge_kutta_integrator), intent(IN)    :: self      !< Actual RK integrator.
  class(euler_1D_openmp),            intent(INOUT) :: U         !< Field to be integrated.
  class(euler_1D_openmp),            intent(INOUT) :: stage(1:) !< Runge-Kutta stages [1:stages].
  real(R_P),                         intent(IN)    :: Dt        !< Time step.
  real(R_P),                         intent(IN)    :: t         !< Time.
  integer(I_P)                                     :: s         !< First stages counter.
  integer(I_P)                                     :: ss        !< Second stages counter.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! computing stages
  do s=1, self%stages
    stage(s) = U ! parallelized
    do ss=1, s - 1
      stage(s)%U = stage(s)%U + stage(ss)%U * (Dt * self%alph(s, ss))
    enddo
    stage(s) = stage(s)%t()  ! parallelized
  enddo
  ! computing new time step
  do s=1, self%stages
    U%U = U%U + stage(s)%U * (Dt * self%beta(s))
  enddo
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate_rk

The stages computing and Euler field new time step loops are not parallelized. Only the defined assignments and the %t() operator are parallel. Now consider the present FOODIE equivalent integrator.

FOODIE integrate

  subroutine integrate(self, U, stage, Dt, t)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with explicit TVD (or SSP) Runge-Kutta scheme.
  !<
  !< @note This method can be used **after** the integrator is created (i.e. the RK coefficients are initialized).
  !---------------------------------------------------------------------------------------------------------------------------------
  class(tvd_runge_kutta_integrator), intent(IN)    :: self      !< Actual RK integrator.
  class(integrand),                  intent(INOUT) :: U         !< Field to be integrated.
  class(integrand),                  intent(INOUT) :: stage(1:) !< Runge-Kutta stages [1:stages].
  real(R_P),                         intent(IN)    :: Dt        !< Time step.
  real(R_P),                         intent(IN)    :: t         !< Time.
  integer(I_P)                                     :: s         !< First stages counter.
  integer(I_P)                                     :: ss        !< Second stages counter.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  select type(stage)
  class is(integrand)
    ! computing stages
    do s=1, self%stages
      stage(s) = U
      do ss=1, s - 1
        stage(s) = stage(s) + stage(ss) * (Dt * self%alph(s, ss))
      enddo
      stage(s) = stage(s)%t(t=t + self%gamm(s) * Dt)
    enddo
    ! computing new time step
    do s=1, self%stages
      U = U +  stage(s) * (Dt * self%beta(s))
    enddo
  endselect
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate

Again there is not explicit parallelism. However, in this FOODIE integrator all part are parallelized because the defined operators + * = are parallel.

Considerations

I think this is the reason why the FOODIE-aware version scales better (I am not considering the absolute cpu time, I am considering only the parallel scaling): in the FOODIE-aware version all ingredients are easily parallelizable (in this sense, having no side effect helps a lot) and they are, whereas in the procedural version it is not trivial (for me) how to parallelize all part (not only the assignments that are already parallelized).

At a first superficial thought, I can modify the procedural version with this:

  !$OMP PARALLEL DO PRIVATE(s, ss) SHARED(U, stage, Dt, self)
  do s=1, self%stages
    stage(s) = U ! parallelized
    do ss=1, s - 1
      stage(s)%U = stage(s)%U + stage(ss)%U * (Dt * self%alph(s, ss))
    enddo
    stage(s) = stage(s)%t()  ! parallelized
  enddo
  ! computing new time step
  !$OMP PARALLEL DO PRIVATE(s) SHARED(U, stage, Dt, self)
  do s=1, self%stages
    U%U = U%U + stage(s)%U * (Dt * self%beta(s))
  enddo

I have not tried, but I guess that this modification does not provide much more parallelism (the stages are only 5...). I can try to tune with collapse or other workarounds, but I feel this is not the best approach.

The question (to specifically @francescosalvadore @muellermichel but also to all): how do you parallelize this procedural integrator? Do I have to restore all the defined operators and parallelized them? In this sense, the abstract FOODIE approach seems to be more natural. Do I have to try to parallelize the inner loops (%U is an array)?

Will spatial differencing algorithms be part of FOODiE?

What the title says - will spatial differencing schemes be part of FOODiE? For example centered, upstream schemes of various orders, MPDATA family of schemes, etc.?

Pure procedures with intent(inout) dummy arguments

Compiling FOODiE with Cray compiler on an XC-30 reveals the following (possible) bug:

ftn -c  type_integrand.f90

module type_integrand
       ^              
ftn-855 crayftn: ERROR TYPE_INTEGRAND, File = type_integrand.f90, Line = 2, Column = 8 
  The compiler has detected errors in module "TYPE_INTEGRAND".  No module information file will be created for this module.

  pure subroutine assignment_integrand(lhs, rhs)
                                       ^         
ftn-1265 crayftn: ERROR TIME_DERIVATIVE, File = type_integrand.f90, Line = 82, Column = 40 
  Non-pointer dummy argument "LHS" to PURE FUNCTION "TIME_DERIVATIVE" must have INTENT(IN) specified for it.

  pure subroutine assignment_real(lhs, rhs)
                                  ^         
ftn-1265 crayftn: ERROR TIME_DERIVATIVE, File = type_integrand.f90, Line = 92, Column = 35 
  Non-pointer dummy argument "LHS" to PURE FUNCTION "TIME_DERIVATIVE" must have INTENT(IN) specified for it.

Cray Fortran : Version 8.3.11 (u83061f83225i83184p83333a83009e83011z83333)
Cray Fortran :                (x8321r83015w83011t8311b83039)
Cray Fortran : Thu Aug 13, 2015  22:37:03
Cray Fortran : Compile time:  0.0080 seconds
Cray Fortran : 103 source lines
Cray Fortran : 3 errors, 0 warnings, 0 other messages, 0 ansi
Cray Fortran : "explain ftn-message number" gives more information about each message.
make: *** [type_integrand.o] Error 1

Apparently, Cray does not like this. Intel and GNU are fine with it. I do not know from the top of my head whether this is allowed by the standard, and it is too late in the evening for me to go look read through the fineprint.

The code compiles without error if the pure attribute is removed for the following procedures: function time_derivative, subroutine assignment_integrand, subroutine assignment_real.

No rule to make target 'foodie_adt_integrand.f90'

Hey, I wanted to say first this project looks super cool and I hope I can get it to work!

I just cloned the repository to my local machine through:
git clone https://github.com/Fortran-FOSS-Programmers/FOODIE

and then after reading the readme here, I tried to:
make foodie

but received the following output and error:

Wills-Freaking-MacBook-Pro:FOODIE willsharpless$ make foodie
Building foodie
/Library/Developer/CommandLineTools/usr/bin/make FC=bgxlf2008_r FCFLAGS="-qsuffix=cpp=f90" --directory=src/lib
/Library/Developer/CommandLineTools/usr/bin/make foodie
make[2]: *** No rule to make target `foodie_adt_integrand.f90', needed by `foodie_adt_integrand.o'.  Stop.
make[1]: *** [all] Error 2
make: *** [foodie] Error 2

I tried running make all but also got the exact same output. Did I miss something in the clone?

Note, I have a Mac with macOS Catalina 10.15.7.

HPC resources: proposal for a grant on European PRACE network

@milancurcic @rouson @zbeekman and all

Dear All,
today I have devoted my break to write a proposal for obtaining HPC resources on the European PRACE network. In few words, PRACE gathers some HPC facilities in the Europe supporting research projects.

The proposal I have written is of Preparatory Access type B type. Such a type is essentially devoted to assessing the performance of codes in order to certify that those codes are allowed for a full PRACE access grant. In few words, the preparatory access grants allow developers to test their codes on HPC facilities. In particular, the type B are grants that do not require support from PRACE people. Consequently, such a grant type is ideal for FOODIE: we just need to assess its performance into a real HPC scenario. Moreover, the limited resources of such a grant should give us more hope that the grant will be obtained. Here at INSEAN we have already obtained one of this grant for another project.

If the grant will be provided, we will obtain access to a small/medium HPC facility. I have requested to use FERMI and Curie thin clusters. Such resources should be allow us to test FOODIE up to thousands of cores in heterogeneous scenario (Fermi is a IBM-BlueGene/Q machine while Curie thin is a Intel® Sandy Bridge EP one).

Attached here there is the draft of our proposal. I would like to have you as collaborators (even if I think this is not strictly required for success of the proposal): if you accept to be a collaborator you will not be asked to do anything, but you will have the possibility to directly access to the clusters for checking/correcting my errors/horrors :-)

In the case you agree to be a collaborator, please check the information I have added about your CV. I think that Milan's data should be correct, while the Damian's and Izaak's ones are very doubtful: in particular I am not sure about your affiliations (for Izaak I have added the University of Maryland, but Izaak has so many others that I am confused about Princeton/Sourcery Institute). The same is for Damian (Sourcery/Stanford). In particular for the proposal I need to know your birth date.

The last deadline for this year is 1st December. If accepted, the grant should start from 1st January 2016.

Let me know what you think about this proposal.

See you soon.

P.S. all other members of the group interested on this are welcome to join us!

PRACE PR - Export on 2015-11-25 at 14.15.22.pdf

Accuracy analysis of the results of tests

Here we can discuss our conclusions about the accuracy of the results obtained from each test.

Passed Tests

List of tests for which the results accuracy is as expected for the formal order of the numerical scheme employed:

• Lorenz' test:
  • explicit:
    • Adams-Bashforth;
    • Euler;
    • Leapfrog;
    • Runge-Kutta:
      • TVD/SSP;
      • Low Storage;
• Oscillation test¹:
  • explicit:
    • Adams-Bashforth;
    • Adams-Bashforth-Moulton;
    • Euler;
    • Leapfrog;
    • Runge-Kutta:
      • TVD/SSP (5 stages solver seems less acurate than the similar low storage 5 stages);
      • Low Storage;
      • embedded³;
• Burgers' test:
  • explicit:
    • Adams-Bashforth;
    • Euler;
    • Leapfrog;
    • Runge-Kutta:
      • TVD/SSP;
      • Low Storage;
• 1D Euler test²:
  • explicit:
    • Adams-Bashforth;
    • Euler;
    • Leapfrog;
    • Runge-Kutta:
      • TVD/SSP;
      • Low Storage;

Issues

1. Analysis completed: all solvers works as expected.
2. Analysis not completed.
3. It is necessary to define a methodology for measuring the accuracy of adaptive time step schemes.

Tests facility

I have added the Lorenz test of Rouson's book. Here we create a list of desired test

• pure ODEs:
  • Lorenz equations, see wiki documentation, test sources;
  • Oscillation equations, see test sources;
• PDEs:
  • Burgers equation, see test sources;
  • 1D Euler equations;
  • 1D Euler equations solved with OpenMP enabled code;

Support for implicit methods

FOODiE currently supports only explicit integration methods, i.e. methods that compute the solution at a future time based on the state at present and/or past times. Implicit methods require a guess of the state at a future time to obtain a more accurate solution. They are useful for some families of problems and in some cases may allow unconditionally stable integration with arbitrarily large time steps, at the expense of accuracy.

I believe that it is essential for FOODiE to include these methods, such as 1st order backward Euler, 2nd order trapezoidal, Adams-Moulton family of methods, Adams-Bashforth-Moulton etc.

Implicit solvers require evaluating the residual function R at time step n+1. This is usually done by a combination of initial guess and iterative procedure like Newton-Raphson iteration or similar. FOODiE should provide such methods for obtaining R at time step n+1. Once they are in place, implementation of the above-mentioned implicit methods should be trivial.

This issue aims to spark discussion on:

• How to implement evaluation of R(n+1) in FOODiE?
• What implicit methods do we want?
• Any other implementation / API suggestions.

Problem with operators overloading

Hi all,

implementing the Burgers' test, I faced with some problems with operators overloading.

Let us consider type_integrand, its definition is

type, abstract :: integrand
  !< Abstract type for building FOODiE ODE integrators.
  contains
    ! public deferred procedures that concrete interpolators must implement
    procedure(time_derivative),      pass(self), deferred, public :: t                       !< Time derivative, residuals function.
    procedure(integrand_op_real),    pass(lhs),  deferred, public :: integrand_multiply_real !< Integrand * real operator.
    procedure(real_op_integrand),    pass(rhs),  deferred, public :: real_multiply_integrand !< Real * integrand operator.
    procedure(symmetric_operator),   pass(lhs),  deferred, public :: add                     !< Integrand + integrand oprator.
    procedure(assignment_integrand), pass(lhs),  deferred, public :: assign_integrand        !< Integrand = integrand.
    procedure(assignment_real),      pass(lhs),  deferred, public :: assign_real             !< Integrand = real.
    ! operators overloading
    generic, public :: operator(+) => add                                              !< Overloading + operator.
    generic, public :: operator(*) => real_multiply_integrand, integrand_multiply_real !< Overloading * operator.
    generic, public :: assignment(=) => assign_integrand, assign_real                  !< Overloading = assignament.
endtype integrand

Now consider the Lorenz concrete extension of the above abstract type

type, extends(integrand) :: lorenz
  !< Lorenz equations field.
  !<
  !< It is a FOODiE integrand class.
  private
  real(R_P), dimension(:), allocatable :: state        !< Solution vector.
  real(R_P)                            :: sigma=0._R_P !< Lorenz \(\sigma\).
  real(R_P)                            :: rho=0._R_P   !< Lorenz \(\rho\).
  real(R_P)                            :: beta=0._R_P  !< Lorenz \(\beta\).
  contains
    procedure, pass(self), public :: output                                          !< Extract Lorenz field.
    procedure, pass(self), public :: t => dLorenz_dt                                 !< Time derivate, resiuduals function.
    procedure, pass(lhs),  public :: integrand_multiply_real => lorenz_multiply_real !< lorenz * real operator.
    procedure, pass(rhs),  public :: real_multiply_integrand => real_multiply_lorenz !< Real * Lorenz operator.
    procedure, pass(lhs),  public :: add => add_lorenz                               !< Lorenz + Lorenz oprator.
    procedure, pass(lhs),  public :: assign_integrand => lorenz_assign_lorenz        !< Lorenz = Lorenz.
    procedure, pass(lhs),  public :: assign_real => lorenz_assign_real               !< Lorenz = real.
endtype lorenz

The overloaded operator * should have both lorenz * real_scalar and real_scalar * lorenz, but only the first work. For example Euler integrator does not compile if I change it

working

  subroutine integrate(field, dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with explicit Euler scheme, 1st order.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(integrand), intent(INOUT) :: field !< Field to be integrated.
  real(R_P),        intent(in)    :: dt    !< Time step.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  field = field + field%t()*dt
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate

not working

changing the position of dt

  subroutine integrate(field, dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate field with explicit Euler scheme, 1st order.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(integrand), intent(INOUT) :: field !< Field to be integrated.
  real(R_P),        intent(in)    :: dt    !< Time step.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  field = field + dt*field%t()
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine integrate

This becomes even more problematic if the high level equation implementation contains a more complex operators combination. For example for the Burgers residual function (involving the space operator) I have had to split the computation into more than on statement (see the delta computation):

  pure function dBurgers_dt(self) result(dState_dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Time derivative of Burgers field, residuals function.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(burgers), intent(IN)    :: self      !< Burgers field.
  class(integrand), allocatable :: dState_dt !< Burgers field time derivative.
  type(burgers),    allocatable :: delta     !< Delta state used as temporary variables.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! preparing temporary variables
  allocate(delta, source=self)
  ! Burgers residuals
  delta = self%xx() * self%nu
  delta = delta - self * self%x()
  call move_alloc (delta, dState_dt)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endfunction dBurgers_dt

Do you have any idea?

P.S. Intel Fortran seems to be more sensible to this problem...

Particle tracer

Dear All,

Can you help me understand if this software is capable of tracing particles in Velocity vector fields? Equation that describes this process is dx_ij / dt = V_i. i = 1,2,3 (3 dimensions), j = 1,..,n (n particles). V_i \in [x0,x1]x[y0,y1]x[z0,z1]. x_ij(0) is given.

`delta%state` initialization missing in dUdt function implementations

Hi Stefano,

I noticed that in the multi-level implementation of R(u,t) function implementations, you did not initialize the state after allocating the array, for example in case of Lorenz:

  pure function dLorenz_dt(self, n) result(dState_dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Time derivative of Lorenz field.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(lorenz),          intent(IN) :: self      !< Lorenz field.
  integer(I_P), optional, intent(IN) :: n         !< Time level.
  class(integrand), allocatable      :: dState_dt !< Lorenz field time derivative.
  type(lorenz),     allocatable      :: delta     !< Delta state used as temporary variables.
  integer                            :: dn        !< Time level, dummy variable.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! preparing temporary delta
  allocate(delta)
  delta%dims = self%dims
  delta%steps = self%steps
  allocate(delta%state(1:self%dims, 1:self%steps))
  !delta%state = 0_R_P ! <<< delta%state should be initialized here
  delta%sigma = self%sigma
  delta%rho = self%rho
  delta%beta = self%beta
  ! Lorenz equations
  dn = self%steps ; if (present(n)) dn = n
  delta%state(1, dn) = self%sigma * (self%state(2, dn) - self%state(1, dn))
  delta%state(2, dn) = self%state(1, dn) * (self%rho - self%state(3, dn)) - self%state(2, dn)
  delta%state(3, dn) = self%state(1, dn) * self%state(2, dn) - self%beta * self%state(3, dn)
  call move_alloc(delta, dState_dt)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endfunction dLorenz_dt

Same is true for the oscillation equations, and I find that initializing to zero now changes the solution from gibberish to just moderately wrong. For example, see the two examples for oscillation equation, first without initialization (notice the y-axis magnitude):

then with initialization of delta%state = 0_R_P:

It is still difficult to test this correctly because of the issue #18.

If you agree with me that delta%state should be initialized to 0 after allocation, I can apply these changes in the current tests that we have.

ODE schemes: a list of schemes to be implemented

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

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

Higher order Adams formula

If there is any interest, I have higher order Adams formulas coded up here.

Object Oriented Design: which pattern(s)?

Hi all,

I have started this new exciting (for me) project.

I started studying the work of Jacob. Then, after reading the Rouson's book, I decided to try a more formal approach.

Presently I have implemented the Abstract Calculus Pattern providing also the same example of Rouson's book the integration of Lorenz equations by means of explicit one step first order Euler scheme.

The pattern is very simple

the abstract integrand

This define the abstract methods necessary for build a FOODiE integrator.

type, abstract :: integrand
  !< Abstract type for building FOODiE ODE integrators.
  contains
    procedure(time_derivative     ), pass(self), deferred :: t
    procedure(asymmetric_operator ), pass(lhs),  deferred :: multiply
    procedure(symmetric_operator  ), pass(lhs),  deferred :: add
    procedure(symmetric_assignment), pass(lhs),  deferred :: assign
    generic :: operator(+) => add
    generic :: operator(*) => multiply
    generic :: assignment(=) => assign
endtype integrand

a FOODiE integrator

Defines in very high level language (by means of operators overloading) the ODE integration scheme
e.g. the explicit Euler is simply

subroutine euler_explicit(field, dt)
  !---------------------------------------------------------------------------------------------------------------------------------
  !< Integrate model with explicit Euler scheme, 1st order.
  !---------------------------------------------------------------------------------------------------------------------------------
  class(integrand), intent(INOUT) :: field !< Field to be integrated.
  real(R_P),        intent(in)    :: dt    !< Time step.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  field = field + field%t()*dt
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine

client side code

The user must extend the abstract integrand defining his/her own equations system. In the case of Lorenz equation it is simply

type, extends(integrand) :: lorenz
  private
  real(R_P), dimension(:), allocatable :: state            !< Solution vector.
  real(R_P)                            :: sigma, rho, beta !< Lorenz parameters.
contains
   procedure, pass(self), public :: t        => dLorenz_dt
   procedure, pass(lhs),  public :: add      => add_lorenz
   procedure, pass(lhs),  public :: multiply => multiply_lorenz
   procedure, pass(lhs),  public :: assign   => assign_lorenz
   procedure, pass(self), public :: output
end type

More details here

First considerations

The approach seems work well. For the lorenz test I have added the generation of a plot of the solution by means of the great pyplot-fortran of Jacob

My main concerns are now related to the efficiency of the approach, but this question will be answered in the future (or maybe @rouson can replay having deeply used such approaches into HPC codes).

How to use HPC resources

Dear All,

it is time to decide how to use our grant.

As I said we have

• 250.000 core-hours on Fermi, CINECA, Italy, a BlueGene/Q system
• 100.000 core-hours on MareNostrum, BSC, Spain, a Intel Sandy Bridge cluster

We will have 6 months to use these resources (I think starting from the date when we activate our accounts).

I am preparing a 2D (hopefully 3D) FOODIE test that I hope will be able to exploit the above resources. I have not a direct experience on MareNostrum, whereas I have used Fermi: the minimum-size job is about thousand of cores (in debug-smallpar queues, that are suitable for our scaling tests, should be smaller), thus the tests could not be too small. In particular we 1GB/core of RAM on Fermi and [2,4,8]GB/core on MareNostrum.

I would like to know if you (all of you, not only the ones directly involved into the paper for example) have some proposal on how to use these resources at best.

In the following a brief summary of the tests I am preparing is reported, I am sorry, but my background is gas-dynamic, I am essentially an inviscid (sub)man 😄

2D Riemann Problem(s)

Classical tests, as reported in the interesting Kurganov and Tadmor paper, e.g.

Shock Diffraction

Very classical test, e.g. Mach 2 shock diffraction on 90° corner, e.g.

Wind Tunnel Step

Even more classical ( 😄 ), the Mach 3 flow into a stepped wind tunnel, e.g.

Unnecessary Finalization

I notice that all of the finalization subroutines in the integrator types are unnecessary. They just set constants to zero or deallocate allocatables. While having an unnecessary finalization subroutine isn't a bit problem in and of itself, the version of gfortran (4.8.4) on the Ubuntu computer in my office does't support finalization and therefore won't compile. Ideally I would just install gfortran5 from the PPA, but the university doesn't trust us mere scientists with administrative rights.

In any case, I've commented out the finalization bindings as a workaround, but given that they aren't necessary and impede compiler support, I would suggest that they shouldn't be present at all. At the very least it would be good if you could add some preprocessor flags around them.

Is necessary select type on dummies class?

Like this

https://github.com/Fortran-FOSS-Programmers/FOODIE/blob/master/src/lib/foodie_integrator_low_storage_runge_kutta.f90#L330

Implementation of 1D Euler's PDEs system test

Here we can discuss the implementation of the 1D Euler's PDEs system test. This is a complex PDEs system (non linear, unstedy, hyperbolic, discontinuos solutions admitting...). It is a good candidate for the benchmark analysis.

I have already implemented it, but not uploaded due to some minor bugs. I should upload it soon in the next week.

Quality of algorithms

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.