What is the time? (AKA Hello world!) is a very simple model. The model consists of a single variable (called ‘time’) and a single differential equation:
d(time)/dt = 1
This way, the value of the variable ‘time’ is equal to the elapsed time in the simulation at any moment.
This tutorial presents the basic structure of daeModel and daeSimulation classes. A typical DAETools simulation requires the following 8 tasks:
The ‘time’ variable plot:
Files
Model report  whats_the_time.xml 
Runtime model report  whats_the_timert.xml 
Source code  whats_the_time.py 
This tutorial introduces several new concepts:
In this example we model a simple heat conduction problem: a conduction through a very thin, rectangular copper plate.
For this problem, we need a twodimensional Cartesian grid (x,y) (here, for simplicity, divided into 10 x 10 segments):
y axis
^

Ly  L T T T T T T T T T R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
 L + + + + + + + + + R
0  L B B B B B B B B B R
> x axis
0 Lx
Points ‘B’ at the bottom edge of the plate (for y = 0), and the points ‘T’ at the top edge of the plate (for y = Ly) represent the points where the heat is applied.
The plate is considered insulated at the left (x = 0) and the right edges (x = Lx) of the plate (points ‘L’ and ‘R’). To model this type of problem, we have to write a heat balance equation for all interior points except the left, right, top and bottom edges, where we need to specify boundary conditions.
In this problem we have to define two distribution domains:
the following parameters:
and a single variable:
The model consists of 5 equations (1 distributed equation + 4 boundary conditions):
Heat balance:
rho * cp * dT(x,y)/dt = k * [d2T(x,y)/dx2 + d2T(x,y)/dy2]; x in (0, Lx), y in (0, Ly)
Neumann boundary conditions at the bottom edge:
k * dT(x,y)/dy = Qb; x in (0, Lx), y = 0
Dirichlet boundary conditions at the top edge:
T(x,y) = Tt; x in (0, Lx), y = Ly
Neumann boundary conditions at the left edge (insulated):
dT(x,y)/dx = 0; y in [0, Ly], x = 0
Neumann boundary conditions at the right edge (insulated):
dT(x,y)/dx = 0; y in [0, Ly], x = Lx
The temperature plot (at t=100s, x=0.5, y=*):
Files
Model report  tutorial1.xml 
Runtime model report  tutorial1rt.xml 
Source code  tutorial1.py 
This tutorial introduces the following concepts:
The model in this example is very similar to the model used in the tutorial 1. The differences are:
The temperature plot (at t=100s, x=0.5, y=*):
Files
Model report  tutorial2.xml 
Runtime model report  tutorial2rt.xml 
Source code  tutorial2.py 
This tutorial introduces the following concepts:
The model in this example is identical to the model used in the tutorial 1. Some additional equations that calculate the total flux at the bottom edge are added to illustrate the array functions.
The temperature plot (at t=100, x=0.5, y=*):
The average temperature plot (considering the whole xy domain):
Files
Model report  tutorial3.xml 
Runtime model report  tutorial3rt.xml 
Source code  tutorial3.py 
This tutorial introduces the following concepts:
In this example we model a very simple heat transfer problem where a small piece of copper is at one side exposed to the source of heat and at the other to the surroundings.
The lumped heat balance is given by the following equation:
rho * cp * dT/dt  Qin = h * A * (T  Tsurr)
where Qin is the power of the heater, h is the heat transfer coefficient, A is the surface area and Tsurr is the temperature of the surrounding air.
The process starts at the temperature of the metal of 283K. The metal is allowed to warm up for 200 seconds, when the heat source is removed and the metal cools down slowly to the ambient temperature.
This can be modelled using the following symmetrical state transition network:
IF t < 200
Qin = 1500 W
ELSE
Qin = 0 W
The temperature plot:
Files
Model report  tutorial4.xml 
Runtime model report  tutorial4rt.xml 
Source code  tutorial4.py 
This tutorial introduces the following concepts:
In this example we use the same heat transfer problem as in the tutorial 4. Again we have a piece of copper which is at one side exposed to the source of heat and at the other to the surroundings.
The process starts at the temperature of 283K. The metal is allowed to warm up, and then its temperature is kept in the interval (320K  340K) for 350 seconds. This is performed by switching the heater on when the temperature drops to 320K and by switching the heater off when the temperature reaches 340K. After 350s the heat source is permanently switched off and the metal is allowed to slowly cool down to the ambient temperature.
This can be modelled using the following nonsymmetrical state transition network:
STN Regulator
case Heating:
Qin = 1500 W
on condition T > 340K switch to Regulator.Cooling
on condition t > 350s switch to Regulator.HeaterOff
case Cooling:
Qin = 0 W
on condition T < 320K switch to Regulator.Heating
on condition t > 350s switch to Regulator.HeaterOff
case HeaterOff:
Qin = 0 W
The temperature plot:
Files
Model report  tutorial5.xml 
Runtime model report  tutorial5rt.xml 
Source code  tutorial5.py 
This tutorial introduces the following concepts:
A simple port type ‘portSimple’ is defined which contains only one variable ‘t’. Two models ‘modPortIn’ and ‘modPortOut’ are defined, each having one port of type ‘portSimple’. The wrapper model ‘modTutorial’ instantiate these two models as its units and connects them by connecting their ports.
Files
Model report  tutorial6.xml 
Runtime model report  tutorial6rt.xml 
Source code  tutorial6.py 
This tutorial introduces the following concepts:
In this example we use the same heat transfer problem as in the tutorial 4. The input power of the heater is defined as a variable. Since there is no equation defined to calculate the value of the input power, the system contains N variables but only N1 equations. To create a wellposed DAE system one of the variable needs to be “fixed”. However the choice of variables is not arbitrary and in this example the only variable that can be fixed is Qin. Thus, the Qin variable represents a degree of freedom (DOF). Its value will be fixed at the beginning of the simulation and later manipulated in the userdefined schedule in the overloaded function daeSimulation.Run().
The plot of the inlet power:
The temperature plot:
Files
Model report  tutorial7.xml 
Runtime model report  tutorial7rt.xml 
Source code  tutorial7.py 
This tutorial introduces the following concepts:
Some time it is not enough to send the results to the DAE Plotter but it is desirable to export them into a specified file format (i.e. for use in other programs). For that purpose, daetools provide a range of data reporters that save the simulation results in various formats. In adddition, daetools allow implementation of custom, userdefined data reporters. As an example, a userdefined data reporter is developed to save the results into a plain text file (after the simulation is finished). Obviously, the data can be processed in any other fashion. Moreover, in certain situation it is required to process the results in more than one way. The daeDelegateDataReporter can be used in those cases. It has the same interface and the functionality like all data reporters. However, it does not do any data processing itself but calls the corresponding functions of data reporters which are added to it using the function AddDataReporter. This way it is possible, at the same time, to send the results to the DAE Plotter and save them into a file (or process the data in some other ways). In this example the results are processed in 10 different ways at the same time.
The model used in this example is very similar to the model in the tutorials 4 and 5.
Files
Model report  tutorial8.xml 
Runtime model report  tutorial8rt.xml 
Source code  tutorial8.py 
This tutorial introduces the following concepts:
Currently there are the following linear equations solvers available:
In this example we use the same conduction problem as in the tutorial 1.
The temperature plot (at t=100s, x=0.5, y=*):
Files
Model report  tutorial9.xml 
Runtime model report  tutorial9rt.xml 
Source code  tutorial9.py 
This tutorial introduces the following concepts:
In this example we use the same conduction problem as in the tutorial 1.
Files
Model report  tutorial10.xml 
Runtime model report  tutorial10rt.xml 
Source code  tutorial10.py 
This tutorial describes the use of iterative linear solvers (AztecOO from the Trilinos project) with different preconditioners (builtin AztecOO, Ifpack or ML) and corresponding solver options. Also, the range of Trilins Amesos solver options are shown.
The model is very similar to the model in tutorial 1, except for the different boundary conditions and that the equations are written in a different way to maximise the number of items around the diagonal (creating the problem with the diagonally dominant matrix). These type of systems can be solved using very simple preconditioners such as Jacobi. To do so, the interoperability with the NumPy package has been exploited and the package itertools used to iterate through the distribution domains in x and y directions.
The equations are distributed in such a way that the following incidence matrix is obtained:
XXX 
 X X X 
 X X X 
 X X X 
 X X X 
 XXX 
 XXX 
 X XXX X 
 X XXX X 
 X XXX X 
 X XXX X 
 XXX 
 XXX 
 X XXX X 
 X XXX X 
 X XXX X 
 X XXX X 
 XXX 
 XXX 
 X XXX X 
 X XXX X 
 X XXX X 
 X XXX X 
 XXX 
 XXX 
 X XXX X 
 X XXX X 
 X XXX X 
 X XXX X 
 XXX 
 XXX 
 X X X 
 X X X 
 X X X 
 X X X 
 XXX
