# Darcy equation- problem in error and convergence rate

I am trying to solve this equation:
$$\mathbf{u} - \nabla p -(-y-\sin^2(\pi x)\sin(2\pi y), -x + \sin(2\pi x)\sin^2(\pi ))=0$$
$$\nabla \cdot \mathbf{u} + 2\pi\sin(\pi x)\sin(2\pi y)\cos(\pi x) - 2\pi\sin(2\pi x)\sin(\pi y)\cos(\pi y)=0$$

here is my code:

import psutil
import time
cpu_start = psutil.cpu_percent(interval=0.1)
start_time = time.time()

from mpi4py import MPI
from petsc4py import PETSc
import ufl
import numpy as np
import matplotlib.pyplot as plt

from basix.ufl import element, mixed_element
from dolfinx import fem, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
div, grad, inner, sin, cos, pi)
mesh_size =[]
error_u = []
error_p = []
for i in range(2,8):
# Create mesh and define function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 2**i, 2**i, mesh.CellType.triangle)
h = 1/2**i
k = 1  # Polynomial degree
U_el = element("RT", domain.basix_cell(), k)
# U_el = element("Lagrange",mesh.topology.cell_name(),2,shape = (mesh.geometry.dim, ))
P_el = element("P", domain.basix_cell(), k)
V_el = mixed_element([U_el, P_el])
V = fem.FunctionSpace(domain, V_el)

(u, p) = TrialFunctions(V)
(v, q) = TestFunctions(V)

# Define the analytical solution
x = SpatialCoordinate(domain)
p_exact_expr = x[0]*x[1] - 0.25
u_exact_expr = ufl.as_vector([-sin(pi*x[0])**2 * sin(2*pi*x[1]), sin(pi*x[1])**2 * sin(2*pi*x[0])])

# Define the source term
f1 = ufl.as_vector([-x[1] - sin(pi*x[0])**2 * sin(2*pi*x[1]), -x[0] + sin(2*pi*x[0]) * sin(pi*x[1])**2])
f2 = 2*pi*sin(pi*x[0])*sin(2*pi*x[1])*cos(pi*x[0]) - 2*pi*sin(2*pi*x[0])*sin(pi*x[1])*cos(pi*x[1])

# Weak form of the equations
dx = Measure("dx", domain)
a = inner(u, v) * dx + inner(p, div(v)) * dx + inner(div(u), q) * dx
L = inner(f1, v) * dx + inner(f2, q) * dx
bcs = []

# Solve the problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
w_h = problem.solve()
u_h, p_h = w_h.split()

# Create function spaces for each component
V_u = fem.FunctionSpace(domain, U_el)
V_p = fem.FunctionSpace(domain, P_el)

# Interpolate the exact velocity vector field
u_exact = fem.Function(V_u)
u_exact_expr = fem.Expression(u_exact_expr, V_u.element.interpolation_points())
u_exact.interpolate(u_exact_expr)

# Interpolate the exact pressure field
p_exact = fem.Function(V_p)
p_exact_expr = fem.Expression(p_exact_expr, V_p.element.interpolation_points())
p_exact.interpolate(p_exact_expr)

# Compute the errors
L2_error_u = fem.assemble_scalar(fem.form(inner(u_h - u_exact, u_h - u_exact) * dx))
L2_error_p = fem.assemble_scalar(fem.form((p_h - p_exact)**2 * dx))

# Take the square root to get the L2 norm
L2_error_u = np.sqrt(domain.comm.allreduce(L2_error_u, op=MPI.SUM))
L2_error_p = np.sqrt(domain.comm.allreduce(L2_error_p, op=MPI.SUM))

print(f"L2 error in velocity: {L2_error_u:.6e}")
print(f"L2 error in pressure: {L2_error_p:.6e}")
error_u.append(L2_error_u)
error_p.append(L2_error_p)
mesh_size.append(h)

plt.figure(figsize = (12,8))
plt.plot(mesh_size, error_u, label = 'Velocity error', marker = 'o')
plt.plot(mesh_size, error_p, label = 'Pressure error', marker = 's')
plt.xlabel('Mesh_size')
plt.ylabel('L2 Error')
plt.title('Error vs Mesh_size')
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.grid(True)
plt.show()

orders_u = []
orders_p = []
for i in range(1, len(mesh_size)):
r_u = np.log(error_u[i-1] / error_u[i]) / np.log(mesh_size[i] / mesh_size[i-1])
r_p = np.log(error_p[i-1] / error_p[i]) / np.log(mesh_size[i] / mesh_size[i-1])
orders_u.append(r_u)
orders_p.append(r_p)

# Print the orders of convergence
print("Order of convergence for u:", orders_u)
print("Order of convergence for p:", orders_p)

cpu_end = psutil.cpu_percent(interval=0.1)
elapsed_time = time.time() - start_time
print(f"Execution time: {elapsed_time} seconds")
print(f"CPU Usage: {cpu_end}%")


The code is running but the errors are very high and convergence rate is almost zero.