Skip to content

Commit

Permalink
correções analysis + mic do rfup
Browse files Browse the repository at this point in the history
  • Loading branch information
@g7fernandes
g7fernandes committed Sep 18, 2019
1 parent 91d132e commit 9b5dcf3
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 37 deletions.
64 changes: 36 additions & 28 deletions analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@ def density_map(x, dimx ,div):
# dimx: dimension of the region in the diraction of the divison
# div: number of divisions
# outut: y: list of arrays that contains the indices of the particles that are in a region that is the position in the array
div = dimx/div
div = int(dimx/div)
xp = x // div # This gives the location of each particle
dmap = [[] for _ in range(div)]
for i in range(len(x)):
dmap[xp].append(i)
dmap[int(xp[i])].append(i)
return dmap

# Find the directory where the files are
Expand Down Expand Up @@ -75,13 +75,20 @@ def density_map(x, dimx ,div):
else:
div = int(div)

KE,tauk,tauxyp,msq,tauyxp = np.zeros((len_list_files,div))
KE = np.zeros((len_list_files,div))
tauk = np.zeros((len_list_files,div))
tauxyp = np.zeros((len_list_files,div))
msq = np.zeros((len_list_files,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"))
esb_len = int(len_list_files)/nesb # length of each ensamble
esb_len = int(len_list_files/nesb) # length of each ensamble
mic = 0
for step in range(len_list_files):
print("Reading files...")
if pbar:
bar = progressbar.ProgressBar(max_value=(len_list_files))
for step in range(1,len_list_files):
esb = step//esb_len
if periodic:
rfup = pd.read_csv(zip_rfup.open('rF_u_P.csv.'+str(step)), header=None, names = ["micx","micy","RxFy","RyFx","u","m","K"])
Expand All @@ -93,24 +100,25 @@ def density_map(x, dimx ,div):
dmap = density_map(pos['x'],dimx,div)

# Depende da região, terá um for
if step%esb_len == 0: # se step == inicio do ensamble
if step%esb_len == 1: # se step == inicio do ensamble
r0 = pos[['x','y']] # um r0 por ensamble
lp = dmap # salva as partículas numa daterminada região
if periodic:
mic = rfup[["micx", "micy"]]*np.array([dimx, dimy])
# if periodic:
# mic = rfup[["micx", "micy"]].to_numpy()*np.array([dimx, dimy])

# Separar aqui por região

for i in range(div):
#mesmas particulas
msq[step,div] = 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.loc[lp[i]] - r0.loc[lp[i]]),axis=1))
# diferentes particulas
KE[step,div] = np.mean(rfup.loc[dmap[i],'K'])
tauk[step,div] = np.sum(vel.loc[dmap[i],'v_x'] * vel.loc[dmap[i],'v_y'] * rfup.loc[dmap[i],'m'])/Vol
tauxyp[step,div] = np.sum(rfup.loc[dmap[i],'RxFy'])/Vol
tauyxp[step,div] = np.sum(rfup.loc[dmap[i],'RyFx'])/Vol

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
tauxyp[step,i] = np.sum(rfup.loc[dmap[i],'RxFy'])/Vol
tauyxp[step,i] = np.sum(rfup.loc[dmap[i],'RyFx'])/Vol

if pbar: #atualiza a barra de progresso
bar.update(step)

kT = KE.mean(axis=0)

Expand Down Expand Up @@ -146,20 +154,20 @@ def density_map(x, dimx ,div):
eta_kp2 = np.polyfit(t,etat_kp2,1)
eta_pp2 = np.polyfit(t,etat_pp2,1)

# Calcula a Difusividade
D = np.polyfit(t,msq,1)

print("Region: {}\n".format(i))
fig, axs = plt.subplot(1,3)
# Plota a difusividade
axs[1,0].figure()
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')
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
string = "D = {}".format(D[0]/4)
print(string + "\n")
axs[1,0].text(0.95,0.05,string, transform=axs[1,0].transAxes, fontsize=10, verticalalignment='bottom', bbox=props)
# # Calcula a Difusividade
# D = np.polyfit(t,msq,1)

# print("Region: {}\n".format(i))
# fig, axs = plt.subplot(1,3)
# # Plota a difusividade
# axs[1,0].figure()
# 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')
# props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
# string = "D = {}".format(D[0]/4)
# print(string + "\n")
# axs[1,0].text(0.95,0.05,string, transform=axs[1,0].transAxes, fontsize=10, verticalalignment='bottom', bbox=props)


# Plota as viscosidades usando taupxy
Expand Down
12 changes: 6 additions & 6 deletions lennard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -476,7 +476,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu)
end do
end do
end do
aux1 = aux1 - 1

! print*, "ID, AUX1", id, aux1
! print*, "FIM"
end subroutine comp_pot
Expand Down Expand Up @@ -1921,7 +1921,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,

cont_db(1) = cont_db(1) + 6
cont_int(1) = cont_int(1) + 4
mic(ptr%p%n, 2) = mic(ptr%p%n, 2) - 1
mic(ptr%p%n, 2) = mic(ptr%p%n, 2) + 1
end if
if (south == 'p' .and. cell(1) == 1 .and. (ids(1) + ids(2)) /= -2) then
! print*, "L 580", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j
Expand All @@ -1931,7 +1931,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
! print*, "> id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo]
cont_db(2) = cont_db(2) + 6
cont_int(2) = cont_int(2) + 4
mic(ptr%p%n, 2) = mic(ptr%p%n, 2) + 1
mic(ptr%p%n, 2) = mic(ptr%p%n, 2) - 1
end if
end if
!!! FIM DA VERIFICAÇÃO SE MUDOU DE DOMÍNIO !!!
Expand Down Expand Up @@ -6144,7 +6144,7 @@ program main
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),kind(0.d0)),real(mic(ii,2),kind(0.d0)),nRfu(ii*6+2), &
[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)
Expand Down Expand Up @@ -6342,7 +6342,7 @@ program main
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),kind(0.d0)),real(mic(ii,2),kind(0.d0)),nRfu(ii*6+2), &
[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)
Expand Down Expand Up @@ -6567,7 +6567,7 @@ program main
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),kind(0.d0)),real(mic(ii,2),kind(0.d0)),nRfu(ii*6+2), &
[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,6,'rF_u_P',j,t,nimpre,start)
Expand Down
6 changes: 3 additions & 3 deletions saida.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,8 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start,concat0)
real(dp) :: time
real(dp), save :: timep, etc, dtimepp
character(LEN=*),parameter :: fmt5 = '(f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )'
character(LEN=*),parameter :: fmt6 = '(f5.0, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )'
character(LEN=*),parameter :: fmt7 = '(f5.0, ", ",f5.0, ", ",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 )'

laux = .false.
if (present(concat0)) laux = concat0
Expand Down Expand Up @@ -52,7 +52,7 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start,concat0)
end do
else if (d == 6) then
do i = 1,n
write(10,fmt6) v(i,1),v(i,2),v(i,3),v(i,4),v(i,5),v(i,6)
write(10,fmt6) int(v(i,1)),int(v(i,2)),v(i,3),v(i,4),v(i,5),v(i,6)
end do
else if (d == 7) then
do i = 1,n
Expand Down

0 comments on commit 9b5dcf3

Please sign in to comment.