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 default daeSimulation.Run() function (reimplemented in Python) is:
def Run(self):
# Python implementation of daeSimulation::Run() C++ function.
import math
while self.CurrentTime < self.TimeHorizon:
# Get the time step (based on the TimeHorizon and the ReportingInterval).
# Do not allow to get past the TimeHorizon.
t = self.NextReportingTime
if t > self.TimeHorizon:
t = self.TimeHorizon
# If the flag is set  a user tries to pause the simulation, therefore return.
if self.ActivityAction == ePauseActivity:
self.Log.Message("Activity paused by the user", 0)
return
# If a discontinuity is found, loop until the end of the integration period.
# The data will be reported around discontinuities!
while t > self.CurrentTime:
self.Log.Message("Integrating from [%f] to [%f] ..." % (self.CurrentTime, t), 0)
self.IntegrateUntilTime(t, eStopAtModelDiscontinuity, True)
# After the integration period, report the data.
self.ReportData(self.CurrentTime)
# Set the simulation progress.
newProgress = math.ceil(100.0 * self.CurrentTime / self.TimeHorizon)
if newProgress > self.Log.Progress:
self.Log.Progress = newProgress
In this example the following schedule is specified:
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 addition, functions defined in shared libraries (.so in GNU/Linux, .dll in Windows and .dylib in macOS) can be used via ctypes Python library and daeCTypesExternalFunction class.
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_ext1’ variable:
The plot of the ‘Heat_ext2’ variable:
The plot of the ‘Value_interp’ variable:
Files
Model report  tutorial14.xml 
Runtime model report  tutorial14rt.xml 
Source code  tutorial14.py 
External function source code  tutorial14_heat_function.c 
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 support variable constraints available in Sundials IDA solver. Benchmarks are available from Matlab documentation.
Absolute Value Function:
dy/dt = fabs(y)
solved on the interval [0,40] with the initial condition y(0) = 1. The solution of this ODE decays to zero. If the solver produces a negative solution value, the computation eventually will fail as the calculated solution diverges to inf. Using the constraint y >= 0 resolves this problem.
The Knee problem:
epsilon * dy/dt = (1t)*y  y**2
solved on the interval [0,2] with the initial condition y(0) = 1. The parameter epsilon is 0 < epsilon << 1 and in this example equal to 1e6. The solution follows the y = 1x isocline for the whole interval of integration which is incorrect. Using the constraint y >= 0 resolves the problem.
In DAE Tools contraints follow the Sundials IDA solver implementation and can be specified using the valueConstraint argument of the daeVariableType class __init__ function:
and changed for individual variables using daeVariable.SetValueConstraint functions.
Absolute Value Function solution plot:
The Knee problem solution plot:
Files
Source code  tutorial20.py 
This tutorial introduces different methods for evaluation of equations in parallel. Equations residuals, Jacobian matrix and sensitivity residuals can be evaluated in parallel using two methods
The Evaluation Tree approach (default)
OpenMP API is used for evaluation in parallel. This method is specified by setting daetools.core.equations.evaluationMode option in daetools.cfg to “evaluationTree_OpenMP” or setting the simulation property:
simulation.EvaluationMode = eEvaluationTree_OpenMP
numThreads controls the number of OpenMP threads in a team. If numThreads is 0 the default number of threads is used (the number of cores in the system). Sequential evaluation is achieved by setting numThreads to 1.
The Compute Stack approach
Equations can be evaluated in parallel using:
OpenMP API for general purpose processors and manycore devices (i.e. Xeon Phi).
This method is specified by setting daetools.core.equations.evaluationMode option in daetools.cfg to “computeStack_OpenMP” or setting the simulation property:
simulation.EvaluationMode = eComputeStack_OpenMP
numThreads controls the number of OpenMP threads in a team. If numThreads is 0 the default number of threads is used (the number of cores in the system). Sequential evaluation is achieved by setting numThreads to 1.
OpenCL framework for streaming processors (GPU, GPGA) and heterogeneous systems.
This type is implemented in an external Python module pyEvaluator_OpenCL. It is up to one order of magnitude faster than the Evaluation Tree approach. However, it does not support external functions nor thermophysical packages.
OpenCL evaluators can use a single or multiple OpenCL devices. It is required to install OpenCL drivers/runtime libraries. InteL: https://software.intel.com/enus/articles/opencldrivers AMD: https://support.amd.com/enus/kbarticles/Pages/OpenCL2Driver.aspx NVidia: https://developer.nvidia.com/opencl
Files
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):
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 
Code verification using the Method of Exact Solutions.
Here, the numerical solution and 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  tutorial_cv_1.xml 
Runtime model report  tutorial_cv_1rt.xml 
Source code  tutorial_cv_1.py 
Code verification using the Method of Manufactured Solutions.
References:
The procedure for the transient convectiondiffusion (Burger’s) equation:
L(f) = df/dt + f*df/dx  D*d2f/dx2 = 0
is the following:
Pick a function (q, the manufactured solution):
q = A + sin(x + Ct)
Compute the new source term (g) for the original problem:
g = dq/dt + q*dq/dx  D*d2q/dx2
Solve the original problem with the new source term:
df/dt + f*df/dx  D*d2f/dx2 = g
Since L(f) = g and g = L(q), consequently we have: f = q. Therefore, the computed numerical solution (f) should be equal to the manufactured one (q).
The terms in the source g term are:
dq/dt = C * cos(x + C*t)
dq/dx = cos(x + C*t)
d2q/dx2 = sin(x + C*t)
The equations are solved for Dirichlet boundary conditions:
f(x=0) = q(x=0) = A + sin(0 + Ct)
f(x=2pi) = q(x=2pi) = A + sin(2pi + Ct)
Numerical vs. manufactured solution plot (no. elements = 60, t = 1.0s):
The normalised global errors and the order of accuracy plots (no. elements = [60, 90, 120, 150], t = 1.0s):
Files
Model report  tutorial_cv_2.xml 
Runtime model report  tutorial_cv_2rt.xml 
Source code  tutorial_cv_2.py 
Code verification using the Method of Manufactured Solutions.
References:
The problem in this tutorial is identical to tutorial_cv_3. The only difference is that the Neumann boundary conditions are applied:
df(x=0)/dx = dq(x=0)/dx = cos(0 + Ct)
df(x=2pi)/dx = dq(x=2pi)/dx = cos(2pi + Ct)
Numerical vs. manufactured solution plot (no. elements = 60, t = 1.0s):
The normalised global errors and the order of accuracy plots (no. elements = [60, 90, 120, 150], t = 1.0s):
Files
Model report  tutorial_cv_3.xml 
Runtime model report  tutorial_cv_3rt.xml 
Source code  tutorial_cv_3.py 
Code verification using the Method of Manufactured Solutions.
Reference:
The problem in this tutorial is the transient convectiondiffusion equation distributed on a rectangular 2D domain with u and v components of velocity:
L(u) = du/dt + (d(uu)/dx + d(uv)/dy)  ni * (d2u/dx2 + d2u/dy2) = 0
L(v) = dv/dt + (d(vu)/dx + d(vv)/dy)  ni * (d2v/dx2 + d2v/dy2) = 0
The manufactured solutions are:
um = u0 * (sin(x**2 + y**2 + w0*t) + eps)
vm = v0 * (cos(x**2 + y**2 + w0*t) + eps)
The terms in the new sources Su and Sv are computed using the daetools derivative functions (dt, d and d2).
Again, the Dirichlet boundary conditions are used:
u(LB, y) = um(LB, y)
u(UB, y) = um(UB, y)
u(x, LB) = um(x, LB)
u(x, UB) = um(x, UB)
v(LB, y) = vm(LB, y)
v(UB, y) = vm(UB, y)
v(x, LB) = vm(x, LB)
v(x, UB) = vm(x, UB)
Steadystate results (w0 = 0.0) for the ucomponent:
Steadystate results (w0 = 0.0) for the vcomponent:
Numerical vs. manufactured solution plot (u velocity component, 40x32 grid):
The normalised global errors and the order of accuracy plots (grids 10x8, 20x16, 40x32, 80x64):
Files
Model report  tutorial_cv_4.xml 
Runtime model report  tutorial_cv_4rt.xml 
Source code  tutorial_cv_4.py 
Code verification using the Method of Manufactured Solutions.
This problem and its solution in COMSOL Multiphysics software is described in the COMSOL blog: Verify Simulations with the Method of Manufactured Solutions (2015).
Here, a 1D transient heat conduction problem in a bar of length L is solved using the FE method:
dT/dt  k/(rho*cp) * d2T/dx2 = 0, x in [0,L]
with the following boundary:
T(0,t) = 500 K
T(L,t) = 500 K
and initial conditions:
T(x,0) = 500 K
The manufactured solution is given by function u(x):
u(x) = 500 + (x/L) * (x/L  1) * (t/tau)
The new source term is:
g(x) = du/dt  k/(rho*cp) * d2u/dx2
The terms in the source g term are:
du_dt = (x/L) * (x/L  1) * (1/tau)
d2u_dx2 = (2/(L**2)) * (t/tau)
Finally, the original problem with the new source term is:
dT/dt  k/(rho*cp) * d2T/dx2 = g(x), x in [0,L]
The mesh is linear (a bar) with a length of 100 m:
The comparison plots for the coarse mesh and linear elements:
The comparison plots for the coarse mesh and quadratic elements:
The comparison plots for the fine mesh and linear elements:
The comparison plots for the fine mesh and quadratic elements:
Files
Model report  tutorial_cv_5.xml 
Runtime model report  tutorial_cv_5rt.xml 
Source code  tutorial_cv_5.py 
Code verification using the Method of Exact Solutions.
Reference (section 3.3):
The problem in this tutorial is 1D transient convectiondiffusion equation:
dc_dt + u*dc/dx  D*d2c/dc2 = 0
The equation is solved using the high resolution cellcentered finite volume upwind scheme with flux limiter described in the article.
Numerical vs. exact solution plots (Nx = [20, 40, 80]):
Files
Model report  tutorial_cv_6.xml 
Runtime model report  tutorial_cv_6rt.xml 
Source code  tutorial_cv_6.py 
Code verification using the Method of Manufactured Solutions.
Reference (section 4.2):
The problem in this tutorial is 1D steadystate convectiondiffusion (Burger’s) equation:
u*dc/dx  D*d2c/dc2 = 0
The manufactured solution is:
c(x) = 0.5 * (1  cos(2*pi*(xa)/(ba))), x in [a,b]
c(x) = 0, otherwise
The new source term is:
s(x) = pi/(ba) * u * sin(2*pi*(xa)/(ba))  2*(pi/(ba))**2 * D * cos(2*pi*(xa)/(ba)), x in [a,b]
s(x) = 0, otherwise
The modified equation:
u*dc/dx  D*d2c/dc2 = s(x)
is solved using the high resolution cellcentered finite volume upwind scheme with flux limiter described in the article.
In order to obtain the consistent discretisation of the convection and the source terms an integral of the source term: S(x) = 1/u * Integral s(x)*dx must be calculated. The result of integration is given as:
S(x) = 0.5 * (1  cos(2*pi*(xa)/(ba))))  pi/(ba) * D/u * sin(2*pi*(xa)/(ba))), x in [a,b]
S(x) = 0, otherwise
Numerical vs. manufactured solution plot (Nx=80):
The normalised global errors and the order of accuracy plots for the Koren flux limiter (grids 40, 60, 80, 120):
The normalised global errors and the order of accuracy plots for the Superbee flux limiter (grids 40, 60, 80, 120):
Files
Model report  tutorial_cv_7.xml 
Runtime model report  tutorial_cv_7rt.xml 
Source code  tutorial_cv_7.py 
Code verification using the Method of Manufactured Solutions.
Reference (page 64):
The problem in this tutorial is 1D transient convectionreaction equation:
dc/dt + dc/dx = c**2
The exact solution is:
c(x,t) = sin(pi*(xt))**2 / (1  t*sin(pi*(xt))**2)
The equation is solved using the high resolution cellcentered finite volume upwind scheme with flux limiter described in the article. The boundary and initial conditions are obtained from the exact solution.
The consistent discretisation of the convection and the source terms cannot be done since the constant C1 in the integral of the source term:
S(x) = 1/u * Integral s(x)*dx = u**/3 + C1
is not known. Therefore, the numerical cell average is used:
Snum(x) = Integral s(x)*dx = s(i) * (x[i]x[i1]).
Numerical vs. manufactured solution plot (Nx=80):
The normalised global errors and the order of accuracy plots for the Koren flux limiter (grids 40, 60, 80, 120):
Files
Model report  tutorial_cv_8.xml 
Runtime model report  tutorial_cv_8rt.xml 
Source code  tutorial_cv_8.py 
Code verification using the Method of Exact Solutions (Solid Body Rotation problem).
Reference (section 4.4.6.1 Solid Body Rotation):
Solid body rotation illustrates the ability of a numerical scheme to transport initial data without distortion. Here, a 2D transient convection problem in a rectangular (0,1)x(0,1) domain is solved using the FE method:
dc/dt + div(u*c) = 0, in Omega = (0,1)x(0,1)
The initial conditions define a conical body which is rotated counterclockwise around the point (0.5, 0.5) using the velocity field u = (0.5  y, x  0.5). The cone is defined by the following equation:
r0 = 0.15
(x0, y0) = (0.5, 0.25)
c(x,y,0) = 1  (1/r0) * sqrt((xx0)**2 + (yy0)**2)
After t = 2pi the body should arrive at the starting position.
Homogeneous Dirichlet boundary conditions are prescribed at all four edges:
c(x,y,t) = 0.0
The mesh is a square (0,1)x(0,1):
The solution plot at t = 0 and t = 2pi (96x96 grid):
Animations for 32x32 and 96x96 grids:
It can be observed that the shape of the cone is preserved. Also, since the above equation is hyperbolic some oscillations in the solution out of the cone appear, which are more pronounced for coarser grids. This problem was resolved in the original example using the flux linearisation technique.
The normalised global errors and the order of accuracy plots (no. elements = [32x32, 64x64, 96x96, 128x128], t = 2pi):
Files
Model report  tutorial_cv_9.xml 
Runtime model report  tutorial_cv_9rt.xml 
Source code  tutorial_cv_9.py 
Code verification using the Method of Exact Solutions (Rotating Gaussian Hill problem).
Reference (section 4.4.6.3 ConvectionDiffusion):
Here, a 2D transient convectiondiffusion problem in a rectangular (1,1)x(1,1) domain is solved using the FE method:
dc/dt + div(u*c)  eps*nabla(c) = 0, in Omega = (1,1)x(1,1)
The exact solution is given by the following function:
(x0, y0) = (0.0, 0.5)
x_bar(t) = x0*cos(t)  y0*sin(t)
y_bar(t) = x0*sin(t) + y0*cos(t)
r2(x,y,t) = (xx_bar(t))**2 + (yy_bar(t))**2
c_exact(x,y,t) = 1.0 / (4*pi*eps*t) * exp(r2(x,y,t) / (4*eps*t))
The initial conditions define a Gaussian hill which is rotated counterclockwise around the point (0.0, 0.0) using the velocity field u = (y, x). Since at t = 0 the value of c_exact is the Dirac delta function it is better to start the simulation at t > 0. Therefore, the simulation is started and t = pi/2 by reinitialising variable c to:
c(x,y,pi/2) = c_exact(x,y,pi/2)
At t = 5/2 pi the peak smeared by the diffusion should arrive at the starting position.
Homogeneous Dirichlet boundary conditions are prescribed at all four edges:
c(x,y,t) = 0.0
The mesh is a rectangle (1,1)x(1,1):
The solution plots at t = pi/2 (the initial peak) and t = 5/2pi (96x96 grid):
Animations for 32x32 and 96x96 grids:
Again, some lowmagnitude oscillations in the solution appear, which are more pronounced for coarser grids. In the original example this problem was resolved using the flux linearisation technique.
The normalised global errors and the order of accuracy plots (no. elements = [32x32, 64x64, 96x96, 128x128], t = 5/2pi):
Files
Model report  tutorial_cv_10.xml 
Runtime model report  tutorial_cv_10rt.xml 
Source code  tutorial_cv_10.py 
Code verification using the Method of Exact Solutions.
The problem is identical to the problem in the tutorial_cv_6. The only difference is that the flow is reversed to test the high resolution scheme for the reversed flow mode.
Numerical vs. exact solution plots (Nx = [20, 40, 80]):
Files
Model report  tutorial_cv_11.xml 
Runtime model report  tutorial_cv_11rt.xml 
Source code  tutorial_cv_11.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):
Animation:
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:
Animation:
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:
Animation:
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:
Animation:
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 magnitude plot:
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 
In this example 2D transient Stokes flow driven by the differences in buoyancy caused by the temperature differences in the fluid is solved (deal.II step31).
The differences to the original problem are that the grid is not adaptive and no stabilisation method is used.
The problem can be described using the Stokes equations:
div(2 * eta * eps(u)) + nabla(p) = rho * beta * g * T in Omega
div(u) = 0 in Omega
dT/dt + div(u * T)  div(k * nabla(T)) = gamma in Omega
The mesh is a simple square (0,1)x(0,1):
The temperature and the velocity vectors plot:
Animation:
Files
Model report  tutorial_dealii_7.xml 
Runtime model report  tutorial_dealii_7rt.xml 
Source code  tutorial_dealii_7.py 
In this example a small parallelplate reactor with an active surface is modelled.
This problem and its solution in COMSOL Multiphysics software is described in the Application Gallery: Transport and Adsorption (id=5).
The transport in the bulk of the reactor is described by a convectiondiffusion equation:
dc/dt  D*nabla^2(c) + div(uc) = 0 in Omega
The material balance for the surface, including surface diffusion and the reaction rate is:
dc_s/dt  Ds*nabla^2(c_s) = k_ads * c * (Gamma_s  c_s) + k_des * c_s in Omega_s
For the bulk, the boundary condition at the active surface couples the rate of the reaction at the surface with the flux of the reacting species and the concentration of the adsorbed species and bulk species:
n⋅(D*nabla(c) + c*u) = k_ads*c*(Gamma_s  c_s) + k_des*c_s
The boundary conditions for the surface species are insulating conditions.
The problem is modelled using two coupled Finite Element systems: 2D for bulk flow and 1D for the active surface. The linear interpolation is used to determine the bulk flow and active surface concentrations at the interface.
The mesh is rectangular with the refined elements close to the left/right ends:
The cs plot at t = 2s:
The c plot at t = 2s:
Animation:
Files
Model report  tutorial_dealii_8.xml 
Runtime model report  tutorial_dealii_8rt.xml 
Source code  tutorial_dealii_8.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 Solubility Solubility
C kg Parac./kg EtOH 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:
Animation:
Files
Source code  tutorial_che_4.py 
Analytical solution  fl_analytical.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 
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 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 
This tutorial illustrates calculation of the sensitivity of the results with respect to the model parameters using forward sensitivity analysis method 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 integration of sensitivities per specified parameters 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  tutorial_sa_1.xml 
Runtime model report  tutorial_sa_1rt.xml 
Source code  tutorial_sa_1.py 
This tutorial illustrates the local derivativebased sensitivity analysis method available in DAE Tools.
The problem is adopted from the section 2.1 of the following article:
The model is very simple and describes a simple reversible chemical reaction A <> B, with reaction rates k1 and k_1 for the direct and inverse reactions, respectively. The reaction rates are uncertain and are described by continuous random variables with known probability density functions. The standard deviation is 0.3 for k1 and 1 for k_1. The standard deviation of the concentration of the species A is approximated using the following expression defined in the article:
stddev(Ca)**2 = stddev(k1)**2 * (dCa/dk1)**2 + stddev(k_1)**2 * (dCa/dk_1)**2
The following derivativebased measures are used in the article:
Derivatives dCa/dk1 and dCa/dk_1 calculated using the forward sensitivity method
Sigma normalised derivatives:
S(k1) = stddev(k1) / stddev(Ca) * dCa/dk1
S(k_1) = stddev(k_1)/ stddev(Ca) * dCa/dk_1
The plot of the concentrations, derivatives and sigma normalised derivatives:
Files
Model report  tutorial_sa_2.xml 
Runtime model report  tutorial_sa_2rt.xml 
Source code  tutorial_sa_2.py 
This tutorial illustrates the global variancebased sensitivity analysis methods available in the SALib python library.
The problem is adopted from the section 2.6 of the following article:
The model describes a thermal analysis of a batch reactor, with exothermic reaction A > B. The model equations are written in dimensionless form.
Three global sensitivity analysis methods available in SALib are applied:
Results from the sensitivity analysis:

