diff --git a/analysis.py b/analysis.py index 1d3b7fe..cca9bf1 100644 --- a/analysis.py +++ b/analysis.py @@ -82,7 +82,7 @@ def density_map(x, dimx ,div): tauyxp = np.zeros((len_list_files,div)) nesb = int(input("Enter the number of assembles in which this simulation will be divided:\n")) -periodic = ('y' == input("Where the boundaries periodic? [y/n]\n")) +periodic = ('y' == input("Were the boundaries periodic? [y/n]\n")) esb_len = int(len_list_files/nesb) # length of each ensamble mic = 0 print("Reading files...") @@ -104,13 +104,13 @@ def density_map(x, dimx ,div): r0 = pos[['x','y']] # um r0 por ensamble lp = dmap # salva as partículas numa daterminada região # if periodic: - # mic = rfup[["micx", "micy"]].to_numpy()*np.array([dimx, dimy]) + # mic = rfup[["micx", "micy"]].to_numpy()*np.array([dimx, dimy]) # Separar aqui por região for i in range(div): #mesmas particulas - # msq[step,i] = np.mean(np.sum(np.square(pos.loc[lp[i],['x','y']] + mic.loc[lp[i]] - r0.loc[lp[i]]),axis=1)) + # msq[step,i] = np.mean(np.sum(np.square(pos.loc[lp[i],['x','y']] + mic[lp[i],:] - r0.loc[lp[i]]),axis=1)) # diferentes particulas KE[step,i] = np.mean(rfup.loc[dmap[i],'K']) tauk[step,i] = np.sum(vel.loc[dmap[i],'v_x'] * vel.loc[dmap[i],'v_y'] * rfup.loc[dmap[i],'m'])/Vol @@ -125,9 +125,21 @@ def density_map(x, dimx ,div): # Calculamos t * viscosidade*2*kT/Vol. kT = KE -etat_kk, etat_kp1, etat_pp1, etat_kp2, etat_pp2 = np.zeros(esb_len) +etat_kk = np.zeros(esb_len) +etat_kp1 = np.zeros(esb_len) +etat_pp1 = np.zeros(esb_len) +etat_kp2 = np.zeros(esb_len) +etat_pp2 = np.zeros(esb_len) + +if (div > 1): + overdiv = np.zeros((div,2)) for i in range(div): + etat_kk = 0*etat_kk + etat_kp1 = 0*etat_kp1 + etat_pp1 = 0*etat_pp1 + etat_kp2 = 0*etat_kp2 + etat_pp2 = 0*etat_pp2 for ii in range(1,esb_len): # ao longo de cada ensemble for j in range(nesb): # ao longo dos ensambles start = j*esb_len @@ -140,27 +152,28 @@ def density_map(x, dimx ,div): # Viscosidade * t - etat_kk = (etat_kk * Vol / (2*kT))/nesb - etat_kp1 = (etat_kp1 * Vol / (kT) )/nesb - etat_pp1 = (etat_pp1 * Vol / (2*kT))/nesb - etat_kp2 = (etat_kp2 * Vol / (kT) )/nesb - etat_pp2 = (etat_pp2 * Vol / (2*kT))/nesb + etat_kk = (etat_kk * Vol / (2*kT[i]))/nesb + etat_kp1 = (etat_kp1 * Vol / (kT[i]) )/nesb + etat_pp1 = (etat_pp1 * Vol / (2*kT[i]))/nesb + etat_kp2 = (etat_kp2 * Vol / (kT[i]) )/nesb + etat_pp2 = (etat_pp2 * Vol / (2*kT[i]))/nesb - # Para descobrir fazemos um curve fit linear + # Para descobrir fazemos um curve fit linear no último terço de pontos t = np.linspace(0,esb_len-1, esb_len) - eta_kk = np.polyfit(t,etat_kk,1) - eta_kp1 = np.polyfit(t,etat_kp1,1) - eta_pp1 = np.polyfit(t,etat_pp1,1) - eta_kp2 = np.polyfit(t,etat_kp2,1) - eta_pp2 = np.polyfit(t,etat_pp2,1) + eta_kk = np.polyfit(t[int(esb_len/2):esb_len],etat_kk[int(esb_len/2):esb_len],1) + eta_kp1 = np.polyfit(t[int(esb_len/2):esb_len],etat_kp1[int(esb_len/2):esb_len],1) + eta_pp1 = np.polyfit(t[int(esb_len/2):esb_len],etat_pp1[int(esb_len/2):esb_len],1) + eta_kp2 = np.polyfit(t[int(esb_len/2):esb_len],etat_kp2[int(esb_len/2):esb_len],1) + eta_pp2 = np.polyfit(t[int(esb_len/2):esb_len],etat_pp2[int(esb_len/2):esb_len],1) # # Calcula a Difusividade # D = np.polyfit(t,msq,1) - # print("Region: {}\n".format(i)) - # fig, axs = plt.subplot(1,3) - # # Plota a difusividade + print("Region: {}\n".format(i)) + fig, axs = plt.subplots(1,3) # axs[1,0].figure() + + # # Plota a difusividade # axs[1,0].scatter(t,msq[:,i]) # axs[1,0].plot(t,t*D[0] + D[1],'k') # axs[1,0].title('Mean square displacement') @@ -173,42 +186,49 @@ def density_map(x, dimx ,div): # Plota as viscosidades usando taupxy #plt.figure() - axs[1,1].scatter(t,etat_kk,c='g',label='t*eta_kk') - axs[1,1].scatter(t,etat_kp1,c='b',label='t*eta_kp_xy') - axs[1,1].scatter(t,etat_pp1,c='r',label='t*eta_pp_xy') + axs[1].scatter(t,etat_kk,c='g',marker='.',label='t*eta_kk') + axs[1].scatter(t,etat_kp1,c='b',marker='.',label='t*eta_kp_xy') + axs[1].scatter(t,etat_pp1,c='r',marker='.',label='t*eta_pp_xy') - axs[1,1].plot(t, t*eta_kk[0] + etat_kk[1],'g') - axs[1,1].plot(t, t*eta_kp1[0] + etat_kp1[1],'b') - axs[1,1].plot(t, t*eta_pp1[0] + etat_pp1[1],'r') + axs[1].plot(t, t*eta_kk[0] + etat_kk[1],'g') + axs[1].plot(t, t*eta_kp1[0] + etat_kp1[1],'b') + axs[1].plot(t, t*eta_pp1[0] + etat_pp1[1],'r') - axs[1,1].plot(t, t*(eta_kk[0] + eta_kp1[0] + eta_pp1[0] ) + etat_kk[1] + etat_kp1[1] + etat_pp1[1],'k',linewidth=2, label='t*eta') - axs[1,1].legend() - axs[1,1].title("Viscosity using tau_yx") + axs[1].plot(t, t*(eta_kk[0] + eta_kp1[0] + eta_pp1[0] ) + etat_kk[1] + etat_kp1[1] + etat_pp1[1],'k',linewidth=2, label='t*eta') + axs[1].legend() + axs[1].set_title("Viscosity using tau_yx") props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) - string = "etaxy = {}".format(eta_kk[0]+eta_kp1[0]+eta_pp1[0]) - axs[1,1].text(0.95,0.05,string, transform=axs[1,1].transAxes, fontsize=10, verticalalignment='bottom', bbox=props) + string = "etaxy = {:.3}".format(eta_kk[0]+eta_kp1[0]+eta_pp1[0]) + overdiv[i,0] = eta_kk[0]+eta_kp1[0]+eta_pp1[0] + axs[1].text(0.05,0.65,string, transform=axs[1].transAxes, fontsize=10, verticalalignment='bottom', bbox=props) # Plota as viscosidades ustando taupyx #plt.figure() - axs[1,2].scatter(t,etat_kk,c='g',label='t*eta_kk') - axs[1,2].scatter(t,etat_kp2,c='b',label='t*eta_kp_yx') - axs[1,2].scatter(t,etat_pp2,c='r',label='t*eta_pp_yx') + axs[2].scatter(t,etat_kk,c='g',marker='.',label='t*eta_kk') + axs[2].scatter(t,etat_kp2,c='b',marker='.',label='t*eta_kp_yx') + axs[2].scatter(t,etat_pp2,c='r',marker='.',label='t*eta_pp_yx') - axs[1,2].plot(t, t*eta_kk[0] + etat_kk[1],'g') - axs[1,2].plot(t, t*eta_kp2[0] + etat_kp2[1],'b') - axs[1,2].plot(t, t*eta_pp2[0] + etat_pp2[1],'r') + axs[2].plot(t, t*eta_kk[0] + etat_kk[1],'g') + axs[2].plot(t, t*eta_kp2[0] + etat_kp2[1],'b') + axs[2].plot(t, t*eta_pp2[0] + etat_pp2[1],'r') - axs[1,2].plot(t, t*(eta_kk[0] + eta_kp2[0] + eta_pp2[0] ) + etat_kk[1] + etat_kp2[1] + etat_pp2[1],'k',linewidth=2, label='t*eta') - axs[1,2].legend() - axs[1,2].title("Viscosity using tau_xy.") + axs[2].plot(t, t*(eta_kk[0] + eta_kp2[0] + eta_pp2[0] ) + etat_kk[1] + etat_kp2[1] + etat_pp2[1],'k',linewidth=2, label='t*eta') + axs[2].legend() + axs[2].set_title("Viscosity using tau_xy.") props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) - string = "etayx = {}".format(eta_kk[0]+eta_kp2[0]+eta_pp2[0]) - axs[1,1].text(0.95,0.05,string, transform=axs[1,2].transAxes, fontsize=10, verticalalignment='bottom', bbox=props) + string = "etayx = {:.3f}".format(eta_kk[0]+eta_kp2[0]+eta_pp2[0]) + overdiv[i,1] = eta_kk[0]+eta_kp2[0]+eta_pp2[0] + axs[2].text(0.05,0.65,string, transform=axs[2].transAxes, fontsize=10, verticalalignment='bottom', bbox=props) print("\nCom tau_xy: \neta_kk = {}\neta_kp = {}\neta_pp = {}\neta = {}\n".format(eta_kk[0],eta_kp1[0],eta_pp1[0], eta_kk[0]+eta_kp1[0]+eta_pp1[0])) print("\nCom tau_x]yx: \neta_kk = {}\neta_kp = {}\neta_pp = {}\neta = {}\n".format(eta_kk[0],eta_kp2[0],eta_pp2[0], eta_kk[0]+eta_kp2[0]+eta_pp2[0])) fig.suptitle("Region {}".format(i)) +plt.figure() +plt.plot(overdiv[:,0],label='eta_xy') +plt.plot(overdiv[:,1],label='eta_yx') +plt.title('Over the positions') + plt.show() \ No newline at end of file diff --git a/lennard.f90 b/lennard.f90 index 983f2c8..ae8a90a 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -6107,6 +6107,7 @@ program main if (id == 1) then tag = 10 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 10 @@ -6119,6 +6120,7 @@ program main if (id == 2) then tag = 20 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 20 @@ -6130,6 +6132,7 @@ program main if (id == 3) then tag = 30 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 30 @@ -6147,7 +6150,7 @@ program main [real(mic(ii+1,1),kind(0.d0)),real(mic(ii+1,2),kind(0.d0)),nRfu(ii*6+2), & nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)] end do - call vec2csv(rFUp,N,7,'rF_u_P',j,t,nimpre,start) + call vec2csv(rFUp,N,6,'rF_u_P',j,t,nimpre,start) else do ii = 0, N-1 rFUp(int(nRfu(ii*6+1)),:) = [nRfu(ii*6+2),nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)] @@ -6195,12 +6198,10 @@ program main ! COMP X call comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,domy,wall,mic,id,np,xih,xic,cold_cells,hot_cells) ! altera posição - !print*, "L 2002", id - ! if (id == 0) read(*,*) - ! call MPI_barrier(MPI_COMM_WORLD, ierr) + call walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ! altera posição e malha - + ! COMP V call comp_v_thermo_pred(malha,mesh,dt,t,np,propriedade,domx,domy, Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells) @@ -6257,8 +6258,7 @@ program main call vec2csv(v,N,2,'velocity',j,t,nimpre,start) end if - ! if (id == 0) read(*,*) - ! call MPI_barrier(MPI_COMM_WORLD, ierr) + deallocate(nxv,nxv_send) @@ -6305,6 +6305,7 @@ program main if (id == 1) then tag = 10 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 10 @@ -6317,6 +6318,7 @@ program main if (id == 2) then tag = 20 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 20 @@ -6328,6 +6330,7 @@ program main if (id == 3) then tag = 30 call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + mic = 0 end if if (id == 0) then tag = 30 @@ -6340,12 +6343,14 @@ program main if (id == 0) then if (wall(1:2) == 'pp' .or. wall(3:4) == 'pp') then + + do ii = 0, N-1 rFUp(int(nRfu(ii*6+1)),:) = & [real(mic(ii+1,1),kind(0.d0)),real(mic(ii+1,2),kind(0.d0)),nRfu(ii*6+2), & nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)] end do - call vec2csv(rFUp,N,7,'rF_u_P',j,t,nimpre,start) + call vec2csv(rFUp,N,6,'rF_u_P',j,t,nimpre,start) else do ii = 0, N-1 rFUp(int(nRfu(ii*6+1)),:) = [nRfu(ii*6+2),nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)] @@ -6353,6 +6358,7 @@ program main call vec2csv(rFUp,N,5,'rF_u_P',j,t,nimpre,start) end if end if + call MPI_barrier(MPI_COMM_WORLD, ierr) ! deallocate(nRfu,nRfu_send,rFUp) end if allocate(LT%lstrdb_N(NMPT*6),LT%lstrdb_S(NMPT*6),LT%lstrdb_E(NMPT*6),LT%lstrdb_W(NMPT*6),LT%lstrdb_D(NMPT*6), & diff --git a/saida.f90 b/saida.f90 index f91701c..c311d88 100644 --- a/saida.f90 +++ b/saida.f90 @@ -18,7 +18,7 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start,concat0) real(dp), save :: timep, etc, dtimepp character(LEN=*),parameter :: fmt5 = '(f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )' character(LEN=*),parameter :: fmt6 = '(i10, ", ",i10, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )' - character(LEN=*),parameter :: fmt7 = '(i10, ", ",i10, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )' + character(LEN=*),parameter :: fmt7 = '(f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )' laux = .false. if (present(concat0)) laux = concat0 diff --git a/tecplot_lines.py b/tecplot_lines.py index a25449d..ce3f813 100644 --- a/tecplot_lines.py +++ b/tecplot_lines.py @@ -1,36 +1,80 @@ import numpy as np import tecplot as tp import pandas as pd -from tecplot.constant import PlotType, Color +from tecplot.constant import PlotType, Color, LinePattern, AxisTitleMode import sys -if '-c' in sys.argv: - tp.session.connect() +# if '-c' in sys.argv: +# tp.session.connect() + +tp.session.connect() # IMPORTA LINHAS DO CSV csv_input = pd.read_csv('Distributions.csv') hl = list(csv_input) #headers list +for i in range(len(hl)): + print("{} | {}\n".format(i,hl[i])) +linhas = input("Entre os numeros dos dados que quer plotar separados por espaço:\n").split() +linhas = [int(x) for x in linhas] + frame = tp.active_frame() frame.activate() frame.height = 8 frame.width = 11 dataset = frame.create_dataset("Data", ['x', 'y']) -for col in hl: +for dado in linhas: + col = hl[dado] x = csv_input[col].to_numpy() zone = dataset.add_ordered_zone(col, len(x)) zone.values('x')[:] = np.linspace(0,1,len(x)) zone.values('y')[:] = x - # Set plot type to XYLine plot = frame.plot(PlotType.XYLine) plot.activate() +############# EIXOS ############################# +x_axis = plot.axes.y_axis(0) +x_axis.title.title_mode = AxisTitleMode.UseText +x_axis.title.text = 'x' +x_axis.fit_range_to_nice() + +y_axis = plot.axes.y_axis(0) +y_axis.title.title_mode = AxisTitleMode.UseText +y_axis.title.text = 'Fraction of Particles' +y_axis.fit_range_to_nice() + +for ax in [x_axis, y_axis]: + ax.title.font.typeface = 'Times' + ax.title.font.bold = False + ax.title.font.italic = False + ax.title.font.size_units = Units.Frame + ax.title.font.size = 7 + +################ LINHAS ##################### + +# Line patters list +# DashDot = 2 DashDotDot = 5 +# Dashed = 1 Dotted = 3 +# LongDash = 4 Solid = 0 + # Show all linemaps and make the lines a bit thicker +cont = 1 for lmap in plot.linemaps(): lmap.show = True lmap.line.line_thickness = 0.6 + lmap.line.line_pattern = LinePattern(cont) + lmap.line.color = Color.Black + cont += 1 + +################# LEGENDA ##################### +legend = plot.legend +legend.font.typeface = 'Times' +legend.show = True + -plot.legend.show = True +################ IMAGEM SALVA ###################### +tp.export.save_png('add_ordered_zones.png', 600, supersample=3) -tp.export.save_png('add_ordered_zones.png', 600, supersample=3) \ No newline at end of file +# # clear plot +# plot.delete_linemaps() \ No newline at end of file