6. User Guide¶

6.2.1. Parameters¶

Parameters are time invariant quantities that do not change during a simulation. Usually a good choice what should be a parameter is a physical constant, number of discretisation points in a domain etc.

Again, parameters are defined in two phases:

• Declaration in the __init__() function
• Initialisation (by setting its value) in the SetUpParametersAndDomains() function

self.myParam = daeParameter("myParam", units, parentModel, "description")

self.myParam = daeParameter("myParam", units, parentModel, "description")
self.myParam.DistributeOnDomain(myDomain)

# Or simply:
self.myParam = daeParameter("myParam", units, parentModel, "description", [myDomain])

# Ordinary parameters:
myParam.SetValue(value)

# Distributed parameters (one-dimensional):
for i in range(myDomain.NumberOfPoints):
    myParam.SetValue(i, value)

myParam.SetValues(values)

6.2.2. Variable types¶

Variable types describe variables and contain the information such as:

• Name: string
• Units: unit object
• LowerBound: float
• UpperBound: float
• InitialGuess: float
• AbsoluteTolerance: float
• ValueConstraint: enumeration daeeVariableValueConstraint

Declaration of variable types is commonly done outside of the model definition (in the module scope):

# Temperature type with units Kelvin, limits 100-1000K, the default value 273K and the absolute tolerance 1E-5
typeTemperature = daeVariableType("Temperature", K, 100, 1000, 273, 1E-5, constraint)

where the argument constraint specifies the value constraint and can be one of:

• eNoConstraint (default)
• eValueGTEQ: imposes >= 0 constraint
• eValueLTEQ: imposes <= 0 constraint
• eValueGT: imposes > 0 constraint
• eValueLT: imposes < 0 constraint

6.2.3. Distribution domains¶

Domains in DAE Tools are used to to create simple arrays of variables, parameters, equations and ports and their distributions in space. They can define uniform (the default) or non-uniform grids (user-specified) and are defined in two phases:

• Declaring a domain in the model
• Initialising it in the simulation

self.myDomain = daeDomain("myDomain", parentModel, units, "description")

6.2.3.2. Initialising domains¶

Domains are initialised in the SetUpParametersAndDomains() function. Simple arrays are defined using the function CreateArray():

# Array of N elements
myDomain.CreateArray(N)

while the domains distributed on a structured grid using the function CreateStructuredGrid():

# Uniform structured grid with N elements and bounds [lowerBound, upperBound]
myDomain.CreateStructuredGrid(N, lowerBound, upperBound)

where the lower and upper bounds are the floating point values or quantity objects. If the floats are used it is assumed that they contain values with the same units as in the domain definition. Using the quantities is advised (to avoid problems when the domain units change).

# Uniform structured grid with 10 elements and bounds [0,1] in centimeters:
myDomain.CreateStructuredGrid(10, 0.0 * cm, 1.0 * cm)

Nota bene

Structured grids with N elements contain N+1 points.

It is also possible to create an unstructured grid (for use in Finite Element models). However, their structure is an implementation detail (i.e. as in deal.II).

In certain situations it is desired to have a non-uniform distribution of the points within the given interval, defined by the lower and upper bounds. In these cases, a non-uniform structured grid can be specified using the attribute Points which contains the list of the points and that can be manipulated by the user:

# First create a structured grid domain
myDomain.CreateStructuredGrid(10, 0.0, 1.0)

# The original 11 points are: [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
# If the system is stiff at the beginning of the domain more points can be placed there
myDomain.Points = [0.0, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.60, 1.00]

The effect of uniform and non-uniform grids is given in Fig. 6.1 (a simple heat conduction problem from the Tutorial 3 has been served as a basis for comparison). Here, there are three cases:

• Black line: the analytic solution
• Blue line (10 intervals): uniform grid - a very rough prediction
• Red line (10 intervals): non-uniform grid - more points at the beginning of the domain

Fig. 6.1 Effect of uniform and non-uniform grids on numerical solution (zoomed to the first 5 points)

More precise results are obtained by using denser grid at the beginning of the interval.

# Notation:
#  - eq is a daeEquation object
#  - myDomain is daeDomain object
#  - myVar is daeVariable object
eq.Residual = myVar() - myDomain[5]

6.2.4. Variables¶

Variables define time varying quantities that change during a simulation and can be:

• Algebraic
• Differential
• Assigned (that is their value is assigned by fixing the number of degrees of freedom - DOF)

They are defined in two phases:

• Declaring a variable in the model
• Initialising it, if required (by assigning its value or setting an initial condition) in the simulation

# Ordinary variables:
self.myVar = daeVariable("myVar", variableType, parentModel, "description")

# Distributed variables:
self.myVar = daeVariable("myVar", variableType, parentModel, "description")
self.myVar.DistributeOnDomain(myDomain)
# or simply:
self.myVar = daeVariable("myVar", variableType, parentModel, "description", [myDomain])

6.2.5. Ports¶

Ports define connection points between model instances for exchange of continuous quantities. Like models, they can contain domains, parameters and variables.

New type of ports are defined by deriving a class from the base daePort:

class myPort(daePort):
    def __init__(self, name, parent = None, description = ""):
        daePort.__init__(self, name, type, parent, description)

        # Declaration/instantiation of domains, parameters and variables
        ...

The port structure (domains, parameters and variables) is declared in the __init__() function. The same rules apply as described in the section Developing models. Two ports are connected using the ConnectPorts() function. Ports are instantiated as inlet or outlet type in the __init__() function:

self.myPort = daePort("myPort", eInletPort, parentModel, "description")

6.2.6. Event ports¶

Event ports define connection points between model instances for exchange of discrete messages/events. Events can be triggered manually (using the SendEvent() function) or when a specified condition is satisfied. The main difference between event and ordinary ports is that the former allow a discrete communication between models while the latter allow a continuous exchange of information.

Messages are floating point values that can be used by a recipient. Upon a reception of an event, a set of actions are executed. The actions are specified in the ON_EVENT() function. The events received by an event port can be recorded by setting the boolean RecordEvents property to True and retrieved using the Events property.

Two event ports are connected using the ConnectEventPorts() function. A single outlet event port can be connected to unlimited number of inlet event ports. Event ports are instantiated in the __init__() function:

self.myEventPort = daeEventPort("myEventPort", eOutletPort, parentModel, "description")

6.2.7. Equations¶

Model equations are specified in an implicit (acausal) form. They can be continuous or discontinuous and distributed on one or more domains. Equations can be distributed on a whole domain, on a portion of it or even on a single point (i.e. equations for boundary conditions).

6.2.7.1. Declaring equations¶

Equations are declared in the DeclareEquations() function:

# Ordinary equation:
eq = model.CreateEquation("MyEquation", "description")

# Distributed equation (on the whole domain, including the boundaries):
eq = model.CreateEquation("MyEquation")
d = eq.DistributeOnDomain(myDomain, eClosedClosed)

Currently, the following options for distributing equations are available:

• Distributed on a closed (whole) domain: \(x \in [x_0, x_n]\)
• Distributed on a left open domain: \(x \in (x_0, x_n]\)
• Distributed on a right open domain: \(x \in [x_0, x_n)\)
• Distributed on a domain open on both sides: \(x \in (x_0, x_n)\)
• Distributed on the lower bound: \(x \in \{ x_0 \}\)
• Distributed on the upper bound: \(x \in \{ x_n \}\)
• Distributed on a given set of points within a domain: i.e. \(x \in \{ x_0, x_3, x_7, x_8 \}\)

