1. Introduction

2. Mathematical and Numerical Models

2.3. The explicit low storage Runge-Kutta class of schemes

As aforementioned, standard Runge-Kutta schemes have the drawback to require $N_S$ auxiliary memory registers to store the necessary stages data. In order to make an efficient use of the available limited computer memory, the class of low storage Runge-Kutta scheme was devised. Essentially, the standard Runge-Kutta class (under some specific conditions) can be reformulated allowing a more efficient memory management. Currently FOODIE provides a class of 2N registers storage Runge-Kutta schemes, meaning that the storage of all stages requires only 2 registers of memory with a word length N (namely the length of the state vector) in contrast to the standard formulation where $N_s$ registers of the same length N are required. This is a dramatic improvement of memory efficiency especially for schemes using a high number of stages ($N_s\ge 4$) where the memory necessary is an half with respect the original formulation. Unfortunately, not all standard Runge-Kutta schemes can be reformulated as a low storage one.

Following the Williamson’s approach (see [13,14,15,16]) the standard coefficients are reformulated to the coefficients vectors A, B and C and the Runge-Kutta algorithm becomes:

$$\begin{array}{c}{K}_1=U\left(t^n\right)\\ K_2=0\\ \left.\begin{array}{c}{K}_2=A_s{K}_2+\Delta t·R\left(t^n+C_s\Delta t,K_1\right)\\ K_1=K_1+B_s{K}_2\end{array}\right\}s=1,2,...N_s\\ U\left(t^{n+1}\right)=K_1\end{array}$$

(5)

Currently FOODIE provides 5/6/7/12/13/14 stages, all 4^th order, 2N registers explicit schemes, the coefficients of which are listed in Table 5.

Table 5. Williamson’s table of low storage Runge-Kutta schemes.

Table 5. Williamson’s table of low storage Runge-Kutta schemes.

Similarly to the TVD/SSP Runge-Kutta class, the low storage class also provides a fail-safe one-stage solver reverting back to the explicit forward Euler solver, that is useful for ROR-like frameworks.

3. Application Program Interface

Listing 1. usage example importing all public entities of FOODIE main module

3.1. The main FOODIE Abstract Data Type: the integrand type

The implemented ACP is based on one main ADT, the integrand type, the definition of which is shown in Listing 2.

Listing 2. integrand type definition

The integrand type does not implement any actual integrand field, it being and abstract type. It only specifies which deferred procedures are necessary for implementing an actual concrete integrand type that can use a FOODIE solver.

As shown in listing 2, the number of the deferred type bound procedures that clients must implement into their own concrete extension of the integrand ADT is very limited: essentially, there are 1 ODE-specific procedure plus some operators definition constituted by symmetric operators between 2 integrand objects, asymmetric operators between integrand and real numbers (and viceversa) and an assignment statement for the creation of new integrand objects. These procedures are analyzed in the following paragraphs.

3.1.1. Time derivative procedure, the residuals function

The abstract interface of the time derivative procedure t is shown in Listing 3.

Listing 3. time derivative procedure interface

This procedure-function takes two arguments, the first passed as a type bounded argument, while the latter is optional, and it returns an integrand object. The passed dummy argument, self, is a polymorphic argument that could be any extensions of the integrand ADT. The optional argument t is the time at which the residuals function must be computed: it can be omitted in the case the residuals function does not depend directly on time.

Commonly, into the concrete implementation of this deferred abstract procedure clients embed the actual ODE equations being solved. As an example, for the Burgers equation, that is a Partial Differential Equations (PDE) system involving also a boundary value problem, this procedure embeds the spatial operator that convert the PDE to a system of algebraic ODE. As a consequence, the eventual concrete implementation of this procedure can be very complex and errors-prone. Nevertheless, the FOODIE solvers are implemented only on the above abstract interface, thus emancipating the solvers implementation from any concrete complexity.

Add citations to Burgers, Adams-Bashfort, leapfrog references.

3.1.2. Symmetric operators procedures

The abstract interface of symmetric procedures is shown in Listing 4.

Listing 4. symmetric operator procedure interface

This interface defines a class of procedures operating on 2 integrand objects, namely it is used for the definition of the operators multiplication, summation and subtraction of integrand objects. These operators are used into the above described ODE solvers, for example see equations (2), (3), (6) or (9). The implementation details of such a procedures class are strictly dependent on the concrete extension of the integrand type. From the FOODIE solvers point of view, we need to known only that first argument passed as bounded one, the left-hand-side of the operator, and the second argument, the right-hand-side of the operator, are two integrand object and the returned object is still an integrand one.

3.1.3. Integrand/real and real/integrand operators procedures

The abstract interfaces of Integrand/real and real/integrand operators procedures are shown in Listing 5.

Listing 5. Integrand/real and real/integrand operators procedure interfaces

These two interfaces are necessary in order to complete the algebra operating on the integrand object class, allowing the multiplication of an integrand object for a real number, circumstance that happens in all solvers, see equations (2), (3), (6) or (9). The implementation details of these procedures are strictly dependent on the concrete extension of the integrand type. From the FOODIE solvers point of view, we need to known only that first argument passed as bounded one, the left-hand-side of the operator, and the second argument, the right-hand-side of the operator, are an integrand object and real number of viceversa and the returned object is still an integrand one.

3.1.4. Integrand assignment procedure

The abstract interface of integrand assignment procedure is shown in Listing 6.

Listing 6. integrand assignment procedure interface

The assignment statement is necessary in order to complete the algebra operating on the integrand object class, allowing the assignment of an integrand object by another one, circumstance that happens in all solvers, see equations (2), (3), (6) or (9). The implementation details of this assignment is strictly dependent on the concrete extension of the integrand type. From the FOODIE solvers point of view, we need to known only that first argument passed as bounded one, the left-hand-side of the assignment, and the second argument, the right-hand-side of the assignment, are two integrand objects.

3.2. The explicit forward Euler solver

The explicit forward Euler solver is exposed (by the FOODIE main module that must imported, see Listing 1) as a single derived type (that is a standard convention for all FOODIE solvers) named euler_explicit_integrator. It provides the type bound procedure (also referred as method) integrate for integrating in time an integrand object, or any of its polymorphic concrete extensions. Consequently, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 7.

Listing 7. definition of an explicit forward Euler integrator

Once an integrator of this type has been instantiated, it can be directly used without any initialization, for example see Listing 8.

Listing 8. example of usage of an explicit forward Euler integrator

where my_integrand is a concrete (valid) extension of integrand ADT.

The complete implementation of the integrate method of the explicit forward Euler solver is reported in Listing 9.

Listing 9. implementation of the integrate method of Euler solver

This method takes three arguments, the first argument is an integrand class, it being the integrand field that must integrated one-step-over in time, the second is the time step used and the third, that is optional, is the current time value that is passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

3.3. The explicit TVD/SSP Runge-Kutta class of solvers

The TVD/SSP Runge-Kutta class of solvers is exposed as a single derived type named tvd_runge_kutta_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• destroy: destroy the integrator previously initialized, eventually freeing the allocated dynamic memory registers;
• integrate: integrate integrand field one-step-over in time.

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 10.

Listing 10: definition of an explicit TVD/SSP Runge-Kutta integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 11.

Listing 11: example of initialization of an explicit TVD/SSP Runge-Kutta integrator

In the Listing 11 a 3-stages solver has been initialized. As a matter of facts, from the equations (3) and (4) a solver belonging to this class is completely defined once the number of stages adopted has been chosen. The complete definition of the tvd_runge_kutta_integrator type is reported into Listing 12. As shown, the Butcher’s coefficients are stored as allocatable arrays the values of which are initialized by the init method.

Listing 12: definition of tvd_runge_kutta_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 13.

Listing 13: example of usage of a TVD/SSP Runge-Kutta integrator

where my_integrand is a concrete (valid) extension of integrand ADT. Listing 13 shows that the memory registers necessary for storing the Runge-Kutta stages must be supplied by the client code.

The complete implementation of the integrate method of the explicit TVD/SSP Runge-Kutta class of solvers is reported in Listing 14.

Listing 14: implementation of the integrate method of explicit TVD/SSP Runge-Kutta class

This method takes five arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third is the stages array for storing the stages computations, the fourth is the time step used and the fifth, that is optional, is the current time value that is passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

It is worth noting that the stages memory registers, namely the array stage, must be passed as argument because it is defined as a not-passed polymorphic argument, thus we are not allowed to define it as an automatic array of the integrate method.

3.4. The explicit low storage Runge-Kutta class of solvers

The low storage variant of Runge-Kutta class of solvers is exposed as a single derived type named ls_runge_kutta_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• destroy: destroy the integrator previously initialized, eventually freeing the allocated dynamic memory registers;
• integrate: integrate integrand field one-step-over in time.

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 15.

Listing 15: definition of an explicit low storage Runge-Kutta integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 16.

Listing 6: example of initialization of an explicit low storage Runge-Kutta integrator

In the Listing 16 a 5-stages solver has been initialized. As a matter of facts, from the equation (5) a solver belonging to this class is completely defined once the number of stages adopted has been chosen. The complete definition of the ls_runge_kutta_integrator type is reported into Listing 17. As shown, the Williamson’s coefficients are stored as allocatable arrays the values of which are initialized by the init method.

Listing 17: definition of ls_runge_kutta_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 18.

Listing 18: example of usage of a low storage Runge-Kutta integrator

where my_integrand is a concrete (valid) extension of integrand ADT. Listing 18 shows that the memory registers necessary for storing the Runge-Kutta stages must be supplied by the client code, as it happens of the TVD/SSP Runge-Kutta class. However, now the registers necessary is always 2, independently on the number of stages used, that in the example considered are 5.

The complete implementation of the integrate method of the explicit low storage Runge-Kutta class of solvers is reported in Listing 19.

Listing 19: implementation of the integrate method of explicit low storage Runge-Kutta class

This method takes five arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third is the stages array for storing the stages computations, the fourth is the time step used and the fifth, that is optional, is the current time value that is passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

It is worth noting that the stages memory registers, namely the array stage, must be passed as argument because it is defined as a not-passed polymorphic argument, thus we are not allowed to define it as an automatic array of the integrate method.

3.5. The explicit Adams-Bashforth class of solvers

