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 :