size is an expression so the result would be just some other expression.

We need more vector buffers for the tail loop!

The old formula contained some 'dim's instead of '2's,
which obviously works for 2D....

Tricky... some grad components depended on the sumfac kernels of the
other directions instead of their input...

... really elegant.

Need to use handler from pymbolic as loopy codegen handlers
are too involved for standalone usage.

until poisson_neumann passes