The explicit Adams-Bashforth class of solvers is exposed as a single derived type named adams_bashforth_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• destroy: destroy the integrator previously initialized, eventually freeing the allocated dynamic memory registers;
• integrate: integrate integrand field one-step-over in time;
• update_previous: auto update (cyclically) previous time steps solutions.

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 20.

Listing 20: definition of an explicit Adams-Bashforth integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 21.

Listing 21: example of initialization of an explicit Adams-Bashforth integrator

In the Listing 21 a 4-steps solver has been initialized. As a matter of facts, from the equation (6) a solver belonging to this class is completely defined once the number of time steps adopted has been chosen. The complete definition of the adams_bashforth_integrator type is reported into Listing 22. As shown, the linear coefficients are stored as allocatable arrays the values of which are initialized by the init method.

Listing 22 definition of adams_bashforth_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 23.

Listing 23: example of usage of an Adams-Bashforth integrator

where my_integrand is a concrete (valid) extension of integrand ADT, times are the time at each 4 steps considered for the current one-step-over integration and previous are the memory registers where previous time steps solutions are saved.

The complete implementation of the integrate method of the explicit Adams-Bashforth class of solvers is reported in Listing 24.

Listing 24: implementation of the integrate method of explicit Adams-Bashforth class

This method takes five arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third are the previous time steps solutions, the fourth is the time step used, the fifth is an array of the time values of the steps considered for the current one-step-over integration that are passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time and the sixth is a logical flag for enabling/disabling the cyclic update of previous time steps solutions. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

It is worth noting that the method also performs the cyclic update of the previous time steps solutions memory registers. This can be disable passing autoupdate=.false.: it is useful in the framework of predictor-corrector solvers.

3.6. The implicit Adams-Moulton class of solvers

The implicit Adams-Moulton class of solvers is exposed as a single derived type named adams_moulton_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• destroy: destroy the integrator previously initialized, eventually freeing the allocated dynamic memory registers;
• integrate: integrate integrand field one-step-over in time;
• update_previous: auto update (cyclically) previous time steps solutions.

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 25.

Listing 25 definition of an implicit Adams-Moulton integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 26.

Listing 26: example of initialization of an implicit Adams-Moulton integrator

In the Listing 26 a 3-steps solver has been initialized. As a matter of facts, from the equation (7) a solver belonging to this class is completely defined once the number of time steps adopted has been chosen. The complete definition of the adams_moulton_integrator type is reported into Listing 27. As shown, the linear coefficients are stored as allocatable arrays the values of which are initialized by the init method.

Listing 27: definition of adams_moulton_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 28.

Listing 28: example of usage of an Adams-Moulton integrator

where my_integrand is a concrete (valid) extension of integrand ADT, times are the time at each 4 steps considered for the current one-step-over integration and previous are the memory registers where previous time steps solutions are saved.

The complete implementation of the integrate method of the implicit Adams-Moulton class of solvers is reported in Listing 29.

Listing 29: implementation of the integrate method of explicit Adams-Moulton class

This method takes six arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third are the previous time steps solutions, the fourth is the time step used, the fifth is an array of the time values of the steps considered for the current one-step-over integration that are passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time and the sixth is a logical flag for enabling/disabling the cyclic update of previous time steps solutions. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

It is worth noting that the method also performs the cyclic update of the previous time steps solutions memory registers. This can be disable passing autoupdate=.false.: it is useful in the framework of predictor-corrector solvers.

3.7. The predictor-corrector Adams-Bashforth-Moulton class of solvers

The predictor-corrector Adams-Bashforth-Moulton class of solvers is exposed as a single derived type named adams_bashforth_moulton_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• destroy: destroy the integrator previously initialized, eventually freeing the allocated dynamic memory registers;
• integrate: integrate integrand field one-step-over in time;

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 30.

Listing 30: definition of an implicit Adams-Moulton integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 31.

Listing 31: example of initialization of an implicit Adams-Moulton integrator

In the Listing 31 a 3-steps solver has been initialized. As a matter of facts, from the equation (8) a solver belonging to this class is completely defined once the number of time steps adopted has been chosen. The complete definition of the adams_moulton_integrator type is reported into Listing 32. As shown, the linear coefficients are stored as allocatable arrays the values of which are initialized by the init method.

Listing 32: definition of adams_bashforth_moulton_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 33.

Listing 33: example of usage of an Adams-Bashforth-Moulton integrator

where my_integrand is a concrete (valid) extension of integrand ADT, times are the time at each 4 steps considered for the current one-step-over integration and previous are the memory registers where previous time steps solutions are saved.

The complete implementation of the integrate method of the predictor-corrector Adams-Bashforth-Moulton class of solvers is reported in Listing 34.

Listing 34: implementation of the integrate method of predictor-corrector Adams-Bashforth-Moulton class

This method takes four arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third is the time step used, and the fourth is the current time. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method.

3.8. The leapfrog solver

The explicit Leapfrog class of solvers is exposed as a single derived type named leapfrog_integrator. This type provides three methods:

• init: initialize the integrator accordingly the possibilities offered by the class of solvers;
• integrate: integrate integrand field one-step-over in time.

As common for FOODIE solvers, for using such a solver it must be previously defined as an instance of the exposed FOODIE integrator type, see Listing 35.

Listing 35: definition of an explicit Leapfrog integrator

