examples/Two_Tanks_System/TwoTanks_System.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Tutorial 2 - Two tanks system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Background information\n",
    "This tutorial shows how to use _numerous_ to create a system of multiple components.\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "The model implementation is developed by instantiating _items_ of different classes, and connecting them by means of _connectors_.\n",
    "\n",
    "The tutorial can be run only importing the extra packages pandas, numpy and plotly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## System governing equations\n",
    "\n",
    "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}$.\n",
    "![Screenshot.png](sketch_tanks.png)\n",
    "\n",
    "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.\n",
    "\n",
    "The mathematical model describing this systme is illusrated below.\n",
    "\n",
    "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:\n",
    "$$\\begin{equation}\n",
    "\\dot{\\mathrm{m}} = \\mathrm{C_v} \\sqrt{\\mathrm{\\Delta p}}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "\n",
    "Moreover, the conservation of mass in the two tanks is expressed by the equations (with the assumption of incompressible fluid):\n",
    "$$\\begin{equation}\n",
    "\\rho A_{0,1} \\mathrm{\\dfrac{dH_1}{dt}} = - \\dot{m}\n",
    "\\end{equation}\n",
    "$$\n",
    "$$\\begin{equation}\n",
    "\\rho A_{0,2}  \\mathrm{\\dfrac{dH_2}{dt}} = + \\dot{m}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "\n",
    "The pressure drop accross the valve can be expressed by the equation:\n",
    "$$\\begin{equation}\n",
    "\\Delta p = \\rho g H_1\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "By substituting this latter expression in the previous equations, the differential system of equations describing the system becomes:\n",
    "$$\\begin{equation}\n",
    "\\rho A_{0,1} \\mathrm{\\dfrac{dH_1}{dt}} = - \\mathrm{C_v} (t)\\sqrt{\\rho g \\mathrm{H_1}}\n",
    "\\end{equation}\n",
    "$$\n",
    "$$\\begin{equation}\n",
    "\\rho A_{0,2} \\mathrm{\\dfrac{dH_2}{dt}} = + \\mathrm{C_v}(t)\\sqrt{\\rho g \\mathrm{H_1}}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Analytical solution\n",
    "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_.\n",
    "\n",
    "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.\n",
    "\n",
    "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:\n",
    "$$\\begin{equation}\n",
    "\\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}\n",
    "\\end{equation}\n",
    "$$\n",
    "$$\\begin{equation}\n",
    "\\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\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "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:\n",
    "$$\\begin{equation}\n",
    "\\mathrm{M} = \\rho \\mathrm{H} \\mathrm{A_0}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "The analytical solution would thus be expressed as:\n",
    "$$\\begin{equation}\n",
    "\\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}\n",
    "\\end{equation}\n",
    "$$\n",
    "$$\\begin{equation}\n",
    "\\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\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Input data\n",
    "The input data used for the analytical solution and the simulation are reported below:\n",
    "* liquid density, $\\rho$ = 1000 kg/m$^3$\n",
    "* cross-sectional area tank 1, $A_{0,1}$ = 0.05 m$^2$\n",
    "* cross-sectional area tank 2, $A_{0,2}$ = 0.1 m$^2$\n",
    "* initial conditions:\n",
    "    * initial height tank 1  $H_{0,1}$ = 0.5 m\n",
    "    * initial height tank 2  $H_{0,2}$ = 0.0 m\n",
    "    \n",
    "\n",
    "The implementation of the model with a time-dependent $\\mathrm{C_v}$ was carried out by using the following formulation:\n",
    "* $\\mathrm{C_v} = \\mathrm{C_{v,0}} \\cdot \\big[sin{(\\omega \\cdot t)}+1\\big]  $\n",
    "and imposing:\n",
    "    * $\\mathrm{C_{v,0}} = 0.1 \\; \\; \\; (\\mathrm{kg}/\\mathrm{s})/\\mathrm{Pa}$\n",
    "    * $\\omega \\; \\; \\;  = 1.5 \\; \\; \\;\\mathrm{rad/s}$\n",
    "\n",
    "A constant value of the valve parameter $\\mathrm{C_v}$ was instead obtained by imposing $\\omega = 0$, so that:\n",
    "* $\\mathrm{C_v} = \\mathrm{C_{v,0}}$ \n",
    "\n",
    "at any time step, and the solution could be compared against the analytical formulation derived above.\n",
    "\n",
    "This is translated into the followinf python code:\n",
    "\n",
    "### Definition of input data - code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "## Input data\n",
    "# Define liquid and geometry for the case in analysis\n",
    "rho_water = 1000           # Water density [kg/m3]\n",
    "A0        = [0.05, 0.1]    # Tanks cross-sectional area [m2]\n",
    "Cv_0      = 0.1            # [(kg s)/Pa] Amplitude of valve characteristic parameter\n",
    "g         = 9.81           # [m/s2] gravitational accelleration\n",
    "\n",
    "## Initial conditions\n",
    "H0        = [0.5, 0]     # Initial condition - liquid height [m]\n",
    "# Estimate the initial value of total mass M [kg] in the two tanks\n",
    "M0        = rho_water * np.multiply(A0, H0)\n",
    "\n",
    "# Caclulate what is the time (analytical solution) to entirely fill the tank (simple case - analytical sol)\n",
    "t_filling = np.sqrt(A0[0]/g) * np.sqrt(M0[0]) *2* 1/Cv_0\n",
    "\n",
    "# Define start and stop time for the analysis and simulation\n",
    "t_start_sim = 0\n",
    "t_stop_sim  = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Analytical soltion: Python code\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "N_t             = 11                           # Number of time steps at which the analytical solution is evaluated\n",
    "time_range      = [t_start_sim ,t_stop_sim]    # Time span in which the solution is plotted\n",
    "\n",
    "def analytical_solution(N_t=N_t, time_range=time_range, g=9.81,  A0=[0.05, 0.1], Cv_0=0.1,\n",
    "                         M0=[25, 0], H0=[0.5, 0] , rho_water=1000):\n",
    "    time_vector       = np.linspace(start = time_range[0], stop = time_range[1], num = N_t)\n",
    "    \n",
    "    # Create a dictionary with time, evolution of liquid heights and mass in the two tanks for the analytical solution:\n",
    "    result_analytical = {'t': [], 'M_1': [], 'M_2': [], 'H_1': [], 'H_2': []}\n",
    "    \n",
    "    # Calculate what is the time after which the first tank is empty and the second tank is full\n",
    "    t_filling = np.sqrt(A0[0]/g) * np.sqrt(M0[0]) *2* 1/Cv_0\n",
    "    \n",
    "    for i in range(N_t):\n",
    "        if time_vector[i]  < t_filling:\n",
    "            M_1_analytical = (np.sqrt(M0[0]) - Cv_0 / 2 * np.sqrt(g / A0[0]) * time_vector[i] ) ** 2\n",
    "            M_2_analytical = np.sqrt(M0[1]) + Cv_0 * np.sqrt(g / A0[0]) * np.sqrt(\n",
    "                M0[0]) * time_vector[i]  - Cv_0 ** 2 * g / A0[0] / 4 * (time_vector[i]  ** 2)    \n",
    "        \n",
    "            H_1_analytical = (np.sqrt(H0[0]) - Cv_0 / 2 /A0[0] * np.sqrt(g/ rho_water) * time_vector[i] ) ** 2\n",
    "            H_2_analytical = np.sqrt(H0[1]) + Cv_0 /A0[1] * np.sqrt(g/ rho_water) * np.sqrt(\n",
    "                H0[0]) * time_vector[i]  - Cv_0 ** 2 * g / A0[0] / 4 /A0[1]/rho_water* (time_vector[i]  ** 2)\n",
    "        else:\n",
    "            M_2_analytical = M0[0] + M0[1]\n",
    "            M_1_analytical = 0\n",
    "            H_2_analytical = (M0[0] + M0[1])/rho_water/A0[1]\n",
    "            H_1_analytical = 0\n",
    "            \n",
    "        result_analytical['t'].append(time_vector[i])\n",
    "        result_analytical['M_1'].append(M_1_analytical)\n",
    "        result_analytical['M_2'].append(M_2_analytical)\n",
    "        result_analytical['H_1'].append(H_1_analytical)\n",
    "        result_analytical['H_2'].append(H_2_analytical)\n",
    "    return result_analytical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "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)\n",
    "\n",
    "data = {'Time, s': result_analytical['t'], 'H1, m':result_analytical['H_1'], 'H2, m':result_analytical['H_2'], \n",
    "        'M1, kg':result_analytical['M_1'],'M2, kg':result_analytical['M_2']}\n",
    "pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Implementation using _numerous_\n",
    "\n",
    "### Preliminary steps\n",
    "The first step for the implementation is to include all the relevant _numerous_ modules. \n",
    "For this tutorial we need:\n",
    "* _Item_ for defining item objects\n",
    "* _Model_ for defining model objects\n",
    "* _ConnectorTwoWay_ for defining a special connectors object\n",
    "* _Subsystems_ for defining the subsystem object\n",
    "* _Simulation_ for defining the simulation\n",
    "* _Equation_ and equation decorator for objects of equation class\n",
    "\n",
    "Moreover, we will need HistoryDataFrame to store results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# We include all the relevant modules from numerous:\n",
    "from numerous.engine.system import Item\n",
    "from numerous.engine.model import Model\n",
    "from numerous.engine.system import Subsystem\n",
    "from numerous.engine.simulation import Simulation\n",
    "from numerous.engine.system import ConnectorTwoWay\n",
    "from numerous.engine.variables import VariableType, VariableDescription, OverloadAction\n",
    "\n",
    "from numerous.multiphysics import Equation\n",
    "from numerous.multiphysics import EquationBase\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Define Tank Equation\n",
    "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.\n",
    "\n",
    "Note that the line of code\n",
    "``` python \n",
    "super().__init__(tag='for info only') \n",
    "```\n",
    "in `__init__` is a necessary line to create any class using _numerous_ in the current release.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "class Tank_Equation(EquationBase):\n",
    "    def __init__(self, tag=\"tank_equation\", H=0.5, rho=1000, A0=0.05, g=9.81):\n",
    "        super().__init__(tag='tank_equation')\n",
    "\n",
    "        self.add_state('H', H)          # [m] Liquid height in the tank\n",
    "        self.add_parameter('rho', rho)  # [kg/m3] Liquid density\n",
    "        self.add_parameter('A0', A0)    # [m2] Tank cross-sectional area\n",
    "        self.add_constant('g', g)       # [m/s2] Gravitational acceleration\n",
    "        self.add_parameter('mdot', 0)   # [kg/s] Mass flow rate \n",
    "        \n",
    "    @Equation()\n",
    "    def eval(self, scope):\n",
    "        # Differential equation for mass conservation in a general tank with a mass flow rate entering or leaving the tan\n",
    "        scope.H_dot = scope.mdot / scope.rho / scope.A0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Define Valve Equation\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%% Define the equation for the valve \n"
    }
   },
   "outputs": [],
   "source": [
    "# We define the equation (using EquationBase class) determining the mass flow rate across the valve\n",
    "class Valve_Equation(EquationBase):    \n",
    "    def __init__(self, Cv_0=0.1, rho=1000, g=9.81, omega=1):\n",
    "        super().__init__(tag='valve_equation')\n",
    "        \n",
    "        self.add_parameter('omega', omega)   # [rad/sec]   Angular frequency of valve characteristic parameter\n",
    "        self.add_parameter('Cv_0', Cv_0)     # [(kg/s)/Pa] Amplitude of valve characteristic parameter\n",
    "        self.add_parameter('mdot1', 0)       # [kg/s]      Mass flow rate in one side of the valve\n",
    "        self.add_parameter('mdot2', 0)       # [kg/s]      Mass flow rate in the other side of the valve\n",
    "        self.add_parameter('H1', 0)          # [m]         Liquid height in the tank 1 connected to the valve (top tank)\n",
    "        self.add_parameter('H2', 0)          # [m]         Liquid height in the tank 2 connected to the valve (bottom tank)\n",
    "        self.add_parameter('rho', rho)       # [kg/m3]     Liquid density\n",
    "        self.add_constant('_g', g)            # [m/s2]      Gravitational acceleration\n",
    "        self.add_parameter('Cv', Cv_0)       # [(kg/s)/Pa] Valve characteristic parameter\n",
    "        self.add_state('time', 0)\n",
    "        \n",
    "    @Equation()\n",
    "    def eval(self,scope):\n",
    "        scope.Cv      = scope.Cv_0 * (np.sin(scope.omega * scope.time)+1)           #[(kg/s)/Pa]\n",
    "        deltaP        = scope.rho * scope._g * (scope.H1)                                  #[Pa]\n",
    "        mdot          = np.sign(deltaP) * np.sqrt(np.absolute(deltaP)) * scope.Cv         #[kg/s]\n",
    "        # The valve will be associated with two mass flow rates (one leaving and one entering the component), \n",
    "        #which - for conservation of mass - have the same magnitude and opposite sign\n",
    "        scope.mdot1   = -mdot                                                             #[kg/s]\n",
    "        scope.mdot2   = mdot                                                              #[kg/s]\n",
    "        scope.time_dot = 1\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Define Tank as Item\n",
    "We define the `Tank` class as an `Item` class. We then create a namespace 'v1' to contain the variables for the `Tank_Equation`.\n",
    "The equation is associated to the namespace using the `add_equations` method, as shown in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "class Tank(Item):\n",
    "    def __init__(self, tag=\"tank\", H=0.5, rho=1000, A0=0.05, g=9.81):\n",
    "        super(Tank, self).__init__(tag)\n",
    "\n",
    "        v1 = self.create_namespace('v1')\n",
    "        v1.add_equations([Tank_Equation(H=H, rho=rho, A0=A0, g=g)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Define Valve as ConnectorTwoWay\n",
    "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.\n",
    "`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. \n",
    "\n",
    "In the code lines\n",
    "``` python \n",
    "super().__init__(tag, side1_name='side1', side2_name='side2')\n",
    "``` \n",
    "we have to specify the names of the two sides.\n",
    "\n",
    "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)\n",
    "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. \n",
    "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. \n",
    "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:\n",
    "    * 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.\n",
    "    * 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%% \n"
    }
   },
   "outputs": [],
   "source": [
    "# Define the valve as a connector item - connecting two tanks\n",
    "class Valve(ConnectorTwoWay):\n",
    "    def __init__(self, tag=\"valve\", Cv_0=0.1, rho=1000, g=9.81, omega=0):\n",
    "        super().__init__(tag, side1_name='side1', side2_name='side2')\n",
    "\n",
    "        #1 Create a namespace for mass flow rate equation and add the valve equation\n",
    "        v1 = self.create_namespace('v1') \n",
    "        v1.add_equations([Valve_Equation(Cv_0=Cv_0, rho=rho, g=g, omega=omega)])\n",
    "   \n",
    "        #2 Create variables H and mdot in side 1 adn 2\n",
    "           #(side 1 represents a connection with one tank, with related liquid height H_1)\n",
    "           #(side 1 represents a connection with the second tank, with related liquid height H_2)\n",
    "        self.side1.v1.create_variable(name='H')\n",
    "        self.side1.v1.create_variable(name='mdot')\n",
    "        self.side2.v1.create_variable(name='H')\n",
    "        self.side2.v1.create_variable(name='mdot')\n",
    "        \n",
    "        # Map variables between binding and internal variables for side 1 and 2\n",
    "        # This is needed to update the values of the variables in the binding according to the equtions of the items\n",
    "        v1.H1              = self.side1.v1.H\n",
    "        v1.H2              = self.side2.v1.H\n",
    "        \n",
    "        self.side1.v1.mdot = v1.mdot1\n",
    "        self.side2.v1.mdot = v1.mdot2\n",
    "       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "###  Create the sub-system of components \n",
    "After defining all the classes for the items that will consitute the system, we are ready for the system assembly.\n",
    "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)\n",
    "1. create the gravitational accelleration constant and assign a value to it\n",
    "2. create two instances of the class Tank called Tank_1 (top tank) and Tank_2 (bottom tank)\n",
    "3. create one instance of the class Valve called Valve_1\n",
    "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.\n",
    "5. register the instanciated items in the Two_Tanks class \n",
    "\n",
    "The inputs needed to the subsystem are:\n",
    "* H0, which is a vector containing the initial state of the system (initial liquid height of Tank_1 and Tank_2)\n",
    "* Cv_0 and omega, which are amplitude and angular frequency of the valve characteristic parameter\n",
    "We assume that geometry (A0) and liquid (rho) are given by the input data as fixed values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Define the subsystem composed by two tanks and one valve connecting them\n",
    "class Two_Tanks(Subsystem):\n",
    "    def __init__(self, tag, H0, Cv_0, omega):\n",
    "        super().__init__(tag)\n",
    "\n",
    "        #1. Gravitational acceleration\n",
    "        g = 9.81\n",
    "        \n",
    "        #2. Instances of Tank class\n",
    "        Tank_1 = Tank('tank_1', H=H0[0], rho=rho_water, A0=A0[0], g=g)\n",
    "        Tank_2 = Tank('tank_2', H=H0[1], rho=rho_water, A0=A0[1], g=g)\n",
    "        \n",
    "        #3. Valve_1 is one instance of valve class\n",
    "        Valve_1 = Valve('valve_1', Cv_0=Cv_0, rho=rho_water, g=g, omega=omega)\n",
    "        \n",
    "        #4. Binding\n",
    "        Valve_1.bind(side1=Tank_1, side2=Tank_2)\n",
    "\n",
    "        #5. Register all the instanciated items in the sub-system\n",
    "        self.register_items([Tank_1, Tank_2, Valve_1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "###  Create the system model and simulation\n",
    "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:\n",
    "* H0, Cv_0, omega, which represent initial conditions of the tanks, and valve characteristic\n",
    "* hdf, which is the historian where the simulation results will be stored\n",
    "* t_start_sim and t_stop_sim determing the time span in which the system will be simulated\n",
    "\n",
    "The steps are the following:\n",
    "1. First a model object m1 is instantiated based on the `Two_Tanks` subsystem with given inputs\n",
    "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\n",
    "3. The `simulation.solve` method is called, and the solution is returned as output of the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def t_1_item_model(H0, Cv_0, omega, t_start_sim, t_stop_sim):\n",
    "    \n",
    "    # 1. Instantiate the model\n",
    "    m1 = Model(Two_Tanks('subsystem', H0, Cv_0, omega))\n",
    "    # 2. Setting up the simulation:\n",
    "    s1  = Simulation(m1, t_start=t_start_sim, t_stop=t_stop_sim, num=5000)\n",
    "    # 3. Call the solve method\n",
    "    sol = s1.solve()\n",
    "    return s1, sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Running the simulations - Case 1: Constant Cv \n",
    "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).\n",
    "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.\n",
    "The progress bar is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "omega = 0\n",
    "\n",
    "# Collect the historical data at each time step\n",
    "s1, sol = t_1_item_model(H0, Cv_0,omega, t_start_sim, t_stop_sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Plotting: comparison vs. analytical solution\n",
    "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.\n",
    "\n",
    "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`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "hdf = s1.model.historian_df\n",
    "x_value =  [x / 500 for x in list(hdf.index)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the comparison\n",
    "# Nb: the package plotly is needed for the plotting \n",
    "from plotly import __version__\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value, y=hdf['subsystem.tank_1.v1.H'],\n",
    "            name='Tank 1 - numerical', mode='lines',line_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=result_analytical['t'], y=result_analytical['H_1'],\n",
    "            name='Tank 1 - analytical', mode='markers', marker_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value, y=hdf['subsystem.tank_2.v1.H'],\n",
    "            name='Tank 2 - numerical', mode='lines',line_color='rgba(152, 0, 0, .8)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=result_analytical['t'], y=result_analytical['H_2'],\n",
    "            name='Tank 2 - analytical',  mode='markers',marker_color='rgba(152, 0, 0, .8)'))\n",
    "              \n",
    "fig.update_layout(title='Liquid height in tanks - analytical vs. numerical solution',\n",
    "                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title=\"time, sec\",\n",
    "                 yaxis_title=\"liquid height H, m\",)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value, y= hdf['subsystem.tank_1.v1.H']*rho_water*A0[0],\n",
    "            name='Tank 1 - numerical', mode='lines',line_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=result_analytical['t'], y=result_analytical['M_1'],\n",
    "            name='Tank 1 - analytical', mode='markers', marker_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value, y=hdf['subsystem.tank_2.v1.H']*rho_water*A0[1],\n",
    "            name='Tank 2 - numerical', mode='lines',line_color='rgba(152, 0, 0, .8)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=result_analytical['t'], y=result_analytical['M_2'],\n",
    "            name='Tank 2 - analytical', mode='markers',marker_color='rgba(152, 0, 0, .8)'))\n",
    "              \n",
    "fig.update_layout(title='Liquid mass in tanks - analytical vs. numerical solution',\n",
    "                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title=\"time, sec\",\n",
    "                 yaxis_title=\"liquid mass M, kg\",)\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Running the simulations - Case 2: Time dependent Cv\n",
    "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.\n",
    "\n",
    "To obtain the solution of the system for Case 2 the procedure is analogous to the previous case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "omega_2 = 1.5 # rad/s\n",
    "    \n",
    "# Solve the model and collect hystorical data\n",
    "s2, sol_2 = t_1_item_model(H0, Cv_0,omega_2, t_start_sim, t_stop_sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Plotting the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "hdf_2  =  s2.model.historian_df\n",
    "x_value_2 =  [x / 500 for x in list(hdf_2.index)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "# Nb: the package plotly is needed for the plotting \n",
    "from plotly import __version__\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value_2, y=hdf_2['subsystem.tank_1.v1.H'],\n",
    "            name='Tank 1', mode='lines',line_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value_2, y=hdf_2['subsystem.tank_2.v1.H'],\n",
    "            name='Tank 2', mode='lines',line_color='rgba(152, 0, 0, .8)'))\n",
    "              \n",
    "fig.update_layout(title='Liquid height in tanks',\n",
    "                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title=\"time, sec\",\n",
    "                 yaxis_title=\"liquid height H, m\",)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value_2, y= hdf_2['subsystem.tank_1.v1.H']*rho_water*A0[0],\n",
    "            name='Tank 1', mode='lines',line_color='rgba(102,51,255,1)'))\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value_2, y=hdf_2['subsystem.tank_2.v1.H']*rho_water*A0[1],\n",
    "            name='Tank 2', mode='lines',line_color='rgba(152, 0, 0, .8)'))\n",
    "              \n",
    "fig.update_layout(title='Liquid mass in tanks',\n",
    "                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title=\"time, sec\",\n",
    "                 yaxis_title=\"liquid mass M, kg\",)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "            x=x_value_2, y= hdf_2['subsystem.valve_1.v1.Cv'],\n",
    "            name='Cv', mode='lines',line_color='rgba(41, 241, 195, 1)'))           \n",
    "fig.update_layout(title='Valve characteristc parameter',\n",
    "                  yaxis_zeroline=False, xaxis_zeroline=False,xaxis_title=\"time, sec\",\n",
    "                 yaxis_title=\"Cv, (kg/s)/Pa\",)\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": [
     "##%% md                       \n",
     "# Two tank system \n"
    ]
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}