Euler-Lagrange equation in curved spacetime: how to calculate $\partial \mathcal{L}/\partial(\nabla_\nu A_\mu)$?

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.

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