Skip to content

Commit

Permalink
correção no lingrad T
Browse files Browse the repository at this point in the history
  • Loading branch information
@g7fernandes
g7fernandes committed Sep 2, 2019
1 parent 3791322 commit c83d1db
Show file tree
Hide file tree
Showing 6 changed files with 135 additions and 109 deletions.
27 changes: 10 additions & 17 deletions anotações.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,20 +49,13 @@ Ideias:
Variar densidade (número)
Variar Temperatura

In [1]: m_C = 12

In [2]: m_Ar = 40

In [3]: fe = (.74*(.335+.127)) # ou 0.64*

In [4]: rs = 10

In [5]: rC = 0.07

In [6]: sigma_Ar = 0.341

In [7]: nC = (sigma_Ar*rs**2*fe/(rC**2) # numero de atomos na região

In [8]: m_Cnd = nC*m_C/m_Ar # massa adimensional do carbono

In [9]: print(nC)
m_C = 12
m_Ar = 40
fe = (.74*(.335+.127)) # ou 0.64*
rs = 5
rC = 0.07
sigma_Ar = 0.341
nC = (sigma_Ar*rs)**2*fe/(rC**2) # numero de atomos na região
m_Cparticle = nC*m_C/m_Ar # massa adimensional do carbono
print(nC)
print(m_Cparticle)
79 changes: 45 additions & 34 deletions lennard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ module fisica
contains
!Calcula velocidades distribuidas de acordo com MaxwellBolzmann
! factor = sqrt(<vd²>) = sqrt(kb*T/m) por componente
subroutine MaxwellBoltzmann(malha,mesh,factor)
subroutine MaxwellBoltzmann(malha,mesh,factor,propriedade)
use linkedlist
use mod1
use randnormal
Expand All @@ -15,15 +15,22 @@ subroutine MaxwellBoltzmann(malha,mesh,factor)
integer, intent(in) :: mesh(2)
type(list_t), pointer :: node
type(data_ptr) :: ptr
type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade

do i = 2,mesh(2)+1
do j = 2,mesh(1)+1
node => list_next(malha(i,j)%list)
do while (associated(node))
ptr = transfer(list_get(node), ptr)
ptr%p%v(1) = ptr%p%v(1) + factor(1)*GaussDeviate()
ptr%p%v(2) = ptr%p%v(2) + factor(2)*GaussDeviate()
ptr%p%F = [0,0]
if (propriedade(ptr%p%grupo)%x_lockdelay == 0) then
ptr%p%v(1) = ptr%p%v(1) + factor(1)*GaussDeviate()
ptr%p%v(2) = ptr%p%v(2) + factor(2)*GaussDeviate()
ptr%p%F = [0,0]
else
ptr%p%v(1) = 0
ptr%p%v(2) = 0
ptr%p%F = [0,0]
end if
node => list_next(node)
end do
end do
Expand Down Expand Up @@ -4278,8 +4285,10 @@ program main

if (termostato_vel_scaling .and. force_lingradT > 1) then
temp_Td(2) = t_fim
print*, "O termostato estará ligado ao longo de toda simulação"
print*, "temp_Td = ", temp_Td
if (id == 0) then
print*, "O termostato estará ligado ao longo de toda simulação"
print*, "temp_Td = ", temp_Td
end if
end if


Expand Down Expand Up @@ -4446,10 +4455,10 @@ program main
deallocate(grupo)

! adiciona às particulas velocidade inicial de acordo com distribuição maxwell boltzmann
if (vd(1) /= 0 .and. vd(2) /= 0) call MaxwellBoltzmann(malha,mesh, vd)
if (vd(1) /= 0 .and. vd(2) /= 0) call MaxwellBoltzmann(malha,mesh, vd,propriedade)
if (vd(1) == 0 .and. vd(2) == 0 .and. Td > 0) then
vd = sqrt([Td, Td])
call MaxwellBoltzmann(malha,mesh,vd)
call MaxwellBoltzmann(malha,mesh,vd,propriedade)
end if
! adicionamos uma celula para corrigir o fato de estarmos usando fantasmas
cold_cells = cold_cells + 1
Expand Down Expand Up @@ -4547,7 +4556,7 @@ program main
! pois eles não vão imprimir x e v pra csv

! id = 0 é o processo raiz

call MPI_barrier(MPI_COMM_WORLD, ierr)
if (id == 0) then
call date_and_time(date,time)
write(*,'("|| Program started at: ", a,":",a,":",a,2x, a,"/",a,"/",a, "||")') &
Expand Down Expand Up @@ -4818,33 +4827,35 @@ program main

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
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

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 = sum(domx_aux)*(Td_hot-Td_cold)/(2*mesh(2))
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.

j2 = j2 + 1
end do
end if
end if

if (i == interv(j)) then
deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, &
Expand Down
12 changes: 6 additions & 6 deletions part_gen_pequeno.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,14 +24,14 @@
import matplotlib.pyplot as plt
from random import randint

ff = 1 # formato da região ff = x/y
Lx = 400
arquivo1 = 'poucas1.csv'
arquivo2 = 'poucas2.csv'
ff = 8 # formato da região ff = x/y
Lx = 2000
arquivo2 = 'particulas8x1.csv'
arquivo1 = 'moleculas8x1.csv'

ratio = 0.8
N1 = 500
N2 = 100
N1 = 5000
N2 = 130

L1 = Lx
L2 = Lx/ff
Expand Down
35 changes: 21 additions & 14 deletions potential_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ def einstein_condutividade(delta_e,na):
a = int(input("Enter the number of the folder\n"))
res_dir = dirlist[a]

zip_rfup = ZipFile(res_dir+'/rFuP.zip','r')
#zip_rfup = ZipFile(res_dir+'/rFuP.zip','r')
zip_positions = ZipFile(res_dir+'/positions.zip','r')
zip_velocities = ZipFile(res_dir+'/velocities.zip','r')

Expand All @@ -91,7 +91,7 @@ def einstein_condutividade(delta_e,na):
quant = []
sigma = []
epsilon = []
rs = [] # raio sólido
rs = [] # raio solido
mass = []
tipo = [0]*N

Expand Down Expand Up @@ -124,18 +124,26 @@ def einstein_condutividade(delta_e,na):
hx = dimx/mesh[0]
hy = dimy/mesh[1]

nsteps = n_files # int(input('Enter the number of steps:\n'))
print("Planned output {} files (steps).\n".format(n_files))
nsteps = int(input('Enter the number of steps (final number):\n'))

density_map = np.zeros((mesh[0],mesh[1], nsteps+1))
Kmap = np.zeros((mesh[0],mesh[1],nsteps+1))
Vmap = np.zeros((mesh[0],mesh[1],nsteps+1))
step = 1

print("There are {} files.".format(n_files-1))
try:
step = int(input("Enter the initial step [0]: "))
except:
step = 0

density_map = np.zeros((mesh[0],mesh[1], nsteps+1 - step))
Kmap = np.zeros((mesh[0],mesh[1],nsteps+1 - step))
Vmap = np.zeros((mesh[0],mesh[1],nsteps+1 - step))
a = nsteps - step
stepini = step
if pbar:
bar = progressbar.ProgressBar(max_value=nsteps)
bar = progressbar.ProgressBar(max_value=(a))
while step <= nsteps:
particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])]
rfup = pd.read_csv(zip_rfup.open('rF_u_P.csv.'+str(step)), header=None, names = ["RxFy","RyFx","u","px","py"])
# rfup = pd.read_csv(zip_rfup.open('rF_u_P.csv.'+str(step)), header=None, names = ["RxFy","RyFx","u","px","py"])
pos = pd.read_csv(zip_positions.open("position.csv."+str(step)), header=None, names = ["x","y"])
vel = pd.read_csv(zip_velocities.open("velocity.csv."+str(step)), header=None, names = ["v_x", "v_y"])
n = [x for x in range(len(pos))]
Expand All @@ -151,24 +159,23 @@ def einstein_condutividade(delta_e,na):
if yp == mesh[1]:
yp = yp - 1
particle_map[xp][yp].append( pos_vel.loc[nn,'n'] )
density_map[xp,yp,step] += 1
density_map[xp,yp,step-a] += 1
# "print C"
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])):
n1 = particle_map[i][j][nn]
m = mass[pos_vel.loc[n1,'tipo']]
Kmap[i,j,step] += (pos_vel.loc[n1,'v_x']**2+pos_vel.loc[n1,'v_y']**2)*m/(2)
Kmap[i,j,step-a] += (pos_vel.loc[n1,'v_x']**2+pos_vel.loc[n1,'v_y']**2)*m/(2)
# Vmap[i,j,step] += pos_vel.loc[n1,'u']
if pbar:
bar.update(step)
bar.update(step-stepini)
step += 1

