com py que permite continuar a partir de simulação

diffusion_coef.py

@@ -110,6 +110,11 @@ def progressbar(step,total,message, stdscr):
 
 mqd = np.zeros(nsteps)
 
+step = input('Enter de initial step:\n')
+try: 
+    step = int(step)
+except:
+    step = 0 
 
 step = 0 
 n1,n2 = 0,0
@@ -150,24 +155,31 @@ def progressbar(step,total,message, stdscr):
         r0 = np.array(r0)
         r1 = np.zeros( (len(sample_list),2) )
         desv = np.zeros((len(particle_map[i][j]),nsteps))
+        desvx = np.zeros((len(particle_map[i][j]),nsteps))
+        desvy = np.zeros((len(particle_map[i][j]),nsteps))
     else:
         for nn in range(len(sample_list)):
             n1 = sample_list[nn]
             r1[nn,:] =  [pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']]
         r1 = np.array(r1) - r0
-        desv[:,step-1] = np.sqrt(r1[:,0]**2 + r1[:,1]**2)
-
+        desvx[:,step-1] = r1[:,0] 
+        desvy[:,step-1] = r1[:,1] 
             # mqd[step-1] += ((r1[0] - r0[nn,0])**2 + (r1[1] - r0[nn,1])**2)
     step += 1
 
-
+desv = np.sqrt(desvx**2 + desvy**2)
 mqd = np.sum(desv**2,axis=0)/len(sample_list)
 
+plt.figure()
 plt.plot(mqd)
 plt.xlabel('time steps x ....')
 plt.ylabel('4D')
 
-# plt.show()
+plt.figure()
+for i in range(10):
+    plt.plot(desvx[np.random.randint(len(desv[:,0])),:])
+
+plt.show()
 
 
 

lennard.f90

@@ -1876,7 +1876,7 @@ program main
     implicit none
 !    Variáveis
     integer :: N,Ntype,i=1, nimpre,j = 1, ii, quant = 0,mesh(2), cont = 1, aux1 = 0,cont2 = 1,domx(2), domy(2)
-    integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4)
+    integer :: subx, suby, NMPT, j2, cold_cells(4), hot_cells(4), nimpre_init
     integer, target :: k
     real(dp), dimension(:,:), allocatable :: v, x, celula !força n e n+1, velocidades e posições
     real(dp), dimension(:), allocatable :: icell,jcell, nxv, nxv_send !dimensões das celulas e vetor de resultado pra imprimir
@@ -1926,7 +1926,9 @@ program main
     call CFG_add(my_cfg, "global%dt",0.0_dp,&
       "time step")
     call CFG_add(my_cfg, "global%t_fim",0.0_dp,&
-      "simulation time")
+      "simulation time") 
+    call CFG_add(my_cfg, "global%nimpre_init",0,&
+      "start nimpre") 
     call CFG_add(my_cfg, "global%nimpre",1,&
       "number of result files")
     call CFG_add(my_cfg, "global%dimX",1.0_dp,&
@@ -1966,6 +1968,7 @@ program main
     call CFG_get(my_cfg,"global%N",N)
     call CFG_get(my_cfg,"global%Ntype",Ntype)
     call CFG_get(my_cfg,"global%t_fim",t_fim)
+    call CFG_get(my_cfg,"global%nimpre_init",nimpre_init)
     call CFG_get(my_cfg,"global%dt",dt)
     call CFG_get(my_cfg,"global%nimpre",nimpre)
     call CFG_get(my_cfg,"global%dimX", dimX)
@@ -2268,6 +2271,12 @@ program main
     j = 1
     j2 = 1
 
+    if (nimpre_init > 0) then
+        j = nimpre_init
+        t = t_fim*nimpre_init/nimpre
+        i = interv(j)
+    end if
+
     do while (t_fim > t)
         ! print*, "L 1990"
         call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t)  !altera Força

part_gen_pequeno.py

@@ -63,9 +63,9 @@
 Ntot = 0 
 arquivo1 = 'bin_peq1.csv'
 
-N = 200
-L1 = 200
-L2 = 200
+N = 1100
+L1 = 400
+L2 = 400
 spac = 0.999*((L1*L2)/(N))**0.5
 d1 = int(L1/spac)
 d2 = int(L2/spac)
@@ -94,9 +94,9 @@
 
 arquivo1 = 'bin_peq2.csv'
 
-N = 200
-L1 = 200
-L2 = 200
+N = 1100
+L1 = 400
+L2 = 400
 spac = 0.999*((L1*L2)/(N))**0.5
 d1 = int(L1/spac)
 d2 = int(L2/spac)

pos_e_vel_inicial_usando_sim_anterior.py

@@ -0,0 +1,95 @@
+# -*- coding: utf-8 -*-
+"""
+Spyder Editor
+
+Este é um arquivo de script temporário.
+"""
+from glob import glob
+import os, shutil
+import configparser
+
+dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__))
+dirlist = glob(dirname + "/*/")
+print("Choose a folder there the results are contained:\nNo | Folder")
+for a in range(len(dirlist)):
+    print("{} | {}\n".format(a,dirlist[a]))
+a = int(input("Enter the number of the folder\n"))
+res_dir = dirlist[a]
+
+config = configparser.ConfigParser()
+list_files = os.listdir(res_dir)
+
+if os.path.isfile(res_dir + 'settings.txt'):
+    config.read(res_dir + 'settings.txt')
+    nimpre =  int(config['global']['nimpre'].split()[0])
+else:
+    print('The still not converted to vtk. Probable incomplete simulation.\n'+\
+          'settings.ini will be used.\n')
+    config.read(dirname + '/settings.ini') 
+    nimpre = len(list_files)/2    
+
+
+N = int(config['global']['N'].split()[0])
+ntype = int(config['global']['Ntype'].split()[0])
+quant = []
+x_files = []
+v_files =[]
+
+for i in range(ntype):
+    quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
+    x_files.append(config['par_'+str(i)]['x'].split()[0])
+    x_files[-1]  = x_files[-1][1:len(x_files[-1])-1] 
+    v_files.append(config['par_'+str(i)]['v_file'].split()[0])
+    v_files[-1]  = v_files[-1][1:len(v_files[-1])-1]
+
+step = input('Extract step no. (max {}) '.format(nimpre-1))
+
+positions = []
+with open(res_dir + 'position.csv.'+step) as file:
+    for line in file:
+        positions.append(line)
+velocities = []
+with open(res_dir + 'velocity.csv.'+step) as file:
+    for line in file:
+        velocities.append(line)
+
+# verificar se não é desejavel fazer backup de arquivos de posição e velocidade
+a == 'n'
+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':
+        try:
+            if i == 0:
+                bkp = 'bkp'
+                os.mkdir(bkp)
+                i = -1
+        except:
+            while i >= 0:
+                try: 
+                    bkp = 'bkp'+str(i)
+                    os.mkdir(bkp)
+                    i = -1
+                except:
+                    i += 1
+        try:                    
+            shutil.move(dirname + '/' + f, dirname  + '/' + bkp + '/' + f)                
+        except:
+            print('file ' +f+ ' not found!')
+
+k = 0            
+for i in range(ntype):
+    with open(x_files[i],'w') as file:
+        for j in range(quant[i]):
+            file.write(positions[k])
+            k += 1
+
+k = 0        
+for i in range(ntype):
+    with open(v_files[i],'w') as file:
+        for j in range(quant[i]):
+            file.write(velocities[k])
+            k += 1        
+
+

settings.ini

@@ -1,60 +1,61 @@
 # PARAMETROS GLOBAIS
 [global]
-    nimpre  = 150 # quntidade de saidas 
-    N = 1940 # número de partículas 
-    Ntype = 3 # número de tipos de partícula
-    t_fim = 1000  # tempo de execução 
+    nimpre  = 200 # quntidade de saidas 
+    N = 1089 # número de partículas 
+    Ntype = 1 # número de tipos de partícula
+    t_fim = 50  # tempo de execução 
+    nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). 
     dt = 0.0001 # passo de tempo
-    dimX = 1200 # Dimensões da região de cálculo
-    dimY = 600 #
-    mesh = 96 48 # malha(elementos)
-    rcut = 2.5 # *sigma
+    dimX = 410 # Dimensões da região de cálculo
+    dimY = 410 #
+    mesh = 82 82 # malha(elementos)
+    rcut = 5 # *sigma
     wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste
     Td = -1 # 2.897189213660259e+24 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E
     temp_Td = 0 1000 100 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] 
     cold_cells = 1 5 1 52 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
     hot_cells = 145 150 1 52 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
-    Td_hot = 100e+23  # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
-    Td_cold = 10e+23 # temperatura nas hot cells constante. Td_cold = -1 desliga termostato
-    vd = 1 1 # velocidade distribuida com Maxwell Boltzmann
-    NMPT = 500 # Número máximo de partículas que poderá mudar de process. Se não souber estimar, usar -1 ou o valor de N.
+    Td_hot = -1 #100e+23  # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
+    Td_cold = -1 #10e+23 # temperatura nas hot cells constante. Td_cold = -1 desliga termostato
+    vd = 2 2 # velocidade distribuida com Maxwell Boltzmann
+    NMPT = 100 # 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
 # PARTICLE PROPERTIES
 # 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. 
 [par_0]
