Functional data analysis enables to exploit the complexity and richness of data measured over continuous domains. When two functional random variables are available, it may be useful to determine their relation by means of a regression model. If the regression function is a linear Hilbert–Schmidt operator between two $L^2$ spaces, we are under the functional linear model with functional response. We propose a novel goodness-of-fit test for the null (composite) hypothesis of this model, against a general, unspecified alternative,leading to an omnibus test. The test statistic is formulated in terms of the quadratic norm over a doubly-projected empirical process, and is easy to compute, interpret and calibrate on its distribution via a wild bootstrap on the residuals.  A flexible hybrid approach involving LASSO regularization and linearly-constrained least-squares is used to perform the selection of the number of dimensions when estimating the residuals. The finite sample behaviour of the test, regarding power and size, is illustrated via a complete simulatio nstudy under varying scenarios. The test is applied to several real datasets to check the validity of the model.