Once an integrator of this type has been instantiated, it must be initialized before used, for example see Listing 36.

Listing 36: example of initialization of an explicit Leapfrog integrator

The complete definition of the leapfrog_integrator type is reported into Listing 37. As shown, the filter coefficients are initialized to zero, suitable values are initialized by the init method.

Listing 37: definition of leapfrog_integrator type

After the solver has been initialized it can be used for integrating an integrand field, as shown in Listing 38.

Listing 38: example of usage of a Leapfrog integrator

where my_integrand is a concrete (valid) extension of integrand ADT, previous are the memory registers where previous time steps solutions are saved, filter_displacement is the register necessary for computing the eventual displacement of the applied filter and times are the time at each 2 steps considered for the current one-step-over integration.

The complete implementation of the integrate method of the explicit Leapfrog class of solvers is reported in Listing 39.

Listing 39: implementation of the integrate method of explicit Leapfrog class

This method takes six arguments, the first argument is passed as bounded argument and it is the solver itself, the second is of an integrand class, it being the integrand field that must integrated one-step-over in time, the third are the previous time steps solutions, the fourth is the optional filter-displacement-register, the fifth is the time step used and the sixth is an array of the time values of the steps considered for the current one-step-over integration that are passed to the residuals function for taking into account the cases where the time derivative explicitly depends on time. The time step is not automatically computed (for example inspecting the passed integrand field), thus its value must be externally computed and passed to the integrate method. It is worth noting that if the filter displacement argument is not passed, the solver reverts back to the standard unfiltered Leapfrog method.

It is worth noting that the method also performs the cyclic update of the previous time steps solutions memory registers. In particular, if the filter displacement argument is passed the method performs the RAW filtering.

4. Tests and Examples

4.1. Oscillation equations test

Let us consider the oscillator problem, it being a simple, yet interesting IVP. Briefly, the oscillator problem is a prototype problem of non dissipative, oscillatory phenomena. For example, let us consider a pendulum subjected to the Coriolis accelerations without dissipation, the motion equations of which can be described by the ODE system (11).

$$\begin{array}{c}{U}_t=R\left(U\right)\\ U=\left[\begin{array}{c}x\\ y\end{array}\right]R\left(U\right)=\left[\begin{array}{c}-fy\\ fx\end{array}\right]\end{array}$$

(11)

where the frequency is chosen as $f={10}^4$. The ODE system (11) is completed by the following initial conditions:

$$\begin{array}{c}x\left(t_0\right)=0\\ y\left(t_0\right)=1\end{array}$$

(12)

where $t_0=0$ is the initial time considered.

The IVP constituted by equations (11) and (12) is (apparently) simple and its exact solution is known:

$$\begin{array}{c}x(t_0+\Delta t)=X_{0}cos(f\Delta t)-y_{0}sin(f\Delta t)\\ y(t_0+\Delta t)=X_{0}sin(f\Delta t)+y_{0}cos(f\Delta t)\end{array}$$

(13)

where $\Delta t$ is an arbitrary time step. This problem is non-stiff meaning that the solution is constituted by only one time-scale, namely the single frequency f.

This problem is only apparently simple. As a matter of facts, in a non dissipative oscillatory problem the eventual errors in the amplitude approximation can rapidly drive the subsequent series of approximations to an unphysical solution. This is of particular relevance if the solution (that is numerically approximated) constitutes a prediction far from the initial conditions, that is the common case in weather forecasting.

Because the Oscillation system (11) posses a closed exact solution, the discussion on this test has twofolds aims: to assess the accuracy of the FOODIE’s built-in solvers comparing the numerical solutions with the exact one and to demonstrate how it is simple to solve this prototypical problem by means of FOODIE.

4.1.1. Errors Analysis

For the analysis of the accuracy of each solver, we have integrated the Oscillation equations (11) with different, decreasing time steps in the range $[5000,2500,1250,625,320,100]$. The error is estimated by the L2 norm of the difference between the exact ($U_e$) and the numerical ($U_{\Delta t}$) solutions for each time step:

$$\epsilon (\Delta t)=\left|\right|U_e-U_{\Delta t}{\left|\right|}_2=\sqrt{\sum _{s=1}^{N_s}{\left(U_e(t_0+s*\Delta t)-U_{\Delta t}(t_0+s*\Delta t)\right)}^2}$$

(14)

where $N_s$ is the total number of time steps performed to reach the final integration time.

Using two pairs of subsequent-decreasing time steps solution is possible to estimate the order of accuracy of the solver employed computing the observed order of accuracy:

$$p=\frac{log10\left(\frac{\epsilon (\Delta t_1)}{\epsilon (\Delta t_2)}\right)}{log10\left(\frac{\Delta t_1}{\Delta t_2}\right)}$$

(15)

where $\frac{\Delta t_1}{\Delta t_2}>1$.

4.1.2. FOODIE aware implementation of an oscillation numerical solver

The IVP (11) can be easily solved by means of FOODIE library. The first block of a FOODIE aware solution consists to define an oscillation integrand field defining a concrete extension of the FOODIE integrand type. Listing 40 reports the implementation of such an integrand field that is contained into the tests suite shipped within the FOODIE library.

Listing 40: implementation of the oscillation integrand type

