How to solve ODE in OpenFOAM

OpenFOAM [1] has a library designed to solve ordinary differential equations (ODEs). It is widely used for solving chemistry in OpenFOAM. This post focuses on how to solve any ODE in OpenFOAM. it is not intended to be a complete review or explaining every ODE solver, however, this post shows step-by-step how to define an ODE and how to write a small program to solve it. Finally, the results will be compared and plotted aginst the extract solution using Gnuplot.

Mathematical background

In order to solve any high order ODE numerically, it should be reduced first to a system of first order ODEs. This process of order reduction produces a number of ODEs equal to the order of the original ODE. Here, the reduction process as described by Ferziger [2] for general high order ODE,

\begin{equation*} y^n = f(x,y,y',\cdots y^{n-1}) \end{equation*}

where \(n\) is the order. Then it can be redefined as (transformation),

\begin{align*} y_1 &= y \\ y_j &= y^{j-1} \quad j=1,2,\cdots n-1 \end{align*}

Finally, the system could be represented as (differentiation),

\begin{align*} f_j = y'_j &= y_{j+1} \quad j=1,2,\cdots n-1\\ f_n = y'_n &= f(x,y_1,y_2,\cdots y_n) \end{align*}

As described above, the process could be summarised into two steps; a transformation of variables and differentiation. In OpenFOAM, there is an extra step required for stiff system solvers only. These solvers require the jacobian matrix and the time partial derivatives.

\begin{equation*} J = \begin{bmatrix} \frac{\partial f_1}{\partial y_1} & \frac{\partial f_1}{\partial y_2} & \cdots &\frac{\partial f_1}{\partial y_n}\\ \frac{\partial f_2}{\partial y_1} & \frac{\partial f_2}{\partial y_2} & \cdots &\frac{\partial f_2}{\partial y_n}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial y_1} & \cdots & \cdots &\frac{\partial f_n}{\partial y_n} \end{bmatrix} \end{equation*}

Where \(f_n\) is the equation \(n\) in the reduced system. Also the partial derivative with respect to \(x\) (independent variable)

\begin{equation*} \frac{\partial f_1}{\partial x} \cdots \frac{\partial f_n}{\partial x} \end{equation*}

Equation of motion

The equation of motion of single degree of freedom system (simple mass-spring system) , which is a second order ODE, will be used here as an example.

\begin{align*} m\ddot{q} + kq &= 0 \\ \ddot{q} &= -\frac{k}{m}q \end{align*}

where \(m\) is the mass, \(k\) is the stiffness and \(q\) is the displacement.

This equation has an exact solution represents the harmonic oscillation of the system as (assuming zero initial velocity);

\begin{equation*} q = q_0 \cos(\omega t) \end{equation*}

where \(q_0\) is the initial displacement and \(\omega\) is the natural frequency of the system \(\omega=\sqrt{k/m}\).

It is the time to apply the methodology described in the previous section (Mathematical background) on the equation of motion.

  • Step A: Transformation from \(q\) to \(y\)

    \begin{align*} y_1 &= q \\ y_2 &= \dot{q} \end{align*}
  • Step B: Differentiation (substitute first)

    \begin{equation*} \ddot{q} = -\frac{k}{m}y_1 \end{equation*}


    \begin{align*} f_1 = \dot{y_1} &= \dot{q} \\ &= y_2\\ f_2 = \dot{y_2} &= \ddot{q} \\ &= -\frac{k}{m}y_1 \end{align*}
  • Step C: Jacobian (optional)

    The Jacobian matrix is required for stiff solvers only. For this case, it is just \((2\times2)\) matrix,

    \begin{equation*} J = \begin{bmatrix} 0.0 & 1.0 \\ -k/m & 0 \end{bmatrix} \end{equation*}

    and time derivatives (independent variable in this case),

    \begin{align*} \frac{\partial f_1}{\partial t} &= 0 \\ \frac{\partial f_2}{\partial t} &= 0 \end{align*}

Where to start?

It is typical in OpenFOAM to search for a similar implementation as a starting point. Fortunately, there is a small test program included in OpenFOAM source files. It is described in details in [3]. So copy it to your run directory then compile and run it. (RKCK45 is one of OpenFOAM ODE solvers).

$ cp -r $FOAM_APP/test/ODE .
$ wmake
$ Test-ODE RKCK45

You should get a lot of numbers on the screen and finally you should see;

Analytical: y(2.0) = 4(0.223891 0.576725 0.352834 0.128943)
Numerical:  y(2.0) = 4(0.223891 0.576725 0.352835 0.128943), dxEst = 0.402302

Please note that since OpenFOAM-2.3 version, the ordinary differential equation solvers have been updated [4]. Therefore Zongyuan [3] report has some outdated parts.


The reduction step is essential to solve the equation and as it is described above. This step is important in OpenFOAM and in any other code like python or Octave/Matlab. The implementation in OpenFOAM is divided into two main classes; ODEsystem and ODEsolver. From its name you can guess the function of each class. ODEsystem is an abstract class defines the system of first order ODEs as explained above. ODEsolver is the base class for all the ODE solvers in OpenFOAM.

Back to Test-ODE.C, open it in any text editor (I prefer atom or qt-creator for big projects) and start examining the code. You will notice that the code is divided into two parts, the first part is the definitionn of class called testODE and the second part is the main function.


