Convergence rate in Poisson equation with fine_mesh solution

dolfinx
1

Dear FEniCS community,

I am trying to solve a simple Poisson equation in \Omega=(0,1)^2 where the left-hand side is constant and u=0 on \partial\Omega. To compute the convergence rates I solve the same problem in a refined mesh. I was expecting to observe a rate O(h^{k+1}) and O(h^k) for the L2 and H1-errors, respectively. However, for k\geq 2 the rate is stuck at 3 and 2. Am I missing something in the code? or is there a regularity estimate of the solution that have to be considered?.

The MWE is below.

import dolfinx.fem as fem
import dolfinx.fem.petsc
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import mesh, default_scalar_type
from mpi4py import MPI
from ufl import dx, grad, inner


def main(N, k):
    msh = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.triangle)

    tdim = msh.topology.dim
    fdim = tdim - 1
    msh.topology.create_connectivity(fdim, tdim)

    Vh = fem.functionspace(msh, element("Lagrange", msh.basix_cell(), k))
    msh.topology.create_connectivity(fdim, tdim)
    boundary_facets = mesh.exterior_facet_indices(msh.topology)
    boundary_dofs = fem.locate_dofs_topological(Vh, fdim, boundary_facets)
    bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, Vh)

    # Next, the variational problem is defined:
    u = ufl.TrialFunction(Vh)
    v = ufl.TestFunction(Vh)
    a = inner(grad(u), grad(v)) * dx
    L = inner(fem.Constant(msh, default_scalar_type(1)), v) * dx

    # +
    problem = fem.petsc.LinearProblem(
        a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )

    uh_lap = problem.solve()

    return uh_lap, Vh, msh


if __name__ == "__main__":
    results = []
    N = [2 ** (k + 2) for k in range(2, 6)]
    degree = 3
    hh, rU, rUH1, eU, eUH1 = [], [], [], [], []

    rU.append(0.0)
    rUH1.append(0.0)

    uhlap_fine, Vh_fine, _ = main(160, k=degree + 3)

    for nn in N:
        uh_lap, Vh, msh = main(nn, degree)
        tdim = msh.topology.dim
        num_cells = msh.topology.index_map(tdim).size_local + msh.topology.index_map(tdim).num_ghosts

        cells = np.arange(num_cells, dtype=np.int32)

        hmax = dolfinx.cpp.mesh.h(msh._cpp_object, tdim, cells).max()

        V_exact = fem.functionspace(msh, ("Lagrange", degree + 3))
        u_exact = fem.Function(V_exact)

        interpolation_data = fem.create_interpolation_data(V_exact, Vh_fine, cells)
        u_exact.interpolate_nonmatching(uhlap_fine, cells, interpolation_data=interpolation_data)
        u_exact.x.scatter_forward()

        # H01 errors
        diff = uh_lap - u_exact
        H1_diff = msh.comm.allreduce(
            fem.assemble_scalar(fem.form(inner(grad(diff), grad(diff)) * dx)),
            op=MPI.SUM,
        )
        H1_diff = np.sqrt(H1_diff)

        # L2 errors
        L2_diff = msh.comm.allreduce(
            fem.assemble_scalar(fem.form(inner(diff, diff) * dx)), op=MPI.SUM
        )

        L2_diff = np.sqrt(L2_diff)
        eU.append(L2_diff)
        eUH1.append(H1_diff)

        nk = N.index(nn)
        hh.append(hmax)

        if nk > 0:
            rU.append(np.log(eU[nk] / eU[nk - 1]) / np.log(hh[nk] / hh[nk - 1]))
            rUH1.append(np.log(eUH1[nk] / eUH1[nk - 1]) / np.log(hh[nk] / hh[nk - 1]))

print('L2-rate')
for k in range(len(N)):
    print(rU[k])
print('H1-rate')
for k in range(len(N)):
    print(rUH1[k])
2

Instead of interpolating the approximate solution to the finest mesh, you are doing the opposite, interpolating onto the coarser mesh.

You should move the coarser solution onto the finer grid.

3

Thanks for the response @dokken. By following your suggestion I have made the following modification:

import dolfinx.fem as fem
import dolfinx.fem.petsc
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import mesh, default_scalar_type
from mpi4py import MPI
from ufl import dx, grad, inner


def main(N, k):
    msh = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.triangle)

    tdim = msh.topology.dim
    fdim = tdim - 1
    msh.topology.create_connectivity(fdim, tdim)

    Vh = fem.functionspace(msh, element("Lagrange", msh.basix_cell(), k))
    msh.topology.create_connectivity(fdim, tdim)
    boundary_facets = mesh.exterior_facet_indices(msh.topology)
    boundary_dofs = fem.locate_dofs_topological(Vh, fdim, boundary_facets)
    bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, Vh)

    # Next, the variational problem is defined:
    u = ufl.TrialFunction(Vh)
    v = ufl.TestFunction(Vh)
    a = inner(grad(u), grad(v)) * dx
    L = inner(fem.Constant(msh, default_scalar_type(1)), v) * dx

    # +
    problem = fem.petsc.LinearProblem(
        a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )

    uhlap = problem.solve()

    return uhlap, Vh, msh


if __name__ == "__main__":
    results = []
    N = [2 ** (k + 2) for k in range(2, 6)]
    degree = 3
    hh, rU, rUH1, eU, eUH1 = [], [], [], [], []

    rU.append(0.0)
    rUH1.append(0.0)

    uhlap_fine, Vh_fine, _ = main(160, k=degree + 3)
    cells_fine_map = Vh_fine.mesh.topology.index_map(Vh_fine.mesh.topology.dim)
    num_cells_on_proc_fine = cells_fine_map.size_local + cells_fine_map.num_ghosts
    cells_fine = np.arange(num_cells_on_proc_fine, dtype=np.int32)

    for nn in N:
        uhlap, Vh, msh = main(nn, degree)
        tdim = msh.topology.dim
        num_cells = msh.topology.index_map(tdim).size_local + msh.topology.index_map(tdim).num_ghosts
        cells = np.arange(num_cells, dtype=np.int32)
        hmax = dolfinx.cpp.mesh.h(msh._cpp_object, tdim, cells).max()

        uh_coarse_to_fine = fem.Function(Vh_fine)
        interpolation_data = fem.create_interpolation_data(Vh_fine,Vh, cells_fine)
        uh_coarse_to_fine.interpolate_nonmatching(uhlap, cells_fine, interpolation_data=interpolation_data)
        uh_coarse_to_fine.x.scatter_forward()

        # H01 errors
        diff = uh_coarse_to_fine - uhlap_fine
        H1_diff = msh.comm.allreduce(
            fem.assemble_scalar(fem.form(inner(grad(diff), grad(diff)) * dx)),
            op=MPI.SUM,
        )
        H1_diff = np.sqrt(H1_diff)

        # L2 errors
        L2_diff = msh.comm.allreduce(
            fem.assemble_scalar(fem.form(inner(diff, diff) * dx)), op=MPI.SUM
        )

        L2_diff = np.sqrt(L2_diff)
        eU.append(L2_diff)
        eUH1.append(H1_diff)

        nk = N.index(nn)
        hh.append(hmax)

        if nk > 0:
            rU.append(np.log(eU[nk] / eU[nk - 1]) / np.log(hh[nk] / hh[nk - 1]))
            rUH1.append(np.log(eUH1[nk] / eUH1[nk - 1]) / np.log(hh[nk] / hh[nk - 1]))

print('L2-rate')
for k in range(len(N)):
    print(rU[k])
print('H1-rate')
for k in range(len(N)):
    print(rUH1[k])

and still observe 3 and 2 for k\geq 2.

4

Interestingly, if I replace your f= with: = 2 * pi**2 * sin(pi * x[0]) * sin(pi * x[1]), which satisfies the equation for u=\sin(\pi x)\sin(\pi y) (and satisfies the boundary condition) I get the expected convergence rates for k=3 and k=4.

For f=1, the solution to this problem is a quadratic function of the form u = ax(1-x) + by(1-y) (to ensure that your bcs hold). There is no point going beyond second order elements, as they resolve the problem to machine precision.

1 Like