diff --git a/codes/Graphics.py b/codes/Graphics.py index 03eb121..91d2685 100644 --- a/codes/Graphics.py +++ b/codes/Graphics.py @@ -132,7 +132,6 @@ def Plot_Histogram(meshnames,paths,colors,outpath): ind = 0 for nm in range(len(meshnames)): - if 'box' in meshnames[nm]: # VOXEL SEQUENCE SIZE Lx = 64 @@ -167,13 +166,11 @@ def Plot_Histogram(meshnames,paths,colors,outpath): heights_method = False for k in range(mesh.num_cells()): - # Coordinates of each cell (A,B,C,D) A = Points[Cells[k]][0] B = Points[Cells[k]][1] C = Points[Cells[k]][2] D = Points[Cells[k]][3] - if heights_method: hs = np.zeros([4]) for l in range(4): @@ -1412,127 +1409,79 @@ def ROUTINE(options): if 'Error-curves' in options: if options['Error-curves']['apply']: - print('--- Error curves ---') - - meas = {} - corrs = {} - times = {} - norm = options['Error-curves']['norm'] - colorset = options['Error-curves']['colors'] - - for resol in options['Error-curves']['resol']: - meas[resol] = {} - corrs[resol] = {} - times[resol] = {} - for dts in options['Error-curves']['times']: - path = options['Error-curves']['folder'] + \ - resol + '/dt' + dts - try: - corrs[resol][dts] = np.loadtxt( - path + '/w_' + str(norm) + 'norm.txt') - except IOError: - print('no curve for ' + resol + '/' + dts) - - try: - meas[resol][dts] = np.loadtxt( - path + '/u_' + str(norm) + 'norm.txt') - times[resol][dts] = np.loadtxt( - path + '/times_' + str(norm) + 'norm.txt') - except IOError: - print('no meas for ' + resol + '/' + dts) - - - if 'meas' in options['Error-curves']['mode']: - for resol in options['Error-curves']['resol']: - for dts in meas[resol].keys(): - if resol == 'H1': - nc = 0 - elif resol == 'H2': - nc = 1 - elif resol == 'H3': - nc = 2 - plt.plot(times[resol][dts], meas[resol] - [dts], color=colorset[nc], label='\$' + resol + '/' + dts + '\$') - - #plt.ylim([0, 1.1]) - plt.ylabel('\$||u||_2\$', fontsize=18) - plt.xlabel('\$time \ \ (s)\$', fontsize=18) - plt.legend(fontsize=16) - plt.savefig(options['Error-curves']['outpath'] + 'meas_' +dts + '.png', - dpi=500, bbox_inches='tight') - plt.show() - - if 'corrs' in options['Error-curves']['mode']: - for resol in options['Error-curves']['resol']: - for dts in corrs[resol].keys(): - if resol == 'H1': - nc = 0 - elif resol == 'H2': - nc = 1 - elif resol == 'H3': - nc = 2 - plt.plot(times[resol][dts], corrs[resol] - [dts], color = colorset[nc],label='\$' + resol + '/' + dts + '\$') - - plt.ylim([0, 16.2]) - plt.ylabel('\$||w||_2\$', fontsize=18) - plt.xlabel('\$time \ \ (s)\$', fontsize=18) - plt.legend(fontsize=16) - plt.savefig(options['Error-curves']['outpath'] + 'corrs_' +dts + '.png', - dpi=500, bbox_inches='tight') - plt.show() - - if 'fracs' in options['Error-curves']['mode']: - for resol in options['Error-curves']['resol']: - for dts in corrs[resol].keys(): - if resol == 'H1': - nc = 0 - elif resol == 'H2': - nc = 1 - elif resol == 'H3': - nc = 2 - plt.plot(times[resol][dts], corrs[resol][dts]/meas[resol][dts], - color=colorset[nc], label='\$' + resol + '/' + dts + '\$') - - plt.ylim([0, 1]) - plt.ylabel('\$||w||_2/||u||_2\$', fontsize=18) - plt.xlabel('\$time \ \ (s)\$', fontsize=18) - plt.legend(fontsize=16) - plt.savefig(options['Error-curves']['outpath'] + 'fracs_' +dts + '.png', - dpi=500, bbox_inches='tight') - plt.show() - - if 'Components' in options: - if options['Components']['apply']: - print('--- Components analysis ---') - - for types in options['Components']['type']: + print('--- Error-curves analysis ---') + ratio = False + for types in options['Error-curves']['type']: + if types=='mean_ratio': + types = 'mean' + ratio = True + if types=='max_ratio': + types = 'max' + ratio = True nc = 0 - for subf in options['Components']['subfolders']: + if len(options['Error-curves']['subfolders'])==0: ucomp = [] wcomp = [] - colorset = options['Components']['colors'] - path = options['Components']['folder'] + subf + '/' + colorset = options['Error-curves']['colors'] + path = options['Error-curves']['folder'] try: ucomp = np.loadtxt(path + '/u'+types+'.txt') wcomp = np.loadtxt(path + '/w'+types+'.txt') times = np.loadtxt(path + '/times.txt') except IOError: - print('no components for ' + subf) - - - plt.plot( - times, ucomp, color=colorset[nc], linestyle='-', label= '\$u'+ subf +'\$' ) - - plt.plot( - times, wcomp, color=colorset[nc], linestyle='--', label='\$w'+subf+'\$') + print('no Error-curves for ' + subf) + if not ratio: + plt.plot( + times, ucomp, color=colorset[nc], linestyle='-', label= '\$u\$' ) + plt.plot( + times, wcomp, color=colorset[nc], linestyle='--', label='\$w\$') + else: + plt.plot( + times, ucomp, color=colorset[nc], linestyle='-', label= '\$u\$' ) + plt.plot( + times, wcomp, color=colorset[nc], linestyle='--', label='\$w\$') nc +=1 + else: + for subf in options['Error-curves']['subfolders']: + ucomp = [] + wcomp = [] + colorset = options['Error-curves']['colors'] + labelset = options['Error-curves']['labels'] + path = options['Error-curves']['folder'] + subf + '/' + try: + ucomp = np.loadtxt(path + '/u'+types+'.txt') + wcomp = np.loadtxt(path + '/w'+types+'.txt') + times = np.loadtxt(path + '/times.txt') + except IOError: + print('no cError-curves for ' + subf) + + if not ratio: + plt.plot( + times, ucomp, color=colorset[nc], linestyle='-', label= '\$u'+ subf +'\$' ) + + plt.plot( + times, wcomp, color=colorset[nc], linestyle='--', label='\$w'+subf+'\$') + else: + wu = wcomp/ucomp + plt.plot( + times, wu, color=colorset[nc], linestyle='-', label= '\$'+ labelset[nc] +'\$' ) + nc +=1 + #plt.ylim([0, 170]) - plt.ylabel('\$velocity \ \ (cm/s)\$', fontsize=18) plt.xlabel('\$time \ \ (s)\$', fontsize=18) plt.legend(fontsize=16) - plt.savefig(options['Components']['outpath'] + types + '.png', dpi=500, bbox_inches='tight') + if options['Error-curves']['title']: + plt.title(options['Error-curves']['title'], fontsize=18) + + if not ratio: + plt.ylabel('\$velocity \ \ (cm/s)\$', fontsize=18) + plt.savefig(options['Error-curves']['outpath'] + types + '.png', dpi=500, bbox_inches='tight') + else: + plt.ylabel('\$w/u\$', fontsize=18) + if 'max' in types: + plt.ylim([0, 1.8]) + plt.savefig(options['Error-curves']['outpath'] + types + '_ratio.png', dpi=500, bbox_inches='tight') plt.show() diff --git a/codes/PostCheck.py b/codes/PostCheck.py index 95fbf73..f3d75e4 100644 --- a/codes/PostCheck.py +++ b/codes/PostCheck.py @@ -175,6 +175,37 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, 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) @@ -216,10 +247,10 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, 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 + B0 = 1.5 # Tesla Magnetic Field Strenght TE = 5e-3 # Echo-time - Phi1 = gamma*10*TE + 0*uvec - Phi2 = gamma*10*TE + np.pi*uvec/VENC + 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) @@ -476,7 +507,10 @@ def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comna 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') + 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() @@ -917,7 +951,6 @@ def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): 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']) @@ -930,11 +963,6 @@ def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): 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' - if indexes[k] < 10: - path_u = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename[0]+'.h5' - path_p = checkpoint_path + '0' + \ - str(indexes[k]) + '/'+filename[1]+'.h5' u = Function(V) p = Function(W) @@ -946,7 +974,7 @@ def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): hdf_p.read(p, 'p/vector_0') hdf_p.close() - print('Computing flows and pressures at timestep number ' + str(k)) + 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))) @@ -1134,10 +1162,21 @@ def ROUTINE(options): MESH = LOADmesh(options['Perturbation']['meshpath']) checkpath = options['Perturbation']['checkpath'] + 'checkpoint/' - out_check = options['Perturbation']['checkpath'] + 'Perturbation/' + 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__':