# Tutorial 2 - Two tanks system

## Background information
This tutorial shows how to use _numerous_ to create a system of multiple components.
The tutorial is aimed at demonstrating the usability and scalability of the model architecture for systems with multiple physical components connected bewteen each other by means of connectors.

The tutorial is built based on a relative simple system, e.g. two tanks are placed on top of each other and connected by a valve. The implementation using _numerous_ might therefore seems tedious and over-complicated, given the limited number of components and thus equations to solve. However, as already pintpointed, we wish to illustrate how systems can be created, and how _numerous_ model architecture would be advantegeous for systems with a high degree of complexity.

The model implementation is developed by instantiating _items_ of different classes, and connecting them by means of _connectors_.

The tutorial can be run only importing the extra packages pandas, numpy and plotly.

## System governing equations

A system of two tanks and a valve is shown in the figure below. H indicates the height of a certain mass M of liquid in each tank. The liquid density is indicated by $\rho$, while the cross-sectional area of each tank is indicated by the parameter $A_0$.The mass flow rate through the valve is determined by a valve characteristic parameter  $\mathrm{C_v}$.


At time 0, the top tank (tank 1) is filled with a certain amount of liquid, determining the height $\mathrm{H_1}$, while the bottom tank (tank 2) is completely empty. If the valve is opened at time 0, a mass flow rate starts flowing through the valve between tank 1 and tank 2. After a given time (that we indicate as $t_\mathrm{filling}$) the top tank will be empty, while the bottom tank will contain the entire mass of liquiq.

The mathematical model describing this systme is illusrated below.

For the sake of simplicity, we assume that the mass flow rate through the valve is proportional to the pressure drop accross the valve by the relation:
$$\begin{equation}
\dot{\mathrm{m}} = \mathrm{C_v} \sqrt{\mathrm{\Delta p}}
\end{equation}
$$


Moreover, the conservation of mass in the two tanks is expressed by the equations (with the assumption of incompressible fluid):
$$\begin{equation}
\rho A_{0,1} \mathrm{\dfrac{dH_1}{dt}} = - \dot{m}
\end{equation}
$$
$$\begin{equation}
\rho A_{0,2}  \mathrm{\dfrac{dH_2}{dt}} = + \dot{m}
\end{equation}
$$


The pressure drop accross the valve can be expressed by the equation:
$$\begin{equation}
\Delta p = \rho g H_1
\end{equation}
$$

By substituting this latter expression in the previous equations, the differential system of equations describing the system becomes:
$$\begin{equation}
\rho A_{0,1} \mathrm{\dfrac{dH_1}{dt}} = - \mathrm{C_v} (t)\sqrt{\rho g \mathrm{H_1}}
\end{equation}
$$
$$\begin{equation}
\rho A_{0,2} \mathrm{\dfrac{dH_2}{dt}} = + \mathrm{C_v}(t)\sqrt{\rho g \mathrm{H_1}}
\end{equation}
$$

## Analytical solution
The system of differential equations above can be quite complex to solve analytically given the time dependence of $\mathrm{C_v}$. Therefore, a numerical solver can be used to solve the system, and we will show in this tutorial how to create a model and to solve it using _numerous_.

However, in order to prove the correct implementation of the solver, we will compare the results of the numerical solution against the analytical solution for the simplyfing case of having a characteristic valve parameter $\mathrm{C_v}$ independent from time.

If we equal the parameter to a constant value over time  $\mathrm{C_v}(t) = \mathrm{C_v}$, it is possible to solve the system analytically. In fact, the integration of the previous equations leads to the solution:
$$\begin{equation}
\mathrm{H_1(t)} = \Bigg[\sqrt{\mathrm{H_{1,0}}} - \dfrac{\mathrm{C_v}}{2 A_{0,1}} \sqrt{\dfrac{g}{\rho}} \cdot t \Bigg]^{2}
\end{equation}
$$
$$\begin{equation}
\mathrm{H_2(t)} = \mathrm{H_{2,0}} + \dfrac{\mathrm{C_v}}{A_{0,2}} \sqrt{ \dfrac{g}{\rho}} \sqrt{\mathrm{H_{1,0}}} \cdot t - \dfrac{\mathrm{C_v^{2}}}{4 A_{0,1} A_{0,2}} \dfrac{g}{\rho} \cdot t^2
\end{equation}
$$

The solution can also be expressed as function of the total mass in the tanks, which is related to the liquid height by the volume as:
$$\begin{equation}
\mathrm{M} = \rho \mathrm{H} \mathrm{A_0}
\end{equation}
$$

The analytical solution would thus be expressed as:
$$\begin{equation}
\mathrm{M_1(t)} = \Bigg[\sqrt{\rho \mathrm{H_{1,0} A_{0,1}}} - \dfrac{\mathrm{C_v}}{2} \sqrt{\dfrac{g}{A_{0,1}}} \cdot  t \Bigg]^{2}
\end{equation}
$$
$$\begin{equation}
\mathrm{M_2(t)} = \mathrm{\rho H_{2,0}} A_{0,2}  + \mathrm{C_v} \sqrt{ \dfrac{g}{A_{0,1}}} \sqrt{\mathrm{\rho \mathrm{H_{1,0} A_{0,1}}}} \cdot t - \dfrac{\mathrm{C_v^{2}}}{4} \dfrac{g}{A_{0,1}} \cdot  t^2
\end{equation}
$$

## Input data
The input data used for the analytical solution and the simulation are reported below:
* liquid density, $\rho$ = 1000 kg/m$^3$
* cross-sectional area tank 1, $A_{0,1}$ = 0.05 m$^2$
* cross-sectional area tank 2, $A_{0,2}$ = 0.1 m$^2$
* initial conditions:
    * initial height tank 1  $H_{0,1}$ = 0.5 m
    * initial height tank 2  $H_{0,2}$ = 0.0 m
    

The implementation of the model with a time-dependent $\mathrm{C_v}$ was carried out by using the following formulation:
* $\mathrm{C_v} = \mathrm{C_{v,0}} \cdot \big[sin{(\omega \cdot t)}+1\big]  $
and imposing:
    * $\mathrm{C_{v,0}} = 0.1 \; \; \; (\mathrm{kg}/\mathrm{s})/\mathrm{Pa}$
    * $\omega \; \; \;  = 1.5 \; \; \;\mathrm{rad/s}$

A constant value of the valve parameter $\mathrm{C_v}$ was instead obtained by imposing $\omega = 0$, so that:
* $\mathrm{C_v} = \mathrm{C_{v,0}}$ 

at any time step, and the solution could be compared against the analytical formulation derived above.

This is translated into the followinf python code:

### Definition of input data - code:


```python
import numpy as np
## Input data
# Define liquid and geometry for the case in analysis
rho_water = 1000           # Water density [kg/m3]
A0        = [0.05, 0.1]    # Tanks cross-sectional area [m2]
Cv_0      = 0.1            # [(kg s)/Pa] Amplitude of valve characteristic parameter
g         = 9.81           # [m/s2] gravitational accelleration

## Initial conditions
H0        = [0.5, 0]     # Initial condition - liquid height [m]
# Estimate the initial value of total mass M [kg] in the two tanks
M0        = rho_water * np.multiply(A0, H0)

# Caclulate what is the time (analytical solution) to entirely fill the tank (simple case - analytical sol)
t_filling = np.sqrt(A0[0]/g) * np.sqrt(M0[0]) *2* 1/Cv_0

# Define start and stop time for the analysis and simulation
t_start_sim = 0
t_stop_sim  = 10
```

### Analytical soltion: Python code
First of all, we define a function which is able to return the analytical solution for a given time span, given the input parameters. The implementation below is based on the analytical solution derived above.


