from dolfin import * import numpy as np import sys import os from mpi4py import MPI from common import inout def LOADmesh(pathtomesh): mesh = Mesh() hdf = HDF5File(mesh.mpi_comm(), pathtomesh, 'r') hdf.read(mesh, '/mesh', False) boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1) hdf.read(boundaries, '/boundaries') #Evel = VectorElement('P',mesh.ufl_cell(),1) #V = FunctionSpace(mesh,Evel) #V1 = VectorElement('P', mesh.ufl_cell(), 1) #V = FunctionSpace(mesh, V1) P1 = FiniteElement('P', mesh.ufl_cell(), 1) V1 = VectorElement('P', mesh.ufl_cell(), 1) V = FunctionSpace(mesh, V1*P1) if rank == 0: print('MESH ' + ' created: H = ' + str(mesh.hmin()) + ' H_max = ' + str(mesh.hmax())) MESH = {} MESH['mesh'] = mesh MESH['FEM'] = V MESH['boundaries'] = boundaries return MESH def LOADsequences(loadpath): # for loading existing sequences alocated in loadpath if rank == 0: print('{reading} ' + loadpath) S = np.load(loadpath) Sqx = S['x'] Sqy = S['y'] Sqz = S['z'] return [Sqx, Sqy, Sqz] def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options): #from dolfin import HDF5File V = MESH['FEM'] W = MESH['FEM'].sub(0).collapse() # Checkpoints folders unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() if mode == 'u': v = Function(W) xdmf_u = XDMFFile(output_path+outname+'.xdmf') dt = options['Velocity']['dt'] for k in range(0,len(indexes),options['Velocity']['undersampling']): te = k*dt path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v.rename('velocity', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') if 'w' in filename: hdf.read(v, 'w/vector_0') else: hdf.read(v, 'u/vector_0') hdf.close() xdmf_u.write(v, te) if mode == 'w': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v = Function(V) v.rename('waux', outname) comm = MPI.COMM_WORLD hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'w/vector_0') hdf.close() wvec = v.vector().get_local() wnorm = wvec - np.mean(wvec) v.vector()[:] = wnorm xdmf_u.write(v, te) te = te + dt numt = numt + 1 if mode == 'gradw': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v = Function(V) gw = Function(W) gw.rename(outname, outname) comm = MPI.COMM_WORLD hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'w/vector_0') hdf.close() wvec = v.vector().get_local() wnorm = wvec - np.mean(wvec) v.vector()[:] = wnorm gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) xdmf_u.write(gw, te) te = te + dt numt = numt + 1 if mode == 'p' or mode == 'p_cib': xdmf_p = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' if filename == 'p': if k < 10 and k > 0: path = checkpoint_path + '0' + \ str(indexes[k]) + '/'+filename+'.h5' p = Function(W) if mode == 'p': barye2mmHg = 1/1333.22387415 if mode == 'p_cib': barye2mmHg = -0.00750062 p.rename('pressure', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(p, 'p/vector_0') hdf.close() p1vec = p.vector().get_local() p1vec = p1vec - np.mean(p1vec) p.vector()[:] = p1vec*barye2mmHg xdmf_p.write(p, te) te = te + dt numt = numt + 1 if mode == 'divu': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' if indexes[k] < 10 and indexes[k] > 0: path = checkpoint_path + '0' + \ str(indexes[k]) + '/'+filename+'.h5' v = Function(V) dv = Function(W) dv.rename('div', outname) comm = MPI.COMM_WORLD hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'u/vector_0') hdf.close() dv.assign(project(div(v), W)) xdmf_u.write(dv, te) te = te + dt numt = numt + 1 def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, options, flow=False, bnds=None): #from dolfin import HDF5File V = MESH['FEM'] W = MESH['FEM'].sub(0).collapse() # Checkpoints folders unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() if flow: QQ = {} ds = Measure('ds', domain=MESH['mesh'], subdomain_data=MESH['boundaries']) n = FacetNormal(MESH['mesh']) for bb in bnds: QQ[bb] = [] if mode == 'Histograms': u = Function(W) if 'w_COR' in filename: h5name = 'w/vector_0' else: h5name = 'u/vector_0' ValuesTime = np.zeros([len(indexes),len(u.vector().get_local())]) ValuesPeak = np.zeros([len(u.vector().get_local())]) ix = int(0) for k in indexes: path = checkpoint_path + str(k) + '/'+filename+'.h5' hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') print('Reading checkpoint ',indexes[k]) hdf.read(u, h5name) time = hdf.attributes(h5name).to_dict()['timestamp'] hdf.close() uvec = u.vector().get_local() ValuesTime[ix,:] = np.abs(uvec) ix +=1 CoordPeak = np.where(ValuesTime==np.max(ValuesTime)) ValuesPeak = np.abs(ValuesTime[CoordPeak[0],:]) freq,edges = np.histogram(ValuesPeak, bins=80, density=True) #Saving the histogram print('Saving at ' + output_path) np.savetxt(output_path + 'hist_frew.txt', freq) np.savetxt(output_path + 'hist_edges.txt', edges) if mode == 'perturbation': u = Function(W) unew = Function(W) u.rename('velocity', outname) unew.rename('velocity', outname) if options['Perturbation']['xdmf']: xdmf_u = XDMFFile(output_path+'u.xdmf') if not options['Perturbation']['type']['SNR']=='inf': Noise = True def Add_Noise(signal,SNR): Psignal = signal**2 Psignal_av = np.mean(Psignal) Psignal_av_db = 10*np.log10(Psignal_av) Pnoise_av_db = Psignal_av_db - SNR Pnoise_av = 10**(Pnoise_av_db/10) noise_std = np.sqrt(Pnoise_av) noise = np.random.normal(0,noise_std,len(signal)) return signal + noise else: Noise = False if not options['Perturbation']['type']['phase_contrast']==0: Phase_Contrast = True else: Phase_Contrast = False noise_in_coil = options['Perturbation']['type']['coil'] for k in indexes: path = checkpoint_path + str(k) + '/'+filename+'.h5' hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(u, 'u/vector_0') time = hdf.attributes('u/vector_0').to_dict()['timestamp'] hdf.close() uvec = u.vector().get_local() if Phase_Contrast: ufactor = options['Perturbation']['type']['phase_contrast']/100 VENC = np.max(np.abs(uvec))*ufactor gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei B0 = 1.5 # Tesla Magnetic Field Strenght TE = 5e-3 # Echo-time Phi1 = 1*B0*TE + 0*uvec Phi2 = 1*B0*TE + np.pi*uvec/VENC M1 = np.cos(Phi1) + 1j*np.sin(Phi1) M2 = np.cos(Phi2) + 1j*np.sin(Phi2) if noise_in_coil: SNR = options['Perturbation']['type']['SNR'] m1x_new = Add_Noise(np.real(M1),SNR) m1y_new = Add_Noise(np.imag(M1),SNR) m2x_new = Add_Noise(np.real(M2),SNR) m2y_new = Add_Noise(np.imag(M2),SNR) M1_new = m1x_new + 1j*m1y_new M2_new = m2x_new + 1j*m2y_new uvec = (np.angle(M2_new) - np.angle(M1_new))*VENC/np.pi unew.vector()[:] = uvec else: uvec = (np.angle(M2) - np.angle(M1))*VENC/np.pi else: if noise_in_coil: raise Exception('In order to perturb in coils some PC should be selected') if not noise_in_coil: if Noise: SNR = options['Perturbation']['type']['SNR'] unew.vector()[:] = Add_Noise(uvec,SNR) else: unew.vector()[:] = uvec print('Writing checkpoint number ',k) write_path = output_path + 'checkpoint/{i}/'.format(i=k) hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') hdf2.write(unew, '/u', time) hdf2.close() if options['Perturbation']['xdmf']: xdmf_u.write(unew, time) if mode == 'interpolation': dt_new = options['Temporal-Interpolation']['dt_new'] dt = options['Temporal-Interpolation']['dt'] Nfolder = len(indexes) Tf = dt*(Nfolder-1) Nfolder_new = np.int(Tf/dt_new + 1) indexes_new = [k for k in range(1,Nfolder_new+1)] u = Function(W) unew = Function(W) u.rename('velocity', outname) unew.rename('velocity', outname) if options['Temporal-Interpolation']['xdmf']: xdmf_u = XDMFFile(output_path+'u.xdmf') uvec = {} # Reading first all the points for every original time-steps for k in indexes: path = checkpoint_path + str(k) + '/'+filename+'.h5' print('Reading timestep number: ' + str(k)) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') if 'w' in filename: hdf.read(u, 'w/vector_0') else: hdf.read(u, 'u/vector_0') hdf.close() uvec[k] = u.vector().get_local() times_old = np.linspace(0,Tf,Nfolder) times_new = np.linspace(0,Tf,Nfolder_new) velnodes = np.zeros(times_old.size) velnodes_new = np.zeros(times_new.size) uvec_new = {} for k in indexes_new: uvec_new[k] = np.zeros(uvec[1].size) from scipy.interpolate import interp1d print('Interpolating in every node across time') # FOR every single node!!! for l in range(len(uvec[1])): for k in range(len(velnodes)): velnodes[k] = uvec[indexes[k]][l] inter_f = interp1d(times_old, velnodes , kind='cubic') velnodes_new = inter_f(times_new) for k in range(len(velnodes_new)): uvec_new[indexes_new[k]][l] = velnodes_new[k] print('Interpolation done') for k in indexes_new: print('Writing timestep number: ' + str(k) ) unew.vector()[:] = uvec_new[k][:] write_path = output_path + 'checkpoint/{i}/'.format(i=k) hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') hdf2.write(unew, '/u', float((k-1)*dt_new)) hdf2.close() if options['Temporal-Interpolation']['xdmf']: xdmf_u.write(unew, (k-1)*dt_new) if mode == 'average': N_av = options['Temporal-Average']['subsampling_rate'] dt = options['Temporal-Average']['dt'] ks = 1 mean_fill = False u = Function(W) usub = Function(W) u.rename('velocity', outname) usub.rename('velocity', outname) if options['Temporal-Average']['xdmf']: xdmf_u = XDMFFile(output_path+'test.xdmf') for k in range(len(indexes)): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' print('Reading timestep number: ' + str(indexes[k])) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') if 'w' in filename: hdf.read(u, 'w/vector_0') else: hdf.read(u, 'u/vector_0') hdf.close() if not mean_fill: mean_u = u.vector().get_local()/N_av mean_fill = True else: mean_u += u.vector().get_local()/N_av if np.mod(k+1,N_av) == 0: mean_fill = False usub.vector()[:] = mean_u print('saving timestep number: ' , ks) write_path = output_path + 'checkpoint/{i}/'.format(i=ks) hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') hdf2.write(usub, '/u', float((k-N_av+1)*dt)) hdf2.close() if options['Temporal-Average']['xdmf']: xdmf_u.write(usub, (k-N_av+1)*dt) ks +=1 if mode == 'u': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' u = Function(V) u.rename('velocity', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(u, 'u/vector_0') hdf.close() if flow: for bb in bnds: QQ[bb].append(assemble(dot(u, n)*ds(bb))) xdmf_u.write(u, k) te = te + dt if flow: for bb in bnds: np.savetxt(output_path+'flowrate'+str(bb)+'.txt', QQ[bb]) if mode == 'w': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v = Function(V) v.rename('waux', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'w/vector_0') hdf.close() wvec = v.vector().get_local() wnorm = wvec - np.mean(wvec) v.vector()[:] = wnorm xdmf_u.write(v, k) if mode == 'gradw': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v = Function(V) gw = Function(W) gw.rename(outname, outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'w/vector_0') hdf.close() wvec = v.vector().get_local() wnorm = wvec - np.mean(wvec) v.vector()[:] = wnorm gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) xdmf_u.write(gw, k) if mode == 'p' or mode == 'p_cib': xdmf_p = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' p = Function(W) if mode == 'p': barye2mmHg = -1/1333.22387415 if mode == 'p_cib': barye2mmHg = -0.00750062 p.rename('pressure', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(p, 'p/vector_0') hdf.close() p1vec = p.vector().get_local()*barye2mmHg p.vector()[:] = p1vec xdmf_p.write(p, k) if mode == 'divu': xdmf_u = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' v = Function(V) dv = Function(W) dv.rename('div', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(v, 'u/vector_0') hdf.close() dv.assign(project(div(v), W)) xdmf_u.write(dv, k) def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comname,options): from dolfin import HDF5File V = MESH['FEM'] W = MESH['FEM'].sub(0).collapse() unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() unsort_indexes0 = os.listdir(reference_path) indexes0 = [int(x) for x in unsort_indexes0] indexes0.sort() if len(indexes)!=len(indexes0): raise Exception('The lengh of the checkpoints are not the same!') if mode =='curves': if options['Error-curves']['undersampling']>1: print('undersampling in the reference measurements!') n_under = options['Error-curves']['undersampling'] indexes0 = indexes0[0::n_under] u = Function(W) w = Function(W) for typs in options['Error-curves']['type']: ucomp = [] wcomp = [] times = [] dt = options['Error-curves']['dt'] for k in range(1,len(indexes)): path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') if 'w' in comname: hdf_w.read(w, 'w/vector_0') else: hdf_w.read(w, 'u/vector_0') hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') hdf_u.close() hdf_w.close() u_vec = u.vector().get_local() w_vec = w.vector().get_local() print('computing error in timestep numer',k) if typs == 'mean': ucomp.append(np.mean(abs(u_vec))) wcomp.append(np.mean(abs(w_vec))) elif typs == 'max': ucomp.append(np.max(abs(u_vec))) wcomp.append(np.max(abs(w_vec))) else: raise Exception('No defined type for curve printing!') times.append(k*dt) np.savetxt(outpath+'u' +typs+'.txt',ucomp) np.savetxt(outpath+'w' +typs+'.txt',wcomp) np.savetxt(outpath+'times.txt',times) if mode == 'histogram': from pathlib import Path import pickle u = Function(V) w = Function(V) errors = {} for k in range(len(indexes)): path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' if indexes0[k] < 10: path_u = reference_path + '0' + \ str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') hdf_w.read(w, 'w/vector_0') hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') hdf_u.close() hdf_w.close() u_vec = u.vector().get_local() w_vec = w.vector().get_local() errors[k] = np.zeros(u_vec.size) for l in range(len(errors[k])): #errors[k][l] = np.nan if u_vec[l]<1e-9: errors[k][l] = -1 else: eta = np.abs(w_vec[l]/u_vec[l]) if np.abs(eta)>50: errors[k][l] = -1 else: errors[k][l] = eta write_path = Path(outpath) fpath = write_path.joinpath('errors.dat') pickle.dump(errors, fpath.open('wb')) if mode == 'h5': xdmf_u = XDMFFile(outpath+'meas.xdmf') #xdmf_ucor = XDMFFile(output_path+'ucor.xdmf') xdmf_w = XDMFFile(outpath+'w.xdmf') ds = Measure("ds", subdomain_data=MESH['boundaries']) u = Function(W) w = Function(W) #ucor = Function(W) for k in range(0, len(indexes), 1): path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' if indexes0[k] < 10: path_u = reference_path + '0' + \ str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') #ucor.rename('ucor', 'ucor') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') hdf_w.read(w, 'w/vector_0') hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') hdf_u.close() hdf_w.close() #u_vec = u.vector().get_local() #w_vec = w.vector().get_local() #ucor.vector()[:] = u_vec + w_vec xdmf_u.write(u, k) #xdmf_ucor.write(ucor, k) xdmf_w.write(w, k) if mode == 'colormap': colormap = XDMFFile(outpath+'colormap.xdmf') #ds = Measure("ds", subdomain_data=MESH['boundaries']) u = Function(W) w = Function(W) cm = Function(W) dt = options['Corrector']['dt'] for k in range(len(indexes)): print('making the colormap in the time',np.round(k*dt,2)) path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' u.rename('meas', 'meas') w.rename('w', 'w') cm.rename('color','color') hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') hdf_w.read(w, 'w/vector_0') hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') hdf_u.close() hdf_w.close() uvec = u.vector().get_local() wvec = w.vector().get_local() cm.vector()[:] = np.sqrt((uvec - wvec)**2) colormap.write(cm, k*dt) if mode == 'error_u': xdmf_u = XDMFFile(output_path+'error_u.xdmf') for k in range(0, len(indexes), 1): path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' if k < 10: path = checkpoint_path + '0'+str(indexes[k]) + '/u.h5' path0 = reference_path + '0'+str(indexes[k]) + '/u.h5' u = Function(V) u0 = Function(V) eu = Function(W) eu.rename('error_u', 'error_u') comm = MPI.COMM_WORLD hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(u, 'u/vector_0') hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') hdf0.read(u0, 'u/vector_0') hdf0.close() hdf.close() eu.assign(project(sqrt(inner(u-u0, u-u0)), W)) xdmf_u.write(eu, te) te = te + dt numt = numt + 1 if mode == 'error_p' or mode == 'error_p_cib': xdmf_p = XDMFFile(output_path+outname+'.xdmf') for k in range(0, len(indexes0), 1): path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' if refname == 'p': if k < 10 and k > 0: path0 = reference_path + '0' + \ str(indexes0[k]) + '/'+refname+'.h5' #zero_point = [12.79353, 14.32866, 6.51101] #zero_point = np.array(zero_point) #ndim = W.mesh().topology().dim() #zero_point = W.tabulate_dof_coordinates().reshape((-1, ndim))[0] if 'cib' in mode: barye2mmHg = 0.00750062 else: barye2mmHg = 1/1333.22387415 p = Function(W) p0 = Function(W) p.rename('error_p', outname) hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(p, 'p/vector_0') hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') hdf0.read(p0, 'p/vector_0') hdf.close() hdf0.close() #p0.vector()[:] -= p0(zero_point) pvec = p.vector().get_local() p0vec = p0.vector().get_local() pvec = pvec - np.mean(pvec) p0vec = p0vec - np.mean(p0vec) errvec = np.sqrt((pvec - p0vec)**2)*barye2mmHg #errvec = np.abs((pvec - p0vec)/(p0vec)) p.vector()[:] = errvec #p.vector()[:] = p0vec xdmf_p.write(p, te) te = te + dt numt = numt + 1 def SEQCIBH5(pathtocib, pathtocibseq, ndata): import csv times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] #times = [8] if ndata == 1: pathtocib += '9_mm_AoCo_phantom_rest_vox0.9mm/' pathtocibseq += '9AoCo_rest_vox0.9/' elif ndata == 2: pathtocib += '9_mm_AoCo_phantom_stress_vox0.9mm/' pathtocibseq += '9AoCo_stress_vox0.9/' elif ndata == 3: pathtocib += '11_mm_AoCo_phantom_rest_vox0.9mm/' pathtocibseq += '11AoCo_rest_vox0.9/' elif ndata == 4: pathtocib += '11_mm_AoCo_phantom_stress_vox0.9mm/' pathtocibseq += '11AoCo_stress_vox0.9/' else: raise Exception('Data file not recognize!') xdmf_v = XDMFFile(pathtocibseq + 'vel.xdmf') mesh = Mesh(pathtocib+'aorta.xml') Evel = VectorElement('Lagrange', mesh.ufl_cell(), 1) V = FunctionSpace(mesh, Evel) boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) VX = V.sub(0).collapse() VY = V.sub(1).collapse() VZ = V.sub(2).collapse() vv_x = Function(VX) vv_y = Function(VY) vv_z = Function(VZ) v = Function(V) v.rename('velocity', 'vel') mx = dof_to_vertex_map(VX) my = dof_to_vertex_map(VY) mz = dof_to_vertex_map(VZ) for ti in times: print('reading aorta '+str(ti) + '.csv') with open(pathtocib + 'aorta'+str(ti)+'.csv', 'r') as csvfile: mylist = [row[0] for row in csv.reader(csvfile, delimiter=';')] Values = np.array(mylist) file_x = np.zeros([len(Values)]) file_y = np.zeros([len(Values)]) file_z = np.zeros([len(Values)]) for l in range(len(Values)): row = Values[l].split(',') file_x[l] = float(row[0]) file_y[l] = float(row[1]) file_z[l] = float(row[2]) # print(np.max(file_z)) #print(np.where(file_z==np.max(file_z)) ) # print(np.max(file_x[np.where(file_z==np.max(file_z))])) # print(np.max(file_y[np.where(file_z==np.max(file_z))])) #S = np.load(pathtocibseq) #ux = S['x'] #uy = S['y'] #uz = S['z'] #ux = ux.transpose((0,2,1,3)) #uy = uy.transpose((0,2,1,3)) #uz = uz.transpose((0,2,1,3)) vv_x.vector()[:] = file_x[mx] vv_y.vector()[:] = file_y[my] vv_z.vector()[:] = file_z[mz] assign(v.sub(0), vv_x) assign(v.sub(1), vv_y) assign(v.sub(2), vv_z) xdmf_v.write(v, ti) # LagrangeInterpolator.interpolate(vv_x,v1_x) # LagrangeInterpolator.interpolate(vv_y,v1_y) # LagrangeInterpolator.interpolate(vv_z,v1_z) xdmf_v.close() def SqtoH5(BOX, MESH, Sqx, Sqy, Sqz, output_path, uname): xdmf_u = XDMFFile(output_path+uname+'.xdmf') Xt = BOX['nodes'] SqXt = BOX['seq'] P1 = BOX['FEM'] V = MESH['FEM'] [Lx, Ly, Lz, Lt] = Sqx.shape VX = V.sub(0).collapse() VY = V.sub(1).collapse() VZ = V.sub(2).collapse() vv_x = Function(VX) vv_y = Function(VY) vv_z = Function(VZ) m = dof_to_vertex_map(P1) if rank == 0: print('{SQtoH5} total number of timesteps: ' + str(Lt)) v = Function(V) v.rename('velocity', uname) dt = (Lt-1) te = 0 for t in range(Lt): if rank == 0: print('timestep number', t) v1_x = Function(P1) v1_y = Function(P1) v1_z = Function(P1) values_x = np.zeros(v1_x.vector().get_local().shape) values_y = np.zeros(v1_y.vector().get_local().shape) values_z = np.zeros(v1_z.vector().get_local().shape) S0x = Sqx[:, :, :, t] S0y = Sqy[:, :, :, t] S0z = Sqz[:, :, :, t] for k in range(Xt.shape[0]): values_x[k] = S0x[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] values_y[k] = S0y[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] values_z[k] = S0z[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] v1_x.vector()[:] = values_x v1_y.vector()[:] = values_y v1_z.vector()[:] = values_z LagrangeInterpolator.interpolate(vv_x, v1_x) LagrangeInterpolator.interpolate(vv_y, v1_y) LagrangeInterpolator.interpolate(vv_z, v1_z) assign(v.sub(0), vv_x) assign(v.sub(1), vv_y) assign(v.sub(2), vv_z) xdmf_u.write(v, te) te = te+dt def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): #from dolfin import HDF5File V = MESH['FEM'] W = MESH['FEM'].sub(1).collapse() # Checkpoints folders unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() from pathlib import Path import pickle # Create the Spheres if options['Estim_Pressure']['method'] == 'spheres': import mshr Npts0 = options['Estim_Pressure']['spheres'][0]['Npts'] Npts1 = options['Estim_Pressure']['spheres'][1]['Npts'] center0 = options['Estim_Pressure']['spheres'][0]['center'] center1 = options['Estim_Pressure']['spheres'][1]['center'] radius0 = options['Estim_Pressure']['spheres'][0]['radius'] radius1 = options['Estim_Pressure']['spheres'][1]['radius'] x0 = Point(np.array(center0, dtype=float)) x1 = Point(np.array(center1, dtype=float)) sphere0 = mshr.Sphere(x0, radius0) sphere1 = mshr.Sphere(x1, radius1) mesh_s1 = mshr.generate_mesh(sphere0, Npts0) mesh_s2 = mshr.generate_mesh(sphere1, Npts1) VS1 = FunctionSpace(mesh_s1, FiniteElement('P', mesh_s1.ufl_cell(), 1)) VS2 = FunctionSpace(mesh_s2, FiniteElement('P', mesh_s2.ufl_cell(), 1)) s1 = Function(VS1) s2 = Function(VS2) p = Function(W) one_mesh = interpolate(Constant(1), W) LagrangeInterpolator.interpolate(s1, one_mesh) LagrangeInterpolator.interpolate(s2, one_mesh) vol1 = assemble(s1*dx) vol2 = assemble(s1*dx) dt = options['Estim_Pressure']['dt'] p_drop_lst = [] time_ = [] for k in range(0, len(indexes)): path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') hdf.read(p, 'p/vector_0') hdf.close() LagrangeInterpolator.interpolate(s1, p) LagrangeInterpolator.interpolate(s2, p) P1 = assemble(s1*dx)/vol1 P2 = assemble(s2*dx)/vol2 p_drop_lst.append(P2-P1) time_.append(k*dt) if rank == 0: print('Pressure drop :', P2-P1) # Saving the Result write_path = Path(outpath) write_path.mkdir(exist_ok=True) methods = filename[2:] data = { 'pdrop': np.array(p_drop_lst), 'time': np.array(time_) } fpath = write_path.joinpath('pdrop_' + methods + '.dat') pickle.dump(data, fpath.open('wb')) def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): #from dolfin import HDF5File V = MESH['FEM'].sub(0).collapse() W = MESH['FEM'].sub(1).collapse() # Checkpoints folders unsort_indexes = os.listdir(checkpoint_path) indexes = [int(x) for x in unsort_indexes] indexes.sort() QQ = {} PP = {} ds = Measure('ds', domain=MESH['mesh'], subdomain_data=MESH['boundaries']) n = FacetNormal(MESH['mesh']) ones = interpolate(Constant(1), W) for bb in bnds: PP[bb] = [] QQ[bb] = [] for k in range(0, len(indexes)): path_u = checkpoint_path + str(indexes[k]) + '/'+filename[0]+'.h5' path_p = checkpoint_path + str(indexes[k]) + '/'+filename[1]+'.h5' u = Function(V) p = Function(W) #comm = MPI.COMM_WORLD hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') hdf_u.read(u, 'u/vector_0') hdf_u.close() hdf_p = HDF5File(MESH['mesh'].mpi_comm(), path_p, 'r') hdf_p.read(p, 'p/vector_0') hdf_p.close() print('Computing flows and pressures at timestep number ' + str(indexes[k])) for bb in bnds: if bb == 2: QQ[bb].append(-assemble(dot(u, n)*ds(bb))) else: QQ[bb].append(assemble(dot(u, n)*ds(bb))) PP[bb].append(assemble(p*ds(bb))/1333.22387415 / assemble(ones*ds(bb))) print('saving the figure at' + output_path) from matplotlib import pyplot as plt from matplotlib import rc rc('text', usetex=True) fig = plt.figure(figsize=(12, 5), dpi=100) t = np.linspace(0, len(indexes), len(indexes)) plt.subplot(1, 2, 2) for bb in bnds: plt.plot(t, QQ[bb], linewidth=1.2, linestyle='-', label='$ Outlet: '+str(bb)+'$') plt.xlabel(r'$frames $', fontsize=20) plt.ylabel(r'$ Q $', fontsize=20) plt.xlim([0, 1.05*max(t)]) plt.legend(fontsize=14) plt.subplot(1, 2, 1) for bb in bnds: plt.plot(t, PP[bb], linewidth=1.2, linestyle='-', label='$ Outlet: '+str(bb)+'$') plt.xlabel(r'$frames $', fontsize=20) plt.ylabel(r'$ P \ \ mmHg$', fontsize=20) plt.xlim([0, 1.05*max(t)]) plt.legend(fontsize=14) fig.savefig(output_path + 'Outlet_Windkessel', dpi=500) def ROUTINE(options): if 'Outlet_Wind' in options: if options['Outlet_Wind']['apply']: if rank == 0: print('--- Outlet Windkessel ---') mesh_path = options['Outlet_Wind']['mesh_path'] out_path = options['Outlet_Wind']['out_path'] filename = options['Outlet_Wind']['filename'] checkpoint = options['Outlet_Wind']['checkpoint'] + 'checkpoint/' bnds = options['Outlet_Wind']['bnds'] MESH = LOADmesh(mesh_path) OUTLETwind(MESH, out_path, checkpoint, filename, bnds) if 'Corrector' in options: if options['Corrector']['apply']: if rank == 0: print('Applying Corrector') MESH = LOADmesh(options['Corrector']['mesh_path']) u_path = options['Corrector']['u_path'] + 'checkpoint/' w_path = options['Corrector']['w_path'] + 'checkpoint/' wname = options['Corrector']['wname'] uname = options['Corrector']['uname'] mode = options['Corrector']['mode'] outpath = options['Corrector']['outpath'] ERRORmap(MESH, mode, outpath, u_path, w_path, uname, wname, options) if 'Velocity' in options: if options['Velocity']['apply']: if rank == 0: print('--- Reading Velocity ---') MESH = LOADmesh(options['Velocity']['mesh_path']) filename = options['Velocity']['filename'] checkpoint_path = options['Velocity']['checkpoint'] + 'checkpoint/' outpath = options['Velocity']['checkpoint'] WORKcheck(MESH, 'u', outpath, checkpoint_path, filename, filename, options) if 'Vel_from_check' in options: if options['Vel_from_check']['apply']: if rank == 0: print('Applying Velocity--MAP from sequence') [Sqx, Sqy, Sqz] = LOADsequences(options['Velocity']['pathtoseq']) BOX = LOADmesh('cib') AORTA = LOADmesh(pathtomesh) SqtoH5(BOX, AORTA, Sqx, Sqy, Sqz, output_path, uname) if 'W_map' in options: if options['W_map']['apply']: if rank == 0: print('Applying W--MAP') MESH = LOADmesh(options['mesh_path']) filename = options['W_map']['filename'] outname = options['W_map']['out_name'] mode = 'w' WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options) if 'GradW_map' in options: if options['GradW_map']['apply']: if rank == 0: print('Applying Grad W--MAP') MESH = LOADmesh(options['mesh_path']) filename = options['GradW_map']['filename'] outname = options['GradW_map']['out_name'] mode = 'gradw' WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options) if 'Pressure' in options: if options['Pressure']['apply']: if rank == 0: print('Applying Pressure-MAP') MESH = LOADmesh(options['mesh_path']) filename = options['Pressure']['filename'] outname = options['Pressure']['out_name'] mode = 'p' WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options) if 'Error_P' in options: if options['Error_P']['apply']: if rank == 0: print('Applying L2 error to Pressure') MESH = LOADmesh(options['mesh_path']) refname = options['Error_P']['refname'] reference_path = options['Error_P']['refpath'] filename = options['Error_P']['filename'] outname = options['Error_P']['out_name'] mode = 'error_p' ERRORmap(MESH, mode, output_path, reference_path, checkpoint_path, refname, filename, outname, options) if 'Estim_Pressure' in options: if options['Estim_Pressure']['apply']: if rank == 0: print('Applying Pressure Estimator') MESH = LOADmesh(options['mesh_path']) filename = options['Estim_Pressure']['filename'] outpath = options['Estim_Pressure']['outpath'] ESTIMpressure(MESH, outpath, checkpoint_path, filename, options) if 'Error-curves' in options: if options['Error-curves']['apply']: if rank == 0: print('--- Error curves ---' ) MESH = LOADmesh(options['Error-curves']['meshpath']) ref_check = options['Error-curves']['ref_check'] + 'checkpoint/' com_check = options['Error-curves']['com_check'] + 'checkpoint/' outpath = options['Error-curves']['outpath'] refname = options['Error-curves']['ref_name'] comname = options['Error-curves']['com_name'] ERRORmap(MESH, 'curves', outpath, ref_check, com_check, refname, comname, options) if 'Temporal-Average' in options: if options['Temporal-Average']['apply']: if rank == 0: print('--- Temporal Average ---') MESH = LOADmesh(options['Temporal-Average']['meshpath']) ref_check = options['Temporal-Average']['original_check'] + 'checkpoint/' out_check = options['Temporal-Average']['out_check'] READcheckpoint(MESH,'average', out_check,ref_check,'u','u',options) if 'Temporal-Interpolation' in options: if options['Temporal-Interpolation']['apply']: if rank == 0: print('--- Temporal Interpolation ---') MESH = LOADmesh(options['Temporal-Interpolation']['meshpath']) ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/' out_check = options['Temporal-Interpolation']['out_check'] READcheckpoint(MESH,'interpolation', out_check,ref_check,'u','u',options) if 'Perturbation' in options: if options['Perturbation']['apply']: if rank==0: print('--- Perturbation in measurements ---') MESH = LOADmesh(options['Perturbation']['meshpath']) checkpath = options['Perturbation']['checkpath'] + 'checkpoint/' name_folder = 'SNR' +options['Perturbation']['type']['SNR'] + 'V' + str(options['Perturbation']['type']['phase_contrast']) out_check = options['Perturbation']['checkpath'] + name_folder + '/' READcheckpoint(MESH,'perturbation', out_check,checkpath,'u','u',options) if 'Histograms' in options: if options['Histograms']['apply']: if rank==0: print('--- Histograms in measurements ---') MESH = LOADmesh(options['Histograms']['meshpath']) checkpath = options['Histograms']['checkpath'] + 'checkpoint/' outpath = options['Histograms']['checkpath'] field_name = options['Histograms']['field_name'] READcheckpoint(MESH,'Histograms', outpath, checkpath,field_name,field_name,options) if __name__ == '__main__': comm = MPI.COMM_WORLD size = comm.Get_size() rank = comm.Get_rank() if len(sys.argv) > 1: if os.path.exists(sys.argv[1]): inputfile = sys.argv[1] if rank == 0: print('Found input file ' + inputfile) else: raise Exception('Command line arg given but input file does not exist:' ' {}'.format(sys.argv[1])) else: raise Exception('An input file is required as argument!') options = inout.read_parameters(inputfile) ROUTINE(options)