The oscillation field extends the integrand ADT making it a concrete type. This derived type is very simple: it has 5 data members for storing the state vector and some auxiliary variables, and it implements all the deferred methods necessary for defining a valid concrete extension of the integrand ADT (plus 2 auxiliary methods that are not relevant for our discussion). The key point is here constituted by the implementation of the deferred methods: the integrand ADT does not impose any structure for the data members, that are consequently free to be customized by the client code. In this example the data members have a very simple, clean and concise structure:

• $dims$ is the number of space dimensions; in the case of equation (11) we have $dims=2$, however the test implementation has been kept more general parametrizing this dimension in order to easily allow future modification of the test-program itself;
• f stores the frequency of the oscillatory problem solved, that is here set to ${10}^4$, but it can be changed at runtime in the test-program;
• U is the state vector corresponding directly to the state vector of equation (11);

As the Listing 40 shows, the FOODIE implementation strictly corresponds to the mathematical formulation embracing all the relevant mathematical aspects into one derived type, a single object. Here we not review the implementation of all deferred methods, this being out of the scope of the present work: the interested reader can see the tests suite sources shipped within the FOODIE library. However, some of these methods are relevant for our discussion, thus they are reviewed.

dOscillation_dt, the oscillation residuals function

Probably, the most important methods for an IVP solver is the residuals function. As a matter of facts, the ODE equations are implemented into the residuals function. However, the FOODIE ADT strongly alleviates the subtle problems that could arise when the ODE solver is hard-implemented within the specific ODE equations. As a matter of facts, the integrand ADT specifies the precise interface the residuals function must have: if the client code implements a compliant interface, the FOODIE solvers will work as expected, reducing the possible errors location into the ODE equations, having designed the solvers on the ADT and not on the concrete type.

Listing 41 reports the implementation of the oscillation residuals function: it is very clear and concise. Moreover, comparing this listing with the equation (11) the close correspondence between the mathematical formulation and Fortran implementation is evident.

Listing 41: implementation of the oscillation integrand residuals function

Add method, an example of oscillation symmetric operator

As a prototype of the operators overloading let us consider the add operator, it being a prototype of symmetric operators, the implementation of which is presented in Listing 42.

Listing 42: implementation of the oscillation integrand add operator

It is very simple and clear: firstly all the auxiliary data are copied into the operator result, then the state vector of the result is populated with the addiction between the state vectors of the left-hand-side and right-hand-side. This is very intuitive from the mathematical point of view and it helps to reduce implementation errors. Similar implementations are possible for all the other operators necessary to define a valid intregrand ADT concrete extension.

assignment of an oscillation object

The assignment overloading of the oscillation type is the last key-method that enforces the conciseness of the FOODIE aware implementation. Listing 43 reports the implementation of the assignment overloading. Essentially, to all the data members of the left-hand-side are assigned the values of the corresponding right-hand-side. Notably, for the assignment of the state vector and of the previous time steps solution array we take advantage of the automatic re-allocation of the left-hand-side variables when they are not allocated or allocated differently from the right-hand-side, that is a Fortran 2003 feature. In spite its simplicity, the assignment overloading is a key-method enabling the usage of FOODIE solver: effectively, the assignment between two integrand ADT variables is ubiquitous into the solvers implementations, see equation (3) for example.

Listing 43: implementation of the oscillation integrand assignment

FOODIE numerical integration

Using the above discussed oscillation type it is very easy to solve IVP (11) by means of FOODIE library. Listing 44 presents the numerical integration of system (11) by means of the Leapfrog RAW-filtered method. In the example, the integration is performed with ${10}^4$ steps with a fixed $\Delta t={10}^2$ until the time $t={10}^6$ is reached. The example shows also that for starting a multi-step scheme such as the Leapfrog one a lower-oder or equivalent order one-scheme is necessary: in the example the first 2 steps are computed by means of one-step TVD/SSP Runge-Kutta 2-stages schemes. Note that the memory registers for storing the Runge-Kutta stages and the RAW filter displacement must be handled by the client code. Listing 44 demonstrates how it is simple, clear and concise to solve a IVP by FOODIE solvers. Moreover, it proves how it is simple and effective to apply different solvers in a coupled algorithm, that greatly simplify the development of new hybrid solvers for self-adaptive time step size.

Listing 44: numerical integration of the oscillation system by means of Leapfrog RAW-filtered method

4.1.3. Adams-Bashforth

Table 9 summarizes the Adams-Bashforth error analysis. As expected, the Adams-Bashforth 1 step solution, that reverts back to the explicit forward Euler one, is unstable for all the $\Delta t$ exercised.

The expected observed orders of accuracy for the Adams-Bashforth solvers using 2, 3 and 4 time steps tend to 1.5, 2.5 and 3.5 that are in agreement with the expected formal order of 2, 3 and 4, respectively. Comparing the errors of the finest time resolution, i.e. $\Delta t=100$, we find that the L2 norm decreases of the 2 orders of magnitude as the solver’s accuracy increases by 1 order. This also means that fixing a tolerance on the errors, the higher is the solver’s accuracy the larger is the time resolution available. As an example, assuming that admissible errors are of $O\left({10}^{-2}\right)$ with the 4-steps solver we can use $\Delta t=625$ performing $N_s=t_{f}inal/625$ numerical integration steps, whereas using a 3-steps solvers we must adopt $\Delta t=100$ performing $6.25\times N_s$ numerical integration steps instead of $N_s$. Considering that the computational costs is only slightly affected by the number of previous time steps considered3, the accuracy order has strong impact on the overall numerical efficiency: to improve the numerical efficiency reducing the computational costs, the usage of high order Adams-Bashforth solvers with larger time steps should be preferred instead of low order solvers with smaller time steps.

