To keep the math simple, we can treat space as a very-fine-but-discrete lattice of very-large-but-finite dimensions, say of size $10^{1000}\times 10^{1000}\times 10^{1000}$ with a lattice spacing of $10^{-500}$ meters.
Let $\hat\psi_a(\vec x)$ denote the components of the fermion field operator, where $a$ is a spinor index, and let $\hat\psi_a^\dagger(\vec x)$ denote its adjoint. For a free field, the Hamiltonian has the form
$$
H\sim\sum_{\vec x} \sum_{a,b}\hat\psi_a^\dagger(\vec x)
D_{ab}\hat\psi_b(\vec x)
\tag{1}
$$
for some differential operator $D_{ab}$ whose details don't affect the essence of the answer. Since we've discretized space, we can streamline the notation by using a single index $n$ for both $a$ and $\vec x$. Then the Hamiltonian is
$$
H\sim\sum_{n,k}\hat\psi_n^\dagger
M_{nk}\hat\psi_k
\tag{2}
$$
where $M$ is an ordinary (but extremely large) matrix.
The goal is to find a representation of these operators as things that operate on a Hilbert space of "functionals" $\Psi[u]$, where $u$ are Grassmann variables, and then to derive an explicit expression for the functional that minimizes the energy defined by $H$.
The Hilbert-space representation
For each value of the index $n$, let $u_n$ be a Grassmann variable. Let $\Psi[u]$ denote a function of all of these variables, with complex coefficients. The Grassmann variables satisfy $u_n^2=0$, and we've replaced continuous space with a finite lattice, so the function $\Psi[u]$ is actually a polynomial with a finite (but very large) number of possible terms. We can think of $\Psi[u]$ as an ordinary many-component vector whose components are the coefficients of the polynomial. We can then define the inner product $\langle\Psi_1|\Psi_2\rangle$ in the obvious way, as the "dot product" of the vector of $\Psi_2$'s coefficients with the complex conjugate of the vector of $\Psi_1$'s coefficients.
The same inner product can also be written as
$$
\langle \Psi_1|\Psi_2\rangle
\sim
\int [dw][du]\
\exp\left(-\sum_n w_n u_n\right)
\Psi_1^*[w] \Psi_2[u]
\tag{3}
$$
using a second set of Grassmann variables $w_n$ that anti-commute with the original variables as well as with each other. (I"m writing "$\sim$" so I don't have to keep track of minus signs.) Here's the intuition behind equation (3): Grassmann integration is defined so that the result is simply the coefficient of the term in the integrand that involves exactly one factor of every integration variable. The exponential factor expands to a sum over all possible products of the factors $w_n u_n$. So, if we think of a polynomial $\Psi[u]$ as a vector of coefficients, then the integral over $[dw][du]$ is equivalent the usual "dot product" of two complex vectors $\Psi_1$ and $\Psi_2$. The exponential factor in (3) is what makes this work.
The field operators are represented by
$$
\hat\psi_n\Psi[u] = u_n\Psi[u]
\hskip2cm
\hat\psi_n^\dagger\Psi[u] \sim \frac{\partial}{\partial u_n}\Psi[u].
\tag{4}
$$
(Again, I'm writing "$\sim$" so that I don't have to keep track of minus signs.) Inside the inner product (3), taking the derivative of $\Psi_2[u]$ with respect to $u_n$ has the same effect as multplying the other wavefunction $\Psi_1[w]$ by $w_n$, so the second operator in (4) really does represent the adjoint of the first operator in (4).
For a more detailed treatment that pays closer attention to minus signs, see chapter 4 in Montvay and Münster's book Quantum Fields on a Lattice.
The ground state
Regardless of what Hilbert-space representation we use, we can choose a basis for the set of field operators $\hat\psi_n$ that diagonalizes the matrix $M$. This is standard textbook material, so suppose that this has already been done. Then the Hamiltonian has the form
$$
H\sim\sum_{n}\omega_n\hat\psi_n^\dagger\hat\psi_n.
\tag{5}
$$
(For example, such a basis exists in which each value of the index $n$ represents some combination of a spinor index $a$ and a momentum $\vec p$.) Then the ground state of the Hamiltonian (2) is annihilated by $\hat\psi_n$ if $\omega_n>0$ and is annihilated by $\hat\psi_n^\dagger$ if $\omega_n<0$.
Finally, in the representation described above, the wavefunction $\Psi[u]$ that satisfies these conditions is just a monomial, namely the product of the $u_n$ for which $\omega_n>0$. Here's why: the vacuum state is the eigenstate of $H$ with minimum eigenvalue. That condition is satisfied by the state that's annihilated by $\hat\psi_n$ whenever $\omega_n>0$ and by $\hat\psi_n^\dagger$ whenever $\omega_n<0$. In the basis that diagonalizes $H$, this implies that the state is just the product of $u_n$s for which $\omega_n>0$.
This monomial is the "wavefunction" of the ground state. Its spatial structure is encoded through the indices $n$, which refer to a basis that diagonalizes the Hamiltonian as in (5). Note that this is a different basis than the one used in (1) where $D$ is the usual Dirac operator. Equations (2)-(4) apply in any basis. Equation (1) corresponds to equation (3.84) in Peskin and Schroeder's An Introduction to Quantum Field Theory, and equation (5) corresponds to (3.90) in the same book. Going from (1) to (5) involves a change of basis, so the $\psi_n$ in equation (5) is a linear combination of the $\psi_n$s in equation (1); I recycled the notation. The details of this basis-change are shown in Peskin and Schroeder.
All of the information about the spatial struction of the "wavefunction" is buried in the relationship between the basis used in (1) and the basis used in (5). I omitted those details because that part is textbook material.
