# Python script for reading and visualising the LBM simulation # results of water retention behaviour in a CT scanned subvolume # of Hamburg Sand # --------------------------------------------------------------- from matplotlib import cm import numpy as np import matplotlib.pyplot as plt from matplotlib import rcParams # Colors: n = 9 my_color = plt.cm.jet(np.linspace(0,1,n)) my_linewidth = 2.0 # Linewidth in plots # ---------------------------------------------------------------------------- SMALL_SIZE = 8 MEDIUM_SIZE = 10 BIGGER_SIZE = 12 plt.rc('font', size=MEDIUM_SIZE) # controls default text sizes plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title # --------------------------- rcParams['font.family'] = 'serif' # ----------- COLOR = 'black' plt.rcParams['grid.color'] = COLOR plt.rcParams['text.color'] = COLOR plt.rcParams['axes.labelcolor'] = COLOR plt.rcParams['xtick.color'] = COLOR plt.rcParams['ytick.color'] = COLOR plt.rcParams['legend.edgecolor'] = COLOR plt.rc('axes',edgecolor=COLOR) print("Reading unsat_results.txt...") my_path = 'unsat_results.txt' time, step_number, Sr, p_mean_gas, p_mean_liquid, pc, gas_cells, vapour_cells, liquid_cells, pore_cells, solid_cells = np.loadtxt(my_path, unpack = True, delimiter = ',', skiprows = 1, usecols = (0,1,2,3,4,5,6,7,8,9,10)) print("Done.") equilibrium_data = True # Show equilibrium data only? savefigs = True # Save figures to file? if equilibrium_data == True: print("----------------------------------------------------------") print("Only equilibrium data at increments of 1000 ts is plotted.") print("----------------------------------------------------------") else: print("--------------------") print("All data is plotted.") print("--------------------") # ------------------------------------ # Select only equilibrium data points: max_time = time[len(time)-1] inc_time = 1000.0 #count_vect = count_vect.astype(int) if max_time >= inc_time: num_incs = int(max_time/inc_time)-1 count_vect = np.zeros(num_incs+1) count_vect[0] = 0 #998 for i in range(0,len(count_vect)-1): count_vect[i+1] = count_vect[i]+1000 time_new = np.array(time)[count_vect.astype(int)] step_number_new = np.array(step_number)[count_vect.astype(int)] Sr_new = np.array(Sr)[count_vect.astype(int)] pc_new = np.array(pc)[count_vect.astype(int)] p_mean_gas_new = np.array(p_mean_gas)[count_vect.astype(int)] p_mean_liquid_new = np.array(p_mean_liquid)[count_vect.astype(int)] gas_cells_new = np.array(gas_cells)[count_vect.astype(int)] liquid_cells_new = np.array(liquid_cells)[count_vect.astype(int)] vapour_cells_new = np.array(vapour_cells)[count_vect.astype(int)] else: time_new = time step_number_new = step_number Sr_new = Sr pc_new = pc p_mean_liquid_new = p_mean_liquid p_mean_gas_new = p_mean_gas gas_cells_new = gas_cells liquid_cells_new = liquid_cells vapour_cells_new = vapour_cells if equilibrium_data == True: time = time_new step_number = step_number_new Sr = Sr_new pc = pc_new p_mean_gas = p_mean_gas_new p_mean_liquid = p_mean_liquid_new gas_cells = gas_cells_new liquid_cells = liquid_cells_new vapour_cells = vapour_cells_new all_cells = pore_cells[0] + solid_cells[0] n = pore_cells[0]/all_cells e = n/(1-n) print("==============") print("n = ", n) print("e = ", e) print("==============") #================================================= step_0_vect = np.append(np.where(step_number==0)[0], np.where(step_number==1)[0][0]) step_1_vect = np.append(np.where(step_number==1)[0], np.where(step_number==2)[0][0]) step_2_vect = np.append(np.where(step_number==2)[0], np.where(step_number==3)[0][0]) step_3_vect = np.append(np.where(step_number==3)[0], np.where(step_number==4)[0][0]) step_4_vect = np.append(np.where(step_number==4)[0], np.where(step_number==5)[0][0]) step_5_vect = np.append(np.where(step_number==5)[0], np.where(step_number==6)[0][0]) step_6_vect = np.append(np.where(step_number==6)[0], np.where(step_number==7)[0][0]) step_7_vect = np.append(np.where(step_number==7)[0], np.where(step_number==8)[0][0]) step_8_vect = np.where(step_number==8)[0] #================================================= plot_1 = plt.figure(1) ## plt.subplot(2,1,1) ax = plt.gca() ax.patch.set_facecolor('white') plt.plot(time[step_0_vect],Sr[step_0_vect], label = 'Step 0', linestyle='-', color = 'k', marker = 'None', markerfacecolor = 'k', markeredgecolor = 'k', markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_1_vect],Sr[step_1_vect], label = 'Step 1', linestyle='-', color = my_color[0], marker = 'None', markerfacecolor = my_color[0], markeredgecolor = my_color[0], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_2_vect],Sr[step_2_vect], label = 'Step 2', linestyle='-', color = my_color[1], marker = 'None', markerfacecolor = my_color[1], markeredgecolor = my_color[1], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_3_vect],Sr[step_3_vect], label = 'Step 3', linestyle='-', color = my_color[2], marker = 'None', markerfacecolor = my_color[2], markeredgecolor = my_color[2], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_4_vect],Sr[step_4_vect], label = 'Step 4', linestyle='-', color = my_color[3], marker = 'None', markerfacecolor = my_color[3], markeredgecolor = my_color[3], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_5_vect],Sr[step_5_vect], label = 'Step 5', linestyle='-', color = my_color[4], marker = 'None', markerfacecolor = my_color[4], markeredgecolor = my_color[4], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_6_vect],Sr[step_6_vect], label = 'Step 6', linestyle='-', color = my_color[5], marker = 'None', markerfacecolor = my_color[5], markeredgecolor = my_color[5], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_7_vect],Sr[step_7_vect], label = 'Step 7', linestyle='-', color = my_color[6], marker = 'None', markerfacecolor = my_color[6], markeredgecolor = my_color[6], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_8_vect],Sr[step_8_vect], label = 'Step 8', linestyle='-', color = my_color[7], marker = 'None', markerfacecolor = my_color[7], markeredgecolor = my_color[7], markersize = 3.0, linewidth = my_linewidth) # plt.ylabel('$S_{\mathrm{r}}$ [-]') # plt.grid(True,color='k',linestyle='--') plt.ylim([0, 1]) # --- plt.subplot(2,1,2) ax = plt.gca() ax.patch.set_facecolor('white') plt.plot(time[step_0_vect],pc[step_0_vect], label = 'Step 0', linestyle='-', color = 'k', marker = 'None', markerfacecolor = 'k', markeredgecolor = 'k', markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_1_vect],pc[step_1_vect], label = 'Step 1', linestyle='-', color = my_color[0], marker = 'None', markerfacecolor = my_color[0], markeredgecolor = my_color[0], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_2_vect],pc[step_2_vect], label = 'Step 2', linestyle='-', color = my_color[1], marker = 'None', markerfacecolor = my_color[1], markeredgecolor = my_color[1], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_3_vect],pc[step_3_vect], label = 'Step 3', linestyle='-', color = my_color[2], marker = 'None', markerfacecolor = my_color[2], markeredgecolor = my_color[2], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_4_vect],pc[step_4_vect], label = 'Step 4', linestyle='-', color = my_color[3], marker = 'None', markerfacecolor = my_color[3], markeredgecolor = my_color[3], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_5_vect],pc[step_5_vect], label = 'Step 5', linestyle='-', color = my_color[4], marker = 'None', markerfacecolor = my_color[4], markeredgecolor = my_color[4], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_6_vect],pc[step_6_vect], label = 'Step 6', linestyle='-', color = my_color[5], marker = 'None', markerfacecolor = my_color[5], markeredgecolor = my_color[5], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_7_vect],pc[step_7_vect], label = 'Step 7', linestyle='-', color = my_color[6], marker = 'None', markerfacecolor = my_color[6], markeredgecolor = my_color[6], markersize = 3.0, linewidth = my_linewidth) # plt.plot(time[step_8_vect],pc[step_8_vect], label = 'Step 8', linestyle='-', color = my_color[7], marker = 'None', markerfacecolor = my_color[7], markeredgecolor = my_color[7], markersize = 3.0, linewidth = my_linewidth) # plt.ylabel('$p_{\mathrm{c}}$ [mu/(lu$\,$ts$^2$)]') plt.xlabel('$t$ [ts]') plt.grid(True,color='k',linestyle='--') plt.ylim([0, 0.01]) plot_1.subplots_adjust(bottom=0.3) handles, labels = ax.get_legend_handles_labels() lgd = ax.legend(handles, labels, loc ='lower center', ncol = 3, bbox_to_anchor = (0.5,-0.85), fancybox=True,shadow=False,facecolor='white', framealpha=1) plot_1.set_size_inches(20/2.54,14/2.54) plot_1.show() # ---------------------------------------------------------------------------- SMALL_SIZE = 12 MEDIUM_SIZE = 16 BIGGER_SIZE = 18 plt.rc('font', size=MEDIUM_SIZE) # controls default text sizes plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title plot_2 = plt.figure(2) ax = plt.gca() ax.patch.set_facecolor('white') plt.plot(pc[step_0_vect],Sr[step_0_vect], label = 'Step 0', linestyle='-', color = 'k', marker = 'None', markerfacecolor = 'k', markeredgecolor = 'k', markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_1_vect],Sr[step_1_vect], label = 'Step 1', linestyle='-', color = my_color[0], marker = 'None', markerfacecolor = my_color[0], markeredgecolor = my_color[0], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_2_vect],Sr[step_2_vect], label = 'Step 2', linestyle='-', color = my_color[1], marker = 'None', markerfacecolor = my_color[1], markeredgecolor = my_color[1], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_3_vect],Sr[step_3_vect], label = 'Step 3', linestyle='-', color = my_color[2], marker = 'None', markerfacecolor = my_color[2], markeredgecolor = my_color[2], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_4_vect],Sr[step_4_vect], label = 'Step 4', linestyle='-', color = my_color[3], marker = 'None', markerfacecolor = my_color[3], markeredgecolor = my_color[3], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_5_vect],Sr[step_5_vect], label = 'Step 5', linestyle='-', color = my_color[4], marker = 'None', markerfacecolor = my_color[4], markeredgecolor = my_color[4], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_6_vect],Sr[step_6_vect], label = 'Step 6', linestyle='-', color = my_color[5], marker = 'None', markerfacecolor = my_color[5], markeredgecolor = my_color[5], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_7_vect],Sr[step_7_vect], label = 'Step 7', linestyle='-', color = my_color[6], marker = 'None', markerfacecolor = my_color[6], markeredgecolor = my_color[6], markersize = 3.0, linewidth = my_linewidth) # plt.plot(pc[step_8_vect],Sr[step_8_vect], label = 'Step 8', linestyle='-', color = my_color[7], marker = 'None', markerfacecolor = my_color[7], markeredgecolor = my_color[7], markersize = 3.0, linewidth = my_linewidth) # plt.ylabel('$S_{\mathrm{r}}$ [-]') plt.xlabel('$p_{\mathrm{c}}$ [mu/(lu$\,$ts$^2$)]') plt.legend(facecolor='white', framealpha=1, ncol = 3) plt.grid(True,color='k',linestyle='--') plt.xlim([0, 0.01]) plt.ylim([0, 1]) plot_2.set_size_inches(18/2.54,16/2.54) plot_2.show() if savefigs == True: print('Saving figures to file...') plot_1.savefig("Sr_and_s_vs_t_LBM_Sand.pdf", bbox_inches='tight') plot_2.savefig("WRC_LBM_Sand.pdf", bbox_inches='tight') print('Done.') else: print('Figures not saved to file.') print("End of program.") raise SystemExit