I am trying to solve the following equation
$$ \nabla\cdot(\sigma(\nabla\boldsymbol{u} + I)) = 0 $$
on a unit square with $ \boldsymbol{u} $ periodic at the boundary and $I$ denoting the identity matrix. I used the method of manufactured solutions to check my code, so chose $\sigma = 1$ and $\boldsymbol{u} = (-\cos(2\pi y),-\sin(2\pi x))$. When comparing my solver with the exact solution, I find that the error varies strangely with mesh size, for example, if I choose a (32,32) mesh, my error is 0.0000204, but when I choose (64,64), my error is 0.002119. Interestingly, when I consider a mesh size of (256,256), the error is down to 0.000235.
Could anybody suggest why this might be? Thank you!
from dolfin import *
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], 0) and near(x[1], 1)) or
(near(x[0], 1) and near(x[1], 0)))) and on_boundary)
def map(self, x, y):
if near(x[0], 1) and near(x[1], 1):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
elif near(x[0], 1):
y[0] = x[0] - 1.
y[1] = x[1]
else:
y[0] = x[0]
y[1] = x[1] - 1.
k = 5
n = 2**k
mesh = UnitSquareMesh(n,n)
V = VectorFunctionSpace(mesh, "CG", 1, constrained_domain = PeriodicBoundary())
u = TrialFunction(V)
v = TestFunction(V)
s = Expression("1.0")
uex = Expression(("-cos(2*DOLFIN_PI*x[1])","-sin(2*DOLFIN_PI*x[0])"), element=V.ufl_element())
t = Expression(("4*DOLFIN_PI*DOLFIN_PI*cos(2*DOLFIN_PI*x[1])","4*DOLFIN_PI*DOLFIN_PI*sin(2*DOLFIN_PI*x[0])"), element = V.ufl_element())
a = -s*inner(grad(u), grad(v))*dx
L = +s*div(v)*dx + inner(t, v)*dx
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs=[])
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "gmres"
solver.parameters["preconditioner"] = "hypre_amg"
solver.parameters["krylov_solver"]["monitor_convergence"] = False
solver.parameters["krylov_solver"]["relative_tolerance"] = 1e-10
solver.solve()
file = File("periodic_trial_%i.pvd" % k)
file << u
file2 = File("periodic_exact.pvd")
u_exact =interpolate(uex, V)
file2 << u_exact
print "error %f" % assemble(inner(u-u_exact, u-u_exact) * dx)