Table 9. Oscillation test: errors analysis of explicit Adams-Bashforth solvers

Table 9. Oscillation test: errors analysis of explicit Adams-Bashforth solvers

Figure 1 shows, for each solver exercised, the $X\left(t\right)$ and $Y\left(t\right)$ solution for $t\in [0,{10}^6]$: the plots into the figure report a global overview of the solution for all the instants considered (left subplots) and a detailed zoom over the last instants of the integration (right subplots) for evaluating the numerical errors accumulation. For the sake of clarity, the strongly unstable solutions are omitted into the subplots concerning the final integration instants, namely the solutions for large $\Delta t$. Figure 1 emphasizes the instability generation for some pairs steps number/$\Delta t$. The 2 and 4 steps solutions are instable for $\Delta t=5000\to f*\Delta t=0.5$. On the contrary, the 3 steps solution is stable, but the amplitude is dumped and the solution vanishes as the integration proceeds. The 2 and 4 steps solutions show a phase error that decreases as the time resolution increases, whereas 3 steps solution has null phase error.

4.1.4. Adams-Bashforth-Moulton

Table 10 summarizes the Adams-Bashforth-Moulton error analysis. The same considerations done for the Adams-Bashforth solutions can repeated for the Adams-Bashforth-Moulton ones, thus they are omitted for the sake of conciseness. An interesting result concerns the observed errors: the $O\left({10}^{-2}\right)$ error is now obtained with $\Delta t=1250$ for the 4-steps solver, thus it is 2 times faster than the corresponding Adams-Bashforth 4-step solver. Considering that the computational costs of a single Adams-Bashforth-Moulton step is only slightly greater than the corresponding Adamas-Bashforth step, the efficiency increasing is not negligible.

Figure 2 shows similar plots of Figure 1 above discussed. Differently from the Adams-bashforth class, the amplitude damping feature is now possessed by the 2-steps solver, see plot Figure 2b, while all solutions show phase errors that decrease as the time resolution increases.

Table 10. Oscillation test: errors analysis of predictor-corrector Adams-Bashforth-Moulton solvers

Table 10. Oscillation test: errors analysis of predictor-corrector Adams-Bashforth-Moulton solvers

4.1.5. Adams-Moulton

Table 11 summarizes the Adams-Moulton error analysis. The implicit Adams-Moulton solvers behave much like the Adams-Bashforth-Moulton ones: they have similar errors and observed orders for the same formal order considered. However, the implicit Adams-Moulton class uses one less step with respect the corresponding Adams-Bashforth-Moulton class: this could lead to the promise of higher computational efficiency. Notwithstanding, for solving the implicit non-linearity embedded into the Adams-Moulton solvers an iterative algorithm must be employed: for the results presented, a 5 iterations of fixed point algorithm have been computed. This strongly reduces the eventual gain of computational efficiency with respect the Adams-Bashforth-Moulton class.

Table 11. Oscillation test: errors analysis of explicit Adams-Moulton solvers; the implicit non-linearity is solved by 5 iterations of fixed point algorithm

Table 11. Oscillation test: errors analysis of explicit Adams-Moulton solvers; the implicit non-linearity is solved by 5 iterations of fixed point algorithm

Figure 3 shows similar plots of Figure 2 above discussed: there are not relevant differences between the 2 classes of solvers.

4.1.6. Leapfrog

The Leapfrog solutions are in agreement with the expected results: both unfiltered and RAW-filtered solutions show an observed order of accuracy that tends to the formal 2^nd order. The two solutions are almost the same.

Table 12. Oscillation test: errors analysis of explicit Leapfrog solvers.

Table 12. Oscillation test: errors analysis of explicit Leapfrog solvers.

Figure 4. Oscillation equations solutions computed by means of Leapfrog solvers.

Figure 4. Oscillation equations solutions computed by means of Leapfrog solvers.

4.1.7. Low Storage Runge-Kutta

Table 13. Oscillationtest: errors analysis of explicit Low Storage Runge-Kutta solvers.

Table 13. Oscillationtest: errors analysis of explicit Low Storage Runge-Kutta solvers.

Figure 5. Oscillation equations solutions computed by means of low storage Runge-Kutta solvers.

Figure 5. Oscillation equations solutions computed by means of low storage Runge-Kutta solvers.

Figure 6. Oscillation equations solutions computed by means of low storage Runge-Kutta solvers.

Figure 6. Oscillation equations solutions computed by means of low storage Runge-Kutta solvers.

4.1.8. TVD/SSP Runge-Kutta

Table 14. Oscillation test: errors analysis of explicit TVD/SSP Runge-Kutta.

Table 14. Oscillation test: errors analysis of explicit TVD/SSP Runge-Kutta.

Figure 7. Oscillation equations solutions computed by means of TVD/SSP Runge-Kutta solvers.

Figure 7. Oscillation equations solutions computed by means of TVD/SSP Runge-Kutta solvers.

