How to use the fbpic.utils.cuda.cuda_tpb_bpg_2d function in fbpic

def erase( self, fieldtype ):
        """
        Sets the field `fieldtype` to zero on the interpolation grid

        Parameter
        ---------
        fieldtype : string
            A string which represents the kind of field to be erased
            (either 'E', 'B', 'J', 'rho')
        """
        if self.use_cuda:
            # Obtain the cuda grid
            dim_grid, dim_block = cuda_tpb_bpg_2d( self.Nz, self.Nr )

            # Erase the arrays on the GPU
            if fieldtype == 'rho':
                cuda_erase_scalar[dim_grid, dim_block](self.rho)
            elif fieldtype == 'J':
                cuda_erase_vector[dim_grid, dim_block](
                      self.Jr, self.Jt, self.Jz)
            elif fieldtype == 'E':
                cuda_erase_vector[dim_grid, dim_block](
                      self.Er, self.Et, self.Ez)
            elif fieldtype == 'B':
                cuda_erase_vector[dim_grid, dim_block](
                      self.Br, self.Bt, self.Bz)
            else:
                raise ValueError('Invalid string for fieldtype: %s'%fieldtype)
        else :

# product of complexs, and the real-complex conversion is negligible.)
        if not self.use_cuda:
            # Initialize real buffer arrays on the CPU
            zero_array = np.zeros((2*Nz, Nr), dtype=np.float64 )
            self.array_in = zero_array.copy()
            self.array_out = zero_array.copy()
        else:
            # Initialize real buffer arrays on the GPU
            # The cuBlas API requires that these arrays be in Fortran order
            zero_array = np.zeros((2*Nz, Nr), dtype=np.float64, order='F')
            self.d_in = cuda.to_device( zero_array )
            self.d_out = cuda.to_device( zero_array )
            # Initialize a cuda stream (required by cublas)
            self.blas = cublas.Blas()
            # Initialize the threads per block and block per grid
            self.dim_grid, self.dim_block = cuda_tpb_bpg_2d(Nz, Nr)

interp[m].Et[:nd,:]*=damp_arr[:,np.newaxis]
                        interp[m].Ez[:nd,:]*=damp_arr[:,np.newaxis]
                        interp[m].Br[:nd,:]*=damp_arr[:,np.newaxis]
                        interp[m].Bt[:nd,:]*=damp_arr[:,np.newaxis]
                        interp[m].Bz[:nd,:]*=damp_arr[:,np.newaxis]
                        if interp[m].use_pml:
                            interp[m].Er_pml[:nd,:]*=damp_arr[:,np.newaxis]
                            interp[m].Et_pml[:nd,:]*=damp_arr[:,np.newaxis]
                            interp[m].Br_pml[:nd,:]*=damp_arr[:,np.newaxis]
                            interp[m].Bt_pml[:nd,:]*=damp_arr[:,np.newaxis]

            if self.right_proc is None:
                # Damp the fields on the CPU or the GPU
                if interp[0].use_cuda:
                    # Damp the fields on the GPU
                    dim_grid, dim_block = cuda_tpb_bpg_2d(
                        nd, interp[0].Nr )
                    for m in range(len(interp)):
                        cuda_damp_EB_right[dim_grid, dim_block](
                            interp[m].Er, interp[m].Et, interp[m].Ez,
                            interp[m].Br, interp[m].Bt, interp[m].Bz,
                            self.d_right_damp, nd)
                        if interp[m].use_pml:
                            cuda_damp_EB_right_pml[dim_grid, dim_block](
                                interp[m].Er_pml, interp[m].Et_pml,
                                interp[m].Br_pml, interp[m].Bt_pml,
                                self.d_right_damp, nd)
                else:
                    # Damp the fields on the CPU
                    damp_arr = self.right_damp
                    for m in range(len(interp)):
                        # Damp the fields in left guard cells

(environment variable NUMBA_NUM_THREADS)
        """
        # Check whether to use cuda
        self.use_cuda = use_cuda
        if (self.use_cuda is True) and (cuda_installed is False) :
            self.use_cuda = False
            print('** Cuda not available for Fourier transform.')
            print('** Performing the Fourier transform on the CPU.')

        # Check whether to use MKL
        self.use_mkl = mkl_installed

        # Initialize the object for calculation on the GPU
        if self.use_cuda:
            # Initialize the dimension of the grid and blocks
            self.dim_grid, self.dim_block = cuda_tpb_bpg_2d( Nz, Nr)

            # Initialize 1d buffer for cufft
            self.buffer1d_in = cuda.device_array(
                (Nz*Nr,), dtype=np.complex128)
            self.buffer1d_out = cuda.device_array(
                (Nz*Nr,), dtype=np.complex128)
            # Initialize the cuda libraries object
            self.fft = cufft.FFTPlan( shape=(Nz,), itype=np.complex128,
                                      otype=np.complex128, batch=Nr )
            self.blas = cublas.Blas()   # For normalization of the iFFT
            self.inv_Nz = 1./Nz         # For normalization of the iFFT

        # Initialize the object for calculation on the CPU
        else:

            # For MKL FFT

"""
        Divide the field `fieldtype` in each cell by the cell volume,
        on the interpolation grid.

        This is typically done for rho and J, after the charge and
        current deposition.

        Parameter
        ---------
        fieldtype :
            A string which represents the kind of field to be divided by
            the volume (either 'rho' or 'J')
        """
        if self.use_cuda :
            # Perform division on the GPU
            dim_grid, dim_block = cuda_tpb_bpg_2d( self.Nz, self.Nr )

            if fieldtype == 'rho':
                cuda_divide_scalar_by_volume[dim_grid, dim_block](
                        self.rho, self.d_invvol )
            elif fieldtype == 'J':
                cuda_divide_vector_by_volume[dim_grid, dim_block](
                        self.Jr, self.Jt, self.Jz, self.d_invvol )
            else:
                raise ValueError('Invalid string for fieldtype: %s'%fieldtype)
        else :
            # Perform division on the CPU
            if fieldtype == 'rho':
                self.rho *= self.invvol[np.newaxis,:]
            elif fieldtype == 'J':
                self.Jr *= self.invvol[np.newaxis,:]
                self.Jt *= self.invvol[np.newaxis,:]

def push_rho(self) :
        """
        Transfer the values of rho_next to rho_prev,
        and set rho_next to zero
        """
        if self.use_cuda :
            # Obtain the cuda grid
            dim_grid, dim_block = cuda_tpb_bpg_2d( self.Nz, self.Nr)
            # Push the fields on the GPU
            cuda_push_rho[dim_grid, dim_block](
                self.rho_prev, self.rho_next, self.Nz, self.Nr )
        else :
            # Push the fields on the CPU
            self.rho_prev[:,:] = self.rho_next[:,:]
            self.rho_next[:,:] = 0.

The number of cells by which the grid should be shifted

        shift_rho: bool, optional
            Whether to also shift the charge density
            Default: True, since rho is only recalculated from
            scratch when the particles are exchanged

        shift_currents: bool, optional
            Whether to also shift the currents
            Default: False, since the currents are recalculated from
            scratch at each PIC cycle
        """
        if grid.use_cuda:
            shift = grid.d_field_shift
            # Get a 2D CUDA grid of the size of the grid
            tpb, bpg = cuda_tpb_bpg_2d( grid.Ep.shape[0], grid.Ep.shape[1] )
            # Shift all the fields on the GPU
            shift_spect_array_gpu[tpb, bpg]( grid.Ep, shift, n_move )
            shift_spect_array_gpu[tpb, bpg]( grid.Em, shift, n_move )
            shift_spect_array_gpu[tpb, bpg]( grid.Ez, shift, n_move )
            shift_spect_array_gpu[tpb, bpg]( grid.Bp, shift, n_move )
            shift_spect_array_gpu[tpb, bpg]( grid.Bm, shift, n_move )
            shift_spect_array_gpu[tpb, bpg]( grid.Bz, shift, n_move )
            if grid.use_pml:
                shift_spect_array_gpu[tpb, bpg]( grid.Ep_pml, shift, n_move )
                shift_spect_array_gpu[tpb, bpg]( grid.Em_pml, shift, n_move )
                shift_spect_array_gpu[tpb, bpg]( grid.Bp_pml, shift, n_move )
                shift_spect_array_gpu[tpb, bpg]( grid.Bm_pml, shift, n_move )
            if shift_rho:
                shift_spect_array_gpu[tpb, bpg]( grid.rho_prev, shift, n_move )
            if shift_currents:
                shift_spect_array_gpu[tpb, bpg]( grid.Jp, shift, n_move )

