Hi everyone,
I have a time dependent temperature solver that fails my mesh convergence test, namely, the L2 error I get (when compared to the exact solution) stays constant with decreasing mesh size. It looks like a simple problem, and I am comparing against an exact solution, so I am really lost as to why my test fails.
My MWE is:
from dolfin import *
import numpy as np
# Choose a mesh. You can also import a previously adapted mesh
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 40, 40)
n = FacetNormal(mesh)
t = Constant(0)
t_old = Constant(0) # used for time stepping the source code
tol = 1E-14
def full_boundary(x, on_boundary):
return on_boundary
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
super().__init__(kwargs)
def eval_cell(self, value, x, ufc_cell):
value[0] = x[0]*x[1] # exact solution evaluated at t=0.0
def value_shape(self):
return()
class GoverningEquations(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A)
for bc in self.bcs:
bc.apply(A)
# Define function space
S = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, S)
# Define Function for the previously converged step and TestFunction(s)
# w0 = Function(W)
dw = TrialFunction(W) # used for the bilinear form
(s1) = TestFunction(W) # used for the linear form
# constant velocity vector
u_vec = Constant((-1.0, 1.0))
theta = 0.5# parameter for the theta scheme in time
dt = 0.001 # timestep size
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
x = SpatialCoordinate(mesh)
facet_f = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # FACET function
# Create the measure
ds = Measure("ds", domain=mesh, subdomain_data=facet_f)
T0 = Function(W)
T = Function(W)
Texact = Expression('x[0]*x[1] + sin(2*pi*t)', t=t, element=S, domain=mesh)
Qexact = Expression('2*pi*cos(2*pi*t) -(x[0] - x[1])', t=t, element=S, domain=mesh)
Qexact_old = Expression('2*pi*cos(2*pi*t) - (x[0] - x[1])', t=t_old, element=S, domain=mesh)
TExact = Function(W)
T_init = InitialConditions(degree=2)
T.interpolate(T_init)
T0.interpolate(T_init)
# T^{k+theta}
T_mid = (1.0 - theta)*T0 + theta*T
Qexact_mid = (1.0 - theta)*Qexact_old + theta*Qexact
L0 = (1/dt) * (T*s1 - T0*s1)*dx() + (dot(u_vec, grad(T_mid)) * s1 + inner(grad(s1), grad(T_mid))) * dx()
L1 = - Qexact_mid*s1 * dx()
L = L0 + L1
J2 = derivative(L, T, dw)
# Step in time
tt = 0.0
T_final = 10*dt
bc = [DirichletBC(W, Texact, full_boundary)]
problem = GoverningEquations(J2, L, bc)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
W2 = FunctionSpace(mesh, S)
while(tt < T_final):
t_old.assign(tt)
tt += dt
print("tt is ", tt)
T0.vector()[:] = T.vector()
t.assign(tt)
solver.solve(problem, T.vector())
T_ex_pr = project(Texact, W2)
error_L2 = errornorm(T_ex_pr, T, 'L2')
# Print errors
print('error_L2 =', error_L2)
The error increases at each timestep, and it stays constant with decreasing mesh size, so I am clearly not solving the equation.
I have tried different manufactured solutions and decreasing dt, always with the same result. I am basing my code on the Cahn Hilliard demo . Any advice would be greatly appreciated!