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_steprepresents the number of steps for time iteration, andCTindicate 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:
- The distribution of particle 1:
