Charge deposit in a pic code

dolfinx
1

Hello everyone,

I have written a beam simulation code using DOLFINx, but I am not sure if the algorithm for calculating charge density from charge deposition is correct. I would appreciate any advice or suggestions. The code is attached below.

def deposit_charge_dg0_parallel(particles, V_rho):

    """

    自洽电荷沉积 (P -> rho) - 采用与对称脚本一致的 DG0 空间

    确保非结构网格单元中心的电荷质心绝对对称

    """

    mesh = V_rho.mesh

    rho = fem.Function(V_rho)

    rho.x.array[:] = 0.0




    if len(particles) > 0:

        points = np.array([p.x for p in particles], dtype=np.float64)

        bb_tree = geometry.bb_tree(mesh, mesh.topology.dim)

        

        # 批量确定点所属单元

        cell_candidates = geometry.compute_collisions_points(bb_tree, points)

        colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points)




        dofmap = V_rho.dofmap

        for i, p in enumerate(particles):

            cells = colliding_cells.links(i)

            if len(cells) > 0:

                p.cell_id = cells[0]  # 缓存单元 ID

                local_dofs = dofmap.cell_dofs(p.cell_id)

                dof = local_dofs[0]

                rho.x.array[dof] += p.q

            else:

                p.cell_id = -1




    # MPI 并行求和汇总

    rho.x.scatter_reverse(la.InsertMode.add)

    

    # 转换为密度 rho = Q/V

    vol_func = fem.Function(V_rho)

    vol_expr = fem.Expression(ufl.CellVolume(mesh), V_rho.element.interpolation_points())

    vol_func.interpolate(vol_expr)

    

    vol_array = vol_func.x.array

    rho_array = rho.x.array

    mask = vol_array > 1e-30

    rho_array[mask] /= vol_array[mask]

    rho_array[~mask] = 0.0

    

    rho.x.scatter_forward()

    return rho