aloma/arnau_op
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

original files

Browse Source
This commit is contained in:
Aloma Blanch 2020-01-08 12:40:10 +01:00
parent 0ce13667eb
commit bd94aa8d65
2 changed files with 211 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
MCS.py Normal file
View File

@ -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')

118
post_Process.py Normal file
View File

@ -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<n[1]):
plt.plot(time[round(N_ts-N_ts/cyc):N_ts],pressure[round(N_ts-N_ts/cyc):N_ts,i]/1333.22, label='ROI-'+str(i+2))
i=i+1
plt.ylabel('Pressure [mmHg]')
plt.xlabel('Time [s]')
#plt.title('5th Cardiac Cylce')
axes = plt.gca()
#xmax = time[-1]
#xmin = time[round(N_ts-N_ts/cyc)]
#axes.set_xlim([xmin,xmax])
#ymin = np.amin(pressure[round(N_ts-N_ts/cyc):N_ts,:]/1333.22)
#ymax = np.amax(pressure[round(N_ts-N_ts/cyc):N_ts,:]/1333.22)+1
#axes.set_ylim([ymin,ymax])
plt.legend(loc='upper right')
plt.show()
# Ploting Pressure over time
i=0
while (i<n[1]):
plt.plot(time,pressure[:,i]/1333.22)
#plt.plot(round(N_ts/cyc), flowRate[round(N_ts/cyc),i], 'ro')
plt.ylabel('Pressure [mmHg]')
plt.xlabel('Time [s]')
plt.title('Pressure waveform ROI-'+str(i+2)+' [cgs]')
#axes = plt.gca()
xmax = time[-1]
xmin = time[0]
axes.set_xlim([xmin,xmax])
ymin = min(pressure[:,0]/1333.22)
ymax = max(pressure[:,0]/1333.22)
axes.set_ylim([ymin,ymax])
plt.show()
i = i + 1
# Plot Flow Rate to analyze that inlet is as expected all outlets all cycles
flowRate_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\MCS_m04\\96-procs_case\\QHistRCR.dat"
#flowRate_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\sim03_RCR_cyc\\QHistRCR_ok.dat"
#flowRate = np.loadtxt(flowRate_filename)
flowRate = np.loadtxt(flowRate_filename, skiprows=1,)
n = flowRate.shape
i=0
while (i<n[1]):
plt.plot(time,flowRate[:,i])
#plt.plot(round(N_ts/cyc), flowRate[round(N_ts/cyc),i], 'ro')
plt.ylabel('Flow Rate Q [mL/s]')
plt.xlabel('Time [s]')
plt.title('Outlet Flow Rate ROI-'+str(i+2)+ ' [cgs]')
#axes = plt.gca()
xmax = time[-1]
xmin = time[0]
axes.set_xlim([xmin,xmax])
ymin = min(flowRate[:,0])
ymax = max(flowRate[:,0])
axes.set_ylim([ymin,ymax])
plt.show()
i = i + 1
# Modifications need to be done in order to Plot only the cycle that user's want
cycle = int(input("Which cycle you whant to plot (velocity plot)? "))
n = flowRate.shape
i=0
while (i<n[1]):
plt.plot(time[round(N_ts-N_ts/cyc):N_ts],flowRate[round(N_ts-N_ts/cyc):N_ts,i])
i = i + 1
plt.ylabel('Flow Rate Q [mL/s]')
plt.xlabel('Time [s]')
plt.title('Out Flow Rate cylce '+str(cycle)+ '- all ROI [cgs]')
#axes = plt.gca()
xmax = time[-1]
xmin = time[0]
axes.set_xlim([xmin,xmax])
ymin = min(flowRate[:,0])
ymax = max(flowRate[:,0])
axes.set_ylim([ymin,ymax])
plt.show()