5. Benchmarks on parallel frameworks

• CoArray Fortran (CAF) model, for shared and distributed memory architectures;
• OpenMP directive-based model, for only shared memory architectures;

5.1. CAF benchmark

This subsection reports the parallel scaling analysis of Euler 1D test programs (with and without FOODIE) being parallelized by means of CoArrays Fortran (CAF) model. This parallel model is based on the concept of coarray introduced into the Fortran 2008 standard: the array syntax is extended introducing the so called codimension that is a new arrays indexing. Essentially, a CAF enabled code is designed to be replicated a certain number of times and all copies, conventionally named images, are executed asynchronously. Each image has its own set of data (memory) and the codimension indexes are used to access to the (remote) memory of the other images. The CAF approach allows the implementation of Partitioned Global Address Space (PGAS) model following the SPMD (single program, multiple data) parallelization paradigm. The programmer must take care of defining the coarray variables and of synchronizing the images when necessary. This approach requires the refactoring of legacy serial codes, but it allows the exploitation of both shared and distributed memory architectures. Moreover, it is a standard feature of Fortran (2008), thus it is not chained to any particular compiler vendors extension.

The benchmarks shows in this section have been done on a dual Intel(R) Xeon(R) CPU X5650 exacores workstation for a total of 12 physical cores, coupled with 24GB of RAM. In order to perform an accurate analysis 4 different codes have considered:

• FOODIE-aware codes:

  -

  serial code;

  -

  CAF-enabled code;

• procedural codes without using FOODIE library:

  -

  serial code;

  -

  CAF-enabled code;

These codes (see A.1 for the implementation details) have been compiled by means of the GNU gfortran compiler v5.2.0 coupled with OpenCoarrays v1.1.04.

The Euler conservation laws are integrated for 30 time steps by means of the TVD RK(5,4) solver: the measured CPU time used for computing the scaling efficiencies is the average of the 30 integrations, thus representing the mean CPU time for computing one time step integration.

For the strong scaling, the benchmark has been conducted with 240000 finite volumes. Figure 8a summarizes the strong scaling analysis: it shows that FOODIE-based code scales similarly to the baseline code without FOODIE.

For the weak scaling the minimum size is 24000 finite volumes and the size is scaled linearly with the CAF images, thus $N_{12}=288000$ cells. Figure 8b summarizes the weak scaling analysis and it essentially confirms that FOODIE-based code scales similarly to the baseline code without FOODIE.

Both strong and weak scaling analysis point out that for the computing architecture considered the parallel scaling is reasonable up to 12 cores, the efficiency being always satisfactory.

To complete the comparison, the absolute CPU-time consumed by the two families of codes (with and without FOODIE) must be considered. Table 15 summarizes the benchmarks results. As shown, procedural and FOODIE-aware codes consume a very similar CPU-time for both the strong and the weak benchmarks. The same results are shown in Figure 9. These results prove that the abstraction of FOODIE environment does not degrade the computational efficiency.

Table 15. caf benchmarks results

Table 15. caf benchmarks results

5.2. OpenMP benchmark

This subsection reports the parallel scaling analysis of Euler 1D test programs (with and without FOODIE) being parallelized by means of OpenMP directives-based paradigm. This parallel model is based on the concept of threads: an OpenMP enabled code start a single (master) threaded program and, at run-time, it is able to generate a team of (many) threads that work concurrently on the parallelized parts of the code, thus reducing the CPU time necessary for completing such parts. The parallelization is made by means of directives explicitly inserted by the programmer: the communications between threads are automatically handled by the compiler (through the provided OpenMP library used as back-end). OpenMP parallel paradigm is not a standard feature of Fortran, rather it is an extension provided by the compiler vendors. This parallel paradigm constitutes an effective and easy approach for parallelizing legacy serial codes, however its usage is limited to shared memory architectures because all threads must have access to the same memory.

The benchmarks shown in this section have been done on a dual Intel(R) Xeon(R) CPU X5650 exacores workstation for a total of 12 physical cores, coupled with 24GB of RAM. In order to perform an accurate analysis 4 different codes have considered:

• FOODIE-aware codes:

  -

  serial code;

  -

  OpenMP-enabled code;

• procedural codes without using FOODIE library:

  -

  serial code;

  -

  OpenMP-enabled code;

These codes (see A.1 for the implementation details) have been compiled by means of the GNU gfortran compiler v5.2.0 with -O2 -fopenmp compilation flags.

The Euler conservation laws are integrated for 30 time steps by means of the TVD RK(5,4) solver: the measured CPU time used for computing the scaling efficiencies is the average of the 30 integrations, thus representing the mean CPU time for computing one time step integration.

For the strong scaling, the benchmark has been conducted with 240000 finite volumes. Figure 10a summarizes the strong scaling analysis: it shows that FOODIE-based code scales similarly to the baseline code without FOODIE.

For the weak scaling the minimum size is 24000 finite volumes and the size is scaled linearly with the OpenMP threads, thus $N_{12}=288000$ cells. Figure 10b summarizes the weak scaling analysis and it essentially confirms that FOODIE-based code scales similarly to the baseline code without FOODIE.

