Stochastic partial differential equations in fenics?


373
views
1
7 months ago by
Alex  
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

1 Answer


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

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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »