From 7394d6c96421631d5d93c0e7fa1c283075dbab7d Mon Sep 17 00:00:00 2001 From: Brantegger Georg Date: Tue, 9 Aug 2022 13:51:57 +0200 Subject: [PATCH] added "vectorized" version of method of characteristics --- .../Druckrohrleitung_class_file.py | 34 ++++++- .../Druckrohrleitung_test_steady_state.ipynb | 98 ++----------------- 2 files changed, 43 insertions(+), 89 deletions(-) diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py index c6ce967..d5f3737 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -240,7 +240,7 @@ class Druckrohrleitung_class: self.v[i] = 0.5*(self.v_old[i+1]+self.v_old[i-1])-0.5/(rho*c)*(self.p_old[i+1]-self.p_old[i-1]) \ +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]+abs(self.v_old[i-1])*self.v_old[i-1]) - self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1]) - 0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \ + self.p[i] = 0.5*(self.p_old[i+1]+self.p_old[i-1])-0.5*rho*c*(self.v_old[i+1]-self.v_old[i-1]) \ +f_D*rho*c*dt/(4*D)*(abs(self.v_old[i+1])*self.v_old[i+1]-abs(self.v_old[i-1])*self.v_old[i-1]) # update overall min and max values for pressure and velocity per node @@ -254,3 +254,35 @@ class Druckrohrleitung_class: # else one can overwrite data by accidient and change two variables at once without noticing self.p_old = self.p.copy() self.v_old = self.v.copy() + + def timestep_characteristic_method_vectorized(self): + # use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones + # they are set with the .set_boundary_conditions_next_timestep() method beforehand + + # constants for cleaner formula + rho = self.density # density of liquid + c = self.c # pressure propagation velocity + f_D = self.f_D # Darcy friction coefficient + dt = self.dt # timestep + D = self.dia # pipeline diameter + g = self.g # graviational acceleration + alpha = self.angle # pipeline angle + + # Vectorized loop + self.v[1:-1] = 0.5*(self.v_old[2:]+self.v_old[:-2])-0.5/(rho*c)*(self.p_old[2:]-self.p_old[:-2]) \ + +dt*g*np.sin(alpha)-f_D*dt/(4*D)*(np.abs(self.v_old[2:])*self.v_old[2:]+np.abs(self.v_old[:-2])*self.v_old[:-2]) + + self.p[1:-1] = 0.5*(self.p_old[2:]+self.p_old[:-2])-0.5*rho*c*(self.v_old[2:]-self.v_old[:-2]) \ + +f_D*rho*c*dt/(4*D)*(np.abs(self.v_old[2:])*self.v_old[2:]-np.abs(self.v_old[:-2])*self.v_old[:-2]) + + # update overall min and max values for pressure and velocity per node + self.p_min = np.minimum(self.p_min,self.p) + self.p_max = np.maximum(self.p_max,self.p) + self.v_min = np.minimum(self.v_min,self.v) + self.v_max = np.maximum(self.v_max,self.v) + + # prepare for next call + # use .copy() to write data to another memory location and avoid the usual python reference pointer + # else one can overwrite data by accidient and change two variables at once without noticing + self.p_old = self.p.copy() + self.v_old = self.v.copy() diff --git a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb index 5e73d33..86d90d3 100644 --- a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb +++ b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -80,48 +80,9 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The pipeline has the following attributes: \n", - "----------------------------- \n", - "Length = 1013.0 m \n", - "Diameter = 0.9 m \n", - "Hydraulic head = 105.0 m \n", - "Number of segments = 50 \n", - "Number of nodes = 51 \n", - "Length per segments = 20.26 m \n", - "Pipeline angle = 0.104 rad \n", - "Pipeline angle = 5.95° \n", - "Darcy friction factor = 0.014 \n", - "Density of liquid = 1000.0 kg/m³ \n", - "Pressure wave vel. = 500.0 m/s \n", - "Simulation timestep = 0.04052 s \n", - "----------------------------- \n", - "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n", - "The pipeline has the following attributes: \n", - "----------------------------- \n", - "Length = 1013.0 m \n", - "Diameter = 0.9 m \n", - "Hydraulic head = 105.0 m \n", - "Number of segments = 50 \n", - "Number of nodes = 51 \n", - "Length per segments = 20.26 m \n", - "Pipeline angle = 0.104 rad \n", - "Pipeline angle = 5.95° \n", - "Darcy friction factor = 0.014 \n", - "Density of liquid = 1000.0 kg/m³ \n", - "Pressure wave vel. = 500.0 m/s \n", - "Simulation timestep = 0.04052 s \n", - "----------------------------- \n", - "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n" - ] - } - ], + "outputs": [], "source": [ "# create objects\n", "\n", @@ -137,7 +98,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -168,7 +129,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -193,48 +154,9 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The cuboid reservoir has the following attributes: \n", - "----------------------------- \n", - "Base area = 74.0 m² \n", - "Outflux area = 0.636 m² \n", - "Current level = 8.0 m\n", - "Critical level low = 0.0 m \n", - "Critical level high = inf m \n", - "Volume in reservoir = 592.0 m³ \n", - "Current influx = 0.773 m³/s \n", - "Current outflux = 0.773 m³/s \n", - "Current outflux vel = 1.215 m/s \n", - "Current pipe pressure = 7.854 mWS \n", - "Simulation timestep = 0.001013 s \n", - "Density of liquid = 1000.0 kg/m³ \n", - "----------------------------- \n", - "\n", - "The pipeline has the following attributes: \n", - "----------------------------- \n", - "Length = 1013.0 m \n", - "Diameter = 0.9 m \n", - "Hydraulic head = 105.0 m \n", - "Number of segments = 50 \n", - "Number of nodes = 51 \n", - "Length per segments = 20.26 m \n", - "Pipeline angle = 0.104 rad \n", - "Pipeline angle = 5.95° \n", - "Darcy friction factor = 0.014 \n", - "Density of liquid = 1000.0 kg/m³ \n", - "Pressure wave vel. = 500.0 m/s \n", - "Simulation timestep = 0.04052 s \n", - "----------------------------- \n", - "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n" - ] - } - ], + "outputs": [], "source": [ "for it_pipe in range(1,nt+1):\n", "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n", @@ -257,7 +179,7 @@ " 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", + " pipe.timestep_characteristic_method_vectorized()\n", "\n", " # prepare for next loop\n", " p_old = pipe.get_current_pressure_distribution()\n", @@ -283,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [