From 9b5dcf3c0f9fabeebb91e90bd1094fa79c1f7057 Mon Sep 17 00:00:00 2001 From: g7fernandes Date: Wed, 18 Sep 2019 12:03:43 -0300 Subject: [PATCH] =?UTF-8?q?corre=C3=A7=C3=B5es=20analysis=20+=20mic=20do?= =?UTF-8?q?=20rfup?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- analysis.py | 64 ++++++++++++++++++++++++++++++----------------------- lennard.f90 | 12 +++++----- saida.f90 | 6 ++--- 3 files changed, 45 insertions(+), 37 deletions(-) diff --git a/analysis.py b/analysis.py index db793ae..1d3b7fe 100644 --- a/analysis.py +++ b/analysis.py @@ -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 @@ -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"]) @@ -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) @@ -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 diff --git a/lennard.f90 b/lennard.f90 index cda11d4..983f2c8 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -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 @@ -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 @@ -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 !!! @@ -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) @@ -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) @@ -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) diff --git a/saida.f90 b/saida.f90 index c025c7b..f91701c 100644 --- a/saida.f90 +++ b/saida.f90 @@ -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 @@ -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