Test collision

This procedure primarily discusses the computation of the collision term within the RBG-Maxwell framework.

For detailed information regarding the specific structure and content of the program, please refer to Collision term

  • Firstly, import the required packages for the program:
import warnings
warnings.filterwarnings("ignore")
import math
# specify the system
from RBG_Maxwell.Collision_database.select_system import which_system

plasma_system = 'QGP_system'
which_system(plasma_system)

from RBG_Maxwell.Collision_database.QGP_system.collision_type import collision_type_for_all_species

from RBG_Maxwell.Unit_conversion.main import determine_coefficient_for_unit_conversion, unit_conversion
import numpy as np

# import the main class Plasma
from RBG_Maxwell.Plasma.main import Plasma

hbar, c, lambdax, epsilon0 = 1., 1., 1.6*10**28, 1.
  • Next, input the model parameters.
##### We work in a 6-dimensional space with grid sizes
# the region corresponds to the centrality of 20%~30%, impact parameter b = 8fm
x_left_bound, y_left_bound, z_left_bound = -6./0.197,-6/0.197,-6/0.197 # GeV^-1
px_left_bound, py_left_bound, pz_left_bound = -2,-2,-2 # GeV
x_grid_size, y_grid_size, z_grid_size = 1, 1, 1 # better be odd numbers
px_grid_size, py_grid_size, pz_grid_size = 35,35,35 # better be odd numbers

# The infinitesimal difference of x and p in the current example
# dx, dy, dz are in unit of fm, dpx, dpy, dpz are in GeV
dx, dy, dz = 2*abs(x_left_bound)/x_grid_size,2*abs(y_left_bound)/y_grid_size,2*abs(z_left_bound)/z_grid_size
dpx, dpy, dpz = 2*abs(px_left_bound)/px_grid_size,2*abs(py_left_bound)/py_grid_size,2*abs(pz_left_bound)/pz_grid_size

# dt is set the same for both BEsolver and EMsolver
dt = 0.00025/0.197 # GeV^-1
dt_upper_limit = float(0.00025/0.197/1000)
dt_lower_limit = float(0.00025/0.197*1000)
# The regions for EM calculations are seperated as the source region and the observation region
# In the current example, these regions aree the same as the BEsolver
x_left_boundary_o, y_left_boundary_o, z_left_boundary_o = x_left_bound, y_left_bound, z_left_bound
x_left_boundary_s, y_left_boundary_s, z_left_boundary_s = x_left_bound, y_left_bound, z_left_bound
x_grid_size_o, y_grid_size_o, z_grid_size_o = x_grid_size, y_grid_size, z_grid_size
x_grid_size_s, y_grid_size_s, z_grid_size_s = x_grid_size, y_grid_size, z_grid_size
dx_o, dy_o, dz_o = dx, dy, dz  
dx_s, dy_s, dz_s = dx, dy, dz

# particle types, 0: classical, 1: fermi, 2: boson
# we have two types of particles
particle_type = [0,0]
# masses
masses = np.array([0.3,0.3],dtype=np.float64) # GeV

# charges correspond to u,d,s,u_bar,d_bar,s_bar and gluon
unit_charge = math.sqrt(4*math.pi/137) # ratinalized Lorentz-Heaviside-QCD
charges = np.array([unit_charge*2/3,-unit_charge/3],dtype=np.float64)
degeneracy = np.array([1,1])

# number of total time steps
n_step = 10001

# momentum grids
npx, npy, npz = px_grid_size, py_grid_size, pz_grid_size

# half_px, half_py, half_pz
# momentum range for x and z direction are not import in this case
half_px, half_py, half_pz = np.array([-px_left_bound]*2), \
np.array([-py_left_bound]*2), np.array([-pz_left_bound]*2)
par_list=[m1**2*c**2, m2**2*c**2, (2*math.pi*hbar)**3, hbar**2*c, d_sigma/(hbar**2)]
# load the collision matrix
flavor, collision_type, particle_order = collision_type_for_all_species()
expected_collision_type = ['2TO2']
  • Then, you can configure the parallelization of the program and specify the number of GPUs to be used.
# number of spatial grids
# must be integers and lists
# odd numbers are recomended
# the maximum spatial gird is limited by CUDA, it's about nx*ny*nz~30000 for each card
nx_o, ny_o, nz_o = [1], [1], [1]

# value of the left boundary
# this is the 
x_left_bound_o, y_left_bound_o, z_left_bound_o = [-6./0.197],[-6/0.197],[-6/0.197]

# number samples gives the number of sample points in MC integration
num_samples = 100

# Only specify one spatial region
number_regions = 1

# each spatial should use the full GPU, this number can be fractional if many regions are chosen
# and only one GPU is available
num_gpus_for_each_region = 1.


# since only one region is specified, this will be empty
sub_region_relations = {'indicator': [[]],\
                        'position': [[]]}

# if np.ones are used, the boundaries are absorbing boundaries
# if np.zeros are used, it is reflection boundary
# numbers in between is also allowed
boundary_configuration = {}
for i_reg in range(number_regions):
    bound_x = np.ones([ny_o[i_reg], nz_o[i_reg]])
    bound_y = np.ones([nz_o[i_reg], nx_o[i_reg]])
    bound_z = np.ones([nx_o[i_reg], ny_o[i_reg]])
    boundary_configuration[i_reg] = (bound_x, bound_y, bound_z)
  • Define the initial distribution of particles, primarily specifying the distribution in position space and momentum space.

    The function gives a step distribution function with the formula f(0,p) = f0*\theta(1-p/Qs), usually Qs ~ 1 GeV, for |p|<Qs, f(0,p) = f0, otherwise f(0,p) = 0. Only specified species will be given the value.

