diff --git "a/anota\303\247\303\265es.txt" "b/anota\303\247\303\265es.txt" index b70ed97..b04086f 100644 --- "a/anota\303\247\303\265es.txt" +++ "b/anota\303\247\303\265es.txt" @@ -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) diff --git a/lennard.f90 b/lennard.f90 index e39bf5d..3c10725 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -2,7 +2,7 @@ module fisica contains !Calcula velocidades distribuidas de acordo com MaxwellBolzmann ! factor = sqrt() = sqrt(kb*T/m) por componente - subroutine MaxwellBoltzmann(malha,mesh,factor) + subroutine MaxwellBoltzmann(malha,mesh,factor,propriedade) use linkedlist use mod1 use randnormal @@ -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 @@ -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 @@ -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 @@ -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, "||")') & @@ -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, & diff --git a/part_gen_pequeno.py b/part_gen_pequeno.py index b0198e7..f1c3ae8 100644 --- a/part_gen_pequeno.py +++ b/part_gen_pequeno.py @@ -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 diff --git a/potential_analysis.py b/potential_analysis.py index a21d543..99c2bec 100644 --- a/potential_analysis.py +++ b/potential_analysis.py @@ -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') @@ -91,7 +91,7 @@ def einstein_condutividade(delta_e,na): quant = [] sigma = [] epsilon = [] -rs = [] # raio sólido +rs = [] # raio solido mass = [] tipo = [0]*N @@ -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))] @@ -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 diff --git a/settings.ini b/settings.ini index 476996b..3b534f7 100755 --- a/settings.ini +++ b/settings.ini @@ -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 @@ -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] diff --git a/teste.f90 b/teste.f90 index f520313..238e1be 100644 --- a/teste.f90 +++ b/teste.f90 @@ -1,5 +1,8 @@ + + program teste - use mpi + + ! use mpi ! call system('mkdir temp') ! open(10,file='temp/test.txt',status="replace") @@ -7,28 +10,40 @@ program teste ! print*,'a' implicit none - real :: sen(2), res(8) - integer ( kind = 4 ) status(MPI_STATUS_SIZE) - integer ( kind = 4 ) ierr - integer ( kind = 4 ) np - integer ( kind = 4 ) tag - integer ( kind = 4 ) id + integer :: a, b + character(len=10) :: time + character(len=10) :: date + + call date_and_time(date,time) + write(*,'("#:Execution date: ", a,":",a,":",a,2x, a,"/",a,"/",a)') & + time(1:2),time(3:4),time(5:8),date(5:6),date(7:8),date(1:4) + print*, date + a = 20 + b = 21 + + + ! real :: sen(2), res(8) + ! integer ( kind = 4 ) status(MPI_STATUS_SIZE) + ! integer ( kind = 4 ) ierr + ! integer ( kind = 4 ) np + ! integer ( kind = 4 ) tag + ! integer ( kind = 4 ) id - call MPI_Init ( ierr ) - call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr ) - call MPI_Comm_size ( MPI_COMM_WORLD, np, ierr ) + ! call MPI_Init ( ierr ) + ! call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr ) + ! call MPI_Comm_size ( MPI_COMM_WORLD, np, ierr ) - res = [0,0,0,0,0,0,0,0] - sen = [1.0,1.0] - sen = sen*id +1 + ! res = (0,0,0,0,0,0,0,0) + ! sen = (1.0,1.0) + ! sen = sen*id +1 - call MPI_ALLGATHER(sen, 2, MPI_REAL, res, 2, MPI_REAL , MPI_COMM_WORLD, ierr) + ! call MPI_ALLGATHER(sen, 2, MPI_REAL, res, 2, MPI_REAL , MPI_COMM_WORLD, ierr) - print '("id ", f3.1, " res ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1)', & - id, res(1),res(2),res(3), res(4) ,res(5) ,res(6), res(7), res(8) + ! print '("id ", f3.1, " res ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1," ", f3.1)', & + ! id, res(1),res(2),res(3), res(4) ,res(5) ,res(6), res(7), res(8) - print*, (2>1 .and. ( 1==2 .or. 3==2 )) - call MPI_Finalize ( ierr ) + ! print*, (2>1 .and. ( 1==2 .or. 3==2 )) + ! call MPI_Finalize ( ierr ) end program teste