Euler-Lagrange equation in curved spacetime: how to calculate $\partial \mathcal{L}/\partial(\nabla_\nu A_\mu)$? My problem is mostly with notation. When we apply the Euler-Lagrange equation to the Lagrangian scalar $$\mathcal{L}= -\frac{1}{4}F^{\mu\nu}F_{\mu\nu}$$
Where $$F_{\mu\nu} = \partial_\mu A_\nu -\partial_\nu A_\mu= \nabla_\mu A_\nu -\nabla_\nu A_\mu$$ and I start with $$\frac{\partial\mathcal{L}}{\partial A_\mu}=\nabla_\nu\left(\frac{\partial \mathcal{L}}{\partial(\nabla_\nu A_\mu)}\right).$$ 
Now I rewrite the lagrangian as $\mathcal{L} = -\frac{1}{4}g^{\mu\alpha}g^{
\nu\beta}(\nabla_\mu A_\nu-\nabla_\nu A_\mu)(\nabla_\alpha A_\beta-\nabla_\beta A_\alpha)$
$$\mathcal{L}= -\frac{1}{4}g^{\mu\alpha}g^{\nu\beta}\left[(\nabla_\mu A_\nu)(\nabla_\alpha A_\beta)-(\nabla_\mu A_\nu)(\nabla_\beta A_\alpha)-(\nabla_\nu A_\mu)(\nabla_\alpha A_\beta)+(\nabla_\nu A_\mu)(\nabla_\beta A_\alpha)
\right]$$
My problem is calculating $\frac{\partial \mathcal{L}}{\partial(\nabla_\nu A_\mu)}$, how are the rules of partial derivatives generalised to covariant ones?
Also, is $\frac{\partial\mathcal{L}}{\partial A_\mu}=0$, is that correct ?
 A: There is an eazier way to derive the equations of motion(without calculating this nasty derivative), using calculus of variations(varying with respect to $A_{ν}$) and discarding surface terms: 
$\delta S=0 \Rightarrow \delta \int d^{4}x\sqrt{-g}F^{μν}F_{μν}= 
\int d^{4}x\sqrt{-g}\delta F^{μν}F_{μν} + \int d^{4}x\sqrt{-g}F^{μν}\delta F_{μν}= 2\int d^{4}x\sqrt{-g}F^{μν}\delta F_{μν}= 
2\int d^{4}x\sqrt{-g}F^{μν}(\nabla_{\mu}δA_{\nu} - \nabla_{\nu}δA_{\mu}) = 
4 \int d^{4}x\sqrt{-g}F^{μν}(\nabla_{\mu}δA_{\nu})$ 
So we have:
$\int d^{4}x\sqrt{-g}\nabla_{\mu}(F^{μν}δA_{\nu}) = \int d^{4}x\sqrt{-g}F^{μν}\nabla_{\mu}(δA_{\nu}) + \int d^{4}x\sqrt{-g}(\nabla_{\mu}F^{μν})δA_{\nu} \Rightarrow
 \int d^{4}x\sqrt{-g}F^{μν}\nabla_{\mu}(δA_{\nu}) = -  \int d^{4}x\sqrt{-g}(\nabla_{\mu}F^{μν})δA_{\nu} =0\Rightarrow $
$$\nabla_{\mu}F^{μν}=0$$
Explanations: 
$\bullet$ $ \int d^{4}x\sqrt{-g}\delta F^{μν}F_{μν} + \int d^{4}x\sqrt{-g}F^{μν}\delta F_{μν}= 2\int d^{4}x\sqrt{-g}F^{μν}\delta F_{μν}$:
$\delta F^{μν}F_{μν} + F^{μν}\delta F_{μν} = F^{μν}\delta F_{μν} + \delta (g^{αμ}g^{βν}F_{αβ})F_{μν} = F^{μν}\delta F_{μν} + δF_{αβ}F^{αβ}$ (since we vary w.r.t $A_{μ} \rightarrow \delta g^{αμ}= \delta g^{βν}=0$
$\bullet$ $2\int d^{4}x\sqrt{-g}F^{μν}(\nabla_{\mu}δA_{\nu} - \nabla_{\nu}δA_{\mu}) = 
4 \int d^{4}x\sqrt{-g}F^{μν}(\nabla_{\mu}δA_{\nu})$ since  $\nabla_{\nu}δA_{\mu} = - \nabla_{\mu}δA_{\nu}$ ($F$ is antisymmetric)
$\bullet$ $\int d^{4}x\sqrt{-g}\nabla_{\mu}(F^{μν}δA_{\nu}) =0$ : Boundary term. The variation of the field ($δΑ_{ν}$) vanishes at the boundary
A: Okay here we go. 
$$\begin{align} \mathcal{L} & = -\frac{1}{4} F^{\mu \nu} F_{\mu \nu}\\ &=-\frac{1}{4} g^{\mu \alpha} g^{v \beta} F_{\alpha \beta} F_{\mu \nu} \\ & = -\frac{1}{4} g^{\mu \alpha} g^{v \beta}\left(\partial_{\alpha} A_{\beta}-\partial_{\beta} A_{\alpha}\right)\left(\partial_{\mu} A_{v}-\partial_{\nu} A_{\mu}\right) \\ & = -\frac{1}{4} g^{\mu \alpha} g^{\nu \beta}\left(\left(\partial_{\alpha} A_\beta\right)\left(\partial_{\mu} A_{\nu}\right)-\left(\partial_{\alpha} A_{\beta}\right)\left(\partial_{\nu} A_{\mu}\right)-\left(\partial_\beta A_{\alpha}\right)\left(\partial_{\mu} A_{\nu}\right)+\left(\partial_\beta A_{\alpha}\right)\left(\partial_{\nu} A_\mu\right)\right). \end{align}$$
Using  $$\frac{\partial\left(\nabla_{\alpha} A_{\beta}\right)}{\partial\left(\nabla_{\mu} A_{\nu}\right)}=\delta_{\alpha}^{\mu} \delta_{\beta}^{\nu}$$ to compute the derivative,
$$\begin{align} \frac{\partial\mathcal{L}}{\partial(\partial_\omega A_\gamma)} & =  -\frac{1}{4}g^{\mu \alpha} g^{\nu \beta} ( \delta_{\alpha}^{\omega} \delta_{\beta}^{\gamma}\left(\partial_{\mu} A_{\nu}\right)+\left(\partial_{\alpha} A_{\beta}\right)\delta_{\mu}^{\omega} \delta_{\nu}^{\gamma} - \delta_{\alpha}^{\omega} \delta_{\beta}^{\gamma}\left(\partial_{\nu} A_{\mu}\right) - \left(\partial_{\alpha} A_{\beta}\right)\delta_{\nu}^{\omega} \delta_{\mu}^{\gamma} \\ & - \delta_{\beta}^{\omega} \delta_{\alpha}^{\gamma}\left(\partial_{\mu} A_{\nu}\right)-\left(\partial_{\beta} A_{\alpha}\right)\delta_{\mu}^{\omega}\delta_{\nu}^{\gamma}+ 
\delta_{\beta}^{\omega} \delta_{\alpha}^{\gamma}\left(\partial_{\nu} A_{\mu}\right) + \left(\partial_{\beta} A_{\alpha}\right)\delta_{\nu}^{\omega} \delta_{\mu}^{\gamma}) \\ & =-\frac{1}{4}( g^{\mu\omega}g^{\nu\gamma}\left(\partial_\mu A_\nu\right) + \left(\partial_\alpha A_\beta\right)g^{\omega\alpha}g^{\gamma\beta} - g^{\mu\omega}g^{\nu\gamma}\left(\partial_\nu A_\mu\right) - \left(\partial_\alpha A_\beta\right)g^{\gamma\alpha}g^{\omega\beta} \\ & -g^{\mu\gamma}g^{\nu\omega}\left(\partial_\mu A_\nu\right) - \left(\partial_\beta A_\alpha\right)g^{\omega\alpha}g^{\gamma\beta} + g^{\mu\gamma}g^{\nu\omega}\left(\partial_\nu A_\mu\right) + \left(\partial_\beta A_\alpha \right)g^{\gamma\alpha}g^{\omega\beta}) \\ & = -\frac{1}{4}\left(\partial^\omega A^\gamma + \partial^\omega A^\gamma - \partial^\gamma A^\omega - \partial^\gamma A^\omega - \partial^\gamma A^\omega - \partial^\gamma A^\omega+ \partial^\omega A^\gamma + \partial^\omega A^\gamma\right) \\ & = -\frac{1}{4}\left(4\partial^\omega A^\gamma - 4\partial^\gamma A^\omega\right) \\ & = \partial^\gamma A^\omega-\partial^\omega A^\gamma \\ & = F^{\gamma\omega}.
\end{align} $$
As you suggest in your question, $\frac{\partial\mathcal{L}}{\partial A_\gamma}=0$ since the field and its derivative are treated as independent coordinates in Lagrangian mechanics. So the Euler-Lagrange Equations read:
$$ 0=\partial_\omega\left(\frac{\partial\mathcal{L}}{\partial(\partial_\omega A_\gamma)}\right)-\frac{\partial\mathcal{L}}{\partial A_\gamma}=\partial_\omega F^{\gamma\omega}=0$$
which is Maxwell's equation in covariant form with zero current density. I'm sure there's an easier way to derive this but that's how I did it... Shoutout to Mathpix.