def step_function_distribution(px_grid_size = 20,
                               py_grid_size = 20,
                               pz_grid_size = 20,
                               x_grid_size = 1,
                               y_grid_size = 1,
                               z_grid_size = 1,
                               half_x = 3,
                               half_y = 3,
                               half_z = 3,
                               half_px = 2,
                               half_py = 2,
                               half_pz = 2,
                               fa0 = 0.2,
                               fb0 = 0.2,
                               Qs = 1.):
    
    # force spatial grid size to be 1
    if x_grid_size != 1 or y_grid_size != 1 or z_grid_size != 1:
        raise AssertionError("make sure that x_grid_size = 1 and y_grid_size == 1 and z_grid_size == 1")
    
    shape_of_distribution_data = [2, x_grid_size, y_grid_size, z_grid_size,
                                     px_grid_size,py_grid_size,pz_grid_size]
      
    initial_distribution = np.zeros(shape_of_distribution_data)
    
    # momentum grid size
    dpx, dpy, dpz = 2*half_px/px_grid_size, 2*half_py/py_grid_size, 2*half_pz/pz_grid_size # GeV
    
    for ipx in range(px_grid_size):
        for ipy in range(py_grid_size):
            for ipz in range(pz_grid_size):
                
                # central value for each grid 
                px = (ipx+0.5)*dpx - half_px
                py = (ipy+0.5)*dpy - half_py
                pz = (ipz+0.5)*dpz - half_pz
                p = math.sqrt(px**2+py**2+pz**2)
                if p<Qs:
                    
                    # loop through particle species
                    for p_type in range(2):
                        if p_type == 0:
                            initial_distribution[p_type,0,0,0,ipx,ipy,ipz] = fa0
                        else:
                            initial_distribution[p_type,0,0,0,ipx,ipy,ipz] = fb0
                
    return initial_distribution
  • Sets the distribution function for particle 0,1.
num_momentum_levels = 1
num_particle_species = 2

# iniital distribution function
f = {}
for i_reg in range(number_regions):
    f[i_reg] = np.zeros([num_momentum_levels, num_particle_species,\
                         nx_o[i_reg], ny_o[i_reg], nz_o[i_reg], npx, npy, npz])
    
initial_distribution = step_function_distribution(px_grid_size,
                                                   py_grid_size,
                                                   pz_grid_size,
                                                   x_grid_size,
                                                   y_grid_size,
                                                   z_grid_size,
                                                   abs(x_left_bound), 
                                                   abs(y_left_bound), 
                                                   abs(z_left_bound),
                                                   abs(px_left_bound), 
                                                   abs(py_left_bound), 
                                                   abs(pz_left_bound),
                                                   fa0 = 0.1,
                                                   fb0 = 0.4,
                                                   Qs = 1.)

f[0][0][0] = initial_distribution[0]/(2*math.pi)**3
f[0][0][1] = initial_distribution[1]/(2*math.pi)**3
# reshape the distribution function in different regions
for i_reg in range(number_regions):
    f[i_reg] = f[i_reg].reshape([num_momentum_levels, num_particle_species,\
                                 nx_o[i_reg]*ny_o[i_reg]*nz_o[i_reg]*npx*npy*npz])
  • Initiate the RBG-Maxwell framework.
BEx, BEy, BEz, BBx, BBy, BBz = [0],[0],[0],[0],[0],[0]

plasma = Plasma(f,par_list, dt, dt_lower_limit, dt_upper_limit,\
                nx_o, ny_o, nz_o, dx, dy, dz, boundary_configuration, \
                x_left_bound_o, y_left_bound_o, z_left_bound_o, \
                int(npx[0]), int(npy[0]), int(npz[0]), half_px, half_py, half_pz,\
                masses, charges, sub_region_relations,\
                flavor, collision_type, particle_type,\
                degeneracy, expected_collision_type,\
                num_gpus_for_each_region,\
                hbar, c, lambdax, epsilon0, time_stride_back,\
                num_samples = 100, drift_order = 2,\
                rho_J_method="raw", GPU_ids_for_each_region = ["1"])
  • Initiate the iterative calculation, where n_step represents the number of steps for time iteration, and CT indicate whether or not to compute the collision term (1 represents computation, 0 represents no computation).
import time
start_time = time.time()
for i_time in range(n_step):  
    
    if i_time%2000 == 0:
        print('Updating the {}-th time step'.format(i_time))
        new_f = plasma.acquire_values("Distribution")[0][0] 
        np.save("/data/sunmingyan/collision/RBG-Maxwell"+str(i_time),new_f)
    plasma.proceed_one_step(i_time, n_step, processes = {'VT':0., 'DT':0., 'CT':1.},\
                            BEx = [0.], BEy = [0.], BEz = [0.], \
                            BBx = [0.], BBy = [0.], BBz = [0.])
end_time = time.time()
print((end_time-start_time)/3600,'hours')
  • The distribution of particle 0:

particle_collision1

  • The distribution of particle 1:

particle_collision2

results matching ""

    No results matching ""