Do contravariant and covariant partial derivatives commute in GR? I'm considering something like this: $\partial_{\mu}\partial^{\nu}A$ .  I feel like we should be able to commute the derivatives so: $\partial_{\mu}\partial^{\nu}A = \partial^{\nu}\partial_{\mu}A$. 
However, surely we can write this as: $\partial_{\mu}\partial^{\nu}A=\partial_{\mu}(g^{\nu \rho}\partial_{\rho}A)$ where immediately it would appear we run into trouble since the partial derivative of the metric is non-zero.
So is it true that we can't actually commute covariant and contravariant derivatives in GR or is there a flaw (quite probably) in my logic?
Note: This same logic also seems to imply that $\partial^{\mu}\partial_{\mu}\ne \partial_{\mu}\partial^{\mu}$, which seems very strange to me.
 A: Covariant derivatives' commutator defines the Riemann tensor, viz. $[\nabla_\mu,\,\nabla_\nu]=R_{\mu\nu\rho\sigma}\nabla^\sigma$. Contracting the left-hand indices obtains $0$, as the Riemann tensor is antisymmetric in its leftmost indices. However, in general raising one index without contraction gives a nonzero commutator for the covariant derivatives. Of course, since covariant derivatives are tensors, so are their commutators.
The problem you've found with partial derivatives with upstairs indices actually provide a motivation for covariant derivatives, but we can explain that with the even simpler question of how these derivatives are defined. Does $\partial^\nu$, whatever that is, act by applying $\partial_\mu$ before or after multiplying by $g^{\mu\nu}$? Clearly, the problem is that partial derivatives are not "metric compatible", which in view of the Leibniz law means they don't annihilate the metric tensor, whereas covariant derivatives are designed to do exactly that.
A: They shouldn't commute, that they seem like they should is mostly a superficial result of the notation.
In the following we rewrite these derivatives in a way such that it's manifest that they're not the same thing. The idea is simply to make sense of these in terms of operations that are a little more global than shifting around terms in expressions on a tangent space. 
Claim
Let $(M, g)$ be a Riemannian manifold, $f \in C^\infty(M)$. On some local coordinate $x^\mu$ defined on open $U \in M$:
$$
\partial_v \partial^\mu f = d\left(~dx^\mu[\nabla f]~\right)[\partial_\nu]
$$
$$
\partial^\mu \partial_v  f = dx^\mu[~\nabla(~df[\partial_v]~)~]
$$
where


*

*$d: C^\infty(M) \to \Omega^1(M) \equiv \Gamma(T^*M)$ ($\Gamma$ denote space of smooth sections in a bundle) is the De Rham differential.

*$\nabla: C^\infty(f) \to \Gamma(TM)$ is the usual gradient, defined by the property


$$
df[v] = \langle\nabla f, v\rangle
$$
Proof of the claim is just a direct computation writing the RHS of the expressions above out in components. A useful identity is the local expression for the gradient: $dx^\mu[\nabla f] = g^{\mu\nu}\partial_v f$ (which can prove as an exercise or look it up in any standard Riemannian geometry text).
On the sum If you're trying to compute a Laplacian, the correct version is $\partial_\mu \partial^\mu$, see, for example, the brief derivation of local formula for the Laplace-Beltrami operator on wikipedia.