use_true_rho : bool, optional
            Whether to use the rho projected on the grid.
            If set to False, this will use div(E) and div(J)
            to evaluate rho and its time evolution.
            In the case use_true_rho==False, the rho projected
            on the grid is used only to correct the currents, and
            the simulation can be run without the neutralizing ions.
        """
        # Check that psatd object passed as argument is the right one
        # (i.e. corresponds to the right mode)
        assert( self.m == ps.m )

        if self.use_cuda :
            # Obtain the cuda grid
            dim_grid, dim_block = cuda_tpb_bpg_2d( self.Nz, self.Nr, 1, 16)
            # Push the fields on the GPU
            if ps.V is None:
                # With the standard PSATD algorithm
                if self.use_pml:
                    # Push the PML split component
                    cuda_push_eb_pml_standard[dim_grid, dim_block](
                        self.Ep_pml, self.Em_pml, self.Bp_pml, self.Bm_pml,
                        self.Ez, self.Bz, ps.d_C, ps.d_S_w,
                        self.d_kr, self.d_kz, self.Nz, self.Nr )
                # Push the regular fields
                cuda_push_eb_standard[dim_grid, dim_block](
                    self.Ep, self.Em, self.Ez, self.Bp, self.Bm, self.Bz,
                    self.Jp, self.Jm, self.Jz, self.rho_prev, self.rho_next,
                    ps.d_rho_prev_coef, ps.d_rho_next_coef, ps.d_j_coef,
                    ps.d_C, ps.d_S_w, self.d_kr, self.d_kz, ps.dt,
                    use_true_rho, self.Nz, self.Nr )

if method == 'replace':
            nz_start = self.n_guard
            nz_end = 2*self.n_guard
        if method == 'add':
            nz_start = 0
            nz_end = 2*self.n_guard
        # Whether or not to send to the left or right neighbor
        copy_left = (self.left_proc is not None)
        copy_right = (self.right_proc is not None)
        Nz = grid_r[0].shape[0]

        # When using the GPU
        if use_cuda:

            # Calculate the number of blocks and threads per block
            dim_grid_2d, dim_block_2d = cuda_tpb_bpg_2d(
                nz_end - nz_start, self.Nr )

            if before_sending:
                # Copy the inner regions of the domain to the buffers
                for m in range(self.Nm):
                    if pml_r is None:
                        # Copy only the regular components
                        copy_vec_to_gpu_buffer[ dim_grid_2d, dim_block_2d ](
                            self.d_send_l[exchange_type],
                            self.d_send_r[exchange_type],
                            grid_r[m], grid_t[m], grid_z[m], m,
                            copy_left, copy_right, nz_start, nz_end )
                    else:
                        # Copy regular components + PML components
                        copy_pml_to_gpu_buffer[ dim_grid_2d, dim_block_2d ](
                            self.d_send_l[exchange_type],

def filter(self, fieldtype) :
        """
        Filter the field `fieldtype`

        Parameter
        ---------
        fieldtype : string
            A string which represents the kind of field to be filtered
            (either 'E', 'B', 'J', 'rho_next' or 'rho_prev')
        """
        if self.use_cuda :
            # Obtain the cuda grid
            dim_grid, dim_block = cuda_tpb_bpg_2d( self.Nz, self.Nr )
            # Filter fields on the GPU
            if fieldtype == 'J' :
                cuda_filter_vector[dim_grid, dim_block](
                        self.Jp, self.Jm, self.Jz, self.Nz, self.Nr,
                        self.d_filter_array_z, self.d_filter_array_r )
            elif fieldtype == 'E' :
                cuda_filter_vector[dim_grid, dim_block](
                        self.Ep, self.Em, self.Ez, self.Nz, self.Nr,
                        self.d_filter_array_z, self.d_filter_array_r )
            elif fieldtype == 'B' :
                cuda_filter_vector[dim_grid, dim_block](
                        self.Bp, self.Bm, self.Bz, self.Nz, self.Nr,
                        self.d_filter_array_z, self.d_filter_array_r )
            elif fieldtype in ['rho_prev', 'rho_next',
                                'rho_next_z', 'rho_next_xy']:
                spectral_rho = getattr( self, fieldtype )