DAE Systems

In Dyssol, you can solve systems of DAE automatically. In this case, the unit should contain one or several additional objects of class CDAEModel. This class is used to describe DAE systems and can be automatically solved with class CDAESolver.

DAE model

CDAEModel()

Basic constructor. Creates an empty DAE model.

Variables

unsigned AddDAEVariable(bool _bIsDifferentiable, double _dVariableInit, double _dDerivativeInit, double _dConstraint)

Adds new algebraic (_bIsDifferentiable = false) or differential (_bIsDifferentiable = true) variable with initial values _dVariableInit and _dDerivativeInit.

Should be called in function Initialize(Time) of the unit. Returns unique index of added variable.

_dConstraint sets the constraint for the variable, different values of _dConstraint are defined as follows:

0.0: no constraint

1.0: Variable ≥ 0.0

−1.0: Variable ≤ 0.0

2.0: Variable > 0.0

−2.0: Variable < 0.0

unsigned GetVariablesNumber()

Returns current number of defined variables.

void ClearVariables()

Removes all defined variables from the model.

Tolerance

void SetTolerance(double _dRTol, double _dATol)

Sets values of relative and absolute tolerances for all variables.

void SetTolerance(double _dRTol, std::vector<double> &_vATol)

Sets value of relative tolerance for all variables and absolute tolerances for each variable separately.

double GetRTol()

Returns current relative tolerance.

double GetATol(unsigned _dIndex)

Returns current absolute tolerance for specified variable.

Virtual functions

virtual void CalculateResiduals(double _dTime, double *_pVars, double *_pDerivs, double *_pRes, void *_pUserData)

Computes the problem residual for given values of the independent variable _dTime, state vector _pVars, and derivatives _pDerivs. Here the DAE system should be specified in implicit form. Function will be called by solver automatically.

_dTime: current value of the independent variable $t$ .

_pVars: pointer to an array of the dependent variables, $y(t)$ .

_pDerivs: pointer to an array of derivatives $y'(t)$ .