where \(x_0\) represents the lower bound and \(x_n\) the upper bound of the domain.

An overview of available options is given in the table below. The examples are given for an equation distributed on two domains: x and y. Green squares represent portions of a domain included in the distributed equation, while white squares represent excluded portions.

x = eClosedClosed; y = eClosedClosed \(x \in [x_0, x_n], y \in [y_0, y_n]\)	x = eOpenOpen; y = eOpenOpen \(x \in ( x_0, x_n ), y \in ( y_0, y_n )\)
x = eClosedClosed; y = eOpenOpen \(x \in [x_0, x_n], y \in ( y_0, y_n )\)	x = eClosedClosed; y = eOpenClosed \(x \in [x_0, x_n], y \in ( y_0, y_n ]\)
x = eLowerBound; y = eClosedOpen \(x = x_0, y \in [ y_0, y_n )\)	x = eLowerBound; y = eClosedClosed \(x = x_0, y \in [y_0, y_n]\)
x = eUpperBound; y = eClosedClosed \(x = x_n, y \in [y_0, y_n]\)	x = eLowerBound; y = eUpperBound \(x = x_0, y = y_n\)

# Notation:
#  - V1, V3, V14 are ordinary variables
eq.Residual = dt(V14()) + V1() / (V14() + 2.5) + sin(3.14 * V3())

# Notation:
#  - V1 is an ordinary variable
#  - V3 and V14 are variables distributed on domains x and y
eq = model.CreateEquation("MyEquation")
dx = eq.DistributeOnDomain(x, eClosedClosed)
dy = eq.DistributeOnDomain(y, eOpenOpen)
eq.Residual = dt(V14(dx,dy)) + V1() / ( V14(dx,dy) + 2.5) + sin(3.14 * V3(dx,dy) )

# Notation:
#  - V1 is an ordinary variable
#  - V3 and V14 are variables distributed on domains x and y
for dx in range(0, x.NumberOfPoints): # x: [x0, xn]
    for dy in range(1, y.NumberOfPoints-1): # y: (y0, yn)
        eq = model.CreateEquation("MyEquation(%d,%d)" % (dx, dy) )
        eq.Residual = dt(V14(dx,dy)) + V1() / ( V14(dx,dy) + 2.5) + sin(3.14 * V3(dx,dy) )

# Notation:
#  - V1 and V2 are variables distributed on the x domain
eq = model.CreateEquation("MyEquation")
dx = eq.DistributeOnDomain(x, eClosedOpen)
eq.Residual = V1(dx) - ( V2(dx) + V2(dx+1) )

# Notation:
#  - Var is an ordinary variable
#  - result is an ordinary variable
eq = self.CreateEquation("...")
eq.Residual = result() - numpy.exp( Var() )

# The above is identical to:
eq.Residual = result() - Exp( Var() )

# Notation:
#  - x is a continuous domain
#  - V1 is a variable distributed on the x domain
#  - V2 is a variable distributed on the x domain
#  - sum is an ordinary variable
#  - ndarr_V1 is one dimensional numpy array with dtype=object
#  - ndarr_V2 is one dimensional numpy array with dtype=object
#  - adarr_V1 is adouble_array object
#  - Nx is the number of points in the domain x

# 1. Create empty numpy arrays as a container for daetools adouble objects
ndarr_V1 = numpy.empty(Nx, dtype=object)
ndarr_V2 = numpy.empty(Nx, dtype=object)

# 2. Fill the created numpy arrays with adouble objects
ndarr_V1[:] = [V1(x) for x in range(Nx)]
ndarr_V2[:] = [V2(x) for x in range(Nx)]

# Now, ndarr_V1 and ndarr_V2 represent arrays of Nx adouble objects each:
#  ndarr_V1 := [V1(0), V1(1), V1(2), ..., V1(Nx-1)]
#  ndarr_V2 := [V2(0), V2(1), V2(2), ..., V2(Nx-1)]

# 3. Create an equation using the common numpy/scipy functions/operators
eq = self.CreateEquation("sum")
eq.Residual = sum() - numpy.sum(ndarr_V1 + 2*ndarr_V2**2)

# If adouble_array is needed after operations on a numpy array, the following two functions can be used:
#   a) static function adouble_array.FromList(python_list)
#   b) static function adouble_array.FromNumpyArray(numpy_array)
# Both return an adouble_array object.
adarr_V1 = adouble_array.FromNumpyArray(ndarr_V1)
print(adarr_V1)

6.2.7.5. Details on autodifferentiation support¶

To calculate a residual and its gradients the operator overloading technique for automatic differentiation is applied (ADOL-C). Model equations are stored in a tree-like data structure called Evaluation Tree. Evaluation Trees consist of nodes representing mathematical operators and functions and their operands. In DAE Tools the basic mathematical operators and functions are re-defined to operate on a modified ADOL-C class adouble extended with the simulator-specific information. In addition, all operators/functions accept the adouble_array objects to support operations on arrays. Once built, Evaluation Trees can be used to calculate equation residuals or derivatives, check units consistency and export equations into the Latex or MathML format. For instance, the equation:

\[\frac{x[0]-x[1]} {x[2]-x[3]} + 1.2 \cdot sin(x[0]) = 0\]

is specified in DAE Tools as:

eq = model.CreateEquation("myEquation")
eq.Residual = (x(0)-x(1))/(x(2)-x(3)) + 1.2 * sin(x(0))

and represented as the evaluation tree in Fig. 6.2.

Fig. 6.2 Evaluation Tree

As it has been described in the previous sections, domains, parameters, and variables contain functions that return adouble/adouble_array objects used to construct the evaluation trees. These functions include functions to get a value of a parameter/variable (__call__()), a time or a partial derivative of a variable (functions dt(), d(), or d2()) or functions to obtain an array of values, time or partial derivatives (array(), dt_array(), d_array(), and d2_array()).

# Notation:
#  - T is a variable distributed on x and y domains
#  - rho, k, and cp are parameters
eq = model.CreateEquation("MyEquation")
dx = eq.DistributeOnDomain(x, eClosedClosed)
dy = eq.DistributeOnDomain(y, eOpenOpen)
eq.Residual = rho() * cp() * dt(T(dx,dy)) - k() * ( d2(T(dx,dy), x) + d2(T(dx,dy), y) )

# The Dirichlet boundary conditions at the "bottom edge" :
bceq = model.CreateEquation("Bottom_BC")
dx = bceq.DistributeOnDomain(x, eClosedClosed)
dy = bceq.DistributeOnDomain(y, eLowerBound)
bceq.Residual = T(dx,dy) - Constant(500 * K)  # Constant temperature (500 K)

# The Neumann boundary conditions at the "top edge":
bceq = model.CreateEquation("Top_BC")
dx = bceq.DistributeOnDomain(x, eClosedClosed)
dy = bceq.DistributeOnDomain(y, eUpperBound)
bceq.Residual = - k() * d(T(dx,dy), y) - Constant(1E6 * W/m**2)  # Constant flux (1E6 W/m2)

