Looking at fenics_pymc3, the package is tiny and the author just uses a few functions from
theano.tensor,theano.gradient,theano.graph.opandfecrto wrap around fenics. Do you think I would easily be able to find equivalent functions in aesara to do the same thing and try to use pymc v4?
I think so. A few import locations might have changed, but the functions and classes are still almost always the same ones. If you can’t find the correct function for something, feel free to ask.
Tightening priors is a great tip, and I don’t understand the reason for the
shape=(1,)thing, but it makes fenics_pymc3 work.
Ah, sounds to me like fenics_pymc3 just expects all parameters to be vectors? You could also call it as
T_preds = theano_fem_solver(k[None], A[None], B[None], C[None], Cs[None])
in that case, so that the strangeness stays close to the source. Won’t make much of a difference though, just keep your final trace object a bit cleaner. 
By the way, I also think changing the likelihood as described by @jessegrabowski here sounds like something that should be worth investigating.