first try at fixing convergence issues at turbine

Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -21,7 +21,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -51,7 +51,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -71,14 +71,14 @@
     "\n",
     "# for while loop\n",
     "total_min_level             = 0.01      # m\n",
-    "total_max_time              = 100     # s\n",
+    "total_max_time              = 1000     # s\n",
     "\n",
     "nt = int(total_max_time//simulation_timestep)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -119,7 +119,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -149,7 +149,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
@@ -158,7 +158,7 @@
        "10.1"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 12,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -170,7 +170,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.8.13 ('DT_Slot_3')",
+   "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
    "language": "python",
    "name": "python3"
   },
@@ -189,7 +189,7 @@
   "orig_nbformat": 4,
   "vscode": {
    "interpreter": {
-    "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
+    "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
    }
   }
  },

Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -22,7 +22,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -36,7 +36,7 @@
     "D       = 0.9                   # pipe diameter [m]\n",
     "h_res   = 10.                   # water level in upstream reservoir [m]\n",
     "n       = 50                  # number of pipe segments in discretization\n",
-    "nt      = 1000                    # number of time steps after initial conditions\n",
+    "nt      = 5000                    # number of time steps after initial conditions\n",
     "f_D     = 0.01                  # Darcy friction factor\n",
     "c       = 400.                  # propagation velocity of the pressure wave [m/s]\n",
     "h_pipe  = 105.                  # hydraulic head without reservoir [m] \n",
@@ -57,7 +57,7 @@
     "# define constants reservoir\n",
     "conversion_pressure_unit    = 'mWS'\n",
     "\n",
-    "area_base                   = 75.                # m²\n",
+    "area_base                   = 75.               # m²\n",
     "area_pipe                   = (D/2)**2*np.pi    # m²\n",
     "critical_level_low          = 0.                # m\n",
     "critical_level_high         = 100.              # m\n",
@@ -69,7 +69,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -84,7 +84,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
@@ -110,7 +110,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -141,7 +141,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -164,7 +164,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -214,7 +214,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -250,7 +250,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.8.13 ('DT_Slot_3')",
+   "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
    "language": "python",
    "name": "python3"
   },
@@ -269,7 +269,7 @@
   "orig_nbformat": 4,
   "vscode": {
    "interpreter": {
-    "hash": "4a28055eb8a3160fa4c7e4fca69770c4e0a1add985300856aa3fcf4ce32a2c48"
+    "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
    }
   }
  },

Turbinen/convergence_turbine.py

