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 stepbystep 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
Contents
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^{n1}) \end{equation*}
where \(n\) is the order. Then it can be redefined as (transformation),
\begin{align*} y_1 &= y \\ y_j &= y^{j1} \quad j=1,2,\cdots n1 \end{align*}
Finally, the system could be represented as (differentiation),
\begin{align*} f_j = y'_j &= y_{j+1} \quad j=1,2,\cdots n1\\ 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 massspring 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*}then,
\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 $ TestODE 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 OpenFOAM2.3 version, the ordinary differential equation solvers have been updated [4]. Therefore Zongyuan [3] report has some outdated parts.
Implementation
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 TestODE.C, open it in any text editor (I prefer atom or qtcreator 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.
myODE
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_; public: myODE(const scalar& mass, const scalar& stiffness) : ODESystem(), m_(mass), k_(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.
Finally, compile the code and run it. Direct the output to a file so it could be used for postprocessing.
$ wmake $ TestODE RKCK45 > log
Validation
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
Comments powered by Disqus