Morris (N = 510)

Param mu mu* mu*_conf Sigma
B 0.367412 0.367412 0.114161 0.546276
gamma 0.040556 0.056616 0.021330 0.111878
psi 0.311563 0.311563 0.103515 0.504398
theta_a 0.326932 0.326932 0.102303 0.490423
theta_0 0.021208 0.023524 0.016015 0.074062


Sobol (N = 6144)

Param S1 S1_conf ST ST_conf
B 0.094110 0.078918 0.581930 0.154737
gamma 0.002416 0.012178 0.044354 0.027352
psi 0.171043 0.087782 0.524579 0.142115
theta_a 0.072535 0.042848 0.523394 0.165736
theta_0 0.002340 0.004848 0.008173 0.005956
Parameter pairs S2 S2_conf
B/gamma 0.180427 0.160979
B/psi 0.260689 0.171367
B/theta_a 0.143261 0.154060
B/theta_0 0.177129 0.156582
gamma/psi 0.000981 0.027443
gamma/theta_a 0.004956 0.036554
gamma/theta_0 0.009392 0.027390
psi/theta_a 0.166155 0.172727
psi/theta_0 0.016434 0.129177
theta_a/theta_0 0.109057 0.127287


FAST (N = 6150)

Param S1 ST
B 0.135984 0.559389
gamma 0.000429 0.029026
psi 0.171291 0.602461
theta_a 0.144617 0.534116
theta_0 0.000248 0.040741

The scatter plot for the Morris method:
The scatter plot for the Sobol method:
The scatter plot for the FAST method:
Files
Model report  tutorial_sa_3.xml 
Runtime model report  tutorial_sa_3rt.xml 
Source code  tutorial_sa_3.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 