From ed1df7cdc5a2ec557ce44dbf7068a6114a42d3dc Mon Sep 17 00:00:00 2001 From: Gabriel Fernandes Date: Wed, 28 Aug 2019 22:05:58 -0300 Subject: [PATCH] =?UTF-8?q?em=20implementa=C3=A7=C3=A3o=20grad=20for=C3=A7?= =?UTF-8?q?ado?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- csv2vtk_particles.py | 2 +- lennard.f90 | 147 ++++++++++++---- lennard_jones_stability_analysis.py | 188 ++++++++++++++++++++ part_gen_pequeno.py | 4 +- potential_analysis.py | 264 ++++++++++++++-------------- saida.f90 | 27 ++- settings.ini | 61 +++---- 7 files changed, 490 insertions(+), 203 deletions(-) create mode 100644 lennard_jones_stability_analysis.py diff --git a/csv2vtk_particles.py b/csv2vtk_particles.py index 2236131..22920d4 100755 --- a/csv2vtk_particles.py +++ b/csv2vtk_particles.py @@ -56,7 +56,7 @@ opt = input('Overwrite? [y/n] ') else: opt = input('Append? (starting on {}) [y/n] '.format(nimpre_init)) - if opt == 'y' or 'opt' == 'Y': + if opt == 'y' or opt == 'Y': aux =False else: folder = input('Enter new folder name:\n') diff --git a/lennard.f90 b/lennard.f90 index 786dcd7..04e2bc6 100755 --- a/lennard.f90 +++ b/lennard.f90 @@ -2681,7 +2681,7 @@ subroutine comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,d end if end if - + ptr%p%flag = .false. @@ -3451,7 +3451,7 @@ subroutine comp_v(malha,mesh,dt,t,propriedade,domx,domy) end do end subroutine comp_v - subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells, id) + subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells) use linkedlist use mod1 use data @@ -3467,7 +3467,7 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic type(data_ptr) :: ptr real(dp):: Kcold, Khot, numpc, numph, aux(4), auxres(16) ! letra grega xi, mas pra não com x confundir será ksi. Estes são os termos de amortecimento para termostato Nosé-Hoover. numpc e nump são o número de partículas nas regiões quentes e frias real(dp), intent(inout) :: xih, xic - integer ( kind = 4 ) :: np, ierr, id + integer ( kind = 4 ) :: ierr !, np, id numph = 0 numpc = 0 @@ -3531,6 +3531,59 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic end do end do + ! do i = hot_cells(3), hot_cells(4) !(domy(1)+bsh(3)), (domy(2)+bsh(4)) + ! do j = hot_cells(1), hot_cells(2) ! (domx(1)+bsh(1)), (domx(2)+bsh(2)) + ! !if (j >= hot_cells(1) .and. j <= hot_cells(2) .and. i >= hot_cells(3) .and. i <= hot_cells(4)) then + ! ! Termostato afeta QUENTE + ! node => list_next(malha(i,j)%list) + ! do while (associated(node)) + ! ptr = transfer(list_get(node), ptr) + ! if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + ! if (propriedade(ptr%p%grupo)%ismolecule) then + ! ! Termostato afeta + ! ! passo preditor + ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(1) + ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(2) + ! Khot = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Khot + ! numph = numph + 1 + ! else + ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) + ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) + + ! end if + ! end if + ! ptr%p%F = [0,0] + ! node => list_next(node) + ! end do + ! end do + ! end do + ! do i = cold_cells(3), cold_cells(4) !(domy(1)+bsh(3)), (domy(2)+bsh(4)) + ! do j = cold_cells(1), cold_cells(2) ! (domx(1)+bsh(1)), (domx(2)+bsh(2)) + ! ! Termostato afeta FRIO + ! ! else if (j >= cold_cells(1) .and. j <= cold_cells(2) .and. i >= cold_cells(3) .and. i <= cold_cells(4)) then + ! node => list_next(malha(i,j)%list) + ! do while (associated(node)) + ! ptr = transfer(list_get(node), ptr) + ! if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + ! if (propriedade(ptr%p%grupo)%ismolecule) then + ! ! Termostato afeta + ! ! passo preditor + ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(1)/2 + ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(2)/2 + ! Kcold = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Kcold + ! numpc = numpc + 1 + ! else + ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) + ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) + ! end if + ! end if + ! ptr%p%F = [0,0] + ! node => list_next(node) + ! end do + ! !end if + ! end do + ! end do + if (np > 1) then aux = [Khot, numph, Kcold, numpc] call MPI_ALLGATHER(aux, 4, MPI_DOUBLE_PRECISION, auxres, 4, MPI_DOUBLE_PRECISION, MPI_COMM_WORLD, ierr) @@ -3562,7 +3615,7 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic ! molecular simulation pag. 52, 101 do pdf end subroutine comp_v_thermo_pred - subroutine comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, hot_cells, cold_cells,id) + subroutine comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, hot_cells, cold_cells) use linkedlist use mod1 use data @@ -3571,7 +3624,7 @@ subroutine comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, h type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha integer :: i,j, bsh(4) - integer, intent(in) :: mesh(2),domx(2),domy(2), cold_cells(4), hot_cells(4), id + integer, intent(in) :: mesh(2),domx(2),domy(2), cold_cells(4), hot_cells(4) !,id type(list_t), pointer :: node real(dp), intent(in) :: dt, t type(data_ptr) :: ptr @@ -4075,12 +4128,12 @@ program main implicit none ! Variáveis integer :: N,Ntype,i=1, nimpre,j = 1, ii, quant = 0,mesh(2), cont = 1, aux1 = 0,cont2 = 1,domx(2), domy(2), aux3 - integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init + integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init, cpu_countrate,ic1, force_lingradT integer, target :: k real(dp), dimension(:,:), allocatable :: v, x, celula, rFUp !força n e n+1, velocidades e posições, [r*F, potencial e momento] real(dp), dimension(:), allocatable :: icell,jcell, nxv, nxv_send, nRfu, nRfu_send !dimensões das celulas e vetor de resultado pra imprimir - real(dp) :: t=0,t_fim,dt,printstep, sigma, epsil, rcut,aux2,start = 0,finish,Td,kb = 1.38064852E-23,vd(2) - real(dp) :: GField(2), temp_Td(3), dimX, dimY, dx_max, Td_hot, Td_cold, xih=0, xic=0, Mc + real(dp) :: t=0,t_fim,dt,printstep, sigma, epsil, rcut,aux2,start,finish,Td,kb = 1.38064852E-23,vd(2) + real(dp) :: GField(2), temp_Td(3), dimX, dimY, dx_max, Td_hot, Td_cold, xih, xic, Mc integer,allocatable :: interv(:), interv_Td(:), grupo(:), rcounts(:), displs(:), mic(:,:), mic_trf(:) ! mic são vetores para contar quantas vezes atravessou a borda periodica type(container), allocatable,dimension(:,:) :: malha type(prop_grupo),allocatable,dimension(:) :: propriedade @@ -4104,6 +4157,10 @@ program main integer ( kind = 4 ) id, ids(8) real(real32) :: nan + start = 0 + xih = 0 + xic = 0 + nan = IEEE_VALUE(nan, IEEE_QUIET_NAN) @@ -4161,6 +4218,8 @@ program main "Thermostat hot wall") call CFG_add(my_cfg,"global%Td_cold",-1.0_dp, & "Thermostat cold wall") + call CFG_add(my_cfg,"global%force_lingradT",0, & + "Force Temperature linear temperature gradient in the X diraction") call CFG_add(my_cfg,"global%cold_cells",(/0, 0, 0, 0/), & "Cold Cells") call CFG_add(my_cfg,"global%hot_cells",(/0, 0, 0, 0/), & @@ -4195,12 +4254,14 @@ program main call CFG_get(my_cfg,"global%hot_cells",hot_cells) call CFG_get(my_cfg,"global%Td_cold",Td_cold) call CFG_get(my_cfg,"global%Td_hot",Td_hot) + call CFG_get(my_cfg,"global%force_lingradT",force_lingradT) call CFG_get(my_cfg,"global%vd",vd) call CFG_get(my_cfg,"global%NMPT",NMPT) call CFG_get(my_cfg,"global%GField",GField) call CFG_get(my_cfg,"global%print_TC",print_TC) + if (NMPT < 1) NMPT = N printstep = ((temp_Td(2)-temp_Td(1))/dt)/temp_Td(3) if (printstep < 1) printstep = 1 @@ -4377,17 +4438,7 @@ program main ! imprime condições iniciais - if (id == 0) then - j = 0 - call system('mkdir temp') !pasta temporária para armazenar os resultados - ! call linked2vec(malha,domx,domy,nxv,aux1) - call vec2csv(x,N,2,'position',j,t,nimpre,start) - call vec2csv(v,N,2,'velocity',j,t,nimpre,start) - ! call vec2csv(celula,N,2,'cell',j) - j = j + 1 - - ! print*,'V_tot0 = ',comp_pot(mesh,malha,propriedade,rcut) - end if + !! AQUI COMEÇA O PARALELISMO @@ -4473,6 +4524,30 @@ program main print '("id = ", I2, " ids = [", I2, " ", I2, " ", I2, " ", I2,"]")', id, ids(1), ids(2), ids(3), ids(4) ! libera a memória utilizada pelos processos não-raiz ! pois eles não vão imprimir x e v pra csv + + ! id = 0 é o processo raiz + + ! if (id == 0) then + ! print*,'Press Return to continue' + ! read(*,*) + ! end if + + if (id == 0) then + call system_clock(ic1,cpu_countrate) + start = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0)) + end if + + if (id == 0) then + j = 0 + call system('mkdir temp') !pasta temporária para armazenar os resultados + ! call linked2vec(malha,domx,domy,nxv,aux1) + call vec2csv(x,N,2,'position',j,t,nimpre,start) + call vec2csv(v,N,2,'velocity',j,t,nimpre,start) + ! call vec2csv(celula,N,2,'cell',j) + j = j + 1 + ! print*,'V_tot0 = ',comp_pot(mesh,malha,propriedade,rcut) + end if + if (id /= 0) then deallocate(v,x) end if @@ -4484,14 +4559,9 @@ program main call clean_mesh(malha, mesh, domx, domy,id,.false.) ! print*, "bbb", id - ! if (id == 0) then - ! print*,'Press Return to continue' - ! read(*,*) - ! end if + call MPI_barrier(MPI_COMM_WORLD, ierr) - ! id = 0 é o processo raiz - if (id == 0) call cpu_time(start) ! Agora dividimos para cada id uma região. i = 0 @@ -4528,18 +4598,18 @@ program main call comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,domy,wall,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 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,propriedade,domx,domy, Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells,id) + call comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy, Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells) ! COMP F call comp_F_thermo(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera força ! COMP V - call comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, hot_cells, cold_cells,id) + call comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, hot_cells, cold_cells) t = t + dt cont2 = cont2+1 @@ -4724,15 +4794,27 @@ program main t = t + dt cont2 = cont2+1 - if (i == interv_Td(j2)) then - if (Td > 0 .and. termostato_vel_scaling) then + if (termostato_vel_scaling .and. i == interv_Td(j2)) then + if (Td > 0) then call corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,t) end if - if (termostato_vel_scaling .and. (Td_hot > 0 .and. Td_cold > 0)) then + if (Td_hot > 0 .and. Td_cold > 0) then call corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedade, np,id,t) end if j2 = j2 + 1 end if + + if (force_lingradT > 0 .and. i == interv_Td(j2)) then + do aux4 = 0, force_lingradT/subx + + domx_aux = [ (domx(2)-domx(1))/(force_lingradT/subx) + + call call corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,t) + + + j2 = j2 + 1 + end do + end if if (i == interv(j)) then deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, & @@ -4894,7 +4976,8 @@ program main ! if (print_TC == 1) deallocate(nRfu,nRfu_send,rFUp) if (id == 0) then - call cpu_time(finish) + call system_clock(ic1,cpu_countrate) + finish = real(ic1)/real(cpu_countrate,kind(0.d0)) print '("Time = ",f10.3," seconds.")',(finish-start) open(unit=22,file='settings.txt',status="old", position="append", action="write") write(22,*) "#:Execution time = ",(finish-start)," seconds." diff --git a/lennard_jones_stability_analysis.py b/lennard_jones_stability_analysis.py new file mode 100644 index 0000000..059e594 --- /dev/null +++ b/lennard_jones_stability_analysis.py @@ -0,0 +1,188 @@ + +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import quad + + +def unprecize(step,p1,p2,sigma,epsil): + r = np.linspace(p1,p2,int((p2-p1)/step)) + f = (1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*r + I = 0 + for i in range(len(f)-1): + I += 0.5*(f[i+1]+f[i])*step + return I + +def trapz(x,y): + I = 0 + for i in range(len(x)-1): + I += 0.5*(y[i+1]+y[i])*(x[i+1]-x[i-1]) + return I + + +def flj(r,sigma, epsil): + return (1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil*r + +print("O que fazer?\n1 - Integrar LJ trapezoidal.\n2 - Distribuição Maxwell-Boltzmann\n3 - Interagir duas partículas.\n") +tarefa = input("Entre 1, 2 ou 3: ") + +if tarefa == "1": + + v_max = 5 + n_points = 100 + dt = np.logspace(-5,-1,n_points) + step = v_max*dt + sigma = 1 + epsil = 1 + m = 1 + p1,p2 = sigma*.98, 3*sigma + desv = np.zeros(n_points) + + r = np.linspace(p1,p2,10000) + f = flj(r,sigma,epsil) + v = 4*epsil*((sigma/r)**6*((sigma/r)**6)-1) + + E_precisa = quad(flj,p1,p2,args=(sigma,epsil)) + + for i in range(n_points): + desv[i] = (unprecize(step[i],p1,p2,sigma,epsil)-E_precisa[0])/E_precisa[0] + + + fig, ax1 = plt.subplots() + ax2 = ax1.twinx() + + dts = [0.00001, 0.0001, 0.001, 0.01, 0.02] + + lns2 = ax2.plot(r,0.00001*f/m,label='v after dt = 0.00001') + lns3 = ax2.plot(r,0.0001*f/m,label='v after dt = 0.0001') + lns4 = ax2.plot(r,0.001*f/m,label='v after dt = 0.001') + lns5 = ax2.plot(r,0.01*f/m,label='v after dt = 0.01') + lns6 = ax2.plot(r,0.01*f/m,label='v after dt = 0.02') + #lns7 = ax2.plot(r,0.05*f/m,label='v after dt = 0.05') + ax2.set_ylabel('speed') + + lns1 = ax1.plot(r,f,'-k',label='force') + ax1.set_xlabel('r') + ax1.set_ylabel('force') + + leg = lns2 + lns3 + lns4 + lns5 + lns6 #+ lns7 + labs = [l.get_label() for l in leg] + ax1.legend(leg, labs, loc=0) + + plt.show() + + plt.semilogx(dt,desv,label='desvio') + plt.legend() + + plt.show() +elif tarefa == "2": + +# pdf maxwell + + plt.figure(7) + ans = "y" + while (ans == "y"): + + prop = input("Enter Temperature and mass and v_max for the graph: ").split() + T = float(prop[0]) + m = float(prop[1]) + + print("mean v = {} at T + {}".format(np.sqrt(2*T),T)) + + v = np.linspace(0,float(prop[2]),1000) + + pdf = (m/(2*np.pi*T))**(3/2)*(4*np.pi*v**2)*np.exp(-(m*v**2/2*T)) + + lab = "T = " + prop[0] + " m = " + prop[1] + plt.plot(v,pdf,label=lab) + plt.xlabel("Velocity") + plt.ylabel("Probabilty density") + ans = input("Draw one more curve? [y/n]") + plt.legend() + plt.show() + +elif tarefa == 3: + +# caso para duas partículas muito rápidas se aproximando + + + t_fim = 2 + print("tfim = {}\n".format(t_fim)) + vini = -float(input('V inicial >0 ? ')) + + p1 = 3*sigma/2 + p2 = -3*sigma/2 + ch = 'c' + w = np.zeros(len(dts)) + j = 0 + while ch != 'cq': + for dt1 in dts: + + v1, v2 = np.zeros(int(t_fim/dt1)), np.zeros(int(t_fim/dt1)) + x1, x2 = np.zeros(int(t_fim/dt1)), np.zeros(int(t_fim/dt1)) + F = np.zeros(int(t_fim/dt1)) + x1[0] = p1 + x2[0] = p2 + v1[0] = vini + v2[0] = -vini + for i in range(len(v1)-1): + f1 = -flj(x1[i]-x2[i],sigma,epsil) + + x1[i+1] = x1[i] + dt1*v1[i] + f1*dt1**2/(2*m) + x2[i+1] = x2[i] + dt1*v2[i] + (-f1)*dt1**2/(2*m) + + f2 = -flj(x1[i+1]-x2[i+1],sigma,epsil) + F[i] = f2 + + v1[i+1] = v1[i] + (f1+f2)*dt1/(2*m) + v2[i+1] = v2[i] -(f1+f2)*dt1/(2*m) + + temp = np.linspace(0,t_fim, int(t_fim/dt1)) + plt.figure(3) + plt.plot(temp,v1,label='V1 dt = '+str(dt1)) + #plt.plot(temp,v2,label='V2 dt = '+str(dt1)) + plt.figure(4) + plt.plot(temp,x1,label='X1 dt = '+str(dt1)) + plt.figure(5) + plt.plot(temp,F,label='F1 dt '+str(dt1)) + w[j] = trapz(x1,F) + j += 1 + + plt.figure(4) + plt.xlabel('Time') + plt.ylabel('Distance') + plt.legend() + + plt.figure(3) + plt.xlabel('Time') + plt.ylabel('Speed') + plt.legend() + + plt.figure(5) + plt.xlabel('temp') + plt.ylabel('force') + plt.legend() + + plt.figure(6) + plt.semilogx(dts,w) + plt.xlabel('dt') + plt.xlabel('work on one particle') + + plt.show() + + ch = input("Change something? \nEnter vini for initial velocity;\nEnter p1 for initial position,\nEnter tfim for final time.\n> Enter c to continue or cq to quit.\n") + while (ch[0] != 'c'): + if ch == 'vini': + vini = -float(input('-V inicial? ')) + elif ch == 'p1': + vini = float(input('p1 inicial? ')) + elif ch == 'tfim': + t_fim = float(input('tfim ? ')) + ch = input("Change something? \nEnter vini for initial velocity;\nEnter p1 for initial position,\nEnter tfim for final time.\n> Enter c to continue or cq to quit.\n") + + + + #plt.plot(temp,x1,label='X1 dt = '+str(dt1)) + #plt.plot(temp,x2,label='X2 dt = '+str(dt1)) + + + diff --git a/part_gen_pequeno.py b/part_gen_pequeno.py index 4c2b1e3..b0198e7 100644 --- a/part_gen_pequeno.py +++ b/part_gen_pequeno.py @@ -26,8 +26,8 @@ ff = 1 # formato da região ff = x/y Lx = 400 -arquivo2 = 'poucas1.csv' -arquivo1 = 'poucas2.csv' +arquivo1 = 'poucas1.csv' +arquivo2 = 'poucas2.csv' ratio = 0.8 N1 = 500 diff --git a/potential_analysis.py b/potential_analysis.py index f601416..a21d543 100644 --- a/potential_analysis.py +++ b/potential_analysis.py @@ -5,11 +5,65 @@ import matplotlib.pyplot as plt from glob import glob import time +from zipfile import ZipFile -a = input("Enter the subdomains mesh dimensions.\n") -a = a.split() -mesh = np.array([int(a[0]), int(a[1])]) +try: + import progressbar + pbar = True +except: + print('progressbar not available. Try one:') + print('conda install -c conda-forge progressbar2') + print('conda install -c conda-forge/label/gcc7 progressbar2') + pbar = False + +kB = 1.38064852*10**-23 + +def integral(F,i,j,dt): + I = np.trapz(F[i,j,:],dx=dt) + return I + +def diff1(F,h): + B = np.zeros((len(F),1)) + for i in range(len(F)): + if i > 1 and i < len(F) - 2: + B[i] = ((1/12)*F[i-2] + (-2/3)*F[i-1] + (2/3)*F[i+1] + (-1/12)*F[i+2])/h + elif i == 0: + B[i] = ((-11/6)*F[0] + 3*F[1] + (-3/2)*F[2] + (1/3)*F[3])/h + elif i == 1: + B[i] = (-2*F[0]-3*F[1]+6*F[2]-1*F[3])/(6*h) + elif i == len(F)-1: + B[i] = (1*F[i-2]-4*F[i-1]+3*F[i+0])/(2*1.0*h**1) + elif i == len(F)-2: + B[i] = (1*F[i-2]-6*F[i-1]+3*F[i+0]+2*F[i+1])/(6*1.0*h**1) + return B +def greenkubo(Tk,Tp,i,j,na,dt): + itv = int(Tp.shape[2]/na) + kk = np.zeros((itv)) + kp = np.zeros((itv)) + pp = np.zeros((itv)) + for a in range(itv): + for b in range(na): + kk[a] += np.trapz(Tk[i,j,b*itv:b*itv+a],dx=dt)**2 + kp[a] += np.trapz(Tp[i,j,b*itv:b*itv+a],dx=dt)*np.trapz(Tk[i,j,b*itv:b*itv+a],dx=dt) + pp[a] += np.trapz(Tp[i,j,b*itv:b*itv+a],dx=dt)**2 + kk = kk/na + pp = pp/na + kp = kp/na + return kk, kp, pp + +def einstein_condutividade(delta_e,na): + lambda_T = np.zeros(delta_e.shape[0], delta_e.shape[1], na) + itv = int(delta_e.shape[2]/na) + for a in range(itv): + for b in range(na): + for i in range(len(delta_e.shape[0])): + for j in range(len(delta_e.shape[1])): + lambda_T[i,j,a] += (delta_e[i,j,b*itv] - delta_e[i,j,b*itv+a])**2 + lambda_T = lambda_T/na + return lambda_T + + dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__)) dirlist = glob(dirname + "/*/") print("Choose a folder there the results are contained:\nNo | Folder") @@ -18,6 +72,12 @@ a = int(input("Enter the number of the folder\n")) res_dir = dirlist[a] +zip_rfup = ZipFile(res_dir+'/rFuP.zip','r') +zip_positions = ZipFile(res_dir+'/positions.zip','r') +zip_velocities = ZipFile(res_dir+'/velocities.zip','r') + +len_list_files = len(zip_positions.namelist()+zip_velocities.namelist()) + config = configparser.ConfigParser() config.read(res_dir + 'settings.txt') N = int(config['global']['N'].split()[0]) @@ -25,19 +85,34 @@ dimy = float(config['global']['dimY'].split()[0]) n_files = int(config['out_files']['out_files'].split()[0]) ntype = int(config['global']['Ntype'].split()[0]) +t_fim = float(config['global']['t_fim'].split()[0]) +#dt = float(config['global']['dt'].split()[0]) +F = np.array([0,0]) quant = [] sigma = [] epsilon = [] rs = [] # raio sólido -m = [] +mass = [] tipo = [0]*N +dt = t_fim/n_files + +a = input("Enter the subdomains mesh dimensions.\n") +a = a.split() +mesh = np.array([int(a[0]), int(a[1])]) +a = input("Enter a location (Starts at 0). xmin xmax ymin ymax\n") +if a == '': + region = [0, mesh[0], 0, mesh[1]] +else: + region = [int(x) for x in a.split()] + +Vol = (dimx/mesh[0])*(dimy/mesh[1]) # volume dos elementos da malha for i in range(ntype): quant.append(int(config['par_'+str(i)]['quantidade'].split()[0])) rs.append(float(config['par_'+str(i)]['rs'].split()[0])) sigma.append(float(config['par_'+str(i)]['sigma'].split()[0])) epsilon.append(float(config['par_'+str(i)]['epsilon'].split()[0])) - m.append(float(config['par_'+str(i)]['m'].split()[0])) + mass.append(float(config['par_'+str(i)]['m'].split()[0])) j,k = 0,0 for i in range(len(quant)): for j in range(quant[i]): @@ -47,148 +122,75 @@ tipo = pd.DataFrame(tipo, columns=["tipo"]) # numero id da partícula hx = dimx/mesh[0] -hy = dimx/mesh[1] -particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])] -density_map = np.zeros((mesh[0],mesh[1])) -kinetic_map = np.zeros((mesh[0],mesh[1])) -potential_map = np.zeros((mesh[0],mesh[1])) -step = input('Choose step:\n') -n1,n2 = 0,0 -while int(step) >= 0: - pos = pd.read_csv(res_dir+"position.csv."+step, header=None, names = ["x","y"]) - vel = pd.read_csv(res_dir+"velocity.csv."+step, header=None, names = ["v_x", "v_y"]) +hy = dimy/mesh[1] + +nsteps = n_files # int(input('Enter the number of steps:\n')) + +density_map = np.zeros((mesh[0],mesh[1], nsteps+1)) +Kmap = np.zeros((mesh[0],mesh[1],nsteps+1)) +Vmap = np.zeros((mesh[0],mesh[1],nsteps+1)) +step = 1 + +if pbar: + bar = progressbar.ProgressBar(max_value=nsteps) +while step <= nsteps: + particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])] + rfup = pd.read_csv(zip_rfup.open('rF_u_P.csv.'+str(step)), header=None, names = ["RxFy","RyFx","u","px","py"]) + pos = pd.read_csv(zip_positions.open("position.csv."+str(step)), header=None, names = ["x","y"]) + vel = pd.read_csv(zip_velocities.open("velocity.csv."+str(step)), header=None, names = ["v_x", "v_y"]) n = [x for x in range(len(pos))] n = pd.DataFrame(n, columns=["n"]) # numero id da partícula + # print("AA") pos_vel = pd.concat([n,pos,vel,tipo],axis=1) - - for nn in range(len(pos_vel)): + # print("B") + for nn in range(quant[0]): # só vai contar as particulas do grupo 1 xp = int(pos_vel.loc[nn,'x']//hx) yp = int(pos_vel.loc[nn,'y']//hy) if xp == mesh[0]: xp = xp - 1 if yp == mesh[1]: yp = yp - 1 - # print('xp = {}, yp = {}\n'.format(xp,yp)) particle_map[xp][yp].append( pos_vel.loc[nn,'n'] ) - # print(particle_map) - # time.sleep(0.5) - density_map[xp,yp] += 1 - kinetic_map[xp,yp] += 0.5*np.sqrt(pos_vel.loc[nn,'v_x']**2 + pos_vel.loc[nn,'v_y']**2)*m[pos_vel.loc[nn,'tipo']] - - for i in range(mesh[0]): - for j in range(mesh[1]): + density_map[xp,yp,step] += 1 + # "print C" + for i in range(region[0],region[1]): + for j in range(region[2],region[3]): for nn in range(len(particle_map[i][j])): n1 = particle_map[i][j][nn] - ep1 = epsilon[pos_vel.loc[n1,'tipo']] - sigma1 = sigma[pos_vel.loc[n1,'tipo']] - rs1 = rs[pos_vel.loc[n1,'tipo']] - # na própria célula - for nn2 in range(nn+1,len(particle_map[i][j])): - n2 = particle_map[i][j][nn2] - rs2 = rs[pos_vel.loc[n2,'tipo']] - ep2 = epsilon[pos_vel.loc[n2,'tipo']] - sigma2 = sigma[pos_vel.loc[n2,'tipo']] - ep = np.sqrt(ep1*ep2) - s = 0.5*(sigma2 + sigma1) - r = np.sqrt((pos_vel.loc[n1,'x'] - pos_vel.loc[n2,'x'])**2 + (pos_vel.loc[n2,'y'] - pos_vel.loc[n1,'y'])**2) - rs1 - rs2 - - potential_map[i,j] += 4*ep*((s/r)**12 - (s/r)**6) - # na celula vizinha i+1, j - if i+1 < mesh[0]: - for nn2 in range(1,len(particle_map[i+1][j])): - n2 = particle_map[i+1][j][nn2] - rs2 = rs[pos_vel.loc[n2,'tipo']] - ep2 = epsilon[pos_vel.loc[n2,'tipo']] - sigma2 = sigma[pos_vel.loc[n2,'tipo']] - ep = np.sqrt(ep1*ep2) - s = 0.5*(sigma2 + sigma1) - r = np.sqrt((pos_vel.loc[n1,'x'] - pos_vel.loc[n2,'x'])**2 + (pos_vel.loc[n2,'y'] - pos_vel.loc[n1,'y'])**2) - rs1 - rs2 - if r == 0: - print('nn {} {}'.format(nn,nn2)) - time.sleep(2) - U = 4*ep*((s/r)**12 - (s/r)**6) - potential_map[i,j] += U/2 - potential_map[i+1,j] += U/2 - # na celula vizinha i+1, j+1 - if i+1 < mesh[0] and j+1 < mesh[1]: - for nn2 in range(1,len(particle_map[i+1][j+1])): - n2 = particle_map[i+1][j+1][nn2] - rs2 = rs[pos_vel.loc[n2,'tipo']] - ep2 = epsilon[pos_vel.loc[n2,'tipo']] - sigma2 = sigma[pos_vel.loc[n2,'tipo']] - ep = np.sqrt(ep1*ep2) - s = 0.5*(sigma2 + sigma1) - r = np.sqrt((pos_vel.loc[n1,'x'] - pos_vel.loc[n2,'x'])**2 + (pos_vel.loc[n2,'y'] - pos_vel.loc[n1,'y'])**2) - rs1 - rs2 - if r == 0: - print('nn {} {}'.format(nn,nn2)) - time.sleep(2) - U = 4*ep*((s/r)**12 - (s/r)**6) - potential_map[i,j] += U/2 - potential_map[i+1,j+1] += U/2 - # na celula vizinha i, j+1 - if j+1 < mesh[1]: - for nn2 in range(1,len(particle_map[i][j+1])): - n2 = particle_map[i][j+1][nn2] - rs2 = rs[pos_vel.loc[n2,'tipo']] - ep2 = epsilon[pos_vel.loc[n2,'tipo']] - sigma2 = sigma[pos_vel.loc[n2,'tipo']] - ep = np.sqrt(ep1*ep2) - s = 0.5*(sigma2 + sigma1) - r = np.sqrt((pos_vel.loc[n1,'x'] - pos_vel.loc[n2,'x'])**2 + (pos_vel.loc[n2,'y'] - pos_vel.loc[n1,'y'])**2) - rs1 - rs2 - if r == 0: - print('nn {} {}'.format(nn,nn2)) - time.sleep(2) - U = 4*ep*((s/r)**12 - (s/r)**6) - potential_map[i,j] += U/2 - potential_map[i,j+1] += U/2 - # na celula vizinha i-1, j+1 - if i-1 >= 0 and j+1 < mesh[1]: - for nn2 in range(1,len(particle_map[i-1][j+1])): - n2 = particle_map[i-1][j+1][nn2] - ep2 = epsilon[pos_vel.loc[n2,'tipo']] - rs2 = rs[pos_vel.loc[n2,'tipo']] - sigma2 = sigma[pos_vel.loc[n2,'tipo']] - ep = np.sqrt(ep1*ep2) - s = 0.5*(sigma2 + sigma1) - r = np.sqrt((pos_vel.loc[n1,'x'] - pos_vel.loc[n2,'x'])**2 + (pos_vel.loc[n2,'y'] - pos_vel.loc[n1,'y'])**2) - rs1 - rs2 - if r == 0: - print('nn {} {}'.format(nn,nn2)) - time.sleep(2) - U = 4*ep*((s/r)**12 - (s/r)**6) - potential_map[i,j] += U/2 - potential_map[i-1,j+1] += U/2 - - mean_Kenergy = (1/2)*kinetic_map**2/density_map - mean_Kenergy[mean_Kenergy == np.inf] = 0 - mean_Kenergy = np.nan_to_num(mean_Kenergy) - - X = np.linspace(0,dimx,mesh[0]) - Y = np.linspace(0,dimy,mesh[1]) - X,Y = np.meshgrid(Y,X) + m = mass[pos_vel.loc[n1,'tipo']] + Kmap[i,j,step] += (pos_vel.loc[n1,'v_x']**2+pos_vel.loc[n1,'v_y']**2)*m/(2) + # Vmap[i,j,step] += pos_vel.loc[n1,'u'] + if pbar: + bar.update(step) + step += 1 + +KE = Kmap/density_map +a = KE.shape[2] +KE = np.nan_to_num(KE) +KE = np.sum(KE,axis=2)/a +print("nsteps {} = a {} ?".format(nsteps,a)) +plt.figure() +if mesh[1] == 1: + T = KE #np.sum(KE,axis=0) # energia cinética ao longo de X + plt.plot(T,label='Temperatura adimensional ao longo de x') +elif mesh[0] == 1: + T = KE #np.sum(KE,axis=1) # energia cinética ao longo de Y + plt.plot(T,label='Temperatura adimensional ao longo de y') +plt.legend() - # gráfico - fig1 = plt.figure(1) - kmap = plt.contourf(X,Y,mean_Kenergy,cmap="coolwarm") - cbar1 = fig1.colorbar(kmap) - cbar1.ax.set_ylabel('Kinetic Energy') +# plt.figure() - fig2 = plt.figure(2) - dmap = plt.contourf(X,Y,density_map) - cbar2 = fig2.colorbar(dmap) - cbar2.ax.set_ylabel('Number of density') +# Vlj = Vmap/a +# if mesh[1] == 1: +# V = np.sum(Vlj,axis=0) # energia cinética ao longo de X +# plt.plot(T,label='Potencial adimensional ao longo de x') +# elif mesh[0] == 1: +# V = np.sum(Vlj,axis=1) # energia cinética ao longo de Y +# plt.plot(T,label='Potencial adimensional ao longo de y') +# plt.legend() - # fig, ax1 = plt.subplots() - # ax1.plot(np.linspace(0,dimx,subspaces),mean_Kenergy,'b') - # ax1.set_xlabel("x - direction") - # ax1.set_ylabel("Mean kinetic energy", color='b') +plt.show() - # ax2 = ax1.twinx() - # ax2.plot(np.linspace(0,dimx,subspaces),concentration_map, 'r') - # ax2.set_ylabel("Number of density",color='r') - # #plt.axis([0,dimx,np.min(mean_Kenergy),np.max(mean_Kenergy)]) - # fig.tight_layout() - plt.show(block=True) - step = input('Choose step or enter -1 to stop:\n') \ No newline at end of file diff --git a/saida.f90 b/saida.f90 index 492548a..bfc911e 100644 --- a/saida.f90 +++ b/saida.f90 @@ -7,7 +7,8 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) implicit none integer, intent(in) :: n,d real(dp), intent(in) :: v(n,d), t, start - integer :: i,step, nimpre + integer :: i,step, nimpre, ic1, cpu_countrate, horas, min + real(dp) :: sec character(*) :: prop character(4) :: extensao = '.csv', passo real(dp) :: time @@ -57,15 +58,26 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) if (prop == "position") then if (step == 0) timep = 0 - call cpu_time(time) + ! call cpu_time(time) + call system_clock(ic1,cpu_countrate) + time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0)) if (step == 1) then - etc = ((time - start)/step+ (time-timep))*0.5*nimpre - (time - start) + call system_clock(ic1,cpu_countrate) + time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0)) + etc = ((time - start)/real(step,kind(0.d0)) + (time-timep))*0.5*real(nimpre,kind(0.d0)) - (time - start) dtimepp = (time-timep) print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc else if (step > 1) then - etc = ((time - start)*2/step + ((time-timep)*4 + dtimepp*4))*(nimpre-step)/10 - print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc - dtimepp = (time-timep) + call system_clock(ic1,cpu_countrate) + time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0)) + etc = ((time - start)*6/real(step,kind(0.d0)) + ((time-timep)*2 + dtimepp*2))* & + (real(nimpre,kind(0.d0))-real(step,kind(0.d0)))/10 + horas = int(etc/3600) + min = (int(etc) - horas*3600)/60 + sec = etc - real(horas*3600 + min*60,kind(0.d0)) + print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", i3,":",i2,":",f4.1 )' & + ,prop//extensao//'.'//trim(passo),t,horas, min, sec + dtimepp = (time-timep) else print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", "unknown" )',prop//extensao//'.'//trim(passo),t end if @@ -77,7 +89,8 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) else print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", "unknown" )',prop//extensao//'.'//trim(passo),t end if - timep = time + ! time = real(ic1,kind(0.d0))/real(cpu_countrate,kind(0.d0)) + ! timep = time end if ! print*, 'Salvo arquivo ',prop//extensao//'.'//trim(passo), "t =", t, "ETC: ", etc diff --git a/settings.ini b/settings.ini index 526fa4f..e13ee66 100755 --- a/settings.ini +++ b/settings.ini @@ -3,27 +3,28 @@ # Os outros parâmetros são adimensionais. Checar anotações.txt T' = T*k_b/epsilon # Argonio: epsilon/kb = 120K, sigma = 0.341nm [global] - nimpre = 50 # quntidade de saidas - N = 400 # número de partículas - Ntype = 1 # número de tipos de partícula - t_fim = 20 # tempo de execução + nimpre = 100 # quntidade de saidas + N = 584 # número de partículas + Ntype = 2 # número de tipos de partícula + t_fim = 50 # tempo de execução nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). - dt = 0.0005 # passo de tempo + dt = 0.0001 # passo de tempo dimX = 400 # Dimensões da região de cálculo dimY = 400 # - mesh = 50 50 # malha(elementos) + mesh = 80 80 # 300 60 # malha(elementos) rcut = 3 # *sigma - wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste - termostato_nose_hoover = .true. #liga termostato Nose Hoover. Este funciona durante toda execução. + wall = 'ppee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste + termostato_nose_hoover = .false. #liga termostato Nose Hoover. Este funciona durante toda execução. termostato_vel_scaling = .false. # liga termostato por velocity scaling Mc = 1 Td = -1 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E/epsilon(literatura) - temp_Td = 0 10 1000 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] - cold_cells = 1 15 1 50 # [cellx_ini, cellx_fin, celly_ini, celly_fin] - hot_cells = 16 50 1 50 # [cellx_ini, cellx_fin, celly_ini, celly_fin] - Td_hot = 1 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato - Td_cold = 1 #0.2 # temperatura nas hot cells constante. Td_cold = -1 desliga termostato - vd = 1 1 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') + temp_Td = 0 10 100 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] + cold_cells = 1 10 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin] banho frio + hot_cells = 190 200 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin] banho quente + Td_hot = 1 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato. + Td_cold = 0.2 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato + force_lingradT = 10 # se > 0 irá dividir região no numero de subregiões entradas na direação X e será forçado um gradiente de temperatura linear. A região dos banhos será ignorada. A temperatura será controlado por um termostato velocity scaling + vd = 0 0 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') NMPT = 1000 # Número máximo de partículas que poderá mudar de process. Se não souber estimar, usar -1 ou o valor de N. GField = 0 0 # Campo de força gravitacional que afeta todas as partículas print_TC = .false. # imprimir dados para calcular coeficiente de transporte @@ -32,31 +33,31 @@ # se v_file é dada, então ela será usada para dar a velocidade às partículas. Deve ser arquivo .csv. Deixar "%%" se não usar. # se ismolecule = .true. então a partícula é levada em conta no termostato [par_0] - x = 'poucas2.csv' + x = 'poucas1.csv' v = 0 0 # velocidade global - v_file = '%%v_file_0.csv' # velocidade de cada partícula + v_file = 'v_file_0.csv' # velocidade de cada partícula m = 1 nome = 'g1' sigma = 1 # Lennard Jones epsilon = 1 # Lennard Jones - quantidade = 400 + quantidade = 484 x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay rs = 0 # raio sólido. Posição de partículas na superfície de um sólido de raio rs fric_term = 0 # fricção artifical não funciona quando o termostato nose hoover está ligado ismolecule = .true. -#[par_1] -# x = 'poucas1.csv' -# v = 0 0 # velocidade global -# v_file = '%%v_file_0.csv' # velocidade de cada partícula -# m = 1 -# nome = 'g2' -# sigma = 1 # Lennard Jones -# epsilon = 1 # Lennard Jones -# quantidade = 20 -# x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay -# rs = 0 # raio sólido. Posição de partículas na superfície de um sólido de raio rs -# fric_term = 0 # fricção artifical -# ismolecule = .false. +[par_1] + x = 'poucas2.csv' + v = 0 0 # velocidade global + v_file = '%%v_file_0.csv' # velocidade de cada partícula + m = 50 + nome = 'g2' + sigma = 1 # Lennard Jones + epsilon = 1 # Lennard Jones + quantidade = 100 + x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay + rs = 1 # raio sólido. Posição de partículas na superfície de um sólido de raio rs + fric_term = 0 # fricção artifical + ismolecule = .false. #[par_2] # x = 'p_g.csv' # v = 0 0