Both strong and weak scaling analysis point out that for the computing architecture considered the parallel scaling is reasonable up to 8 cores: using 12 cores the measured efficiencies become unsatisfactory, reducing below the 60%.

To complete the comparison, the absolute CPU-time consumed by the two families of codes (with and without FOODIE) must be considered. Table 16 summarizes the benchmarks results. As shown, procedural and FOODIE-aware codes consume a very similar CPU-time for both the strong and the weak benchmarks. The same results are shown in Figure 11. These results prove that the abstraction of FOODIE environment does not degrade the computational efficiency.

Table 16. OpenMP benchmarks results

Table 16. OpenMP benchmarks results

6. Concluding Remarks and Perspectives

References

1. Ince, E.L. Ordinary differential equations; Courier Corporation, 1956.
2. Jose Luis Galan–García and Pedro Rodriguez–Cielos and Maria Ángeles Galan–Garcia and Morgan le Goff and Yolanda Padilla–Dominguez and Pablo Rodríguez–Padilla and Ivan Atencia and Gabriel Aguilera–Venegas. SODES: Solving ordinary differential equations step by step. Journal of Computational and Applied Mathematics 2023, 428, 115127. [Google Scholar] [CrossRef]
3. Owoyele, O.; Pal, P. ChemNODE: A neural ordinary differential equations framework for efficient chemical kinetic solvers. Energy and AI 2022, 7, 100118. [Google Scholar] [CrossRef]
4. Nagy, D.; Plavecz, L.; Hegedus, F. The art of solving a large number of non-stiff, low-dimensional ordinary differential equation systems on GPUs and CPUs. Communications in Nonlinear Science and Numerical Simulation 2022, 112, 106521. [Google Scholar] [CrossRef]
5. Nascimento, R.G.; Fricke, K.; Viana, F.A. A tutorial on solving ordinary differential equations using Python and hybrid physics-informed neural network. Engineering Applications of Artificial Intelligence 2020, 96, 103996. [Google Scholar] [CrossRef]
6. Backus, J.W.; Beeber, R.J.; Best, S.; Goldberg, R.; Haibt, L.M.; Herrick, H.L.; Nelson, R.A.; Sayre, D.; Sheridan, P.B.; Stern, H.; Ziller, I.; Hughes, R.A.; Nutt, R. The FORTRAN automatic coding system. Papers presented at the February 26-28, 1957, western joint computer conference: Techniques for reliability. ACM: New York, NY, USA, 1957; IRE-AIEE-ACM ’57 (Western), pp. 188–198. [CrossRef]
7. Draft International Standard for Fortran 2003. Published by the International Fortran standards committee, ISO/IEC JTC1/SC22, 2004.
8. Rouson, D.; Xia, J.; Xu, X. Scientific software design: the object-oriented way; Cambridge University Press, 2011.
9. Butcher, J. Numerical Methods for Ordinary Differential Equations, 2 ed.; Wiley, 2008; p. 482.
10. Gottlieb, S.; Ketcheson, D.I.; Shu, C.W. High order strong stability preserving time discretizations. Journal of Scientific Computing 2009, 38, 251–289. [Google Scholar] [CrossRef]
11. Maury Jr, J.; Segal, G. Cowell type numerical integration as applied to satellite orbit computation. Technical report, 1969.
12. Shu, C.W. Total-variation-diminishing time discretizations. SIAM Journal on Scientific and Statistical Computing 1988, 9, 1073–1084. [Google Scholar] [CrossRef]
13. Williamson, J.H. Low-storage runge-kutta schemes. Journal of computational physics 1980, 35, 48–56. [Google Scholar] [CrossRef]
14. Carpenter, M.H.; Kennedy, C.A. Fourth-order 2N-storage Runge-Kutta schemes. Technical report, 1994.
15. Allampalli, V.; Hixon, R.; Nallasamy, M.; Sawyer, S.D. High-accuracy large-step explicit Runge–Kutta (HALE-RK) schemes for computational aeroacoustics. Journal of Computational Physics 2009, 228, 3837–3850. [Google Scholar] [CrossRef]
16. Niegemann, J.; Diehl, R.; Busch, K. Efficient low-storage Runge–Kutta schemes with optimized stability regions. Journal of Computational Physics 2012, 231, 364–372. [Google Scholar] [CrossRef]
17. Robert, A.J. The integration of a low order spectral form of the primitive meteorological equations. Journal of the Meteorological Society of Japan. Ser. II 1966, 44, 237–245. [Google Scholar] [CrossRef]
18. Asselin, R. Frequency filter for time integrations. Monthly Weather Review 1972, 100, 487–490. [Google Scholar] [CrossRef]
19. Williams, P.D. The RAW filter: An improvement to the Robert–Asselin filter in semi-implicit integrations. Monthly Weather Review 2011, 139, 1996–2007. [Google Scholar] [CrossRef]
20. Zaghi, S.; Di Mascio, A.; Favini, B. Application of WENO-Positivity-Preserving Schemes to Highly Under-Expanded Jets. Journal of Scientific Computing 2016, 69, 1033–1057. [Google Scholar] [CrossRef]
21. Zaghi, S. OFF, Open source Finite volume Fluid dynamics code: A free, high-order solver based on parallel, modular, object-oriented Fortran API. Computer Physics Communications 2014, pp. [Google Scholar] [CrossRef]