Any ideas why we don't get quadratic convergence for this uncoupled Heat Equation system?

Question

The MMS works great for 2D but fails to show quadratic convergence in 1D - wondering why?

#Testing MMS on Heat Equation system - works well (in 2-D)!

from dolfin import *

from math import log

def Heat(n):
#Define mesh and geometry
#mesh = Rectangle(0,0,5,5,n,n)
mesh = Interval(n,0,1)

# Define function spaces
V = FunctionSpace(mesh, "CG" , 1) 
Y = MixedFunctionSpace([V, V])

# Time variables
dt = Constant(0.02); t = float(dt); T = 0.09


#a_expr = "1 + x[0]*x[0] + A*x[1]*x[1] + B*t" #exact solution 2D
#b_expr = "1 + x[0]*x[0] + A*x[1]*x[1] + B*t"

a_expr = "1 + A*x[1]*x[1] + B*t" #exact solution 1D
b_expr = "1 + A*x[1]*x[1] + B*t"

aexact = Expression(a_expr, A=3.0, B=1.2, t=0)
bexact = Expression(b_expr, A=3.0, B=1.2, t=0)

abexact = Expression((a_expr,b_expr), A=3.0, B=1.2, t=0)

#f = Constant(1.2 - 2. - 2*3.0) #MMS Source term 2D

f = Constant(1.2 - 2*3.0) #MMS Source term 1D


# Previous solution
gamma0 = interpolate(abexact,Y)
(alpha0,beta0) = split(gamma0)


#Define Function and TestFunctions
gamma = Function(Y); (alpha, beta) = split(gamma) #current solution
(v, w) = split(TestFunction(Y))


# Define variational form
F = (alpha*v + dt*inner(grad(alpha),grad(v)) - alpha0*v - dt*f*v + beta*w + dt*inner(grad(beta),grad(w)) - beta0*w - dt*f*w)*dx #uncoupled system 


#Boundary Conditions
bcs_a = DirichletBC(Y.sub(0), aexact, DomainBoundary())
bcs_b = DirichletBC(Y.sub(1), bexact, DomainBoundary())
bcs = [bcs_a, bcs_b]


# Solve problem
while t <= T:
    aexact.t = t
    bexact.t = t
    solve(F == 0 , gamma, bcs=bcs)
    gamma0.assign(gamma)
    t += float(dt)

alpha, beta = gamma.split()
error1  = errornorm(alpha,aexact,"L2",mesh = mesh)
error2  = errornorm(beta,bexact,"L2",mesh = mesh)


return error1

Err = [0,0]
Err[0] = Heat(50)
Err[1] = Heat(100)

print "errors:(%s,%s)" % (Err[0],Err[1])

print (log(Err[0]/Err[1],2))#, log(Err[1]/Err[2],2))

Comments

Please fix the code formatting, and post complete (but simple) code.

---

Hi Garth. Thanks for the quick response. It is greatly appreciated. Managed to fix the Heat equations system but here is the actual problem we were looking to solve - however it still does not converge quadratically. We are sure the variational form is correct so it seems something is wrong updating the solutions in time/imposing BCs/splitting functions...

#Solving Non-linear Schrodinger equation
from dolfin import *
from math import log


def GPE(n):
  # Define mesh and geometry
  mesh = Interval(n,0,1)

  # Define function spaces
  V = FunctionSpace(mesh, "CG" , 1) #where alpha,beta lives
  Y = MixedFunctionSpace([V, V])

  # Time variables
  dt = Constant(0.02); T = 0.09; t = float(dt)

  g_real = "(C*cos(B*x[0] + (A*A - B*B)*t))/cosh(A*x[0] - 2*A*B*t)"
  g_imag = "(C*sin(B*x[0] + (A*A - B*B)*t))/cosh(A*x[0] - 2*A*B*t)"

  gr = Expression(g_real, t=0, A = 1.0, B = 0.5, C = sqrt(2.0))
  gi = Expression(g_imag, t=0, A = 1.0, B = 0.5, C = sqrt(2.0))

  gri = Expression((g_real,g_imag), t=0, A = 1.0, B = 0.5, C = sqrt(2.0))


  #Previous solution
  gamma0 = interpolate(gri,Y)
  (alpha0,beta0) = split(gamma0)


  #Define Function and TestFunctions
  gamma = Function(Y); (alpha, beta) = split(gamma) #current solution
  (v, w) = split(TestFunction(Y))


  #'Practice' Constant
  K = Constant(1.0)

  #Define Trap Potential as harmonic
  P = Constant(0.0)


  # Define variational form
  F = (alpha*v + beta*w + dt*(inner(grad(alpha),grad(w)) - inner(grad(beta),grad(v))) + dt*P*(alpha*w - beta*v) - dt*K*(alpha**2 + beta**2)*(alpha*w - beta*v) - (beta0*w + alpha0*v))*dx


  #Boundary Conditions
  bcs_a = DirichletBC(Y.sub(0), gr, DomainBoundary())
  bcs_b = DirichletBC(Y.sub(1), gi, DomainBoundary())
  bcs = [bcs_a, bcs_b]


  # Solve problem
  while t <= T:
  gr.t = t 
  gi.t = t 
  solve(F == 0 , gamma, bcs=bcs)
  gamma0.assign(gamma)
  t += float(dt)

  alpha, beta = gamma.split()
  error1  = errornorm(alpha,gr,"L2",mesh = mesh)
  error2  = errornorm(beta,gi,"L2",mesh = mesh)

return error1


Err = [0,0,0]
Err[0] = GPE(50)
Err[1] = GPE(100)
Err[2] = GPE(200)


print "errors:(%s,%s,%s)" % (Err[0],Err[1],Err[2])

print (log(Err[0]/Err[1],2), log(Err[1]/Err[2],2))