Suppose I have an arbitrary PDE solver I’ve built. This solver solves simple linear problems that look like Ax=b.

Now suppose I specify a *distribution* of source terms (the b vector). This could be a straightforward probability distribution, or maybe a distribution described by a data-driven model. Has anyone tried to leverage Julia’s type system to cleanly propagate that source distribution through the PDE solver, such that you also get a solution distribution?

To be clear, I’m not looking for a way to Monte Carlo *sample* this problem; I’m interested in running a single solve (if possible). I’d also like to avoid a *symbolic* approach (I’m using a discretized PDE solver).