We can mimic the proof of this fact for open sets in $\mathbb{R}^n$. Here's a sketch of one possible way to do it.
Step 1
Instead of considering the operator $-\Delta$ we perturb it a bit my considering $L_\mu = -\Delta + \mu I$, where $I$ denotes the identity and $\mu >0$. The weak formulation of the problem $L_\mu u = f$ is then given by integrating by parts: for any $v \in H^1(M)$ we have
$$
\int_M f v = \int_M L_\mu u v = \int_M \nabla_g u \cdot \nabla_g v + \mu uv,
$$
where $\nabla_g$ is the covariant derivative determined by the metric $g$. The reason why we introduce $\mu$ is evident here: the right hand side of this now determines an inner-product on $H^1(M)$. If we had set $\mu=0$ then we would have to worry about removing the functions that are constant on each connected component of $M$.
Step 2
The above allows us to use the Lax-Milgram lemma to establish the solvability of the weak problem $L_\mu u =f$ for any $f \in (H^1(M))^\ast$. Actually, L-M is overkill and we could just use the Riesz representation. In other words, for any given $f \in (H^1)^\ast$ we can find a unique $u \in H^1$ such that
$$
\int_M \nabla_g u \cdot \nabla_g u + \mu uv = \langle f,v \rangle \text{ for all } v \in H^1.
$$
Moreover,
$$
\Vert u \Vert_{H^1} \le C \Vert f \Vert_{(H^1)^\ast}.
$$
Step 3
The above establishes that the map $L_\mu : H^1 \to (H^1)^\ast$ is an isomorphism. Now we consider the map $L_\mu^{-1}$, but restricted to $L^2(M) \hookrightarrow (H^1(M))^\ast$, i.e. we consider only $f \in L^2(M)$, which is more restrictive than $f \in (H^1)^\ast$. Doing so, we find that $L_\mu^{-1} : L^2 \to H^1$ is a bounded linear map. But we have the compact embedding $H^1 \subset \subset L^2$ due to Rellich's theorem, and consequently the map $L_\mu^{-1} : L^2 \to L^2$ is compact.
Step 4
Next we prove that $L_\mu^{-1}$ is self-adjoint on $L^2$. This follows directly, and I'll leave it as an exercise. A similar argument shows that $L_\mu^{-1}$ is a positive operator.
Step 5
We now know that $L_\mu^{-1}: L^2 \to L^2$ is a compact self-adjoint positive operator. The spectral theory of such operators now tells us that the spectrum of $L_\mu^{-1}$ consists of $0$ and a countable sequence $\rho_n$ of positive real eigenvalues such that $\rho_n \to 0$ as $n \to \infty$. Moreover, the associated eigenfunctions form an orthonormal basis of $L^2$.
Step 6
I leave it as an exercise to verify that $L_\mu^{-1} u = \rho_n u$ if and only if $L_\mu u = \rho_n^{-1} u$. This means that the spectrum of $L_\mu$ consists of a sequence $0 < \rho_n^{-1} \to \infty$. Finally, we now note that
$$
L_\mu u = \rho_n^{-1} u \Leftrightarrow -\Delta u = (\rho_n^{-1} - \mu) u.
$$
Thus the spectrum of $-\Delta$ is the sequence of real eigenvalues $\lambda_n = (\rho_n^{-1} - \mu)$, which satisfy $\lambda_n \to \infty$ as $n \to \infty$. In particular, the eigenvalues are discrete.
Given a point $p\in M$, you can recover the metric at $p$ as follows. Choose any local coordinates $(x^1,\dots,x^n)$ such that $p$ has the coordinate representation $(0,\dots,0)$. For any indices $k,l$, let $f_{kl}(x) = \tfrac 1 2 x^k x^l$. If you expand out the expression for $\Delta (f_{kl})$ and use the fact that $x^k=x^l=0$ at $p$, you'll find that $\Delta(f_{kl})(p) = g^{kl}$. Then you can use matrix inversion to recover $g_{kl}$.
If your oracle is restricted to acting only on globally defined functions, you can use a bump function to extend the coordinate functions to global smooth functions $(u^1,\dots,u^n)$ that agree with $(x^1,\dots,x^n)$ in a neighborhood of $p$, and let $f_{kl}(x) = u^k u^l$.
This is a special case of a much more general construction: if $P$ is an $m$th order scalar linear partial differential operator, the coefficients of its highest-order terms determine a coordinate-independent symmetric contravariant $m$-tensor field called the principal symbol, which can be evaluated at any point by applying $P$ to $m$-fold products of coordinate functions. The principal symbol of $\Delta$ is the associated cometric, which is the induced metric on $T^*M$.
Best Answer
Your mistake is in the computation of $\ast \, df$: The formula for the Hodge star operator is $\alpha \wedge \ast \beta = \langle \alpha , \beta \rangle \, \mathrm{vol}$, which also involves the inverse metric. The relevant consequence is
$$ \ast \, dx^i = \sum_{j=1}^n (-1)^{j-1} \sqrt{\lvert g \rvert} g^{ij} dx^1 \wedge \dotsm \wedge \widehat{dx^j} \wedge \dotsm \wedge dx^n . $$
The apparent sign error is because your expected formula is incorrect. The Hodge Laplacian as you define it is a nonnegative operator, while your expected operator is nonpositive (check by integrating by parts).