Test point expansion of pure electron plasma

This example is primarily designed to simulate the drifting process of a circle pure electron plasma in the YOZ plane.

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

  • Firstly, import the required packages for the program:
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

# import the main class Plasma
from RBG_Maxwell.Plasma.main import Plasma
  • Next, input the model parameters to create a unit conversion table.

    Here, the parameters are given in units of the International System of Units (SI).

# here we use a pure electron system
# the relavant PIC code is given by Jian-Nan Chen

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

# velocity
v_max = 5*1.87683*10**6

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

# maximum momentum
momentum = 10**(-30)*v_max

# the momentum grid is set to be 
# half_px=half_pz=half_py=momentum
npy=npz=201
npx=1
dpz = dpy = 2*momentum/npy
dpx = 2*momentum
dp_volume = dpx*dpy*dpz
dp = (dpy+dpz)/2

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

# number of maximum particles in each phase grid
n_max = 5*10**(-15)/(1.6*10**(-19))/(npx*npy*npz)

# number of averaged particles in each spatial grid
nx = 1
ny = nz = 101
n_average = 5*10**(-15)/(1.6*10**(-19))/(nx*ny*nz)


E = 10**5
B = 10**(-5)

# Now find the coefficient
hbar, c, lambdax, epsilon0 = determine_coefficient_for_unit_conversion(dt, dx, dx*dy*dz, dp, dp_volume,\
                                                                       n_max, n_average, v_max, E, B)
  • Perform conversions between the International System of Units (SI) and Natural units based on the given parameters.

    The specific content of the program can be referred to in Unit_Conversion. In addition, we have provided an example to help users better understand the unit conversion program, Test-Conversion.

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)
  • Perform unit conversion on the input parameters and specify the type of ion collision.
# time step, and spatial infinitesimals
# dt is 10**(-13) s, dx = dy = dz = 10**(-5) m
dt, dx, dy, dz = dt*conversion_table['second'], \
                 dx*conversion_table['meter'], \
                 dy*conversion_table['meter'], \
                 dz*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 = npx, npy, npz

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

dpx, dpy, dpz = 2*half_px/npx, 2*half_py/npy, 2*half_pz/npz
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 = [nx], [ny], [nz]

# value of the left boundary
# this is the 
x_left_bound_o, y_left_bound_o, z_left_bound_o = [-nx/2*dx],\
                                                 [-ny/2*dy],\
                                                 [-nz/2*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 = 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.

  • Momentum distribution of particles
import math
v_collect = np.zeros([npx,npy,npz])
for ipx in range(npx):
    px = dpx[0]*(ipx+0.5) - half_px[0]
    for ipy in range(npy):
        py = dpy[0]*(ipy+0.5) - half_py[0]
        for ipz in range(npz):
            pz = dpz[0]*(ipz+0.5) - half_pz[0]
            # current velocity
            pnorm = math.sqrt((px**2+py**2+pz**2))
            Enorm = math.sqrt(px**2+py**2+pz**2+masses[0]**2*c**2)
            v = pnorm/Enorm
            v_collect[ipx, ipy, ipz] = v
v_average = masses[0]*\
            1.87683*10**6*conversion_table['meter']/conversion_table['second']*\
            1/math.sqrt(1-(1.87683*10**6/c)**2)
sigma = 0.2*v_average
p_v = 1/math.sqrt(2*math.pi*sigma**2)*np.exp(-(np.array(v_collect-v_average))**2/(2*sigma**2))
p_v = p_v/p_v.sum()*5*10**(-15)/(1.6*10**(-19))
  • Position distribution of particles
r_collect = np.zeros([nx_o[0],ny_o[0],nz_o[0]])
for ix in range(nx_o[0]):
    x = ix - int(nx_o[0]/2)
    for iy in range(ny_o[0]):
        y = iy - int(ny_o[0]/2)
        for iz in range(nz_o[0]):
            z = iz - int(nz_o[0]/2)
            r = math.sqrt(x**2+y**2+z**2)
            r_collect[ix, iy, iz] = r
sigma = 3
r_v = 1/math.sqrt(2*math.pi*sigma**2)*np.exp(-(np.array(r_collect))**2/(2*sigma**2))
r_v = r_v/r_v.sum()
num_momentum_levels = 1
  • The total particle distribution setting
# 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 momentum distribution is a Gaussian function
phase_space_volume = dx*dy*dz*dpx*dpy*dpz
@jit
def setup(npx,npy,npz,nx_o,ny_o,nz_o,f,r_v,p_v):
    for ipx in range(npx):
        for ipy in range(npy):
            for ipz in range(npz):
                for ix in range(nx_o[0]):
                    for iy in range(ny_o[0]):
                        for iz in range(nz_o[0]):                       
                            # current velocity
                            f[ix,iy,iz,ipx,ipy,ipz] = r_v[ix,iy,iz]*p_v[ipx,ipy,ipz]
    return f

f[0][0][0] = setup(npx,npy,npz,nx_o,ny_o,nz_o,f[0][0][0],r_v,p_v)
# 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])

The schematic diagram below illustrates the initial position distribution of the particles in the yoz plane:

particle_dis

  • Define the electromagnetic field being used.

    The parameter โ€˜drifit_orderโ€™ can control the order of differentiation in the calculation. When โ€˜drifit_orderโ€™ is set to 1, the program performs first-order differentiation. When โ€˜drifit_orderโ€™ is set to 2, the program performs second-order differentiation.

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 VT and DT indicate whether or not to compute the Vlasov term and Drift term (1 represents computation, 0 represents no computation).
n_step = 5001
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':0., 'DT':1., 'CT':0.},\
                            BEx = BEx, BEy = BEy, BEz = BEz, BBx = BBx, BBy = BBy, BBz = BBz)
    if i_time%500 == 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()
  • Plotting the results using matplotlib :

    # spatial distribution
    import matplotlib.pyplot as plt
    xi, yi = np.mgrid[1:102:1,1:102: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)][0][0].reshape([nx_o[0],ny_o[0],nz_o[0]])[0])
    

The calculated result of drift term using second-order differentiation is shown below:

particle_Drifit2

results matching ""

    No results matching ""