Sorry for the delay in answering. Since we answered question 1 in another thread, I will just focus on question 2 here. I will keep this fairly general, rather than focusing on the particular model in the question at hand, but I will assume that the coarse-graining method we are interested in is integrating out degrees of freedom with momenta greater than some cutoff. i.e., if our degrees of freedom are $m(\mathbf{q})$, we coarse grain by integrating out all modes with $|\mathbf{q}| > b\Lambda$, for $b < 1$ and $\Lambda$ the maximum wavelength.
I'll split this answer into a short version and a long version.
The short version is:
One normally doesn't attempt to compute the cumulant generating functional (CGF) for non-Gaussian field theories (at least, not that I have seen). Instead, once one has the "renormalized" field theory (i.e., obtained the coarse-grained and rescaled action retaining only the relevant interactions), one can identify the Feynman diagram rules and use those to systematically improve the mean-field estimates of the statistical moments. i.e., to calculate loop-corrections to the tree-level (Gaussian) approximations. Since this is normally done as a perturbative series expansion for each statistical moment, there's usually no useful $F[J]$ to write down, since it would basically just be written as a series
$$F'[J'] = \int d\mathbf{q}'~ J'(\mathbf{q}') \cdot \kappa_1(\mathbf{q}') + \int d\mathbf{q}_1' d\mathbf{q}_2'~ J'(\mathbf{q}_1')^T \kappa_2(\mathbf{q}_1',\mathbf{q}_2') J'(\mathbf{q}_2') + \dots,$$
where $\kappa_1(\mathbf{q}')$ and $\kappa_2(\mathbf{q}_1',\mathbf{q}_2')$ are the first and second order cumulants, respectively (the mean and the covariance) that would be calculated using Feynman diagrams, with $\dots$ indicating the higher order cumulants. I have used primes to indicate that these are the rescaled momenta (to connect to the notation in the longer answer). If all of the cumulants are only estimated up to tree-level, then in principle this series should sum up to $F_{\rm gaussian}$, and you could in principle try to organize the remaining terms (from the loop corrections) into an $F_{\rm corrections}$, but it would still be in the form of the series in $J'$ and up to whatever order of the loop approximation you've calculated things up to. I am not aware of a systematic perturbative method for calculating a global approximation to $F'[J']$. So, for the action in the question in particular, to calculate the moments you would typically work out the Feynman diagram rules for your coarse-grained action $S'[m'] = \beta H'[m']$ and use those to calculate the statistical moments.
This said, there is an approach called the "non-perturbative renormalization group" which could, in principle, be used to renormalize the CGF, although normally it focuses on the Legendre transform of the CGF, which is the average effective action $\Gamma[M(\mathbf{q}')]$, where $M(\mathbf{q}') \equiv \frac{\delta F'[J'(\mathbf{q}')]}{\delta J'(\mathbf{q}')}$ is the Legendre field conjugate to the source field $J'(\mathbf{q}')$. The average effective action also contains all of the information about the statistical moments. Even in this method, however, the goal is not usually to calculate an approximation for $\Gamma[M(\mathbf{q}')]$, but rather to calculate critical exponents or sometimes the "vertex functions" $\Gamma^{(n)}[\mathbf{q}'_1,\dots,\mathbf{q}'_n] \equiv \frac{\delta^n \Gamma[M(\mathbf{q}')]}{\delta M(\mathbf{q}'_1) \dots \delta M(\mathbf{q}'_n)}$, usually only for up to small $n$ and usually for the purposes of estimating the scaling form of the correlation functions near critical points. (Statistical moments can be obtained from the vertex functions). Both approaches have been applied to the $O(N)$ model, which is basically the model in the question. This paper reports some results using the non-perturbative methods, although it's rather technical. The key figure relevant to this discussion is Fig. 4, which plots the $\Gamma^{(2)}(p)/p^{2-\eta}$, where $\Gamma^{(2)}(p)$ is obtained from the 2-point vertex $\Gamma^{(2)}[\mathbf{q}_1,\mathbf{q}_2] = \Gamma^{(2)}(|\mathbf{q}'_1|) \delta(\mathbf{q}'_1-\mathbf{q}'_2)$ (the detlta function is due to translation invariance) and $p^{2-\eta}$ is the expected scaling of the function as $p \rightarrow 0$. Fig. 6 also plots the scaling function $g(x)$ obtained by this method, obtained from the 2-point correlation function $G^{(2)}(p) = \Gamma^{(2)}(p)^{-1} g(p\xi)$, with $\xi$ the correlation length near (but not at) criticality.
The longer version is:
All I really want to do here is add on some detail to support a couple of the statements made above and try to elucidate the relationship between the renormalization group calculation and the cumulant generating functional (CGF).
To start, consider the CGF $F[J(\mathbf{q})]$ for the fine-grained model, before we do any coarse graining:
$$e^{F[J(\mathbf{q})]} \equiv \int \mathcal D m(\mathbf{q})~e^{-S[m(\mathbf{q})] + \int d\mathbf{q}~J(\mathbf{q}) \cdot m(\mathbf{q})},$$
where $S[m(\mathbf{q})]$ is the action (equal to $\beta H[m(\mathbf{q}]$ in the question) and the sources $J(\mathbf{q})$ have been written explicitly in momentum space already. Suppose we could evaluate this integral exactly to obtain $F[J(\mathbf{q})]$. As you know, from this quantity we could obtain all statistical moments for the original fine-grained degrees of freedom $m(\mathbf{q})$ by functional differentiation.
Now, consider the coarse-grained action $S_b[m(\mathbf{q})]$ defined by integrating out modes with $|\mathbf{q}| > b \Lambda$:
$$e^{-S_b[\mathbf{q}]} \equiv \int \mathcal D m(|\mathbf{q}| > b\Lambda)~e^{-S[m(\mathbf{q})]},$$
where $\mathcal D m(|\mathbf{q}| > b\Lambda)$ is a shorthand to denote that we are only integrating out the high-momentum modes. Note that I have not performed the rescaling step yet.
Now, we can also write down the CGF $F_b[J(\mathbf{q})]$ for this coarse-grained action:
$$e^{F_b[J(\mathbf{q})]} \equiv \int \mathcal D m(|\mathbf{q}| \leq b \Lambda)~e^{-S_b[m(\mathbf{q})] + \int d\mathbf{q}~J(\mathbf{q}) \cdot m(\mathbf{q})},$$
where we integrate over the remaining modes with $|\mathbf{q}| \leq b \Lambda$.
Now we can ask: how is $F[J(\mathbf{q})]$ related to $F_b[J(\mathbf{q})]$? The answer is that
$$F_b[J(\mathbf{q})] = F[J(|\mathbf{q}| \leq b\Lambda),J(|\mathbf{q}| > b\Lambda) = 0].$$ That is, the coarse-grained cumulant generating functional is obtained (in this case) by simply setting the source terms $J(\mathbf{q})$ to zero for all sources with momenta $|\mathbf{q}| > b \Lambda$. The significance of this is that if you could calculate the full CGF $F[J]$ you can trivially obtain the coarse-grained CGF.
But what about the CGF after we also rescale the momenta and degrees of freedom? i.e., if we change variables to $\mathbf{q} \rightarrow b^{-1} \mathbf{q}'$ and $m(\mathbf{q}) \rightarrow z m'(\mathbf{q}')$, what is the corresponding CGF, $F'[J'(\mathbf{q}')]$? If we make this change of variables in our definition of $F_b[J(\mathbf{q})]$ above the action term will just become the rescaled action $S'[m'(\mathbf{q}')]$ (plus constant factors from the Jacobian, which we could alternatively absorb into the implicit normalization). So, we can focus on the source term, which becomes $\int d\mathbf{q}'~b^{-d} z J(b^{-1} \mathbf{q}') \cdot m'(\mathbf{q}')$. If we want our CGF to be the moments of the rescaled variables $m'(\mathbf{q}')$, then we expect the source-term in $F'[J'(\mathbf{q}')]$ ought to look like $\int d\mathbf{q}'~J'(\mathbf{q}') \cdot m'(\mathbf{q}')$, which motivates us to define
$$J'(\mathbf{q}') \equiv b^{-d} z J(b^{-1} \mathbf{q}').$$
i.e., the sources for the rescaled action are just rescalings of the sources for the coarse-grained action. Thus, if you could calculate the full CGF $F[J(\mathbf{q})]$ the CGF for the coarse-grained and rescaled theory is obtained by setting the appropriate source-terms to $0$, rescaling the remaining source terms, and then taking the limit of infinitely many iterations of the coarse-graining+rescaling step.
Ok, so what is the takeaway from this explanation so far? It's this: if we could calculate the full CGF for a theory we wouldn't need to renormalize it.
So what is renormalization doing for us here? Well, in the fully rescaled+coarse-grained action we expect that if we set the original bare parameters of the action to appropriate values (essentially, we tune the theory to the critical manifold), then as we repeatedly coarse-grain+rescale the model the so-called "irrelevant" interactions will be driven to zero and the relevant interactions will flow toward a scale-free fixed point that has lost memory of the initial fine-grained action. (It is very important to note that this is true only of the rescaled theory: if we were to only perform the coarse-graining steps and not the rescaling the irrelevant terms would not be suppressed and while we could still end up with a well-defined action it will depend on details of the original fine-grained theory and would include all interactions generated by coarse-graining, precluding writing down Feynman rules for calculating moments.).
The result is that we typically end up with a much simpler looking action for which we can write down Feynman rules and use perturbative methods for computing the statistical moments of the theory. For example, for the model in the question, the potential $U[m(\mathbf{q})]$ may be relatively complicated, but it turns out that only the quadratic and quartic interactions are relevant near the Gaussian fixed point in the renormalization group flow. We can write down Feynman rules for this simpler action and attempt to calculate the statistical moments (2-point functions, etc.), treating the quartic interaction as a perturbation. Up to artifacts introduced by approximations we make so that we can actually perform the renormalization group calculation, the results should in principle match what we would obtain by calculating the moments from the original fine-grained theory (after making the rescalings and taking the limit of infinitely many coarse-graining+rescaling steps).
This general picture doesn't change much if we use the non-perturbative renormalization group approach I mentioned in the short answer, except that it gives an alternate way to try and calculate the correlation functions for the renormalized action.
