Counting number of solutions to a specific SAT formula
This is not an answer, just a long comment.

First, you say that you are unsure if this problem in in #P.  It depends on how you define your input, which can either be the $n$-by-$n$ graph or just $n$.  I would say that it is more natural to take the graph as input, in which case, an NP machine can guess $n^2$ bits (which is a polynomial, linear even, in the size of the input) and can check in polynomial time, linear even, if all your constraints are satisfied.  If the input is just $n$, then $n^2$ is exponentially larger than the input size $\log(n)$, so we can't say your problem is in #P, but it is in the counting class related to NEXP (in the same way that #P is the counting class related to NP).

Second, this input issue also affects what types of hardness we can expect to prove. If the input is the grid graph, then you have at most one input for each input length.  Such problems are said to be sparse (which means there is at most a polynomial number of inputs for each input length) and are not expected to be hard by (a generalization of) Mahaney's theorem.  But if you pick $n$ to be your input so that set of instances is no longer sparse, then your output is exponentially larger than your input.  What I am trying to say here is that the typical notions of hardness don't apply very well to your problem.

Third, I can say something about the possibility of this problem being tractable.  My research is in counting problems defined by local constraints, which we formally define as Holant problems. The difference between Holant problems and your problem is that we don't put such tight restrictions on the inputs graphs.  We allow restrictions like planar, bipartie, and prescribed degree (such at $k$-regular or at most degree $k$), none of which defines a sparse set of graphs.

The best I can do to cast your problem as a Holant problem is as follows.  For omitted definitions, see the preliminary section in a recent paper of mine.  In the Holant framework, edges take values from a finite domain and vertices impose constraints on the values assigned to its incident edges.  To (partially) model your problem, we consider the Boolean domain $\{0,1\}$.  As a first attempt to capture your problem, think of an $n$-by-$n$ grid graph as the input.  To associate a bit with each vertex $v$, we add an edge $e_v$ incident $v$ and to a new vertex $u_v$.  To communicate this bit value to the neighboring internal vertices, we duplicate all the edges between internal vertices (i.e. vertices not on the boundary) in the original $n$-by-$n$ grid graph.

The vertices on the boundary have no constraint, which is represented by the signature $[1,1]$ tensored with itself twice for the corner vertices and three times for the vertices on the sides.  The remaining interior vertices pick an edge in each direction and require these four edges and the other edge added above (with its associated bit) to take the same value. Then we impose your constraint by comparing this value with the four values assigned to the remaining four edges (one edge in each direction).  Of course, between each pair of internal vertices, they have to agree on which edge each uses to communicate their bit value.  One way to decide this, in a planar embedding of this graph, is that each vertex should use the edge that is encountered first in a clockwise traversal of its incident edges.

The paper of mine that I linked to above is a dichotomy for Holant problems that use symmetric signatures (which means that the output is invariant under any permutation of the input).  Unfortunately, the constraint that I used to capture your problem at internal vertices is not symmetric and it is still an open problem to prove the generalization of that dichotomy for general (i.e. not necessarily symmetric) signatures.  Worse yet, we currently do not even know how to decide if this Holant problem is easy or hard according to what we conjecture the dichotomy to be.  Furthermore, we can impose restrictions on the input graphs like being planar and bipartite that are satisfied by your grid graph, the hope being that this Holant problem may become tractable under these restrictions. However, a dichotomy for such Holant problems is years, maybe even a decade away.

In closing, let me emphasize that this Holant problem that I have described is "more difficult" than your problem (when the input is the grid graph).  The input to this (planar, bipartite) Holant problem can only use the complicated arity 9 constraint (that was assigned to the internal vertices) and the (unconstraining) unary constraint $[1,1]$ but the graph is not restricted to grid graphs. Formally speaking, your problem (when the input is the grid graph) reduces to this Holant problem. So, if this Holant problem were easy, then your problem (when the input is the grid graph) is easy too, but if this Holant problem is hard, then I want to say that your problem is hard, but I don't know how to do so (as I discussed above).