def DeclareEquations(self):
    daeModel.DeclareEquations(self)

    # Auxiliary functions and objects
    rho     = self.rho()
    Q       = lambda i:      self.Q(i)
    cp      = lambda x,y:    self.cp(x,y)
    k       = lambda x,y:    self.k(x,y)
    T       = lambda x,y:    self.T(x,y)
    dT_dt   = lambda x,y: dt(self.T(x,y))
    dT_dx   = lambda x,y:  d(self.T(x,y), self.x, eCFDM)
    dT_dy   = lambda x,y:  d(self.T(x,y), self.y, eCFDM)
    d2T_dx2 = lambda x,y: d2(self.T(x,y), self.x, eCFDM)
    d2T_dy2 = lambda x,y: d2(self.T(x,y), self.y, eCFDM)

    eq = self.CreateEquation("HeatBalance", "Heat balance equation valid on the open x and y domains")
    x = eq.DistributeOnDomain(self.x, eOpenOpen)
    y = eq.DistributeOnDomain(self.y, eOpenOpen)
    eq.Residual = rho * cp(x,y) * dT_dt(x,y) - k(x,y) * (d2T_dx2(x,y) + d2T_dy2(x,y))

    eq = self.CreateEquation("BC_bottom", "Neumann boundary conditions at the bottom edge (constant flux)")
    x = eq.DistributeOnDomain(self.x, eOpenOpen)
    y = eq.DistributeOnDomain(self.y, eLowerBound)
    eq.Residual = -k(x,y) * dT_dy(x,y) - Q(0)

    eq = self.CreateEquation("BC_top", "Neumann boundary conditions at the top edge (constant flux)")
    x = eq.DistributeOnDomain(self.x, eOpenOpen)
    y = eq.DistributeOnDomain(self.y, eUpperBound)
    eq.Residual = -k(x,y) * dT_dy(x,y) - Q(1)

    eq = self.CreateEquation("BC_left", "Neumann boundary conditions at the left edge (insulated)")
    x = eq.DistributeOnDomain(self.x, eLowerBound)
    y = eq.DistributeOnDomain(self.y, eClosedClosed)
    eq.Residual = dT_dx(x,y)

    eq = self.CreateEquation("BC_right", " Neumann boundary conditions at the right edge (insulated)")
    x = eq.DistributeOnDomain(self.x, eUpperBound)
    y = eq.DistributeOnDomain(self.y, eClosedClosed)
    eq.Residual = dT_dx(x,y)

...

eq.Residual = rho * cp(x,y) * dT_dt(x,y) - k(x,y) * (d2T_dx2(x,y) + d2T_dy2(x,y))

...

eq.Residual = self.rho() * self.cp() * dt(self.T(x,y)) - \
              self.k() * (d2(self.T(x,y), self.x, eCFDM) + d2(self.T(x,y), self.y, eCFDM))

6.2.8. State Transition Networks¶

Discontinuous equations are equations that take different forms subject to certain conditions. In DAE Tools they are modelled using the concept of State Transition Networks (STN). State Transition Networks can be:

1. Reversible (internally described by daeIF class)
2. Irreversible (internally described by daeSTN class)

Every STN consists of two or more states (daeState class). Each state is described by a set of equations and conditions for transitions between states. Irreversible STNs can optionally contain a set of actions to be performed when a specified logical condition is satisfied. At any given moment, there is one active state.

For example, a flow of fluid through a pipe is described by three different regimes:

• Laminar: if Reynolds number is less than 2,100
• Transient: if Reynolds number is greater than 2,100 and less than 10,000
• Turbulent: if Reynolds number is greater than 10,000

From any of these three states the system can go to any other state. This type of discontinuities is called a reversible discontinuity and constructed using the IF(), ELSE_IF(), ELSE() and END_IF() functions:

IF(Re() <= 2100)                    # Laminar flow
#... (equations go here)

ELSE_IF(Re() > 2100 & Re() < 10000) # Transient flow
#... (equations go here)

ELSE()                              # Turbulent flow
#... (equations go here)

END_IF()

The comparison operators operate on adouble objects and floating point values. Units consistency is strictly checked and expressions including plain numbers values are allowed only if a variable or parameter is dimensionless. The following expressions are valid:

# Notation:
#  - T is a variable with units: K
#  - m is a variable with units: kg
#  - p is a dimensionless parameter

# T < 0.5 K
T() < Constant(0.5 * K)

# (T >= 300 K) or (m < 1 kg)
(T() >= Constant(300 * K)) | (m < Constant(0.5 * kg))

# p <= 25.3 (use of the Constant function not necessary)
p() <= 25.3

Reversible discontinuities can be symmetrical and non-symmetrical. The above example is symmetrical. Non-symmetrical STNs are used to describe the concept of hysteresis.

Discontinuities can also be irreversible. For instance, a CPU and its power dissipation can operate in three operating modes:

• Normal mode
  • switch to Power saving mode if CPU load is below 5%
  • switch to Fried mode if the temperature is above 110 degrees
• Power saving mode
  • switch to Normal mode if CPU load is above 5%
  • switch to Fried mode if the temperature is above 110 degrees
• Fried mode
  • Damn, no escape from here... go to the nearest shop and buy a new one! Or, donate some money to DAE Tools project :-)

From the Normal mode the system can either go to the Power saving mode or to the Fried mode. The same stands for the Power saving mode: the system can either go to the Normal mode or to the Fried mode. However, once the temperature exceeds 110 degrees the CPU dies (let’s say it is heavily overclocked) and the system remains in this state permanently. This type of discontinuities is called an irreversible discontinuity and can be described using STN(), STATE(), END_STN() functions:

STN("CPU")

STATE("Normal")
#... (equations go here)
ON_CONDITION( CPULoad() < 0.05,       switchToStates = [ ("CPU", "PowerSaving") ] )
ON_CONDITION( T() > Constant(110*K),  switchToStates = [ ("CPU", "Fried") ] )

STATE("PowerSaving")
#... (equations go here)
ON_CONDITION( CPULoad() >= 0.05,      switchToStates = [ ("CPU", "Normal") ] )
ON_CONDITION( T() > Constant(110*K),  switchToStates = [ ("CPU", "Fried") ] )

STATE("Fried")
#... (equations go here)

END_STN()

The function ON_CONDITION() is used to define actions to be performed when the specified condition is satisfied. In addition, the function ON_EVENT() can be used to define actions to be performed when an event is triggered on a specified event port. Details on how to use ON_CONDITION() and ON_EVENT() functions can be found in the OnCondition actions and OnEvent actions sections, respectively.

More information about state transition networks can be found in daeSTN, daeIF and in Tutorials.

6.2.9. OnCondition actions¶

The function ON_CONDITION() is used to define actions to be performed when a specified condition is satisfied. The available actions include:

• Changing the active state in specified State Transition Networks (argument switchToStates)
• Re-assigning or re-initialising specified variables (argument setVariableValues)
• Triggering an event on the specified event ports (argument triggerEvents)
• Executing user-defined actions (argument userDefinedActions)

For instance, to execute some actions when the temperature becomes greater than 340 K the following can be used:

def DeclareEquations(self):
    ...

    self.ON_CONDITION( T() > Constant(340*K), switchToStates     = [ ('STN', 'StateName'), ... ],
                                              setVariableValues  = [ (variable, newValue), ... ],
                                              triggerEvents      = [ (eventPort, eventMessage), ... ],
                                              userDefinedActions = [ userDefinedAction, ... ] )

where the first argument of the ON_CONDITION() function is a condition specifying when the actions will be executed and:

• switchToStates is a list of tuples (string ‘STN Name‘, string ‘State name to activate‘)
• setVariableValues is a list of tuples (daeVariable object, adouble object)
• triggerEvents is a list of tuples (daeEventPort object, adouble object)
• userDefinedActions is a list of user defined objects derived from the base daeAction class

More details on how to use ON_CONDITION() function can be found in Tutorial 13.

6.2.10. OnEvent actions¶

The function ON_EVENT() is used to define actions to be performed when an event is triggered on the specified event port. The same type of actions are available as in the ON_CONDITION() function.

For instance, to execute some actions when an event is triggered on an event port the following can be used:

def DeclareEquations(self):
    ...

    self.ON_EVENT( eventPort, switchToStates     = [ ('STN', 'StateName'), ... ],
                              setVariableValues  = [ (variable, newValue), ... ],
                              triggerEvents      = [ (eventPort, eventMessage), ... ],
                              userDefinedActions = [ userDefinedAction, ... ] )

where the first argument of the ON_EVENT() function is the daeEventPort object to be monitored for events, while the rest of the arguments is the same as in the ON_CONDITION() function.

More details on how to use ON_EVENT() function can be found in Tutorial 13.

6.2.11. User-defined actions¶

User-defined actions are executed in a response to a specified condition in OnCondition handlers or in a response to an event triggered in OnEvent handlers.

They are created by deriving a class from the daeAction base and implementing the Execute() function. The Execute() function takes no arguments. If some information from the model is required they should be specified in the constructor.

User-defined actions do not return a value and should not change the values of variables (other types of actions must be used for that purpose), but perform some user-defined operations. The source code for a simple action that prints a message with the data sent to a specified event port is given below:

# User-defined action executed when an event is triggered on a specified event port.
class simpleUserAction(daeAction):
    def __init__(self, eventPort):
        daeAction.__init__(self)

        # Store the daeEventPort object for later use.
        self.eventPort = eventPort

    def Execute(self):
        # The floating point value of the data sent when the event is triggered
        # can be retrieved using the daeEventPort.EventData property.
        msg = 'simpleUserAction executed; input data = %f' % self.eventPort.EventData

        print('********************************************************')
        print(msg)
        print('********************************************************')

def DeclareEquations(self):
    ...

    # User-defined action objects should be stored in the model, otherwise
    # they will get destroyed when they go out of scope.
    self.action = simpleUserAction(self.eventPort)

    # The actions executed when the event on the inlet 'eventPort' event port is received.
    # daeEventPort defines the operator() which returns adouble object that can be used
    # at the moment when the action is executed to get the value of the event data.
    self.ON_EVENT(self.eventPort, userDefinedActions = [self.action])

More details on how to use user-defined actions can be found in Tutorial 13.

6.2.12. External functions¶

The concept of external functions in DAE Tools is used to handle and evaluate user-defined functions or functions from external libraries. External functions can return scalar (daeScalarExternalFunction) or vector (daeVectorExternalFunction) values.

External functions are created by deriving a class from the daeScalarExternalFunction base, specifying its arguments in the constructor and implementing the Calculate() function. The source code for a simple \(F(x) = x ^ 2\) external function is given below:

class F(daeScalarExternalFunction):
    def __init__(self, Name, parentModel, units, x):
        # Instantiate the scalar external function by specifying
        # the dictionary with arguments {'name' : adouble-object}
        arguments = {}
        arguments["x"]  = x

        daeScalarExternalFunction.__init__(self, Name, parentModel, units, arguments)

    def Calculate(self, values):
        # Calculate function is used to calculate a value and a derivative of the external
        # function per given argument (if requested). Here, a simple function is given by:
        #    F(x) = x**2

        # Procedure:
        # 1. Get the arguments from the dictionary 'values': {'arg-name' : adouble-object}.
        #    Every adouble object has two properties: Value and Derivative that can be
        #    used to evaluate function or its partial derivatives per arguments
        #    (partial derivatives are used to fill in a Jacobian matrix necessary to solve
        #    a system of non-linear equations using the Newton method).
        x = values["x"]

        # 2. Always calculate the value of a function (the derivative is zero by default).
        res = adouble(x.Value ** 2)

        # 3. If a function derivative per one of its arguments is requested,
        #    the derivative part of that argument will be non-zero.
        #    In that case, investigate which derivative is requested and calculate it
        #    using the chain rule: f'(x) = x' * df(x)/dx
        if x.Derivative != 0:
            # A derivative per 'x' was requested; its value is: x' * 2x
            res.Derivative = x.Derivative * (2 * x.Value)

        # 4. Return the result as a adouble object (contains both a value and a derivative)
        return res

def DeclareEquations(self):
    ...

    # Create external function (it has to be created in DeclareEquations!),
    # specify its units (here for simplicity dimensionless) and
    # arguments (here only a single argument: x)
    # External function objects should be stored in the model, otherwise
    # they will get destroyed when they go out of scope.
    self.F = F("F", self, unit(), self.x())

    # External function can now be used in daetools equations.
    # Its value can be obtained using the operator() (python special function __call__)
    eq = self.CreateEquation("...", "...")
    eq.Residual = ... self.F() ...

A more complex example is given in the Tutorial 14. There, the external function concept is used to interpolate a set of values using the scipy.interpolate.interp1d object.

class extfn_interp1d(daeScalarExternalFunction):
    def __init__(self, Name, parentModel, units, times, values, Time):
        arguments = {}
        arguments["t"] = Time

        # Instantiate interp1d object and initialise interpolation using supplied (time,y) values.
        self.interp = scipy.interpolate.interp1d(times, values)

        # During the solver iterations, the function is called very often with the same arguments.
        # Therefore, cache the last interpolated value to speed up a simulation.
        self.cache = None

        daeScalarExternalFunction.__init__(self, Name, parentModel, units, arguments)

    def Calculate(self, values):
        # Get the argument from the dictionary of arguments' values.
        time = values["t"].Value

        # Here we do not need to return a derivative for it is not a function of variables.

        # First check if an interpolated value was already calculated during the previous call.
        # If it was, return the cached value (the derivative part is always equal to zero in this case).
        if self.cache:
            if self.cache[0] == time:
                return adouble(self.cache[1])

        # The time received is not in the cache and has to be interpolated.
        # Convert the result to float datatype since daetools can't accept
        # numpy.float64 types as arguments at the moment.
        interp_value = float(self.interp(time))
        res = adouble(interp_value, 0)

        # Save it in the cache for later use.
        self.cache = (time, res.Value)

        return res

The extfn_interp1d class is used here to approximate some function f:

using its t ad y values:

def DeclareEquations(self):
    ...

    # Create scipy.interp1d interpolation external function.
    # Create 'times' and 'values' arrays to be used for interpolation:
    times  = numpy.arange(0.0, 1000.0)
    values = 2*times
    # The external function accepts only a single argument: the current time in the simulation
    # that can be obtained using the Time() daetools function.
    # The external function units are seconds.
    self.interp1d = extfn_interp1d("interp1d", self, s, times, values, Time())

Alternatively, DAE Tools can utilise functions defined in shared libraries via daeCTypesExternalFunction class. As an argument it accepts a function pointer from libraries loaded using Python ctypes. Sample usage can be found in the Tutorial 14:

def DeclareEquations(self):
    ...

    # Load the library:
    self.ext_lib = ctypes.CDLL("libheat_function.so")

    # Function arguments:
    arguments = {}
    arguments['m']     = self.m()
    arguments['cp']    = self.cp()
    arguments['dT/dt'] = dt(self.T())

    # Function pointer ('calculate' function is used):
    function_ptr = self.ext_lib.calculate

    self.exfnHeat2 = daeCTypesExternalFunction("heat_function", self, W, function_ptr, arguments)

The calculate function is defined in the heat_function c shared library:

#include <string.h>

typedef struct
{
    double Value;
    double Derivative;
}
adouble_c;

#if defined(_WIN32) || defined(WIN32) || defined(WIN64) || defined(_WIN64)
#define DLLEXPORT  extern "C" __declspec(dllexport)
#else
#define DLLEXPORT
#endif

DLLEXPORT adouble_c calculate(const adouble_c values[], const char* names[], int no_arguments);

adouble_c calculate(const adouble_c values[], const char* names[], int no_arguments)
{
    adouble_c result;
    memset(&result, 0, sizeof(adouble_c));

    /* Get the arguments' values. */
    adouble_c m, cp, dT_dt;
    for(int i = 0; i < no_arguments; i++)
    {
        if(strcmp(names[i], "m") == 0)
            m = values[i];
        else if(strcmp(names[i], "cp") == 0)
            cp = values[i];
        else if(strcmp(names[i], "dT/dt") == 0)
            dT_dt = values[i];
    }

    /* Calculate the value. */
    result.Value = m.Value * cp.Value * dT_dt.Value;

    /* Calculate the derivative. */
    if(m.Derivative != 0) /* A derivative per 'm' was requested */
        result.Derivative = m.Derivative * (cp.Value * dT_dt.Value);
    else if(cp.Derivative != 0) /* A derivative per 'cp' was requested */
        result.Derivative = cp.Derivative * (m.Value * dT_dt.Value);
    else if(dT_dt.Derivative != 0) /* A derivative per 'dT_dt' was requested */
        result.Derivative = dT_dt.Derivative * (m.Value * cp.Value);

    return result;
}

The library can be compiled using the following commands:

# GNU/Linux gcc:
gcc -fPIC -shared -o libheat_function.so tutorial4_heat_function.c

# macOS gcc:
gcc -fPIC -dynamiclib -o libheat_function.dylib tutorial14_heat_function.c

# Windows vc++:
cl /LD tutorial14_heat_function.c /link /dll /out:heat_function.dll

6.3.1. The Finite Difference Method¶

DAE Tools support numerical simulation of partial differential equations on structured grids using the Finite Difference Method. Three different methods are provided:

• Backward Finite Difference method (eBFDM)
• Forward Finite Difference method (eFFDM)
• Center Finite Difference method (eCFDM)

The partial derivatives of the first and second order can be specified using the functions d() and d2().

As an illustration, the 1D convection-diffusion-reaction equation:

can be specified using the Center Finite Difference Method:

class modTutorial(daeModel):
    ...
    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        # Notation:
        #  - c is a state variable
        #  - x is a domain object
        #  - u is velocity

        # Declare some auxiliary functions to make equations more readable
        c       = lambda i: self.c(i)
        dc_dt   = lambda x: dt(c(x))
        dc_dx   = lambda x: d (c(x), self.x, eCFDM)
        d2c_dx2 = lambda x: d2(c(x), self.x, eCFDM)
        s       = lambda i: c(i)**2

        # Declare the Convection-Diffusion-Reaction equation distributed on (0, L]:
        eq = self.CreateEquation("c")
        eq.DistributeOnDomain(self.x, eOpenClosed)
        eq.Residual = dc_dt(x) + u * dc_dx(x) - D * d2c_dx2(x) - s(x)

        # Boundary conditions at x = 0:
        eq = self.CreateEquation("c(0)")
        eq.Residual = c(0) - 0.0

6.3.2. The Finite Volume Method¶

DAE Tools support numerical simulation of partial differential equations on 1D structured grids using the Finite Volume Method (high-resolution upwind scheme with flux limiter).

Consider the 1D convection-diffusion-reaction equation:

A cell-centered finite-volume discretisation yields the semi-discrete equation [1] [2]:

where the half-integer indices refer to cell faces \(\delta\Omega_{i-{1 \over 2}}\) and \(\delta\Omega_{i+{1 \over 2}}\) between cell centers \(\Omega_{i-1}\) and \(\Omega_i\) as presented in the figure below:

The accuracy of the above finite volume discretisation is determined by the way in which the cell-face fluxes are computed. Applying the high-resolution upwind scheme with flux limiter [1] [2] for the cell-face state \(c_{i+{1 \over 2}}\) results in the following equation:

where \(\phi\) is the flux limiter function and \(r_{i + {1 \over 2}}\) the upwind ratio of consecutive solution gradients:

There is a large number of flux limiters [3] implemented in DAE Tools:

• CHARM [not 2nd order TVD] (Zhou, 1995):

  \(\phi(r)= \begin{cases} \frac{r\left(3r+1\right)}{\left(r+1\right)^{2}} & r>0, \lim_{r \rightarrow \infty} \phi(r)=3 \\ 0 & otherwise \end{cases}\)

• HCUS (not 2nd order TVD) (Waterson and Deconinck, 1995):

  \(\phi(r) = \frac{ 1.5 \left(r+\left| r \right| \right)}{ \left(r+2 \right)} , \lim_{r \rightarrow \infty}\phi_{hc}(r) = 3\)

• HQUICK (not 2nd order TVD) (Waterson and Deconinck, 1995):

  \(\phi(r) = \frac{2 \left(r + \left|r \right| \right)}{ \left(r+3 \right)}, \lim_{r \rightarrow \infty}\phi_{hq}(r) = 4\)

• Koren (Koren, 1993):

  \(\phi(r) = \max \left[ 0, \min \left(2 r, \left(2 + r \right)/3, 2 \right) \right], \lim_{r \rightarrow \infty}\phi(r) = 2\)

• minmod - symmetric (Philip and Roe, 1986):

  \(\phi (r) = \max \left[ 0 , \min \left( 1 , r \right) \right], \lim_{r \rightarrow \infty}\phi(r) = 1\)

• monotonized central (MC) – symmetric (van Leer, 1977):

  \(\phi (r) = \max \left[ 0 , \min \left( 2 r, 0.5 (1+r), 2 \right) \right] , \lim_{r \rightarrow \infty}\phi(r) = 2\)

• Osher (Chatkravathy and Osher, 1983):

  \(\phi (r) = \max \left[ 0 , \min \left( r, \beta \right) \right], \left(1 \leq \beta \leq 2 \right), \lim_{r \rightarrow \infty}\phi (r) = \beta\)

