Solving PDEs over distributions

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).

There are a few ways to do it.

Measurements.jl is the simplest but only linear i.e. propagates normal distribution approximates, which may be too crude for many applications. MonteCarloMeasurements is particle based and no faster than solving N times, but is a nice interface for it. Polynomial chaos expansions are quite a good approach for PDEs and strike a nice balance between nonlinearity and performance, but are a bit more advanced in API as well.