Solving stochastic PDEs in Fipy

Question

I would like to simulate coupled PDEs with a Gaussian white noise field, and was unable to find any examples or documentation that suggests how it should be done. In particular, I am interested in Cahn-Hilliard-like systems with noise:

d/dt(phi) = div(grad(psi)) + div(noise)

psi = f(phi) + div(grad(phi))

Is there a way to implement this in Fipy?

Answer 1

You can add a GaussianNoiseVariable as a source to your equations.

For a non-conserved field, you would do, e.g.,

noise = fp.GaussianNoiseVariable(mesh=..., mean=..., variance=...)
  :
  :
eq = fp.TransientTerm(...) == fp.DiffusionTerm(...) + ... + noise

for step in steps:
   noise.scramble()
      :
      :
   eq.solve(...)

For a conserved field, you would use:

eq = fp.TransientTerm(...) == fp.DiffusionTerm(...) + ... + noise.faceGrad.divergence