• ospre - symmetric (Waterson and Deconinck, 1995):

  \(\phi (r) = \frac{1.5 \left(r^2 + r \right) }{\left(r^2 + r +1 \right)} , \lim_{r \rightarrow \infty}\phi (r) = 1.5\)

• smart (not 2nd order TVD) (Gaskell and Lau, 1988):

  \(\phi(r) = \max \left[ 0, \min \left(2 r, \left(0.25 + 0.75 r \right), 4 \right) \right], \lim_{r \rightarrow \infty}\phi(r) = 4\)

• superbee – symmetric (Roe, 1986):

  \(\phi (r) = \max \left[ 0, \min \left( 2 r , 1 \right), \min \left( r, 2 \right) \right] , \lim_{r \rightarrow \infty}\phi (r) = 2\)

• Sweby – symmetric (Sweby, 1984):

  \(\phi (r) = \max \left[ 0 , \min \left( \beta r, 1 \right), \min \left( r, \beta \right) \right], \left(1 \leq \beta \leq 2 \right), \lim_{r \rightarrow \infty}\phi (r) = \beta\)

• UMIST (Lien and Leschziner, 1994):

  \(\phi(r) = \max \left[ 0, \min \left(2 r, \left(0.25 + 0.75 r \right), \left(0.75 + 0.25 r \right), 2 \right) \right] , \lim_{r \rightarrow \infty}\phi(r) = 2\)

• van Albada 1 - symmetric (van Albada, et al., 1982):

  \(\phi (r) = \frac{r^2 + r}{r^2 + 1 } , \lim_{r \rightarrow \infty}\phi (r) = 1\)

• van Albada 2 : alternative form (not 2nd order TVD; Kermani, 2003)

  \(\phi (r) = \frac{2 r}{r^2 + 1}, \lim_{r \rightarrow \infty}\phi(r) = 0\)

• van Leer - symmetric (van Leer, 1974):

  \(\phi (r) = \frac{r + \left| r \right| }{1 + \left| r \right| } , \lim_{r \rightarrow \infty}\phi (r) = 2\)

For the diffusive flux, the gradient \(\left( \partial c \over \partial x \right)_{i + {1 \over 2}}\) is evaluated using the standard second-order accurate central difference formula:

except at the inflow and outflow boundaries where:

The above convection-diffusion-reaction equation can be specified using the daeHRUpwindSchemeEquation class with the following functions:

• Accumulation term in the cell-centered finite-volume discretisation: dc_dt():

  \(dc\_dt(i) = \int_{\Omega_i} {\partial c_i \over \partial t} dx\)

• Convection term in the cell-centered finite-volume discretisation: dc_dx() (may contain the \(\mathbf{S} = {1 \over u} \int_{\Omega_i} s(x) dx\) integral for the consistent discretisation of the convection and the source terms):

  \(dc\_dx(i) = c_{i + {1 \over 2}} - c_{i - {1 \over 2}}\)

  or (if the source integral \(\mathbf{S}\) has been specified):

  \(dc\_dx(i) = \left( c_{i + {1 \over 2}} - \mathbf{S}_{i + {1 \over 2}} \right) - \left( c_{i - {1 \over 2}} - \mathbf{S}_{i - {1 \over 2}} \right)\)

• Diffusion term in the cell-centered finite-volume discretisation: d2c_dx2():

  \(d2c\_dx2(i) = \left( \partial c \over \partial x \right)_{i + {1 \over 2}} - \left( \partial c \over \partial x \right)_{i - {1 \over 2}}\)

• Source term in the cell-centered finite-volume discretisation: source():

  \(source(i) = \int_{\Omega_i} s_i dx\)

as given in the example below:

class modTutorial(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)

    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        xp = self.x.Points
        Nx = self.x.NumberOfPoints

        c = lambda i: self.c(i)

        # 1. Declare the HR upwind scheme object:
        #     - c is a state variable
        #     - x is a domain object
        #     - Phi_Koren is a flux limiter function
        hr = daeHRUpwindSchemeEquation(self.c, self.x, daeHRUpwindSchemeEquation.Phi_Koren, 1e-10)

        # 2. Define the source term function
        def s(i):
            return self.c(i)**2

        # 3. Declare the Convection-Diffusion-Reaction equation distributed on (0, L]:
        for i in range(1, Nx):
            eq = self.CreateEquation("c(%d)" % i)
            eq.Residual = hr.dc_dt(i) + u * hr.dc_dx(i) - D * hr.d2c_dx2(i) - hr.source(s,i)

        # Boundary conditions at x=0:
        eq = self.CreateEquation("c(0)")
        eq.Residual = c(0) - ...

It is desired that the discretisation of the source term is consistent with that of the advection operator. For this purpose, if the source term integral: \(S(x) = {1 \over u} \int_{\Omega_i} s(x) dx\) can be calculated analytically, the convection term can be rewritten as:

and the semi-discrete equation becomes:

The example above now becomes:

class modTutorial(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)

    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        xp = self.x.Points
        Nx = self.x.NumberOfPoints

        c = lambda i: self.c(i)

        # 1. Declare the HR upwind scheme object:
        #     - c is a state variable
        #     - x is a domain object
        #     - u is velocity
        #     - Phi_Koren is a flux limiter function
        hr = daeHRUpwindSchemeEquation(self.c, self.x, daeHRUpwindSchemeEquation.Phi_Koren, 1e-10)

        # 2. Consistent discretisation of convection and source terms:
        #    Calculate an analytical integral of the source term S = 1/u * Integral(s(x)*dx)
        def S(i):
            C1 = 0.0
            return c(i)**2 / 3 + C1

        # 3. Declare the Convection-Diffusion-Reaction equation distributed on (0, L]
        #    (with the S integral in the convection term):
        for i in range(1, Nx):
            eq = self.CreateEquation("c(%d)" % i)
            eq.Residual = hr.dc_dt(i) + u * hr.dc_dx(i, S) - D * hr.d2c_dx2(i)

        # Boundary Conditions:
        eq = self.CreateEquation("c(0)")
        eq.Residual = c(0) - ...

Footnotes

6.3.3. The Finite Element Method¶

DAE Tools support numerical simulation of partial differential equations on unstructured grids using the Finite Element method. The main idea is to utilise available state-of-the-art Finite Element libraries (i.e. deal.II) for low-level tasks such as mesh loading, management of finite element spaces, degrees of freedom, assembly of the system stiffness and mass matrices and the system load vector, setting the boundary conditions etc.. After the assembly phase in an external library the matrices are used to generate a set of equations in the following form:

where \(x_j\) is a state variable, \(M_{ij}\) and \(A_{ij}\) are mass and stiffness matrices and \(F_i\) is the load vector. This system is in a general case a DAE system, although it can also be a linear or a non-linear (if the mass matrix is zero). The generated set of equations are solved together with the rest of equations in the model.

The unique feature of this approach is a capability to use DAE Tools objects (parameters and variables) as a native data type in deal.II functions to specify boundary conditions, time varying coefficients and source terms. This way, the non-linear finite element systems are automatically supported and the equations resulting from the finite element discretisation are fully integrated with the rest of the model equations. Moreover, multiple FE systems can be created and coupled together.

The information required (from the modeller’s perspective):

