Custom stencil operations

Let me start by saying this is a help request; not as issue as such. If there is a more appropriate place for those, I am all ears.

On a high level; I want to solve a 3d laplacian with non-constant coefficients. That means I cannot use any off-the-shelf optimised convolutional libraries since they all assume constant stencils per definition. I currently have this implemented using first order differential operators as div(non_constant_field*grad(T)); but thats quite inefficient, and I want to fuse that loop into a single stencil operation. Without reinventing the wheel of how to do efficient stencil operations on the GPU.

It seems to me that loopy is the right tool for the job. But I am not quite confident as to the best way to go about it. I believe the most optimal implementation is a wavefront of threads looping through the array, to make the best use of local memory; whereby each thread processes multiple array elements. But perhaps there are competitive methods that map better to the loopy paradigm?

Having examples of this nature in the project strikes me as very helpful; and id be happy to add a PR with such examples if I manage to figure out something worthwhile. That is, if I am not reinventing the wheel in the first place, since I could imagine people have used loopy before with exactly these type of use cases in mind.