_pRes: output residual vector $F(t, y, y')$ .

_pUserData: pointer to user’s data. Is used to provide access from this function to unit’s data.

virtual void ResultsHandler(double _dTime, double *_pVars, double *_pDerivs, void *_pUserData)

Processing the results returned by the solver at each calculated step. Called by solver every time when the solution in new time point is ready.

_dTime: current value of the independent variable $t$ .

_pVars: current values of the dependent variables, $y(t)$ .

_pDerivs: current values of derivatives $y'(t)$ .

_pUserData: pointer to user’s data. Is used to provide access from this function to unit’s data.

Other functions

void Clear()

Removes all data from the model.

void SetUserData(void *_pUserData)

Set pointer to user’s data. This data will be returned in overloaded functions CalculateResiduals and ResultsHandler. Usually is used to provide access from these functions to unit’s data.

DAE solver

CDAESolver()

Basic constructor. Creates an empty solver.

Model

bool SetModel(CDAEModel *_pModel)

Sets model to a solver. Should be called in function Initialize of the unit.

Returns false on error.

bool SetMaxStep()

Sets maximum time step for solver. Should be used in Initialize before the function SetModel.

bool Calculate(double _dTime, unsigned _nModel = KIN_NONE)

Solves problem on a given time point _dTime. Should be called in function Simulate of the unit.

Default value of _nModel is KIN_NONE, a Newton-based iteration method. Other available models can be found in declaration for KINSOL parameters in the *.cpp file, including: KIN_FP (fixed point), KIN_LINESEARCH and KIN_PICARD.

Returns false on error.

bool Calculate(double _dStartTime, double _dEndTime)

Solves problem on a given time interval. Should be called in function Simulate of the unit. Returns false on error.

Other functions

void SaveState()

Saves current state of the solver. Should be called in function SaveState of the unit.

void LoadState()

Loads last saved state of the solver. Should be called in function LoadState of the unit.

std::string GetError()

Returns error’s description. Can be called if function SetModel or Calculate returns false.

Application example

Assume that you are going to solve a dynamic DAE system,

$\begin{cases} \dfrac{dx}{dt} = -0.04 \cdot x + 10^4 \cdot y \cdot z &x(0) = 0.04 \dfrac{dy}{dt} = 0.04\cdot x - 10^4 \cdot y \cdot z - 3 \cdot 10^7 \cdot y^2 &y(0) = -0.04 x + y + z = 1 \end{cases}$

where $x$ , $y$ , $z$ are fractions of solid, liquid and vapor phases of the output stream of the unit.

Now you want to develope a new unit for for automatic calculation of this DAE by using built-in solvers of Dyssol, just do the following steps:

Add a new template unit to your solution and name it DynamicUnitWithSolver, please refer to Unit development.

In file Unit.h, there is already a class of DAE model CMyDAEModel defined with:

two functions, which must be overridden (CalculateResiduals and ResultsHandler) and
a class of unit CUnit with two additional variables (for DAE model CMyDAEModel m_Model, for DAE solver CDAESolver m_Solver).

Add three variables in class CMyDAEModel to store indices for $x$ , $y$ and $z$ variables, name them m_nS, m_nL and m_nV respectively.

After adding description of class CMyDAEModel, your code should look like this:

class CMyDAEModel : public CDAEModel
{
public:
        unsigned m_nS; //solid fraction
        unsigned m_nL; //liquid fraction
        unsigned m_nV; //vapor fraction

public:
        void CalculateResiduals(double _dTime, double* _pVars, double* _pDers, double* _pRes, void* _pUserData);
        void ResultsHandler(double _dTime, double* _pVars, double* _pDers, void *_pUserData);
};

Setup unit’s basic info (name, author’s name, unique ID, ports) as described in Unit development.Then provide model with pointer to your unit, in order to have an access to the unit from functions of the model. Now your unit’s constructor should look like this:

CUnit::CUnit()
{
        // Unit basic information
        m_sUnitName = "DummyUnit4";
        m_sAuthorName = "Author";
        m_sUniqueID = "344BCC0048AA4c3a9117F20A9F8AF9A8"; //an example ID

        // Add ports for in- and outlets
        AddPort("InPort", INPUT_PORT);
        AddPort("OutPort", OUTPUT_PORT);

        // Set this unit as user data of the model applied
        m_Model.SetUserData(this);
}

Implement function Initialize of the unit:

Since function Initialize is called every time when simulation starts, all variables must be previously removed from the model by calling ClearVariables.

m_Model.ClearVariables();

Now you can add 3 variables with specified initial conditions to the model (using function AddDAEVariable) according to the equation system. Use phase fractions of the input stream as initials.

Set tolerances to the model using function SetTolerance. As tolerances for the model, global tolerances of the system can be used.

m_Model.SetTolerance( GetRelTolerance(), GetAbsTolerance() );

Now, you can connect your model to the solver by calling function SetModel. To receive errors from solver, connect it to the global errors handling procedure.

if( !m_Solver.SetModel(&m_Model) )
        RaiseError(m_Solver.GetError());

Your code should look like this after following all steps above:

void CUnit::Initialize(double _dTime)
{
        // Get pointer to inlet stream
        CMaterialStream* pInpStream = GetPortStream("InPort");

        // Clear previous state variables in model
        m_Model.ClearVariables();

        // Add state variables to model
        double dSFr = pInpStream->GetSinglePhaseProp(_dTime, FRACTION, SOA_SOLID);
        double dLFr = pInpStream->GetSinglePhaseProp(_dTime, FRACTION, SOA_LIQUID);
        double dVFr = pInpStream->GetSinglePhaseProp(_dTime, FRACTION, SOA_VAPOR);
        m_Model.m_nS = m_Model.AddDAEVariable(true, dSFr, 0.04, 1.0);
        m_Model.m_nL = m_Model.AddDAEVariable(true, dLFr, -0.04, 1.0);
        m_Model.m_nV = m_Model.AddDAEVariable(false, dVFr, 0, 1.0);

        // Set tolerance to model
        m_Model.SetTolerance( GetRelTolerance(), GetAbsTolerance() );

        // Set model to solver
        if( !m_Solver.SetModel(&m_Model) )
                RaiseError(m_Solver.GetError());
}

Connect solver to a system saving / loading procedure in functions SaveState() and LoadState() of the unit:

void CUnit::SaveState()
{
        m_Solver.SaveState();
}

void CUnit::LoadState()
{
        m_Solver.LoadState();
}

In function Simulate of the unit, calculation procedure should be run by calling function Calculate of the solver. Additionally, solver can be connected to the system’s errors handling procedure to receive possible errors during the calculation. Unit’s Simulate function after that must look as follows:

void CUnit::Simulate(double _dStartTime, double _dEndTime)
{
        // Get pointers to inlet and outlet streams
        CMaterialStream *pInputStream = GetPortStream("InPort");
        CMaterialStream *pOutputStream = GetPortStream("OutPort");

        pOutputStream->RemoveTimePointsAfter(_dStartTime, true);
        pOutputStream->CopyFromStream(pInputStream, _dStartTime); // Copy the inlet stream information to outlet stream

        // Run solver and check if errors take place
        if( !m_Solver.Calculate(_dStartTime, _dEndTime) )
                RaiseError( m_Solver.GetError() );
}

In function CalculateResiduals, DAE system in implicit form should be described. According to the given equation system

$\begin{cases} \dfrac{dx}{dt} = -0.04 \cdot x + 10^4 \cdot y \cdot z &x(0) = 0.04 \dfrac{dy}{dt} = 0.04\cdot x - 10^4 \cdot y \cdot z - 3 \cdot 10^7 \cdot y^2 &y(0) = -0.04 x + y + z = 1 \end{cases}$

and the definition of residual _pRes (difference between new and old values from last calculation), the equaitons are expressed below. Here the values are arranged in vector _pVars according to sequence $[x, y, z]$ , in _pDers according to $[dx, dy]$ .

void CMyDAEModel::CalculateResiduals(double _dTime, double* _pVars, double* _pDers, double* _pRes, void* _pUserData)
{
        _pRes[0] = _pDers[0] - (-0.04*_pVars[0] + 1.0e4*_pVars[1]*_pVars[2]);
        _pRes[1] = _pDers[1] - ( 0.04*_pVars[0] - 1.0e4*_pVars[1]*_pVars[2] - 3.0e7*_pVars[1]*_pVars[1] );
        _pRes[2] = _pVars[0] + _pVars[1] + _pVars[2] - 1;
}

Last step is handling of results from the solver (CMyDAEModel::ResultsHandler). Calculated fractions can be set here to the output stream of the unit. To access to the unit’s data, use the pointer _pUserData defined previously:

void CMyDAEModel::ResultsHandler(double _dTime, double* _pVars, double* _pDerivs, void *_pUserData)
{
        // Get pointers to streams
        CUnit *unit = static_cast<CUnit*>(_pUserData);
        CMaterialStream *pStream = static_cast<CMaterialStream*>(unit->GetPortStream("OutPort"));

        // Add time point to outlet stream
        pStream->AddTimePoint(_dTime);

        // Set calculated results to corresponding phases
        pStream->SetSinglePhaseProp(_dTime, FRACTION, SOA_SOLID, _pVars[0]);
        pStream->SetSinglePhaseProp(_dTime, FRACTION, SOA_LIQUID, _pVars[1]);
        pStream->SetSinglePhaseProp(_dTime, FRACTION, SOA_VAPOR, _pVars[2]);
}