• Mesh file with the specified boundary indicators (integers)
• Variables (degrees of freedom in deal.II) and their Finite Element spaces
• Quadrature formulas for elements and their faces
• The weak form of the problem which contains expressions for the cells’ and boundary faces’ contributions to the system mass and stiffness matrix and the load vector.

The weak form expressions are specified using the DAE Tools API that wraps deal.II concepts used to assembly the matrices/vectors. The weak forms in daetools contain expressions as they would appear in typical nested for loops. In deal.II a typical cell assembly loop (in C++) would look like (i.e. a very simple example given in step-7):

FEValues<dim>             fe_values(...);
std::vector<unsigned int> local_dof_indices;

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                               endc = dof_handler.end();
for(; cell != endc; ++cell)
{
    fe_values.reinit(cell);
    cell->get_dof_indices(local_dof_indices);

    for(unsigned int q_point = 0; q_point < n_q_points; ++q_point)
    {
        for(unsigned int i = 0; i < dofs_per_cell; ++i)
        {
            for(unsigned int j = 0; j < dofs_per_cell; ++j)
            {
                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
                                      fe_values.shape_grad(j,q_point)
                                      +
                                      fe_values.shape_value(i,q_point) *
                                      fe_values.shape_value(j,q_point)) *
                                    fe_values.JxW(q_point));
            }

            cell_rhs(i) += (fe_values.shape_value(i,q_point) *
                            rhs_values [q_point] * fe_values.JxW(q_point));
        }
    }
}

An equivalent problem in DAE Tools creates an execution context used in a generic loop to evaluate the cell/face contributions:

FEValues<dim>             fe_values(...);
std::vector<unsigned int> local_dof_indices;

// Create evaluation context objects where the expressions will be evaluated:
feCellContext<dim> cellContext(fe_values, local_dof_indices, ...);

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                               endc = dof_handler.end();
for(; cell != endc; ++cell)
{
    // Update the evaluation context with the current cell.
    // It will call fe_values.reinit(cell), cell->get_dof_indices(local_dof_indices) etc.
    update(cellContext, cell);

    for(unsigned int q_point = 0; q_point < n_q_points; ++q_point)
    {
        // update cellContext with the current quadrature point index

        for(unsigned int i = 0; i < dofs_per_cell; ++i)
        {
            // update cellContext with the current i loop dof index

            for(unsigned int j = 0; j < dofs_per_cell; ++j)
            {
                // update cellContext with the current j loop dof index

                cell_matrix(i,j) += evaluate(cell_contribution_to_stiffness_matrix, cellContext);
            }

            // update cellContext with the current i, j, q indices

            cell_rhs(i) += evaluate(cell_contribution_to_load_vector, cellContext);
        }
    }
}

Apart from specifying the weak formulation using the single expressions evaluated in a generic loop as shown above, tuples of expressions representing independent terms evaluated in the q, i and j loops can be also used. This is useful for complex weak form expressions involving non-linear terms (i.e. dof approximations) and results in simpler expressions and a faster evaluation. The usage of separate items for different loop iterations is illustrated below:

FEValues<dim>             fe_values(...);
std::vector<unsigned int> local_dof_indices;

// Create evaluation context objects where the expressions will be evaluated:
feCellContext<dim> cellContext(fe_values, local_dof_indices, ...);

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                               endc = dof_handler.end();
for(; cell != endc; ++cell)
{
    // Update the evaluation context with the current cell.
    // It will call fe_values.reinit(cell), cell->get_dof_indices(local_dof_indices) etc.
    update(cellContext, cell);

    for(unsigned int q_point = 0; q_point < n_q_points; ++q_point)
    {
        // Temporary storage
        Vector<adouble>     cell_rhs_temp(dofs_per_cell);
        FullMatrix<adouble> cell_matrix_temp(dofs_per_cell, dofs_per_cell);

        // update cellContext with the current quadrature point index

        q_loop_term = evaluate(q_loop_expression, cellContext);

        for(unsigned int i = 0; i < dofs_per_cell; ++i)
        {
            // update cellContext with the current i loop dof index

            i_loop_term = evaluate(i_loop_expression, cellContext);

            for(unsigned int j = 0; j < dofs_per_cell; ++j)
            {
                // update cellContext with the current j loop dof index

                j_loop_term = evaluate(j_loop_expression, cellContext);

                cell_matrix_temp(i,j) += i_loop_term * j_loop_term;
            }

            cell_rhs_temp(i) += i_loop_term;
        }

        cell_rhs_temp    *= q_loop_term;
        cell_matrix_temp *= q_loop_term;
    }

    // Add the temporary data to the cell matrix/vector
    cell_rhs    += cell_rhs_temp;
    cell_matrix += cell_matrix_temp;
}

Obviously, the generic loops can be used to solve many FE problems but not all. However, they can support a large number of problems at the moment. In the future they will be expanded to support a broader class of problems.

Equations produced by the Finite Element matrix assembly can be extremely long/complex. Under the hood, they will be simplified as much as possible. Some examples are given below.

• From Tutorial deal.II 2:

  

• From Tutorial deal.II 3:

  

• From Tutorial deal.II 7:

  

DAE Tools provide four main classes for Finite Element models using deal.II library:

1. dealiiFiniteElementDOF_1D, dealiiFiniteElementDOF_2D and dealiiFiniteElementDOF_3D

   In deal.II it represents a degree of freedom distributed on a finite element domain. In DAE Tools it represents a variable distributed on a finite element domain.

2. dealiiFiniteElementSystem_1D, dealiiFiniteElementSystem_2D and dealiiFiniteElementSystem_3D (implements daeFiniteElementObject)

   It is a wrapper around deal.II FESystem<dim> class and handles all finite element related details. It uses information about the mesh, quadrature and face quadrature formulas, degrees of freedom and the FE weak formulation to assemble the system’s mass matrix (Mij), stiffness matrix (Aij) and the load vector (Fi).

3. dealiiFiniteElementWeakForm_1D, dealiiFiniteElementWeakForm_2D and dealiiFiniteElementWeakForm_3D

   Contains weak form expressions for the contribution of FE cells to the system/stiffness matrices, the load vector, boundary conditions and (optionally) surface/volume integrals (as an output).

4. daeFiniteElementModel

   daeModel-derived class that use system matrices/vectors from the dealiiFiniteElementSystem_nD object to generate a system of equations.

A typical use-case scenario consists of the following steps:

1. The starting point is a definition of the dealiiFiniteElementSystem_nD class (where nD can be 1D, 2D or 3D). That includes specification of:

   • Mesh file in one of the formats supported by deal.II (GridIn)
   • Degrees of freedom (as a python list of dealiiFiniteElementDOF_nD objects). Every dof has a name which will be also used to declare DAE Tools variable with the same name, description, finite element space FE (deal.II FiniteElement<dim> instance) and the multiplicity
   • Quadrature formulas for elements and their faces

2. Creation of daeFiniteElementModel object (similarly to the ordinary DAE Tools model) with the finite element system object as the last argument.

