Skip to content
Snippets Groups Projects
Commit eff68266 authored by Sasiri Juliana Vargas Urbano's avatar Sasiri Juliana Vargas Urbano
Browse files

Upload New File

parent bfa1acf9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 483 additions and 0 deletions
Numba.ipynb 0 → 100644
+ 483
− 0
View file @ eff68266
%% Cell type:code id: tags:
``` python
import random as rm
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import time
from numba import jit, njit, prange
import multiprocessing as mp
```
%% Cell type:code id: tags:
``` python
#------------------------
# FUNCIONES ADICIONALES
#------------------------
#------------------------------------------
# Distribución de espaciamientos discreta
#------------------------------------------
#La distribución de espaciamiento discreta se encuentra con la función P_i_xt(r_,l_,L_),
#la cual utiliza P_i(r,l,i).
#Devuelve un arreglo por ejemplo: [0,0,1,0,0,0]
#Indicando que la distancia entre rma y lmb es 2
@njit
def P_i(r,l,i):
salida=0.
if (abs(r-l)==i):
salida=1.
return salida
@njit
def P_i_xt(r_,l_,L_):
salida = np.zeros(2*L_+1)
for j in range(0,2*L_+1):
salida[j] = P_i(r_,l_,j)
return salida
```
%% Cell type:code id: tags:
``` python
#------------------------------------
# Separación de concentraciones
#------------------------------------
#Entra un arreglo del tipo: [-1,-1, 0, 1, 1]
# Y devuelve dos arreglos: [ 1, 1, 0, 0, 0] --> Para las partículas A
# [ 0, 0, 0, 1, 1] --> Para las partículas B
@njit
def concentracionA_B(entrada):
salida_A = np.zeros(len(entrada))
salida_B = np.zeros(len(entrada))
for i in range(len(entrada)):
if entrada[i] == -1:
salida_A[i] = 1
elif entrada[i] == 1:
salida_B[i] = 1
return salida_A, salida_B
```
%% Cell type:code id: tags:
``` python
#------------
# DINÁMICA
#------------
@njit
def Dinamica(vec, l, Pr, Pj, Ps):
#Vector que almacena las concentraciones de las partículas A y B
concentracion=np.copy(vec)
#Vector que almacena la distribución de espaciamientos
Pesp=np.zeros(len(concentracion))
# Números aleatorios:
#--------------------
#Se tienen 4 números aleatorios. 3 de ellos representan la posibilidad de tener o no,
#reacción, incerción y salto. El restante es el sitio en el que ocurrirá la dinámica
Escogida = rm.randint(0,2*l) #Num aleatorio que escoge el sitio de red para efectuar la dinpamica
particula = concentracion[Escogida] #Tipo de partícula escogida
p_r = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de reacción
p_j = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de incerción
p_s = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de salto
# Saltos
# -------
#Para los saltos se tiene en cuenta lo siguiente:
#Si 0 < p_s < Ps -> ocurre un salto hacia la derecha
#Si Ps < p_s < 2*Ps -> ocurre un salto hacia la izquierda
#(Siempre y cuando el sitio al que saltará la partícula este vacío)
#Esto se realiza para que podamos tener información tanto de la posibilidad de salto,
#como de la posibilidad de no saltar
#Ejemplo: Si Ps=0.5, entonces:
#Salto a la derecha : 0 < p_s < 0.5 --> 50% de las veces
#Salto a la izquierda: 0.5 < p_s < 1 --> 50% de las veces
#No salta --> 0% de las veces
if p_s < 2*Ps:
#Salto general (Se pueden escoger todos los sitios, excepto los extremos)
if Escogida != 0 and Escogida != 2*l and concentracion[Escogida] != 0:
#Salto derecha
if p_s < Ps and concentracion[Escogida + 1] == 0:
concentracion[Escogida + 1] = particula
concentracion[Escogida] = 0
#Salto izquierda
elif Ps < p_s < 2*Ps and concentracion[Escogida - 1] == 0:
concentracion[Escogida - 1] = particula
concentracion[Escogida] = 0
#Extremo izquierdo (dominio de la partícula A)
elif Escogida == 0:
#Incercion: Si p_j < Pj, entonces se efectúa la incerción
if p_j < Pj and concentracion[Escogida] == 0:
concentracion[Escogida] = -1
#Salto: Se considera que hay salto a la derecha
elif p_s < Ps and concentracion[0] != 0 and concentracion[1] == 0:
concentracion[1] = particula
concentracion[0] = 0
#Extremo derecho (dominio de la partícula B)
elif Escogida == 2*l:
#Incersion: Si p_j < Pj, entonces se efectúa la incerción
if p_j < Pj and concentracion[Escogida] == 0:
concentracion[Escogida] = 1
#Salto: Se considera que hay salto a la izquierda
elif p_s < Ps and concentracion[2*l] != 0 and concentracion[2*l-1] == 0:
concentracion[2*l-1] = particula
concentracion[2*l] = 0
# Reaccion
#----------
#La reacción se efectúa si p_r < Pr y si la distancia entre rma y lmb es igual a 1 (si A y B son vecinos)
try:
#El try y el except es puesto porque hay ocasiones en las que se tiene por ejemplo:
#[-1,-1,-1,0,0,0]
#Entonces en este caso no es posible encontrar un lmb
#rma
where1=np.where(concentracion == -1) #devuelve una tupla
rmav=where1[0] #Arreglo de las ubicaciones de las partículas A
rma=rmav[-1] #Partícula más a la derecha
#lmb
where2=np.where(concentracion == 1)
lmbv=where2[0] #Arreglo de las ubicaciones de las partículas B
lmb=lmbv[0] #Partícula más a la izquierda
if (p_r < Pr) and (abs(rma-lmb) == 1): #Las partículas se aniquilan
concentracion[rma] = 0
concentracion[lmb] = 0
#Se encuentran los nuevos rma y lmb para encontrar la distr. de espaciamientos
#rma
where1=np.where(concentracion == -1)
rmav=where1[0]
rma=rmav[-1]
#lmb
where2=np.where(concentracion == 1)
lmbv=where2[0]
lmb=lmbv[0]
#Función que devuelve un arreglo de la distribución de espaciamientos
Pesp= P_i_xt(rma,lmb,l)
except:
pass
return concentracion, Pesp # Devuelve 2 arreglos
```
%% Cell type:code id: tags:
``` python
# ---------
# CICLOS
# ---------
#Repite la función 'Dinamica' una cantidad M de veces, determinada por 'time'
@njit
def ciclos(c_ini, l, Pr, Pj, Ps, time):
resultado=np.copy(c_ini) # copia de la condición inicial
salida_conc_A=np.zeros(len(resultado)) # arreglo de salida para la concentración de A
salida_conc_B=np.zeros(len(resultado)) # arreglo de salida para la concentración de B
salida_esp=np.zeros(len(resultado)) # arreglo de salida para la dist. de espaciamientos
for i in range(time+1):
resultado, espaciamiento = Dinamica(resultado, l, Pr, Pj, Ps) #Dinamica para el tiempo i
A, B = concentracionA_B(resultado) # Separación de arreglos
return A, B, espaciamiento # Devuelve 3 arreglos
```
%% Cell type:code id: tags:
``` python
# ----------------
# REPETICIONES
# ----------------
#Repite la función 'ciclos' una cantidad determinada (rep) de repeticiones
pool = mp.Pool(mp.cpu_count())
global results
results=[]
def get_result(result):
results.append(result)
#Uso de paralelización
def repeticiones(rep, c_ini, l, Pr, Pj, Ps, time):
for i in prange(rep):
pool.apply_async(ciclos, args=(c_ini, l, Pr, Pj, Ps, time), callback=get_result)
# Cada 1000 iteraciones sacamos un print en pantalla para conocer el estado de la simulacion
if not(i % 100):
print('Working in iteration number {} out of {}.'.format(i, rep))
```
%% Cell type:code id: tags:
``` python
#--------------
# ANALISIS
#--------------
def ANALISIS(L=5, dx=1., dt=1., D=0.5, j=1., K=1., tiempo=2000, rep=10, opcion='graficos', nombre = 'Data'):
#PARÁMETROS
beta= K/dx # Tasa de reacción
lamb = D/dx**2 #Tasa de salto
P_S = lamb * dt # Probabilidad de salto
P_J= j #Probabilidad de incersion
P_R= beta*dt #Probabilidad de reacción
#CONDICION INICIAL
# 0 : Vacío
# -1 : Partícula A
# 1 : Particula B
condicion_inicial = np.zeros(2*L+1)
for i in range(2*L+1): #Ilustracion para L=3: [-1,-1,-1,1,1,1]
if i < L:
condicion_inicial[i] = -1
else:
condicion_inicial[i] = 1
#-----------------------------------OPCIONES------------------------------------
if opcion=='graficos':
x = np.linspace(-L, L, 2*L+1) #Posiciones
z = np.linspace(0, 2*L, 2*L+1) #Distancias dist. de espaciamientos
salida_A=np.zeros(len(x))
salida_B=np.zeros(len(x))
salida_Desp=np.zeros(len(z))
# Se almacena la información de TODA la dinámica en 3 arreglos de salida
# Se suman cada iteración de tiempo para posteriormente hacer su respectiva normalización y promedio
# y encontrar la densidad local contínua y dist. de espaciamientos contínua
#Función que aplica la función ciclos para 'rep' repeticiones
repeticiones(rep, condicion_inicial, L, P_R, P_J, P_S, tiempo)
#ESTO FALTARIA POR OPTIMIZAR!!!!!!!! ( HACE LENTO EL CÓDIGO )
pool.close()
pool.join()
#La variable global results es la que almacena los resultados
for item in results:
A, B, D = item
salida_A += A
salida_B += B
salida_Desp += D
results.clear() #Se limpia la variable global
#Normalización y cálculo de la densidad de concentraciones
salida_A = salida_A/(rep*2*L)
salida_B = salida_B/(rep*2*L)
salida_Desp = salida_Desp/salida_Desp.sum()
#--------------------------------------------------------
#------------------ARCHIVOS------------------------------
#--------------------------------------------------------
#Escribe 3 archivos txt en la carpeta Archi:
# - Info para A
# - Info para B
# - Info para dist. de espaciamientos
print('TERMINA LA SIMULACIÓN')
#-----------------------------------------------
info = 'Archi/'+str(nombre)+'_info.txt'
File_info = open(info, 'w')
File_info.write('L: '+str(L)+'\n')
File_info.write('D: '+str(D)+'\n')
File_info.write('j: '+str(j)+'\n')
File_info.write('K: '+str(K)+'\n')
File_info.write('tiempo: '+str(tiempo)+'\n')
File_info.write('Repetición: '+str(rep)+'\n')
#-----------------------------------------------
nombre1 = 'Archi/'+str(nombre)+'A.txt'
nombre2 = 'Archi/'+str(nombre)+'B.txt'
nombre3 = 'Archi/'+str(nombre)+'Desp.txt'
File_A = open(nombre1, 'w')
File_B = open(nombre2, 'w')
File_P_i = open(nombre3, 'w')
for i in range(len(x)):
File_A.write(str(x[i])+' '+str(salida_A[i])+'\n')
File_B.write(str(x[i])+' '+str(salida_B[i])+'\n')
File_P_i.write(str(z[i])+' '+str(salida_Desp[i])+'\n')
File_A.close()
File_B.close()
File_P_i.close()
#----------------------------------------------------------
#----------------------GRAFICOS----------------------------
#----------------------------------------------------------
#-------------------PERFILES DE DENSIDAD-------------------------
fig1=plt.figure(1)
plt.plot(x,salida_A,c='b')
plt.plot(x,salida_B,c='r')
plt.legend(['Particulas A', 'Particulas B'],loc='upper right',bbox_to_anchor=(1.35, 1))
plt.xlabel('x')
plt.ylabel('Perfil de densidad')
plt.show()
#---------------DISTRIBUCIÓN DE ESPACIAMIENTOS--------------------
fig2=plt.figure(2)
plt.plot(z,salida_Desp,c='b', alpha=0.5)
plt.xlabel('x')
plt.ylabel('Distribución de espaciamientos')
plt.show()
#-------------------------------------------------------------------
#---------------
# Dinamica
#---------------
#Imprime en cada paso de tiempo, el vector concentración
#Con esta opción se puede visualizar lo que ocurre en cada paso de tiempo
elif opcion=='dinamica':
salida_din=np.copy(condicion_inicial)
for i in range(0,tiempo):
print('tiempo ', i)
salida_din, otro = Dinamica(salida_din, L, P_R, P_J, P_S)
print(salida_din)
print(' ')
#--------------------------------------------------------------------
#---------------
# Estacionario
#---------------
#Comprueba si se llega al estado estacionario en el tiempo ingresado. Esto se hace comparando
#el sistema en el tiempo = 'tiempo' y tiempo = 'tiempo' + 1000
elif opcion=='estacionario':
pruebas=[]
#Evolución hasta el tiempo = 'tiempo'
A_est, B_est, Desp_est = ciclos(condicion_inicial, L, P_R, P_J, P_S, tiempo)
estacionario = B_est - A_est
pruebas.append(estacionario)
print('Estado en la iteración '+ str(tiempo)+':')
print(estacionario)
print(' ')
#Evolución hasta el tiempo = 'tiempo + 1000'
A_est2, B_est2, Desp_est2 = ciclos(estacionario, L, P_R, P_J, P_S, 1000)
estacionario2 = B_est2 - A_est2
pruebas.append(estacionario2)
print('Estado en la iteración '+ str(tiempo+1000)+':')
print(estacionario2)
print(' ')
salida_estacionario = pruebas[0] == pruebas[1]
num_verdadero = list(salida_estacionario).count(True)
print('Estado estacionario: ' + str(len(salida_estacionario) == num_verdadero))
return len(salida_estacionario) == num_verdadero
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment