diff --git a/csv2vtk_particles.py b/csv2vtk_particles.py index 9f14cb3..2236131 100755 --- a/csv2vtk_particles.py +++ b/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,14 +113,21 @@ 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 @@ -119,7 +135,12 @@ 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)) \ No newline at end of file + #shutil.move(via+'/cell.csv.'+str(fnum),via+'/'+folder+'/cell.csv.'+str(fnum)) diff --git a/debug_lennard.sh b/debug_lennard.sh index 15c7839..7d51201 100755 --- a/debug_lennard.sh +++ b/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 + + + + diff --git a/lennard.f90 b/lennard.f90 index b542215..4e72507 100755 --- a/lennard.f90 +++ b/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,7 +1473,8 @@ 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 @@ -1480,7 +1482,8 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, 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 diff --git a/pack.sh b/pack.sh index 8c96458..02386b5 100755 --- a/pack.sh +++ b/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 diff --git a/part_gen_pequeno.py b/part_gen_pequeno.py index 1301391..0912bc8 100644 --- a/part_gen_pequeno.py +++ b/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' diff --git a/pos_e_vel_inicial_usando_sim_anterior.py b/pos_e_vel_inicial_usando_sim_anterior.py index 11fa672..51c523e 100644 --- a/pos_e_vel_inicial_usando_sim_anterior.py +++ b/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' diff --git a/scale_positions.py b/scale_positions.py index 1707aad..17994bc 100644 --- a/scale_positions.py +++ b/scale_positions.py @@ -15,7 +15,7 @@ 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]) )) @@ -23,11 +23,11 @@ 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) diff --git a/settings.ini b/settings.ini index 9e4acf4..ee5ad7b 100755 --- a/settings.ini +++ b/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