```python
import pandas as pd

N_t             = 11                           # Number of time steps at which the analytical solution is evaluated
time_range      = [t_start_sim ,t_stop_sim]    # Time span in which the solution is plotted

def analytical_solution(N_t=N_t, time_range=time_range, g=9.81,  A0=[0.05, 0.1], Cv_0=0.1,
                         M0=[25, 0], H0=[0.5, 0] , rho_water=1000):
    time_vector       = np.linspace(start = time_range[0], stop = time_range[1], num = N_t)
    
    # Create a dictionary with time, evolution of liquid heights and mass in the two tanks for the analytical solution:
    result_analytical = {'t': [], 'M_1': [], 'M_2': [], 'H_1': [], 'H_2': []}
    
    # Calculate what is the time after which the first tank is empty and the second tank is full
    t_filling = np.sqrt(A0[0]/g) * np.sqrt(M0[0]) *2* 1/Cv_0
    
    for i in range(N_t):
        if time_vector[i]  < t_filling:
            M_1_analytical = (np.sqrt(M0[0]) - Cv_0 / 2 * np.sqrt(g / A0[0]) * time_vector[i] ) ** 2
            M_2_analytical = np.sqrt(M0[1]) + Cv_0 * np.sqrt(g / A0[0]) * np.sqrt(
                M0[0]) * time_vector[i]  - Cv_0 ** 2 * g / A0[0] / 4 * (time_vector[i]  ** 2)    
        
            H_1_analytical = (np.sqrt(H0[0]) - Cv_0 / 2 /A0[0] * np.sqrt(g/ rho_water) * time_vector[i] ) ** 2
            H_2_analytical = np.sqrt(H0[1]) + Cv_0 /A0[1] * np.sqrt(g/ rho_water) * np.sqrt(
                H0[0]) * time_vector[i]  - Cv_0 ** 2 * g / A0[0] / 4 /A0[1]/rho_water* (time_vector[i]  ** 2)
        else:
            M_2_analytical = M0[0] + M0[1]
            M_1_analytical = 0
            H_2_analytical = (M0[0] + M0[1])/rho_water/A0[1]
            H_1_analytical = 0
            
        result_analytical['t'].append(time_vector[i])
        result_analytical['M_1'].append(M_1_analytical)
        result_analytical['M_2'].append(M_2_analytical)
        result_analytical['H_1'].append(H_1_analytical)
        result_analytical['H_2'].append(H_2_analytical)
    return result_analytical
```

We can run the equation to get the analytical soluton for the input data defined above, and we can create a table (using pandas data frame) containing the analytical solution at each evaluated time. The code is shown below.


```python
result_analytical = analytical_solution(g=g, A0=A0, M0=M0, rho_water = rho_water, N_t=N_t, time_range=time_range, Cv_0=Cv_0)

data = {'Time, s': result_analytical['t'], 'H1, m':result_analytical['H_1'], 'H2, m':result_analytical['H_2'], 
        'M1, kg':result_analytical['M_1'],'M2, kg':result_analytical['M_2']}
pd.DataFrame(data)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Time, s</th>
      <th>H1, m</th>
      <th>H2, m</th>
      <th>M1, kg</th>
      <th>M2, kg</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.0</td>
      <td>0.500000</td>
      <td>0.000000</td>
      <td>25.000000</td>
      <td>0.000000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1.0</td>
      <td>0.369739</td>
      <td>0.065131</td>
      <td>18.486929</td>
      <td>6.513071</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2.0</td>
      <td>0.259097</td>
      <td>0.120451</td>
      <td>12.954859</td>
      <td>12.045141</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3.0</td>
      <td>0.168076</td>
      <td>0.165962</td>
      <td>8.403788</td>
      <td>16.596212</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4.0</td>
      <td>0.096674</td>
      <td>0.201663</td>
      <td>4.833718</td>
      <td>20.166282</td>
    </tr>
    <tr>
      <th>5</th>
      <td>5.0</td>
      <td>0.044893</td>
      <td>0.227554</td>
      <td>2.244647</td>
      <td>22.755353</td>
    </tr>
    <tr>
      <th>6</th>
      <td>6.0</td>
      <td>0.012732</td>
      <td>0.243634</td>
      <td>0.636577</td>
      <td>24.363423</td>
    </tr>
    <tr>
      <th>7</th>
      <td>7.0</td>
      <td>0.000190</td>
      <td>0.249905</td>
      <td>0.009506</td>
      <td>24.990494</td>
    </tr>
    <tr>
      <th>8</th>
      <td>8.0</td>
      <td>0.000000</td>
      <td>0.250000</td>
      <td>0.000000</td>
      <td>25.000000</td>
    </tr>
    <tr>
      <th>9</th>
      <td>9.0</td>
      <td>0.000000</td>
      <td>0.250000</td>
      <td>0.000000</td>
      <td>25.000000</td>
    </tr>
    <tr>
      <th>10</th>
      <td>10.0</td>
      <td>0.000000</td>
      <td>0.250000</td>
      <td>0.000000</td>
      <td>25.000000</td>
    </tr>
  </tbody>
</table>
</div>



## Implementation using _numerous_

### Preliminary steps
The first step for the implementation is to include all the relevant _numerous_ modules. 
For this tutorial we need:
* _Item_ for defining item objects
* _Model_ for defining model objects
* _ConnectorTwoWay_ for defining a special connectors object
* _Subsystems_ for defining the subsystem object
* _Simulation_ for defining the simulation
* _Equation_ and equation decorator for objects of equation class

Moreover, we will need HistoryDataFrame to store results.


```python
# We include all the relevant modules from numerous:
from numerous.engine.system import Item
from numerous.engine.model import Model
from numerous.engine.system import Subsystem
from numerous.engine.simulation import Simulation
from numerous.engine.system import ConnectorTwoWay
from numerous.engine.variables import VariableType, VariableDescription, OverloadAction

from numerous.multiphysics import Equation
from numerous.multiphysics import EquationBase

```

### Define Tank Equation
The first item that we will model is the tank. Before creating the tank item, we need to define the equation to apply, and we thus create a Tank_Equation item using Equation class. We need to define all the parameters and constants of the equation, by using the method ` Equation.add_parameter ` .  _g_ (gravitational accelleration) is the only variable defined as a constant, as its value cannot be modified.

Note that the line of code
``` python 
super().__init__(tag='for info only') 
```
in `__init__` is a necessary line to create any class using _numerous_ in the current release.

The only _state_ variable in the tank equation is given by the tank height, which is determined by the differential equation presented in the description above and reported in the equation definition.


```python
class Tank_Equation(EquationBase):
    def __init__(self, tag="tank_equation", H=0.5, rho=1000, A0=0.05, g=9.81):
        super().__init__(tag='tank_equation')

        self.add_state('H', H)          # [m] Liquid height in the tank
        self.add_parameter('rho', rho)  # [kg/m3] Liquid density
        self.add_parameter('A0', A0)    # [m2] Tank cross-sectional area
        self.add_constant('g', g)       # [m/s2] Gravitational acceleration
        self.add_parameter('mdot', 0)   # [kg/s] Mass flow rate 
        
    @Equation()
    def eval(self, scope):
        # Differential equation for mass conservation in a general tank with a mass flow rate entering or leaving the tan
        scope.H_dot = scope.mdot / scope.rho / scope.A0
```

### Define Valve Equation
The second item that we will model is the valve. Before creating the valve item, we need to define the equation to apply, and we thus create a `Valve_Equation` item using `EquationBase` class.

Please, note the use of the **global variable** time, which is defined as ` scope.globals.time ` in _numerous_ , and thus it has not to be specified as the other parameters and constants.

Note that in the `Valve_Equation` no state is defined, since the valve characteristic parameter is a function of time, but we have an explicit formulation for it. 


```python
# We define the equation (using EquationBase class) determining the mass flow rate across the valve
class Valve_Equation(EquationBase):    
    def __init__(self, Cv_0=0.1, rho=1000, g=9.81, omega=1):
        super().__init__(tag='valve_equation')
        
        self.add_parameter('omega', omega)   # [rad/sec]   Angular frequency of valve characteristic parameter
        self.add_parameter('Cv_0', Cv_0)     # [(kg/s)/Pa] Amplitude of valve characteristic parameter
        self.add_parameter('mdot1', 0)       # [kg/s]      Mass flow rate in one side of the valve
        self.add_parameter('mdot2', 0)       # [kg/s]      Mass flow rate in the other side of the valve
        self.add_parameter('H1', 0)          # [m]         Liquid height in the tank 1 connected to the valve (top tank)
        self.add_parameter('H2', 0)          # [m]         Liquid height in the tank 2 connected to the valve (bottom tank)
        self.add_parameter('rho', rho)       # [kg/m3]     Liquid density
        self.add_constant('g', g)            # [m/s2]      Gravitational acceleration
        self.add_parameter('Cv', Cv_0)       # [(kg/s)/Pa] Valve characteristic parameter
        
    @Equation()
    def eval(self,scope,global_variables):
        scope.Cv      = scope.Cv_0 * (np.sin(scope.omega * global_variables.time)+1)           #[(kg/s)/Pa]
        deltaP        = scope.rho * scope.g * (scope.H1)                                  #[Pa]
        mdot          = np.sign(deltaP) * np.sqrt(np.absolute(deltaP)) * scope.Cv         #[kg/s]
        # The valve will be associated with two mass flow rates (one leaving and one entering the component), 
        #which - for conservation of mass - have the same magnitude and opposite sign
        scope.mdot1   = -mdot                                                             #[kg/s]
        scope.mdot2   = mdot                                                              #[kg/s]
        
