Plasma

This program mainly deals with the temporal evolution of multi-particle systems, including：

  • Parallel computation of single-particle distribution function；
  • Collision term Calculation；
  • Drifit term Calculation；
  • Charge density and current density calculation；
  • Based on Jefimenko’s equations, the electric and magnetic fields can be solved；
  • External Electromagnetic Force Calculation.

Plasma

To understand the application background and principles of the program, you can refer to the following article.

Towards a full solution of the relativistic Boltzmann equation for quark-gluon matter on GPUs

The advantages of the program

  • Parallelize the program using CUDA；
  • Realized the conversion between flexible unit system and international unit system；
  • Use Ray framework to efficiently call multiple servers.；
  • Create a clear framework to integrate and call each sub-module in an orderly manner；
  • Ensure momentum conservation in the simulation process by using the “symmetric sampling” method；
  • Reduce the diffusion phenomenon in the simulation process by using the Unconditionally Positivity Preserving Finite Difference Scheme (UPFD) and Flux Limiter method；
  • Utilize GPU for efficient integration.

The structure of the program

  • Init
    • collision_type_for_all_species（）
    • Plasma（）
    • plasma_single_GPU（）
    • check_input_legacy（）
    • EMsolver（）
  • Main
    • proceed_one_step（）
      • PF_PT（）
        • external_forces（）
        • Vlasov（）
        • Collision（）
        • Drift（）
      • update_f_single_GPU（）
    • return_self（）
  • other_calculation
    • find_largest_time_steps（）
    • boundary_coding（）

The operational logic of the entire program is depicted in the ensuing diagram:

flowchart

collision_type_for_all_species（）

The main purpose of this function is to establish and configure collision parameters.：

  • 1、Provides an index of the particle types used in the program；
particle_order = 'u (0), d (1), s (2), ubar (3), dbar (4), sbar (5), fluon (6)'
  • 2、Provides an index of the type of reaction performed by the particle；
