The Riemannian metric extends to a metric on all tensors of rank $(k,l)$, and in particular it extends to $p$-forms. The very definition of the Hodge star operator is that
$$
\omega \wedge \ast \eta = \langle \omega, \eta \rangle d\mathrm{vol}
$$
where $\langle \omega, \eta \rangle$ is the pairing on $p$-forms induced by the metric, and $d\mathrm{vol}$ is the volume form (or density) induced by the metric. So
$$
\int_M \omega \wedge \ast \eta = \int_M \langle \omega, \eta \rangle d\mathrm{vol}
$$
which is very obviously symmetric and positive definite.
If you are using some other definition of the Hodge star (for example using te $\epsilon$ symbol), then all you have to do is to check that it is equivalent to the condition $\omega \wedge \ast \eta = \langle \omega, \eta \rangle d\mathrm{vol}$.
Disclaimer: the comment section is overgrowing so here is an answer that I hope will erase all your doubts.
First part:
As said in the comment section, $(1)$ is a inner product on the vector space $\Lambda^p_xM$ while $(2)$ is an inner product on the vector space $\Omega^p(M)$. The link between them is that $(2)$ is obtained by integrating $(1)$ over the whole manifold $M$.
An analog is this: consider the set $M=[0,1]$. Then for each $x\in [0,1], T_x[0,1] = \mathbb{R}$ and one can define an inner product on $T_x[0,1] = \mathbb{R}$ by $\langle a,b \rangle_x = a\times b$. This is $(1)$.
A vector field on $[0,1]$ is just a smooth function $f:[0,1] \to \mathbb{R}$, and $(2)$ is here
$$
\langle f,g\rangle = \int_0^1 f(x) g(x) \mathrm{d}x = \int_0^1 \langle f, g\rangle_x \mathrm{d}x.
$$
Second part:
If $V$ is a vector space and $\langle\cdot,\cdot\rangle$ is an inner product, one can create a Riemannian metric on $V$, thought as a manifold the following way. As a vector space, the tangent bundle of $V$ is trivial:
$$TV = V\times V$$
and one can define the Riemannian metric $g_v = \langle\cdot,\cdot\rangle$ for $v\in V$. It is a constant Riemannian metric because the canonical trivialization makes $g_v$ to be a function independant of $v \in V$.
Take $V = \Omega^p(M)$ and $\langle \alpha,\beta\rangle = \int_M \alpha \wedge \star \beta$. Now, forget that $V$ and $\langle\cdot,\cdot\rangle$ is defined thanks to a Riemannian manifold $(M,g)$ and just look at its structure: it is a vector space with an inner product. Hence, for this inner product $\|\alpha\|$ is a number.
If you really want to think of this construction as a Riemannian manifold, like in the first paragraph, then $\|\alpha\|$ will be a function:
$$
\|\alpha\| : \beta \in \Omega(M)^p \mapsto \|\alpha\|(\beta) = \|\alpha\|\in \mathbb{R}
$$
which is constant and does not take points of $M$ as entries.
Comment: if you really do not understand what I said, here is just a question for you: for $x \in M$, how would you define $\left(\int_M \alpha\wedge \star \beta\right)(x)$?
This is the exact same thing as this question: how would you define $\left(\int_0^1 t^2 \mathrm{d}t\right)\left(\frac{1}{2}\right)$?
Best Answer
To cut a long story short, you don't use a Fréchet derivative but a Gateaux derivative, which doesn't run into any of the issues you bring up.
Anyhow, if $A$ is a fixed $1$-form on $N$, then for any $1$-form $\delta A$ of compact support in the interior of $N^4$, you can still define the variation $\delta S(A)(\delta A)$ of $S$ about $A$ in the direction of $\delta A$ by the Gateaux derivative $$ \delta S(A)(\delta A) = \left.\frac{d}{d\epsilon}\right\rvert_{\epsilon=0}S(A+\epsilon\delta A) = \lim_{\epsilon \to 0}\frac{S(A+\epsilon\delta A) - S(A)}{\epsilon}. $$ So, in this case, since $$ L(A+\epsilon\delta A) = -\frac{1}{2}d(A+\epsilon \delta A) \wedge \star d(A+\epsilon \delta A) -(A + \epsilon \delta A) \wedge \star J\\ = L(A) + \epsilon\left( -\tfrac{1}{2}dA \wedge \star d(\delta A) - \tfrac{1}{2}d(\delta A) \wedge \star dA -\delta A \wedge \star J\right) + O(\epsilon^2)\\ = L(A) + \epsilon\left(-d(\delta A) \wedge \star dA - \delta A \wedge \star J \right) + O(\epsilon^2)\\ = L(A) + \epsilon\left( -\delta A \wedge \left(d\!\star\!dA - \star J \right) -d\left(\delta A \wedge \star dA\right) \right) + O(\epsilon^2) $$ it follows that $$ \delta S(A)(\delta A) = \int_{N^4} \delta A \wedge \left(-d\!\star\!dA - \star J\right) - \int_{N^4} d(\delta A \wedge \star dA)\\ = \int_{N^4} \delta A \wedge \left(-d\!\star\!dA - \star J\right) = \int_{N^4} g\left(\delta A,-\Delta A - J\right)\,dvol, $$ where the Laplacian $\Delta$ on $1$-forms is defined (in Lorentzian signature) by $\Delta = \star d \! \star \! d$. Since the bilinear form $\langle A,B \rangle := \int_{N^4} g(A,B)\,dvol$ is still non-degenerate even in the Lorentzian case, it still follows that $\delta S(A)(\delta A) = 0$ for all suitable variations $\delta A$ if and only if the Euler–Lagrange equation $\Delta A = -J$ holds.