testODE is basically the ODE system definition which is inherited from the base abstract class ODEsystem. Let's modify it to represent our equation. Before we start just rename it myODE. Also, modify the class constructor to allow passing \(m\) and \(k\) values from the main application.

class myODE
public ODESystem
    //- Mass of the system
    const scalar m_;
    //- Stiffness of the system
    const scalar k_;


    myODE(const scalar& mass, const scalar& stiffness)

This class has only three functions, namely, nEqns, derivatives and jacobian. nEqns represents the number of equation of this system which is essentially the order of the original equation (for this case is \(2\)).

label nEqns() const
    return 2;

The second function is derivatives which is the definition of the system of the first order ODEs (Step B);

void derivatives
    const scalar x,
    const scalarField& y,
    scalarField& dydx
) const
    dydx[0] = y[1];
    dydx[1] = -(k_/m_+VSMALL)*y[0]

The third and final function in this class is jacobian and it represents the jacobian matrix (Step C) and the time derivatives of the system.

void jacobian
    const scalar x,
    const scalarField& y,
    scalarField& dfdx,
    scalarSquareMatrix& dfdy
) const
    dfdx[0] = 0.0;
    dfdx[1] = 0.0;

    dfdy[0][0] = 0.0;
    dfdy[0][1] = 1.0;

    dfdy[1][0] = -(k_/m_+VSMALL);
    dfdy[1][1] = 0;

main program

Now we have myODE and it is time to solve the system. The first line inside the main function is basically to allow the program to read the ODEsolver name from the command line. We will keep it and add extra two parameters for \(m\) and \(k\) (lines 7 & 8).

Please read the comments inside the code carefully. One point could cause confusion here is yStart (line 31). It is scalarField has the same size as the system (order). yStart has two functions here, the fist one is applying the initial conditions. In our case and according to (Step A), these initial conditions are the system initial displacement and velocity. The second function of yStart is storing the output results of the ODEsolver.

int main(int argc, char *argv[])
    argList args(argc, argv);

    // system properties (mass and stiffness)
    const scalar mass = 2.0; // Kg
    const scalar stiff = 2048.0; // N/m
    // Initial displacement
    const scalar y0 = 0.2; // m
    const label n = 200;
    const scalar endTime = 1.0; // s

    // Create the ODE system as object of type myODE
    myODE ode(mass,stiff);

    // Create dictionary and add the odeSolver name
    dictionary dict;
    dict.add("solver", args[1]);

    // Create the selected ODE system solver
    autoPtr<ODESolver> odeSolver = ODESolver::New(ode, dict);

    // Initialise the ODE system fields

    // Start time
    scalar xStart = 0.0;
    // time step
    const scalar dx = endTime/n;
    // Initial displacement and velocity
    scalarField yStart(ode.nEqns());
    yStart[0] = y0;
    yStart[1] = 0.0;
    // integration initial step
    scalar dxEst = 0.1;
    scalar xEnd = 0;

    // Required to store dydx
    scalarField dyStart(ode.nEqns());

    // Integration loop
    for (label i=0; i<n; i++)
        xEnd = xStart + dx;
        ode.derivatives(xStart, yStart, dyStart);
        odeSolver->solve(xStart, xEnd, yStart, dxEst);
        xStart = xEnd;
        Info << xStart << "   " << yStart[0] << endl;

    Info<< "\nEnd\n" << endl;

    return 0;

Finally, compile the code and run it. Direct the output to a file so it could be used for postprocessing.

$ wmake
$ Test-ODE RKCK45 > log


As mentioned before (Equation of motion), this ODE has an exact solution. It is basically a harmonic oscillation which could be represented as (after applying the initial conditions and system properties),

\begin{equation*} q(t) = 0.2 \cos(32t) \end{equation*}

A quick and easy way to compare the results is by using gnuplot to visualise the results.

Gnuplot Script:
# system properties
mass = 2.0
stiff = 2048.0
x0 = 0.2

set grid
set samples 200
set xlabel 'Time (s)'
set ylabel 'Displacement (m)'
set yrange [-0.3:0.3]

f(x) = x0*cos(sqrt(stiff/mass)*x)
plot 'log' u 1:2 w p lw 2 t "OpenFOAM", f(x) w l lw 2 t "Exact"

Final remarks

This post is meant to be minimal as possible and to be focused on solving ODE. Therefore, there are many possibilities to improve the main program. The system properties can be read from the command line directly or , even better, all the inputs could be included in one dictionary (OpenFoam style). Also, the output result could be formatted in a better elegant way. Moreover, the results can be written to a file. Finally, the solve function is an overloaded function, so it can be used with different input argument.

Please feel free to comment below. Your feedback will be highly appreciated.

  • Update: 23 June 2016, scalar type added to myODE constructor as pointed out by Mohamed Ouda
  • Update: 28 June 2016, variables names are revised to be consistent and more clear according to Francisco Angel comments
[1] OpenFOAM® and OpenCFD® are registered trademarks of OpenCFD Limited, the producer OpenFOAM software. All registered trademarks are property of their respective owners. This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks. Hassan Kassem is not associated to OpenCFD.
[2] Ferziger, Joel H. 1998. Numerical Methods for Engineering Applications. Wiley.
[3] Zongyuan, Gu 2009. Introduction to ODE solvers and their application in OpenFOAM. Chalmers.
[4] Numerical methods. OpenFOAM 2.3.0 relase notes 2014.


Comments powered by Disqus