# 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:

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

2. 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);
       };

3. 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);
       }

4. 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());
       }

5. 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();
       }

6. 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() );
       }

7. 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;
       }

8. 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]);
       }