KE = Kmap/density_map
a = KE.shape[2]

KE = np.nan_to_num(KE)
KE = np.sum(KE,axis=2)/a
print("nsteps {} = a {} ?".format(nsteps,a))
plt.figure()
if mesh[1] == 1:
T = KE #np.sum(KE,axis=0) # energia cinética ao longo de X
Expand Down
40 changes: 20 additions & 20 deletions settings.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,27 +4,27 @@
# Argonio: epsilon/kb = 120K, sigma = 0.341nm
[global]
nimpre = 100 # quntidade de saidas
N = 5100 # número de partículas
N = 4100 # número de partículas
Ntype = 2 # número de tipos de partícula
t_fim = 10000 # tempo de execução
t_fim = 35000 # 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 = 20000 # Dimensões da região de cálculo
dimY = 5000 #
mesh = 200 40 # 300 60 # malha(elementos)
dimX = 10000 # Dimensões da região de cálculo
dimY = 1250 #
mesh = 200 30 # 300 60 # malha(elementos)
rcut = 3 # *sigma
wall = 'ppee' # 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 = .true. # liga termostato por velocity scaling
termostato_nose_hoover = .true. #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)
temp_Td = 0 35000 250 # Tempo de velocity scaling [inicial, final, iterações_por_scaling]
cold_cells = 1 10 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
hot_cells = 190 200 1 40 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
temp_Td = 0 10000 500 # Tempo de velocity scaling [inicial, final, iterações_por_scaling]
cold_cells = 1 12 1 30 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
hot_cells = 190 200 1 30 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
Td_hot = 2 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
Td_cold = 1 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato
force_lingradT = 20 # (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 = 0 0 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td')
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')
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
Expand All @@ -33,29 +33,29 @@
# 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 = 'moleculas_sc.csv'
x = 'moleculas8x1_sc.csv'
v = 0 0 # velocidade global
v_file = 'v_file_0.csv' # velocidade de cada partícula
v_file = '%%v_file_0.csv' # velocidade de cada partícula
m = 1
nome = 'g1'
sigma = 1 # Lennard Jones
epsilon = 1 # Lennard Jones
quantidade = 5000
quantidade = 4000
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 = 'particulas_sc.csv'
x = 'particulas8x1_sc.csv'
v = 0 0 # velocidade global
v_file = '%%v_file_0.csv' # velocidade de cada partícula
m = 100
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
x_lockdelay = 1000 # só vai poder mudar de posição a partir de t = x_lockdelay
rs = 10 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
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.
#[par_2]
Expand Down
Loading

0 comments on commit c83d1db

Please sign in to comment.