pos_e_vel corrigido e outros prob

csv2vtk_particles.py

@@ -30,14 +30,32 @@
 folder = 'result'
 i = 1
 
+#Ler a quantidade de arquivos de settings.
+
+config = configparser.ConfigParser()
+config.read('settings.ini')
+print('Reading settings.ini')
+N = int(config['global']['N'].split()[0])
+nimpre =  int(config['global']['nimpre'].split()[0])
+ntype = int(config['global']['Ntype'].split()[0])
+nimpre_init = int(config['global']['nimpre_init'].split()[0])
 
+if config['global']['print_TC'].split()[0] == '.false.':
+    print_TC = False
+else:
+    print_TC = True
+
+#Checar se existe arquivo
 while aux:
     try:
         os.mkdir(folder)
         aux = False
     except FileExistsError:
         print('Folder {} exists'.format(folder))
-        opt = input('Overwrite? [y/n] ')
+        if nimpre_init == 0:
+            opt = input('Overwrite? [y/n] ')
+        else:
+            opt = input('Append? (starting on {}) [y/n] '.format(nimpre_init))
         if opt == 'y' or 'opt' == 'Y':
             aux =False 
         else:
@@ -49,18 +67,7 @@
 
 print('\nThe results will be saved at {}\n'.format(folder))
 
-#Ler a quantidade de arquivos de settings.
 
-config = configparser.ConfigParser()
-config.read('settings.ini')
-print('Reading settings.ini')
-N = int(config['global']['N'].split()[0])
-nimpre =  int(config['global']['nimpre'].split()[0])
-ntype = int(config['global']['Ntype'].split()[0])
-if config['global']['print_TC'].split()[0] == '.false.':
-    print_TC = False
-else:
-    print_TC = True
 quant = []
 rs = [] # raio sólido
 sigma = []
@@ -74,10 +81,12 @@
 rs = rs + sigma*(2**(1/6))
 a = os.listdir('temp')
 no_out_files = nimpre
-if len(a)/2 < nimpre:
+if len(a)/2 < nimpre-nimpre_init:
     no_out_files = int(len(a)/2)
-    print("Propable incomplete execution. Processing {} files\n.".format(len(a)/2))
-    nimpre = int(len(a)/2)-1
+    print("Propable incomplete execution. Processing {} files".format(len(a)/2))
+    nimpre = int(len(a)/2)-1 + nimpre_init
+    input("nimpre = {}, nimpre_init = {}. OK?\n".format(nimpre,nimpre_init))
+if nimpre_init > 0: print("Starting on step {}".format(nimpre_init))
 
 tipo = np.zeros(N)
 rsol = np.zeros(N)
@@ -104,22 +113,34 @@
 cz = np.zeros(N)
 
 nID = np.linspace(1,N,N)
-
-zip_positions = ZipFile(via+'/'+folder+'/positions.zip','w')
-zip_velocities = ZipFile(via+'/'+folder+'/velocities.zip','w')
-
+if nimpre_init == 0:
+    zip_positions = ZipFile(via+'/'+folder+'/positions.zip','w')
+    zip_velocities = ZipFile(via+'/'+folder+'/velocities.zip','w')
+else:
+    zip_positions = ZipFile(via+'/'+folder+'/positions.zip','a')
+    zip_velocities = ZipFile(via+'/'+folder+'/velocities.zip','a')    
+
+
 print('Converting...')
+position_list = zip_positions.namelist()
+aux1 = 0
+flag = True
 if pbar:
-    bar = progressbar.ProgressBar(max_value=nimpre)
-for fnum in range(0,nimpre+1):
+    bar = progressbar.ProgressBar(max_value=(nimpre-nimpre_init))
+for fnum in range(nimpre_init,nimpre+1):
     with open('temp/position.csv.'+str(fnum),encoding='utf-8') as file_locus:
         csv_lector = csv.reader(file_locus,delimiter = ',')
         i = 0
         for linea in csv_lector:
             x[i] = linea[0]
             y[i] = linea[1]
             i = i+1
-    zip_positions.write(via+'/temp/position.csv.'+str(fnum), 'position.csv.'+str(fnum))
+    while flag:
+        if any('position.csv.'+str(fnum+aux1) in s for s in position_list):
+            aux1 += 1
+        else:
+            flag = False
+    zip_positions.write(via+'/temp/position.csv.'+str(fnum), 'position.csv.'+str(fnum+aux1))
     # shutil.move(via+'/temp/position.csv.'+str(fnum),via+'/'+folder+'/position.csv.'+str(fnum)) 
     with open('temp/velocity.csv.'+str(fnum),encoding='utf-8') as file_velocitas:
         csv_lector = csv.reader(file_velocitas,delimiter = ',')
