79 lines
2.8 KiB
Python
79 lines
2.8 KiB
Python
|
#!/usr/bin/env python
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
import matplotlib.pyplot as plt
|
||
|
|
import matplotlib.ticker as mtick
|
||
|
|
from tkinter import Tk
|
||
|
|
from tkinter.filedialog import askopenfilename, asksaveasfile, askdirectory
|
||
|
|
import pandas as pd
|
||
|
|
import tkinter as tk
|
||
|
|
from scipy.signal import find_peaks
|
||
|
|
|
||
|
|
|
||
|
|
def error_plot(folder,t_step,r_criteria,save):
|
||
|
|
# Load iterations and residual error
|
||
|
|
histor = folder + '/histor.dat'
|
||
|
|
input = open(histor, 'r')
|
||
|
|
output = open(folder + "/histor_better.dat", 'w')
|
||
|
|
output.writelines(line.strip() +'\n' for line in input)
|
||
|
|
input.close()
|
||
|
|
output.close()
|
||
|
|
error_info = pd.read_csv(folder + "/histor_better.dat", sep=' ', header=None, usecols=(0,1,2))
|
||
|
|
|
||
|
|
# Select only last iteration of residual error
|
||
|
|
error=[]
|
||
|
|
for n in range(1,error_info.shape[0]):
|
||
|
|
if (error_info[0][n]>error_info[0][n-1]):
|
||
|
|
error.append(error_info[2][n-1])
|
||
|
|
|
||
|
|
time = np.linspace(start = t_step, stop = len(error)*t_step, num = len(error))
|
||
|
|
|
||
|
|
# Plots of interest
|
||
|
|
# Liniar Scale
|
||
|
|
plt.figure()
|
||
|
|
plt.plot(time,error)
|
||
|
|
plt.plot(time,r_criteria*np.ones(len(error)),'r')
|
||
|
|
plt.ylabel('Residual error')
|
||
|
|
plt.xlabel('Time steps')
|
||
|
|
plt.title('Last nonlinear residual error for each time step')
|
||
|
|
plt.grid(True)
|
||
|
|
if save: plt.savefig(plt_folder + case + '_Last_nonlin_res_error.pdf')
|
||
|
|
|
||
|
|
# Semilog scale
|
||
|
|
plt.figure()
|
||
|
|
plt.semilogy(time,error)
|
||
|
|
plt.semilogy(time,r_criteria*np.ones(len(error)),'r:')
|
||
|
|
plt.ylabel('Residual error')
|
||
|
|
plt.xlabel('Time steps')
|
||
|
|
plt.title('Log - Last nonlinear residual error for each time step')
|
||
|
|
plt.grid(True)
|
||
|
|
if save: plt.savefig(plt_folder + case + '_Log_Last_nonlin_res_error.pdf')
|
||
|
|
|
||
|
|
|
||
|
|
def periodicity(project,folder,dt,T_cyc,n_cyc):
|
||
|
|
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
|
||
|
|
time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
|
||
|
|
peak_P = []
|
||
|
|
peak_P_pos = []
|
||
|
|
Nc = round(T_cyc/dt)
|
||
|
|
for i in range(0,n_cyc):
|
||
|
|
peak_P.append(np.amax(pressure[i*Nc:Nc*(i+1),-1])/1333.22)
|
||
|
|
peak_P_pos.append(np.where(pressure[i*Nc:Nc*(i+1),-1] == np.amax(pressure[i*Nc:Nc*(i+1),-1]))[0][0]+Nc*i)
|
||
|
|
|
||
|
|
peak_Pdiff = [peak_P[n]-peak_P[n-1] for n in range(1,len(peak_P))]
|
||
|
|
peak_Pdiff = list(map(abs, peak_Pdiff))
|
||
|
|
|
||
|
|
fig, ax = plt.subplots()
|
||
|
|
ax.plot(time,pressure[:,-1]/1333.22)
|
||
|
|
ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike')
|
||
|
|
ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
|
||
|
|
title='Pressure @ last outlet')
|
||
|
|
ax.spines['right'].set_visible(False)
|
||
|
|
ax.spines['top'].set_visible(False)
|
||
|
|
ax.legend(loc=0)
|
||
|
|
plt.show()
|
||
|
|
|
||
|
|
if (peak_Pdiff[-1]<=1): print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))
|
||
|
|
|
||
|
|
|