Help on Monte Carlo Simulation of temperature relaxation of He-Ar mixture

[jorvz]

Hey, so i`ve been working on this code for a while for my research but it just doesn't seem to work well, i think it is misscalculating in some step, probaly on the relative velocities, because it keeps exploding the value for some really high values, but i am not shure. Can someone help me finding the source of the problem here?

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import random as rd
from matplotlib import animation
from matplotlib.animation import PillowWriter

# Constants
kB = 1.38e-23  # Boltzmamn constant
amu = 1.66e-27  # Atomic mass unity
m_He = 4 * amu  # He mass
m_Ar = 38.9 * amu  # Ar mass

# Initial conditions
n = 100000  # Total number os particles
L = 10.  # Length of the box
T_He = 300 # initial temperature of He
T_Ar = 300
#T_mistura = 500 #Temperature of the mixture
conc_He = 0.90  # molar fraction of He in the mixture
#T_Ar = (T_mistura - conc_He*T_He)/(1.-conc_He) #initial temperature of Ar
T_mistura = conc_He*T_He + (1. - conc_He)*T_Ar

n_He = int(n * conc_He)  # Number of He particles
n_Ar = n - n_He  # Number of Ar particles
sigma_He = np.sqrt((kB * T_He) / m_He)  # standard deviation for He speeds
sigma_Ar = np.sqrt((kB * T_Ar) / m_Ar)  #standard deviation for Ar speeds

T = 100 #Total time for the simulations
dt = 0.00002 # time interval

va = np.sqrt(2.*kB*T_mistura/(m_He + m_Ar)) #termal speed of the mixture
def generate():
    # Generating the initial speed and position of the particles
    x = L * np.random.rand(n, 3)  # random positions
    v = np.zeros((n, 3),dtype='float64')  # velocities array
    v[:n_He,:] = np.random.normal(0, sigma_He, size=(n_He, 3))  # He velocities

    v[n_He:,:] = np.random.normal(0, sigma_Ar, size=(n_Ar, 3))  # Ar velocities
    
    return v/va, x #va to let it dimentionless

def plot(v):
    fig, ax = plt.subplots(1,2, figsize=(25,10))
    
    a = ax[0]
    v_mag_He = np.linalg.norm(v[:n_He,:], axis=1)  # Magnitude das velocidades do hélio
    a.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He')

    # Plot the histogram for the initial velocities of Argon
    v_mag_Ar = np.linalg.norm(v[n_He:,:], axis=1)  # Magnitude das velocidades do argônio
    a.hist(v_mag_Ar, bins=50, density=True, alpha=0.5, label='Ar')

    # Configuring the graph
    a.set_title('Distribuição de velocidades')
    a.set_xlabel('Velocidade (m/s)')
    a.set_ylabel('Densidade de probabilidade')
    a.legend(fontsize=20)
    
    a.grid()
    
    a = ax[1]
    v_mag_He = v[:n_He,0]  # Magnitude das velocidades do hélio
    a.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He')

    # Plot the histogram for the initial velocities of Argon
    v_mag_Ar = v[n_He:,0]  # Magnitude das velocidades do argônio
    a.hist(v_mag_Ar, bins=50, density=True, alpha=0.5, label='Ar')

    # Configuring the graph
    a.set_title('Distribuição de velocidades')
    a.set_xlabel('Velocidade (m/s)')
    a.set_ylabel('Densidade de probabilidade')
    a.legend(fontsize=20)
    
    a.grid()
    
    
    plt.show()
#Importing the matrix for viscosity
muHA = pd.read_csv('tables/muHA.dat', sep="\s+",index_col=False).set_index(['0'])

#Importing the matrix for cosx and the cross sections
xiHe4He4 = pd.read_csv('tables/xiHe4He4.dat', sep="\s+", skiprows=4, header=None)
xiHe4Ar = pd.read_csv('tables/xiHe4Ar.dat', sep="\s+", skiprows=4, header=None)
xiArAr = pd.read_csv('tables/xiArAr.dat', sep="\s+", skiprows=4, header=None)

def mover_particulas(v, x):
    # Move the particles in every Δt of the simulation
    #x = x + v*dt
    x,v = bordas(x,v)
    return x

def temp_cinetica(vel):
    
    u = (1./n)*np.sum(v,axis=0)
    
    T_He = (m_He*va**2)/(3. * n_He * kB)* np.sum(v[:n_He]**2)
    T_Ar = (m_Ar*va**2)/(3. * n_Ar * kB)* np.sum(v[n_He:]**2)
    T = conc_He*T_He + (1.-conc_He)*T_Ar
    return T_He, T_Ar, T


def bordas(pos, vel):
    """
    Mantém as partículas dentro de uma caixa de tamanho L,
    fazendo com que elas refletam nas paredes da caixa quando 
    colidem com elas.
    
    Args:
        pos (numpy array): Matriz Nx3 com as posições das N partículas.
        vel (numpy array): Matriz Nx3 com as velocidades das N partículas.
        L (float): Tamanho da caixa quadrada.
        
    Returns:
        numpy array: Matriz Nx3 com as novas posições das N partículas.
        numpy array: Matriz Nx3 com as novas velocidades das N partículas.
    """
    # Verifica se alguma partícula está fora da caixa.
    new_pos = pos.copy()
    new_vel = vel.copy()

    # Verifica se as partículas estão dentro da caixa e corrige as posições e velocidades caso contrário
    for i in range(pos.shape[0]):
        for j in range(pos.shape[1]):
            if new_pos[i,j] < 0:
                new_pos[i,j] = abs(new_pos[i,j])
                new_vel[i,j] = -new_vel[i,j]
            elif new_pos[i,j] > L:
                new_pos[i,j] = 2*L - new_pos[i,j]
                new_vel[i,j] = -new_vel[i,j]

    return new_pos, new_vel


def colisao(v, tipo,sg):
    # type 1 = He-He
    # type 2 = Ar-Ar
    # type 3 = He-Ar
    p1 = 0
    p2 = 0
    
    while p1==p2:
        if tipo==1:
            #Selecting particles in the interval [0, n_He]
            p1 = int(np.floor((n_He)*rd.random()))
            p2 = int(np.floor((n_He)*rd.random()))
            
            #Particles masses
            m1 = m_He
            m2 = m_He
            
            #G constant for the type of collision
            G = 400.
            
            #Cross section matrix
            σ = xiHe4He4[100].to_numpy()
            
            #matrix for the values of cosX
            ξ = xiHe4He4.to_numpy() 
        
        elif tipo == 2:
            #Selecting particles in the interval (n_He, n_Ar]
            p1 = int(n_He + np.floor(rd.random()*n_Ar))
            p2 = int(n_He + np.floor(rd.random()*n_Ar))
            
            #Particles masses
            m1 = m_Ar
            m2 = m_Ar
            
            #G constant for the type of collision
            G = 100.
            
            #Cross section matrix
            σ = xiArAr[100].to_numpy()
            ξ = xiArAr.to_numpy() 
        
        else:
            #Selecting particles in the interval [0, n_He] and (n_He, n_Ar]
            p1 = int(np.floor((n_He)*rd.random()))
            p2 = int(n_He + np.floor(rd.random()*n_Ar))
            
            #Particle masses
            m1 = m_He
            m2 = m_Ar
            
            #G constant for the type of collision
            G = 300.
            
            #Cross section matrix
            σ = xiHe4Ar[100].to_numpy()
            
            #matrix for the values of cosX
            ξ = xiHe4Ar.to_numpy() 
    
    
        
    v1 = v[p1,:] #speed of the particle 1
    v2 = v[p2,:] #speed of the particle 2
    
    #Calculating relative speed
    g = v1 - v2
    gr = np.linalg.norm(g) #norm of g
    
    #Post-collisional relative speed
    gl = np.zeros(3,dtype='float64')
    
    
    #Center of mass velocity
    G_cm = (m1*v1 + m2*v2)/(m1 + m2)
    
    #Calcuting the index of cross section array
    k = int(np.floor((np.log(1. + gr*va/G)/np.log(1.005)) + 0.5))
   
    if k>899:k=899
        
    vr = np.sqrt(g[1]*g[1] + g[2]*g[2])  

        
    if ((σ[k]*gr)/(sg) > rd.random()): #Accept-rejection condition
        
        print(σ[k], gr)  
        n = int(np.floor(rd.random()*100))
        
        #Monte Carlo parameters fot the post collisional velocities
        e =  2. * rd.random() * np.pi
        cosx = ξ[k][n]
        sinx = np.sqrt(1. - cosx*cosx)
        
        if vr>1e-15:
        #Components of the post collisional relative speed
            gl[0] = g[0]*cosx + vr*sinx*np.sin(e)
            gl[1] = g[1]*cosx + (gr*g[2]*np.cos(e) - g[0]*g[1]*np.sin(e))*sinx / vr
            gl[2] = g[2]*cosx - (gr*g[1]*np.cos(e) + g[0]*g[2]*np.sin(e))*sinx / vr
        else:
            gl[0] = g[0]*cosx
            gl[1] = g[0]*sinx*np.cos(e)
            gl[2]= g[0]*sinx*np.sin(e)
               
        
        #calculating the post collision velocities
        v1 = G_cm + 0.5*gl
        v2 = G_cm - 0.5*gl
        
        #setting the new velocities
        v[p1,:] = v1
        v[p2,:] = v2
        return v, σ[k]*gr
    else:
        return v, sg
    
    
t=0.
μ = muHA[f"{conc_He}"][np.floor(T_mistura)]
σg11 = 15.
σg12 = 15.
σg22 = 15.
v, x = generate()
v_t = np.zeros((n, 3, T))
Temp = np.zeros((T,3))
for t in range(0,T,1):
    
    #plot(v)
    #x = mover_particulas(v,x)
    
    
    np.append(Temp[t,:],temp_cinetica(v)) 
    N_col_HH = int((1./n)*n_He*(n_He - 1.)*μ*dt*σg11)
    N_col_AA = int((1./n)*n_Ar*(n_Ar - 1.)*μ*dt*σg22)
    N_col_HA = int((2./n)*n_He*n_Ar*μ*dt*σg12)
    
    print(f'N11 = {N_col_HH}, N12 = {N_col_HA}, N22 ={N_col_AA}')
    print(f'σ11 = {σg11}, σ12 = {σg12}, σ22 ={σg22}')

    for i in range(0,N_col_HH):
        v, sg11 = colisao(v, 1,σg11)
        if (sg11 > σg11): σg11 = sg11
    for i in range(0,N_col_AA):
        v, sg22 = colisao(v, 2,σg22)
        if (sg22 > σg22): σg22 = sg22
    for i in range(0,N_col_HA):
        v, sg12 = colisao(v, 3,σg12)
        if (sg12 > σg12): σg12 = sg12
    
    print(t)
plt.plot(np.arange(0,T,1), Temp[:,0],label='He')
plt.plot(np.arange(0,T,1), Temp[:,1],label='Ar')
plt.plot(np.arange(0,T,1), Temp[:,2], label='T')
plt.legend()
fig, ax = plt.subplots(1,1, figsize=(10,10))

def animate(i):
    ax.clear()
    v_mag_He = np.linalg.norm(v_t[:n_He,:,i], axis=1)  # Magnitude das velocidades do hélio
    ax.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He')

    # Plot the histogram for the initial velocities of Argon
    v_mag_Ar = np.linalg.norm(v_t[n_He:,:,i], axis=1)  # Magnitude das velocidades do argônio
    ax.hist(v_mag_Ar, bins=50, density=True, alpha=0.5,label='Ar')

    # Configuring the graph
    #ax.title('Distribuição de velocidades')
    ax.set_xlabel('Velocidade (m/s)',fontsize=20)
    ax.set_ylabel('Densidade de probabilidade',fontsize=20)
    
    #ax.set_xlim(0,10)
    #ax.set_ylim(0,2)
    ax.grid()
    
    
ani = animation.FuncAnimation(fig, animate,frames=T,interval=10)
ani.save('ani.gif',writer='pillow')

[farshid]

there are some steps that you can take to debug the code and narrow down the potential causes of the issue:

Check the inputs and outputs of the functions to make sure that they are consistent and make sense.

Print out the intermediate results of the calculations to see where things might be going wrong.

Use a debugger to step through the code and see where the program might be diverging from what you expect.

Use unit tests to validate the output of the code against known inputs and outputs.

Based on the code, it looks like you are simulating the motion of particles in a box, and calculating their velocities and positions. Some potential sources of the problem could be the calculation of relative velocities or the handling of the boundary conditions.

One thing you can do is to print out the values of the velocities and positions at each step of the simulation to see if there are any extreme or unexpected values. You can also plot the distributions of the velocities to see if they look reasonable.

Another thing you can try is to use smaller time steps and see if the problem persists. Sometimes, if the time step is too large, the simulation can become unstable and produce unexpected results.