@@ -129,7 +150,7 @@
             vy[i] = linea[1]
 
             i = i+1
-    zip_velocities.write(via+'/temp/velocity.csv.'+str(fnum),'velocity.csv.'+str(fnum))
+    zip_velocities.write(via+'/temp/velocity.csv.'+str(fnum),'velocity.csv.'+str(fnum+aux1))
     # shutil.move(via+'/temp/velocity.csv.'+str(fnum),via+'/'+folder+'/velocity.csv.'+str(fnum))
 
     fin = 0
@@ -145,19 +166,18 @@
         tipos = tipo[ini:fin]
         rsols = rsol[ini:fin]
         nIDs = nID[ini:fin]
-        pointsToVTK(via+'/'+folder +'/grupo'+ str(grupo) + '_' +str(fnum), xs, ys, zs, data = {"Vx" : vxs, "Vy" : vys, "Tipo" : tipos, "raio_solido" : rsols, "nID" : nIDs })       
+        pointsToVTK(via+'/'+folder +'/grupo'+ str(grupo) + '_' +str(fnum+aux1), xs, ys, zs, data = {"Vx" : vxs, "Vy" : vys, "Tipo" : tipos, "raio_solido" : rsols, "nID" : nIDs })       
         grupo += 1
 
     if pbar:
-        bar.update(fnum)
+        bar.update(fnum-nimpre_init)
     if fnum == nimpre+1:
         time.sleep(0.1)
 
 
 zip_positions.close()
 zip_velocities.close()
-
-shutil.rmtree('temp')
+if nimpre > nimpre_init: shutil.rmtree('temp')
 
 if print_TC:
     a = os.listdir('temp2')
@@ -199,4 +219,4 @@
 
             #i = i+1
 
-    #shutil.move(via+'/cell.csv.'+str(fnum),via+'/'+folder+'/cell.csv.'+str(fnum))  
+    #shutil.move(via+'/cell.csv.'+str(fnum),via+'/'+folder+'/cell.csv.'+str(fnum))  

debug_lennard.sh

@@ -1,2 +1,22 @@
 #!/bin/bash
-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
+read -p "[d]ebug with mpiifort or [c]ompile with OpenMPI+gfortran: " option
+
+if [option = 'c']
+then
+    echo "compiling with mpifort"
+    /usr/lib64/mpi/gcc/openmpi/bin/mpifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -o lennard.out
+    echo "output: lennard.out"
+    read -p "Run lennard.out? [y/n]" option2
+    if [option2 = 'y']
+    then
+        /usr/lib64/mpi/gcc/openmpi/bin/mpirun -n 1 ./lennard.out
+    fi
+else
+    echo "debuging with intel mpiifort"
+    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"
+fi
+
+
+
+

lennard.f90

@@ -47,9 +47,9 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu)
         real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut, fric_term !fric_term = força de ficção
         real(dp) :: x1(2),v1(2),x2(2),p1(2), rs1, rs2, coss, sine, u   
         real(dp), intent(inout), allocatable, dimension(:) :: nRfu
-        integer :: i,j,ct = 0, n1, n2 !,ptr, ptrn
+        integer :: i,j,ct = 0, n1, n2, dox(2), doy(2)
         integer, intent(in) :: mesh(:),domx(2),domy(2), id
-        real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2, dox(2), doy(2)
+        real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2
         type(list_t), pointer :: node, next_node
         type(data_ptr) :: ptr,ptrn
         integer ( kind = 4 ), intent(in) :: ids(8)
@@ -356,9 +356,9 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu)
                                 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] 
-                                    nRfu(n1*6-4:n1*6-2) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                    nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
                                         + nRfu(n1*6-4:n1*6-2)
-                                    nRfu(n2*6-4:n2*6-2) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                    nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
                                         + nRfu(n2*6-4:n2*6-2)
                                     if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
                                         nRfu(n2*6-5) = real(n2,kind(0.d0))
@@ -448,7 +448,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu)
                                     [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
                                 nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
                                     + nRfu(n1*6-4:n1*6-2)
-                                nRfu(n2*6-4:n2*6-2) = [-aux2(2)*r*coss, -aux2(1)*r*sine, u] & 
+                                nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
                                     + nRfu(n2*6-4:n2*6-2)
                                 if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
                                     nRfu(n2*6-5) = real(n2,kind(0.d0))
@@ -679,9 +679,9 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t)
         real(dp), intent(in) :: GField(2), t
         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 !,ptr, ptrn
+        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, dox(2), doy(2)
+        real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2
         type(list_t), pointer :: node, next_node
         type(data_ptr) :: ptr,ptrn
         integer ( kind = 4 ), intent(in) :: ids(8)
@@ -1463,7 +1463,8 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
                         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%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
@@ -1472,15 +1473,17 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
                         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%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%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
@@ -1573,7 +1576,8 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
 
                         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%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

pack.sh

@@ -3,8 +3,8 @@ read -e -p "Enter particle generator file: " VAR1
 
 if [ "$VAR1" = "" ]; then
     echo "No particle generator given."
-    tar -cvf lennard_files.tar csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py
+    tar -cvf lennard_files.tar 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 csv2vtk_particles.py lennard $VAR1 settings.ini verify_settings.py visualizar.py
+    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
 

part_gen_pequeno.py

@@ -24,7 +24,7 @@
 import matplotlib.pyplot as plt
 from random import randint
 
-ff = 1 # formato da região ff = x/y
+ff = 5 # formato da região ff = x/y
 Lx = 400
 arquivo2 = 'poucas1.csv'
 arquivo1 = 'poucas2.csv'

pos_e_vel_inicial_usando_sim_anterior.py

@@ -46,11 +46,15 @@
     try:
         v_files.append(config['par_'+str(i)]['v_file'].split()[0])
         v_files[-1]  = v_files[-1][1:len(v_files[-1])-1]
+        if v_files[-1][0] == '%':
+            print("no velocity file used")
+            v_files[-1] = 'v_file_'+str(i)+'.csv'
     except:
         print("no velocity file used")
         v_files.append('v_file_'+str(i)+'.csv')
         pass 
-
+
+
 step = input('Extract step no. (max {}) '.format(nimpre-1))
 
 positions = []
@@ -77,8 +81,9 @@
 i = 0
 for f in x_files + v_files:
     if os.path.isfile(f) and a == 'n':
-        a = input('Old positions files exist. Backup? y/n ')
-    if a != 'n' or a != 'N': 
+        a = input('Old positions/velocities files exist.\nBackup? {} y/n '.format(f))
+        print(a)
+    if a != 'n': 
         try:
             if i == 0:
                 bkp = 'bkp'

scale_positions.py

@@ -15,19 +15,19 @@
 s = input('Entre os fatores de escala x e y seguidos separados por espaço:\n').split()
 
 data[:,0] = data[:,0]*float(s[0])
-data[:,0] = data[:,0]*float(s[1])
+data[:,1] = data[:,1]*float(s[1])
 
 print('menores posições: [{}, {}] '.format( np.amin(data[:,0]), np.amin(data[:,1]) ))
 print('maiores posições: [{}, {}] '.format( np.amax(data[:,0]), np.amax(data[:,1]) ))
 
 le = len(data)
 for i in range(le-1):
     pos1 = data[i,:]
-    for j in range(i+1,le+1):
+    for j in range(i+1,le):
         pos2 = data[j,:]
         dx = abs(pos1[0]-pos2[0])
         dy = abs(pos1[1]-pos2[1])
-        res = sqrt(dx**2+dy**2)
+        res = np.sqrt(dx**2+dy**2)
         if j == 1 and i == 0:
             mdx = dx
             mdy = dy 
@@ -40,6 +40,6 @@
 print('Menores distâncias entre dois pontos.\ndx = {}, dy = {}, absoluta = {}'.format(mdx,mdy,mres))
 
 
-np.savetxt(file_name[0:len(file_name-4)] + '_sc.csv', data)
+np.savetxt(file_name[0:len(file_name)-4] + '_sc.csv', data)
 
 

settings.ini

@@ -29,7 +29,7 @@
 # se v_file é dada, então ela será usada para dar a velocidade às partículas. Deve ser arquivo .csv. Deixar "%%" se não usar. 
 # se ismolecule = .true. então a partícula é levada em conta no termostato 
 [par_0]
-    x = 'gas.csv' 
+    x = 'poucas1.csv' 
     v = 0 0 # velocidade global
     v_file = '%%v_file_0.csv' # velocidade de cada partícula
     m = 1