Macro_quantities
The program primarily defines the charge density (ρ) and current density (J).
charged_density（）
The program primarily initializes the charge density and current density for electrons ：
Since the particle type is already specified as electrons, the initialization process focuses solely on defining the position space and momentum space.
def charged_density(f, num_particle_species, total_spatial_grids, \
masses, charges, degeneracy, \
total_phase_grids, momentum_volume_element, npx, npy, npz, \
nx, ny, nz, half_px, half_py, half_pz, \
dx, dy, dz, dpx, dpy, dpz, num_momentum_levels,\
blockspergrid_total_phase, threadsperblock):
electric_rho = cupy.zeros([total_spatial_grids])
electric_Jx = cupy.zeros([total_spatial_grids])
electric_Jy = cupy.zeros([total_spatial_grids])
electric_Jz = cupy.zeros([total_spatial_grids])
electric_rho_J_kernel[blockspergrid_total_phase, threadsperblock]\
(f, electric_rho, electric_Jx, electric_Jy, electric_Jz, \
masses, charges, degeneracy,\
total_phase_grids, num_particle_species, total_spatial_grids,\
momentum_volume_element, npx, npy, npz, nx, ny, nz, \
half_px, half_py, half_pz, \
dx, dy, dz, dpx, dpy, dpz, num_momentum_levels)
return electric_rho, electric_Jx, electric_Jy, electric_Jz
electric_rho_J_kernel()
The function primarily calculates the charge density and current density for electrons.
First, it associates the GPU block index with the position and momentum of the particles：
i_grid = cuda.grid(1)
if i_grid < total_grid:
# convert one-d index into six-d
ipz = i_grid%npz
ipz_rest = i_grid//npz
ipy = ipz_rest%npy
ipy_rest = ipz_rest//npy
ipx = ipy_rest%npx
ipx_rest = ipy_rest//npx
iz = ipx_rest%nz
iz_rest = ipx_rest//nz
iy = iz_rest%ny
iy_rest = iz_rest//ny
ix = iy_rest%nx
i_s = iz + iy*nz + ix*nz*ny
Then, it calculates the charge density ：
The particle type loop is necessary because there are multiple electrons with different masses, so their contributions are summed.
for p_type in range(num_of_particle_types):
# distribution function multiplied by volume
modified_f = f_x_p_t[0, p_type, i_grid]*momentum_volume_element[p_type]*charges[p_type]
# if modified_f is less than the numerical accuracy, turn it into 0
if abs(modified_f) < 10**(-19):
modified_f = 0.
# rho at each phase space point
# rho only considers the 0-th momentum level
cuda.atomic.add(rho, i_s, modified_f)
Finally, it calculates the current density：
for p_type in range(num_of_particle_types):
# distribution function multiplied by volume
modified_f = f_x_p_t[0, p_type, i_grid]*momentum_volume_element[p_type]*charges[p_type]
# if modified_f is less than the numerical accuracy, turn it into 0
if abs(modified_f) < 10**(-19):
modified_f = 0.
cuda.atomic.add(rho, i_s, modified_f)
px = ((ipx+0.5)*dpx[p_type] - half_px[p_type])
py = ((ipy+0.5)*dpy[p_type] - half_py[p_type])
pz = ((ipz+0.5)*dpz[p_type] - half_pz[p_type])
# energy for current grid
mp_squared = masses[p_type]**2
# p0 for current grid
p0 = math.sqrt(mp_squared*c**2+px**2+py**2+pz**2)
# vx, vy, vz
vx = 0.
vy = 0.
vz = 0.
if p0 > 10**(-19):
vx = c*px/p0
vy = c*py/p0
vz = c*pz/p0
if abs(px) < 0.5*dpx[p_type]:
vx = 0.
if abs(py) < 0.5*dpy[p_type]:
vy = 0.
if abs(pz) < 0.5*dpz[p_type]:
vz = 0.
fvx = modified_f*vx
fvy = modified_f*vy
fvz = modified_f*vz
# if fv is less than the numerical accuracy, turn it into 0
if abs(fvx) < 10**(-19):
fvx = 0.
if abs(fvy) < 10**(-19):
fvy = 0.
if abs(fvz) < 10**(-19):
fvz = 0.
cuda.atomic.add(Jx, i_s, fvx)
cuda.atomic.add(Jy, i_s, fvy)
cuda.atomic.add(Jz, i_s, fvz)