```

### Define Tank as Item
We define the `Tank` class as an `Item` class. We then create a namespace 'v1' to contain the variables for the `Tank_Equation`.
The equation is associated to the namespace using the `add_equations` method, as shown in the code below.


```python
class Tank(Item):
    def __init__(self, tag="tank", H=0.5, rho=1000, A0=0.05, g=9.81):
        super(Tank, self).__init__(tag)

        v1 = self.create_namespace('v1')
        v1.add_equations([Tank_Equation(H=H, rho=rho, A0=A0, g=g)])
```

### Define Valve as ConnectorTwoWay
Once that we have defined the equation describing the mass flow rate flowing through the valve, we need to create the Valve as a class `ConnectorTwoWay` and to assign an equation to it.
`ConnectorTwoWay` is a special case of a `Connector` class, and the reader is referred to _numerous_ documentation for an exhuastive explanation. The peculiarity of this connector is the possibility of defining two sides, i.e. variables can be binded to the connectors by specifying to different items as sides. 

In the code lines
``` python 
super().__init__(tag, side1_name='side1', side2_name='side2')
``` 
we have to specify the names of the two sides.

The steps that we have to take are the following (refer to the numbering in the code comments # to see which lines of codes belong to the different steps)
1. We create a namespace 'v1' to contain the variables for the valve equation. This is done using the Item method `Item.create_namespace `. The namespace is then associated to an equation using the `add_equations` method. 
2. We create variables at each side of the connector item, and we associated them to the same namespace containing the valve equation. The variables must be created because when we first instantiate the ConnectorTwoWay object no information on side1 and side2 is passed. 
3. The binding between the ConnectorTwoWay and the two items at each side is done, using the variables previosuly created in the name space. In this particular example:
    * the value of v1.H1 and v1.H2 (liquid heights of the tanks connected to the valve, stored inside the valve object) must point to the respective tank heights in the two side objects. This implies that the value of H is determined by the tank equation and not by the valve equation.
    * the value of the mass flow rate entering or leaving each tank (for example the value self.side1.v1.mdot stored inside the side1 object (tank 1)) must point to the mass flow rate flowing through the valve (in this case determined by the valve equation)


```python
# Define the valve as a connector item - connecting two tanks
class Valve(ConnectorTwoWay):
    def __init__(self, tag="valve", Cv_0=0.1, rho=1000, g=9.81, omega=0):
        super().__init__(tag, side1_name='side1', side2_name='side2')

        #1 Create a namespace for mass flow rate equation and add the valve equation
        v1 = self.create_namespace('v1') 
        v1.add_equations([Valve_Equation(Cv_0=Cv_0, rho=rho, g=g, omega=omega)])
   
        #2 Create variables H and mdot in side 1 adn 2
           #(side 1 represents a connection with one tank, with related liquid height H_1)
           #(side 1 represents a connection with the second tank, with related liquid height H_2)
        self.side1.v1.create_variable(name='H')
        self.side1.v1.create_variable(name='mdot')
        self.side2.v1.create_variable(name='H')
        self.side2.v1.create_variable(name='mdot')
        
        # Map variables between binding and internal variables for side 1 and 2
        # This is needed to update the values of the variables in the binding according to the equtions of the items
        v1.H1              = self.side1.v1.H
        v1.H2              = self.side2.v1.H
        
        self.side1.v1.mdot = v1.mdot1
        self.side2.v1.mdot = v1.mdot2
       
```

###  Create the sub-system of components 
After defining all the classes for the items that will consitute the system, we are ready for the system assembly.
We create a special class of `Subsystem`, inside which we: (refer to the numbering in the code comments # to see which lines of codes belong to the different steps)
1. create the gravitational accelleration constant and assign a value to it
2. create two instances of the class Tank called Tank_1 (top tank) and Tank_2 (bottom tank)
3. create one instance of the class Valve called Valve_1
4. bind Tank_1 and Tank_2 by assigning each of them to the two sides of Valve_1. we use the `ConnectorTwoWay.bind` method for this.
5. register the instanciated items in the Two_Tanks class 

The inputs needed to the subsystem are:
* H0, which is a vector containing the initial state of the system (initial liquid height of Tank_1 and Tank_2)
* Cv_0 and omega, which are amplitude and angular frequency of the valve characteristic parameter
We assume that geometry (A0) and liquid (rho) are given by the input data as fixed values


```python
# Define the subsystem composed by two tanks and one valve connecting them
class Two_Tanks(Subsystem):
    def __init__(self, tag, H0, Cv_0, omega):
        super().__init__(tag)

        #1. Gravitational acceleration
        g = 9.81
        
        #2. Instances of Tank class
        Tank_1 = Tank('tank_1', H=H0[0], rho=rho_water, A0=A0[0], g=g)
        Tank_2 = Tank('tank_2', H=H0[1], rho=rho_water, A0=A0[1], g=g)
        
        #3. Valve_1 is one instance of valve class
        Valve_1 = Valve('valve_1', Cv_0=Cv_0, rho=rho_water, g=g, omega=omega)
        
        #4. Binding
        Valve_1.bind(side1=Tank_1, side2=Tank_2)

        #5. Register all the instanciated items in the sub-system
        self.register_items([Tank_1, Tank_2, Valve_1])
```

###  Create the system model and simulation
Finally we are ready to define the model and the simulation of the implemented system, and we do it by creating a function named `t_1_item_model`. The inputs to the function are given by:
* H0, Cv_0, omega, which represent initial conditions of the tanks, and valve characteristic
* hdf, which is the historian where the simulation results will be stored
* t_start_sim and t_stop_sim determing the time span in which the system will be simulated

The steps are the following:
1. First a model object m1 is instantiated based on the `Two_Tanks` subsystem with given inputs
2. A simulation s1 is connected to the model object, and some solver settings are chosen. The parameter 'num' is used to specify the number of steps at which the solution is evaluated between t_start_sim and t_stop_sim
3. The `simulation.solve` method is called, and the solution is returned as output of the function


```python
def t_1_item_model(H0, Cv_0, omega, t_start_sim, t_stop_sim):
    
    # 1. Instantiate the model
    m1 = Model(Two_Tanks('subsystem', H0, Cv_0, omega))
    # 2. Setting up the simulation:
    s1  = Simulation(m1, t_start=t_start_sim, t_stop=t_stop_sim, num=5000)
    # 3. Call the solve method
    sol = s1.solve()
    return s1, sol
```

## Running the simulations - Case 1: Constant Cv 
We impose a Cv to be constant in time by setting the angular frequency $\omega$ to be equal to 0 (i.e. the cosine of 0 will be equal to 1 at each time step).
To obtain the solution of the system for Case 1, we simply need to call the `t_1_item_model` function with the inputs defined by the input data, and omega = 0.
The progress bar is shown below.


```python
omega = 0

# Collect the historical data at each time step
s1, sol = t_1_item_model(H0, Cv_0,omega, t_start_sim, t_stop_sim)
```

### Plotting: comparison vs. analytical solution
To plot the solution, we add the following code which uses `plotly`. We shall not go into detail on the code, but simply include it here as it's used for the purpose illustrating the results of this tutorial only.

We can plot the solution by accessing the `Model` object `historian`, which contains a log of all variables as a `Pandas` dataframe. The `Simulation` object, `sim1` contains the `model` object and the time logged variables are accessible through `sim1.model.historian`. 


```python

```


```python
# Plot the comparison
# Nb: the package plotly is needed for the plotting 
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import plotly.graph_objs as go

fig = go.Figure()
fig.add_trace(go.Scatter(
            x=hdf.df.index.total_seconds(), y=hdf.df['subsystem.tank_1.v1.H'],
            name='Tank 1 - numerical', mode='lines',line_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=result_analytical['t'], y=result_analytical['H_1'],
            name='Tank 1 - analytical', mode='markers', marker_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=hdf.df.index.total_seconds(), y=hdf.df['subsystem.tank_2.v1.H'],
            name='Tank 2 - numerical', mode='lines',line_color='rgba(152, 0, 0, .8)'))
fig.add_trace(go.Scatter(
            x=result_analytical['t'], y=result_analytical['H_2'],
            name='Tank 2 - analytical',  mode='markers',marker_color='rgba(152, 0, 0, .8)'))
              
