Let's take a look at the Feynman propagator for a massive scalar field:
$$D_F(x-y)=\int\frac{dp^3}{(2\pi)^3}\int\frac{dp^0}{2\pi}\frac{ie^{-ip \cdot (x-y)}}{p^2-m^2}$$
We can use this as the Green's function for a function $u(x)$ with a corresponding $f(y)$:
$$u(x)=\int f(y)G(x-y)dy$$
My question is about the $u(x)$ and the $f(y)$. How are these related to the fields? ARE they the fields? I mean, can we insert, in this case, the scalar field $\phi(x)$?
The Feynman propagator is the time-ordered two point correlation function of the field \begin{equation} \langle 0 | T\phi(x)\phi(y) | 0 \rangle = D_F(x,y) \end{equation}
Because $D_F$ obeys the equation for a Green's function for the Klein-Gordon equation (please double check signs and factors of $2$ and $i$) \begin{equation} \left(\square_x+m^2\right) D_F(x,y) = i \delta^{(4)}(x-y) \end{equation} we can use it to find solutions to the sourced Klein-Gordon equation. This is most useful in classical field theory.
In other words, given a source $J(x)$, we can solve the equation \begin{equation} \left(\square+m^2\right)\phi = J \end{equation} by \begin{equation} \phi(x)=-i\int d^4 y G_F(x,y) J(y) \end{equation} However, note that the Feynman propagator will generically yield a complex solution to the equations of motion, even for a real scalar field $\phi$. In fact within classical field theory, the Feynman propagator is a strange object; it requires a future boundary condition, and it will lead to solutions that have support outside the light cone of local disturbances $J$. So in classical field theory, we are almost always interested in the retarded propagator $D_R(x,y)$, which is causal, only requires past boundary conditions, and vanishes outside the light cone of $J$.
The reason that the Feynman propagator is useful in quantum field theory is because when solving the Dyson expansion for the $S$ matrix, we find we need to compute time ordered correlation functions of fields.
You were right in your comment that what I had suggested, i.e. choosing $u(x)=\phi(x)$ was trivial and could not be used to infer any useful results.
A more insightful choice would be $f(y)=\phi(y)$. Then, $$u(x)=\int dy \phi(y) G(x-y)$$ and by applying the Klein-Gordon operator on $u(x)$, one gets $$(\Box_x+m^2)u(x)=\int dy \delta^{(4)}(x-y)\phi(y)=\phi(x)$$ which is also not telling us much. It is just telling us that the KG operator is invertible.
However, you asked how are $u(x)$ and $f(y)$ associated with fields. The answer is simply "they do not", as the expression you write in your answer holds for any 2 functions of spacetime points. However, since $$u(x)=\int dy f(y) G(x-y)$$ holds for any two functions of spacetime points, this means that it holds if one chooses $f(y)$ to be the scalar field etc. I do not think that there is more to that, apart from the conclusion that the Green's function creates a particle associated with the scalar field at spacetime point $x$, after having destroyed the same kind of particle at spacetime point $y$.
I hope this answers your question in a satisfactory manner. If anything is still unclear, then let me know in the comments.
