How to use RBG-Maxwell

RBG_Maxwell

The main purpose of this part is to let users quickly master the use of RBG-Maxwell framework.

RBG-Maxwell RBG-Maxwell CUDA: v10.0 (shields.io) Ray: v1.0 (shields.io) Plasma: simulation (shields.io)

1、Usage via an example

The following codes demonstrate an example of how to use RBG-Maxwell.

1.1、 Set the initial conditions

  • 1、First, we need to invoke the following package:
import warnings
warnings.filterwarnings("ignore")

# specify the system
from RBG_Maxwell.Collision_database.select_system import which_system

plasma_system = 'Fusion_system'
which_system(plasma_system)

from RBG_Maxwell.Collision_database.Fusion_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
from RBG_Maxwell.Plasma.main import Plasma
  • 2、Then, we need to specify the unit conversion factor:
    • Determine the conversion factors for the International System of Units (IS) and the Flexible System of Units (FS) by configuring the spatial grid, velocity, and charge parameters.
# here we use a pure electron system

# give the quantities in SI
# the spatial grid is chosen to be dx=dy=10**(-5) m
dx = dy = 10**(-5)
dz = 1.

# velocity is roughly 5*10**(6) m/s
v = 5*10**6

# charge
Q = 1.6*10**(-19) 

# momentum is roughly 10**(-30)kg*10**7m/s
momentum = 10**(-23)
dp = (10**(-25)*10**(-23)*10**(-23))**(1/3)

# the total number of particles are 5*10**(-13)/(1.6*10**(-19))
# put these particles in 71 spatial grids in z direction
# in 201 spatial grids in y direction
# and 100 momentum grids
# in each phase grid we have dn = 21.89755448111554
# the average value of distribution is roughly 
dn = 0.2189755448111554
f = dn/(dp**3*dx*dy*dz)
df = f

# time scale
dt = 10**(-13)

# Now find the coefficient
hbar, c, lambdax, epsilon0 = determine_coefficient_for_unit_conversion(df, dx, dp, Q, f, v, dt)

conversion_table = \
unit_conversion('SI_to_LHQCD', coef_J_to_E=lambdax, hbar=hbar, c=c, k=1., epsilon0=epsilon0)
conversion_table_reverse = \
unit_conversion('LHQCD_to_SI', coef_J_to_E=lambdax, hbar=hbar, c=c, k=1., epsilon0=epsilon0)

For a detailed procedure of unit conversion you can refer to Conversion.

  • 3、Next, we need to initialize the plasma system:
    • The primary dataset comprises the parameters for time-space discretization, grid quantities, particle classification, and collision classification.
    • In the present demonstration, a two-dimensional system consisting of pure electron plasma is utilized, and the initial spatial arrangement of the system is depicted in below.

Plasma_initialization

# time step, and spatial infinitesimals
# dt is 10**(-13) s, dx = dy = dz = 10**(-5) m
dt, dx, dy, dz = 10**(-13)*conversion_table['second'], \
                 10**(-5)*conversion_table['meter'], \
                 10**(-5)*conversion_table['meter'], \
                 10**(-5)*conversion_table['meter']
dt_upper_limit = float(10**(-1)*conversion_table['second'])
dt_lower_limit = float(10**(-9)*conversion_table['second'])

# we have only one type of particle e-
num_particle_species = 1

# treat the electron as classical particles
particle_type = np.array([0])

# masses, charges and degenericies are
masses, charges, degeneracy = np.array([9.11*10**(-31)*conversion_table['kilogram']]), \
                              np.array([-1.6*10**(-19)*conversion_table['Coulomb']]),\
                              np.array([1.])

# momentum grids
npx, npy, npz = 1, 201, 1

# 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([9.11*10**(-31)*5*10**6*conversion_table['momentum']]), \
                            np.array([9.11*10**(-31)*5*10**6*conversion_table['momentum']]),\
                            np.array([9.11*10**(-31)*5*10**6*conversion_table['momentum']])

par_list=[m1**2*c**2, m2**2*c**2, (2*math.pi*hbar)**3, hbar**2*c, d_sigma/(hbar**2)]


dpx, dpy, dpz = 2*half_px/npx, 2*half_py/npy, 2*half_pz/npz

# load the collision matrix
flavor, collision_type, particle_order = collision_type_for_all_species()
expected_collision_type = ['2TO2']

The parameters related to collisions can be found in Collision_database. The program describes in detail how to set up different colliding plasmas. We also set up the quark-gluon plasma system and the fusion system in this program.

  • 4、Set parallel calculation parameters for the plasma system:
    • Including the number of Monte Carlo particles, the number of regions, the number of GPUs in the regions, etc.
# 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], [251], [111]

# value of the left boundary
# this is the 
x_left_bound_o, y_left_bound_o, z_left_bound_o = [-0.5*dx],\
                                                 [-125.5*dy],\
                                                 [-55.5*dz]

# 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 = 0.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)

For details, please refer to Plasma.

  • 5、Set the distribution function and boundary conditions of the plasma system
# 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])


# The initial velocity of the electrons is 1.87683*10**6 m/s, corresponds to the momentum value
# 9.11*10**(-31)*1.87683*10**6*conversion_table['momentum'] ~ 408.770512.
# The following code specifies the momentum grid index
dpy = 2*half_py/npy
a = 9.11*10**(-31)*1.87683*10**6*conversion_table['momentum']
ipy = [i for i in range(npy) if (-half_py+dpy*(i-0.5))<=a<=(-half_py+dpy*(i+1))][0]

'''
The total number of particles is 5*10**(-13)/(1.6*10**(-19)) ~ 31249999.999999996.
Put these particles in 101 grids, the number density of the particles is 
31249999.999999996/(101*dx*dy*dz) ~ 756837.0957070973 -> Delta_N/Delta_V.
The particles only possess the momentum region of size dpx*dpy*dpz,
hence the distribution function at each phase space grid is
756837.0957070973/(dpx*dpy*dpz) ~ 756837.0957070973/(2*half_px/npx*2*half_py/npy*2*half_pz/npz)
~ 1234.63197049
'''
dn_dv = 5*10**(-14)/(1.6*10**(-19))/(101*dx*dy*dz*dpx*dpy*dpz)
f[0][0, 0, 0,9,5:106,0,ipy,0] = dn_dv


# 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])

'''
We add an external magnetic field of 10 T in the +y direction
'''
BBy = [10*conversion_table['Tesla']*np.ones(nx_o[0]*ny_o[0]*nz_o[0])]
BEx, BEy, BEz, BBx, BBz = [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"])

1.2、System evolution and results output

  • Set the time step and perform the plasma system evolution.
n_step = 10001
number_rho = []
EM = []
charged_rho = []
dis = []
VT= []
DT = []
import time
start_time = time.time()
for i_time in range(n_step):  
    
    # if i_time%1000 == 0:
    #     dis.append(plasma.acquire_values("Distribution"))            
    plasma.proceed_one_step(i_time, n_step, processes = {'VT':1., 'DT':1., 'CT':0.},\
                            BEx = BEx, BEy = BEy, BEz = BEz, BBx = BBx, BBy = BBy, BBz = BBz)
    if i_time%1000 == 0:     
        print('Updating the {}-th time step'.format(i_time))
        number_rho.append(plasma.acquire_values("number_rho/J"))
        charged_rho.append(plasma.acquire_values("Electric rho/J"))
    EM.append(plasma.acquire_values('EM fields on current region'))
end_time = time.time()
  • Using pictures to show the evolution of the system.
# spatial distribution
import matplotlib.pyplot as plt
xi, yi = np.mgrid[1:252:1,1:112:1]
fig, axes = plt.subplots(ncols=5, nrows=2, figsize = (15,5))
for jj in range(2):
    for kk in range(5):
        axes[jj,kk].pcolormesh(xi, yi, number_rho[(jj*5+kk+1)][0][0].reshape([nx_o[0],ny_o[0],nz_o[0]])[0])

result

Users can easily build own plasma system by learning this example.

For more details in the framework, please refer to

Plasma

For more examples, please refer to

Example

results matching ""

    No results matching ""