Not getting expected rate of convergence for galerkin method w.r.t L2 error

General
1

Hello dr @dokken I am working on plain galerkin method and while computing rate of convergence w.r.t L2 error I am not getting order 2 convergence as expected for linear degree polynomials. Please let me know if I am doing something wrong. I am using dolfinx 0.8.0.0 version.
Here is the entire code:

import numpy as np
import matplotlib.pyplot as plt

from dolfinx.fem import (Constant, Function, functionspace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem, io, mesh
from dolfinx.mesh import locate_entities_boundary
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, inner
from petsc4py.PETSc import ScalarType

def solve_adv_diff(N=100, eps= 1, degree=1):

    # define mesh and function space
    msh = create_unit_square(MPI.COMM_WORLD,N,N)
    V = functionspace(msh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    x = SpatialCoordinate(msh)

    # define boundary conditions

    # left facet
    facet_2 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.isclose(x[0], 0))
    dof_2 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_2)
    bc_2 = dirichletbc(ScalarType(1), dof_2,V)

    # bottom left_facet
    facet_1 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]<=0.5))
    dof_1 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_1)
    bc_1 = dirichletbc(ScalarType(1), dof_1,V)

    # bottom right_facet
    facet_3 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]>0.5))
    dof_3 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_3)
    bc_3 = dirichletbc(ScalarType(0), dof_3,V)
    bcs = [bc_1,bc_2, bc_3]

    # define variational form
    vec_f = Constant(msh,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
    a = eps*inner(grad(u), grad(v))*dx + v*inner(vec_f,grad(u))*dx
    f = Constant(msh, ScalarType(0))
    L = f * v * dx

    # solving problem
    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    return problem.solve()

def error_L2(uh, u_ex, degree_raise = 3):
    # Create higher order function space
    tol = 1e-9
    degree = u_ex.function_space.ufl_element().degree
    family = "CG"
    mesh = u_ex.function_space.mesh
    W = functionspace(mesh, (family, degree + degree_raise))
    u_fine = Function(W)
    u_fine.interpolate(u_ex)
    u_coarse = Function(W)
    # Interpolate approx. solution in fine solution function space

    if uh.function_space.mesh != mesh:
        from dolfinx.fem import create_nonmatching_meshes_interpolation_data
        interpolation_data = create_nonmatching_meshes_interpolation_data(
            u_coarse.function_space.mesh._cpp_object,
            u_coarse.function_space.element,
            uh.function_space.mesh._cpp_object,tol)
        u_coarse.interpolate(uh, nmm_interpolation_data=interpolation_data)
    else:
        u_coarse.interpolate(uh)
    # writing fine and coarse solution for visualization
    from dolfinx import io
    with io.XDMFFile(mesh.comm, "u_fine.xdmf", "w") as file:
        file.write_mesh(mesh)
        file.write_function(u_fine)
    with io.XDMFFile(mesh.comm, "u_coarse.xdmf", "w") as file:
        file.write_mesh(mesh)
        file.write_function(u_coarse)

    # Compute the error in fine solution function space
    e_W = Function(W)
    e_W.x.array[:] = u_fine.x.array - u_coarse.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

degrees = [1]
Ns = [8,16,32,64]

for degree in degrees:
    Es = np.zeros(len(Ns))
    hs = np.zeros(len(Ns), dtype=np.float64)
    u_fine = solve_adv_diff(128,1,degree)
    for i, N in enumerate(Ns):
        u_approx = solve_adv_diff(N,1,degree)
        Es[i] = error_L2(u_approx, u_fine,degree_raise=3)
        hs[i] = 1. / Ns[i]
        print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}")
    rates = np.log(Es[1:] / Es[:-1]) / np.log(hs[1:] / hs[:-1])
    print(f"Polynomial degree {degree:d}, Rates {rates}")

here is the output :

image
image739×142 5.71 KB

2

Please don’t call out people by name demanding help. Those who provide help here are all volunteers and provide guidance as a service to the community. We are not on call for personalised help with esoteric research projects. If you require immediate advice the quickest path is to consult with your advisor or collaborators.

Regarding your problem. It looks like you are committing the variational crime that your boundary conditions do not conform with the regularity demanded by the H^1 CG space you’re employing (i.e. "Lagrange").

See for example that you impose u = 1 on the bottom left facets y=0, x \leq 0.5 and on the bottom right facets u = 0 on y = 0, x > 0.5. As a consequence I’d expect your convergence rates to be suboptimal. And this is indeed what you observe.

3

Dr. @nate I apologize for the misunderstanding and any offense caused by my post. It was not my intention to call anyone out or demand help. I am still learning to navigate things in fenicsx. I was just requesting help for the issue I have been facing from past several weeks. I truly appreciate your response and in the future I will refrain from tagging particular expert for requesting help.

1 Like