In this paper, Wald presents a quite general construction of a sympletic form for classical fields. If I understood (which I might have not, and in that case corrections are highly appreciated), the procedure is somewhat along these lines:
Let a spacetime $(M,g)$ be given and let $\mathcal{L}$ be a Lagrangian density specifying a field theory in $(M,g)$. The action is $$S[\phi]=\int_M \mathcal{L}[\phi](x) \sqrt{g} \ d^4x.$$
The equations of motion of the theory follow from $\delta S =0$. It turns out that when we compute $\delta S$ we get a boundary term. So the full variation is $$\delta S = \int \bigg[\bigg(\dfrac{\partial \mathcal{L}}{\partial \phi}+\nabla_a \dfrac{\partial\mathcal{L}}{\partial (\nabla_a\phi)} \delta\phi\bigg)+\nabla_a\bigg(\dfrac{\partial \mathcal{L}}{\partial(\nabla_a \phi)}\delta \phi\bigg) \bigg]\sqrt{g} d^4x $$
The boundary term is usually discarded when computing the equations of motion out of the argument that everything vanishes at infinity.
Nonetheless, Wald calls this thing the presympletic current potential density $$\Theta^a= \dfrac{\partial \mathcal{L}}{\partial (\nabla_a \phi)}\delta\phi$$
The phase space $\Gamma$ of the theory is the manifold of solutions to the equations of motion. A tangent vector is a variation $\delta \phi$ of a solution $\phi$, which is itself a function.
In the above, $\Theta^a$ is a map defined at each $\phi\in \Gamma$ and is a linear function of variations $\delta \phi \in T_\phi \Gamma$. It takes values on vector fields defined on spacetime.
One can then define the presympletic current to be a map $\omega^a$ defined on $T_\phi\Gamma\times T_\phi\Gamma$ taking values on vector fields on spacetime by $$\omega^aa[\phi,\delta_1 \phi,\delta_2\phi]=\delta_1 \Theta^a[\phi,\delta_2 \phi]-\delta_2\Theta^a[\phi,\delta_1\phi]$$
What we mean is: compute $\Theta^a[\phi,\delta_2\phi]$ and then vary this considering that the field varies by $\delta \phi_1$, and the same for the other term.
Define the sympletic form by integrating the pre-sympletic current over a Cauchy surface $$\Omega_\phi(\delta_1 \phi,\delta_2\phi)=\int_{\Sigma} \omega^a[\phi,\delta_1\phi,\delta_2\phi] d\Sigma_a$$ By construction this is, at each $\phi$ a bilinear function of two variations, which is skew symmetric and takes values on $\mathbb{R}$. Thus it is a two-form on $\Gamma$ as a sympletic form would be.
As a recipe this is fine. Indeed this yields the expected sympletic forms.
But how on Earth would anyone think of doing this? I mean, the main reasons I'm unconfortable with it, is that I look at it trying to see how anyone would think of doing all this to get a sympletic form for a classical field theory and I simply can't get it.
So what is the intuition behind this construction? What is the motivation to construct the sympletic form for classical field theories like that? How would one get the idea of doing this?