Problem with solving poisson in primal form

Legacy DOLFIN
1

Hello,

I’ve been trying to test the convergence of a solver that I wrote based on the post here: https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/poisson/python/documentation.html#:~:text=The%20Poisson%20equation%20is%20the,gon%20ΓN.

Im using an exact solution that is polynomial, thus with a sufficiently high degree FEM space I should be able to recover it exactly.

Im posting the code I’m using

from dolfin import *
import numpy as np


def solver_online(A, mesh, degree, f):
    # Create mesh and define function space

    V = FunctionSpace(mesh, "Lagrange", degree)

    def boundary(x):
        return (
            x[0] < DOLFIN_EPS
            or x[0] > 1.0 - DOLFIN_EPS
            or x[1] < DOLFIN_EPS
            or x[1] > 1.0 - DOLFIN_EPS
        )

    # Define boundary condition
    u0 = Constant(0.0)
    bc = DirichletBC(V, u0, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    df = Expression(f, degree=degree)
    a = inner(Constant(A) * grad(u), grad(v)) * dx
    L = df * v * dx

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)
    return u, V


def test_primal_convergence(hs, degree, f, p_expression) :
    for h in hs:
        mesh = UnitSquareMesh(
            round(1 / (h * np.sqrt(2))),
            round(1 / (h * np.sqrt(2))),
        )

        print(f"{h = }")

        p, V = solver_online(
            np.array([[5, 1], [1, 5]]),
            mesh,
            degree,
            f,
        )

        p_exact = interpolate(
            Expression(p_expression, degree=degree),
            V,
        )
        set_log_level(50)

        p_err = np.sqrt(assemble((p - p_exact) ** 2 * dx(domain=mesh)))

        print(f"{p_err = }\n")


if __name__ == "__main__":
    hs = [0.25, 0.125, 0.0625]
    p_expression = "x[0]*(1-x[0])*x[1]*(1-x[1])"
    degree = 3
    f = "2*(1-2*x[0])*(1-2*x[1])-10*x[1]*(1-x[1])-10*x[0]*(1-x[0])"

    print("running primal formulation convergence test\n\n")
    test_primal_convergence(hs, degree, f, p_expression)

This is what I’m getting printed on the terminal:

running primal formulation convergence test


h = 0.25
Solving linear variational problem.
p_err = 0.06668253828972745

h = 0.125
p_err = 0.06666780805961044

h = 0.0625
p_err = 0.06666677073449483`
2

Looks like you haven’t implemented the correct source term (as the error remains constant). I am not at a computer, and cannot have a look at your code in detail, But i would suggest having a look at Error control: Computing convergence rates — FEniCSx tutorial

1 Like
3

I was forgetting about the minus sign:(

  • div A grad p = f

thank you very much!