@@ -0,0 +1,169 @@
+from time import time
+import numpy as np
+#importing pressure conversion function
+import sys
+import os
+current = os.path.dirname(os.path.realpath(__file__))
+parent = os.path.dirname(current)
+sys.path.append(parent)
+from functions.pressure_conversion import pressure_conversion
+
+class Francis_Turbine_test:
+ # units 
+    # make sure that units and print units are the same
+    # units are used to label graphs and print units are used to have a bearable format when using pythons print()
+    density_unit        = r'$\mathrm{kg}/\mathrm{m}^3$'
+    flux_unit           = r'$\mathrm{m}^3/\mathrm{s}$'
+    LA_unit             = '%'
+    pressure_unit       = 'Pa'
+    time_unit           = 's'
+    velocity_unit       = r'$\mathrm{m}/\mathrm{s}$'
+    volume_unit         = r'$\mathrm{m}^3$'
+
+    density_unit_print      = 'kg/m³'
+    flux_unit_print         = 'm³/s'
+    LA_unit_print           = '%'
+    pressure_unit_print     = 'mWS' 
+    time_unit_print         = 's'
+    velocity_unit_print     = 'm/s'
+    volume_unit_print       = 'm³'
+
+    g = 9.81    # m/s²  gravitational acceleration
+
+
+ # init
+    def __init__(self, Q_nenn,p_nenn,t_closing=-1.,timestep=-1.):
+        self.Q_n    = Q_nenn                                    # nominal flux 
+        self.p_n    = p_nenn                                    # nominal pressure
+        self.LA_n   = 1. # 100%                                 # nominal Leitapparatöffnung
+        h           = pressure_conversion(p_nenn,'Pa','MWs')    # nominal pressure in terms of hydraulic head
+        self.A      = Q_nenn/(np.sqrt(2*self.g*h)*0.98)         # Ersatzfläche
+
+        self.dt     = timestep                                  # simulation timestep
+        self.t_c = t_closing                                    # closing time
+        self.d_LA_max_dt = 1/t_closing                          # maximal change of LA per second
+
+        # initialize for get_info() - parameters will be converted to display -1 if not overwritten
+        self.p  = pressure_conversion(-1,self.pressure_unit_print,self.pressure_unit)
+        self.Q  = -1.
+        self.LA = -0.01 
+
+
+# setter
+    def set_LA(self,LA,display_warning=True):
+        # set Leitapparatöffnung
+        self.LA = LA
+        # warn user, that the .set_LA() method should not be used ot set LA manually
+        if display_warning == True:
+            print('Consider using the .update_LA() method instead of setting LA manually')
+
+    def set_timestep(self,timestep,display_warning=True):
+        # set Leitapparatöffnung
+        self.dt = time
+        # warn user, that the .set_LA() method should not be used ot set LA manually
+        if display_warning == True:
+            print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!')
+
+    def set_pressure(self,pressure):
+        # set pressure in front of the turbine
+        self.p = pressure
+
+#getter
+    def get_current_Q(self):
+        # return the flux through the turbine, based on the current pressure in front
+            #  of the turbine and the Leitapparatöffnung
+        if self.p < 0:
+            self.Q = 0
+        else:
+            self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(self.p/self.p_n)
+        return self.Q
+
+    def get_current_pressure(self):
+        return self.p
+
+    def get_current_LA(self):
+        return self.LA
+
+    def get_info(self, full = False):
+        new_line = '\n'
+        p   = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print)
+        p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print)
+
+        
+        if full == True:
+            # :<10 pads the self.value to be 10 characters wide
+            print_str = (f"Turbine has the following attributes: {new_line}" 
+                f"----------------------------- {new_line}"
+                f"Type                  =       Francis {new_line}"
+                f"Nominal flux          =       {self.Q_n:<10} {self.flux_unit_print} {new_line}"
+                f"Nominal pressure      =       {round(p_n,3):<10} {self.pressure_unit_print}{new_line}"
+                f"Nominal LA            =       {self.LA_n*100:<10} {self.LA_unit_print} {new_line}"
+                f"Closing time          =       {self.t_c:<10} {self.time_unit_print} {new_line}"
+                f"Current flux          =       {self.Q:<10} {self.flux_unit_print} {new_line}" 
+                f"Current pipe pressure =       {round(p,3):<10} {self.pressure_unit_print} {new_line}"
+                f"Current LA            =       {self.LA*100:<10} {self.LA_unit_print} {new_line}"
+                f"Simulation timestep   =       {self.dt:<10} {self.time_unit_print} {new_line}"
+                f"----------------------------- {new_line}")
+        else:
+            # :<10 pads the self.value to be 10 characters wide
+            print_str = (f"The current attributes are: {new_line}" 
+                f"----------------------------- {new_line}"
+                f"Current flux          =       {self.Q:<10} {self.flux_unit_print} {new_line}" 
+                f"Current pipe pressure =       {round(p,3):<10} {self.pressure_unit_print} {new_line}"
+                f"Current LA            =       {self.LA*100:<10} {self.LA_unit_print} {new_line}"
+                f"----------------------------- {new_line}")
+
+        print(print_str)
+
+# methods
+    def update_LA(self,LA_soll):
+        # update the Leitappartöffnung and consider the restrictions of the closing time of the turbine
+        LA_diff = self.LA-LA_soll                   # calculate the difference to the target LA
+        LA_diff_max = self.d_LA_max_dt*self.dt     # calculate the maximum change in LA based on the given timestep  
+        LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max]))    # calulate the correct change in LA
+        
+        LA_new = self.LA-LA_diff
+        if LA_new < 0.:
+            LA_new = 0.
+        elif LA_new > 1.:
+            LA_new = 1.
+        self.set_LA(LA_new,display_warning=False)
+
+    def set_steady_state(self,ss_flux,ss_pressure):
+        # calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the
+            # turbine at the steady state LA
+        ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure)
+        if ss_LA < 0 or ss_LA > 1:
+            raise Exception('LA out of range [0;1]')
+        self.set_LA(ss_LA,display_warning=False)
+
+    def converge(self,area_pipe,pressure_s2l_node,velocity_s2l_node,alpha,f_D,dt):
+        eps = 1e-9
+        error = 1.
+        i = 0
+        p = pressure_s2l_node
+        v = velocity_s2l_node
+        rho = 1000
+        g = self.g
+        c = 400
+        D = area_pipe
+        p_old = self.get_current_pressure()
+        Q_old = self.get_current_Q()
+        v_old = Q_old/area_pipe
+
+
+        while error > eps:
+            self.set_pressure(p_old)
+            Q_new = self.get_current_Q()
+            v_new = Q_new/area_pipe
+            p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v
+
+            error = abs(Q_old-Q_new)
+            Q_old = Q_new.copy()
+            p_old = p_new.copy()
+            v_old = v_new.copy()
+            i = i+1
+            if i == 1e6:
+                print('did not converge')
+                break
+        self.Q = Q_new