The temperature plot (at t=100s, x=0.5, y=*):
Files
Model report  tutorial11.xml 
Runtime model report  tutorial11rt.xml 
Source code  tutorial11.py 
This tutorial describes the use and available options for superLU direct linear solvers:
The model is the same as the model in tutorial 1, except for the different boundary conditions.
The temperature plot (at t=100s, x=0.5, y=*):
Files
Model report  tutorial12.xml 
Runtime model report  tutorial12rt.xml 
Source code  tutorial12.py 
This tutorial introduces the following concepts:
In this example we use the very similar model as in the tutorial 5.
The simulation output should show the following messages at t=100s and t=350s:
...
********************************************************
simpleUserAction2 message:
This message should be fired when the time is 100s.
********************************************************
...
********************************************************
simpleUserAction executed; input data = 427.464093129832
********************************************************
...
The plot of the ‘event’ variable:
The temperature plot:
Files
Model report  tutorial13.xml 
Runtime model report  tutorial13rt.xml 
Source code  tutorial13.py 
In this tutorial we introduce the external functions concept that can handle and execute functions in external libraries. The daeScalarExternalFunctionderived external function object is used to calculate the heat transferred and to interpolate a set of values using the scipy.interpolate.interp1d object.
In this example we use the same model as in the tutorial 5 with few additional equations.
The simulation output should show the following messages at the end of simulation:
...
scipy.interp1d statistics:
interp1d called 1703 times (cache value used 770 times)
The plot of the ‘Heat_ext’ variable:
The plot of the ‘Value_interp’ variable:
Files
Model report  tutorial14.xml 
Runtime model report  tutorial14rt.xml 
Source code  tutorial14.py 
This tutorial introduces the following concepts:
In this example we use the same model as in the tutorial 4 with the more complex STN:
IF t < 200
IF 0 <= t < 100
IF 0 <= t < 50
Qin = 1600 W
ELSE
Qin = 1500 W
ELSE
Qin = 1400 W
ELSE IF 200 <= t < 300
Qin = 1300 W
ELSE
Qin = 0 W
The plot of the ‘Qin’ variable:
The temperature plot:
Files
Model report  tutorial15.xml 
Runtime model report  tutorial15rt.xml 
Source code  tutorial15.py 
This tutorial shows how to use DAE Tools objects with NumPy arrays to solve a simple stationary heat conduction in one dimension using the Finite Elements method with linear elements and two ways of manually assembling a stiffness matrix/load vector:
d2T(x)/dx2 = F(x); x in (0, Lx)
Linear finite elements discretisation and simple FE matrix assembly:
phi phi
(k1) (k)
* *
*  * *  *
*  * *  *
*  * *  *
*  * *  *
*  *  *
*  * *  *
*  * *  *
*  * *  *
*  * element (k) *  *
**+++++++++++++++++++**
x x
(ki (k)
\_________ _________/

dx
The comparison of the analytical solution and two ways of assembling the system is given in the following plot:
Files
Model report  tutorial16.xml 
Runtime model report  tutorial16rt.xml 
Source code  tutorial16.py 
This tutorial introduces the following concepts:
In this example we use the same heat transfer problem as in the tutorial 7.
The screenshot of the TCP/IP log server:
The temperature plot:
Files
Model report  tutorial17.xml 
Runtime model report  tutorial17rt.xml 
Source code  tutorial17.py 
This tutorial shows one more problem solved using the NumPy arrays that operate on DAE Tools variables. The model is taken from the Sundials ARKODE (ark_analytic_sys.cpp). The ODE system is defined by the following system of equations:
dy/dt = A*y
where:
A = V * D * Vi
V = [1 1 1; 1 2 1; 0 1 2];
Vi = 0.25 * [5 1 3; 2 2 2; 1 1 1];
D = [0.5 0 0; 0 0.1 0; 0 0 lam];
lam is a large negative number.
The analytical solution to this problem is:
Y(t) = V*exp(D*t)*Vi*Y0
for t in the interval [0.0, 0.05], with initial condition y(0) = [1,1,1]’.
The stiffness of the problem is directly proportional to the value of “lamda”. The value of lamda should be negative to result in a wellposed ODE; for values with magnitude larger than 100 the problem becomes quite stiff.
In this example, we choose lamda = 100.
The solution:
lamda = 100
reltol = 1e06
abstol = 1e10

t y0 y1 y2

0.0050 0.70327 0.70627 0.41004
0.0100 0.52267 0.52865 0.05231
0.0150 0.41249 0.42145 0.16456
0.0200 0.34504 0.35696 0.29600
0.0250 0.30349 0.31838 0.37563
0.0300 0.27767 0.29551 0.42383
0.0350 0.26138 0.28216 0.45296
0.0400 0.25088 0.27459 0.47053
0.0450 0.24389 0.27053 0.48109
0.0500 0.23903 0.26858 0.48740

The plot of the ‘y0’, ‘y1’, ‘y2’ variables:
Files
Model report  tutorial18.xml 
Runtime model report  tutorial18rt.xml 
Source code  tutorial18.py 
This tutorial introduces the thermo physical property packages.
Since there are many thermo packages with a very different API the CapeOpen standard has been adopted in daetools. This way, all thermo packages implementing the CapeOpen thermo interfaces are automatically vailable to daetools. Those which do not are wrapped by the class with the CapeOpen conforming API. At the moment, two types of thermophysical property packages are implemented:
The central point is the daeThermoPhysicalPropertyPackage class. It can load any COM component that implements CapeOpen 1.1 ICapeThermoPropertyPackageManager interface or the CoolProp thermo package.
The framework provides lowlevel functions (specified in the CapeOpen standard) in the daeThermoPhysicalPropertyPackage class and the higherlevel functions in the auxiliary daeThermoPackage class defined in the daetools/pyDAE/thermo_packages.py file. The lowlevel functions are defined in the ICapeThermoCoumpounds and ICapeThermoPropertyRoutine CapeOpen interfaces. These functions come in two flavours:
The daeThermoPackage auxiliary class offers functions to calculate specified properties, for instance:
All functions return properties in the SI units (as specified in the CapeOpen 1.1 standard).
Known issues:
In this tutorial, we use a very simple model: a quantity of liquid (water + ethanol mixture) is heated using the constant input power. The model uses a thermo package to calculate the commonly used transport properties such as specific heat capacity, thermal conductivity, dynamic viscosity and binary diffusion coefficients. First, the lowlevel functions are tested for CapeOpen and CoolProp packages in the test_single_phase, test_two_phase, test_coolprop_single_phase functions. The results depend on the options selected in the CapeOpen package (equation of state, etc.). Then, the model that uses a thermo package is simulated.
The plot of the specific heat capacity as a function of temperature:
Files
Model report  tutorial19.xml 
Runtime model report  tutorial19rt.xml 
Source code  tutorial19.py 
This tutorial illustrates the sensitivity analysis features in DAE Tools.
This model has one state variable (T) and one degree of freedom (Qin). Qin is set as a parameter for sensitivity analysis.
The sensitivity analysis is enabled and the sensitivities can be reported to the data reporter like any ordinary variable by setting the boolean property simulation.ReportSensitivities to True.
Raw sensitivity matrices can be saved into a specified directory using the simulation.SensitivityDataDirectory property (before a call to Initialize). The sensitivity matrics are saved in .mmx coordinate format where the first dimensions is Nparameters and second Nvariables: S[Np, Nvars].
The plot of the sensitivity of T per Qin:
Files
Model report  tutorial20.xml 
Runtime model report  tutorial20rt.xml 
Source code  tutorial20.py 
This tutorial illustrates the additional sensitivity analysis features. Here, the numerical sensitivities for the Constant coefficient first order equations are compared to the available analytical solution.
The sensitivity analysis is enabled and the sensitivities are reported to the data reporter. The sensitivity data can be obtained in two ways:
The comparison between the numerical and the analytical sensitivities:
Files
Model report  tutorial21.xml 
Runtime model report  tutorial21rt.xml 
Source code  tutorial21.py 
This tutorial presents a userdefined simulation which instead of simply integrating the system shows the pyQt graphical user interface (GUI) where the simulation can be manipulated (a sort of interactive operating procedure).
The model in this example is the same as in the tutorial 4.
The simulation.Run() function is modifed to show the graphical user interface (GUI) that allows to specify the input power of the heater (degree of freedom), a time period for integration, and a reporting interval. The GUI also contains the temperature plot updated in real time, as the simulation progresses.
The screenshot of the pyQt GUI:
Files
Model report  tutorial_adv_1.xml 
Runtime model report  tutorial_adv_1rt.xml 
Source code  tutorial_adv_1.py 
This tutorial demonstrates a solution of a discretized population balance using high resolution upwind schemes with flux limiter.
Reference: Qamar S., Elsner M.P., Angelov I.A., Warnecke G., SeidelMorgenstern A. (2006) A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering 30(67):11191131. doi:10.1016/j.compchemeng.2006.02.012
It shows a comparison between the analytical results and various discretization schemes:
The problem is from the section 3.1 Sizeindependent growth.
Could be also found in: Motz S., Mitrović A., Gilles E.D. (2002) Comparison of numerical methods for the simulation of dispersed phase systems. Chemical Engineering Science 57(20):43294344. doi:10.1016/S00092509(02)003494
The comparison of number density functions between the analytical solution and the highresolution scheme:
The comparison of number density functions between the analytical solution and the I order upwind, II order upwind and II order central difference schemes:
Files
Model report  tutorial_adv_2.xml 
Runtime model report  tutorial_adv_2rt.xml 
Source code  tutorial_adv_2.py 
This tutorial introduces the following concepts:
The model represent a simple multiplier block. It contains two inlet and two outlet ports. The outlets values are equal to inputs values multiplied by a multiplier “m”:
out1.y = m1 x in1.y
out2.y[] = m2[] x in2.y[]
where multipliers m1 and m2[] are:
STN Multipliers
case variableMultipliers:
dm1/dt = p1
dm2[]/dt = p2
case constantMultipliers:
dm1/dt = 0
dm2[]/dt = 0
(that is the multipliers can be constant or variable).
The ports in1 and out1 are scalar (width = 1). The ports in2 and out2 are vectors (width = 1).
Achtung, Achtung!! Notate bene:
The plot of the inlet ‘y’ variable and the multiplied outlet ‘y’ variable for the constant multipliers (m1 = 2):
The plot of the inlet ‘y’ variable and the multiplied outlet ‘y’ variable for the variable multipliers (dm1/dt = 10, m1(t=0) = 2):
FMI CrossCheck results:
Start the DAEPlotter (or change the data reporter below).
Execute:
./fmuCheck.linux64 n 10 tutorial_adv_3.fmu
Results:
[INFO][FMUCHK] FMI compliance checker 2.0 [FMILibrary: 2.0] build date: Aug 22 2014
[INFO][FMUCHK] Will process FMU tutorial_adv_3.fmu
[INFO][FMILIB] XML specifies FMI standard version 2.0
[INFO][FMUCHK] Model name: tutorial_adv_3
[INFO][FMUCHK] Model GUID: e9654532099811e6957b9cb70d5dfdfc
[INFO][FMUCHK] Model version:
[INFO][FMUCHK] FMU kind: CoSimulation
[INFO][FMUCHK] The FMU contains:
0 constants
3 parameters
0 discrete variables
4 continuous variables
2 inputs
2 outputs
0 local variables
0 independent variables
0 calculated parameters
6 real variables
0 integer variables
0 enumeration variables
0 boolean variables
1 string variables
[INFO][FMUCHK] Printing output file header
time,out_1.y,out_2.y
[INFO][FMUCHK] Model identifier for CoSimulation: tutorial_adv_3
[INFO][FMILIB] Loading 'linux64' binary with 'default' platform types
[INFO][FMUCHK] Version returned from CS FMU: 2.0
***********************************************************************
Version: 1.5.0
Copyright: Dragan Nikolic
Homepage: http://www.daetools.com
@ @
@ @@@@@ @@@@@ @@@@@ @@@@@ @@@@@ @ @@@@@
@@@@@@ @ @ @ @ @ @ @ @ @ @
@ @ @@@@@@ @@@@@@ @ @ @ @ @ @ @@@@@
@ @ @ @ @ @ @ @ @ @ @ @
@@@@@@ @@@@@@ @@@@@ @@@ @@@@@ @@@@@ @@@@ @@@@@
***********************************************************************
DAE Tools is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License version 3
as published by the Free Software Foundation.
DAE Tools is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses>.
***********************************************************************
Creating the system...
The system created successfully in: 0.002 s
Starting the initialization of the system... Done.
[INFO][FMUCHK] Initialized FMU for simulation starting at time 0
0,1.0000000010000001E+00,2.0000000039999999E+00
10,1.0100000000000001E+02,4.0200000000000006E+02
20,2.0099999999999994E+02,8.0199999999999977E+02
30,3.0100000000000000E+02,1.2020000000000000E+03
40,4.0100000000000000E+02,1.6020000000000000E+03
50,5.0099999999999989E+02,2.0019999999999995E+03
60,6.0100000000000000E+02,2.4020000000000000E+03
70,7.0100000000000000E+02,2.8020000000000000E+03
80,8.0100000000000000E+02,3.2020000000000000E+03
90,9.0099999999999977E+02,3.6019999999999991E+03
100,1.0009999999999998E+03,4.0019999999999991E+03
[INFO][FMUCHK] Simulation finished successfully at time 100
FMU check summary:
FMU reported:
0 warning(s) and error(s)
Checker reported:
0 Warning(s)
0 Error(s)
Files
Model report  tutorial_adv_3.xml 
Runtime model report  tutorial_adv_3rt.xml 
Source code  tutorial_adv_3.py 
This tutorial illustrates the C++ MPI code generator. The model is identical to the model in the tutorial 11.
The temperature plot (at t=100s, x=0.5128, y=*):
Files
Model report  tutorial_adv_4.xml 
Runtime model report  tutorial_adv_4rt.xml 
Source code  tutorial_adv_4.py 
An introductory example of the support for Finite Elements in daetools. The basic idea is to use an external library to perform all lowlevel tasks such as management of mesh elements, degrees of freedom, matrix assembly, management of boundary conditions etc. deal.II library (www.dealii.org) is employed for these tasks. The mass and stiffness matrices and the load vector assembled in deal.II library are used to generate a set of algebraic/differential equations in the following form: [Mij]{dx/dt} + [Aij]{x} = {Fi}. Specification of additional equations such as surface/volume integrals are also available. The numerical solution of the resulting ODA/DAE system is performed in daetools together with the rest of the model equations.
The unique feature of this approach is a capability to use daetools variables to specify boundary conditions, time varying coefficients and nonlinear terms, and evaluate quantities such as surface/volume integrals. This way, the finite element model is fully integrated with the rest of the model and multiple FE systems can be created and coupled together. In addition, nonlinear and DAE finite element systems are automatically supported.
In this tutorial the simple transient heat conduction problem is solved using the finite element method:
dT/dt  kappa/(rho*cp)*nabla^2(T) = g(T) in Omega
The mesh is rectangular with two holes, similar to the mesh in step49 deal.II example:
Dirichlet boundary conditions are set to 300 K on the outer rectangle, 350 K on the inner ellipse and 250 K on the inner diamond.
The temperature plot at t = 500s (generated in VisIt):
Files
Model report  tutorial_dealii_1.xml 
Runtime model report  tutorial_dealii_1rt.xml 
Source code  tutorial_dealii_1.py 
In this example a simple transient heat convectiondiffusion equation is solved.
dT/dt  kappa/(rho*cp)*nabla^2(T) + nabla.(uT) = g(T) in Omega
The fluid flows from the left side to the right with constant velocity of 0.01 m/s. The inlet temperature for 0.2 <= y <= 0.3 is iven by the following expression:
T_left = T_base + T_offset*sin(pi*t/25) on dOmega
creating a bubblelike regions of higher temperature that flow towards the right end and slowly diffuse into the bulk flow of the fluid due to the heat conduction.
The mesh is rectangular with the refined elements close to the left/right ends:
The temperature plot at t = 500s:
Files
Model report  tutorial_dealii_2.xml 
Runtime model report  tutorial_dealii_2rt.xml 
Source code  tutorial_dealii_2.py 
In this example the CahnHilliard equation is solved using the finite element method. This equation describes the process of phase separation, where two components of a binary mixture separate and form domains pure in each component.
dc/dt  D*nabla^2(mu) = 0, in Omega
mu = c^3  c  gamma*nabla^2(c)
The mesh is a simple square (0100)x(0100):
The concentration plot at t = 500s:
Files
Model report  tutorial_dealii_3.xml 
Runtime model report  tutorial_dealii_3rt.xml 
Source code  tutorial_dealii_3.py 
In this tutorial the transient heat conduction problem is solved using the finite element method:
dT/dt  kappa * nabla^2(Τ) = g in Omega
The mesh is a pipe submerged into water receiving the heat of the sun at the 45 degrees from the topleft direction, the heat from the suroundings and having the constant temperature of the inner tube. The boundary conditions are thus:
The pipe mesh is given below:
The temperature plot at t = 3600s:
Files
Model report  tutorial_dealii_4.xml 
Runtime model report  tutorial_dealii_4rt.xml 
Source code  tutorial_dealii_4.py 
In this example a simple flow through porous media is solved (deal.II step20).
K^{1} u + nabla(p) = 0 in Omega
div(u) = f in Omega
p = g on dOmega
The mesh is a simple square:
The velocity plot at t = 500s:
Files
Model report  tutorial_dealii_5.xml 
Runtime model report  tutorial_dealii_5rt.xml 
Source code  tutorial_dealii_5.py 
A simple steadystate diffusion and firstorder reaction in an irregular catalyst shape (Proc. 6th Int. Conf. on Mathematical Modelling, Math. Comput. Modelling, Vol. 11, 375319, 1988) applying Dirichlet and Robin type of boundary conditions.
D_eA * nabla^2(C_A)  k_r * C_A = 0 in Omega
D_eA * nabla(C_A) = k_m * (C_A  C_Ab) on dOmega1
C_A = C_Ab on dOmega2
The catalyst pellet mesh:
The concentration plot:
The concentration plot for Ca=Cab on all boundaries:
Files
Model report  tutorial_dealii_6.xml 
Runtime model report  tutorial_dealii_6rt.xml 
Source code  tutorial_dealii_6.py 
Continuously Stirred Tank Reactor with energy balance and Van de Vusse reactions:
A > B > C
2A > D
Reference: G.A. Ridlehoover, R.C. Seagrave. Optimization of Van de Vusse Reaction Kinetics Using Semibatch Reactor Operation, Ind. Eng. Chem. Fundamen. 1973;12(4):444447. doi:10.1021/i160048a700
The concentrations plot:
The temperatures plot:
Files
Model report  tutorial_che_1.xml 
Runtime model report  tutorial_che_1rt.xml 
Source code  tutorial_che_1.py 
Binary distillation column model.
Reference: J. Hahn, T.F. Edgar. An improved method for nonlinear model reduction using balancing of empirical gramians. Computers and Chemical Engineering 2002; 26:13791397. doi:10.1016/S00981354(02)001205
The liquid fraction after 120 min (x(reboiler)=0.935420, x(condenser)=0.064581):
The liquid fraction in the reboiler (tray 1) and in the condenser (tray 32):
Files
Model report  tutorial_che_2.xml 
Runtime model report  tutorial_che_2rt.xml 
Source code  tutorial_che_2.py 
Batch reactor seeded crystallisation using the method of moments.
References (model equations and input parameters):
The main assumptions:
Solubility of Paracetamol in ethanol:

Temperature, C Solubility, kg Parac./kg EtOH Solubility, mol Parac./m3 EtOH

0 0.11362 593.0387
10 0.14128 737.4215
20 0.17568 916.9562
30 0.21845 1140.2008
40 0.27163 1417.7972
50 0.33777 1762.9779
60 0.42000 2192.1973

The supersaturation plot:
The concentration plot:
The recovery plot:
The yield plot:
The total number of crystals plot:
Files
Model report  tutorial_che_3.xml 
Runtime model report  tutorial_che_3rt.xml 
Source code  tutorial_che_3.py 
This example shows a comparison between the analytical results and the discretised population balance equations results solved using the cell centered finite volume method employing the high resolution upwind scheme (Van Leer kinterpolation with k = 1/3) and a range of flux limiters.
This tutorial can be run from the console only.
The problem is from the section 4.1.1 Sizeindependent growth I of the following article:
and also from the section 3.1 Sizeindependent growth of the following article:
The growthonly crystallisation process was considered with the constant growth rate of 1μm/s and the following initial number density function:
n(L,0): 1E10, if 10μm < L < 20μm
0, otherwise
The crystal size in the range of [0, 100]μm was discretised into 100 elements. The analytical solution in this case is equal to the initial profile translated right in time by a distance Gt (the growth rate multiplied by the time elapsed in the process).
The flux limiters used in the model are:
Comparison of L1 and L2norms (ni_HR  ni_analytical):

Scheme L1 L2

superbee 1.786e+10 7.016e+09
Sweby 2.817e+10 8.614e+09
Koren 3.015e+10 9.293e+09
smart 2.961e+10 9.326e+09
MC 3.258e+10 9.807e+09
HCUS 3.638e+10 1.001e+10
HQUICK 3.622e+10 1.005e+10
vanLeerMinmod 3.581e+10 1.011e+10
vanLeer 3.874e+10 1.059e+10
ospre 4.139e+10 1.094e+10
UMIST 4.363e+10 1.136e+10
Osher 4.579e+10 1.156e+10
vanAlbada1 4.574e+10 1.157e+10
minmod 5.653e+10 1.325e+10
vanAlbada2 5.456e+10 1.331e+10

The comparison of number density functions between the analytical solution and the solution obtained using highresolution scheme with the Superbee flux limiter at t=60s:
The comparison of number density functions between the analytical solution and the solution obtained using highresolution scheme with the Koren flux limiter at t=60s:
Files
Source code  tutorial_che_4.py 
Analytical solution  fl_analytical.py 
Flux limiters  flux_limiters.py 
Similar to the chem. eng. example 4, this example shows a comparison between the analytical results and the discretised population balance equations results solved using the cell centered finite volume method employing the high resolution upwind scheme (Van Leer kinterpolation with k = 1/3) and a range of flux limiters.
This tutorial can be run from the console only.
The problem is from the section 4.1.2 Sizeindependent growth II of the following article:
and also from the section 3.2 Sizeindependent growth of the following article:
Again, the growthonly crystallisation process was considered with the constant growth rate of 0.1μm/s and with the different initial number density function:
n(L,0): 0, if L <= 2.0μm
1E10, if 2μm < L <= 10μm (region I)
0, if 10μm < L <= 18μm
1E10*cos^2(pi*(L26)/64), if 18μm < L <= 34μm (region II)
0, if 34μm < L <= 42μm
1E10*sqrt(1(L50)^2/64), if 42μm < L <= 58μm (region III)
0, if 58μm < L <= 66μm
1E10*exp((L70)^2/(2sigma^2)), if 66μm < L <= 74μm (region IV)
0, if 74μm < L
The crystal size in the range of [0, 100]μm was discretised into 200 elements. The analytical solution in this case is equal to the initial profile translated right in time by a distance Gt (the growth rate multiplied by the time elapsed in the process).
Comparison of L1 and L2norms (ni_HR  ni_analytical):

Scheme L1 L2

superbee 4.464e+10 1.015e+10
smart 4.727e+10 1.120e+10
Koren 4.861e+10 1.141e+10
Sweby 5.435e+10 1.142e+10
MC 5.129e+10 1.162e+10
HQUICK 5.531e+10 1.194e+10
HCUS 5.528e+10 1.194e+10
vanLeerMinmod 5.600e+10 1.202e+10
vanLeer 5.814e+10 1.225e+10
ospre 6.131e+10 1.252e+10
UMIST 6.181e+10 1.259e+10
Osher 6.690e+10 1.275e+10
vanAlbada1 6.600e+10 1.281e+10
minmod 7.751e+10 1.360e+10
vanAlbada2 7.901e+10 1.413e+10

The comparison of number density functions between the analytical solution and the solution obtained using highresolution scheme with the Superbee flux limiter at t=100s:
The comparison of number density functions between the analytical solution and the solution obtained using highresolution scheme with the Koren flux limiter at t=100s:
Files
Source code  tutorial_che_5.py 
Analytical solution  fl_analytical.py 
Flux limiters  flux_limiters.py 
Model of a lithiumion battery based on porous electrode theory as developed by John Newman and coworkers. In particular, the equations here are based on a summary of the methodology by Karen E. Thomas, John Newman, and Robert M. Darling,
Thomas K., Newman J., Darling R. (2002). Mathematical Modeling of Lithium Batteries in Advances in Lithiumion Batteries. Springer US. 345392. doi:10.1007/0306475081_13
A few simplifications have been made rather than implementing the more complete model described there. For example, the following assumptions have (currently) been made:
The up to date version of the model is available at Raymond’s GitHub repository: https://github.com/raybsmith/daetoolsexamplebattery.
The voltage plot:
The current plot:
Files
Model report  tutorial_che_6.xml 
Runtime model report  tutorial_che_6rt.xml 
Source code  tutorial_che_6.py 
Steadystate Plug Flow Reactor with energy balance and first order reaction:
A > B
The problem is example 9.4.3 from the section 9.4 Nonisothermal Plug Flow Reactor from the following book:
The dimensionless concentration plot:
The dimensionless temperature plot (adiabatic and nonisothermal cases):
Files
Model report  tutorial_che_7.xml 
Runtime model report  tutorial_che_7rt.xml 
Source code  tutorial_che_7.py 
Model of a gas separation on a porous membrane with a metal support. The model employs the Generalised MaxwellStefan (GMS) equations to predict fluxes and selectivities. The membrane unit model represents a generic twodimensonal model of a porous membrane and consists of four models:
The retentate compartment, the porous membrane, the support layer and the permeate compartment are coupled via molar flux, temperature, pressure and gas composition at the interfaces. The model is described in the section 2.2 Membrane modelling of the following article:
and in the original Krishna article:
This version is somewhat simplified for it only offers an extended Langmuir isotherm. The Ideal Adsorption Solution theory (IAS) and the Real Adsorption Solution theory (RAS) described in the articles are not implemented here.
The problem modelled is separation of hydrocarbons (CH4+C2H6) mixture on a zeolite (silicalite1) membrane with a metal support from the section ‘Binary mixture permeation’ of the following article:
The CH4 and C2H6 fluxes, and CH4/C2H6 selectivity plots for two cases: GMS and GMS(Dij=∞), 1:1 mixture, and T = 303 K:
Files
Model report  tutorial_che_8.xml 
Runtime model report  tutorial_che_8rt.xml 
Source code  tutorial_che_8.py 
Membrane unit  membrane_unit.py 
Variable types  membrane_variable_types.py 
Membrane model  membrane.py 
Support model  support.py 
In/out compartment  compartment.py 
Chemical reaction network from the Dow Chemical Company described in the following article:
The sensitivity analysis is enabled and the sensitivities are reported to the data reporter. The sensitivity data can be obtained in two ways:
The concentrations plot (u1, u3, u4):
The concentrations plot (u6, u8):
The sensitivities plot (k2*du1/dk2, k2*du2/dk2, k2*du3/dk2, k2*du4/dk2, k2*du5/dk2):
Files
Model report  tutorial_che_9.xml 
Runtime model report  tutorial_che_9rt.xml 
Source code  tutorial_che_9.py 
The implementations of the COPS tests differ from the original ones in following:
As a consequence, the results slightly differ from the published results. In addition, the solver options should be tuned to achieve faster convergence.
Optimisation of the CSTR model and Van de Vusse reactions given in tutorial_che_1:
Not fully implemented yet.
Reference: G.A. Ridlehoover, R.C. Seagrave. Optimization of Van de Vusse Reaction Kinetics Using Semibatch Reactor Operation, Ind. Eng. Chem. Fundamen. 1973;12(4):444447. doi:10.1021/i160048a700
Files
Model report  tutorial_che_opt_1.xml 
Runtime model report  tutorial_che_opt_1rt.xml 
Source code  tutorial_che_opt_1.py 
COPS test 5: Isomerization of αpinene (parameter estimation of a dynamic system).
Very slow convergence.
Determine the reaction coefficients in the thermal isometrization of αpinene (y1) to dipentene (y2) and alloocimen (y3) which in turn produces α and βpyronene (y4) and a dimer (y5).
Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS273, 2004. PDF
Experimental data taken from: Rocha A.M.A.C., Martins M.C., Costa M.F.P., Fernandes, E.M.G.P. (2016) Direct sequential based firefly algorithm for the αpinene isomerization problem. 16th International Conference on Computational Science and Its Applications, ICCSA 2016, Beijing, China. doi:10.1007/9783319420851_30
Run options:
Currently, the parameter estimation results are (solver options/scaling should be tuned):
Fobj 57.83097
p1 5.63514e05
p2 2.89711e05
p3 1.39979e05
p4 18.67874e05
p5 2.23770e05
The concentration plots (for optimal ‘p’ from the literature):
Files
Model report  tutorial_che_opt_2.xml 
Runtime model report  tutorial_che_opt_2rt.xml 
Source code  tutorial_che_opt_2.py 
COPS test 6: Marine Population Dynamics. (Not working properly)
Given estimates of the abundance of the population of a marine species at each stage (for example, nauplius, juvenile, adult) as a function of time, determine stage specific growth and mortality rates.
Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS273, 2004. PDF
Experimental data generated following the procedure described in the COPS test.
Run options:
Currently, the parameter estimation results are (suboptimal results obtained, solver options/scaling should be tuned):
Fobj = 1.920139e+8
m(0) 3.358765e01
m(1) 4.711709e01
m(2) 1.120524e01
m(3) 8.509170e02
m(4) 9.683579e02
m(5) 1.919142e01
m(6) 2.418778e01
m(7) 2.421000e01
g(0) 1.152995e+00
g(1) 7.529383e01
g(2) 5.024174e01
g(3) 5.704327e01
g(4) 4.180333e01
g(5) 3.185407e01
g(6) 2.250250e01
The distribution moments 1,2,5,6 plots (for optimal results from the literature):
The distribution moments 3,4,7,8 plots (for optimal results from the literature):
Files
Model report  tutorial_che_opt_3.xml 
Runtime model report  tutorial_che_opt_3rt.xml 
Source code  tutorial_che_opt_3.py 
COPS test 12: Catalytic Cracking of Gas Oil.
Determine the reaction coefficients for the catalytic cracking of gas oil into gas and other byproducts.
Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS273, 2004. PDF
Experimental data generated following the procedure described in the COPS test.
Run options:
Currently, the parameter estimation results are (solver options/scaling should be tuned):
Fobj = 4.841995e3
p1 = 10.95289
p2 = 7.70601
p3 = 2.89625
The concentration plots (for optimal ‘p’ from the literature):
Files
Model report  tutorial_che_opt_4.xml 
Runtime model report  tutorial_che_opt_4rt.xml 
Source code  tutorial_che_opt_4.py 
COPS test 13: Methanol to Hydrocarbons.
Determine the reaction coefficients for the conversion of methanol into various hydrocarbons.
Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS273, 2004. PDF
Experimental data generated following the procedure described in the COPS test.
Run options:
Currently, the parameter estimation results are (solver options/scaling should be tuned):
Fobj = 1.274997e2
p1 = 2.641769
p2 = 1.466245
p3 = 1.884254
p4 = 1.023885
p5 = 0.471067
The concentration plots (for optimal ‘p’ from the literature):
The concentration plots (for optimal ‘p’ from this optimisation):
Files
Model report  tutorial_che_opt_5.xml 
Runtime model report  tutorial_che_opt_5rt.xml 
Source code  tutorial_che_opt_5.py 
COPS optimisation test 14: Catalyst Mixing.
Determine the optimal mixing policy of two catalysts along the length of a tubular plug flow reactor involving several reactions.
Reference: Benchmarking Optimization Software with COPS 3.0, Mathematics and Computer Science Division, Argonne National Laboratory, Technical Report ANL/MCS273, 2004. PDF
In DAE Tools numerical solution of dynamic optimisation problems is obtained using the Direct Sequential Approach where, given a set of values for the decision variables, the system of ODEs are accurately integrated over the entire time interval using specific numerical integration formulae so that the objective functional can be evaluated. Therefore, the differential equations are satisfied at each iteration of the optimisation procedure.
In the COPS test, the problem is solved using the Direct Simultaneous Approach where the equations that result from a discretisation of an ODE model using orthogonal collocation on finite elements (OCFE), are incorporated directly into the optimisation problem, and the combined problem is then solved using a largescale optimisation strategy.
The results: fobj = 4.79479E2 (for Nh = 100) and 4.78676E02 (for Nh = 200).
The control variables plot (for Nh = 100):
The control variables plot (for Nh = 200):
Files
Model report  tutorial_che_opt_6.xml 
Runtime model report  tutorial_che_opt_6rt.xml 
Source code  tutorial_che_opt_6.py 
This tutorial introduces IPOPT NLP solver, its setup and options.
Files
Model report  opt_tutorial1.xml 
Runtime model report  opt_tutorial1rt.xml 
Source code  opt_tutorial1.py 
This tutorial introduces Bonmin MINLP solver, its setup and options.
Files
Model report  opt_tutorial2.xml 
Runtime model report  opt_tutorial2rt.xml 
Source code  opt_tutorial2.py 
This tutorial introduces NLOPT NLP solver, its setup and options.
Files
Model report  opt_tutorial3.xml 
Runtime model report  opt_tutorial3rt.xml 
Source code  opt_tutorial3.py 
This tutorial shows the interoperability between DAE Tools and 3rd party optimization software (scipy.optimize) used to minimize the Rosenbrock function.
DAE Tools simulation is used to calculate the objective function and its gradients, while scipy.optimize.fmin function (NelderMead Simplex algorithm) to find the minimum of the Rosenbrock function.
Files
Model report  opt_tutorial4.xml 
Runtime model report  opt_tutorial4rt.xml 
Source code  opt_tutorial4.py 
This tutorial shows the interoperability between DAE Tools and 3rd party optimization software (scipy.optimize) used to fit the simple function with experimental data.
DAE Tools simulation object is used to calculate the objective function and its gradients, while scipy.optimize.leastsq function (a wrapper around MINPACK’s lmdif and lmder) implementing LevenbergMarquardt algorithm is used to estimate the parameters.
Files
Model report  opt_tutorial5.xml 
Runtime model report  opt_tutorial5rt.xml 
Source code  opt_tutorial5.py 
daeMinpackLeastSq module test.
Files
Model report  opt_tutorial6.xml 
Runtime model report  opt_tutorial6rt.xml 
Source code  opt_tutorial6.py 
This tutorial introduces monitoring optimization progress.
Files
Model report  opt_tutorial7.xml 
Runtime model report  opt_tutorial7rt.xml 
Source code  opt_tutorial7.py 