calculo coeficientes em py

diffusion_coef.py

@@ -0,0 +1,173 @@
+import numpy as np 
+import pandas as pd 
+import os
+import configparser
+import matplotlib.pyplot as plt
+from glob import glob
+import time
+import curses # para progressbar
+
+def integral(F,i,j,dt):
+    # sh = np.shape(F)
+    I = np.trapz(F[i,j,:],dx=dt)                
+    # I = np.zeros((sh[0],sh[1]))
+    # xi = [-np.sqrt(3/5), 0, np.sqrt(3/5)]
+    #xp = np.linspace(0, dt*sh[2], sh[2])
+    #w = [5/9, 8/9, 5/9]
+    # for i in range(sh[0]):
+    #     for j in range(sh[1]):
+            # for s in range(sh[2]-1):
+                # I[i,j] += 0.5*(dt)*( \
+                #     w[0]*np.interp( dt*xi[0]/2 + ((s+1)*dt+s*dt)/2 , xp, F[i,j,:]) + \
+                #     w[1]*np.interp( dt*xi[0]/2 + ((s+1)*dt+s*dt)/2  , xp, F[i,j,:]) + \
+                #     w[2]*np.interp( dt*xi[0]/2 + ((s+1)*dt+s*dt)/2  , xp, F[i,j,:]) \
+                #     )
+    return I
+
+def progressbar(step,total,message, stdscr):
+    # usar import curses
+    # inicializar uma strig vazia fora do loop para stdcr
+    if step == 0:  
+        stdscr = curses.initscr() 
+        stdscr.addstr(0, 0, message + "|" + " "*20 + "| {} %".format(0))
+        stdscr.refresh()
+        return stdscr
+    elif step < total-1:
+        a = int(20*step/total)
+        b = 20-a
+        stdscr.addstr(0, 0, message + "|" + "#"*a + " "*b + "| {} %".format(5*a))
+        stdscr.refresh()
+        return stdscr
+        #print("#",end='')
+    else:
+        a = int(20*step/total)
+        b = 20-a        
+        stdscr.addstr(0, 0, message + "|" + "#"*a + " "*b + "| {} %".format(5*a))
+        curses.endwin()
+        print("\n")    
+        return "0"
+
+
+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()
+config.read(res_dir + 'settings.txt')
+N = int(config['global']['N'].split()[0])
+dimx = float(config['global']['dimX'].split()[0])
+dimy = float(config['global']['dimY'].split()[0])
+n_files = int(config['out_files']['out_files'].split()[0])
+
+ntype = int(config['global']['Ntype'].split()[0])
+t_fim = float(config['global']['t_fim'].split()[0])
+dt = float(config['global']['dt'].split()[0])
+F = np.array([0,0])
+quant = []
+sigma = []
+epsilon = []
+rs = [] # raio sólido
+mass = []
+tipo = [0]*N
+
+dt = t_fim/n_files
+
+a = input("Enter the subdomains mesh dimensions.\n")
+a = a.split()
+mesh = np.array([int(a[0]), int(a[1])])
+a = input("Enter a location. xmin xmax ymin ymax\n")
+if a == '':
+    region = [0, mesh[0], 0, mesh[1]]
+else:
+    region = [int(x) for x in a.split()]
+
+Vol = (dimx/mesh[0])*(dimy/mesh[1]) # volume dos elementos da malha
+for i in range(ntype):
+    quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
+    rs.append(float(config['par_'+str(i)]['rs'].split()[0]))
+    sigma.append(float(config['par_'+str(i)]['sigma'].split()[0]))
+    epsilon.append(float(config['par_'+str(i)]['epsilon'].split()[0]))
+    mass.append(float(config['par_'+str(i)]['m'].split()[0]))
+j,k = 0,0
+for i in range(len(quant)):
+    for j in range(quant[i]):
+        tipo[j+k] = i
+
+    k = sum(quant[0:i+1])
+tipo = pd.DataFrame(tipo, columns=["tipo"]) # numero id da partícula
+
+hx = dimx/mesh[0]
+hy = dimy/mesh[1]
+
+nsteps = n_files # int(input('Enter the number of steps:\n'))
+
+density_map = np.zeros((mesh[0],mesh[1], nsteps+1))
+r = np.zeros((mesh[0],mesh[1],nsteps+1))
+
+mqd = np.zeros(nsteps)
+
+
+step = 0 
+n1,n2 = 0,0
+stdscr = 's' #para barra de progresso
+sample_list = []
+r0 = []
+r1 = []
+r = np
+
+while step <= nsteps: 
+
+    particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])]
+    #stdscr = progressbar(step,nsteps,'Computing:',stdscr)
+    print(step)
+    pos = pd.read_csv(res_dir+"position.csv."+str(step), header=None, names = ["x","y"])
+    vel = pd.read_csv(res_dir+"velocity.csv."+str(step), header=None, names = ["v_x", "v_y"])
+    n = [x for x in range(len(pos))]
+    n = pd.DataFrame(n, columns=["n"]) # numero id da partícula
+    pos_vel = pd.concat([n,pos,vel,tipo],axis=1)
+
+    for nn in range(len(pos_vel)):
+        xp = int(pos_vel.loc[nn,'x']//hx)
+        yp = int(pos_vel.loc[nn,'y']//hy)
+        if xp == mesh[0]:
+            xp = xp - 1
+        if yp == mesh[1]:
+            yp = yp - 1
+        particle_map[xp][yp].append( pos_vel.loc[nn,'n'] )
+        density_map[xp,yp,step] += 1
+
+    if (step == 0):
+        for i in range(region[0],region[1]):
+            for j in range(region[2],region[3]):
+                for nn in range(len(particle_map[i][j])):
+                    sample_list.append(particle_map[i][j][nn])
+                    n1 = particle_map[i][j][nn]
+                    r0.append([pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']])
+        r0 = np.array(r0)
+        r1 = np.zeros( (len(sample_list),2) )
+        desv = 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)
+
+            # mqd[step-1] += ((r1[0] - r0[nn,0])**2 + (r1[1] - r0[nn,1])**2)
+    step += 1
+
+
+mqd = np.sum(desv**2,axis=0)/len(sample_list)
+
+plt.plot(mqd)
+plt.xlabel('time steps x ....')
+plt.ylabel('4D')
+
+# plt.show()
+
+
+

lennard.f90

@@ -1890,7 +1890,7 @@ program main
     type(CFG_t) :: my_cfg
     character(5) :: particle
     character(4) :: wall
-    character(11) :: nome,arquivox, arquivov = '%' !nome de partícula deve ter até 11 caracteres
+    character(20) :: nome,arquivox, arquivov = '%' !nome de partícula deve ter até 20 caracteres
     type(string) :: part_nomes(10) ! vetor de strings
     character(1) :: optio
     type(data_ptr) :: ptr !integer,pointer :: ptr,ptrn

part_gen_pequeno.py

@@ -0,0 +1,126 @@
+# #!/usr/bin/env python3
+# # -*- coding: utf-8 -*-
+# """
+# Created on Sun Jan 27 22:47:50 2019
+
+# Este programa é para colocar as partículas em seus lugares iniciais
+
+# @author: gabriel
+# """
+
+# import numpy as np
+# import matplotlib.pyplot as plt
+# from random import randint
+
+# arquivo1 = 'pequeno.csv'
+
+# N = 400
+# L1 = 200
+# L2 = 200
+# spac = 0.999*((L1*L2)/(N))**0.5
+# d1 = int(L1/spac)
+# d2 = int(L2/spac)
+# N = d1*d2
+# N1 = N*2
+# # Posição das partículas
+# p1 = np.zeros([N1,2])
+
+# #posição de referência
+# refe = np.array([0.01,0.01])
+# print("d1*spac = {}, L1 = {}, d2*spac = {}, L2 = {}\n".format(d1*spac+refe[0],L1, d2*spac+refe[1], L2))
+# print("N = {}, spac = {}\n".format(N*2,spac))
+# cont = 0
+
+
+# for i in range(d1):
+#     for j in range(d2):
+#         p1[cont,:] = [i*spac+refe[0],j*spac+refe[1]]
+#         cont = cont+1
+# refe = np.array([1199.9,599.9])     
+
+# with open(arquivo1,'w') as file:
+#     for i in range(N1):
+#         file.write('{},{}\n'.format(p1[i,0],p1[i,1]))
+
+# plt.scatter(p1[:,0],p1[:,1], label='pequeno')
+# plt.show()
+
+
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Jan 27 22:47:50 2019
+
+Este programa é para colocar as partículas em seus lugares iniciais
+
+@author: gabriel
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from random import randint
+
+Ntot = 0 
+arquivo1 = 'bin_peq1.csv'
+
+N = 200
+L1 = 200
+L2 = 200
+spac = 0.999*((L1*L2)/(N))**0.5
+d1 = int(L1/spac)
+d2 = int(L2/spac)
+N = d1*d2
+Ntot = N + Ntot
+# Posição das partículas
+p1 = np.zeros([N,2])
+
+#posição de referência
+refe = np.array([0.01,0.01])
+print("d1*spac = {}, L1 = {}, d2*spac = {}, L2 = {}\n".format(d1*spac+refe[0],L1, d2*spac+refe[1], L2))
+print("N1 = {}, spac = {}\n".format(N,spac))
+cont = 0
+
+
+for i in range(d1):
+    for j in range(d2):
+        p1[cont,:] = [i*spac+refe[0],j*spac+refe[1]]
+        cont = cont+1   
+
+with open(arquivo1,'w') as file:
+    for i in range(N):
+        file.write('{},{}\n'.format(p1[i,0],p1[i,1]))
+
+plt.scatter(p1[:,0],p1[:,1], label='pequeno1')
+
+arquivo1 = 'bin_peq2.csv'
+
+N = 200
+L1 = 200
+L2 = 200
+spac = 0.999*((L1*L2)/(N))**0.5
+d1 = int(L1/spac)
+d2 = int(L2/spac)
+N = d1*d2
+Ntot = N + Ntot
+# Posição das partículas
+p1 = np.zeros([N,2])
+
+#posição de referência
+refe = np.array([0.01+spac/2,0.01+spac/2])
+print("d1*spac = {}, L1 = {}, d2*spac = {}, L2 = {}\n".format(d1*spac+refe[0],L1, d2*spac+refe[1], L2))
+print("N2 = {}, spac = {}\n".format(N,spac))
+cont = 0
+
+for i in range(d1):
+    for j in range(d2):
+        p1[cont,:] = [i*spac+refe[0],j*spac+refe[1]]
+        cont = cont+1  
+
+with open(arquivo1,'w') as file:
+    for i in range(N):
+        file.write('{},{}\n'.format(p1[i,0],p1[i,1]))
+
+print('Ntot = {}'.format(Ntot))
+plt.scatter(p1[:,0],p1[:,1], label='pequeno2')
+
+plt.show()