Turbinen/turbine_convergence_test.ipynb

@@ -0,0 +1,286 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from convergence_turbine import Francis_Turbine_test\n",
+    "\n",
+    "#importing pressure conversion function\n",
+    "import sys\n",
+    "import os\n",
+    "current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
+    "parent = os.path.dirname(current)\n",
+    "sys.path.append(parent)\n",
+    "from functions.pressure_conversion import pressure_conversion\n",
+    "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class\n",
+    "from Druckrohrleitung.Druckrohrleitung_class_file import Druckrohrleitung_class"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib qt5\n",
+    "\n",
+    "#Turbine\n",
+    "Q_nenn          = 0.85  # m³/s\n",
+    "p_nenn          = pressure_conversion(10.6,'bar','Pa')\n",
+    "closing_time    = 5    #s\n",
+    "\n",
+    "#define constants pipe\n",
+    "\n",
+    "g       = 9.81                  # gravitational acceleration [m/s²]\n",
+    "rho     = 1000.                 # density of water [kg/m³]\n",
+    "\n",
+    "L       = 1000.                 # length of pipeline [m]\n",
+    "D       = 0.9                   # pipe diameter [m]\n",
+    "h_res   = 10.                   # water level in upstream reservoir [m]\n",
+    "n       = 50                  # number of pipe segments in discretization\n",
+    "nt      = 10000                    # number of time steps after initial conditions\n",
+    "f_D     = 0.01                  # Darcy friction factor\n",
+    "c       = 400.                  # propagation velocity of the pressure wave [m/s]\n",
+    "h_pipe  = 105.                  # hydraulic head without reservoir [m] \n",
+    "alpha   = np.arcsin(h_pipe/L)   # Höhenwinkel der Druckrohrleitung \n",
+    "\n",
+    "\n",
+    "# preparing the discretization and initial conditions\n",
+    "initial_flux    = 0.8                                    # m³/s\n",
+    "initial_level   = h_res                                   # m\n",
+    "dx              = L/n                                   # length of each pipe segment\n",
+    "dt              = dx/c                                  # timestep according to method of characterisitics\n",
+    "nn              = n+1                                   # number of nodes\n",
+    "pl_vec          = np.arange(0,nn,1)*dx                  # pl = pipe-length. position of the nodes on the pipeline\n",
+    "t_vec           = np.arange(0,nt,1)*dt                  # time vector\n",
+    "h_vec           = np.arange(0,nn,1)*h_pipe/n            # hydraulic head of pipeline at each node\n",
+    "\n",
+    "\n",
+    "# define constants reservoir\n",
+    "conversion_pressure_unit    = 'mWS'\n",
+    "\n",
+    "area_base                   = 75.               # m²\n",
+    "area_pipe                   = (D/2)**2*np.pi    # m²\n",
+    "critical_level_low          = 0.                # m\n",
+    "critical_level_high         = 100.              # m\n",
+    "\n",
+    "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
+    "nt_eRK4                     = 1              # number of simulation steps of reservoir in between timesteps of pipeline              \n",
+    "simulation_timestep         = dt/nt_eRK4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "V = Ausgleichsbecken_class(area_base,area_pipe,critical_level_low,critical_level_high,simulation_timestep)\n",
+    "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
+    "\n",
+    "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
+    "pipe.set_pressure_propagation_velocity(c)\n",
+    "pipe.set_number_of_timesteps(nt)\n",
+    "pipe.set_steady_state(initial_flux,initial_level,area_base,pl_vec,h_vec)\n",
+    "\n",
+    "\n",
+    "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
+    "T1 = Francis_Turbine_test(Q_nenn,p_nenn,closing_time,timestep=dt)\n",
+    "T1.set_steady_state(initial_flux,initial_pressure_turbine)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# initialization for timeloop\n",
+    "\n",
+    "level_vec       = np.zeros_like(t_vec)\n",
+    "level_vec[0]    = V.get_current_level()\n",
+    "\n",
+    "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
+    "v_old = pipe.get_current_velocity_distribution()\n",
+    "p_old = pipe.get_current_pressure_distribution()\n",
+    "\n",
+    "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
+    "    # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
+    "        # through the time evolution of the reservoir respectively  \n",
+    "    # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
+    "v_boundary_res  = np.zeros_like(t_vec)\n",
+    "v_boundary_tur  = np.zeros_like(t_vec)\n",
+    "p_boundary_res  = np.zeros_like(t_vec)\n",
+    "p_boundary_tur  = np.zeros_like(t_vec)\n",
+    "\n",
+    "# set the boundary conditions for the first timestep\n",
+    "v_boundary_res[0]   = v_old[0]\n",
+    "v_boundary_tur[0]   = v_old[-1] \n",
+    "p_boundary_res[0]   = p_old[0]\n",
+    "p_boundary_tur[0]   = p_old[-1]\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1,axs1 = plt.subplots(2,1)\n",
+    "axs1[0].set_title('Pressure distribution in pipeline')\n",
+    "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
+    "axs1[0].set_ylabel(r'$p$ [mWS]')\n",
+    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa',conversion_pressure_unit),marker='.')\n",
+    "axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_old,'Pa',conversion_pressure_unit))])\n",
+    "\n",
+    "axs1[1].set_title('Velocity distribution in pipeline')\n",
+    "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
+    "axs1[1].set_ylabel(r'$v$ [m/s]')\n",
+    "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
+    "# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
+    "\n",
+    "fig1.tight_layout()\n",
+    "plt.pause(1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "for it_pipe in range(1,nt):\n",
+    "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
+    "    # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
+    "    V.set_pressure(p_old[0])\n",
+    "    V.set_outflux(v_old[0]*area_pipe)\n",
+    "    # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
+    "    for it_res in range(nt_eRK4):\n",
+    "        V.timestep_reservoir_evolution()                                                             \n",
+    "    level_vec[it_pipe] = V.get_current_level() \n",
+    "\n",
+    "    \n",
+    "    # set boundary conditions for the next timestep of the characteristic method\n",
+    "    p_boundary_res[it_pipe] = V.get_current_pressure()\n",
+    "    T1.set_pressure(p_old[-1])\n",
+    "    T1.converge(area_pipe,p_old[-2],v_old[-2],alpha,f_D,dt)\n",
+    "    v_boundary_tur[it_pipe] = T1.get_current_Q()/area_pipe\n",
+    "\n",
+    "    # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
+    "    pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
+    "    pipe.v[0] = (pipe.v[0]+V.get_current_outflux()/area_pipe)/2\n",
+    "    p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
+    "    v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
+    "\n",
+    "    # perform the next timestep via the characteristic method\n",
+    "    pipe.timestep_characteristic_method()\n",
+    "\n",
+    "    # prepare for next loop\n",
+    "    p_old = pipe.get_current_pressure_distribution()\n",
+    "    v_old = pipe.get_current_velocity_distribution()\n",
+    "\n",
+    "    # plot some stuff\n",
+    "        # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
+    "    lo_00.remove()\n",
+    "    lo_01.remove()\n",
+    "    # lo_02.remove()\n",
+    "        # plot new pressure and velocity distribution in the pipeline\n",
+    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa', conversion_pressure_unit),marker='.',c='blue')\n",
+    "    lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
+    "    \n",
+    "    fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(round(t_vec[-1],2)))\n",
+    "    fig1.canvas.draw()\n",
+    "    fig1.tight_layout()\n",
+    "    plt.pause(0.000001)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig2,axs2 = plt.subplots(2,2)\n",
+    "axs2[0,0].set_title('Pressure Reservoir')\n",
+    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))\n",
+    "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[0,0].set_ylabel(r'$p$ [mWS]')\n",
+    "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))])\n",
+    "\n",
+    "axs2[0,1].set_title('Velocity Reservoir')\n",
+    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
+    "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
+    "\n",
+    "axs2[1,0].set_title('Pressure Turbine')\n",
+    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))\n",
+    "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[1,0].set_ylabel(r'$p$ [mWS]')\n",
+    "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))])\n",
+    "\n",
+    "axs2[1,1].set_title('Velocity Turbine')\n",
+    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
+    "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "axs2[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n",
+    "\n",
+    "fig2.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x1783803c0d0>]"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "plt.plot(level_vec)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
+   "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.8.13"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Untertweng_mit_Pegelregler.ipynb

@@ -42,7 +42,7 @@
     "\n",
     "\n",
     "# pipeline\n",
-    "L               = 535.+478.                     # length of pipeline [m]\n",
+    "L               = (535.+478.)                     # length of pipeline [m]\n",
     "D               = 0.9                           # pipe diameter [m]\n",
     "A_pipe          = D**2/4*np.pi                  # pipeline area\n",
     "h_pipe          = 105                           # hydraulic head without reservoir [m] \n",
@@ -51,7 +51,7 @@
     "f_D             = 0.014                           # Darcy friction factor\n",
     "c               = 500.                          # propagation velocity of the pressure wave [m/s]\n",
     "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
-    "nt              = 4500                          # number of time steps after initial conditions\n",
+    "nt              = 9000                          # number of time steps after initial conditions\n",
     "\n",
     "# derivatives of the pipeline constants\n",
     "dx              = L/n                           # length of each pipe segment\n",
@@ -184,6 +184,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "error_vec = np.zeros_like(t_vec)\n",
     "# loop through time steps of the pipeline\n",
     "for it_pipe in range(1,pipe.nt+1):\n",
     "\n",
@@ -218,6 +219,7 @@
     "    p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
     "    v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
     "\n",
+    "    error_vec[it_pipe] = abs(v_boundary_res[it_pipe]-V.get_current_outflux()/A_pipe)\n",
     "\n",
     "    # perform the next timestep via the characteristic method\n",
     "    pipe.timestep_characteristic_method()\n",
@@ -288,6 +290,26 @@
     "fig2.tight_layout()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x1ac81d70af0>]"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "plt.semilogy(t_vec,error_vec)"
+   ]
   }
  ],
  "metadata": {