collision_type22=np.array([[1,0,0,0,0,5,2,3,3,4], ...)
  • 3、Provides an index of the particles used in the reaction type.
flavor22=np.array([[[0,0,0,0],[1,0,1,0],[2,0,2,0],[4,0,4,0],[5,0,5,0],[6,0,6,0], \
                    [3,0,3,0],[4,1,3,0],[5,2,3,0],[6,6,3,0]], ...)

check_input_legacy(f, dt, \
                   nx_o, ny_o, nz_o, dx, dy, dz, boundary_configuration, \
                   x_left_bound_o, y_left_bound_o, z_left_bound_o, \
                   npx, npy, npz, half_px, half_py, half_pz,\
                   masses, charges,\
                   sub_region_relations,\
                   num_gpus_for_each_region,\
                   num_samples,\
                   flavor, collision_type, particle_type,\
                   degeneracy, expected_collision_type)

For a specific and detailed description, please refer to collision_type_for_all_species(), where the function is thoroughly explained.

Plasma（）

This function primarily constructs the class, initializes the pertinent parameters employed by the program, and subsequently transfers said parameters to the Graphics Processing Unit (GPU).

  • The determination of the input parameter type is conducted prior to the initialization process.：
check_input_legacy(f, dt, \
                   nx_o, ny_o, nz_o, dx, dy, dz, boundary_configuration, \
                   x_left_bound_o, y_left_bound_o, z_left_bound_o, \
                   npx, npy, npz, half_px, half_py, half_pz,\
                   masses, charges,\
                   sub_region_relations,\
                   num_gpus_for_each_region,\
                   num_samples,\
                   flavor, collision_type, particle_type,\
                   degeneracy, expected_collision_type)

For a specific and detailed description, please refer to Input parameters, where the function is thoroughly explained.

  • The maximal temporal increment necessary for the electromagnetic wave propagation from position A to position B prior to the initialization of EMsolver() is determined:

len_time

def find_largest_time_steps(dx_o, dy_o, dz_o, x_left_boundary_o, y_left_boundary_o, z_left_boundary_o, \
                            dx_s, dy_s, dz_s, x_left_boundary_s, y_left_boundary_s, z_left_boundary_s, \
                            x_grid_size_o, y_grid_size_o, z_grid_size_o, \
                            x_grid_size_s, y_grid_size_s, z_grid_size_s,\
                            dt, c):

plasma_single_GPU（）

The primary function of the program involves receiving the parameters communicated by the Central Processing Unit (CPU) and transmitting them to the Graphics Processing Unit (GPU). Moreover, the boundary conditions are handled in accordance with distinct blocks.

Where bound_indicator is: indicator, indexes represent different surfaces of the block：

0→ back; 1→ front; 2→ left; 3→ right side; 4→ below; 5→ Above；

sub_region_relations = {'indicator': [[0,3,4],[0,2,4],[0,3,5],[0,2,5],\
                                      [1,3,4],[1,2,4],[1,3,5],[1,2,5]]}

To facilitate the interaction of block boundaries, the ghost boundaries method is employed. This method entails the following steps:

  • Initially, the bound_indicator is utilized to identify potential overlaps between the current block and other blocks;
  • Subsequently, the boundary conditions along the negative axis of the coordinate system are adjusted;
  • Finally, the size of the calculation area is modified by subtracting the portion that represents the “coincidence surface.
for i_reg in range(number_regions):
    # index for rho/J in generating EM
    self.id_xl, self.id_xr, self.id_yl, self.id_yr, self.id_zl, self.id_zr = 0, self.nx, 0, self.ny, 0, self.nz
    # the source region needs to be shrinked, according to the ghost boundaries
    self.x_left_bound_adj, self.y_left_bound_adj, self.z_left_bound_adj = self.x_left_bound, self.y_left_bound, self.z_left_bound
    self.nx_adj, self.ny_adj, self.nz_adj = self.nx, self.ny, self.nz

    # check the -x, -y, -z, +x, +y, +z directions           
    if 0 in bound_indicator:
       self.nx_adj -= 1
       self.x_left_bound_adj += self.dx
       self.id_xl += 1
    if 1 in bound_indicator:
       self.nx_adj -= 1
       self.id_xr -= 1
    if 2 in bound_indicator:
       self.ny_adj -= 1
       self.y_left_bound_adj += self.dy
       self.id_yl += 1
    if 3 in bound_indicator:
       self.ny_adj -= 1
       self.id_yr -= 1
    if 4 in bound_indicator:
       self.nz_adj -= 1
       self.z_left_bound_adj += self.dz
       self.id_zl += 1
    if 5 in bound_indicator:
       self.nz_adj -= 1  
       self.id_zr -= 1
  • The initialization process precedes the electromagnetic (EM) calculation and involves the transmission of parameters:
self.EMsolver.append(EMsolver(self.len_time_snapshots[i_reg], \
                                          self.nx_o[i_reg], self.ny_o[i_reg], self.nz_o[i_reg], \
                                          self.nx_adj, self.ny_adj, self.nz_adj,\
                                          self.dx_o[i_reg], self.dy_o[i_reg], self.dz_o[i_reg], \
                                          self.x_left_bound_o[i_reg], self.y_left_bound_o[i_reg], \
                                          self.z_left_bound_o[i_reg], \
                                          self.dx, self.dy, self.dz, \
                                          self.x_left_bound_adj, self.y_left_bound_adj, self.z_left_bound_adj, \
                                          self.dt, self.epsilon0, self.c))

EMsolver()

The primary function of the program involves the initialization of electromagnetic (EM) calculation parameters. Additionally, it performs simultaneous computations of Coff_E and Coff_B。

For a specific and detailed description, please refer to EMsolver, where the function is thoroughly explained.

  • Calculation Coff_EandCoff_B。

coeffEB

self.coeff_E = 1/(4*math.pi*epsilon0)*dx_s*dy_s*dz_s
self.coeff_B = -1/(4*math.pi*epsilon0*c**2)*dx_s*dy_s*dz_s

boundary_coding()

boundary

  • surface_index_taken in the code indicates a boundary that will be used elsewhere（The inner part of the picture, where the arrow starts）；
  • The surface index in the code represents the boundary that will be replaced（In the outer part of the image, the arrow points to the position）；
def boundary_coding(nx, ny, nz, npx, npy, npz, \
                    number_momentum_levels, num_particle_species, \
                    bound_indicator):

    # these are the boundaries that need to be taken
    surface_index_taken = np.array([[1,    2,    0,    ny,   0,    nz],\
                                    [nx-2, nx-1, 0,    ny,   0,    nz],\
                                    [0,    nx,   1,    2,    0,    nz],\
                                    [0,    nx,   ny-2, ny-1, 0,    nz],\
                                    [0,    nx,   0,    ny,   1,    2],\
                                    [0,    nx,   0,    ny,   nz-2, nz-1]])

    # these are the boundaries that need to be exchanged
    surface_index = np.array([[0,    1,    0,    ny, 0,    nz],\
                              [nx-1, nx,   0,    ny, 0,    nz],\
                              [0,    nx,   0,    1,  0,    nz],\
                              [0,    nx,   ny-1, ny, 0,    nz],\
                              [0,    nx,   0,    ny, 0,    1],\
                              [0,    nx,   0,    ny, nz-1, nz]])
    
    return surface_index_taken, surface_index

Main()——The main part of the program

proceed_one_step（）——Plasma

The primary functions are delineated as follows:

  • Evaluation of all blocks and determination of the corresponding time increment (dt) deemed suitable；
  • Time evolution in accordance with the aforementioned dt:
    • Computation of charge density and current density via the distribution function.；
    • Calculation of electromagnetic fields employing the Jefimenko’s equation.；
    • Estimation of remote electromagnetic forces；
    • Computation of Drifit terms；
    • Computation of Collision terms；
    • Updating the distribution function——f(t+1) = dt*(f(t) - Vt - Dt + Ct)
  • Exchange coordinate boundaries, electromagnetic fields.
def proceed_one_step(self, i_time, total_time_steps, processes = {'VT':1., 'DT':1., 'CT':1.}, \
                         BEx = 0., BEy = 0., BEz = 0., BBx = 0., BBy = 0., BBz = 0.):

        whether_terminate = []
        # evaluate the distribution function in all regions
        for i_reg in range(self.number_regions):
            whether_terminate.append(self.plasmas[i_reg].proceed_one_step.remote(i_time, total_time_steps, self.global_dt, try_evaluate = True, processes = processes, BEx = BEx[i_reg], BEy = BEy[i_reg], BEz = BEz[i_reg], BBx = BBx[i_reg], BBy = BBy[i_reg], BBz = BBz[i_reg]))
            
        for i_reg in range(self.number_regions):
            whether_terminate[i_reg] = ray.get(whether_terminate[i_reg])

        # terminate if illlegal results have been obtained
        if any(whether_terminate) == True:
            raise AssertionError("The code has been terminated since the obtained distribution is unphysical!") 
        
        # exchange the boundaries after evaluation
        if self.number_regions > 1:
            self.exchange_boundaries()
            
        # exchange EM fields for Vlasov term
        if processes['VT']>0.5:
            self.exchange_EM_fields()

        # load stream_ticker
        # stream_ticker transfers a small data from GPU to CPU
        # this gives a signal to the GPU that whether each time step is finished 
        self.acquire_values("stream_ticker")

proceed_one_step（）——Plasma_Single_GPU

This function mainly passes the electromagnetic fields to the GPU and also makes distribution function update program calls；

def proceed_one_step(self, current_time_step, total_time_steps, \
                     dt, try_evaluate, \
                     processes = {'VT':1., 'DT':1., 'CT':1.},\
                     BEx = 0., BEy = 0., BEz = 0., BBx = 0., BBy = 0., BBz = 0.):

    # sending the background EM fields to GPU
    BEx, BEy, BEz, BBx, BBy, BBz = cupy.asarray(BEx), cupy.asarray(BEy), cupy.asarray(BEz),\
    cupy.asarray(BBx), cupy.asarray(BBy), cupy.asarray(BBz)

    # evalute new distribution
    whether_terminate = self.PF_PT(processes, BEx, BEy, BEz, BBx, BBy, BBz, \
                                   dt, current_time_step, total_time_steps)
    return whether_terminate

PF_PT（）

The fundamental principle governing this function primarily relies on input parameters to discern the requisite terms, which predominantly encompass the drift term, Vlasov term, external force, and collision term.

DT----Drift()      
VT----Vlasov()
VT and eval_External_Forces----external_forces()
CT----Collision()

Through parameter selection, the program can solve for partial results, thus helping the user to refine the scope and focus of the study.

def PF_PT(self, processes, BEx, BEy, BEz, BBx, BBy, BBz, dt, current_time_step, total_time_steps):

    # calculate the external forces
    if processes['VT']>0.5:
        BEx, BEy, BEz, BBx, BBy, BBz = cupy.asarray(BEx), cupy.asarray(BEy), cupy.asarray(BEz),\
        cupy.asarray(BBx), cupy.asarray(BBy), cupy.asarray(BBz)
        self.external_forces(BEx, BEy, BEz, BBx, BBy, BBz)
        self.rho_J_to_EB()
    else:
        self.Fx, self.Fy, self.Fz = cupy.zeros_like(self.f), cupy.zeros_like(self.f), cupy.zeros_like(self.f)

    # if 'CT' is instantiated, calculate the collision terms at all levels
    if processes['CT'] > 0.5:
        # calculate collision term
        self.CT = Collision_term(self.flavor, self.collision_type, \
                                 self.particle_type, self.degeneracy, self.num_samples,\
                                 self.f, self.masses, self.total_phase_grids, self.rng_states,\
                                 self.num_particle_species, \
                                 self.npx, self.npy, self.npz, self.nx, self.ny, self.nz, self.dp, \
                                 self.half_px, self.half_py, self.half_pz, \
                                 self.dpx, self.dpy, self.dpz, \
                                 self.blockspergrid_total_phase, self.threadsperblock, \
                                 self.middle_npx, self.middle_npy, self.middle_npz, \
                                 self.number_momentum_levels,\
                                 self.hbar,self.c,self.lambdax,self.allowed_collision_type)

    else:
    	# otherwise set collision to be zero
        self.CT = cupy.zeros_like(self.f)

    # evaluate Drift and Vlasov terms
    whether_terminate = self.DT_VT(self.f, dt, current_time_step, processes)
    return whether_terminate

rho_J_to_EB（）

Initially, the charge density and current density are computed for the given context：

self.electric_rho, self.electric_Jx, self.electric_Jy, self.electric_Jz = \
charged_density(self.f, self.num_particle_species, self.total_spatial_grids, \
                self.masses, self.charges, \
                self.total_phase_grids, self.momentum_volume_element, self.npx, self.npy, self.npz, \
                self.nx, self.ny, self.nz, self.half_px, self.half_py, self.half_pz, \
                self.dx, self.dy, self.dz, self.dpx, self.dpy, self.dpz, self.number_momentum_levels,\
                self.blockspergrid_total_phase, self.threadsperblock, self.c)

For a specific and detailed description, please refer to Macro_quantities, where the function is thoroughly explained.

Then, based on the charge density and current density, the electric and magnetic fields are solved using the Jefimenko’s equations.

if self.rho_J_method == "raw":            

    for i_reg in range(self.number_regions):
        self.EMsolver[i_reg].Jefimenko_solver(self.electric_rho.reshape([self.nx, self.ny, self.nz])[self.id_xl:self.id_xr,  \
                                                                           self.id_yl:self.id_yr,self.id_zl:self.id_zr].flatten(),self.electric_Jx.reshape([self.nx, self.ny, self.nz])[self.id_xl:self.id_xr,  \
                                                                           self.id_yl:self.id_yr,self.id_zl:self.id_zr].flatten(),self.electric_Jy.reshape([self.nx, self.ny, self.nz])[self.id_xl:self.id_xr,   \
                                                                           self.id_yl:self.id_yr,self.id_zl:self.id_zr].flatten(),self.electric_Jz.reshape([self.nx, self.ny, self.nz])[self.id_xl:self.id_xr, self.id_yl:self.id_yr,self.id_zl:self.id_zr].flatten(), self.quasi_neutral)
        self.Ex_dis[i_reg], self.Ey_dis[i_reg], self.Ez_dis[i_reg], \
        self.Bx_dis[i_reg], self.By_dis[i_reg], self.Bz_dis[i_reg] = \
        self.EMsolver[i_reg].GEx.get(), self.EMsolver[i_reg].GEy.get(), self.EMsolver[i_reg].GEz.get(),\
        self.EMsolver[i_reg].GBx.get(), self.EMsolver[i_reg].GBy.get(), self.EMsolver[i_reg].GBz.get()

    elif self.rho_J_method == 'quasi_neutral':
        pass
    else:
        raise AssertionError("Plasma method must be raw or quasi_neutral!")

For a specific and detailed description, please refer to EMsolver, where the function is thoroughly explained.

external_forces（）

  • remote electromagnetic force calculation.
def external_forces(self, BEx, BEy, BEz, BBx, BBy, BBz):

    self.Ex_total, self.Ey_total, self.Ez_total, self.Bx_total, self.By_total, self.Bz_total = self.Ex.sum(axis=0)+BEx, self.Ey.sum(axis=0)+BEy, self.Ez.sum(axis=0)+BEz, \
    self.Bx.sum(axis=0)+BBx, self.By.sum(axis=0)+BBy, self.Bz.sum(axis=0)+BBz
    self.Fx, self.Fy, self.Fz = \
    External_forces(self.masses, self.charges,\
                    self.total_phase_grids, self.num_particle_species, \
                    self.npx,self.npy,self.npz,self.nx,self.ny,self.nz, \
                    self.half_px, self.half_py, self.half_pz, \
                    self.dx, self.dy, self.dz, self.dpx, self.dpy, self.dpz, \
                    self.blockspergrid_total_phase, self.threadsperblock, self.number_momentum_levels,\
                    self.Ex_total, self.Ey_total, self.Ez_total, self.Bx_total, self.By_total, self.Bz_total, self.c)
    print(self.By_total.sum(), self.Bz_total.sum(),cupy.abs(self.Fx).sum(),cupy.abs(self.Fy).sum(),cupy.abs(self.Fz).sum())

For a specific and detailed description, please refer to Externel-forces, where the function is thoroughly explained.

Collision()

  • Collision term calculation:
self.CT = Collision_term(self.flavor, self.collision_type, \
                         self.particle_type, self.degeneracy, self.num_samples,\
                         self.f, self.masses, self.total_phase_grids, self.rng_states,\
                         self.num_particle_species, \
                         self.npx, self.npy, self.npz, self.nx, self.ny, self.nz, self.dp, \
                         self.half_px, self.half_py, self.half_pz, \
                         self.dpx, self.dpy, self.dpz, \
                         self.blockspergrid_total_phase, self.threadsperblock, \
                         self.middle_npx, self.middle_npy, self.middle_npz, \
                         self.number_momentum_levels,\
                         self.hbar,self.c,self.lambdax,self.allowed_collision_type)

For a specific and detailed description, please refer to Collision_term(), where the function is thoroughly explained.

Vlasov_Drift()

  • Calculation of the Vlasov and Drift terms：
def DT_VT(self, f_local, dt, current_time_step, processes):

    self.f = Drift_Vlasov_terms(f_local, self.Fx, self.Fy, self.Fz, \
                                self.masses, self.total_phase_grids, self.num_particle_species, \
                                self.npx, self.npy, self.npz, self.nx, self.ny, self.nz,\
                                self.half_px, self.half_py, self.half_pz, \
                                self.dx, self.dy, self.dz, self.dpx, self.dpy, self.dpz, self.number_momentum_levels,\
                                self.x_bound_config, self.y_bound_config, self.z_bound_config, \
                                self.blockspergrid_total_phase, self.threadsperblock,\
                                self.CT, dt, self.c, current_time_step, \
                                processes['DT'], processes['VT'], self.drift_order)

    if current_time_step%time_stride_back==0 and current_time_step != 0 :
            # check is the obtained self.fp is legal.
        whether_legal, illegal_type, value, which_species, which_position = check_legacy_of_distributions(self.f, self.particle_type, self.hbar, self.dt, self.degeneracy)
        if whether_legal == 1:
            self.last_f = self.f 
        else:
            self.f = self.last_f  
            return 0        
    return 1

For a specific and detailed description, please refer to Vlasov_Drifit_terms, where the function is thoroughly explained.

acquire_values()

The function extracts the result of the calculation by setting the parameters

def acquire_values(self, quantity):
        
        quantities_required = {}
        for i_reg in range(self.number_regions):   
            quantities_required[i_reg] = ray.get(self.plasmas[i_reg].return_self.remote(quantity))
          
        return quantities_required

return_self()

The function is primarily responsible for data transmission and acquisition, transferring the calculation results from Plasma_single_GPU to Plasma.

def return_self(self, process):

    if 'Collision term' == process:
        return self.CT.get()

    if 'Distribution' == process:
        return self.f.get()

    if 'Forces' == process:
        return self.Fx.get(), self.Fy.get(), self.Fz.get()

    if "Electric rho/J" == process:
        return self.electric_rho.get(), self.electric_Jx.get(), self.electric_Jy.get(), self.electric_Jz.get()

    if "number_rho/J" == process:
        self.get_rho_J()
        return self.number_rho.get(), self.number_Jx.get(), self.number_Jy.get(), self.number_Jz.get()

    if "EM fields on current region" == process:
        return self.Ex_total.get(), self.Ey_total.get(), self.Ez_total.get(), self.Bx_total.get(), self.By_total.get(), self.Bz_total.get()

    if "EM fields to other regions" == process:
        return self.Ex_dis, self.Ey_dis, self.Ez_dis, self.Bx_dis, self.By_dis, self.Bz_dis
    if "stream_ticker" == process:
        return self.stream_ticker.get()

results matching ""

    No results matching ""