Externel-forces
The program primarily calculates the long-range electromagnetic force in the Vlasov term of the Boltzmann equation.
External_forces（）
The program primarily initializes the long-range electromagnetic force
def External_forces(masses, charges,\
total_phase_grids, num_particle_species, \
npx,npy,npz,nx,ny,nz, \
half_px, half_py, half_pz, \
dx, dy, dz, dpx, dpy, dpz, \
blockspergrid_total_phase, threadsperblock, number_momentum_levels,\
Ex, Ey, Ez, Bx, By, Bz):
# allocate a space in GPU to save the electromagnetic forces
Fx = cupy.zeros([number_momentum_levels, num_particle_species, total_phase_grids])
Fy = cupy.zeros([number_momentum_levels, num_particle_species, total_phase_grids])
Fz = cupy.zeros([number_momentum_levels, num_particle_species, total_phase_grids])
EM_force_kernel[blockspergrid_total_phase, threadsperblock]\
(Fx, Fy, Fz, masses, charges,\
total_phase_grids, num_particle_species, \
npx,npy,npz,nx,ny,nz, \
half_px, half_py, half_pz, \
dx, dy, dz, dpx, dpy, dpz, \
Ex, Ey, Ez, Bx, By, Bz, number_momentum_levels, c)
return Fx, Fy, Fz
EM_force_kernel()
The function primarily calculates the long-range electromagnetic force；
First, we establish a one-to-one correspondence between the GPU block indices and the position space and momentum space；
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
Next, we perform calculations for the long-range electromagnetic force ：
for p_type in range(num_of_particle_types):
mp_squared = masses[p_type]**2
# loop through all momentum levels
for i_level in range(number_momentum_levels):
px = ((ipx+0.5)*dpx[p_type] - half_px[p_type])/(npx**i_level)
py = ((ipy+0.5)*dpy[p_type] - half_py[p_type])/(npy**i_level)
pz = ((ipz+0.5)*dpz[p_type] - half_pz[p_type])/(npz**i_level)
charge = charges[p_type]
# 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.
Fx[i_level,p_type,i_grid] = charge*(Efx + (vy*Bfz - vz*Bfy))
Fy[i_level,p_type,i_grid] = charge*(Efy + (vz*Bfx - vx*Bfz))
Fz[i_level,p_type,i_grid] = charge*(Efz + (vx*Bfy - vy*Bfx))