From bd94aa8d6559a2013ae88fd5fa69feb4b3427372 Mon Sep 17 00:00:00 2001 From: Aloma Blanch Date: Wed, 8 Jan 2020 12:40:10 +0100 Subject: [PATCH] original files --- MCS.py | 93 ++++++++++++++++++++++++++++++++++++++ post_Process.py | 118 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 211 insertions(+) create mode 100644 MCS.py create mode 100644 post_Process.py diff --git a/MCS.py b/MCS.py new file mode 100644 index 0000000..5a7e409 --- /dev/null +++ b/MCS.py @@ -0,0 +1,93 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Dec 11 07:27:15 2019 + +@author: 6035blancha +""" + +import numpy as np +import matplotlib.pyplot as plt + +T_cyc = 0.477 +#cyc = int(input("Number of cardiac cycles simulated on all meshes? ")) +cyc = 5 + + +path = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\" +file = "-procs_case\\PHistRCR.dat" + + + + +# Mesh convergence study pressure evaluation +# Load 1st Mesh results +N_ts_m1 = 4760 +nCPUs_m1 = "48" +folder_project_m1 = "MCS_m01" +#N_ts_m1 = int(input("Number of Time Steps used for this mesh? ")) +#nCPUs_m1 = str(input("Number of processors used for this mesh: ")) +#folder_project_m1 = str(input("Type folder where project data is located: ")) +pressure_filename_m1 = path + folder_project_m1 + "\\" + nCPUs_m1 + file +pressure_m1 = np.loadtxt(pressure_filename_m1, skiprows=1,) +time_m1 = np.linspace(0,T_cyc*cyc,pressure_m1.shape[0]) + +# Load 2nd Mesh results +N_ts_m2 = 7930 +nCPUs_m2 = "48" +folder_project_m2 = "MCS_m02" +#N_ts_m2 = int(input("Number of Time Steps used for this mesh? ")) +#nCPUs_m2 = str(input("Number of processors used for this mesh: ")) +#folder_project_m2 = str(input("Type folder where project data is located: ")) +pressure_filename_m2 = path + folder_project_m2 + "\\" + nCPUs_m2 + file +pressure_m2 = np.loadtxt(pressure_filename_m2, skiprows=1,) +time_m2 = np.linspace(0,T_cyc*cyc,pressure_m2.shape[0]) + +# Load 3rd Mesh results +N_ts_m3 = 10200 +nCPUs_m3 = "48" +folder_project_m3 = "MCS_m03" +#N_ts_m3 = int(input("Number of Time Steps used for this mesh? ")) +#nCPUs_m3 = str(input("Number of processors used for this mesh: ")) +#folder_project_m3 = str(input("Type folder where project data is located: ")) +pressure_filename_m3 = path + folder_project_m3 + "\\" + nCPUs_m3 + file +pressure_m3 = np.loadtxt(pressure_filename_m3, skiprows=1,) +time_m3 = np.linspace(0,T_cyc*cyc,pressure_m3.shape[0]) + +# Load 4rd Mesh results +N_ts_m4 = 14280 +nCPUs_m4 = "96" +folder_project_m4 = "MCS_m06" +#N_ts_m4 = int(input("Number of Time Steps used for this mesh? ")) +#nCPUs_m4 = str(input("Number of processors used for this mesh: ")) +#folder_project_m4 = str(input("Type folder where project data is located: ")) +pressure_filename_m4 = path + folder_project_m4 + "\\" + nCPUs_m4 + file +pressure_m4 = np.loadtxt(pressure_filename_m4, skiprows=1,) +time_m4 = np.linspace(0,T_cyc*cyc,pressure_m4.shape[0]) + +# Load 5th Mesh results +N_ts_m5 = 15870 +nCPUs_m5 = "96" +folder_project_m5 = "MCS_m05" +#N_ts_m5 = int(input("Number of Time Steps used for this mesh? ")) +#nCPUs_m5 = str(input("Number of processors used for this mesh: ")) +#folder_project_m5 = str(input("Type folder where project data is located: ")) +pressure_filename_m5 = path + folder_project_m5 + "\\" + nCPUs_m5 + file +pressure_m5 = np.loadtxt(pressure_filename_m5, skiprows=1,) +time_m5 = np.linspace(0,T_cyc*cyc,pressure_m5.shape[0]) + + +for n in range(0,pressure_m4.shape[1]-1): + plt.figure() +# l1, = plt.plot(time_m1,pressure_m1[:,n]/1333.22,'b', label='0.2M') + l2, = plt.plot(time_m2,pressure_m2[:,n]/1333.22,'r', label='0.4M') + l3, = plt.plot(time_m3,pressure_m3[:,n]/1333.22,'g', label='0.8M') + l4, = plt.plot(time_m4,pressure_m4[:,n]/1333.22,'y', label='0.8M') + l5, = plt.plot(time_m5,pressure_m5[:,n]/1333.22,'k', label='0.8M') + plt.show() + + +#n=2 +#l2, = plt.plot(time_m2[round(N_ts_m2-N_ts_m2/cyc):N_ts_m2],pressure_m2[round(N_ts_m2-N_ts_m2/cyc):N_ts_m2,n]/1333.22,'r', label='0.4M') +#l3, = plt.plot(time_m3[round(N_ts_m3-N_ts_m3/cyc):N_ts_m3],pressure_m3[round(N_ts_m3-N_ts_m3/cyc):N_ts_m3,n]/1333.22,'g', label='0.8M') +#l4, = plt.plot(time_m4[round(N_ts_m4-N_ts_m4/cyc):N_ts_m4],pressure_m4[round(N_ts_m4-N_ts_m4/cyc):N_ts_m4,n]/1333.22,'y', label='0.8M') +#l5, = plt.plot(time_m5[round(N_ts_m5-N_ts_m5/cyc):N_ts_m5],pressure_m5[round(N_ts_m4-N_ts_m5/cyc):N_ts_m5,n]/1333.22,'k', label='0.8M') \ No newline at end of file diff --git a/post_Process.py b/post_Process.py new file mode 100644 index 0000000..17b4a42 --- /dev/null +++ b/post_Process.py @@ -0,0 +1,118 @@ + #!/usr/bin/env python +import numpy as np +import matplotlib.pyplot as plt + +T_cyc = 0.477 +N_ts = 14280 +cyc = 5 + +pressure_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\MCS_m04\\96-procs_case\\PHistRCR.dat" +#pressure_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\sim03_RCR_2cyc\\PHistRCR_ok.dat" +#pressure = np.loadtxt(pressure_filename) +pressure = np.loadtxt(pressure_filename, skiprows=1,) + +time = np.linspace(0,T_cyc*cyc,pressure.shape[0]) + +# Ploting Pressure over iterations +# use this to have the plot in different window %matplotlib qt +n=2 +plt.figure() +plt.plot(time,pressure[:,n]/1333.22) +plt.ylabel('Pressure [mmHg]') +plt.xlabel('Time [s]') +plt.title('Presure evolution ROI-2') +axes = plt.gca() +xmax = time[-1] +xmin = time[0] +axes.set_xlim([xmin,xmax]) +ymin = min(pressure[:,n]/1333.22) +ymax = max(pressure[:,n]/1333.22+1) +axes.set_ylim([ymin,ymax]) +plt.show() + +# Ploting Pressure in all the vessels same plot over time +n = pressure.shape +i=0 +plt.figure() +while (i