Externel-forces

The program primarily calculates the long-range electromagnetic force in the Vlasov term of the Boltzmann equation.

long_rang_electromagnetic_force

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 :

long-rang-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))

results matching ""

    No results matching ""