Skip to content

Commit

Permalink
corrigido mic + melhorado analise
Browse files Browse the repository at this point in the history
  • Loading branch information
g7fernandes committed Sep 19, 2019
1 parent 9b5dcf3 commit 7562d1d
Show file tree
Hide file tree
Showing 4 changed files with 126 additions and 56 deletions.
100 changes: 60 additions & 40 deletions analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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...")
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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')
Expand All @@ -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()
22 changes: 14 additions & 8 deletions lennard.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)]
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -6340,19 +6343,22 @@ 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)]
end do
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), &
Expand Down
2 changes: 1 addition & 1 deletion saida.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
58 changes: 51 additions & 7 deletions tecplot_lines.py
Original file line number Diff line number Diff line change
@@ -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)
# # clear plot
# plot.delete_linemaps()

0 comments on commit 7562d1d

Please sign in to comment.