#import tkinter #import matplotlib #matplotlib.use('TkAgg') import numpy as np import matplotlib.pyplot as plt import sys import os import math # # from matplotlib import rc # # #rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) # # rc('font',**{'family':'serif','serif':['Helvetica']}) # # rc('text', usetex=True) # import matplotlib.pyplot as plt # #plt.rcParams["font.family"] = "Arial" # plt.rcParams.update({ # "text.usetex": True, # "font.family": 'Helvetica' # }) from pathlib import Path def save_two_vectors_to_csv(vec1, vec2, out_path, filename): """ Save two vectors as a 2-column CSV file (no header). Parameters ---------- vec1 : array-like First vector (column 1) vec2 : array-like Second vector (column 2) out_path : str or Path Folder where the CSV will be saved filename : str Name of the CSV file (e.g., "output.csv") """ vec1 = np.asarray(vec1) vec2 = np.asarray(vec2) if vec1.shape != vec2.shape: raise ValueError("Vectors must have the same length.") out_dir = Path(out_path) out_dir.mkdir(parents=True, exist_ok=True) file_path = out_dir / filename data = np.column_stack((vec1, vec2)) np.savetxt(file_path, data, delimiter=",", fmt="%.8f") return file_path caseName = sys.argv[1] #Specificy case that you are plotting vec = np.loadtxt("driver.out") file_key = open("gpvar2", "r") #gpvar2 content = file_key.readlines() file_key.close() codeV = [] indexV = [] for i in range(len(content)): value = content[i].split("=") code = str(value[0]) index = int(value[1]) codeV.append(code) indexV.append(index) def lookup(valtext): indF = codeV.index(valtext) index = indexV[indF] return index-1 def fc_crack(cg, lg): p_int = 0.25*(math.pi**2)*(cg**3)/(lg**3) #Crack intersection probability for Bethe Lattice Z = 4 if ((1/p_int) - 0.75) < 0: #Usually this implies percolation is reacher return 1 else: pintf = np.sqrt( (1/p_int) - 0.75) - 0.5 fc = 1 - 4*(pintf**3) + 3*(pintf**4) return min(fc,1) output_vals_poro = ["poro"] output_vals_p = ["p"] output_vals_energy = ["e"] output_vals_strainEv = ["EvT11","EvT22","EvT33","EvT12", "EvT31", "EvT23"] output_vals_main_stress = ["TT11","TT22","TT33"] output_vals_dev_stress = ["TT12","TT31","TT23"] output_vals_permeability = ["K11","K22","K33"] output_vals_permeability2 = ["K12", "K31", "K23"] output_vals_c_length = ["c11","c22","c33"] output_vals_c_dens = ["l11","l22","l33"] output_vals_w = ["w11","w22","w33"] output_vals_dK = ["dKI1","dKII1", "dKI2","dKII2", "dKI3","dKII3"] output_vals_pstrain = ["plastic_strain"] out_wildcard = ["tjump"] #Use here for c33 inc due to dynamic loutput_vals_poro = [r"\phi"] loutput_vals_p = [r"p"] loutput_vals_energy = [r"e"] loutput_vals_strainEv = [r"Ev$_{11}$", r"Ev$_{22}$", r"Ev$_{33}$", r"Ev$_{12}$", r"Ev$_{31}$", r"Ev$_{23}$"] loutput_vals_main_stress = [r"$\sigma_{11}$", r"$\sigma_{22}$", r"$\sigma_{33}$"] loutput_vals_dev_stress = [r"$\sigma_{12}$", r"$\sigma_{31}$", r"$\sigma_{23}$"] loutput_vals_permeability = [r"k$_{11}$", r"k$_{22}$", r"k$_{33}$"] loutput_vals_permeability2 = [r"k$_{12}$", r"k$_{31}$", r"k$_{23}$"] loutput_vals_c_length = [r"c$_{1}$", r"c$_{2}$", r"c$_{3}$"] loutput_vals_c_dens = [r"l$_{1}$", r"l$_{2}$", r"l$_{3}$"] loutput_vals_w = [r"w$_{1}$", r"w$_{2}$", r"w$_{3}$"] loutput_vals_dK = ["dKI1","dKII1", "dKI2","dKII2", "dKI3","dKII3"] time = vec[:,0] c1v = vec[:,lookup(output_vals_c_length[0])] c2v = vec[:,lookup(output_vals_c_length[1])] c3v = vec[:,lookup(output_vals_c_length[2])] l1v = vec[:,lookup(output_vals_c_dens[0])] l2v = vec[:,lookup(output_vals_c_dens[1])] l3v = vec[:,lookup(output_vals_c_dens[2])] w1v = vec[:,lookup(output_vals_w[0])] w2v = vec[:,lookup(output_vals_w[1])] w3v = vec[:,lookup(output_vals_w[2])] porov = vec[:,lookup(output_vals_poro[0])] fc1 = np.zeros(len(c1v)) fc2 = np.zeros(len(c1v)) fc3 = np.zeros(len(c1v)) Sig11 = vec[:,lookup(output_vals_main_stress[0])] Sig22 = vec[:,lookup(output_vals_main_stress[1])] Sig33 = vec[:,lookup(output_vals_main_stress[2])] Sig12 = vec[:,lookup(output_vals_dev_stress[0])] Sig23 = vec[:,lookup(output_vals_dev_stress[1])] Sig31 = vec[:,lookup(output_vals_dev_stress[2])] SigDev11 = np.zeros(len(Sig11)) SigDev22 = np.zeros(len(Sig11)) SigDev33 = np.zeros(len(Sig11)) SigVM = np.zeros(len(Sig11)) kmag1v = np.zeros(len(Sig11)) kmag2v = np.zeros(len(Sig11)) kmag3v = np.zeros(len(Sig11)) crack_porov = np.zeros(len(Sig11)) for i in range(len(c1v)): p = (1/3)*(Sig11[i]+Sig22[i]+Sig33[i]) SigDev11[i] = Sig11[i] - p SigDev22[i] = Sig22[i] - p SigDev33[i] = Sig33[i] - p #SigVM[i] = ( (3/2) * ( Sig11[i]**2 + Sig22[i]**2 + Sig33[i]**2 + 2 * ( Sig12[i]**2 + Sig23[i]**2 + Sig31[i]**2) ) - 0.5 * (Sig11[i]+Sig22[i]+Sig33[i])**2 )**0.5 SigVM[i] = (0.5 * ( (Sig11[i])**2 + (Sig22[i])**2 + (Sig33[i])**2 + 6*(Sig12[i]**2 + Sig23[i]**2 + Sig31[i]**2) ) )**0.5 for i in range(len(c1v)): fc1[i] = fc_crack(c1v[i], l1v[i]) fc2[i] = fc_crack(c2v[i], l2v[i]) fc3[i] = fc_crack(c3v[i], l3v[i]) image2_1D = True shearXZ = True def crack_poro(c11, c22, c33, w11, w22, w33, l11, l22, l33): if image2_1D: return (2 * math.pi * ( ( 1 * (c33))**2 )* 1*(w33)) / ( 1*(l33) )**3 else: return (2 * math.pi * ( (0.33333 * (c11+c22+c33))**2 )* 0.33333*(w11+w22+w33)) / ( 0.333333*(l11+l22+l33) )**3 crack_poro_limit = 0.5 for i in range(len(c1v)): cporo = crack_poro(c1v[i], c2v[i], c3v[i], w1v[i], w2v[i], w3v[i], l1v[i], l2v[i], l3v[i]) crack_porov[i] = min(cporo, crack_poro_limit) kmag1v[i] = (2/15) * (porov[i] + crack_porov[i]) * fc1[i] * w1v[i]**2 kmag2v[i] = (2/15) * (porov[i] + crack_porov[i]) * fc2[i] * w2v[i]**2 kmag3v[i] = (2/15) * (porov[i] + crack_porov[i]) * fc3[i] * w3v[i]**2 #Now let's go beyond time plotting if image2_1D: #################################################################################################### #1D Results output V2 plt.figure(figsize=(23, 20)) plt.subplots_adjust(wspace=0.52, hspace=0.4) # plt.subplot(3,5,1) # index = lookup(output_vals_p[0]) # plt.plot( time[1:] , vec[:,index][1:], ".--", linewidth=0.5 ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"p", size=20) # plt.title("Pressure", size=20) # #plt.legend(fontsize = 18) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,3,3) # index = lookup(output_vals_energy[0]) # plt.plot( time , vec[:,index] ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$e$", size=20) # plt.title("Energy", size=20) # #plt.legend(fontsize = 18) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,3,3) # plt.plot( time , vec[:,lookup(output_vals_strainEv[0])], label= output_vals_strainEv[0]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[1])], label= output_vals_strainEv[1]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[2])], label= output_vals_strainEv[2]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[3])], label= output_vals_strainEv[3]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[4])], label= output_vals_strainEv[4]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[5])], label= output_vals_strainEv[5]) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Ev$", size=20) # plt.title("Elastic Strain", size=20) # plt.legend(fontsize = 16, framealpha=0.7) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) plt.subplot(3,4,1) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[0])][1:], "rs" , label= loutput_vals_main_stress[0], linewidth=0.5 ) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[1])][1:], "g.-", label= loutput_vals_main_stress[1], linewidth=0.5 ) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[2])][1:], "b*--", label= loutput_vals_main_stress[2], linewidth=2 ) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[0])][1:], "ms" , label= loutput_vals_dev_stress[0], linewidth=0.5 ) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[1])][1:], "y.-", label= loutput_vals_dev_stress[1], linewidth=0.5 ) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[2])][1:], "c*--", label= loutput_vals_dev_stress[2], linewidth=2 ) plt.plot( time[1:] , vec[:,lookup(output_vals_p[0])][1:], "k*--", label= loutput_vals_p[0], linewidth=2 ) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma$", size=20) plt.title("Main Stresses and Pressure", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) save_two_vectors_to_csv( time[1:] , c3v[1:]/ c3v[1] , "./output_results", "time_c_.dat") # plt.subplot(3,4,2) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[0])][1:], "rs" , label= loutput_vals_strainEv[0], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[1])][1:], "g.-", label= loutput_vals_strainEv[1], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[2])][1:], "b*--", label= loutput_vals_strainEv[2], linewidth=2 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[3])][1:], "ms" , label= loutput_vals_strainEv[3], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[4])][1:], "y.-", label= loutput_vals_strainEv[4], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_strainEv[5])][1:], "c*--", label= loutput_vals_strainEv[5], linewidth=2 ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\epsilon$", size=20) # plt.title("Strain Components", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) plt.subplot(3,4,2) if shearXZ: plt.plot( vec[:,lookup(output_vals_dev_stress[1])][1:], vec[:,lookup(output_vals_strainEv[4])][1:] ) plt.xlabel(r"Stress ($\sigma_{xz}$)", size=20) plt.ylabel(r"Strain ($\epsilon_{xz}$)", size=20) else: plt.plot( vec[:,lookup(output_vals_main_stress[2])][1:], vec[:,lookup(output_vals_strainEv[2])][1:] ) plt.xlabel(r"Stress ($\sigma_{zz}$)", size=20) plt.ylabel(r"Strain ($\epsilon_{zz}$)", size=20) plt.title("Stress vs. Strain", size=20) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) # plt.subplot(3,4,2) # plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[0])][1:], "ms" , label= loutput_vals_dev_stress[0], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[1])][1:], "y.-", label= loutput_vals_dev_stress[1], linewidth=0.5 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[2])][1:], "c*--", label= loutput_vals_dev_stress[2], linewidth=2 ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\sigma$", size=20) # plt.title("Off-Diagonal Stresses", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) plt.subplot(3,4,3) index = lookup(out_wildcard[0]) if shearXZ: plt.plot( vec[:,lookup(output_vals_dev_stress[1])][1:] , vec[:,index][1:], "g*--", label= 'gf_dyn_I', linewidth=2 ) plt.plot( vec[:,lookup(output_vals_dev_stress[1])][1:] , vec[:,lookup(output_vals_dK[1])][1:], "m.--", label='gf_dyn_II', linewidth=2 ) plt.xlabel(r"Stress ($\sigma_{xz}$)", size=20) else: plt.plot( vec[:,lookup(output_vals_main_stress[2])][1:] , vec[:,index][1:], "g*--", label= 'gf_dyn_I', linewidth=2 ) plt.plot( vec[:,lookup(output_vals_main_stress[2])][1:] , vec[:,lookup(output_vals_dK[1])][1:], "m.--", label='gf_dyn_II', linewidth=2 ) plt.xlabel(r"Stress ($\sigma_{zz}$)", size=20) plt.ylabel(r"$\Delta{K}_{I} / \Delta{K}_{II}$", size=20) plt.title("Stress vs. Dynamic S.I.F. (Mode I and II)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,4,4) if shearXZ: plt.plot( vec[:,lookup(output_vals_dev_stress[1])][1:] , vec[:,lookup(output_vals_pstrain[0])][1:], "g*--", linewidth=2 ) plt.xlabel(r"Stress ($\sigma_{xz}$)", size=20) else: plt.plot( vec[:,lookup(output_vals_main_stress[2])][1:] , vec[:,lookup(output_vals_pstrain[0])][1:], "g*--", linewidth=2 ) plt.xlabel(r"Stress ($\sigma_{zz}$)", size=20) plt.ylabel(r"$\epsilon_{p}$", size=20) plt.title("Stress vs. Plastic Strain", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,4,5) if shearXZ: plt.plot(vec[:,lookup(output_vals_dev_stress[1])][1:], vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1]) plt.xlabel(r"Stress ($\sigma_{xz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_dev_stress[1])][1:], vec[:,lookup(output_vals_w[2])][1:]/ vec[:,lookup(output_vals_w[2])][1], "./output_results", "Sigma_w_.dat") else: plt.plot(vec[:,lookup(output_vals_main_stress[2])][1:], vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1]) plt.xlabel(r"Stress ($\sigma_{zz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_main_stress[2])][1:], vec[:,lookup(output_vals_w[2])][1:]/ vec[:,lookup(output_vals_w[2])][1], "./output_results", "Sigma_w_.dat") plt.ylabel(r"$w$", size=20) plt.title("Stress vs. Aperture", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,4,6) if shearXZ: plt.plot(np.abs(vec[:,lookup(output_vals_strainEv[4])][1:]), vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1]) plt.xlabel(r"Strain ($\epsilon_{xz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_strainEv[4])][1:], vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1] , "./output_results", "Ev_w_.dat") else: plt.plot(np.abs(vec[:,lookup(output_vals_strainEv[2])][1:]), vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1]) plt.xlabel(r"Strain ($\epsilon_{zz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_strainEv[2])][1:], vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1], "./output_results","Ev_w_.dat") plt.ylabel(r"$w$", size=20) plt.title("Strain vs. Aperture", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,4,7) if shearXZ: plt.plot(np.abs(vec[:,lookup(output_vals_dev_stress[1])][1:]), vec[:,lookup(output_vals_c_length[0])][1:] / vec[:,lookup(output_vals_c_length[0])][1], label= 'n1', linewidth=9 ) plt.plot(np.abs(vec[:,lookup(output_vals_dev_stress[1])][1:]), vec[:,lookup(output_vals_c_length[1])][1:] / vec[:,lookup(output_vals_c_length[1])][1], "m-", label= 'n2', linewidth=6 ) plt.plot(np.abs(vec[:,lookup(output_vals_dev_stress[1])][1:]), vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], "y.-", label= 'n3', linewidth=3 ) plt.xlabel(r"Stress ($\sigma_{xz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_dev_stress[1])][1:], vec[:,lookup(output_vals_c_length[2])][1:]/ vec[:,lookup(output_vals_c_length[2])][1] , "./output_results", "Sigma_c_.dat") else: plt.plot(np.abs(vec[:,lookup(output_vals_main_stress[2])][1:]), vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1]) plt.xlabel(r"Stress ($\sigma_{zz}$)", size=20) save_two_vectors_to_csv(vec[:,lookup(output_vals_main_stress[2])][1:], vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1] , "./output_results", "Sigma_c_.dat") plt.ylabel(r"$c$", size=20) plt.title("Stress vs. Length", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,4,8) plt.plot(vec[:,lookup(output_vals_pstrain[0])][1:], vec[:,lookup(output_vals_c_dens[2])][1:] / vec[:,lookup(output_vals_c_dens[2])][1]) plt.xlabel(r"Plastic strain ($\epsilon_{p}$)", size=20) plt.ylabel(r"$l$", size=20) plt.title("Plastic strain vs. Distance", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) save_two_vectors_to_csv(vec[:,lookup(output_vals_pstrain[0])][1:], vec[:,lookup(output_vals_c_dens[2])][1:] / vec[:,lookup(output_vals_c_dens[2])][1] , "./output_results", "pstrain_l_.dat") plt.subplot(3,4,9) plt.plot(vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], fc3[1:]) plt.xlabel(r"$c$", size=20) plt.ylabel(r"$f$", size=20) plt.title("Length vs. Connectivity", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) save_two_vectors_to_csv(vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], fc3[1:], "./output_results", "c_f_.dat") plt.subplot(3,4,10) plt.plot(vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1], kmag3v[1:]/kmag3v[1] ) plt.xlabel(r"$w$", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) plt.title("Aperture vs. Rel. Permeability", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) save_two_vectors_to_csv(vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1], kmag3v[1:]/kmag3v[1], "./output_results", "w_k_.dat") plt.subplot(3,4,11) plt.plot(vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], kmag3v[1:]/kmag3v[1] ) plt.xlabel(r"$c$", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) plt.title("Length vs. Rel. Permeability", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) save_two_vectors_to_csv(vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1] , kmag3v[1:]/kmag3v[1], "./output_results", "c_k_.dat") plt.subplot(3,4,12) plt.plot(vec[:,lookup(output_vals_c_dens[2])][1:] / vec[:,lookup(output_vals_c_dens[2])][1], kmag3v[1:]/kmag3v[1] ) plt.xlabel(r"$l$", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) plt.title("Distance vs. Rel. Permeability", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) # plt.subplot(3,4,3) # plt.plot( time[1:] , SigDev11[1:], "ms" , label= ' TT1\' ', linewidth=0.5 ) # plt.plot( time[1:] , SigDev22[1:], "y.-", label= ' TT2\' ', linewidth=0.5 ) # plt.plot( time[1:] , SigDev33[1:], "c*--", label= ' TT3\' ', linewidth=2 ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\sigma'$", size=20) # plt.title("Main Dev. Stresses", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,5) # # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[0])][1:], "rs" ,label= loutput_vals_dK[0]) # # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[2])][1:], "g.-", label= loutput_vals_dK[2]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[4])][1:], "rs--", label= 'dK_I') # # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[1])][1:], "ms" ,label= loutput_vals_dK[1]) # # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[3])][1:], "y.-", label= loutput_vals_dK[3]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[5])][1:], "b.--", label='dK_II') # index = lookup(out_wildcard[0]) # plt.plot( time[1:] , vec[:,index][1:], "g*--", label= 'gf_dyn_I', linewidth=2 ) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[1])][1:], "m.--", label='gf_dyn_II', linewidth=2 ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\Delta{K}_{I} / \Delta{K}_{II}$", size=20) # plt.title("S.I.F. (Mode I and II)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,4) # plt.plot(time[1:] , SigVM[1:]) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\sigma_{VM}$", size=20) # plt.title("Von Mises Stress", size=20) # #plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,5,9) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[0])][1:] / vec[:,lookup(output_vals_permeability[0])][1], "rs" , label= loutput_vals_permeability[0]) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[1])][1:] / vec[:,lookup(output_vals_permeability[1])][1], "g.-", label= loutput_vals_permeability[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[2])][1:] / vec[:,lookup(output_vals_permeability[2])][1], "b.--", label= loutput_vals_permeability[2]) # plt.xlabel(r"Time ($\mu$s)", size=20) # #plt.ylabel(r"$K$ (mm$^2$)", size=20) # plt.ylabel(r"$Rel. \bar{k}$", size=20) # #plt.title("Permeab. (Main)", size=20) # plt.title("Rel. Permeab. (Main)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,5,10) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[0])][1:] / vec[:,lookup(output_vals_permeability2[0])][1] , "ms" ,label= loutput_vals_permeability2[0]) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[1])][1:] / vec[:,lookup(output_vals_permeability2[1])][1] , "y.-", label= loutput_vals_permeability2[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[2])][1:] / vec[:,lookup(output_vals_permeability2[2])][1] , "c.--", label= loutput_vals_permeability2[2]) # plt.xlabel(r"Time ($\mu$s)", size=20) # #plt.ylabel(r"$K$ (mm$^2$)", size=20) # plt.ylabel(r"$Rel. \bar{k}$", size=20) # #plt.title("Permeab. (Offd)", size=20) # plt.title("Rel. Permeab. (Offd)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,5) # # plt.plot( time[1:] , vec[:,lookup(output_vals_w[0])][1:] / vec[:,lookup(output_vals_w[0])][1] , "rs", label= loutput_vals_w[0]) # # plt.plot( time[1:] , vec[:,lookup(output_vals_w[1])][1:] / vec[:,lookup(output_vals_w[1])][1], "g.-", label= loutput_vals_w[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1], "b.-", label= r'$\bar{w}$') # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Rel. \bar{w}$", size=20) # #plt.ylabel(r"$\bar{w}$ (mm)", size=20) # plt.title("Crack aperture", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,6) # # plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[0])][1:] / vec[:,lookup(output_vals_c_length[0])][1], "rs", label= loutput_vals_c_length[0]) # # plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[1])][1:] / vec[:,lookup(output_vals_c_length[1])][1], "g.-", label= loutput_vals_c_length[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], "b.-", label= r'$\bar{c}$') # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Rel. \bar{c}$", size=20) # #plt.ylabel(r"$\bar{c}$ (mm)", size=20) # plt.title("Crack length", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,7) # # plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[0])][1:] / vec[:,lookup(output_vals_c_dens[0])][1], "rs", label= output_vals_c_dens[0]) # # plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[1])][1:] / vec[:,lookup(output_vals_c_dens[1])][1], "g.-", label= output_vals_c_dens[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[2])][1:] / vec[:,lookup(output_vals_c_dens[2])][1], "b.-", label= r'$\bar{l}$') # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Rel. \bar{l}$", size=20) # #plt.ylabel(r"$\bar{c}$ (mm)", size=20) # plt.title("Crack distance", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,8) # # plt.plot( time[1:] , fc1[1:], "rs" , label= 'fc1') # # plt.plot( time[1:] , fc2[1:], "g.-", label= 'fc2') # plt.plot( time[1:] , fc3[1:], "b.--", label= r'$f_c$') # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$f_c$", size=20) # plt.title("Connectivity Factor", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,9) # index = lookup(output_vals_poro[0]) # plt.plot( time[1:] , vec[:,index][1:], label=r"$\phi$" ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\phi$", size=20) # plt.title("Intrinsic Porosity", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,10) # index = lookup(output_vals_poro[0]) # plt.plot( time[1:] , crack_porov[1:], "rs", label=r"$\phi_c$" ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\phi$", size=20) # plt.title("Crack Porosity", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,11) # # plt.plot( time[1:] , kmag1v[1:]/kmag1v[1], "rs" , label= r'$\bar{\alpha}_{1}$') # # plt.plot( time[1:] , kmag2v[1:]/kmag2v[1], "g.-", label= r'$\bar{\alpha}_{2}$') # plt.plot( time[1:] , kmag3v[1:]/kmag3v[1], "b.--", label= r'$\bar{k} / \bar{k}_o$') # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Rel. \bar{k}$", size=20) # #plt.ylabel(r"$\bar{k}$ (mm$^2$)", size=20) # plt.title("K. Mag. (main)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,4,12) # # plt.plot( time[1:] , kmag1v[1:], "rs" , label= r'$\bar{\alpha}_{1}$') # # plt.plot( time[1:] , kmag2v[1:], "g.-", label= r'$\bar{\alpha}_{2}$') # plt.plot( time[1:] , kmag3v[1:], "b.--", label= r'$\bar{k}$') # plt.xlabel(r"Time ($\mu$s)", size=20) # #plt.ylabel(r"$Rel. \bar{k}$", size=20) # plt.ylabel(r"$\bar{k}$ (mm$^2$)", size=20) # plt.title("K. Mag. (main)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,5,15) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[1])][1:], "ms" ,label= loutput_vals_dK[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[3])][1:], "y.-", label= loutput_vals_dK[3]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[5])][1:], "c.--", label= loutput_vals_dK[5]) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\Delta{K}_{II}$", size=20) # plt.title("S.I.F. (Mode II)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) #plt.savefig("tests_plot.jpeg") plt.savefig("./plot_tests/"+ str(caseName) + "_outputv2.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_outputv2.png") else: ########################################################################################################################################################################### #ORIGINAL OUTPUT V2 FOR 3D plt.figure(figsize=(28, 17)) plt.subplots_adjust(wspace=0.52, hspace=0.4) plt.subplot(3,5,1) index = lookup(output_vals_p[0]) plt.plot( time[1:] , vec[:,index][1:], ".--", linewidth=0.5 ) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"p", size=20) plt.title("Pressure", size=20) #plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) # plt.subplot(3,3,3) # index = lookup(output_vals_energy[0]) # plt.plot( time , vec[:,index] ) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$e$", size=20) # plt.title("Energy", size=20) # #plt.legend(fontsize = 18) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) # plt.subplot(3,3,3) # plt.plot( time , vec[:,lookup(output_vals_strainEv[0])], label= output_vals_strainEv[0]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[1])], label= output_vals_strainEv[1]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[2])], label= output_vals_strainEv[2]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[3])], label= output_vals_strainEv[3]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[4])], label= output_vals_strainEv[4]) # plt.plot( time , vec[:,lookup(output_vals_strainEv[5])], label= output_vals_strainEv[5]) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$Ev$", size=20) # plt.title("Elastic Strain", size=20) # plt.legend(fontsize = 16, framealpha=0.7) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) plt.subplot(3,5,2) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[0])][1:], "rs" , label= loutput_vals_main_stress[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[1])][1:], "g.-", label= loutput_vals_main_stress[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_main_stress[2])][1:], "b.--", label= loutput_vals_main_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma$", size=20) plt.title("Main Stresses", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,3) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[0])][1:], "ms" , label= loutput_vals_dev_stress[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[1])][1:], "y.-", label= loutput_vals_dev_stress[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_dev_stress[2])][1:], "c.--", label= loutput_vals_dev_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma$", size=20) plt.title("Off-Diagonal Stresses", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,4) plt.plot( time[1:] , SigDev11[1:], "ms" , label= ' TT1\' ') plt.plot( time[1:] , SigDev22[1:], "y.-", label= ' TT2\' ') plt.plot( time[1:] , SigDev33[1:], "c.--", label= ' TT3\' ') plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma'$", size=20) plt.title("Main Dev. Stresses", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,5) plt.plot(time[1:] , SigVM[1:]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma_{VM}$", size=20) plt.title("Von Mises Stress", size=20) #plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,6) index = lookup(output_vals_poro[0]) plt.plot( time[1:] , vec[:,index][1:], label=r"$\phi$" ) plt.plot( time[1:] , crack_porov[1:], "rs", label=r"$\phi_c$" ) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\phi$", size=20) plt.title("Porosity", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,7) plt.plot( time[1:] , kmag1v[1:]/kmag1v[1], "rs" , label= r'$\bar{\alpha}_{1}$') plt.plot( time[1:] , kmag2v[1:]/kmag2v[1], "g.-", label= r'$\bar{\alpha}_{2}$') plt.plot( time[1:] , kmag3v[1:]/kmag3v[1], "b.--", label= r'$\bar{\alpha}_{3}$') plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) #plt.ylabel(r"$\bar{k}$ (mm$^2$)", size=20) plt.title("K. Mag. (main)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,8) plt.plot( time[1:] , fc1[1:], "rs" , label= 'fc1') plt.plot( time[1:] , fc2[1:], "g.-", label= 'fc2') plt.plot( time[1:] , fc3[1:], "b.--", label= 'fc3') plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$f_c$", size=20) plt.title("Connectivity Factor", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,9) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[0])][1:] / vec[:,lookup(output_vals_permeability[0])][1], "rs" , label= loutput_vals_permeability[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[1])][1:] / vec[:,lookup(output_vals_permeability[1])][1], "g.-", label= loutput_vals_permeability[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability[2])][1:] / vec[:,lookup(output_vals_permeability[2])][1], "b.--", label= loutput_vals_permeability[2]) plt.xlabel(r"Time ($\mu$s)", size=20) #plt.ylabel(r"$K$ (mm$^2$)", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) #plt.title("Permeab. (Main)", size=20) plt.title("Rel. Permeab. (Main)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,10) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[0])][1:] / vec[:,lookup(output_vals_permeability2[0])][1] , "ms" ,label= loutput_vals_permeability2[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[1])][1:] / vec[:,lookup(output_vals_permeability2[1])][1] , "y.-", label= loutput_vals_permeability2[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_permeability2[2])][1:] / vec[:,lookup(output_vals_permeability2[2])][1] , "c.--", label= loutput_vals_permeability2[2]) plt.xlabel(r"Time ($\mu$s)", size=20) #plt.ylabel(r"$K$ (mm$^2$)", size=20) plt.ylabel(r"$Rel. \bar{k}$", size=20) #plt.title("Permeab. (Offd)", size=20) plt.title("Rel. Permeab. (Offd)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,11) plt.plot( time[1:] , vec[:,lookup(output_vals_w[0])][1:] / vec[:,lookup(output_vals_w[0])][1] , "rs", label= loutput_vals_w[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_w[1])][1:] / vec[:,lookup(output_vals_w[1])][1], "g.-", label= loutput_vals_w[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_w[2])][1:] / vec[:,lookup(output_vals_w[2])][1], "b.--", label= loutput_vals_w[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$Rel. \bar{w}$", size=20) #plt.ylabel(r"$\bar{w}$ (mm)", size=20) plt.title("Crack aperture", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,12) plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[0])][1:] / vec[:,lookup(output_vals_c_length[0])][1], "rs", label= loutput_vals_c_length[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[1])][1:] / vec[:,lookup(output_vals_c_length[1])][1], "g.-", label= loutput_vals_c_length[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_c_length[2])][1:] / vec[:,lookup(output_vals_c_length[2])][1], "b.--", label= loutput_vals_c_length[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$Rel. \bar{c}$", size=20) #plt.ylabel(r"$\bar{c}$ (mm)", size=20) plt.title("Crack length", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,13) plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[0])][1:] / vec[:,lookup(output_vals_c_dens[0])][1], "rs", label= output_vals_c_dens[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[1])][1:] / vec[:,lookup(output_vals_c_dens[1])][1], "g.-", label= output_vals_c_dens[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_c_dens[2])][1:] / vec[:,lookup(output_vals_c_dens[2])][1], "b.--", label= output_vals_c_dens[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$Rel. \bar{l}$", size=20) #plt.ylabel(r"$\bar{c}$ (mm)", size=20) plt.title("Crack distance", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,14) plt.plot( time[1:] , kmag1v[1:], "rs" , label= r'$\bar{\alpha}_{1}$') plt.plot( time[1:] , kmag2v[1:], "g.-", label= r'$\bar{\alpha}_{2}$') plt.plot( time[1:] , kmag3v[1:], "b.--", label= r'$\bar{\alpha}_{3}$') plt.xlabel(r"Time ($\mu$s)", size=20) #plt.ylabel(r"$Rel. \bar{k}$", size=20) plt.ylabel(r"$\bar{k}$ (mm$^2$)", size=20) plt.title("K. Mag. (main)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(3,5,15) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[0])][1:], "rs" ,label= loutput_vals_dK[0]) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[2])][1:], "g.-", label= loutput_vals_dK[2]) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[4])][1:], "b.--", label= loutput_vals_dK[4]) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[1])][1:], "ms" ,label= loutput_vals_dK[1]) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[3])][1:], "y.-", label= loutput_vals_dK[3]) plt.plot( time[1:] , vec[:,lookup(output_vals_dK[5])][1:], "c.--", label= loutput_vals_dK[5]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\Delta{K}_{I} / \Delta{K}_{II}$", size=20) plt.title("S.I.F. (Mode I and II)", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) # plt.subplot(3,5,15) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[1])][1:], "ms" ,label= loutput_vals_dK[1]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[3])][1:], "y.-", label= loutput_vals_dK[3]) # plt.plot( time[1:] , vec[:,lookup(output_vals_dK[5])][1:], "c.--", label= loutput_vals_dK[5]) # plt.xlabel(r"Time ($\mu$s)", size=20) # plt.ylabel(r"$\Delta{K}_{II}$", size=20) # plt.title("S.I.F. (Mode II)", size=20) # plt.legend(fontsize = 16, framealpha=0.2) # plt.tick_params(labelsize=18) # plt.xticks(fontsize=18) # plt.yticks(fontsize=18) #plt.savefig("tests_plot.jpeg") plt.savefig("./plot_tests/"+ str(caseName) + "_outputv2.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_outputv2.png") shear_flag = True #Stress display (3) plt.figure(figsize=(12, 16)) plt.subplots_adjust(wspace=0.6, hspace=0.25) plt.tight_layout(pad=3) plt.subplot(2,1,1) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[0])], "rs" , label= loutput_vals_main_stress[0]) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[1])], "g.-", label= loutput_vals_main_stress[1]) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[2])], "b.--", label= loutput_vals_main_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma$ (GPa)", size=20) plt.title("Main Stresses", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(2,1,2) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[0])], "ms" , label= loutput_vals_dev_stress[0]) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[1])], "y.-", label= loutput_vals_dev_stress[1]) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[2])], "c.--", label= loutput_vals_dev_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=20) plt.ylabel(r"$\sigma$ (GPa)", size=20) plt.title("Off-Diagonal Stresses", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.savefig("./plot_tests/"+ str(caseName) + "_outputv3.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_outputv3.png") #Microvariable change display (4) plt.figure(figsize=(12, 32)) plt.subplots_adjust(wspace=0.6, hspace=0.25) plt.tight_layout(pad=3) plt.subplot(4,1,1) index = lookup(output_vals_poro[0]) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , vec[:,index] ) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=20) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=20) plt.ylabel(r"$\phi$", size=20) plt.title("Porosity", size=20) #plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,2) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , vec[:,lookup(output_vals_w[0])], "rs", label= loutput_vals_w[0]) plt.plot( xcomp , vec[:,lookup(output_vals_w[1])], "g.-", label= loutput_vals_w[1]) plt.plot( xcomp , vec[:,lookup(output_vals_w[2])], "b.--", label= loutput_vals_w[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=20) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=20) plt.ylabel(r"$\bar{w}_{1}$ (mm)", size=20) plt.title("Crack aperture", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,3) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , vec[:,lookup(output_vals_c_length[0])], "rs", label= loutput_vals_c_length[0]) plt.plot( xcomp , vec[:,lookup(output_vals_c_length[1])], "g.-", label= loutput_vals_c_length[1]) plt.plot( xcomp , vec[:,lookup(output_vals_c_length[2])], "b.--", label= loutput_vals_c_length[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=20) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=20) plt.ylabel(r"$\bar{c}_{1}$ (mm)", size=20) plt.title("Crack length", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,4) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , vec[:,lookup(output_vals_c_dens[0])], "ms" ,label= loutput_vals_c_dens[0]) plt.plot( xcomp , vec[:,lookup(output_vals_c_dens[1])], "y.-", label= loutput_vals_c_dens[1]) plt.plot( xcomp , vec[:,lookup(output_vals_c_dens[2])], "c.--", label= loutput_vals_c_dens[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=20) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=20) plt.ylabel(r"$\bar{l}_{1}$ (mm)", size=20) plt.title("Crack spacing", size=20) plt.legend(fontsize = 16, framealpha=0.2) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.savefig("./plot_tests/"+ str(caseName) + "_outputv4.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_outputv4.png") #Microvariable change display (5) plt.figure(figsize=(12, 32)) plt.subplots_adjust(wspace=0.6, hspace=0.25) plt.tight_layout(pad=3) plt.subplot(4,1,1) index = lookup(output_vals_poro[0]) xcomp = vec[:,index] plt.plot( xcomp , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mm$^2$)", size=20) plt.xlabel(r"$\phi$", size=20) plt.title("Permeab. (Main) vs. Porosity", size=20) #plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,2) xcomp = vec[:,lookup(output_vals_w[0])] plt.plot( xcomp , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mm$^2$)", size=20) plt.xlabel(r"$\bar{w}_{1}$ (mm)", size=20) plt.title("Permeab. (Main) vs. Aperture", size=20) plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,3) xcomp = vec[:,lookup(output_vals_c_length[0])] plt.plot( xcomp , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mm$^2$)", size=20) plt.xlabel(r"$\bar{c}_{1}$ (mm)", size=20) plt.title("Permeab. (Main) vs. Length", size=20) plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.subplot(4,1,4) xcomp = vec[:,lookup(output_vals_c_dens[0])] plt.plot( xcomp , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mm$^2$)", size=20) plt.xlabel(r"$\bar{l}_{1}$ (mm)", size=20) plt.title("Permeab. (Main) vs. Distance", size=20) plt.legend(fontsize = 18) plt.tick_params(labelsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.savefig("./plot_tests/"+ str(caseName) + "_outputv5.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_outputv5.png") #Individual saved plots #TIME STRESSES #Stress Main plt.figure(figsize=(15, 11)) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[0])], "rs" , label= loutput_vals_main_stress[0]) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[1])], "g.-", label= loutput_vals_main_stress[1]) plt.plot( time , 1000*vec[:,lookup(output_vals_main_stress[2])], "b.--", label= loutput_vals_main_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=34) plt.ylabel(r"$\sigma$ (GPa)", size=34) plt.title("Main Stresses", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_stress_m.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_stress_m.png") plt.figure(figsize=(15, 11)) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[0])], "ms" , label= loutput_vals_dev_stress[0]) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[1])], "y.-", label= loutput_vals_dev_stress[1]) plt.plot( time , 1000*vec[:,lookup(output_vals_dev_stress[2])], "c.--", label= loutput_vals_dev_stress[2]) plt.xlabel(r"Time ($\mu$s)", size=34) plt.ylabel(r"$\sigma$ (GPa)", size=34) plt.title("Off-Diagonal Stresses", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_stress_o.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_stress_o.png") #STRESSES GEO plt.figure(figsize=(15, 11)) index = lookup(output_vals_poro[0]) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , vec[:,index] ) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=34) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=34) plt.ylabel(r"$\phi$", size=34) plt.title("Porosity", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_geo_phi.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_geo_phi.png") plt.figure(figsize=(15, 11)) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , 1000000*vec[:,lookup(output_vals_w[0])], "rs", label= loutput_vals_w[0]) plt.plot( xcomp , 1000000*vec[:,lookup(output_vals_w[1])], "g.-", label= loutput_vals_w[1]) plt.plot( xcomp , 1000000*vec[:,lookup(output_vals_w[2])], "b.--", label= loutput_vals_w[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=34) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=34) plt.ylabel(r"$\bar{w}_{1}$ (nm)", size=34) plt.title("Crack aperture", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_geo_w.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_geo_w.png") plt.figure(figsize=(15, 11)) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_length[0])], "rs", label= loutput_vals_c_length[0]) plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_length[1])], "g.-", label= loutput_vals_c_length[1]) plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_length[2])], "b.--", label= loutput_vals_c_length[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=34) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=34) plt.ylabel(r"$\bar{c}_{1}$ ($\mu$m)", size=34) plt.title("Crack length", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_geo_c.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_geo_c.png") plt.figure(figsize=(15, 11)) xcomp = 1000*abs(vec[:,lookup(output_vals_main_stress[2])]) if int(shear_flag) == 1: xcomp = 1000*vec[:,lookup(output_vals_dev_stress[0])] plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_dens[0])], "ms" ,label= loutput_vals_c_dens[0]) plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_dens[1])], "y.-", label= loutput_vals_c_dens[1]) plt.plot( xcomp , 1000*vec[:,lookup(output_vals_c_dens[2])], "c.--", label= loutput_vals_c_dens[2]) plt.xlabel(r"$\sigma_{33}$ (GPa)", size=34) if int(shear_flag) == 1: plt.xlabel(r"$\sigma_{12}$ (GPa)", size=34) plt.ylabel(r"$\bar{l}_{1}$ ($\mu$m)", size=34) plt.title("Crack distance", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_geo_l.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_geo_l.png") #GEO PERM kfac = 1000*(1e-6)*(1/9.869233e-13) plt.figure(figsize=(15, 11)) index = lookup(output_vals_poro[0]) xcomp = vec[:,index] plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mD)", size=34) plt.xlabel(r"$\phi$", size=34) plt.title("Permeab. (Main) vs. Porosity", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_k_phi.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_k_phi.png") plt.figure(figsize=(15, 11)) xcomp = 1000000*vec[:,lookup(output_vals_w[0])] plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mD)", size=34) plt.xlabel(r"$\bar{w}_{1}$ (nm)", size=34) plt.title("Permeab. (Main) vs. Aperture", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_k_w.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_k_w.png") plt.figure(figsize=(15, 11)) xcomp = 1000*vec[:,lookup(output_vals_c_length[0])] plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.ylabel(r"$k$ (mD)", size=34) plt.xlabel(r"$\bar{c}_{1}$ ($\mu$m)", size=34) plt.title("Permeab. (Main) vs. Length", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_k_c.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_k_c.png") plt.figure(figsize=(15, 11)) xcomp = 1000*vec[:,lookup(output_vals_c_dens[0])] plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[1])], "g.-", label= loutput_vals_permeability[1]) plt.plot( xcomp , kfac*vec[:,lookup(output_vals_permeability[2])], "b.--", label= loutput_vals_permeability[2]) plt.gca().invert_xaxis() plt.ylabel(r"$k$ (mD)", size=34) plt.xlabel(r"$\bar{l}_{1}$ ($\mu$m)", size=34) plt.title("Permeab. (Main) vs. Distance", size=34) plt.legend(fontsize = 34, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.savefig("./plot_tests/"+ str(caseName) + "_k_l.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "_k_l.png") #Compare with analytical relations #Useful info for models k0 = vec[:,lookup(output_vals_permeability[0])][0] index = lookup(output_vals_poro[0]) poroV = vec[:,index] poro0 = vec[:,index][0] index = lookup(output_vals_p[0]) pres = vec[:,index] #Other constants dp = 0.1 alpha = 3 sigma_o = 1 kc = k0 sigmaB = 0.2 sigmaR = 2 gamma = 1 epstrain = 0.0001 #KC original (only uses porosity) def KC_original(poro): k = (poro**2 * dp**2) / (180 * (1 - poro)**2) return k #KC simple (only uses porosity) def KC_simple(poro): k = k0 * (poro/poro0)**alpha return k #Pressure related (only uses effective stress) def K_pres(pres): k = k0 * math.exp(-pres/sigma_o) return k #Pressure and KC (uses both porosity and effective stress) def KC_pres(poro, pres): k = k0 * (poro/poro0)**alpha * math.exp(-pres/sigma_o) return k #Pressure, KC and percolation (sets thresholds based on material properties based on strain) def KC_perc(pres, strain, sigma1, sigma3): if sigma1 - sigma3 < sigmaB: k = kc elif (sigma1 - sigma3 >= sigmaB) and (sigma1 - sigma3 <= sigmaR): k = k0 * math.exp(-gamma*pres) * (1 - math.exp(-strain/epstrain)) else: k = D*strain**3 return k KC_simpleV = np.zeros(len(time)) KC_originalV = np.zeros(len(time)) K_presV = np.zeros(len(time)) KC_presV = np.zeros(len(time)) KC_percV = np.zeros(len(time)) for i in range(len(time)): KC_simpleV[i] = KC_original(poroV[i]) KC_originalV[i] = KC_simple(poroV[i]) K_presV[i] = K_pres(pres[i]) KC_presV[i] = KC_pres(poroV[i], pres[i]) KC_percV[i] = KC_perc(pres[i], vec[:,lookup(output_vals_strainEv[0])][i], vec[:,lookup(output_vals_dev_stress[0])][i], pres[i]) plt.figure(figsize=(14, 11)) #1 mm2 = 9.86923e-10 mD #1 m2 = 9.86923e-16 mD mDconv = 9.86923e-10 plt.plot( time , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0] + " Sim.") plt.plot( time , KC_simpleV, "b-" , label= "KC Simple") plt.plot( time , KC_originalV, "g-" , label= "KC Original") plt.plot( time , K_presV, "k-" , label= "Pressure related") plt.plot( time , KC_presV, "m--" , label= "KC and Pressure related") plt.plot( time , KC_percV, "c--" , label= "KC, pres. and percolation related") plt.ylabel(r"$k$ ($mm^2$)", size=34) plt.xlabel(r"Time ($\mu$s)", size=34) plt.title("Permeab. Simulation vs. Analytical", size=34) plt.legend(fontsize = 20, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.yscale("log") #plt.show() plt.savefig("./plot_tests/"+ str(caseName) + "__SimvsAnalysis.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "__SimvsAnalysis.png") #Compare with data for rock salt plt.figure(figsize=(14, 11)) strainP40T24V = [1.45,2.91,3.92,5.01,5.79,7.11,7.54,9.19,9.7,11.83] strainP40T48V = [0.73,1.15,1.55,2.4,3.27,3.65,4.7,5.7,6.06,7.27,8.56,9.84,10.2,10.66,11.69,12.66,13.04] P40T24V = [0.0000000000183,0.0000000000547,0.0000000000998,0.000000000143,0.000000000141,0.000000000143,0.000000000144,0.000000000139,0.000000000131,0.000000000115] P40T48V = [0.000000000000248,0.0000000000155,0.0000000000224,0.000000000108,0.000000000277,0.000000000309,0.000000000386,0.000000000403,0.000000000396,0.000000000375,0.000000000339,0.000000000299,0.000000000293,0.000000000276,0.000000000249,0.000000000223,0.000000000217] plt.plot( strainP40T24V , P40T24V, "r-*" , label= "Sample P40T24") plt.plot( strainP40T48V , P40T48V, "b-*" , label= "Sample P40T48") plt.plot( abs(vec[:,lookup(output_vals_strainEv[0])]*100), vec[:,lookup(output_vals_permeability[0])], "ks-" , label= output_vals_strainEv[0] +" strain effect sim") #plt.plot( time , vec[:,lookup(output_vals_permeability[0])], "rs" , label= loutput_vals_permeability[0] + " Sim.") # plt.plot( time , KC_simpleV, "b-" , label= "KC Simple") # plt.plot( time , KC_originalV, "g-" , label= "KC Original") # plt.plot( time , K_presV, "k-" , label= "Pressure related") # plt.plot( time , KC_presV, "p-" , label= "KC and Pressure related") #plt.plot( time , KC_percV, "c-" , label= "KC, pres. and percolation related") plt.ylabel(r"$k$ ($mm^2$)", size=34) plt.xlabel(r"% strain", size=34) plt.title("Permeab. vs. Strain for Samples (Peach, 1996)", size=34) plt.legend(fontsize = 20, framealpha=0.2) plt.tick_params(labelsize=26) plt.xticks(fontsize=26) plt.yticks(fontsize=26) plt.yscale("log") #plt.xscale("log") #plt.show() plt.savefig("./plot_tests/"+ str(caseName) + "__rocksalt_data.png") plt.savefig("/g/g91/moncadalopez1/ktests/"+ str(caseName) + "__rocksalt_data.png")