### Stochastic partial differential equations in fenics?

373

views

1

Has anyone worked on linear/nonlinear stochastic PDEs in fenics? Is this possible to do with the current code or does it require new functions/elements?

Thank you,

Alex

Thank you,

Alex

Community: FEniCS Project

### 1 Answer

4

Several years ago I amused myself by solving a stochastic Poisson equation using the standard polynomial chaos expansion. I completely forget now the details and in general we add now a stochastic dimension \(p\) to the original problem: the 1d Poisson becomes \(p\)-dimensional.

The code should still work in the current stable version

Thank you. That's what I was looking for. Do you have any references/books on polynomial chaos expansion and how they are related with other methods like Euler-Maruyama?

The code should still work in the current stable version

```
from fenics import *
import pylab as pl
from scipy.integrate import quad
# Create mesh and define function space
mesh = UnitSquareMesh(100, 100)
element = FiniteElement("CG", mesh.ufl_cell(), 1)
_V = FunctionSpace(mesh, element)
def calculate_p(p):
element_V = MixedElement([element,]*p)
V = FunctionSpace(mesh, element_V)
# Define boundary condition
bc = DirichletBC(V, Constant([0.0, ]*p), "on_boundary")
# PCE matrix and vector
psi = Legendre()
P = pl.zeros([p, p])
Pv = pl.zeros([p, p])
for i in range(p):
for j in range(p):
def f(xi):
return psi.eval(i, xi)*(xi+2)*psi.eval(j, xi)*0.5
P[i, j] = quad(f, -1, 1)[0]
def f(xi):
return psi.eval(i, xi)*0.5
Pv[i, i] = quad(f, -1, 1)[0]
Pv = as_tensor(Pv)
P = as_tensor(P)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant([1.0, ]*p)
a = inner(P*grad(u), grad(v))*dx
L = inner(Pv*f, v)*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc, solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"})
u0 = Function(_V, name="Expectation")
assign(u0, u.sub(0))
sv = Function(_V, name="Std. variation")
for i in range(1, p):
def f(xi):
return psi.eval(i, xi)*psi.eval(i, xi)*0.5
c = quad(f, -1, 1)[0]
ui = Function(_V)
assign(ui, u.sub(i))
sv.vector()[:] += c*ui.vector()*ui.vector()
sv.vector()[:] = pl.sqrt(sv.vector().array())
sv.vector().apply("insert")
# Save solution in XDMF format
f = XDMFFile("p" + str(p) + "_mean.xdmf")
f.write(u0)
f = XDMFFile("p" + str(p) + "_sv.xdmf")
f.write(sv)
return u0, sv
def calculate_error(p):
u0, sv = calculate_p(p)
err_u0 = errornorm(u0, u0_ref, degree_rise=0)/u0_ref.vector().norm("l2")
err_sv = errornorm(sv, sv_ref, degree_rise=0)/sv_ref.vector().norm("l2")
return err_u0, err_sv
# Reference solution
u0_ref, sv_ref = calculate_p(10)
# Convergence rate
err_u0_list = []
err_sv_list = []
for p in range(2, 8):
err_u0, err_sv = calculate_error(p)
err_u0_list.append(err_u0)
err_sv_list.append(err_sv)
pl.semilogy(range(2, 8), err_u0_list, label=r"Error on $\mathrm{E}(u)$")
pl.semilogy(range(2, 8), err_sv_list, label=r"Error on $\mathrm{Std}(u)$")
pl.legend(loc="best")
pl.xlabel("PCE degree $p$")
pl.ylabel("Relative error / $p=10$")
pl.grid()
pl.savefig("conv_p.pdf")
```

written
7 months ago by
Alex

Please login to add an answer/comment or follow this question.