From ffd21376e8ef3b2628f008c3cb3fff3a99816946 Mon Sep 17 00:00:00 2001 From: g7fernandes Date: Wed, 4 Sep 2019 17:42:05 -0300 Subject: [PATCH] =?UTF-8?q?implementado=20rota=C3=A7=C3=A3o=20-=20Em=20tes?= =?UTF-8?q?tes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- debug_lennard.sh | 6 +- lennard.f90 | 2257 ++++++++++++++++++++++++++++++++++++---------- mod0.f90 | 1 + pack.sh | 2 +- rot_par.py | 80 ++ saida.f90 | 2 +- settings.ini | 36 +- 7 files changed, 1902 insertions(+), 482 deletions(-) create mode 100644 rot_par.py diff --git a/debug_lennard.sh b/debug_lennard.sh index 879295a..6cb8520 100755 --- a/debug_lennard.sh +++ b/debug_lennard.sh @@ -18,9 +18,9 @@ else mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -O0 -shared-intel -nocheck -fno-omit-frame-pointer -fasynchronous-unwind-tables -fexceptions -o lennard.out echo "output: lennard.out" else - echo "compiling with intel mpiifort without options" - mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -o lennard.out - echo "output: lennard" + echo "compiling with intel mpiifort without optimization" + mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -O0 -g -traceback -check all -o lennard.out + echo "output: lennard.out" fi fi diff --git a/lennard.f90 b/lennard.f90 index dee7dc6..a3fe754 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -635,7 +635,7 @@ subroutine corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedad use mpi type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade - real(dp),intent(in) :: Td_hot, Td_cold!energia cinética + real(dp),intent(in) :: Td_hot, Td_cold !energia cinética real(dp) :: T_hot, T_cold,K,kb = 1.38064852E-23, betac, betah integer, intent(in) :: domx(2),domy(2), cold_cells(4), hot_cells(4) type(container), allocatable,dimension(:,:),intent(in) :: malha @@ -689,7 +689,7 @@ function comp_fric(r,v,fric_term) result(f) end function comp_fric ! atualiza forças - subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + subroutine comp_FT(GField,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,domy, ids, id, t,dt,partlst) use linkedlist use mod1 use data @@ -698,15 +698,42 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha - real(dp), intent(in) :: GField(2), t + real(dp), intent(in) :: dt, GField(2), pr(5), t ! outras propriedades da particula que gira [A, B, alpha, beta, phase] + real(dp), intent(inout), dimension(:) :: Tor ! torque + real(dp), intent(inout), dimension(:) :: theta ! ângulo das partículas que giram, tangencial real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção - real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine - integer :: i,j,ct = 0, dox(2), doy(2) !,ptr, ptrn + real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine, A, B, alpha,beta,ph + integer :: i,j, ct = 0, dox(2), doy(2) !,ptr, ptrn integer, intent(in) :: mesh(:),domx(2),domy(2), id - real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2 + real(dp) :: Fi(2)=0,r, aux2(2),aux3,fR(2), fric_term1, fric_term2 + integer, save :: pa type(list_t), pointer :: node, next_node type(data_ptr) :: ptr,ptrn + integer, intent(out), dimension(:) :: partlst integer ( kind = 4 ), intent(in) :: ids(8) + + A = pr(1) + B = pr(2) + alpha = pr(3) + beta = pr(4) + ph = pr(5) + + ! quanta a quantidade de partículas antes das que giram + print*, "t,2*dt", t, 2*dt + if (t < 2*dt) then + rs = 0 + i = 0 + pa = 0 + do while (rs == 0) + i = i+1 + rs = propriedade(i)%rs + pa = pa + propriedade(i)%quant + end do + pa = pa - propriedade(i)%quant + print*, "pa calculado!",id,pa + end if + print*, "pa = ", pa + !Lennard Jones fR = [0,0] @@ -731,9 +758,12 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) ! ! print*, "L 253", ptr%p%n ! end if do while (associated(node)) + ! print*, "Encontrada", ct, "celulas, id=", id ptr = transfer(list_get(node), ptr) !particula selecionada + ! PRINT*, "L 759" ptr%p%flag = .true. ! indica ao comp_x que a partícula precisa ser calculada + ! PRINT*, "L 764" x1 = ptr%p%x v1 = ptr%p%v m1 = propriedade(ptr%p%grupo)%m @@ -745,12 +775,14 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then ptr%p%F = ptr%p%F + GField*m1 end if + ! print '("x1 =", f10.6, " ", f10.6, " n ", i2, " cell ",i2," ",i2)',x1(1),x1(2), ptr%p%n,i,j ! if (id == 0) read(*,*) !calcular a força desta com todas as outras partículas next_node => list_next(node) ! próxima partícula da célula node => list_next(node) ! a ser computado com a particula selecionada + ! PRINT*, "L 780" ! NA PRÓPRIA CELULA do while (associated(node)) ptrn = transfer(list_get(node), ptrn) !outra particula selecionada @@ -780,22 +812,78 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) r = r - rs1 - rs2 !raio ! print*, "L 389 r", r, "id",id if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ! print*, "L 395 r", r, "id",id - !print '("(x1(1)-x2(1)) ", f7.3 ," (rs1+rs2)*coss ", f7.3, " (x1(2)-x2(2)) ", f7.3, " (rs1+rs2)*sine ", f7.3 )', (x1(1)-x2(1)), & - !(rs1+rs2)*coss, (x1(2)-x2(2)), (rs1+rs2)*sine - fric_term = (fric_term1+fric_term2)/2 - if (fric_term > 0) then - fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + if (rs1 == 0 .and. rs2 == 0) then + aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & + ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + ! print*, "r1 > 0 gamma", gamma + ! print*, "aux2, aux3",aux2, aux3 + ! print*, "x1 =", x1 + ! print*, "x2 =", x2 + ! read(*,*) + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + ! print*, "r1 > 0 gamma", gamma + ! print*, "aux2, aux3",aux2, aux3 + ! print*, "x1 =", x1 + ! print*, "x2 =", x2 + ! read(*,*) + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR - end if + end if node => list_next(node) ! próxima partícula da célula end do - + ! print*, "L 871" !Células ao redor !i é linha, j é coluna if (i /= mesh(2)+2) then ! se não for a última linha if (j == mesh(1)+2) then ! se for a última coluna @@ -828,17 +916,66 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 666 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - fric_term = (fric_term1+fric_term2)/2 + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1))*(sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)* & + (B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 & + * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula end do @@ -872,14 +1009,64 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 710 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 & + * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & + (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)* & + (B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & + ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & + (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 & + - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR end if @@ -919,14 +1106,60 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 757 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & + ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) & + * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)*( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR end if @@ -961,14 +1194,61 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 798 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & + ((sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR end if @@ -1004,14 +1284,60 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 841 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & + ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR end if @@ -1048,14 +1374,60 @@ subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) sine = (x1(2)-x2(2))/r r = r - rs1 - rs2 !raio ! print*, "L 885 r", r, "id",id - if (r <= rcut) then + if (rs1 == 0 .and. rs2 == 0) then aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + else if (rs1 > 0 .and. rs2 == 0) then + + gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 & + * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptr%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptr%p%n - pa) = 1 + else if (rs2 > 0 .and. rs1 == 0) then + gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * & + (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & + ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & + 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & + ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) + + Tor(-pa +ptrn%p%n) = aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) + ptrn%p%F = aux2 + ptr%p%F + ptr%p%F = -aux2 + ptrn%p%F + aux3*[-sine,-coss] + partlst(ptrn%p%n - pa) = 1 + else + ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos + aux2 = -(1/r**2)*(sigma*(1+B)/r)**6*(1-2*(sigma*(1+B)/r)**6)*24*epsil*(1+A)* & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + ! print*, "L 395 r", r, "id",id + + ! talvez o fricterm seja interessante em algum momento, vou deixar aqui + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) end if + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR end if @@ -1923,383 +2295,996 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id ! call MPI_barrier(MPI_COMM_WORLD, ierr) end subroutine comp_F_thermo - subroutine comp_FT_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + + + ! atualiza posições + subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) use linkedlist use mod1 use data use mod0 use mpi - + use matprint + + character(1) :: north, south, east, west type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha - real(dp), intent(in) :: GField(2), t - real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut - real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine - integer :: i,j,ct = 0, dox(2), doy(2) !,ptr, ptrn - integer, intent(in) :: mesh(:),domx(2),domy(2), id - real(dp) :: Fi(2)=0,r, aux2(2) - type(list_t), pointer :: node, next_node - type(data_ptr) :: ptr,ptrn - integer ( kind = 4 ), intent(in) :: ids(8) - - !Lennard Jones - dox = domx - doy = domy - if (sum(ids(1:4)) > -4) then - !caso paralelo - if (domx(1) > 1) dox(1) = domx(1) - 1 - if (domy(1) > 1) doy(1) = domy(1) - 1 - if (domx(2) < mesh(1)+2) dox(2) = domx(2) + 1 - if (domy(2) < mesh(2)+2) doy(2) = domy(2) + 1 - else - dox = domx - doy = domy + ! integer, intent(inout) :: mic(:) + integer :: i,j,k, cell(2), status(MPI_STATUS_SIZE), count, destD,dx1,dx2,dy1,dy2 + integer, intent(in) :: N,mesh(2),domx(2),domy(2) + type(list_t), pointer :: node, previous_node + real(dp), intent(in) :: dt, t, dx_max + character(4), intent(in) :: wall + type(data_ptr) :: ptr + real(dp) :: x(2),m, dx(2) + real(dp), intent(in) :: icell(:), jcell(:) + integer ( kind = 4 ), intent(in) :: ids(8), id, np + integer ( kind = 4 ) :: tag,cont_db(8),cont_int(8) + type(lstr) :: LT + ! real(dp),intent(inout) :: celula(:,:) + ! IDS correspondem às celulas [N,S,E,W] + cont_db = 0 + cont_int = 0 + + north = wall(1:1) + south = wall(2:2) + east = wall(3:3) + west = wall(4:4) + ! if (id == 0) read(*,*) + + ! Esvazia as celulas emprestadas + if (domy(1) > 1) then + i = domy(1)-1 + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "A", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + else if (domy(1) == 1) then !esvazia celulas fantasma + i = domy(1) + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "A", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do end if - - do i = doy(1),doy(2) ! i é linha - do j = dox(1),dox(2) + if (domy(2) < mesh(2)+2) then + i = domy(2)+1 + do j = domx(1),domx(2) + previous_node => malha(i,j)%list node => list_next(malha(i,j)%list) - ! if (associated(node)) then - ! ptr = transfer(list_get(node), ptr) - ! ! print*, "L 253", ptr%p%n - ! end if - do while (associated(node)) - ! print*, "Encontrada", ct, "celulas, id=", id - ptr = transfer(list_get(node), ptr) !particula selecionada - ptr%p%flag = .true. ! indica ao comp_x que a partícula precisa ser calculada - x1 = ptr%p%x - v1 = ptr%p%v - m1 = propriedade(ptr%p%grupo)%m - rs1 = propriedade(ptr%p%grupo)%rs !raio sólido - sigma_a = propriedade(ptr%p%grupo)%sigma - epsil_a = propriedade(ptr%p%grupo)%epsilon - if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then - ptr%p%F = ptr%p%F + GField*m1 - end if - !calcular a força desta com todas as outras partículas - next_node => list_next(node) ! próxima partícula da célula - node => list_next(node) ! a ser computado com a particula selecionada - - ! NA PRÓPRIA CELULA - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + rs1 + rs2 - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a *epsil_b ) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + rs1 + rs2 - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a *epsil_b ) - else - rcut = r_cut*sigma_a + rs1 + rs2 - sigma = sigma_a - epsil = epsil_a + ! print*, associated(node), "B", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + else if (domy(2) == mesh(2)+2) then + i = domy(2) + do j = domx(1),domx(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "B", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + end if + if (domx(2) < mesh(1)+2) then + j = domx(2)+1 + do i = domy(1),domy(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "D", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + else if (domx(2) == mesh(1)+2) then + j = domx(2) + do i = domy(1),domy(2) + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! print*, associated(node), "D", id + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + end if + if (domx(1) > 1) then + j = domx(1)-1 + do i = domy(1),domy(2) + node => list_next(malha(i,j)%list) + previous_node => malha(i,j)%list + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + else if (domx(1) == 1) then + j = domx(1) + do i = domy(1),domy(2) + node => list_next(malha(i,j)%list) + previous_node => malha(i,j)%list + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end do + end if + ! PARA 4 processos. Limpa o que está na diagonal adjacente + if (np > 1 .and. sum(ids(5:8)) > -4) then + if (domy(1) /= 1) i = domy(1) -1 + if (domy(2) /= mesh(2)+2) i = domy(2) +1 + if (domx(1) /= 1) j = domx(1) -1 + if (domx(2) /= mesh(1) + 2) j = domx(2) +1 + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + deallocate(ptr%p) + call list_remove(previous_node) + node => list_next(previous_node) + ! node => list_next(node) + end do + end if + ! print*, "cont_int", cont_int + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + dx1 = domx(1) + dx2 = domx(2) + dy1 = domy(1) + dy2 = domy(2) + if (dy1 == 1) dy1 = 2 + if (dx1 == 1) dx1 = 2 + if (dy2 == mesh(2)+2) dy2 = mesh(2)+1 + if (dx2 == mesh(1)+2) dx2 = mesh(1)+1 + + do i = dy1, dy2 + do j = dx1,dx2 + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! ptr%p%flag = true significa que ainda não foi computado para evitar loop infinito + + if (associated(node)) then + ptr = transfer(list_get(node), ptr) + end if + + do while (associated(node) .and. ptr%p%flag) + ! print*, "L PART", i,j, "id", id + m = propriedade(ptr%p%grupo)%m + ! aqui vemos se já passou o tempo que as posições ficaram travadas. + if (propriedade(ptr%p%grupo)%x_lockdelay > t) ptr%p%flag = .false. + + if (ptr%p%flag) then + dx(1) = dt*ptr%p%v(1) + ptr%p%F(1)*dt**2/(2*m) + dx(2) = dt*ptr%p%v(2) + ptr%p%F(2)*dt**2/(2*m) + else + dx = [0,0] + end if + ptr%p%flag = .false. + if ((dx(1)**2 + dx(2)**2) >= dx_max) then + print '("Particulas rápidas demais! dx =", f18.5, " ", f18.5, " | n ", i4)', dx(1), dx(2), ptr%p%n + print*, "F",ptr%p%F, "v", ptr%p%v + print*, "x", ptr%p%x(1), ptr%p%x(2) + call system('killall lennard.out') + dx = [dx(1)/dx(1),dx(1)/dx(1)]*dx_max + read(*,*) + end if + if (isnan(dx(1)) .or. isnan(dx(2))) then + print '("NaN NaN Nan Batman! dx =", f18.5, " ", f18.5, " | n ", i4)', dx(1), dx(2), ptr%p%n + print*, "F",ptr%p%F, "v", ptr%p%v + print*, "x", ptr%p%x(1), ptr%p%x(2) + call system('killall lennard.out') + end if + + ptr%p%x(1) = ptr%p%x(1) + dx(1) !dt*ptr%p%v(1) + ptr%p%F(1)*dt**2/(2*m) + ptr%p%x(2) = ptr%p%x(2) + dx(2) !dt*ptr%p%v(2) + ptr%p%F(2)*dt**2/(2*m) + ! print*,'F_0 = ',ptr%p%F(1),ptr%p%F(2), "n", ptr%p%n, "id", id + ! print '("x_0 = [",f10.3,", ",f10.3, "] ", "i, j =", i2, " ", i2, " n ",i2 )',ptr%p%x(1),ptr%p%x(2),i,j, ptr%p%n!"id", id, "n", ptr%p%n + ! print '("v_0 = [",f10.3,", ",f10.3, "] ", "i, j =", i2, " ", i2, " n ",i2 )',ptr%p%v(1),ptr%p%v(2),i,j, ptr%p%n!"id", id, "n", ptr%p%n + ! Ordena as celulas + x = ptr%p%x + ! teste se a partícula pode pular celulas + ! Teste se escapou para a esquerda + if (x(1) <= jcell(j-1)) then + cell = [i,j-1] + !para esquerda para cima + if (x(2) <= icell(i-1)) then + cell = [i-1,j-1] + !para esquerda e para baixo + else if (x(2) > icell(i)) then + cell = [i+1,j-1] + end if + + ! Teste se escapou para a direita + else if (x(1) > jcell(j)) then + cell = [i,j+1] + !para direita para cima + if (x(2) <= icell(i-1)) then + cell = [i-1,j+1] + !para direita e para baixo + else if (x(2) > icell(i)) then + cell = [i+1,j+1] + ! print*, "icell(i)", icell(i), "jcell(j)", jcell(j), "x=", x + ! ! print*, "L PUTAMERDA" + end if + + else if (x(2) <= icell(i-1)) then + ! print*,'!Teste se escapou para cima' + ! print*,x(2) ,'<', icell(i-1) + cell = [i-1, j] + !Teste se escapou para baixo + else if (x(2) > icell(i)) then + ! print*,'!Teste se escapou para baixo' + ! print*,x(2), '>', icell(i) + cell = [i+1,j] + else + cell = [i,j] + end if + + ! print*, "Contint(1)", cont_int(1) + + ! VERIFICAÇÃO SE MUDOUD DE DOMÍNIO + if (cell(1) /= i .or. cell(2) /= j) then !Mudou de celula + ! print*, "MUDOU" + ! if (id == 0) read(*,*) + ! Se a partícula chegar na fronteira do domínio que um processador + ! cuida, então esta partícula é colocada na lista linkada referente + ! a este célula com list_change. Além disso ela precisa ser + ! transferida para o outro processo vizinho para que ele possa + ! calcular sua influência no domínio. + ! Se a partícula ultrapassar o domínio, então ela é transferida + ! para o próximo processo e dealocada da lista com list_remove + ! ! ! print*, "L FORÇA", ptr%p%F, id + if (cell(2) >= domx(2) .and. cell(2) < mesh(1)+2 ) then + ! print*, "L 538", cell(1), cell(2),"part", ptr%p%n,"i,j", i, j + ! 6 elementos + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = [x(1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + ! print*, "id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "id",id,"transferindo para o leste", LT%lstrint_E, "cont", cont_int(3) + ! if (cell(2) > domx(2)) then + ! ! print*, 'affs' + ! ! read(*,*) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_E + + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 end if - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 389 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + if (cell(2) <= domx(1) .and. cell(2) > 1) then + ! print*, "L 554", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "domx", domx + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1), x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "id",id,"transferindo para o oeste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(2) < domx(1)) then + ! ! print*, 'affs' + ! ! read(*,*) + ! ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_W + + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + end if + if (cell(1) >= domy(2) .and. cell(1) < mesh(2)+2) then + ! print*, "L 567", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(1) > domy(1)) then + ! ! read(*,*) + ! ! print*, 'affs' + ! ! call list_remove(previous_node) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_N(cont_int(1)+1:cont_int(1)+4), "cont", cont_int(1) + ! ! node => list_next(previous_node) + ! ! else + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! end if + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + end if + if (cell(1) <= domy(1) .and. cell(1) > 1) then + ! print*, "L 580", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L 583" + ! print*, "id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! if (cell(1) < domy(1)) then + ! ! read(*,*) + ! ! print*, 'affs' + ! ! call list_remove(previous_node) + ! print*, "Removido", ptr%p%n, "de malha", i,j + ! print*, "Irá para", LT%lstrint_S(cont_int(2)+1:cont_int(2)+4), "cont", cont_int(1) + ! ! node => list_next(previous_node) + ! ! else + ! ! !!!! ! print*, "L 590" + ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) + ! ! node => list_next(previous_node) + ! ! !!!! ! print*, "L 624" + ! end if + cont_db(2) = cont_db(2) + 6 + cont_int(2) = cont_int(2) + 4 + end if + + ! CASO PERIODICO + + if (east == 'p' .and. cell(2) == mesh(1)+2 ) then + ! print*, "L 538", cell(1), cell(2),"part", ptr%p%n,"i,j", i, j + ! 6 elementos + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1)- jcell(mesh(1)+1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + ! print*, "id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1), 2,ptr%p%n,ptr%p%grupo] + + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + end if + if (west == 'p' .and. cell(2) == 1) then + ! print*, "L 554", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "domx", domx + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = & + [x(1)+ jcell(mesh(1)+1), x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),mesh(1)+1,ptr%p%n,ptr%p%grupo] + + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 + end if + if (north == 'p' .and. cell(1) == mesh(2)+2) then + ! print*, "L 567", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + ! print*, "id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = & + [x(1),x(2) - icell(mesh(2)+1), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [2,cell(2),ptr%p%n,ptr%p%grupo] + + cont_db(2) = cont_db(2) + 6 + cont_int(2) = cont_int(2) + 4 + end if + if (south == 'p' .and. cell(1) == 1) then + ! print*, "L 580", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j + LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = & + [x(1),icell(mesh(2)+1) + x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [mesh(2)+1,cell(2),ptr%p%n,ptr%p%grupo] + + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + end if + + !!! FIM DA VERIFICAÇÃO SE MUDOU DE DOMÍNIO !!! + + ! Antes aqui tinha um else para o caso de não mudar de processo. + ! Removi todos os call list_remove e deixei que a partícula fosse para + ! a área emprestada que será limpa na próxima execução de call comp_x. Isto pois a célula precisa + ! sentir a presença da partícula que ela perdeu mas ficou ao lado + ! print*, "mudou2" + ! if (id == 0) read(*,*) + + call list_change(previous_node,malha(cell(1),cell(2))%list) + node => list_next(previous_node) + + else !Não mudou de celula + ! Mas tem que avisar quais ainda estão na + ! vizinhança + + ! print*, "NAO mudou cell",cell, "domx", domx, "domy", domy + if (cell(2) == domx(2)) then + !!! ! print*, "L 538", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + ! 6 elementos + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = [x(1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 + + else if (cell(2) == domx(1)) then + !!! ! print*, "L 554", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo para o oeste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + end if + + if (cell(1) == domy(1)) then + !!! ! print*, "L 567", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L 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 + + elseif (cell(1) == domy(2)) then + !!! ! print*, "L 580", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + end if + + ! CASO PERIODICO + + if (east == 'p' .and. cell(2) == mesh(1)+1) then + !!! ! print*, "L 538", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + ! 6 elementos + LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1)- jcell(mesh(1)+1),x(2), & + ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + ! 4 elementos + LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),1,ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo para o oeste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(4) = cont_db(4) + 6 + cont_int(4) = cont_int(4) + 4 + + else if (west == 'p' .and. cell(2) == 2) then + !!! ! print*, "L 554", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_E(cont_db(3)+1:cont_db(3)+6) = & + [x(1)+ jcell(mesh(1)+1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),mesh(1)+2,ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(3) = cont_db(3) + 6 + cont_int(3) = cont_int(3) + 4 + end if + + if (south == 'p' .and. cell(1) == mesh(2)+1) then + !!! ! print*, "L 567", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_N(cont_db(2)+1:cont_db(2)+6) = & + [x(1),icell(mesh(2)+1) + x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_N(cont_int(2)+1:cont_int(2)+4) = [1,cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + cont_db(1) = cont_db(1) + 6 + cont_int(1) = cont_int(1) + 4 + + elseif (north == 'p' .and. cell(1) == 2) then + !!! ! print*, "L 580", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_S(cont_db(1)+1:cont_db(1)+6) = & + [x(1),x(2) - icell(mesh(2)+1), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_S(cont_int(1)+1:cont_int(1)+4) = [mesh(2)+2,cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L 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 + end if + + previous_node => node + node => list_next(node) + ! print*, "L740 FORÇA", ptr%p%F, id + end if + if (associated(node)) ptr = transfer(list_get(node), ptr) + end do + end do + end do + + ! transferir os termos que estão no encontro de 4 processos + ! print*, "L 794", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + if (np > 1 .and. sum(ids(5:8)) > -4) then + do i = 5,8 + !identifica qual o processo na diagonal + if (ids(i) > -1) then + destD = i + end if + end do + ! ! print*, "L 787", ids, "D", destD + ! Aqui está com suporte a 4 processadores por enquanto + ! Primeiro transfere para as celulas emprestadas + + if (domy(1) /= 1) i = domy(1) + if (domy(2) /= mesh(2)+2) i = domy(2) + if (domx(1) /= 1) j = domx(1) + if (domx(2) /= mesh(1) + 2) j = domx(2) + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + ! if (associated(node)) ptr = transfer(list_get(node), ptr) !para que isto? + ! if (associated(node)) print*, "transferencia diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + ! Segundo transfere para as celulas do domínio caso a partícula tenha mudado de processo + if (domy(1) /= 1) i = domy(1) -1 + if (domy(2) /= mesh(2)+2) i = domy(2) +1 + if (domx(1) /= 1) j = domx(1) -1 + if (domx(2) /= mesh(1) + 2) j = domx(2) +1 + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "mudança diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"mudando diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do + !Caso Periódico. Só será necessário se os domínios dos processos formarem região 2x2 + if (wall == 'pppp') then + ! Primeiro transfere para as celulas emprestadas (periodico) + ! O x aqui tem outra função + !extremidades arestas norte e sul + if (domy(1) == 1) then + i = domy(1) + 1 + x = [0.0_dp, icell(mesh(2)+1)] + cell = [mesh(2)+2, 1] + end if + if (domy(2) == mesh(2)+2) then + i = domy(2) - 1 + x = [0.0_dp, -icell(mesh(2)+1)] + cell = [1, 1] + end if + if (domx(1) /= 1) j = domx(1) + if (domx(2) /= mesh(1) + 2) j = domx(2) + cell(2) = j - !Células ao redor !i é linha, j é coluna - if (i /= mesh(2)+2) then ! se não for a última linha - if (j == mesh(1)+2) then ! se for a última coluna - node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha apenas - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "transferencia diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + !extremidades leste oeste + if (domy(1) /= 1) i = domy(1) + if (domy(2) /= mesh(2)+2) i = domy(2) + if (domx(1) == 1) then + j = 2 + x = [jcell(mesh(1)+2), 0.0_dp] + cell = [i, mesh(1) + 2] + end if + if (domx(2) == mesh(1) + 2) then + j = domx(2) - 1 + x = [-jcell(mesh(1)+1), 0.0_dp] + cell = [i, 1] + end if + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "transferencia diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + !extremidades nas pontas + if (domy(1) == 1) then + i = domy(1) +1 + cell(1) = mesh(2) + 2 + x(2) = icell(mesh(2)+1) + end if + if (domy(2) == mesh(2)+2) then + i = domy(2) - 1 + cell(1) = 1 + x(2) = - icell(mesh(2)+1) + end if + if (domx(1) == 1) then + j = domx(1)+1 + cell(2) = mesh(1) + 2 + x(1) = jcell(mesh(1)+1) + end if + if (domx(2) == mesh(1) + 2) then + j = domx(2)-1 + cell(2) = 1 + x(1) = - jcell(mesh(1)+2) + end if + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "transferencia diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"transferindo diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + ! Segundo transfere para as celulas do domínio caso a partícula tenha mudado de processo + ! (periodico) + + ! norte e sul + + if (domy(1) == 1) then + i = domy(1) + x = [0.0_dp, icell(mesh(2)+1)] + cell = [mesh(2)+1, 1] + end if + if (domy(2) == mesh(2)+2) then + i = domy(2) + x = [0.0_dp, -icell(mesh(2)+1)] + cell = [1, 0] + end if + if (domx(1) /= 1) j = domx(1) + if (domx(2) /= mesh(1) + 2) j = domx(2) + cell(2) = j + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "mudança diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"mudando diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + !extremidades leste oeste + if (domy(1) /= 1) i = domy(1) + if (domy(2) /= mesh(2)+2) i = domy(2) + if (domx(1) == 1) then + j = 1 + x = [jcell(mesh(1)+1), 0.0_dp] + cell = [i, mesh(1) + 2] + end if + if (domx(2) == mesh(1) + 2) then + j = domx(2) + x = [-jcell(mesh(1)+1), 0.0_dp] + cell = [i, 2] + end if + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "mudança diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"mudando diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + !extremidades nas pontas + if (domy(1) == 1) then + i = domy(1) + cell(1) = mesh(2) + 2 + x(2) = icell(mesh(2)+1) + end if + if (domy(2) == mesh(2)+2) i = domy(2) + i = domy(2) + cell(1) = 1 + x(2) = - icell(mesh(2)+1) + end if + if (domx(1) == 1) then + j = domx(1) + cell(2) = mesh(1) + 2 + x(1) = jcell(mesh(1)+1) + end if + if (domx(2) == mesh(1) + 2) then + j = domx(2) + cell(2) = 1 + x(1) = - jcell(mesh(1)+2) + end if + + previous_node => malha(i,j)%list + node => list_next(malha(i,j)%list) + ! if (associated(node)) ptr = transfer(list_get(node), ptr) + ! if (associated(node)) print*, "mudança diagonal", i,j, "ID", id + do while(associated(node)) + ptr = transfer(list_get(node), ptr) + x = ptr%p%x ! não sei pq essas linhas não estavam aqui antes. Atentar + cell = [i, j] + ! ! print*, "L 806", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n + LT%lstrdb_D(cont_db(destD)+1:cont_db(destD)+6) = [x(1)+ptr%p%x(1),x(2)+ptr%p%x(2), ptr%p%v(1),ptr%p%v(2), & + ptr%p%F(1),ptr%p%F(2)] + LT%lstrint_D(cont_int(destD)+1:cont_int(destD)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! print*, "L id",id,"mudando diagonal", [cell(1),cell(2),ptr%p%n], "p/ id", ids(destD) + cont_db(destD) = cont_db(destD) + 6 + cont_int(destD) = cont_int(destD) + 4 + previous_node => node + node => list_next(node) + end do + + ! Fim da transfência para diagonal + end if + + ! print*, "L 854", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, "L SOMA", (ids(1) + ids(2) + ids(3) +ids(4)), "id", id + ! print*, "enviando" + if (np > 1) then !se for paralelo + ! print*, 'L 862 >> id', id + ! call MPI_SEND(BUF, COUNT, DATATYPE, DEST, TAG, COMM, IERROR) + tag = 1 + if (ids(1) >= 0) call MPI_SEND(LT%lstrdb_N, cont_db(1), MPI_DOUBLE_PRECISION, ids(1), tag,MPI_COMM_WORLD, ierr) + if (ids(2) >= 0) call MPI_SEND(LT%lstrdb_S, cont_db(2), MPI_DOUBLE_PRECISION, ids(2), tag,MPI_COMM_WORLD, ierr) + if (ids(3) >= 0) call MPI_SEND(LT%lstrdb_E, cont_db(3), MPI_DOUBLE_PRECISION, ids(3), tag,MPI_COMM_WORLD, ierr) + if (ids(4) >= 0) call MPI_SEND(LT%lstrdb_W, cont_db(4), MPI_DOUBLE_PRECISION, ids(4), tag,MPI_COMM_WORLD, ierr) + if (sum(ids(5:8)) > -4) then + call MPI_SEND(LT%lstrdb_D, cont_db(destD), MPI_DOUBLE_PRECISION, ids(destD), tag,MPI_COMM_WORLD, ierr) + end if + tag = 2 + if (ids(1) >= 0) call MPI_SEND(LT%lstrint_N, cont_int(1), MPI_integer, ids(1), tag,MPI_COMM_WORLD, ierr) + if (ids(2) >= 0) call MPI_SEND(LT%lstrint_S, cont_int(2), MPI_integer, ids(2), tag,MPI_COMM_WORLD, ierr) + if (ids(3) >= 0) call MPI_SEND(LT%lstrint_E, cont_int(3), MPI_integer, ids(3), tag,MPI_COMM_WORLD, ierr) + if (ids(4) >= 0) call MPI_SEND(LT%lstrint_W, cont_int(4), MPI_integer, ids(4), tag,MPI_COMM_WORLD, ierr) + if (sum(ids(5:8)) > -4) then + call MPI_SEND(LT%lstrint_D, cont_int(destD), MPI_integer, ids(destD), tag,MPI_COMM_WORLD, ierr) + end if - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 666 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do - - if (j /= 1) then !se não for a primeira coluna - node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if + cont_db = 0 + cont_int = 0 + ! print*, "L 935 TUDO ENVIADO", id + tag = 1 + if (ids(1) >= 0) then + call MPI_probe(ids(1),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(1) , ierr ) + call MPI_RECV (LT%lstrdb_N, cont_db(1), MPI_DOUBLE_PRECISION,ids(1), tag, MPI_COMM_WORLD, status, ierr) + end if + + if (ids(2) >= 0) then + call MPI_probe(ids(2),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(2) , ierr ) + call MPI_RECV (LT%lstrdb_S, cont_db(2), MPI_DOUBLE_PRECISION,ids(2), tag, MPI_COMM_WORLD, status, ierr) + + end if + if (ids(3) >= 0) then + call MPI_probe(ids(3),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(3) , ierr ) + call MPI_RECV (LT%lstrdb_E, cont_db(3), MPI_DOUBLE_PRECISION,ids(3), tag, MPI_COMM_WORLD, status, ierr) + end if - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 710 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do - end if - else - !interagirá com a próxima linha e coluna, e na diagonal - node => list_next(malha(i,j+1)%list) - - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if + if (ids(4) >= 0) then + call MPI_probe(ids(4),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(4) , ierr ) + call MPI_RECV (LT%lstrdb_W, cont_db(4), MPI_DOUBLE_PRECISION,ids(4), tag, MPI_COMM_WORLD, status, ierr) + end if - x2 = ptrn%p%x - v2 = ptrn%p%v - - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 757 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do + tag = 2 + if (ids(1) >= 0) then + call MPI_probe(ids(1),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(1) , ierr ) + call MPI_RECV (LT%lstrint_N, cont_int(1), MPI_integer,ids(1), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(2) >= 0) then + call MPI_probe(ids(2),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(2) , ierr ) + call MPI_RECV (LT%lstrint_S, cont_int(2), MPI_integer,ids(2), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(3) >= 0) then + call MPI_probe(ids(3),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(3) , ierr ) + call MPI_RECV (LT%lstrint_E, cont_int(3), MPI_integer,ids(3), tag, MPI_COMM_WORLD, status, ierr) + end if + if (ids(4) >= 0) then + call MPI_probe(ids(4),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_integer , cont_int(4) , ierr ) + call MPI_RECV (LT%lstrint_W, cont_int(4), MPI_integer,ids(4), tag, MPI_COMM_WORLD, status, ierr) + end if + if (sum(ids(5:8)) > -4) then + tag = 1 + ! RECEBE DA DIAGONAL (suporte para 4 processos) + ! call MPI_probe(,tag, MPI_COMM_WORLD, status, ierr) + call MPI_probe(ids(destD),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_db(destD) , ierr ) + call MPI_RECV (LT%lstrdb_D, cont_db(destD), MPI_DOUBLE_PRECISION,ids(destD), tag, MPI_COMM_WORLD, status, ierr) + tag = 2 + ! Recebe INT da DIAGONAL. Aqui está com suporte para 4 processadores + call MPI_probe(ids(destD),tag, MPI_COMM_WORLD, status, ierr) + CALL MPI_Get_count ( status , MPI_DOUBLE_PRECISION , cont_int(destD) , ierr ) + call MPI_RECV (LT%lstrint_D, cont_int(destD), MPI_DOUBLE_PRECISION,ids(destD), tag, MPI_COMM_WORLD, status, ierr) - node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 798 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do - - node => list_next(malha(i+1,j+1)%list) !interagirá com a próxima linha e coluna - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if + ! DIAGONAL, tem que alterar para o caso com >4 processos + j = destD + do i = 0, cont_db(j)/6 + if (cont_db(j) > 0) then + ! print*, "LA 1108" + ! print*, "D al", id, LT%lstrdb_D(i*6+1:i*6+2) + allocate(ptr%p) + ptr%p%x = LT%lstrdb_D(i*6+1:i*6+2) + ! print*, ptr%p%x + ptr%p%v = LT%lstrdb_D(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_D(i*4+4) + ptr%p%F = LT%lstrdb_D(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_D(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_D(i*4+1), & + LT%lstrint_D(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "D allocado em", LT%lstrint_D(i*4+1), LT%lstrint_D(i*4+2), "id =", id, "posição",ptr%p%x + ! print*, "W allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 + end if + end do + end if - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 841 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do - - if (j /= 1) then !se não for a primeira coluna - node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil_a - end if + ! print*, "cont db", cont_db - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 885 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do - end if - end if + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, "Tudo recebido", id + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, 'ID APOS BARREIRA', id + ! print*, "L 1033" + do j = 1,4 + do i = 0, cont_db(j)/6 + !! ! print*, "L 840" , i, j, "id", id, "cont_db/6", cont_db(j)/6 + if (j == 1 .and. cont_db(1) > 0 ) then + ! print*, "LA 1133" + allocate(ptr%p) + + ptr%p%x = LT%lstrdb_N(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_N(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_N(i*4+4) + ptr%p%F = LT%lstrdb_N(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_N(i*4+3) !identidade da partícula importante para imprimir + ! !print*, 'LLLL', LT%lstrint_N(i*4+1), LT%lstrint_N(i*4+2) + call list_insert(malha(LT%lstrint_N(i*4+1), & + LT%lstrint_N(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "N allocado em", LT%lstrint_N(i*4+1), LT%lstrint_N(i*4+2), "id =", id + ! print*, "N allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 - else ! se for a última lina, só interage com a celua ao lado - node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna - do while (associated(node)) - ptrn = transfer(list_get(node), ptrn) !outra particula selecionada - sigma_b = propriedade(ptrn%p%grupo)%sigma - epsil_b = propriedade(ptrn%p%grupo)%epsilon - rs2 = propriedade(ptrn%p%grupo)%rs - ! Lorenz-Betherlot rule for mixing epsilon sigma - if (sigma_a > sigma_b) then - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else if (sigma_a < sigma_b) then - rcut = r_cut*sigma_b + (rs1 + rs2) - sigma = 0.5*(sigma_a + sigma_b) - epsil = sqrt(epsil_a*epsil_b) - else - rcut = r_cut*sigma_a + (rs1 + rs2) - sigma = sigma_a - epsil = epsil__a - end if - - x2 = ptrn%p%x - v2 = ptrn%p%v - r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) - coss = (x1(1)-x2(1))/r - sine = (x1(2)-x2(2))/r - r = r - rs1 - rs2 !raio - ! print*, "L 931 r", r, "id",id - if (r <= rcut) then - aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & - [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - end if - node => list_next(node) ! próxima partícula da célula - end do end if - node => next_node + if (j == 2 .and. cont_db(2) > 0) then + ! print*, "LA 1151", id + allocate(ptr%p) + ! print*, "L 859 <<<", id + ptr%p%x = LT%lstrdb_S(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_S(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_S(i*4+4) + ptr%p%F = LT%lstrdb_S(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_S(i*4+3) !identidade da partícula importante para imprimir + ! print*, "L 865 <<<", LT%lstrint_S(i*6+1:i*6+6) + call list_insert(malha(LT%lstrint_S(i*4+1), & + LT%lstrint_S(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "S allocado em", LT%lstrint_S(i*4+1), LT%lstrint_S(i*4+2), "id =", id + ! print*, "S allocado F = ",ptr%p%F + !! ! print*, "L 869 <<<" + cont_db(j) = cont_db(j) - 6 + end if + if (j == 3 .and. cont_db(3) > 0) then + ! print*, "LA 1170" + allocate(ptr%p) + ptr%p%x = LT%lstrdb_E(i*6+1:i*6+2) + ptr%p%v = LT%lstrdb_E(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_E(i*4+4) + ptr%p%F = LT%lstrdb_E(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_E(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_E(i*4+1), & + LT%lstrint_E(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "E allocado em", LT%lstrint_E(i*4+1), LT%lstrint_E(i*4+2), "x=",ptr%p%x + ! print*, "E allocado", LT%lstrint_E + cont_db(j) = cont_db(j) - 6 + end if + if (j == 4 .and. cont_db(4) > 0 ) then + ! print*, "LA 1185" + allocate(ptr%p) + ptr%p%x = LT%lstrdb_W(i*6+1:i*6+2) + ! print*, ptr%p%x + ptr%p%v = LT%lstrdb_W(i*6+3:i*6+4) + ptr%p%grupo = LT%lstrint_W(i*4+4) + ptr%p%F = LT%lstrdb_W(i*6+5:i*6+6) + ! ! print*, "L FORÇÆ", ptr%p%F + ptr%p%n = LT%lstrint_W(i*4+3) !identidade da partícula importante para imprimir + call list_insert(malha(LT%lstrint_W(i*4+1), & + LT%lstrint_W(i*4+2))%list, data=transfer(ptr, list_data)) + ! print*, ptr%p%n, "W allocado em", LT%lstrint_W(i*4+1), LT%lstrint_W(i*4+2), "id =", id + ! print*, "W allocado F = ",ptr%p%F + cont_db(j) = cont_db(j) - 6 + end if end do - end do - end do - !print*, 'Linha 143' - ! if (id == 0) read(*,*) - ! call MPI_barrier(MPI_COMM_WORLD, ierr) - end subroutine comp_FT_thermo - - ! atualiza posições - subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) + end do + ! print*, "Tudo Alocado" + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end if + + ! print*, "L 1115 fim comp_x", id + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end subroutine comp_x + + subroutine comp_xT(icell,jcell,malha,N,mesh,propriedade,Tor,theta,omega, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) use linkedlist use mod1 use data @@ -2310,20 +3295,42 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, character(1) :: north, south, east, west type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha - ! integer, intent(inout) :: mic(:) integer :: i,j,k, cell(2), status(MPI_STATUS_SIZE), count, destD,dx1,dx2,dy1,dy2 integer, intent(in) :: N,mesh(2),domx(2),domy(2) type(list_t), pointer :: node, previous_node real(dp), intent(in) :: dt, t, dx_max + real(dp), intent(inout), dimension(:) :: theta + real(dp), intent(in), dimension(:) :: Tor, omega character(4), intent(in) :: wall type(data_ptr) :: ptr - real(dp) :: x(2),m, dx(2) + logical :: auxl ! variavel auxiliar lógica + real(dp) :: x(2),m, dx(2),rs + real(dp), save :: MI !Mi momento de inercia (1/2)m*r^2 para discos real(dp), intent(in) :: icell(:), jcell(:) integer ( kind = 4 ), intent(in) :: ids(8), id, np integer ( kind = 4 ) :: tag,cont_db(8),cont_int(8) type(lstr) :: LT ! real(dp),intent(inout) :: celula(:,:) - ! IDS correspondem às celulas [N,S,E,W] + ! IDS correspondem às celulas [N,S,E,W] + + ! o raio sólido será o da primeira partícula com raio sólido encontrado + ! de acordo com a ordem no settings.ini + if (t < 2*dt) then + rs = 0 + i = 0 + do while (rs == 0) + i = i+1 + rs = propriedade(i)%rs + end do + MI = 0.5*propriedade(i)%m*rs**2 + print*, "Mi calculado!",id, mi + end if + + if (t > 2*dt .and. t < 6*dt) print*, "MI,id", mi, id + + !Computamos a variação de theta + theta = theta + dt*omega + Tor*dt**2/(2*MI) + cont_db = 0 cont_int = 0 @@ -2464,6 +3471,8 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, ! node => list_next(node) end do end if + + ! print*, "L 3464" ! print*, "cont_int", cont_int ! if (id == 0) read(*,*) ! call MPI_barrier(MPI_COMM_WORLD, ierr) @@ -2476,19 +3485,27 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, if (dy2 == mesh(2)+2) dy2 = mesh(2)+1 if (dx2 == mesh(1)+2) dx2 = mesh(1)+1 + + do i = dy1, dy2 do j = dx1,dx2 previous_node => malha(i,j)%list node => list_next(malha(i,j)%list) ! ptr%p%flag = true significa que ainda não foi computado para evitar loop infinito - + if (associated(node)) then ptr = transfer(list_get(node), ptr) + auxl = ptr%p%flag + ! print*, "L 3481 associado" end if + + ! print*, "L 3489", associated(node) - do while (associated(node) .and. ptr%p%flag) + do while (associated(node) .and. auxl) + ! print*, "L PART", i,j, "id", id m = propriedade(ptr%p%grupo)%m + ! print*, "L 3496", m ! aqui vemos se já passou o tempo que as posições ficaram travadas. if (propriedade(ptr%p%grupo)%x_lockdelay > t) ptr%p%flag = .false. @@ -2514,11 +3531,9 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, call system('killall lennard.out') end if + ptr%p%x(1) = ptr%p%x(1) + dx(1) !dt*ptr%p%v(1) + ptr%p%F(1)*dt**2/(2*m) ptr%p%x(2) = ptr%p%x(2) + dx(2) !dt*ptr%p%v(2) + ptr%p%F(2)*dt**2/(2*m) - ! print*,'F_0 = ',ptr%p%F(1),ptr%p%F(2), "n", ptr%p%n, "id", id - ! print '("x_0 = [",f10.3,", ",f10.3, "] ", "i, j =", i2, " ", i2, " n ",i2 )',ptr%p%x(1),ptr%p%x(2),i,j, ptr%p%n!"id", id, "n", ptr%p%n - ! print '("v_0 = [",f10.3,", ",f10.3, "] ", "i, j =", i2, " ", i2, " n ",i2 )',ptr%p%v(1),ptr%p%v(2),i,j, ptr%p%n!"id", id, "n", ptr%p%n ! Ordena as celulas x = ptr%p%x ! teste se a partícula pode pular celulas @@ -2581,48 +3596,20 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, ! 4 elementos ! print*, "id",id,"transferindo para o leste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] LT%lstrint_E(cont_int(3)+1:cont_int(3)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "id",id,"transferindo para o leste", LT%lstrint_E, "cont", cont_int(3) - ! if (cell(2) > domx(2)) then - ! ! print*, 'affs' - ! ! read(*,*) - ! print*, "Removido", ptr%p%n, "de malha", i,j - ! print*, "Irá para", LT%lstrint_E - cont_db(3) = cont_db(3) + 6 cont_int(3) = cont_int(3) + 4 end if if (cell(2) <= domx(1) .and. cell(2) > 1) then - ! print*, "L 554", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j - ! print*, "domx", domx LT%lstrdb_W(cont_db(4)+1:cont_db(4)+6) = [x(1), x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_W(cont_int(4)+1:cont_int(4)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "id",id,"transferindo para o oeste", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! if (cell(2) < domx(1)) then - ! ! print*, 'affs' - ! ! read(*,*) - ! ! print*, "Removido", ptr%p%n, "de malha", i,j - ! print*, "Removido", ptr%p%n, "de malha", i,j - ! print*, "Irá para", LT%lstrint_W cont_db(4) = cont_db(4) + 6 cont_int(4) = cont_int(4) + 4 end if if (cell(1) >= domy(2) .and. cell(1) < mesh(2)+2) then - ! print*, "L 567", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j - ! print*, "id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! if (cell(1) > domy(1)) then - ! ! read(*,*) - ! ! print*, 'affs' - ! ! call list_remove(previous_node) - ! print*, "Removido", ptr%p%n, "de malha", i,j - ! print*, "Irá para", LT%lstrint_N(cont_int(1)+1:cont_int(1)+4), "cont", cont_int(1) - ! ! node => list_next(previous_node) - ! ! else - ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) - ! ! node => list_next(previous_node) - ! end if + cont_db(1) = cont_db(1) + 6 cont_int(1) = cont_int(1) + 4 end if @@ -2630,21 +3617,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, ! print*, "L 580", cell(1), cell(2), "part", ptr%p%n, "i,j", i, j LT%lstrdb_S(cont_db(2)+1:cont_db(2)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_S(cont_int(2)+1:cont_int(2)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! ! print*, "L 583" - ! print*, "id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! if (cell(1) < domy(1)) then - ! ! read(*,*) - ! ! print*, 'affs' - ! ! call list_remove(previous_node) - ! print*, "Removido", ptr%p%n, "de malha", i,j - ! print*, "Irá para", LT%lstrint_S(cont_int(2)+1:cont_int(2)+4), "cont", cont_int(1) - ! ! node => list_next(previous_node) - ! ! else - ! ! !!!! ! print*, "L 590" - ! ! call list_change(previous_node,malha(cell(1),cell(2))%list) - ! ! node => list_next(previous_node) - ! ! !!!! ! print*, "L 624" - ! end if + cont_db(2) = cont_db(2) + 6 cont_int(2) = cont_int(2) + 4 end if @@ -2742,7 +3715,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, !!! ! print*, "L 580", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] cont_db(1) = cont_db(1) + 6 cont_int(1) = cont_int(1) + 4 end if @@ -2784,7 +3757,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, LT%lstrdb_S(cont_db(1)+1:cont_db(1)+6) = & [x(1),x(2) - icell(mesh(2)+1), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_S(cont_int(1)+1:cont_int(1)+4) = [mesh(2)+2,cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "L id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L 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 end if @@ -2793,7 +3766,10 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, node => list_next(node) ! print*, "L740 FORÇA", ptr%p%F, id end if - if (associated(node)) ptr = transfer(list_get(node), ptr) + if (associated(node)) then + ptr = transfer(list_get(node), ptr) + auxl = ptr%p%flag + end if end do end do end do @@ -3283,7 +4259,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, ! print*, "L 1115 fim comp_x", id ! if (id == 0) read(*,*) ! call MPI_barrier(MPI_COMM_WORLD, ierr) - end subroutine comp_x + end subroutine comp_xT ! atualiza posições subroutine comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,domy,wall,id,np,xih,xic,cold_cells,hot_cells) @@ -3695,7 +4671,7 @@ subroutine comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,d !!! ! print*, "L 580", cell(1), cell(2), domy(1), domy(2), "part", ptr%p%n LT%lstrdb_N(cont_db(1)+1:cont_db(1)+6) = [x(1),x(2), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_N(cont_int(1)+1:cont_int(1)+4) = [cell(1),cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L id",id,"transferindo para o norte", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] cont_db(1) = cont_db(1) + 6 cont_int(1) = cont_int(1) + 4 end if @@ -3737,7 +4713,7 @@ subroutine comp_x_thermo(icell,jcell,malha,N,mesh,propriedade,t,dt,ids,LT,domx,d LT%lstrdb_S(cont_db(1)+1:cont_db(1)+6) = & [x(1),x(2) - icell(mesh(2)+1), ptr%p%v(1),ptr%p%v(2), ptr%p%F(1),ptr%p%F(2)] LT%lstrint_S(cont_int(1)+1:cont_int(1)+4) = [mesh(2)+2,cell(2),ptr%p%n,ptr%p%grupo] - ! print*, "L id",id,"transferindo para o sul", [cell(1),cell(2),ptr%p%n,ptr%p%grupo] + ! ! print*, "L 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 end if @@ -4276,6 +5252,63 @@ subroutine comp_v(malha,mesh,dt,t,propriedade,domx,domy) end do end subroutine comp_v + subroutine comp_vT(malha,mesh,dt,t,omega,Tor,propriedade,domx,domy,id) + use linkedlist + use mod1 + use data + use mod0 + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + integer :: i,j, bsh(4),id + integer, intent(in) :: mesh(2),domx(2),domy(2) + type(list_t), pointer :: node + real(dp), save :: MI + real(dp) :: rs + real(dp), intent(in) :: dt, t + real(dp),dimension(:), intent(inout) :: omega, Tor(:) + type(data_ptr) :: ptr + + ! velocidade de rotação das partículas que giram + if (t < 2*dt) then + rs = 0 + i = 0 + do while (rs == 0) + i = i+1 + rs = propriedade(i)%rs + end do + MI = 0.5*propriedade(i)%m*rs**2 + print*, "Mi calculado compV!",id, mi + end if + + if (t > 2*dt .and. t < 4*dt) print*, "MI,id compV", mi, id + + omega = omega + Tor*dt/(2*MI) + Tor = 0 ! reinicia os torques para a próxima iteração + + bsh = [0,0,0,0] + if (domx(1) == 1) bsh(1) = 1 + if (domx(2) == mesh(1)+2) bsh(2) = -1 + if (domy(1) == 1) bsh(3) = 1 + if (domy(2) == mesh(2)+2) bsh(4) = -1 + + do i = (domy(1)+bsh(3)), (domy(2)+bsh(4)) + do j = (domx(1)+bsh(1)), (domx(2)+bsh(2)) + node => list_next(malha(i,j)%list) + do while (associated(node)) + ptr = transfer(list_get(node), ptr) + ! print*, "n", ptr%p%n, "F =",ptr%p%F, "v =",ptr%p%v + if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + 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 + ptr%p%F = [0,0] + node => list_next(node) + end do + end do + end do + end subroutine comp_vT + 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 @@ -4899,15 +5932,17 @@ program main use, intrinsic :: iso_fortran_env, only: real32 implicit none + + real(dp),parameter :: PI = 3.1415927, kb = 1.38064852E-23 ! 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 :: N,Ntype,i=1, nimpre,j = 1, ii, quant = 0,mesh(2), cont = 1, aux1 = 0,cont2 = 1,domx(2), domy(2), aux3, num_rot, pa integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init, cpu_countrate,ic1, force_lingradT, aux4, domx_aux(2) 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,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 + real(dp), dimension(:), allocatable :: aux5,omega,theta,Tor,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,finish,Td,vd(2) + real(dp) :: GField(2), temp_Td(3), dimX, dimY, dx_max, Td_hot, Td_cold, xih, xic, Mc, rs, pr(5) + integer,allocatable :: partlst(:),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 type(list_t), pointer :: node @@ -4935,6 +5970,7 @@ program main start = 0 xih = 0 xic = 0 + pr = [0,0,0,0,0] nan = IEEE_VALUE(nan, IEEE_QUIET_NAN) @@ -5048,6 +6084,7 @@ program main printstep = ((temp_Td(2)-temp_Td(1))/dt)/temp_Td(3) if (printstep < 1) printstep = 1 !aloca as Variáveis com base no número de partículas + allocate(v(N,2),x(N,2),grupo(N),interv(nimpre),interv_Td(nint(printstep)), & jcell(mesh(2)+1),icell(mesh(1)+1),celula(N,2),propriedade(Ntype)) 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), & @@ -5057,7 +6094,7 @@ program main allocate(mic(N,2),mic_trf(2*N)) mic = 0 end if - + ! allocate(nxv(N*5),nxv_send(N*5)) aux2 = dimX/mesh(1) @@ -5093,13 +6130,16 @@ program main call CFG_add(my_cfg,particle//"%fric_term",1.0_dp, & "Friction term") call CFG_add(my_cfg,particle//"%ismolecule",.true., & - "Is Molecule") + "Is Molecule") + call CFG_add(my_cfg,particle//"%pr",(/0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), & + "Constants to modify the lennard-jones potential of a particle (not molecule).") end do dx_max = 10*dimx !pra definir critério de estabilidade no uso de malha do i = 0,(Ntype-1) write(particle,'(a,i0)') 'par_',i call CFG_get(my_cfg, particle//"%quantidade", quant) + propriedade(i+1)%quant = quant call CFG_get(my_cfg, particle//"%nome", nome) part_nomes(i+1)%str = nome call CFG_get(my_cfg, particle//"%x", arquivox) @@ -5111,7 +6151,9 @@ program main call CFG_get(my_cfg, particle//"%rs", propriedade(i+1)%rs) call CFG_get(my_cfg, particle//"%fric_term",propriedade(i+1)%fric_term) call CFG_get(my_cfg, particle//"%ismolecule",propriedade(i+1)%ismolecule) + if (.not. propriedade(i+1)%ismolecule) call CFG_get(my_cfg, particle//"%pr",pr) + ! le o arquivo com posições if (dx_max < propriedade(i+1)%sigma*rcut/2) then dx_max = propriedade(i+1)%sigma*rcut/2 @@ -5136,6 +6178,7 @@ program main end if close(20) end do + dx_max = dx_max**2 ! mostra informações lidas if (id == 0) then @@ -5368,8 +6411,300 @@ program main end if if (id == 0) call system('mkdir temp2') end if + print*, "INICIANDO" + !! CASE OF RUGGED LENNARD JONES THAT CAN MAKE THE PARTICLE ROTATE + if (pr(1)+pr(2) > 0) then ! PARTICLE ROTATES + if (id == 0) then + deallocate(x,v) + allocate(x(N,3),v(N,3)) + end if + rs = 0 + ii = 0 + pa = 0 + do while (rs == 0) + ii = ii+1 + rs = propriedade(ii)%rs + pa = pa + propriedade(ii)%quant + end do + pa = pa - propriedade(ii)%quant + num_rot = propriedade(ii)%quant + allocate(omega(num_rot), Tor(num_rot), theta(num_rot),partlst(num_rot), aux5(num_rot*3)) + Tor = 0 + theta = 0 + partlst = 0 + ! Calcula a rotação da partícula. Velocity scaling. + do while (t_fim > t) + ! COMPF + print*, "L 6404" + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + call comp_FT(GField,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,t,dt,partlst) !altera Força + ! IDS são os ids das regiões vizinhas, vetor(4) int posições, + ! DOMX e DOMY são vetores(2) int com o domínio (quat. de celulas) + ! que o processo vai cuidar (subdivisões) + print*, "L 6408" + theta = theta*partlst + omega = omega*partlst + + ! COMP X + call comp_xT(icell,jcell,malha,N,mesh,propriedade,Tor,theta,omega, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) ! altera posição + print*, "L 6414" + do ii = 1,num_rot + if ( partlst(ii) > -1 ) then ! particula não está no domínio + theta(ii) = 0 + omega(ii) = 0 + end if + end do + ! compartilha valores de theta, torque e omega (pos, torque, vel angular) + ! num_rot, numero de partícuals que giram + ! print*, "L 6426" + + if (np > 1) then + if (id /= 0) then + aux5 = [theta,omega,Tor] + ! aux5(1:num_rot) = theta ! vetor auxiliar pra o mpi, theta + ! aux5(num_rot+1:2*num_rot) = omega ! " " ", omega + ! aux5(2*num_rot+1:3*num_rot) = Tor ! " " ", Tor + tag = id + call MPI_SEND(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD, ierr) + else ! processo raiz + ! recebe dos processos servos. Aqui vamos somar. Garantiremos que os valores de theta e omega + ! nos processos onde a partícula não está sejam = 0 + call MPI_RECV(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 1, 1, MPI_COMM_WORLD, status, ierr) + theta = theta + aux5(1:num_rot) + omega = omega + aux5(num_rot+1:2*num_rot) + Tor = Tor + aux5(2*num_rot+1:3*num_rot) + call MPI_RECV(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 2, 2, MPI_COMM_WORLD, status, ierr) + theta = theta + aux5(1:num_rot) + omega = omega + aux5(num_rot+1:2*num_rot) + Tor = Tor + aux5(2*num_rot+1:3*num_rot) + call MPI_RECV(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 3, 3, MPI_COMM_WORLD, status, ierr) + theta = theta + aux5(1:num_rot) + omega = omega + aux5(num_rot+1:2*num_rot) + Tor = Tor + aux5(2*num_rot+1:3*num_rot) + + ! omega é angulo e portanto periódico, para não estourar vamos pegar o resto da divisão por 2*pi + omega = modulo(omega,2*PI) + ! agora vamos distribuir Tor, omega e theta para todos os processos, todos terão os mesmos valores + aux5 = [theta,omega,Tor] + ! aux5(1:num_rot) = theta! vetor auxiliar pra o mpi, theta + ! aux5(num_rot+1:2*num_rot) = omega ! " " ", omega + ! aux5(2*num_rot+1:3*num_rot) = Tor ! " " ", Tor + + end if + end if + print*, "L 6455", np + if (np > 1) then + call MPI_BCAST(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (id /= 0) then + theta = aux5(1:num_rot) + omega = aux5(num_rot+1:2*num_rot) + Tor = aux5(2*num_rot+1:3*num_rot) + end if + end if + + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + print*, "L 6484" + call walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ! altera posição e malha + print*, "L 6486" + ! COMP F + call comp_FT(GField,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,t,dt,partlst) !altera força + print*, "L 6470" + + ! COMP V + call comp_vT(malha,mesh,dt,t,omega,Tor,propriedade,domx,domy,id) !altera velocidade + print*, "L 6474" + + t = t + dt + cont2 = cont2+1 + if (i == interv_Td(j2)) then + if (force_lingradT > 0 .and. i == interv_Td(j2)) then + ! ao longo da direção x, o gradiente terá [force_lingradT] passos. Teremos (domx(2)-domx(1))/(force_lingradT/subx) celulas por passo + do aux4 = 0, force_lingradT/subx-1 + !cellx por processo ! passos por processo + domx_aux = [1+aux4*(domx(2)-domx(1)) / (force_lingradT/subx), & + (aux4+1)*(domx(2)-domx(1)) / (force_lingradT/subx) ] + if (domx_aux(1) == 1) domx_aux(1) = 2 + if (domx_aux(2) == mesh(1)+2) domx_aux(2) = mesh(1) + 1 + Td = Td_cold + sum(domx_aux)*(Td_hot-Td_cold)/(2*mesh(1)) + ! print*, Td, Td_cold, Td_hot,sum(domx_aux),2*mesh(1) + + call corr_Kglobal(malha,domx_aux,domy,Td,propriedade, np,id,t,.true.) + ! Este .true. indica para ignorar mpi, ou seja, pedaços de outras regiões + ! não serão levados em conta. + end do + end if + + 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 (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 + j2 = j2 + 1 + end if + print*, "L 6535 " + if (i == interv(j)) then + deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, & + LT%lstrint_N, LT%lstrint_S, LT%lstrint_E, LT%lstrint_W, LT%lstrint_D) + allocate(nxv(N*5),nxv_send(N*5)) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, 'MANDANDO PRINTAR', id, i + + call linked2vec(malha,mesh,domx,domy,nxv_send,aux1) + ! ! ! print*, "L nxv_send", nxv_send + !! print*, "LINKED 2 VEC, id", id, "aux1", aux1 + if (np > 1) then ! paralelo + ! call MPI_GATHER(sbuf, scount, MPI_integer, rbuf, rcount, MPI_integer, root, MPI_COMM_WORLD, ierr) + call MPI_GATHER(aux1, 1, MPI_integer, rcounts, 1, MPI_integer, 0, MPI_COMM_WORLD, ierr) + + ! print*, "Aux1 =", aux1, "id =", id + ! if (id == 0) print*, "rcounts =", rcounts + ! !print*, 'L IMPRE >', id + displs = 0 + + if (id == 0) then + do ii = 2,np + displs(ii) = displs(ii-1) + rcounts(ii-1) + end do + ! print*, "displs", displs + end if + ! print*,'E_tot0 = ',(comp_pot(mesh,malha,propriedade,rcut) +comp_K(malha,mesh,propriedade)) + ! MPI_GATHERV (sbuf, scount, stype, rbuf, rcounts, displs, rtype, root, comm, ierr) + call MPI_GATHERV(nxv_send, aux1, MPI_DOUBLE_PRECISION, nxv, rcounts, & + displs, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + else + nxv = nxv_send + end if + print*, "L 6568" + if (id == 0) then + do ii = 0, N-1 + if (ii <= pa) then + x(int(nxv(ii*5+1)),:) = [nxv(ii*5+2),nxv(ii*5+3),0.0_dp] + v(int(nxv(ii*5+1)),:) = [nxv(ii*5+4),nxv(ii*5+5),0.0_dp] + else + x(int(nxv(ii*5+1)),:) = [nxv(ii*5+2),nxv(ii*5+3),theta] + v(int(nxv(ii*5+1)),:) = [nxv(ii*5+4),nxv(ii*5+5),omega] + end if + end do + + call vec2csv(x,N,3,'position',j,t,nimpre,start) + call vec2csv(v,N,3,'velocity',j,t,nimpre,start) + + end if + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + + deallocate(nxv,nxv_send) + + if (print_TC) then + if ((wall(1:2) == 'pp' .or. wall(3:4) == 'pp') .and. id /= 0) then + mic_trf(1:N) = mic(:,1) + mic_trf(N+1:2*N) = mic(:,2) + end if + + nRfu = 0.0_dp + + call comp_pot(mesh,malha,propriedade,rcut,domx,domy, ids, id, t, nRfu) + + if (np > 1) then ! paralelo + aux3 = 1 + do ii = 1,N + if (nRfu(ii*6-5) /= 0) then + nRfu_send(aux3:aux3+5) = nRfu(ii*6-5:ii*6) + aux3 = aux3 + 6 + end if + end do + aux1 = aux3 -1 + + ! call MPI_GATHER(sbuf, scount, MPI_integer, rbuf, rcount, MPI_integer, root, MPI_COMM_WORLD, ierr) + call MPI_GATHER(aux1, 1, MPI_integer, rcounts, 1, MPI_integer, 0, MPI_COMM_WORLD, ierr) + + displs = 0 + + if (id == 0) then + do ii = 2,np + displs(ii) = displs(ii-1) + rcounts(ii-1) + end do + ! print*, "displs", displs + end if + + ! MPI_GATHERV (sbuf, scount, stype, rbuf, rcounts, displs, rtype, root, comm, ierr) + call MPI_GATHERV(nRfu_send, aux1, MPI_DOUBLE_PRECISION, nRfu, rcounts, & + displs, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + end if + + if ((wall(1:2) == 'pp' .or. wall(3:4) == 'pp') .and. np > 1) then + ! tag = 1 + if (id == 1) then + tag = 10 + call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + end if + if (id == 0) then + tag = 10 + ! MPI_RECV(BUF, COUNT, DATATYPE, SOURCE, TAG, COMM, STATUS, IERROR) + call MPI_RECV(mic_trf,2*N,MPI_integer,1,tag,MPI_COMM_WORLD,status,ierr) + mic(:,1) = mic(:,1) + mic_trf(1:N) + mic(:,2) = mic(:,2) + mic_trf(N+1:2*N) + end if + call MPI_barrier(MPI_COMM_WORLD, ierr) + if (id == 2) then + tag = 20 + call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + end if + if (id == 0) then + tag = 20 + call MPI_RECV(mic_trf,2*N,MPI_integer,2,tag,MPI_COMM_WORLD,status,ierr) + mic(:,1) = mic(:,1) + mic_trf(1:N) + mic(:,2) = mic(:,2) + mic_trf(N+1:2*N) + end if + call MPI_barrier(MPI_COMM_WORLD, ierr) + if (id == 3) then + tag = 30 + call MPI_SEND(mic_trf,2*N,MPI_integer,0,tag,MPI_COMM_WORLD,ierr) + end if + if (id == 0) then + tag = 30 + call MPI_RECV(mic_trf,2*N,MPI_integer,3,tag,MPI_COMM_WORLD,status,ierr) + mic(:,1) = mic(:,1) + mic_trf(1:N) + mic(:,2) = mic(:,2) + mic_trf(N+1:2*N) + end if + call MPI_barrier(MPI_COMM_WORLD, ierr) + end if + + 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),kind(0.d0)),real(mic(ii,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) + 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 + ! 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), & + LT%lstrint_N(NMPT*4),LT%lstrint_S(NMPT*4),LT%lstrint_E(NMPT*4),LT%lstrint_W(NMPT*4), LT%lstrint_D(NMPT*4)) + j = j+1 + end if + i = i+1 + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, t + end do - if (termostato_nose_hoover) then + else if (termostato_nose_hoover) then ! TERMOSTATO NOSE HOOVER do while (t_fim > t) !print*, "L 1990" @@ -5384,7 +6719,7 @@ 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 diff --git a/mod0.f90 b/mod0.f90 index d52161b..e4dc22e 100644 --- a/mod0.f90 +++ b/mod0.f90 @@ -15,6 +15,7 @@ module mod0 real(dp) :: x_lockdelay ! só vai poder mudar de posição a partir de t = x_lockdelay real(dp) :: fric_term logical :: ismolecule + integer :: quant end type prop_grupo ! LSTR lista de transferência pro MPI diff --git a/pack.sh b/pack.sh index 02386b5..f266802 100755 --- a/pack.sh +++ b/pack.sh @@ -3,7 +3,7 @@ read -e -p "Enter particle generator file: " VAR1 if [ "$VAR1" = "" ]; then echo "No particle generator given." - tar -cvf lennard_files.tar pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py + tar -cvf lennard_files.tar potential_analysis.py pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py else tar -cvf lennard_files.tar pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard $VAR1 settings.ini verify_settings.py visualizar.py fi diff --git a/rot_par.py b/rot_par.py new file mode 100644 index 0000000..24da5c8 --- /dev/null +++ b/rot_par.py @@ -0,0 +1,80 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.ticker import MaxNLocator + +sigma_ref = 0.341 + +delta = 0.005 +sig0 = 1 +ep0 = 1 +rs = 5 # raio solido +rm = .245 / sigma_ref #raio do atomo +# .245 é a distância entra atomos em estrutura hexagonal rombica + +x = y = np.arange(-(3*sig0+rs), 3*sig0+rs, delta) +X, Y = np.meshgrid(x, y) + +r = np.sqrt(X**2 + Y**2) +np.place(r,r 20, np.nan) + +plt.figure() +plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:]) +plt.title("potential") + + +fig1, ax1 = plt.subplots() +cs1 = ax1.contourf(X, Y, V,levels=20) +cbar = fig1.colorbar(cs1) +ax1.set_title("potential") + +#force + +np.place(Fr,abs(Fr) > 5000, np.nan) + +# plt.figure() +# plt.plot(X[int(len(X)/2),:] ,Fr[int(len(X)/2),:]) +# plt.title("force") + + +cmap = plt.get_cmap('PiYG') +fig2, ax2 = plt.subplots() +cs2 = ax2.contourf(X, Y, Fr,levels=20, cmap=cmap) +cbar = fig2.colorbar(cs2) +ax2.set_title('radial force') + + +np.place(Fg,abs(Fg) > 5000, np.nan) +fig3, ax3 = plt.subplots() +cs3 = ax3.contourf(X, Y, Fg, levels=20, cmap=cmap) +cbar = fig3.colorbar(cs3) +ax3.set_title('tangential force') + + +plt.show() \ No newline at end of file diff --git a/saida.f90 b/saida.f90 index dec7af6..710b2fe 100644 --- a/saida.f90 +++ b/saida.f90 @@ -49,7 +49,7 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) end do else if (d == 7) 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,fmt7) v(i,1),v(i,2),v(i,3),v(i,4),v(i,5),v(i,6),v(i,7) end do end if diff --git a/settings.ini b/settings.ini index 08b77cf..ef7193b 100755 --- a/settings.ini +++ b/settings.ini @@ -3,18 +3,18 @@ # 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 = 100 # quntidade de saidas - N = 4100 # número de partículas + nimpre = 200 # quntidade de saidas + N = 2 # número de partículas Ntype = 2 # número de tipos de partícula - t_fim = 40000 # tempo de execução + t_fim = 20 # tempo de execução nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). dt = 0.001 # passo de tempo - dimX = 10000 # Dimensões da região de cálculo - dimY = 1250 # - mesh = 200 30 # 300 60 # malha(elementos) + dimX = 600 # Dimensões da região de cálculo + dimY = 600 # + mesh = 8 8 # 300 60 # malha(elementos) rcut = 3 # *sigma - wall = 'ppee' # 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 = 'eeee' # 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) @@ -24,7 +24,7 @@ Td_hot = 2 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato Td_cold = 0.8 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato force_lingradT = -1 # (NUMERO MULTIPLO DE 4 RECOMENDADO!) 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 = 1.4 1.4 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') + 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,32 +32,36 @@ # v é velocidade global, entrar dois número correspondentes ao vetor vx vy que serão aplicadas em todas as partículas # 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 +# pr = constantes [A, B, alpha, beta] do lennard jones rugoso que simulará um cristal de atomos. Se o valor for dado, a particula poderá girar. +# o grupo com pr deve ser o último e deverá haver apenas um. Ver em rot_par.py como será o potencial [par_0] - x = 'moleculas8x1_sc.csv' - v = 0 0 # velocidade global + x = 'mole.csv' + v = 5 0 # velocidade global v_file = '%%v_file_0.csv' # velocidade de cada partícula m = 1 nome = 'g1' sigma = 1 # Lennard Jones epsilon = 1 # Lennard Jones - quantidade = 4000 + quantidade = 1 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 = 'particulas8x1_sc.csv' - v = 0 0 # velocidade global + x = 'parti.csv' + v = 5 0 # velocidade global v_file = '%%v_file_1.csv' # velocidade de cada partícula m = 60 nome = 'g2' sigma = 0.94 # Lennard Jones epsilon = 0.428 # Lennard Jones - quantidade = 100 + quantidade = 1 x_lockdelay = 1000 # só vai poder mudar de posição a partir de t = x_lockdelay rs = 5 # 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. + ismolecule = .false. + pr = 0.104105 0.104105 43.72584 43.72584 0.071847 # constantes [A, B, alpha, beta, ph] do lennard jones rugoso que simulará um cristal de atomos + #[par_2] # x = 'p_g.csv' # v = 0 0