### Stochastic partial differential equations in fenics?

373
views
1
7 months ago by
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
Community: FEniCS Project

4
7 months ago by
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
from fenics import *
import pylab as pl

# 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)
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
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")
​

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?
written 7 months ago by Alex