This paper develops a novel approach to the problem of source current identification for the diffusion equation in connection with geophysical self-potential measurements. The problem is split into two subproblems: (a) the scalar source identification, and (b) solution of the divergence equation. For subproblem (a), we design an algorithm for reconstructing the scalar source function, which does not require solving the Fredholm integral equation of the first kind. Instead, the problem is reformulated as a linear operator equation, which is solved by a projection method. The dimension of the subspace, in which the source function is sought, is independent of the dimension of the forward problem, leading to reduction of the size of the inverse problem. Numerical experiments with exact and noisy data are presented. For subproblem (b), the divergence equation was posed as a minimization problem, which, by means of Lagrangian formalism, was reduced to a system of partial differential equation of Stokes type with a unique solution. To demonstrate how this framework can be used in a practical application, we implemented the algorithm in the two-dimensional physical space, using a finite-different discretization on staggered grids.