Suppose I want to solve $u_{xy} = xy$ via Green's Function. This will correspond to the associated PDE $G_{xy} = \delta(x - x_G,\ y - y_G)$, and I want my boundary conditions for this Green problem to be $G(0,\ y) = 0,\ G_x(x,\ 0) = 0$.
The first step logically would actually be to find the adjoint of the differential operator and boundary conditions for $u$ to arrive at the boundary value problem for $G$, but I'll just assume everything is self-adjoint for now, and return to that step at the end. I solve $G_{xy} = \delta(x - x_G,\ y - y_G),\ G(0,\ y) = 0,\ G_x(x,\ 0) = 0$ using a double Laplace Transform. My double-transformed equation is $\hat{\hat G} = \frac{1}{rs}e^{-y_G s - x_G r}$, and unwinding the transforms, my Green's Function turns out to be $G = H[x - x_G]H[y - y_G]$, where $H[\ ]$ is the Heaviside Function.
The solution to the original PDE thus should be $u = \int_0^L \int_0^W x_G y_G H[x - x_G]H[y - y_G]\ dy_G\ dx_G$. Making use of the identity $H[a - b] = 1 - H[b - a]$, I can eliminate the Heaviside Functions by modifying the integration limits, and then straightforward calc computations yield $u = \frac{x^2 y^2}{4}$ (the domain dimensions $L$ and $W$ cancel).
If the operator/boundary conditions are indeed self-adjoint, and if I did everything correctly, then $u_{xy} = xy,\ u(0,\ y) = 0,\ u_x(x,\ 0) = 0$ should yield $u = \frac{x^2 y^2}{4}$. Solving this problem by simple integration and application of the boundary conditions confirms this answer, which in turn strongly confirms the self-adjointness hypothesis as well as my computations so far.
Returning to the beginning, all that remains is to find the adjoint to verify this. I have $u_{xy} = xy$, and I take the inner product of both sides with $G$. So I have to use integration by parts on $\int_0^L \int_0^W u_{xy}G\ dy\ dx$ with respect to $y$. This gives $\int_0^L (u_x G |_{y = 0}^{y = W} - \int_0^W u_x G_y\ dy) dx$. After distributing the outer integral, swapping the order of integration in the double integral, and applying integration by parts again, this time with respect to $x$, I get $\int_0^L u_x G |_{y = 0}^{y = W}\ dx - \int_0^W u G_y |_{x = 0}^{x = L}\ dy + \int_0^L \int_0^W uG_{xy}\ dy\ dx$.
The final term is indeed the original inner product with the differential operator moved from $u$ to $G$. I just have to check the boundary conditions now; I plug $u(0,\ y) = 0$ and $u_x(x,\ 0) = 0$ into $\int_0^L u_x G |_{y = 0}^{y = W}\ dx - \int_0^W u G_y |_{x = 0}^{x = L}\ dy$, and am left with $\int_0^L u_x(x,\ W)G(x,\ W)\ dx - \int_0^W u(L,\ y)G_y(L,\ y)\ dy$.
However, my understanding is that I now have to pick two boundary conditions in $G$ that make this expression $0$, and all I have left to use are $G(x,\ W)$ and $G_y(L,\ y)$. Setting these to $0$ won't give me the $G(0,\ y) = 0$ and $G_x(x,\ 0) = 0$ boundary conditions I expected. This seems to suggest that the operator $\frac{\partial^2}{\partial y\ \partial x}$ is actually not self-adjoint. But this doesn't make sense either, because when I assumed self-adjointness, I got the correct answer, and when I use simple integration to solve $u_{xy} = xy$ and apply the boundary conditions in $u$ which would produce the needed boundary conditions in $G$ (they are $u(x,\ W) = 0$ and $u_y(L,\ y) = 0$, which I can find by repeating the entire adjoint computation swapping the initial order of integration), I get $u = \frac{x^2 y^2 + L^2 W^2 - x^2 W^2 - L^2 y^2}{4}$, which no longer matches my Green's Function method solution $u = \frac{x^2 y^2}{4}$.
What am I doing wrong here? Is there some simple mistake in my computations, or do I understand the procedure of finding the adjoint operator/boundary conditions incorrectly?
You have $G(x,y; x_G, y_G) = H[x - x_G]\cdot H[y - y_G]$ and $$ \int_I x_G y_G\cdot G(x,y;x_G, y_G) \,\text dx_G \text dy_G = \frac14 x^2y^2, \quad (x,y)\in I = [0, L]\times[0, W] $$ You are supposed to multiply $G(s, t; x, y)$ at two sides of the original equation, and intergrate over $(x,y) \in I$ $$ \text{RHS} = \int_I x y\cdot G(s, t; x, y)\,\text dx \text dy = \frac14 s^2t^2, \quad (s,t)\in I $$ And for $\text{LHS}$, $$ \begin{align} \text{LHS} &= \int_I \partial_x\partial_y u(x,y) \cdot G(s, t; x,y) \,\text dx \text dy \\ &= \cdots \\ &= \int_0^L \partial_x u(x,W)\cdot G(s, t; x, W)\,\text dx - \int_0^W u(L,y)\cdot \partial_yG(s,t;L,y) \,\text dy\, + u(s, t) \end{align} $$ in which, $$ G(s,t; x, W) = H[s - x]\cdot H[t - W] = 0, \quad \text{since } t < W\\ \partial_y G(s,t; L, y) = H[s-L]\cdot (-\delta(t-y)) = 0, \quad\text{since } s < L $$ So, $$ \text{LHS} = u(s,t) = \frac14 s^2t^2 = \text{RHS} $$