3. Definition of the weak form of the problem using the following functions (a version exist for 1D, 2D and 3D):

   • phi_1D(): corresponds to shape_value in deal.II
   • dphi_1D(): corresponds to shape_grad in deal.II
   • d2phi_1D(): corresponds to shape_hessian in deal.II
   • phi_vector_1D(): corresponds to shape_value of vector dofs in deal.II
   • dphi_vector_1D(): corresponds to shape_grad of vector dofs in deal.II
   • d2phi_vector_1D(): corresponds to shape_hessian of vector dofs in deal.II
   • div_phi_1D(): corresponds to divergence in deal.II
   • JxW_1D(): corresponds to the mapped quadrature weight in deal.II
   • xyz_1D(): returns the point for the specified quadrature point in deal.II
   • normal_1D(): corresponds to the normal_vector in deal.II
   • function_value_1D(): wraps Function<dim> object that returns a value
   • function_gradient_1D(): wraps Function<dim> object that returns a gradient
   • function_adouble_value_1D(): wraps Function<dim,adouble> object that returns adouble value
   • function_adouble_gradient_1D(): wraps Function<dim,adouble> object that returns adouble gradient
   • dof_1D(): returns daetools variable at the given index (adouble object)
   • dof_approximation_1D(): returns FE approximation of a quantity as a daetools variable (adouble object)
   • dof_gradient_approximation_1D(): returns FE gradient approximation of a quantity as a daetools variable (adouble object)
   • dof_hessian_approximation_1D(): returns FE hessian approximation of a quantity as a daetools variable (adouble object)
   • vector_dof_approximation_1D(): returns FE approximation of a vector quantity as a daetools variable (adouble object)
   • vector_dof_gradient_approximation_1D(): returns FE approximation of a vector quantity as a daetools variable (adouble object)
   • adouble_1D(): wraps any daetools expression to be used in matrix assembly
   • tensor1_1D(): wraps deal.II Tensor<rank=1>
   • tensor2_1D(): wraps deal.II Tensor<rank=2>
   • tensor3_1D(): wraps deal.II Tensor<rank=3>
   • adouble_tensor1_1D(): wraps deal.II Tensor<rank=1,adouble>
   • adouble_tensor2_1D(): wraps deal.II Tensor<rank=2,adouble>
   • adouble_tensor3_1D(): wraps deal.II Tensor<rank=3,adouble>

   and constants: fe_i, fe_j and fe_q used to access the current indexes in the DOF and quadrature points loops.

The weak form contains the following contributions:

• Aij - cell contribution to the system stiffness matrix
• Mij - cell contribution to the system mass matrix
• Fi - cell contribution to the system load vector
• boundaryFaceAij - boundary face contribution to the system stiffness matrix
• boundaryFaceFi - boundary face contribution to the load vector
• innerCellFaceAij - inner cell face contribution to the system stiffness matrix
• innerCellFaceFi - inner cell face contribution to the load vector
• functionsDirichletBC - Dirichlet boundary conditions
• surfaceIntegrals - surface integrals
• volumeIntegrals - volume integrals

Example for the Helmholtz problem (step-7 tutorial from deal.II). The strong form of the Helmholtz equation is given by (with Neumann boundary conditions):

The weak form of the above equation is:

Cell contributions to the stiffness matrix and the load vector are:

The c++ implementation in deal.II is given in the step-7 example and in the previous c++ listing. In DAE Tools the above equation is specified in the following way (\(\delta \Omega\) boundary is marked with id = 0 and for simplicity, \(f\) and \(g\) are assumed constant):

class modTutorial(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)

        dofs = [dealiiFiniteElementDOF_2D(name='u',
                                          description='u description',
                                          fe = FE_Q_2D(1),
                                          multiplicity=1)]

        self.fe_system = dealiiFiniteElementSystem_2D(meshFilename    = '...',        # path to the .msh file
                                                      quadrature      = QGauss_2D(3), # quadrature formula
                                                      faceQuadrature  = QGauss_1D(3), # face quadrature formula
                                                      dofs            = dofs)         # degrees of freedom

        self.fe_model = daeFiniteElementModel('Helmholtz', self, 'Helholtz problem', self.fe_system)

    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        # Create some auxiliary objects for readability
        phi_i  = phi_2D ('u', fe_i, fe_q)
        phi_j  = phi_2D ('u', fe_j, fe_q)
        dphi_i = dphi_2D('u', fe_i, fe_q)
        dphi_j = dphi_2D('u', fe_j, fe_q)
        JxW    = JxW_2D(fe_q)
        f      = 1.0
        g      = 1.0

        weakForm = dealiiFiniteElementWeakForm_2D(Aij = (dphi_i*dphi_j + phi_i*phi_j) * JxW,
                                                  Mij = 0.0,
                                                  Fi  = phi_i * f * JxW,
                                                  boundaryFaceFi = {0 : phi_i * g * JxW} )

        self.fe_system.WeakForm = weakForm

In DAE Tools v1.7.1 an additional way of specifying the weak formulation has been introduced. Now, weak form contributions can also be python lists containing feExpression objects or tuples of feExpression objects. Tuples can contain two or three items. Matrix contributions contain three items representing q-loop, i-loop and j-loop expressions in the deal.ii matrix assembly. Load vector contributions contain two items representing q-loop and i-loop expressions in the deal.ii vector assembly. This way the expressions produced by the finite element system assembly can be much simpler and the simulations faster. The new way of specifying weak formulations is useful mostly for specification of non-linear terms (i.e. those including scalar or vector DOF approximations).

For instance, the convection term in the heat convection-diffusion equation (tutorial_dealii_7.py) can be specified in a more efficient way:

class modTutorial(daeModel):
    ...
    def DeclareEquations(self):
        ...

        # Contributions from the Stokes equation:
        Aij_u_gradient     = ...
        Aij_p_gradient     = ...
        # Using the new way (q-loop, i-loop tuple for the load vector contributions):
        Fi_buoyancy        = (T_dof, -rho * beta * (gravity * phi_vector_u_i) * JxW)

        # Contributions from the continuity equation:
        Aij_continuity     = ...

        # Contributions from the heat convection-diffusion equation:
        Mij_T_accumulation = ...
        # Using the new way (q-loop, i-loop, j-loop tuple for matrix contributions):
        Aij_T_convection   = (u_dof, phi_T_i, dphi_T_j * JxW)
        Aij_T_diffusion    = ...
        Fi_T_source        = ...

        # Total contributions (using the new way - python lists of expressions or tuples):
        Mij = [Mij_T_accumulation]
        Aij = [Aij_u_gradient + Aij_p_gradient + Aij_continuity + Aij_T_diffusion,  Aij_T_convection]
        Fi  = [Fi_T_source, Fi_buoyancy]

        weakForm = dealiiFiniteElementWeakForm_2D(Aij = Aij,
                                                  Mij = Mij,
                                                  Fi  = Fi)

Using the old way the vector DOF approximation (\(u_{dof}\)) is evaluated \(N_{quadrature\,points} \cdot N_{dofs\,per\,cell} \cdot N_{dofs\,per\,cell}\) times while only \(N_{quadrature\,points}\) evaluations are required. A vector dof approximation is an expensive operation calculated in the following way:

where \(\phi^u(j,q)\) is a rank=1 tensor. Reduction in the number of evaluations by \(N_{dofs\,per\,cell}^2\) leads to the huge amounts of memory and computation time saved.

More information about the finite element method in DAE Tools can be found in the API reference and in Finite Element Tutorials (in particular Tutorial deal.II 1).