-    x = 'grupo0.csv' 
+    x = 'bin_peq1.csv' 
     v = 0 0 # velocidade global
     v_file = '%vgrupo0.csv' # velocidade de cada partícula
     m = 1
     nome = 'g1'
-    sigma = 1 # Lennard Jones
+    sigma = 2 # Lennard Jones
     epsilon = 1 # Lennard Jones
-    quantidade = 1900
+    quantidade = 1089
     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				
-[par_1]
-    x = 'p_p.csv'
-    v = 0 0
-    v_file = "%vp_p.csv"
-    m = 20
-    nome = 'g3'
-    sigma = 1 # Lennard Jones
-    epsilon = 1 # Lennard Jones
-    quantidade = 20
-    x_lockdelay = 300 # só vai poder mudar de posição a partir de t = x_lockdelay
-    rs = 5
-    fric_term = 0 # fricção artifical			
-[par_2]
-    x = 'p_g.csv'
-    v = 0 0
-    v_file = "%vp_g.csv"
-    m = 20
-    nome = 'g3'
-    sigma = 1 # Lennard Jones
-    epsilon = 1 # Lennard Jones
-    quantidade = 20
-    x_lockdelay = 300 # só vai poder mudar de posição a partir de t = x_lockdelay
-    rs = 5
-    fric_term = 0 # fricção artifical		
+#[par_1]
+#    x = 'p_p.csv'
+#    v = 0 0
+#    v_file = "%vp_p.csv"
+#    m = 20
+#    nome = 'g3'
+#    sigma = 1 # Lennard Jones
+#    epsilon = 1 # Lennard Jones
+#    quantidade = 20
+#    x_lockdelay = 300 # só vai poder mudar de posição a partir de t = x_lockdelay
+#    rs = 5
+#    fric_term = 0 # fricção artifical			
+#[par_2]
+#    x = 'p_g.csv'
+#    v = 0 0
+#    v_file = "%vp_g.csv"
+#    m = 20
+#    nome = 'g3'
+#    sigma = 1 # Lennard Jones
+#    epsilon = 1 # Lennard Jones
+#    quantidade = 20
+#    x_lockdelay = 300 # só vai poder mudar de posição a partir de t = x_lockdelay
+#    rs = 5
+#    fric_term = 0 # fricção artifical		

viscosity.py

@@ -264,8 +264,8 @@ def progressbar(step,total,message, stdscr):
             etakp[i,j,step] = Vol*Ik[i,j,step]*Ip[i,j,step]/density_map[i,j,step]
             etapp[i,j,step] = Vol*Ip[i,j,step]**2/density_map[i,j,step]
 
-# plt.plot(etakk[0,0,:]+etakp[0,0,:]+etapp[0,0,:])
-# plt.show()
+plt.plot(etakk[0,0,:]+etakp[0,0,:]+etapp[0,0,:])
+plt.show()
 
 
 

viscosity_allen.py

@@ -117,7 +117,9 @@ def progressbar(step,total,message, stdscr):
 sample_list = []
 r0 = []
 p0 = []
-Q = 
+Qxy = np.zeros((nsteps,1))
+Qyx = np.zeros((nsteps,1))
+#eta = np.zeros((nsteps,1))
 
 while step <= nsteps: 
 
@@ -150,17 +152,25 @@ def progressbar(step,total,message, stdscr):
                     m = mass[pos_vel.loc[n1,'tipo']]
                     p0.append([pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m])
         r0 = np.array(r0)
+        r1 = np.zeros( (len(sample_list),2) )
+        p0 = np.array(p0)
+        p1 = np.zeros( (len(sample_list),2) )
     else:
         for nn in range(len(sample_list)):
             n1 = sample_list[nn]
-            r1.append([pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']])
+            r1[nn,:] =  [pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']]
             m = mass[pos_vel.loc[n1,'tipo']]
-            p1.append([pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m])
+            p1[nn,:] = [pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m]
+
+        Qxy[step-1] = np.sum(r1[:,0]*p1[:,1], axis=0)/Vol
+        Qyx[step-1] = np.sum(r1[:,1]*p1[:,0], axis=0)/Vol
     step += 1
 
-    mqd = mqd/len(sample_list)
+eta = (Qxy - np.sum(r0[:,0]*p0[:,1]))/len(sample_list) #completar 
+
 
-# plt.show()
+plt.plot(eta)
+plt.show()