Yeah, that’s the jist of it. You get the output of the PDE, transform them into probabilities, and pass it into Bernoulli. There might be a logit_p parameter in pm.Bernoulli you could use instead of of p, to avoid having to call tt.nnet.sigmoid yourself (I think logit_p does some other stability tricks under the hood, but I don’t remember if it’s available in v3).
All I meant by “centered on the PDE output” was that you can take T_preds and do additional modeling with it. For example you could say that a fire at location (i, j) depends on 1) the predicted temperature at (i,j), and 2) the predicted temperatures in the 8 adjacent locations. So then probs could be a Multivariate Normal with mu=T_pred (that the “centering” I’m refering to), and the covariance matrix has the appropriate structure, maybe with just one diagonal variance and one off-diagonal adjacent variance parameter to estimate.
You could also combine T_pred with other features of locations (i, j) you might have.
Anyway, I guess none of that is helpful given your current troubles. Can you be more specific about the types of errors you are getting? I don’t know anything about fenics_pymc3 specifically, but checking the github repo it looks like its not actively maintained anymore. Have you looked into SunODE, which is compatible with PyMC 4.0 and actively maintained?