fig.update_layout(title='Liquid height in tanks - analytical vs. numerical solution',
                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title="time, sec",
                 yaxis_title="liquid height H, m",)
fig.show()
```


```python
fig = go.Figure()
fig.add_trace(go.Scatter(
            x=hdf.df.index.total_seconds(), y= hdf.df['subsystem.tank_1.v1.H']*rho_water*A0[0],
            name='Tank 1 - numerical', mode='lines',line_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=result_analytical['t'], y=result_analytical['M_1'],
            name='Tank 1 - analytical', mode='markers', marker_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=hdf.df.index.total_seconds(), y=hdf.df['subsystem.tank_2.v1.H']*rho_water*A0[1],
            name='Tank 2 - numerical', mode='lines',line_color='rgba(152, 0, 0, .8)'))
fig.add_trace(go.Scatter(
            x=result_analytical['t'], y=result_analytical['M_2'],
            name='Tank 2 - analytical', mode='markers',marker_color='rgba(152, 0, 0, .8)'))
              
fig.update_layout(title='Liquid mass in tanks - analytical vs. numerical solution',
                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title="time, sec",
                 yaxis_title="liquid mass M, kg",)
fig.show()

```

## Running the simulations - Case 2: Time dependent Cv
We assign now a value to the parameter $\omega$, so that the valve has a characteristic parameter that is dependent from time, with a sinusoidal behaviour.

To obtain the solution of the system for Case 2 the procedure is analogous to the previous case.


```python
omega_2 = 1.5 # rad/s
    
# Solve the model and collect hystorical data
hdf_2     = SimpleHistoryDataFrame()
s2, sol_2 = t_1_item_model(H0, Cv_0,omega_2, hdf_2, t_start_sim, t_stop_sim)
```

### Plotting the results


```python
# Plot the results
# Nb: the package plotly is needed for the plotting 
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import plotly.graph_objs as go

fig = go.Figure()
fig.add_trace(go.Scatter(
            x=hdf_2.df.index.total_seconds(), y=hdf_2.df['subsystem.tank_1.v1.H'],
            name='Tank 1', mode='lines',line_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=hdf_2.df.index.total_seconds(), y=hdf_2.df['subsystem.tank_2.v1.H'],
            name='Tank 2', mode='lines',line_color='rgba(152, 0, 0, .8)'))
              
fig.update_layout(title='Liquid height in tanks',
                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title="time, sec",
                 yaxis_title="liquid height H, m",)
fig.show()
```


```python
fig = go.Figure()
fig.add_trace(go.Scatter(
            x=hdf_2.df.index.total_seconds(), y= hdf_2.df['subsystem.tank_1.v1.H']*rho_water*A0[0],
            name='Tank 1', mode='lines',line_color='rgba(102,51,255,1)'))
fig.add_trace(go.Scatter(
            x=hdf_2.df.index.total_seconds(), y=hdf_2.df['subsystem.tank_2.v1.H']*rho_water*A0[1],
            name='Tank 2', mode='lines',line_color='rgba(152, 0, 0, .8)'))
              
fig.update_layout(title='Liquid mass in tanks',
                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title="time, sec",
                 yaxis_title="liquid mass M, kg",)
fig.show()
```


```python
fig = go.Figure()
fig.add_trace(go.Scatter(
            x=hdf_2.df.index.total_seconds(), y= hdf_2.df['subsystem.valve_1.v1.Cv'],
            name='Cv', mode='lines',line_color='rgba(41, 241, 195, 1)'))           
fig.update_layout(title='Valve characteristc parameter',
                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title="time, sec",
                 yaxis_title="Cv, (kg/s)/Pa",)
fig.show()
```

5409925aa87cbac27d7ef3b2d6f1894a4ea5aebc
ipynb
Jupyter Notebook
examples/Two_Tanks_System/TwoTanks_System.ipynb
fossilfree/Numerous
121cd00456f263cab9def10cb99a21655f5832c9
BSD-3-Clause
examples/Two_Tanks_System/TwoTanks_System.ipynb
fossilfree/Numerous
121cd00456f263cab9def10cb99a21655f5832c9
BSD-3-Clause
examples/Two_Tanks_System/TwoTanks_System.ipynb
fossilfree/Numerous
121cd00456f263cab9def10cb99a21655f5832c9
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES